text:
#################################################################################################
# AIM 1: ROC
#################################################################################################
## Deps
using Distributed
using CodecZlib
using StatsPlots
using Distributions
using DelimitedFiles
using Plots.PlotMeasures

## Constants
const noise_levels = ["$(s)" for s in collect(2.0:0.5:3.5)]
const data_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Simulations/Aim-1"

# Thresholds
const nano_thrs = vcat(collect(-500.0:10.0:-20.0), collect(-20.0:0.1:20.0), collect(20.0:10.0:500.0))

# Colors
const call_thrs = Dict("Nanopolish" => nano_thrs)
const cllr_color_code = Dict("Nanopolish" => "#D55E00")
const blnd_col = ["#999999","#E69F00","#56B4E9","#009E73","#F0E442","#0072B2","#D55E00","#CC79A7"]
const s_color_cde = Dict("0.5" => blnd_col[1], "1.0" => blnd_col[2], "1.5" => blnd_col[3], "2.0" => blnd_col[4],
                         "2.5" => blnd_col[5], "3.0" => blnd_col[6], "3.5" => blnd_col[7], "4.0" => blnd_col[8])

## Default attributes
default(titlefont=(14, "arial"),guidefont=(16, "arial"),tickfont=(12, "arial"))

# Function to threshold calls
@everywhere function thresh_calls(conf_x, thrs)

    # Init output vectors
    out_conf_x = []

    # Compare confidence score to threshold to determine call
    @inbounds for i = 1:length(conf_x)

        if conf_x[i] == "n/a"
            # Fails to detect
            push!(out_conf_x, 0)
        elseif conf_x[i] < thrs
            # Negative
            push!(out_conf_x, -1)
        elseif conf_x[i] >= thrs
            # Positive
            push!(out_conf_x, 1)
        else
            println("Something went wrong when thresholding calls ...")
        end

    end

    return out_conf_x

end

# Function to compute sensitivity (recall, or true positive rate)
@everywhere function comp_sens(truth, pred)

    # Get positives
    pos_ind = pred .== 1

    # Get true positives
    if length(truth[pos_ind]) != 0
        true_pos = sum(truth[pos_ind] .== pred[pos_ind])
    else
        return 0.0
    end

    # Return sensitivity
    return true_pos / sum(truth .== 1)

end

# Function to compute specificity (selectivity or true negative rate)
@everywhere function comp_spec(truth, pred)

    # Get negatives
    neg_ind = pred .== -1

    # Get prediction & truth
    if length(truth[neg_ind]) != 0
        true_neg = sum(truth[neg_ind] .== pred[neg_ind])
    else
        return 0.0
    end

    # Return specificity
    return true_neg / sum(truth .== -1)

end

@everywhere function pmap_thresh(true_x, conf_x, t)

    # Threshold calls
    thresh_x = thresh_calls(conf_x, t)

    # Return
    return comp_sens(true_x, thresh_x), 1.0 - comp_spec(true_x, thresh_x)

end

function get_tpr_fpr(caller, s)

    # Read in data
    stream = GzipDecompressorStream(open("$(data_dir)/$(caller)/gm12878_chr22_sigma_$(s)_$(caller)_tuples_sample.tsv.gz"))
    in_data = readdlm(stream)
    close(stream)

    # Remove NaN
    kp_ind = .! map(isequal("n/a"), in_data[:,2])
    in_data = in_data[kp_ind,:]

    # Get true calls and confidence scores for predictions
    true_x = in_data[:,1]
    conf_x = in_data[:,2]

    # Get tpr and fpr
    pmap_out = pmap(t -> pmap_thresh(true_x, conf_x, t), call_thrs[caller])

    # Disentangle data
    tpr_vec = [x[1] for x in pmap_out]
    fpr_vec = [x[2] for x in pmap_out]

    # Return tuple
    return tpr_vec, fpr_vec

end

function get_auc(fpr_vec, tpr_vec)

    # Integrate ROC curve
    auc = 0.0
    @inbounds for i in 1:(length(fpr_vec) - 1)
        x = fpr_vec[i] - fpr_vec[i + 1]
        y = tpr_vec[i]
        auc += x * y
    end

    # Return AUC
    return auc

end

#################################################################################################
# ROC & AUC
#################################################################################################

# Init plots
p = plot(ylabel="sensitivity", xlabel="1-specificity", title="ROC - Nanopolish", xlim=(0, 1));

# Init AUROC array
auc_nano = []

# Loop over noise levels
for s in noise_levels

    # Print noise level
    println("Working on sigma=$(s)")

    # Get TPR & FPR
    tpr_vec, fpr_vec = get_tpr_fpr("Nanopolish", s)

    # Plot curve
    col = s_color_cde[s]
    plot!(p, fpr_vec, tpr_vec, seriestype=:line, alpha=0.5, color=col, label="sd=$(s)")

    # Calculate AUC
    push!(auc_nano, get_auc(fpr_vec, tpr_vec))

end

# Plot ROC ROC_Curves
plot(p,size=(700, 600))
savefig("$(data_dir)/ROC-Curves.pdf")

# Plot AUROC
col = cllr_color_code["Nanopolish"]
p4 = plot(xlabel="signal noise level (sd)", ylabel="AUC-ROC", ylim=(0.5, 1.0));
plot!(p4, noise_levels, auc_nano, seriestype=:scatter, markershape=:circle, color=col, label="")
plot!(p4, noise_levels, auc_nano, seriestype=:line, alpha=0.5, color=col, label="")
plot(p4,size=(700, 600))
savefig("$(data_dir)/AUC-ROC.pdf")


text:
import StatsBase: fit, predict
using SimilaritySearch:
    Sequential, KnnResult, search, empty!, Item, Index

#using TextModel

export kmap, centroid!, partition, knearestreferences, sequence, fit, predict

include("criterions.jl")
include("fftraversal.jl")
include("dnet.jl")
include("kclass.jl")
include("oneclass.jl")

"""
Transforms `objects` to a new representation space induced by ``(refs, dist, kernel)``
- `refs` a list of references
- `kernel` a kernel function (and an embedded distance) with signature ``(T, T) \\rightarrow Float64``
"""
function kmap(objects::AbstractVector{T}, kernel, refs::AbstractVector{T}) where {T}
    # X = Vector{T}(length(objects))
    m = Vector{Vector{Float64}}(undef, length(objects))
    @inbounds for i in 1:length(objects)
        u = Vector{Float64}(undef, length(refs))
        obj = objects[i]
        for j in 1:length(refs)
            u[j] = kernel(obj, refs[j])
        end

        m[i] = u
    end

    return m
end

"""
Groups items in `objects` using a nearest neighbor rule over `refs`.
The output is controlled using a callback function. The call is performed in `objects` order.

- `callback` is a function that is called for each `(objID, refItem)`
- `objects` is the input dataset
- `dist` a distance function ``(T, T) \\rightarrow \\mathbb{R}``
- `refs` the list of references
- `k` specifies the number of nearest neighbors to use
- `indexclass` specifies the kind of index to be used, a function receiving `(refs, dist)` as arguments,
    and returning the new metric index

Please note that each object can be related to more than one group ``k > 1`` (default ``k=1``)
"""
function partition(callback::Function, dist::Function, objects::AbstractVector{T}, refs::Index; k::Int=1) where T
    res = KnnResult(k)
    for i in 1:length(objects)
        empty!(res)
        callback(i, search(refs, dist, objects[i], res))
    end
end

"""
Creates an inverted index from references to objects.
So, an object ``u`` is in ``r``'s posting list iff ``r``
is among the ``k`` nearest references of ``u``.

"""
function invindex(dist::Function, objects::AbstractVector{T}, refs::Index; k::Int=1) where T
    π = [Vector{Int}() for i in 1:length(refs.db)]
    # partition((i, p) -> push!(π[p.objID], i), dist, objects, refs, k=k)
    partition(dist, objects, refs, k=k) do i, res
        for p in res
            push!(π[p.objID], i)
        end
    end
    π
end

"""
Returns the nearest reference (identifier) of each item in the dataset
"""
function sequence(dist::Function, objects::AbstractVector{T}, refs::Index) where T
    s = Vector{Int}(length(objects))
    partition(dist, objects, refs) do i, res
        s[i] = first(res).objID
    end
    s
end

"""
Returns an array of k-nearest neighbors for `objects`
"""
function knearestreferences(dist::Function, objects::AbstractVector{T}, refs::Index) where T
    s = Vector{Vector{Int}}(length(objects))
    partition(dist, objects, refs) do i, res
        s[i] = [p.objID for p in res]
    end
    s
end

"""
Computes the centroid of the list of objects

- Use the dot operator (broadcast) to convert several groups of objects
"""
function centroid!(objects::AbstractVector{Vector{F}})::Vector{F} where {F <: AbstractFloat}
    u = copy(objects[1])
    @inbounds for i in 2:length(objects)
        w = objects[i]
        @simd for j in 1:length(u)
            u[j] += w[j]
        end
    end

    f = 1.0 / length(objects)
    @inbounds @simd for j in 1:length(u)
        u[j] *= f
    end

    return u
end

ext: jl


text:
"""
$(TYPEDEF)
"""
abstract type BasicEnsembleAlgorithm <: EnsembleAlgorithm end

"""
$(TYPEDEF)
"""
struct EnsembleThreads <: BasicEnsembleAlgorithm end

"""
$(TYPEDEF)
"""
struct EnsembleDistributed <: BasicEnsembleAlgorithm end

"""
$(TYPEDEF)
"""
struct EnsembleSplitThreads <: BasicEnsembleAlgorithm end

"""
$(TYPEDEF)
"""
struct EnsembleSerial <: BasicEnsembleAlgorithm end

#=
if (kwargs[:parallel_type] == != :none && kwargs[:parallel_type] == != :threads)
  error("Distributed arrays cannot be generated via none or threads")
end
(batch_size != trajectories) && warn("batch_size and reductions are ignored when !collect_result")

elapsed_time = @elapsed u = DArray((trajectories,)) do I
    solve_batch(prob,alg,kwargs[:parallel_type] ==,I[1],pmap_batch_size,kwargs...)
end
return EnsembleSolution(u,elapsed_time,false)
=#

function __solve(prob::AbstractEnsembleProblem,
                 alg::Union{DEAlgorithm,Nothing};
                 kwargs...)
    #=
    if alg isa EnsembleAlgorithm
      @error "You forgot to pass a DE solver algorithm! Only a EnsembleAlgorithm has been supplied. Exiting"
    end
    =#
    if :parallel_type ∈ keys(kwargs)
      if kwargs[:parallel_type] == :none
        ensemblealg = EnsembleSerial()
      elseif kwargs[:parallel_type] == :pmap || kwargs[:parallel_type] == :parfor
        ensemblealg = EnsembleDistributed()
      elseif kwargs[:parallel_type] == :threads
        ensemblealg = EnsembleThreads()
      elseif kwargs[:parallel_type] == :split_threads
        ensemblealg = EnsembleSplitThreads()
      else
        @error "parallel_type value not recognized"
      end
    else
      ensemblealg = EnsembleThreads()
    end
    if :num_monte ∈ keys(kwargs)
      @warn "num_monte has been replaced by trajectories"
      trajectories = kwargs[:num_monte]
    else
      @assert :trajectories ∈ keys(kwargs)
      trajectories = kwargs[:trajectories]
    end
    if :parallel_type ∈ keys(kwargs)
      @warn "parallel_type has been deprecated. Please refer to the docs for the new dispatch-based system."
    end
    __solve(prob,alg,ensemblealg;trajectories=trajectories,kwargs...)
end

function __solve(prob::AbstractEnsembleProblem,
                 alg::Union{DEAlgorithm,Nothing},
                 ensemblealg::BasicEnsembleAlgorithm;
                 trajectories, batch_size = trajectories,
                 pmap_batch_size = batch_size÷100 > 0 ? batch_size÷100 : 1, kwargs...)

  num_batches = trajectories ÷ batch_size
  num_batches * batch_size != trajectories && (num_batches += 1)

  u = deepcopy(prob.u_init)
  converged = false
  elapsed_time = @elapsed for i in 1:num_batches
    if i == num_batches
      I = (batch_size*(i-1)+1):trajectories
    else
      I = (batch_size*(i-1)+1):batch_size*i
    end
    batch_data = solve_batch(prob,alg,ensemblealg,I,pmap_batch_size,kwargs...)
    u,converged = prob.reduction(u,batch_data,I)
    converged && break
  end
  if typeof(u) <: Vector{Any}
    _u = map(i->u[i],1:length(u))
  else
    _u = u
  end
  return EnsembleSolution(_u,elapsed_time,converged)
end

function batch_func(i,prob,alg,I,kwargs...)
  iter = 1
  new_prob = prob.prob_func(deepcopy(prob.prob),i,iter)
  rerun = true
  x = prob.output_func(solve(new_prob,alg;kwargs...),i)
  if !(typeof(x) <: Tuple)
      @warn("output_func should return (out,rerun). See docs for updated details")
      _x = (x,false)
  else
    _x = x
  end
  rerun = _x[2]
  while rerun
      iter += 1
      new_prob = prob.prob_func(deepcopy(prob.prob),i,iter)
      x = prob.output_func(solve(new_prob,alg;kwargs...),i)
      if !(typeof(x) <: Tuple)
          @warn("output_func should return (out,rerun). See docs for updated details")
          _x = (x,false)
      else
        _x = x
      end
      rerun = _x[2]
  end
  _x[1]
end

function solve_batch(prob,alg,::EnsembleDistributed,I,pmap_batch_size,kwargs...)
  wp=CachingPool(workers())
  batch_data = let
    pmap(wp,I,batch_size=pmap_batch_size) do i
      batch_func(i,prob,alg,I,kwargs...)
    end
  end
  map(i->batch_data[i],1:length(batch_data))
end

function solve_batch(prob,alg,::EnsembleSerial,I,pmap_batch_size,kwargs...)
  batch_data = let
    map(I) do i
      batch_func(i,prob,alg,I,kwargs...)
    end
  end
  map(i->batch_data[i],1:length(batch_data))
end

function solve_batch(prob,alg,::EnsembleThreads,I,pmap_batch_size,kwargs...)
  batch_data = Vector{Any}(undef,length(I))
  let
    Threads.@threads for batch_idx in axes(batch_data, 1)
        i = I[batch_idx]
        iter = 1
        new_prob = prob.prob_func(deepcopy(prob.prob),i,iter)
        x = prob.output_func(solve(new_prob,alg;kwargs...),i)
        if !(typeof(x) <: Tuple)
            warn("output_func should return (out,rerun). See docs for updated details")
            _x = (x,false)
        else
          _x = x
        end
        rerun = _x[2]

        while rerun
            iter += 1
            new_prob = prob.prob_func(deepcopy(prob.prob),i,iter)
            x = prob.output_func(solve(new_prob,alg;kwargs...),i)
            if !(typeof(x) <: Tuple)
                warn("output_func should return (out,rerun). See docs for updated details")
                _x = (x,false)
            else
              _x = x
            end
            rerun = _x[2]
        end
        batch_data[batch_idx] = _x[1]
    end
  end
  map(i->batch_data[i],1:length(batch_data))
end

function solve_batch(prob,alg,::EnsembleSplitThreads,I,pmap_batch_size,kwargs...)
  wp=CachingPool(workers())
  batch_data = let
    pmap(wp,1:nprocs(),batch_size=pmap_batch_size) do i
      thread_monte(prob,I,alg,i,kwargs...)
    end
  end
  _batch_data = vector_batch_data_to_arr(batch_data)
end

function thread_monte(prob,I,alg,procid,kwargs...)
  start = I[1]+(procid-1)*length(I)
  stop = I[1]+procid*length(I)-1
  portion = start:stop
  batch_data = Vector{Any}(undef,length(portion))
  let
    Threads.@threads for i in portion
      iter = 1
      new_prob = prob.prob_func(deepcopy(prob.prob),i,iter)
      rerun = true
      x = prob.output_func(solve(new_prob,alg;kwargs...),i)
      if !(typeof(x) <: Tuple)
          warn("output_func should return (out,rerun). See docs for updated details")
          _x = (x,false)
      else
        _x = x
      end
      rerun = _x[2]
      while rerun
          iter += 1
          new_prob = prob.prob_func(deepcopy(prob.prob),i,iter)
          x = prob.output_func(solve(new_prob,alg;kwargs...),i)
          if !(typeof(x) <: Tuple)
              warn("output_func should return (out,rerun). See docs for updated details")
              _x = (x,false)
          else
            _x = x
          end
          rerun = _x[2]
      end
      batch_data[i - start + 1] = _x[1]
    end
  end
  batch_data
end

function vector_batch_data_to_arr(batch_data)
  _batch_data = Vector{Any}(undef,sum((length(x) for x in batch_data)))
  idx = 0
  @inbounds for a in batch_data
    for x in a
      idx += 1
      _batch_data[idx] = x
    end
  end
  map(i->_batch_data[i],1:length(_batch_data))
end


text:
"""
    monte_carlo_return

Perform a monte_carlo_return using the provided environment, gvf, and start_state. For an environment to be compatible it must
implement `reset!(env, state)`

The interface is limited in that it can't do rollouts on compositional GVFs.
"""
function monte_carlo_return(env,
                            gvf,
                            start_state,
                            num_returns,
                            γ_thresh=1e-6,
                            max_steps=Int(1e6),
                            rng=Random.GLOBAL_RNG)

    returns = zeros(num_returns)

    for ret in 1:num_returns
        step = 0
        term = false
        cumulative_gamma = 1.0

        cur_state = start!(env, cur_state)
        next_action = StatsBase.sample(rng, policy(gvf), cur_state)

        while cumulative_gamma > γ_thresh &&
            step < max_steps &&
            term == false

            # Take action
            action = next_action
            next_state, r, term = step!(env, cur_state, action; rng=rng)

            # Get next action for GVFs
            next_action = StatsBase.sample(rng, policy(gvf), cur_state)

            # Update Return
            c, γ, pi_prob = get(gvf, cur_state, action, next_state, next_action, nothing)
            returns[ret] += cumulative_gamma*(1 - term)*c
            cumulative_gamma *= γ

            cur_state = next_state
            step += 1
        end
    end

    return returns
end

monte_carlo_returns(env, gvf, state_states, num_returns, γ_thresh, max_steps=Int(1e6), rng=Random.GLOBAL_RNG) =
    [monte_carlo_return(env, gvf, st, num_returns, γ_thres, max_steps=Int(1e6), rng=Random.GLOBAL_RNG) for st in start_states]


text:
#=
  Euler #21 in Julia.

  Problem 21
  """
  Let d(n) be defined as the sum of proper divisors of n (numbers less
  than n which divide evenly into n).
  If d(a) = b and d(b) = a, where a /= b, then a and b are an amicable
  pair and each of a and b are called amicable numbers.

  For example, the proper divisors of 220 are
  1, 2, 4, 5, 10, 11, 20, 22, 44, 55 and 110; therefore d(220) = 284.
  The proper divisors of 284 are 1, 2, 4, 71 and 142; so d(284) = 220.

  Evaluate the sum of all the amicable numbers under 10000.
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.12387843s
function euler21a()
    s = Dict()
    for a in 1:9999
        b = sum(all_divisors3(a));
        c = sum(all_divisors3(b));
        if a != b && a == c
            s[a] = 1;
            s[b] = 1;
        end
    end
    # return sum(keys(s).|>i->parse(Int,i))
    return keys(s)|>sum
end

#  Little more functional
# 0.12656805s
function euler21b()
    return sum(filter(v->v[1]!==v[2] && v[1] === v[3],(1:9999).|>a-> (
                  b = sum(all_divisors3(a));
                  c = sum(all_divisors3(b));
                  [a,b,c]
         )
         ).|>v->v[1])
end

# 0.07385759s
function euler21c()
    n = 9999
    s = zeros(Int,n)
    for i in 2:n
        s[i] = sum(all_divisors3(i))
    end
    a = []
    for i in 2:n
        # ignore perfect numbers...
        if s[i] <= n && i==s[s[i]] && i!= s[i]
            push!(a,i)
        end
    end
    return sum(a)
end

# run_euler(euler21a)
# run_euler(euler21b);
run_euler(euler21c);

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
const BUILTINS = Pair{String,Builtin}[
    "len"=>Builtin(
        function (args::Vararg{Object}; env::Environment = Environment())
            if length(args) != 1
                return ErrorObj(
                    "argument error: wrong number of arguments. Got $(length(args)) instead of 1",
                )
            end

            arg = args[1]
            if isa(arg, StringObj)
                return IntegerObj(length(arg.value))
            elseif isa(arg, ArrayObj)
                return IntegerObj(length(arg.elements))
            elseif isa(arg, HashObj)
                return IntegerObj(length(arg.pairs))
            end

            return ErrorObj(
                "argument error: argument to `len` is not supported, got $(type_of(arg))",
            )
        end,
    ),
    "first"=>Builtin(
        function (args::Vararg{Object}; env::Environment = Environment())
            if length(args) != 1
                return ErrorObj(
                    "argument error: wrong number of arguments. Got $(length(args)) instead of 1",
                )
            end

            arg = args[1]
            if isa(arg, StringObj)
                return length(arg.value) >= 1 ? StringObj(string(first(arg.value))) : _NULL
            elseif isa(arg, ArrayObj)
                return length(arg.elements) >= 1 ? first(arg.elements) : _NULL
            end

            return ErrorObj(
                "argument error: argument to `first` is not supported, got $(type_of(arg))",
            )
        end,
    ),
    "last"=>Builtin(
        function (args::Vararg{Object}; env::Environment = Environment())
            if length(args) != 1
                return ErrorObj(
                    "argument error: wrong number of arguments. Got $(length(args)) instead of 1",
                )
            end

            arg = args[1]
            if isa(arg, StringObj)
                return length(arg.value) >= 1 ? StringObj(string(last(arg.value))) : _NULL
            elseif isa(arg, ArrayObj)
                return length(arg.elements) >= 1 ? last(arg.elements) : _NULL
            end

            return ErrorObj(
                "argument error: argument to `last` is not supported, got $(type_of(arg))",
            )
        end,
    ),
    "rest"=>Builtin(
        function (args::Vararg{Object}; env::Environment = Environment())
            if length(args) != 1
                return ErrorObj(
                    "argument error: wrong number of arguments. Got $(length(args)) instead of 1",
                )
            end

            arg = args[1]
            if isa(arg, StringObj)
                if length(arg.value) >= 1
                    _, start = iterate(arg.value)
                    return StringObj(arg.value[start:end])
                else
                    return _NULL
                end
            elseif isa(arg, ArrayObj)
                l = length(arg.elements)
                return l >= 1 ? ArrayObj(arg.elements[2:end]) : _NULL
            end

            return ErrorObj(
                "argument error: argument to `rest` is not supported, got $(type_of(arg))",
            )
        end,
    ),
    "push"=>Builtin(
        function (args::Vararg{Object}; env::Environment = Environment())
            if length(args) != 2 && length(args) != 3
                return ErrorObj(
                    "argument error: wrong number of arguments. Got $(length(args)) instead of 2 or 3",
                )
            end

            if length(args) == 2
                if !isa(args[1], ArrayObj)
                    return ErrorObj(
                        "argument error: argument to `push` is not supported, got $(type_of(args[1]))",
                    )
                end

                arr = args[1]
                elements = copy(arr.elements)
                push!(elements, args[2])

                return ArrayObj(elements)
            else
                if !isa(args[1], HashObj)
                    return ErrorObj(
                        "argument error: argument to `push` is not supported, got $(type_of(args[1]))",
                    )
                end

                hash = args[1]
                pairs = copy(hash.pairs)
                push!(pairs, args[2] => args[3])

                return HashObj(pairs)
            end
        end,
    ),
    "puts"=>Builtin(function (args::Vararg{Object}; env::Environment = Environment())
        for arg in args
            if isa(arg, StringObj)
                println(env.output, arg.value)
            else
                println(env.output, arg)
            end
        end

        return _NULL
    end),
    "type"=>Builtin(
        function (args::Vararg{Object}; env::Environment = Environment())
            if length(args) != 1
                return ErrorObj(
                    "argument error: wrong number of arguments. Got $(length(args)) instead of 1",
                )
            end

            return StringObj(type_of(args[1]))
        end,
    ),
]

get_builtin_by_name(name::String) = begin
    for (_name, builtin) in BUILTINS
        if name == _name
            return builtin
        end
    end
end


text:
struct PdeBSM{B <: AbstractBlackScholesProcess}
  process::B
end

diffusion(pde::PdeBSM, t::Float64, x::Float64) = diffusion(pde.process, t, x)
drift(pde::PdeBSM, t::Float64, x::Float64) = drift(pde.process, t, x)
function discount(pde::PdeBSM, t::Float64, ::Float64)
  if abs(t) < 1e-8
    t = 0
  end

  return forward_rate(pde.process.riskFreeRate, t, t, ContinuousCompounding(), NoFrequency()).rate
end

mutable struct PdeConstantCoeff
  diffusion::Float64
  drift::Float64
  discount::Float64
end

function PdeConstantCoeff(process::AbstractBlackScholesProcess, t::Float64, x::Float64)
  pde = PdeBSM(process)
  diffusion_ = diffusion(pde, t, x)
  drift_ = drift(pde, t, x)
  discount_ = discount(pde, t, x)

  return PdeConstantCoeff(diffusion_, drift_, discount_)
end

# diffusion(pdecc::PdeConstantCoeff, ::Float64, ::Float64) = pdecc.diffusion
# drift(pdecc::PdeConstantCoeff, ::Float64, ::Float64) = pdecc.drift
# discount

function generate_operator!(pdecc::PdeConstantCoeff, t::Float64, tg::TransformedGrid, L::TridiagonalOperator)
  @inbounds @simd for i = 2:length(tg.grid) - 1
    sigma = pdecc.diffusion
    nu = pdecc.drift
    r = pdecc.discount

    sigma2 = sigma * sigma

    pd = -(sigma2 / tg.dxm[i] - nu) / tg.dx[i]
    pu = -(sigma2 / tg.dxp[i] + nu) / tg.dx[i]
    pm = sigma2 / (tg.dxm[i] * tg.dxp[i]) + r
    set_mid_row!(L, i, pd, pm, pu)
  end

  return L
end


text:
using BenchmarkTools

function func1(a::Array{Float64,4})

   #div = foo(a)
   @time div = zeros(size(a))

   @time div .= foo(a)

   a .= a + div

end

function func2(a::Array{Float64,4})

   let div = ones(10,10,10,8)
      #div = foo(a)
      a .= a + div
      #div[1] += 1
      #@show div[1]
   end

end

function foo(a::Array{Float64,4})
   div = a .+ 1
end

function foo(a::Array{Float64,1})
   div = a .+ 1
end

const func3 = let y = zeros(100)
    function inner_foo(x)
        for i in 1 : 100
            y[i] = 2*x[i] - 1
        end
        return sum(y)
    end
end

function test()
   a = zeros(10,10,10,8)
   #=
   @time for i=1:10
      func1(a)
   end
   =#
   #=
   @time for i=1:10
      func2(a)
   end
   =#
   b = zeros(100)
   @time for i=1:10
      func3(b)
   end
end


text:
function jaccard(s1, s2)::Float64
  return length(intersect(s1, s2)) / length(union(s1, s2))
end

function simpson(s1, s2)::Float64
  return length(intersect(s1, s2)) / min(length(s1), length(s2))
end

function dice(s1, s2)::Float64
  return (2 * length(intersect(s1, s2))) / (length(s1) + length(s2))
end

function okapi(h::Hypergraph, he::Int, node::Int, params=Dict("k1"=>2.0, "b"=>0.75))::Float64
  tf = 1.0 / length(getvertices(h, he))
  idf = log(((nhe(h)-length(gethyperedges(h, node)))+0.5) / (length(gethyperedges(h, node))+0.5))
  dl = length(getvertices(h, he))
  dl_ave = params["dl_ave"]
  k1 = haskey(params, "k1") ? params["k1"] : 2.0
  b = haskey(params, "b") ? params["b"] : 0.75
  v = (idf * tf * (k1+1)) / ((tf * k1) * (1 - b + b * (dl/dl_ave)))
  return v
end

function tf(h::Hypergraph, he::Int, node::Int, params=Dict())
  return 1.0 / length(getvertices(h, he))
end

function idf(h::Hypergraph, he::Int, node::Int, params=Dict())
  return log(nhe(h) / length(gethyperedges(h, node))) + 1
end

function tfidf(h::Hypergraph, he::Int, node::Int, params=Dict())::Float64
  return tf(h, he, node) * idf(h, he, node)
end

function transpose_tf(h::Hypergraph, he::Int, node::Int, params=Dict())
  # println(length(gethyperedges(h, node)))
  return 1.0 / length(gethyperedges(h, node))
end

function transpose_idf(h::Hypergraph, he::Int, node::Int, params=Dict())
  return log(nhv(h) / length(getvertices(h, he))) + 1
end

function transpose_tfidf(h::Hypergraph, he::Int, node::Int, params=Dict())
  return transpose_tf(h, he, node) * transpose_idf(h, he, node)
end

function tfidf_dfitf(h::Hypergraph, he::Int, node::Int, params=Dict())
  return tfidf(h, he, node) * transpose_tfidf(h, he, node)
end

function general_weight(h::Hypergraph, he::Int, node::Int, dl_ave)::Float64
  return 1.0 / length(getvertices(h, he))
end

function random_weight(h::Hypergraph=1, he=1, node=1, dl_ave=1)::Float64
  return rand(1)[1]
end


text:
#
# @description Defines the Cryptography module.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2015 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
module Cryptography

include("hashers.jl")
include("padders.jl")
include("ciphers.jl")

end


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

# jitter used to not work for categorical data

plot(x=rand(['a','b','c','d'],100), y=rand(100),
      Geom.point, Scale.x_discrete, Stat.x_jitter(range=0.2))


text:
include("../../AF.jl")

module MNIST

using AF

include("load.jl")

end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
using PyCall: @pycheckn, pyincref_, __pycall!

#########################################################################
struct PyWrapFn{N, RT}
    o::PyPtr
    pyargsptr::PyPtr
    ret::PyObject
end

function PyWrapFn(o::Union{PyObject, PyPtr}, nargs::Int, returntype::Type=PyObject)
    pyargsptr = ccall((@pysym :PyTuple_New), PyPtr, (Int,), nargs)
    ret = PyNULL()
    optr = o isa PyPtr ? o : o.o
    pyincref_(optr)
    return PyWrapFn{nargs, returntype}(optr, pyargsptr, ret)
end

(pf::PyWrapFn{N, RT})(args...) where {N, RT} =
    convert(RT, _pycall!(pf.ret, pf.pyargsptr, pf.o, args, N, C_NULL))

(pf::PyWrapFn{N, RT})() where {N, RT} =
    convert(RT, __pycall!(pf.ret, pf.pyargsptr, pf.o, C_NULL))

"""
```
pywrapfn(o::PyObject, nargs::Int, returntype::Type{T}=PyObject) where T
```
Wrap a callable PyObject/PyPtr possibly making calling it more performant. The
wrapped version (of type `PyWrapFn`) reduces the number of allocations made for
passing its arguments, and re-uses the same PyObject as its return value each
time it is called.

Mainly useful for functions called in a tight loop. After wrapping, arguments
should be passed in a tuple, rather than directly, e.g. `wrappedfn((a,b))` rather
than `wrappedfn(a,b)`.
Example
```
@pyimport numpy as np

# wrap a 2-arg version of np.random.rand for creating random matrices
randmatfn = pywrapfn(np.random["rand"], 2, PyArray)

# n.b. rand would normally take multiple arguments, like so:
a_random_matrix = np.random["rand"](7, 7)

# but we call the wrapped version with a tuple instead, i.e.
# rand22fn((7, 7)) not
# rand22fn(7, 7)
for i in 1:10^9
    arr = rand22fn((7,7))
    ...
end
```
"""
pywrapfn(o::PyObject, nargs::Int, returntype::Type=PyObject) =
    PyWrapFn(o, nargs, returntype)

"""
```
pysetargs!(w::PyWrapFn{N, RT}, args)
```
Set the arguments with which to call a Python function wrapped using
`w = pywrapfn(pyfun, ...)`
"""
function pysetargs!(pf::PyWrapFn{N, RT}, args) where {N, RT}
    check_pyargsptr(pf)
    pysetargs!(pf.pyargsptr, args, N)
end

"""
```
pysetarg!(w::PyWrapFn{N, RT}, arg, i::Integer=1)
```
Set the `i`th argument to be passed to a Python function previously
wrapped with a call to `w = pywrapfn(pyfun, ...)`
"""
function pysetarg!(pf::PyWrapFn{N, RT}, arg, i::Integer=1) where {N, RT}
    check_pyargsptr(pf)
    pysetarg!(pf.pyargsptr, arg, i)
end

"""
See check_pyargsptr(nargs::Int) above
"""
function check_pyargsptr(pf::PyWrapFn{N, RT}) where {N, RT}
    if unsafe_load(pf.pyargsptr).ob_refcnt > 1
        pydecref_(pf.pyargsptr)
        pf.pyargsptr =
            @pycheckn ccall((@pysym :PyTuple_New), PyPtr, (Int,), nargs)
    end
end


text:
# FIX THIS. NEITHER ELEGANT NOT PERFORMANT
export Material

mutable struct Material
    name::String
    color::RGBA{N0f8}
    mesh_size::Float64
end

function Material(;name::String = "", color = RGBA{N0f8}(0, 0, 1, 1), mesh_size = 1)
    return Material(name, parse(RGBA{N0f8}, color), Float64(mesh_size))
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
@testset "Test components" begin
    data = IS.SystemData()

    name = "component1"
    component = IS.TestComponent(name, 5)
    IS.add_component!(data, component)

    component = IS.get_component(IS.TestComponent, data, name)
    @test component.name == name

    components = IS.get_components(IS.TestComponent, data)
    @test length(components) == 1

    i = 0
    for component in IS.iterate_components(data)
        i += 1
    end
    @test i == 1

    io = IOBuffer()
    show(io, "text/plain", components)
    output = String(take!(io))
    expected = "TestComponent: $i"
    @test occursin(expected, output)

    IS.remove_component!(data, collect(components)[1])
    components = IS.get_components(IS.TestComponent, data)
    @test length(components) == 0

    IS.add_component!(data, component)
    components = IS.get_components_by_name(IS.InfrastructureSystemsComponent, data, name)
    @test length(components) == 1
    @test components[1].name == name

    IS.remove_components!(IS.TestComponent, data)
    components = IS.get_components(IS.TestComponent, data)
    @test length(components) == 0

    @test_throws ArgumentError IS.remove_components!(IS.TestComponent, data)
end

@testset "Test masked components" begin
    data = IS.SystemData()
    initial_time = Dates.DateTime("2020-09-01")
    resolution = Dates.Hour(1)
    ta = TimeSeries.TimeArray(range(initial_time; length = 24, step = resolution), ones(24))
    ts = IS.SingleTimeSeries(data = ta, name = "test")

    for i in 1:3
        name = "component_$(i)"
        component = IS.TestComponent(name, 5)
        IS.add_component!(data, component)
        IS.add_time_series!(data, component, ts)
    end

    component = IS.get_component(IS.TestComponent, data, "component_2")
    @test component isa IS.InfrastructureSystemsComponent
    IS.mask_component!(data, component)
    @test IS.get_component(IS.TestComponent, data, "component_2") === nothing
    @test IS.get_masked_component(IS.TestComponent, data, "component_2") isa
          IS.TestComponent
    @test collect(IS.get_masked_components(IS.TestComponent, data)) == [component]
    @test IS.get_masked_components_by_name(
        IS.InfrastructureSystemsComponent,
        data,
        "component_2",
    ) == [component]
    @test IS.get_time_series(IS.SingleTimeSeries, component, "test") isa IS.SingleTimeSeries
    @test IS.is_attached(component, data.masked_components)

    # This needs to return time series for masked components.
    @test length(collect(IS.get_time_series_multiple(data, type = IS.SingleTimeSeries))) ==
          3
end


text:
using Test
using DifferenceEquations
using Distributions
using LinearAlgebra
using Random

include("matrix_vector_of_vectors.jl")
include("linear_likelihood.jl")
include("linear_simulations.jl")
include("quadratic_likelihood.jl")
include("quadratic_simulations.jl")
include("sciml_interfaces.jl")


text:
using Meshes

# Produce a level set function that is a noisy version of the distance from
# the origin (such that level sets are noisy spheres).
#
# The noise should exercise marching tetrahedra's ability to produce a water-
# tight surface in all cases (unlike standard marching cubes).
#
N = 10
sigma = 1.0
distance = Float32[ sqrt(float32(i*i+j*j+k*k)) for i = -N:N, j = -N:N, k = -N:N ]
distance = distance + sigma*rand(2*N+1,2*N+1,2*N+1)

# Extract an isosurface.
#
lambda = N-2*sigma # isovalue

tic()
msh = isosurface(distance,lambda)
toc()

# Simplify the mesh
#
#msh2 = simplify(msh,0.1)

# Export the mesh to a ply.
#
# The mesh can be visualized, e.g., in MeshLab (http://meshlab.sourceforge.net/).
#
exportToPly(msh,"noisy_sphere.ply")

# test contatenation
msh3 = merge(msh,msh)

function testUnionNotInterection()
	# creates Union, Not, Intersection of cube and sphere
	#
	# http://en.wikipedia.org/wiki/Constructive_solid_geometry

	tic()
	# volume of interest
	x_min, x_max = -1, 15
	y_min, y_max = -1, 5
	z_min, z_max = -1, 5
	scale = 8

	b1(x,y,z) = box(   x,y,z, 0,0,0,3,3,3)
	s1(x,y,z) = sphere(x,y,z, 3,3,3,sqrt(3))
	f1(x,y,z) = min(b1(x,y,z), s1(x,y,z))  # UNION
	b2(x,y,z) = box(   x,y,z, 5,0,0,8,3,3)
	s2(x,y,z) = sphere(x,y,z, 8,3,3,sqrt(3))
	f2(x,y,z) = max(b2(x,y,z), -s2(x,y,z)) # NOT
	b3(x,y,z) = box(   x,y,z, 10,0,0,13,3,3)
	s3(x,y,z) = sphere(x,y,z, 13,3,3,sqrt(3))
	f3(x,y,z) = max(b3(x,y,z), s3(x,y,z))  # INTERSECTION
	f(x,y,z) = min(f1(x,y,z), f2(x,y,z), f3(x,y,z))

	vol = volume(f, x_min,y_min,z_min,x_max,y_max,z_max, scale)
	msh = isosurface(vol, 0.0)
	exportToStl(msh, "wiki_csg.stl")
	toc()
end
testUnionNotInterection()

function testCylinders()
	tic()
	c1(x,y,z) = cylinderX(x,y,z, 0,0,1,-2,4)
	c2(x,y,z) = cylinderY(x,y,z, 0,0,1,-2,4)
	c3(x,y,z) = cylinderZ(x,y,z, 0,0,1,-2,4)
	f(x,y,z) = min(c1(x,y,z), c2(x,y,z), c3(x,y,z))

	vol = volume(f, -3,-3,-3, 5,5,5, 2)
	msh = isosurface(vol, 0.0)
	exportToStl(msh, "cylinders_test.stl")
	toc()
end
testCylinders()


text:
using Pkg
using SafeTestsets

# Start Test Script

@time begin
@time @safetestset "Tests básicos" begin include("cinta_test.jl") end
end  # @time

max_stars_repo_name: AndresParraSilva/julia-start

max_issues_repo_name: AndresParraSilva/julia-start

max_forks_repo_name: AndresParraSilva/julia-start


text:
"""
    maximum_weight_maximal_matching{T<:Real}(g, w::Dict{Edge,T})
    maximum_weight_maximal_matching{T<:Real}(g, w::Dict{Edge,T}, cutoff)

Given a bipartite graph `g` and an edgemap `w` containing weights associated to edges,
returns a matching with the maximum total weight among the ones containing the
greatest number of edges.

Edges in `g` not present in `w` will not be considered for the matching.

The algorithm relies on a linear relaxation on of the matching problem, which is
guaranteed to have integer solution on bipartite graps.

Eventually a `cutoff` argument can be given, to reduce computational times
excluding edges with weights lower than the cutoff.

The package JuMP.jl and one of its supported solvers is required.

The returned object is of type `MatchingResult`.
"""
function maximum_weight_maximal_matching end

function maximum_weight_maximal_matching(g::Graph, solver::AbstractMathProgSolver, w::AbstractMatrix{T}, cutoff::R) where {T<:Real, R<:Real}
    return maximum_weight_maximal_matching(g, solver, cutoff_weights(w, cutoff))
end

function maximum_weight_maximal_matching(g::Graph, solver::AbstractMathProgSolver, w::AbstractMatrix{T}) where {T<:Real}
# TODO support for graphs with zero degree nodes
# TODO apply separately on each connected component
    bpmap = bipartite_map(g)
    length(bpmap) != nv(g) && error("Graph is not bipartite")
    v1 = findin(bpmap, 1)
    v2 = findin(bpmap, 2)
    if length(v1) > length(v2)
        v1, v2 = v2, v1
    end

    nedg = 0
    edgemap = Dict{Edge,Int}()

    for j in 1:size(w,2)
        for i in 1:size(w,1)
            if w[i,j] > 0.0
                nedg += 1
                edgemap[Edge(i,j)] = nedg
                edgemap[Edge(j,i)] = nedg
            end
        end
    end

    model = Model(solver=solver)
    @variable(model, x[1:length(w)] >= 0)

    for i in v1
        idx = Vector{Int}()
        for j in neighbors(g, i)
            if haskey(edgemap, Edge(i,j))
                push!(idx, edgemap[Edge(i,j)])
            end
        end
        if length(idx) > 0
            @constraint(model, sum(x[id] for id = idx) == 1)
        end
    end

    for j in v2
        idx = Vector{Int}()
        for i in neighbors(g, j)
            if haskey(edgemap, Edge(i,j))
                push!(idx, edgemap[Edge(i,j)])
            end
        end

        if length(idx) > 0
            @constraint(model, sum(x[id] for id = idx) <= 1)
        end
    end

    @objective(model, Max, sum(w[src(e),dst(e)] * x[edgemap[e]] for e in keys(edgemap)))

    status = solve(model)
    status != :Optimal && error("JuMP solver failed to find optimal solution.")
    sol = getvalue(x)

    all(Bool[s == 1 || s == 0 for s in sol]) || error("Found non-integer solution.")

    cost = getobjectivevalue(model)

    mate = fill(-1, nv(g))
    for e in edges(g)
        if w[src(e),dst(e)] > zero(T)
            inmatch = convert(Bool, sol[edgemap[e]])
            if inmatch
                mate[src(e)] = dst(e)
                mate[dst(e)] = src(e)
            end
        end
    end

    return MatchingResult(cost, mate)
end

"""
    cutoff_weights copies the weight matrix with all elements below cutoff set to 0
"""
function cutoff_weights(w::AbstractMatrix{T}, cutoff::R) where {T<:Real, R<:Real}
    wnew = copy(w)
    for j in 1:size(w,2)
        for i in 1:size(w,1)
            if wnew[i,j] < cutoff
                wnew[i,j] = zero(T)
            end
        end
    end
    wnew
end


text:
const USE_GPU = false  # Use GPU? If this is set false, then no GPU needs to be available
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
    macro pow(args...)  esc(:(CUDA.pow($(args...)))) end
else
    @init_parallel_stencil(Threads, Float64, 2)
    pow(x,y) = x^y
    macro pow(args...)  esc(:(pow($(args...)))) end
end
using Plots, Printf, Statistics

@parallel function compute_nonlin!(K_mu::Data.Array, EtaR::Data.Array, Phi::Data.Array, Eta::Data.Array, k_μ0::Data.Number, n::Data.Number)
    @all(K_mu) = k_μ0*@pow(@all(Phi), n)
    @all(EtaR) = @all(Eta)
    return
end

@parallel_indices (ix,iy) function compute_EtaR!(EtaR::Data.Array, Pe::Data.Array, ηR::Data.Number)
    if (ix<size(EtaR,1) && iy<size(EtaR,2)) if (Pe[ix+1,iy+1]<0.0) EtaR[ix,iy] = ηR; end;  end
    return
end

@parallel function compute_flux!(qx::Data.Array, qy::Data.Array, K_mu::Data.Array, Pe::Data.Array, Δρg::Data.Number, dx::Data.Number, dy::Data.Number)
    @all(qx) = -@av_xi(K_mu)* @d_xi(Pe)/dx
    @all(qy) = -@av_yi(K_mu)*(@d_yi(Pe)/dy + Δρg)
    return
end

@parallel function compute_Pe_Phi!(dPedt::Data.Array, Pe::Data.Array, Phi::Data.Array, dPhidt::Data.Array, qx::Data.Array, qy::Data.Array, EtaR::Data.Array, β::Data.Number, dt::Data.Number, dx::Data.Number, dy::Data.Number)
    @all(dPedt)  = 1.0/β*(-@d_xa(qx)/dx -@d_ya(qy)/dy -1.0/@all(EtaR)*@inn(Pe))
    @inn(Pe)     = @inn(Pe) + dt*@all(dPedt)
    @all(dPhidt) = -β*@all(dPedt) - 1.0/@all(EtaR)*@inn(Pe)
    @inn(Phi)    = max(0.0, @inn(Phi) + dt*@all(dPhidt))
    return
end

##################################################
@views function poroVEP2D()
    chan      = 0                    # chan=1: channel setup; chan=0: blob setup
    # Physics
    sc        = 1.0+chan*30          # domain extend scaling
    lx, ly    = 250.0/sc, 500.0/sc   # domain extend
    k_μ0      = 1.0                  # permeability/fuild viscosity
    Δρg       = 1.0                  # density*gravity (total-fluid)
    n         = 3.0                  # Karman-Cozeny exponant
    η         = 1.0                  # compaction bulk viscosity
    ηR        = η/(1.0+chan*1000)    # decompaction bulk viscosity
    β         = 0.2                  # bulk compressibility
    λ         = 20.0/sc              # wave length
    ϕ0        = 1.0                  # background porosity
    ϕA        = 2.0                  # perturbation amplitude
    t         = 0.0                  # physical time initialisation
    # Numerics
    nx, ny    = 256, 512             # numerical grid resolution; should be a mulitple of 32 for optimal GPU perf
    nt        = 4e4                  # if chan=0 setup - number of timesteps
    nout      = 5e2                  # if chan=0 setup - plotting frequency
    # nt        = 2e6                  # if chan=1 setup - number of timesteps
    # nout      = 2e4                  # if chan=1 setup - plotting frequency
    # Derived numerics
    dx, dy    = lx/(nx-1), ly/(ny-1) # cell sizes
    # Array allocations
    Eta       = η*@ones(nx-2,ny-2)
    EtaR      =  @zeros(nx-2,ny-2)
    Pe        =  @zeros(nx  ,ny  )
    dPedt     =  @zeros(nx-2,ny-2)
    dPhidt    =  @zeros(nx-2,ny-2)
    K_mu      =  @zeros(nx  ,ny  )
    qx        =  @zeros(nx-1,ny-2)
    qy        =  @zeros(nx-2,ny-1)
    # Initial conditions
    Phi       = ϕ0*ones(nx  ,ny  )
    Radc      = zeros(nx  ,ny  )
    Radc     .= [ (((ix-1)*dx-0.5*lx)/λ/4.)^2 + (((iy-1)*dy-0.2*ly)/λ)^2 for ix=1:size(Radc,1), iy=1:size(Radc,2)]
    Phi[Radc.<1.0] .= Phi[Radc.<1.0] .+ ϕA
    Phi[2:end-1,2:end-1] .= Phi[2:end-1,2:end-1] .+ 1.0/4.1*(diff(diff(Phi[:,2:end-1],dims=1),dims=1).+diff(diff(Phi[2:end-1,:],dims=2),dims=2))
    Phi       = Data.Array(Phi)
    dt        = 1.0
    # Preparation of visualisation
    ENV["GKSwstype"]="nul"; if isdir("viz2D_out")==false mkdir("viz2D_out") end; loadpath = "./viz2D_out/"; anim = Animation(loadpath,String[])
    println("Animation directory: $(anim.dir)")
    X, Y      = -lx/2:dx:lx/2, -ly/2:dy:ly/2
    # Time loop
    for it = 1:nt
        if (it==11)  global wtime0 = Base.time()  end
        @parallel compute_nonlin!(K_mu, EtaR, Phi, Eta, k_μ0, n)
        @parallel compute_EtaR!(EtaR, Pe, ηR)
        if mod(it,100)==1  dt = dx^2/(maximum(K_mu)/β)/5.0  end
        @parallel compute_flux!(qx, qy, K_mu, Pe, Δρg, dx, dy)
        @parallel compute_Pe_Phi!(dPedt, Pe, Phi, dPhidt, qx, qy, EtaR, β, dt, dx, dy)
        t = t + dt
        # Visualisation
        if mod(it,nout)==0
            @printf("it=%d, time=%1.3e \n", it, t)
            p1 = heatmap(X, Y, Array(Phi)', aspect_ratio=1, xlims=(X[1],X[end]), ylims=(Y[1],Y[end]), c=:inferno, title="Porosity")
            p2 = heatmap(X, Y, Array(Pe)',  aspect_ratio=1, xlims=(X[1],X[end]), ylims=(Y[1],Y[end]), c=:inferno, title="Effective pressure")
            # display(plot(p1, p2))
            plot(p1, p2); frame(anim)
        end
    end
    # Performance
    wtime    = Base.time() - wtime0
    A_eff    = (2*2+2)/1e9*nx*ny*sizeof(Data.Number)  # Effective main memory access per iteration [GB] (Lower bound of required memory access: H and dHdτ have to be read and written (dHdτ for damping): 4 whole-array memaccess; B has to be read: 1 whole-array memaccess)
    wtime_it = wtime/(nt-10)                          # Execution time per iteration [s]
    T_eff    = A_eff/wtime_it                         # Effective memory throughput [GB/s]
    @printf("Total steps=%d, time=%1.3e sec (@ T_eff = %1.2f GB/s) \n", nt, wtime, round(T_eff, sigdigits=2))
    return
    gif(anim, "poroVEP2D.gif", fps = 15)
end

poroVEP2D()


text:
function DBMIMMPSolver(A,b,C,d,VariableDef;XLB=[],XUB=[],MU = .9,EnhancedUB = 0,EnhancedCut = 0,SelectVariableMode = 1,SelectNodeMode = 0, TimeLim = 2e+10,
                rEPS = .000001, aEPS = .000001)
    ### Initialization
    t0 = time_ns();
    nVar = size(A)[2];
    if XLB == []
        XLB = -1e9*ones(nVar)
    end
    if XUB== []
        XUB = 1e9*ones(nVar)
    end

    # Algorithm's setting ##################################################

    # Initializing the Type 2 node selection
    Type2Iter = 0;
    Type2LastLB = 0;
    Type2Flag = 0;
    Type2MaxIter = 1000;
    Type2MaxImprovement = 0.05;

    # Initializing the Type 3 or 4 variable selection
    if (SelectVariableMode == 3) || (SelectVariableMode == 4)
        TempSelectVariableMode = SelectVariableMode
        SelectVariableMode = 1
    end
    # Initializing the Type 3 or 4 node selection
    if (SelectNodeMode == 3) || (SelectNodeMode == 4)
        TempSelectNodeMode = SelectNodeMode
        SelectNodeMode = 0
    end
    ########################################################################

    const  inf = 9999999;
    const  EPS = 0.000001
    nObj = size(C)[1];
    nVar = size(A)[2];
    nCons = length(b);

    #  Initial Global Bounds
    GLB = 0;
    GUB = inf;

    # Initialization of Temp Variables
    NodeCounter = 1;
    RO1 = zeros(sum(VariableDef[2:3]));
    n1 = zeros(sum(VariableDef[2:3]));
    RO2 = zeros(sum(VariableDef[2:3]));
    n2 = zeros(sum(VariableDef[2:3]));

    ## Main Algorithm
    Tree = []
    ## Preprocessing
    Cuts = [];
    if EnhancedUB == 1 | EnhancedCut == 1
        nVars = size(C,2)
        Xs, Ys, GUB, Cuts = Preprocess2(A,b,C,d,XLB,XUB,VariableDef,TimeLim = .01*nVars)
        if Xs != []
            OptX = Xs;
        end
        if EnhancedUB == 0
            UB = inf;
        end
        if EnhancedCut == 0
            Cut = []
        end
    else
        OptX = [];
    end
    # ---------------------------------------------
    Time = (time_ns() - t0)/1e9
    if TimeLim-Time < 0.001;
        println("Time Limit reached! The solution is not a proved optimal")
        return GUB, OptX
    end
    # ---------------------------------------------

    GLB, X_Tilde,NewUB,NewX = ShaoRSolver(A,b,C,d,XLB,XUB,VariableDef,CUTs = Cuts,TimeLim = TimeLim - Time);

    if GLB == -2
        Time = (time_ns() - t0)/1e9
        println("Something went wrong, the problem seems to be infeasible or has numerical issues!")
        return -2, []
    end
    iter = 1;
    if isempty(X_Tilde)
        if GUB > NewUB
            GUB = NewUB
            OptX = NewX
            GLB = GUB
            return GUB, OptX
        else
            println("The problem is infeasible!")
            return GUB, []
        end
    elseif VariableDef[1] == nVar
        OptX = X_Tilde;
        Time = (time_ns() - t0)/1e9
    elseif all(isinteger.(X_Tilde[Int(VariableDef[1]+1):end]))
        # We get the optimal solution
        println("The relaxed solution was optimal")
        OptX = X_Tilde;
        GUB = GLB
        Time = (time_ns() - t0)/1e9
    else
        delta1 = (ceil.(X_Tilde[VariableDef[1]+1:end])-X_Tilde[VariableDef[1]+1:end]).*RO1./(EPS+n1)
        delta2 = (X_Tilde[VariableDef[1]+1:end]-floor.(X_Tilde[VariableDef[1]+1:end])).*RO2./(EPS+n2)
        SelectedVariable = SelectVariable(VariableDef,X_Tilde,delta1,delta2,MU,
        n1 = n1, n2 = n2, TAV =100, ZZ = GLB ,xlb = XLB,xub = XUB, Type = SelectVariableMode)

        # Node 1
        NewXUB = copy(XUB);
        NewXUB[SelectedVariable] = floor(X_Tilde[SelectedVariable])
        # Time commands
        # ---------------------------------------------
        Time = (time_ns() - t0)/1e9
        if TimeLim-Time < 0.001;
            println("Time Limit reached! The solution is not a proved optimal")
            return GUB, OptX
        end
        # ---------------------------------------------
        Z, X_Tilde1,NewUB = ShaoRSolver(A,b,C,d,XLB,NewXUB,VariableDef,CUTs = Cuts,TimeLim = TimeLim - Time);
        if GUB > NewUB
            GUB = NewUB
        end
        # Time commands
        # ---------------------------------------------
        Time = (time_ns() - t0)/1e9
        if TimeLim-Time < 0.001;
            println("Time Limit reached! The solution is not a proved optimal")
            return GUB, OptX
        end
        # ---------------------------------------------
        if Z == -2
            println("Something went wrong, the problem seems to be infeasible or has numerical issues!")
            return -2, []
        end
        if !isempty(Z) && Z < GUB
            if all(isinteger.(X_Tilde1[VariableDef[1]+1:end]))
                GUB = Z;
                OptX = copy(X_Tilde1);
            else
                Tree = [Node(copy(X_Tilde1),copy(XLB),NewXUB,Z,GUB)];
                RO2[SelectedVariable-VariableDef[1]] = RO2[SelectedVariable-VariableDef[1]] + (Z-GLB)/(X_Tilde[SelectedVariable] - NewXUB[SelectedVariable])
                n2[SelectedVariable-VariableDef[1]] = n2[SelectedVariable-VariableDef[1]] + 1
            end
        end
        NodeCounter = NodeCounter + 1;

        # Node 2

        NewXLB = copy(XLB);
        NewXLB[SelectedVariable] = ceil(X_Tilde[SelectedVariable])
        # Time commands
        # ---------------------------------------------
        Time = (time_ns() - t0)/1e9
        if TimeLim-Time < 0.001;
            println("Time Limit reached! The solution is not a proved optimal")
            return GUB, OptX
        end
        # ---------------------------------------------
        Z, X_Tilde2,NewUB = ShaoRSolver(A,b,C,d,NewXLB,XUB,VariableDef,CUTs = Cuts,TimeLim = TimeLim - Time);
        if GUB > NewUB
            GUB = NewUB
        end
        # Time commands
        # ---------------------------------------------
        Time = (time_ns() - t0)/1e9
        if TimeLim-Time < 0.001;
            println("Time Limit reached! The solution is not a proved optimal")
            return GUB, OptX
        end
        # ---------------------------------------------
        if Z == -2
            println("Something went wrong, the problem seems to be infeasible or has numerical issues!")
            return -2, []
        end
        if !isempty(Z) && Z < GUB
            if all(isinteger.(X_Tilde2[VariableDef[1]+1:end]))
                GUB = Z;
                OptX = copy(X_Tilde2);
            else
                push!(Tree,Node(copy(X_Tilde2),NewXLB,copy(XUB),Z,GUB))
                RO1[SelectedVariable-VariableDef[1]] = RO1[SelectedVariable-VariableDef[1]] + (Z-GLB)/(NewXLB[SelectedVariable] - X_Tilde[SelectedVariable] )
                n1[SelectedVariable-VariableDef[1]] = n1[SelectedVariable-VariableDef[1]] + 1
            end
        end
        NodeCounter = NodeCounter + 1;
        if !isempty(Tree)
            GLB = GLBUpdater(Tree)
        end
        # Main Loop
        while (length(Tree) > 0) && (abs(GUB-GLB) >= aEPS) && ((GUB-GLB)/GUB >= rEPS)
            # Type 2 Node Selection modification
            # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            Type2Iter = Type2Iter + 1;
            if Type2Iter > Type2MaxIter
                Type2Iter = 0;
                Type2LastLB = GLB;
                Type2Flag = 1 - Type2Flag;
            elseif Type2Flag == 1 && (abs(Type2LastLB-GLB)/Type2LastLB >Type2MaxImprovement)
                Type2Flag = 0;
                Type2Iter = 0;
            end
            # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

            # Type 3 or 4 Variable Selection modification
            # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            if (NodeCounter > 10) && ((SelectVariableMode == 3) || (SelectVariableMode == 4))
                SelectVariableMode = TempSelectVariableMode;
            end
            # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

            # Type 3 or 4 Variable Selection modification
            # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            if (NodeCounter > 10) && ((SelectNodeMode == 3) || (SelectNodeMode == 4))
                SelectNodeMode = TempSelectNodeMode
            end

            # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

            # Select node to branch
            SelectedNode = SelectNode(Tree,VariableDef = VariableDef, MU = MU, RO1 = RO1,
                                     RO2 = RO2, Type = SelectNodeMode, Type2Flag = Type2Flag);
            node = deepcopy(Tree[SelectedNode])
            deleteat!(Tree,SelectedNode);

            # Select variable to branch
            delta1 = (ceil.(node.X[VariableDef[1]+1:end])-node.X[VariableDef[1]+1:end]).*RO1./(EPS+n1)
            delta2 = (node.X[VariableDef[1]+1:end]-floor.(node.X[VariableDef[1]+1:end])).*RO2./(EPS+n2)
            SelectedVariable = SelectVariable(VariableDef,node.X,delta1,delta2,MU,n1 = n1, n2 = n2,
                                    TAV =100,ZZ = node.LB, xlb = node.XLB, xub = node.XUB, Type = SelectVariableMode)

            # Node 1,TimeLim = TimeLim - Time
            NewXUB = copy(node.XUB);
            NewXUB[SelectedVariable] = floor(node.X[SelectedVariable])
            # Time commands
            # ---------------------------------------------
            Time = (time_ns() - t0)/1e9
            if TimeLim-Time < 0.001;
                println("Time Limit reached! The solution is not a proved optimal")
                return GUB, OptX
            end
            # ---------------------------------------------
            Z, X_Tilde,NewUB = ShaoRSolver(A,b,C,d,copy(node.XLB),NewXUB,VariableDef,CUTs = Cuts,TimeLim = TimeLim - Time);
            if GUB > NewUB
                GUB = NewUB
            end
            # Time commands
            # ---------------------------------------------
            Time = (time_ns() - t0)/1e9
            if TimeLim-Time < 0.001;
                println("Time Limit reached! The solution is not a proved optimal")
                return GUB, OptX
            end
            # ---------------------------------------------
            if Z == -2
                println("Something went wrong, the problem seems to be infeasible or has numerical issues!")
                return -2, []
            end
            if !isempty(Z) && Z < GUB
                if all(isinteger.(X_Tilde[VariableDef[1]+1:end]))
                    GUB = Z;
                    OptX = copy(X_Tilde);
                    # Type 2 Node Selection modification
                    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
                    if Type2Flag == 0
                        Type2Flag = 1;
                        Type2Iter = 0;
                        Type2LastLB = GLB;
                    end
                    # Type 2 Node Selection modification
                    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
                else
                    push!(Tree,Node(copy(X_Tilde),copy(node.XLB),NewXUB,Z,GUB));
                    RO2[SelectedVariable-VariableDef[1]] = RO2[SelectedVariable-VariableDef[1]] +
                        (Z-node.LB)/(node.X[SelectedVariable] - NewXUB[SelectedVariable])
                    n2[SelectedVariable-VariableDef[1]] = n2[SelectedVariable-VariableDef[1]] + 1
                end
            end
            NodeCounter = NodeCounter + 1;
            # Node 2
            NewXLB = copy(node.XLB);
            NewXLB[SelectedVariable] = ceil(node.X[SelectedVariable])
            # Time commands
            # ---------------------------------------------
            Time = (time_ns() - t0)/1e9
            if TimeLim-Time < 0.001;
                println("Time Limit reached! The solution is not a proved optimal")
                return GUB, OptX
            end
            # ---------------------------------------------
            Z, X_Tilde,NewUB = ShaoRSolver(A,b,C,d,NewXLB,copy(node.XUB),VariableDef,CUTs = Cuts,TimeLim = TimeLim - Time);
            if GUB > NewUB
                GUB = NewUB
            end
            # Time commands
            # ---------------------------------------------
            Time = (time_ns() - t0)/1e9
            if TimeLim-Time < 0.001;
                println("Time Limit reached! The solution is not a proved optimal")
                return GUB, OptX
            end
            # ---------------------------------------------
            if Z == -2
                println("Something went wrong, the problem seems to be infeasible or has numerical issues!")
                return -2, []
            end
            if !isempty(Z) && Z < GUB
                if all(isinteger.(X_Tilde[VariableDef[1]+1:end]))
                    GUB = Z;
                    OptX = copy(X_Tilde);
                    # Type 2 Node Selection modification
                    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
                    if Type2Flag == 0
                        Type2Flag = 1;
                        Type2Iter = 0;
                        Type2LastLB = GLB;
                    end
                    # Type 2 Node Selection modification
                    # +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
                else
                    push!(Tree,Node(copy(X_Tilde),NewXLB,copy(node.XUB),Z,GUB))
                    RO1[SelectedVariable-VariableDef[1]] = RO1[SelectedVariable-VariableDef[1]] +
                        (Z-node.LB)/(NewXLB[SelectedVariable] - node.X[SelectedVariable] )
                    n1[SelectedVariable-VariableDef[1]] = n1[SelectedVariable-VariableDef[1]] + 1
                end
            end
            NodeCounter = NodeCounter + 1;

            if !isempty(Tree)
                GLB = GLBUpdater(Tree)
            end
            # println("Iteration = $iter  |Time = $Time| Node = $NodeCounter  |  GLB = $GLB  |   GUB =  $GUB" )
            iter = iter + 1;

        end

    end

    Time = (time_ns() - t0)/1e9
    return GUB, OptX
end


text:
writing_output = false
if VERSION < v"1.5"
    ver = "111"
else
    ver = "150"
end

@everywhere Random.seed!(42)

weights = rand(400)
weights = weights ./ sum(weights)

test_sys_resample    = SMC.resample(weights, method = :systematic)
test_multi_resample  = SMC.resample(weights, method = :multinomial)
test_poly_resample   = SMC.resample(weights, method = :polyalgo)

if writing_output
    jldopen("reference/resample_version=" * ver * ".jld2",
            true, true, true, IOStream) do file
        write(file, "sys", test_sys_resample)
        write(file, "multi", test_multi_resample)
        write(file, "poly", test_poly_resample)
    end
end

saved_sys_resample   = load("reference/resample_version=" * ver * ".jld2", "sys")
saved_multi_resample = load("reference/resample_version=" * ver * ".jld2", "multi")
saved_poly_resample  = load("reference/resample_version=" * ver * ".jld2", "poly")

####################################################################

@testset "Resampling methods" begin
    @test test_sys_resample   == saved_sys_resample
    @test test_multi_resample == saved_multi_resample
    @test test_poly_resample  == saved_poly_resample
end


text:
using DrWatson
@quickactivate "AccuracyAtTop_aaai"

include(srcdir("utilities.jl"))
include(srcdir("tfco.jl"))

# ------------------------------------------------------------------------------------------
# Settings
# ------------------------------------------------------------------------------------------
Dataset_Settings = Dict(
    :dataset => [FashionMNIST, CIFAR100, Molecules],
    :posclass => [0, 1, 1],
)

Train_Settings = Dict(
    :batchsize => 32,
    :epochs => 200,
    :optimiser => Descent,
    :steplength => 0.0001,
    :seed => collect(1:10),
)

Model_Settings = Dict(
    :type => [TFCO],
    :arg => [0.01],
    :surrogate => missing,
    :reg => 1e-3,
    :buffer => missing,
)

run_simulations_tfco(Dataset_Settings, Train_Settings, Model_Settings)
run_evaluation(Dataset_Settings)


text:
# # [Ross-Macdonald Malaria Modeling](@id ross-macdonald)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/Ross-Macdonald.ipynb)

# Authors: Sean Wu and Sophie Libkind

using AlgebraicDynamics.DWDDynam
using Catlab.WiringDiagrams
using LabelledArrays

using Catlab.Graphics
using DelayDiffEq, DifferentialEquations
using Plots

# The Ross-Macdonald (RM) model is the canonical model of malaria transmission, first expressed by
# Sir Ronald Ross in 1911. The simplest form of the model is given by a system of 2 nonlinear ordinary
# differential equations describing $Z$, parasite prevalence in the mosquito vector population, and $X$,
# parasite prevalence in the vertebrate host population:
# $$ \dot X = mazbZ(1 - X) - rX, \quad \dot Z = acX (e^{-gn} - Z) - gZ  \quad (1)$$

#-

# The RM model is of supreme importance to modeling of malaria and other mosquito transmitted diseases
# because it is the simplest model that retains the fundamental processes involved in transmission: recovery
# of infected hosts, death of infected mosquitoes, and parasite transmission between discordant pairs of
# vectors and hosts. In this sense it represents the irreducible complexity to modeling mosquito borne
# diseases, such that any simpler model would be insufficient.

# The model assumes that mosquitoes take bloodmeals (bites) on vertebrate hosts with rate $a$.
# Therefore, the per-capita rate at which susceptible mosquitoes become infected is $acX$, where $c$
# is the transmission efficiency from infectious humans to susceptible mosquitoes (probability of
# parasites successfully invading the mosquito given a bite occurs) and $X$ is the probability the
#  bite lands on an infectious human. $1-Z$ is the proportion of mosquitoes who are susceptible but
# given an incubation period of $n$ days during which mosquitoes suffer per-capita mortality at rate
# $g$, only a proportion $e^{-gn}$ survive incubation to become infectious. Therefore the classic
# RM model will have statics that are consistent with the incorporation of an incubation period of
# fixed duration, but dynamics will not, as the delay is not incorporated into the equations.

#-

# Susceptible humans acquire infection at a rate $mabZ$, where $m = M/H$, the ratio of mosquitoes to
# humans, $a$ is as before, and $b$ is the transmission efficiency from infectious mosquitoes to
# susceptible humans. Therefore $mabZ$ gives the per-capita rate at which susceptible humans recieve
# parasites from the mosquito population. Infectious humans recover at a rate $r$ (the rate at which the
# immune system clears parasites). Because $m$ is considered a constant parameter, the model ignores
# seasonal fluctuations in mosquito (or host) population sizes. However, because the timescale of host
# dynamcis is typically slow in relation to mosquito dynamics, the approximation is often suitable
# over a transmission season.

#-

# The model's threshold criterion between the trivial equilibrium where both $Z,X$ are equal to zero
# (no disease) and the endemic equilibrium can be expressed by the basic reproductive number, giving
# the expected number of secondary infections arising from the introduction of a single infectious human,
# $$R_0 = \frac{ma^2bce^{-gn}}{rg}.$$ The endemic equilibrium is therefore:
# $$\bar{X}= \frac{R_{0}-1}{R_{0} + \frac{ac}{g}}, \quad \bar{Z} = \frac{ac\bar{X}}{g + ac\bar{X}} e^{-gn}.$$
#-
# From these equations, the nonlinear relationships that characterize observed malaria transmission can be
# understood. We implement them as follows:

# ## Diagram of Systems
# First we must construct a diagram of systems which describes the interaction between the mosquito and
# host populations. The arrows between the two subsystems represents the bidirectional infection during bloodmeals.

bloodmeal = WiringDiagram([], [:mosquitos, :humans])
mosq_box   = add_box!(bloodmeal, Box(:mosquitos, [:x], [:z]))
human_box  = add_box!(bloodmeal, Box(:humans, [:z], [:x]))
output_box = output_id(bloodmeal)

add_wires!(bloodmeal, Pair[
    (mosq_box, 1)  => (human_box, 1),
    (human_box, 1) => (mosq_box, 1),
    (mosq_box, 1)  => (output_box, 1),
    (human_box, 1) => (output_box, 2)]
)

to_graphviz(bloodmeal)

# ## ODE Model
# Next we implement the concrete mosquito and host dynamics given in Equation (1), and apply them to the diagram
# of systems. This composition is the complete Ross-Macdonald model which we can solve and plot.

dZdt = function(u,x,p,t)
    Z = u[1]
    X = x[1]
    [p.a*p.c*X*(exp(-p.g*p.n) - Z) - p.g*Z]
end

dXdt = function(u,x,p,t)
    X = u[1]
    Z = x[1]
    [p.m*p.a*p.b*Z*(1 - X) - p.r*X]
end

#-

mosquito_model = ContinuousMachine{Float64}(1, 1, 1, dZdt, (u,p,t) -> u)
human_model    = ContinuousMachine{Float64}(1, 1, 1, dXdt, (u,p,t) ->  u)

malaria_model = oapply(bloodmeal,
    Dict(:humans => human_model, :mosquitos => mosquito_model)
)

#-
params = LVector(a = 0.3, b = 0.55, c = 0.15,
    g = 0.1, n = 10, r = 1.0/200, m = 0.5)

u0 = [0.1, 0.3]
tspan = (0.0, 365.0*2)

prob = ODEProblem(malaria_model, u0, tspan, params)
sol = solve(prob, Tsit5());
#-

plot(sol, malaria_model,
    lw=2, title = "Ross-Macdonald Malaria model",
    xlabel = "time", ylabel = "proportion infectious",
    color = ["magenta" "blue"]
)

## Plot the equilibrium behavior as well
a, b, c, g, n, r, m = params
R0 = (m*a^2*b*c*exp(-g*n))/(r*g)
X̄ = (R0 - 1)/(R0 + (a*c)/g)
Z̄ = (a*c*X̄)/(g + a*c*X̄)*exp(-g*n)

N = length(sol)
plot!(sol.t, fill(X̄, N), label = "human equilibrium", ls = :dash, lw = 2, color = "blue")
plot!(sol.t, fill(Z̄, N), label = "mosquito equilibrium", ls = :dash, lw = 2, color = "magenta")

# ## Delay Model
# The previous models did not capture the incubation period for the disease in the
# mosquito population. To do so we can replace the models with delay differential equations
# and apply them to the same diagram of systems representing the bloodmeal.

dzdt_delay = function(u,x,h,p,t)
    Y, Z = u
    Y_delay, Z_delay = h(p, t - p.n)
    X, X_delay = x[1]

    [p.a*p.c*X*(1 - Y - Z) -
        p.a*p.c*X_delay*(1 - Y_delay - Z_delay)*exp(-p.g*p.n) -
        p.g*Y,
    p.a*p.c*X_delay*(1 - Y_delay - Z_delay)*exp(-p.g*p.n) -
        p.g*Z]
end

dxdt_delay = function(u,x,h,p,t)
    X, = u
    Z, _ = x[1]
    [p.m*p.a*p.b*Z*(1 - X) - p.r*X]
end

#-

mosquito_delay_model = DelayMachine{Float64, 2}(
    1, 2, 1, dzdt_delay, (u,h,p,t) -> [[u[2], h(p,t - p.n)[2]]])

human_delay_model = DelayMachine{Float64, 2}(
    1, 1, 1, dxdt_delay, (u,h,p,t) -> [[u[1], h(p, t - p.n)[1]]])

malaria_delay_model = oapply(bloodmeal,
    Dict(:humans => human_delay_model, :mosquitos => mosquito_delay_model)
)
#-
params = LVector(a = 0.3, b = 0.55, c = 0.15,
    g = 0.1, n = 10, r = 1.0/200, m = 0.5)

u0_delay = [0.09, .01, 0.3]
tspan = (0.0, 365.0*3)
hist(p,t) = u0_delay;

prob = DDEProblem(malaria_delay_model, u0_delay, [], hist, tspan, params)
alg = MethodOfSteps(Tsit5())
sol = solve(prob, alg)

plot(sol, label=["non-infectious mosquito population" "infectious mosquito population" "host population"],
    lw=2, title = "Ross-Macdonald malaria model",
    xlabel = "time", ylabel = "proportion infectious",
    color = ["magenta" "red" "blue"]
)

# While the equilibrium points of the two models are identical, they exhibit different dynamical behavior
# before settling down to equilibrium. Because models are often used to examine how the system may respond
# to intervention, incorporating additional biological realism can produce more plausible results regarding
# the expected time for interventions to affect the system. In any case, such simple models are best used as
# tools to explain basic concepts in malaria epidemiology rather than descriptions of real systems, which are
# likely far from equilibrium and affected by weather, climate, policy, and other external forces.


text:
using Catlab.WiringDiagrams

f = singleton_diagram(Box(:f, [:A], [:A]))
g = singleton_diagram(Box(:g, [:B], [:B]))
h = singleton_diagram(Box(:h, [:A,:B], [:C]))
diagram = compose(otimes(f,g),h)


text:
using Test
using ColorTypes: RGB, Gray, N0f8
using FileIO, ImageCore, ImageMagick, Dates, Statistics
using Statistics, StatsBase

import VideoIO

createmode = false

testdir = dirname(@__FILE__)
videodir = joinpath(testdir, "..", "videos")

VideoIO.TestVideos.available()
VideoIO.TestVideos.download_all()

swapext(f, new_ext) = "$(splitext(f)[1])$new_ext"

isarm() = Base.Sys.ARCH in (:arm,:arm32,:arm7l,:armv7l,:arm8l,:armv8l,:aarch64,:arm64)

#@show Base.Sys.ARCH

@noinline function isblank(img)
    all(c->green(c) == 0, img) || all(c->blue(c) == 0, img) || all(c->red(c) == 0, img) || maximum(rawview(channelview(img))) < 0xcf
end

@testset "Reading of various example file formats" begin
    for name in VideoIO.TestVideos.names()
        @testset "Reading $name" begin
            first_frame_file = joinpath(testdir, swapext(name, ".png"))
            !createmode && (first_frame = load(first_frame_file))

            f = VideoIO.testvideo(name)
            v = VideoIO.openvideo(f)

            time_seconds = VideoIO.gettime(v)
            @test time_seconds == 0

            if !createmode && (size(first_frame, 1) > v.height)
                first_frame = first_frame[1+size(first_frame,1)-v.height:end,:]
            end

            # Find the first non-trivial image
            img = read(v)
            i=1
            while isblank(img)
                read!(v, img)
                i += 1
            end
            # println("$name vs. $first_frame_file - First non-blank frame: $i") # for debugging
            createmode && save(first_frame_file,img)
            if isarm()
                !createmode && (@test_skip img == first_frame)
            else
                !createmode && (@test img == first_frame)
            end

            for i in 1:50
                read!(v,img)
            end
            fiftieth_frame = img
            timebase = v.avin.video_info[1].stream.time_base
            tstamp = v.aVideoFrame[1].pkt_dts
            video_tstamp = v.avin.video_info[1].stream.first_dts
            fiftytime = (tstamp-video_tstamp)/(convert(Float64,timebase.den)/convert(Float64,timebase.num))

            while !eof(v)
                read!(v, img)
            end

            seek(v,float(fiftytime))
            read!(v,img)

            @test img == fiftieth_frame

            # read first frames again, and compare
            seekstart(v)

            read!(v, img)

            while isblank(img)
                read!(v, img)
            end

            if isarm()
                !createmode && (@test_skip img == first_frame)
            else
                !createmode && (@test img == first_frame)
            end

            close(v)
        end
    end
end

@testset "IO reading of various example file formats" begin
    for name in VideoIO.TestVideos.names()
        # TODO: fix me?
        (startswith(name, "ladybird") || startswith(name, "NPS")) && continue
        @testset "Testing $name" begin
            first_frame_file = joinpath(testdir, swapext(name, ".png"))
            first_frame = load(first_frame_file)

            filename = joinpath(videodir, name)
            v = VideoIO.openvideo(open(filename))

            if size(first_frame, 1) > v.height
                first_frame = first_frame[1+size(first_frame,1)-v.height:end,:]
            end
            img = read(v)
            # Find the first non-trivial image
            while isblank(img)
                read!(v, img)
            end

            if isarm()
                @test_skip img == first_frame
            else
                @test img == first_frame
            end
            while !eof(v)
                read!(v, img)
            end
        end
    end

    VideoIO.testvideo("ladybird") # coverage testing
    @test_throws ErrorException VideoIO.testvideo("rickroll")
    @test_throws ErrorException VideoIO.testvideo("")
end
@testset "Reading video metadata" begin
    @testset "Reading Storage Aspect Ratio: SAR" begin
        # currently, the SAR of all the test videos is 1, we should get another video with a valid SAR that is not equal to 1
        vids = Dict("ladybird.mp4" => 1, "black_hole.webm" => 1, "crescent-moon.ogv" => 1, "annie_oakley.ogg" => 1)
        @test all(VideoIO.aspect_ratio(VideoIO.openvideo(joinpath(videodir, k))) == v for (k,v) in vids)
    end
    @testset "Reading video duration, start date, and duration" begin
        # tesing the duration and date & time functions:
        file = joinpath(videodir, "annie_oakley.ogg")
        @test VideoIO.get_duration(file) == 24224200/1e6
        @test VideoIO.get_start_time(file) == DateTime(1970, 1, 1)
        @test VideoIO.get_time_duration(file) == (DateTime(1970, 1, 1), 24224200/1e6)
    end
end

@testset "Encoding video across all supported colortypes" begin
    for el in [UInt8, RGB{N0f8}]
        @testset "Encoding $el imagestack" begin
            imgstack = map(x->rand(el,100,100),1:100)
            props = [:priv_data => ("crf"=>"22","preset"=>"medium")]
            encodedvideopath = VideoIO.encodevideo("testvideo.mp4",imgstack,framerate=30,AVCodecContextProperties=props, silent=true)
            @test stat(encodedvideopath).size > 100
            rm(encodedvideopath)
        end
    end
end

@testset "Video encode/decode accuracy (read, encode, read, compare)" begin
    file = joinpath(videodir, "annie_oakley.ogg")
    f = VideoIO.openvideo(file)
    imgstack_rgb = []
    imgstack_gray = []
    while !eof(f)
        img = collect(read(f))
        img_gray = convert(Array{Gray{N0f8}},img)
        push!(imgstack_rgb,img)
        push!(imgstack_gray,img_gray)
    end
    @testset "Lossless Grayscale encoding" begin
        file_lossless_gray_copy = joinpath(videodir, "annie_oakley_lossless_gray.mp4")
        prop = [:color_range=>2, :priv_data => ("crf"=>"0","preset"=>"medium")]
        codec_name="libx264"
        VideoIO.encodevideo(file_lossless_gray_copy,imgstack_gray,codec_name=codec_name,AVCodecContextProperties=prop, silent=true)

        fcopy = VideoIO.openvideo(file_lossless_gray_copy,target_format=VideoIO.AV_PIX_FMT_GRAY8)
        imgstack_gray_copy = []
        while !eof(fcopy)
            push!(imgstack_gray_copy,collect(read(fcopy)))
        end
        close(f)
        @test eltype(imgstack_gray) == eltype(imgstack_gray_copy)
        @test length(imgstack_gray) == length(imgstack_gray_copy)
        @test size(imgstack_gray[1]) == size(imgstack_gray_copy[1])
        @test !any(.!(imgstack_gray .== imgstack_gray_copy))
    end

    @testset "Lossless RGB encoding" begin
        file_lossless_rgb_copy = joinpath(videodir, "annie_oakley_lossless_rgb.mp4")
        prop = [:priv_data => ("crf"=>"0","preset"=>"medium")]
        codec_name="libx264rgb"
        VideoIO.encodevideo(file_lossless_rgb_copy,imgstack_rgb,codec_name=codec_name,AVCodecContextProperties=prop, silent=true)

        fcopy = VideoIO.openvideo(file_lossless_rgb_copy)
        imgstack_rgb_copy = []
        while !eof(fcopy)
            img = collect(read(fcopy))
            push!(imgstack_rgb_copy,img)
        end
        close(f)
        @test eltype(imgstack_rgb) == eltype(imgstack_rgb_copy)
        @test length(imgstack_rgb) == length(imgstack_rgb_copy)
        @test size(imgstack_rgb[1]) == size(imgstack_rgb_copy[1])
        @test !any(.!(imgstack_rgb .== imgstack_rgb_copy))
    end

    @testset "UInt8 accuracy during read & lossless encode" begin
        # Test that reading truth video has one of each UInt8 value pixels (16x16 frames = 256 pixels)
        f = VideoIO.openvideo(joinpath(testdir,"precisiontest_gray_truth.mp4"),target_format=VideoIO.AV_PIX_FMT_GRAY8)
        frame_truth = collect(rawview(channelview(read(f))))
        h_truth = fit(Histogram, frame_truth[:], 0:256)
        @test h_truth.weights == fill(1,256) #Test that reading is precise

        # Test that encoding new test video has one of each UInt8 value pixels (16x16 frames = 256 pixels)
        img = Array{UInt8}(undef,16,16)
        for i in 1:256
            img[i] = UInt8(i-1)
        end
        imgstack = []
        for i=1:24
            push!(imgstack,img)
        end
        props = [:color_range=>2, :priv_data => ("crf"=>"0","preset"=>"medium")]
        VideoIO.encodevideo(joinpath(testdir,"precisiontest_gray_test.mp4"), imgstack,
            AVCodecContextProperties = props,silent=true)
        f = VideoIO.openvideo(joinpath(testdir,"precisiontest_gray_test.mp4"),
            target_format=VideoIO.AV_PIX_FMT_GRAY8)
        frame_test = collect(rawview(channelview(read(f))))
        h_test = fit(Histogram, frame_test[:], 0:256)
        @test h_test.weights == fill(1,256) #Test that encoding is precise (if above passes)
    end

    @testset "Correct frame order when reading & encoding" begin
        @testset "Frame order when reading ground truth video" begin
            # Test that reading a video with frame-incremental pixel values is read in in-order
            f = VideoIO.openvideo(joinpath(testdir,"ordertest_gray_truth.mp4"),target_format=VideoIO.AV_PIX_FMT_GRAY8)
            frame_ids_truth = []
            while !eof(f)
                img = collect(rawview(channelview(read(f))))
                push!(frame_ids_truth,img[1,1])
            end
            @test frame_ids_truth == collect(0:255) #Test that reading is in correct frame order
        end
        @testset "Frame order when encoding, then reading video" begin
            # Test that writing and reading a video with frame-incremental pixel values is read in in-order
            imgstack = []
            img = Array{UInt8}(undef,16,16)
            for i in 0:255
                push!(imgstack,fill(UInt8(i),(16,16)))
            end
            props = [:color_range=>2, :priv_data => ("crf"=>"0","preset"=>"medium")]
            VideoIO.encodevideo(joinpath(testdir,"ordertest_gray_test.mp4"), imgstack,
                AVCodecContextProperties = props,silent=true)
            f = VideoIO.openvideo(joinpath(testdir,"ordertest_gray_test.mp4"),
                target_format=VideoIO.AV_PIX_FMT_GRAY8)
            frame_ids_test = []
            while !eof(f)
                img = collect(rawview(channelview(read(f))))
                push!(frame_ids_test,img[1,1])
            end
            @test frame_ids_test == collect(0:255) #Test that reading is in correct frame order
        end
    end
end

#VideoIO.TestVideos.remove_all()