text:
using Clang.Generators
using Clang.JLLEnvs
using Clang_jll

cd(@__DIR__)

include_dir = joinpath(Clang_jll.artifact_dir, "include") |> normpath
clang_dir = joinpath(include_dir, "clang-c")

for target in JLLEnvs.JLL_ENV_TRIPLES
    @info "processing $target"

    options = load_options(joinpath(@__DIR__, "generator.toml"))
    options["general"]["output_file_path"] = joinpath(@__DIR__, "..", "lib", "$target.jl")

    args = get_default_args(target)
    push!(args, "-I$include_dir")

    headers = detect_headers(clang_dir, args)

    ctx = create_context(headers, args, options)

    build!(ctx)
end


text:
using Documenter, TrigPolys

makedocs(
    sitename="TrigPolys",

    format = Documenter.HTML(
        # See https://github.com/JuliaDocs/Documenter.jl/issues/868
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),

    pages = [
        "Introduction" => "index.md",
    ]
)

deploydocs(
    repo   = "github.com/JuliaAlgebra/TrigPolys.jl.git",
)

ext: jl


text:
using Base.Meta

export @objc, @classes
export objcm # # TODO: remove

callerror() = error("ObjectiveC call: use [obj method] or [obj method:param ...]")

# TODO: Remove? Don't know when this gets used
function flatvcat(ex::Expr)
  any(ex->isexpr(ex, :row), ex.args) || return ex
  flat = Expr(:hcat)
  for row in ex.args
    isexpr(row, :row) ?
      push!(flat.args, row.args...) :
      push!(flat.args, row)
  end
  return calltransform(flat)
end

"""
    calltransform(ex::Expr)

Transform and expression such as `[NSString new]` into

    message(NSString, Selector("new"))

For this to work you would have had to have loaded the class `NSString` first.
This is done with:

    @classes NSString

That line basically translates into the following Julia code:

    NSString = Class("NSString")
"""
function calltransform(ex::Expr)
  obj = objcm(ex.args[1])
  params = ex.args[2:end]
  isempty(params) && callerror()

  # Simple case such as [NSString new]
  # params[1] would be symbol :new
  if typeof(params[1]) == Symbol
    if length(params) > 1
        callerror()
    end
    return :($message($obj, $(Selector(params[1]))))
  end

  # Make sure each parameter is valid
  # Assume params = [setObject:12 forKey: foo]
  ok = all(params) do param
      # assume param = [call:, :, setObject, 12]
      isexpr(param, :call, 3) &&
      param.args[1] == :(:)   &&
      typeof(param.args[2]) == Symbol
  end
  ok || callerror()

  param_names = [param.args[2] for param in params]

  # Turn e.g. [setObject:12 forKey: foo] into "setObject:forKey:"
  selector = map(param_names) do param
    string(param, ':')
  end |> join |> Selector

  values = [objcm(param.args[3]) for param in params]
  :($message($obj, $selector, $(values...)))
end

function objcm(ex::Expr)
    # check if we got an objective-c call such as [NSSound soundNamed:name]
    # an objc call is a hcat of symbols and expressions
    if isexpr(ex, :hcat)
        calltransform(ex)
    elseif isexpr(ex, :vcat) # TODO: Remove? Don't know when this gets used
        flatvcat(ex)
    elseif isexpr(ex, [:block, :let])
        Expr(:block, map(objcm, ex.args)...)
    else
        # This expression did not have any ObjC-like expression which
        # needed transormation, but subexpression may.
        Expr(ex.head, map(objcm, ex.args)...)
    end
end

objcm(ex) = ex

"""
    @objc expr
Interprets `expr` as Objective-C code to call. This creates a new string:

    str = @objc [NSString new]
"""
macro objc(ex)
  esc(objcm(ex))
end

# Import Classes

"""
    @classes klass

Imports the class `klass`. It basically assigns a class object to the
given Julia variable `klass`. So this would lookup the `NSString` class and assign
the corresponding `NSString` class object to the a variable named `NSString`

    @classes NSString

You can load multiple classes:

    @classes NSString, NSBundle, NSArray, NSObject
"""
macro classes(names)
    if typeof(names) == Symbol
        names = [names]
    else
        names = names.args
    end

    classdefs = map(names) do name
        :(const $(esc(name)) = Class($(Expr(:quote, name))))
    end

    Expr(:block, classdefs..., nothing)
end


text:
# Part of submodule Utils of BetaML - the Beta Machine Learning Toolkit
# Various helper/ utility functions concerning stochastiticy management

"""
    FIXEDSEED

Fixed seed to allow reproducible results.
This is the seed used to obtain the same results under unit tests.

Use it with:
- `myAlgorithm(;rng=FIXEDRNG)`             # always produce the same sequence of results on each run of the script ("pulling" from the same rng object on different calls)
- `myAlgorithm(;rng=copy(FIXEDRNG)`        # always produce the same result (new rng object on each call)
"""
const FIXEDSEED = 123

"""
    FIXEDRNG

Fixed ring to allow reproducible results

Use it with:
- `myAlgorithm(;rng=FIXEDRNG)`         # always produce the same sequence of results on each run of the script ("pulling" from the same rng object on different calls)
- `myAlgorithm(;rng=copy(FIXEDRNG))`   # always produce the same result (new rng object on each function call)

"""
const FIXEDRNG  = StableRNG(FIXEDSEED) #StableRNG(FIXEDSEED) Random.default_rng() #MersenneTwister(FIXEDSEED)
#const FIXEDRNG  = MersenneTwister(FIXEDSEED) #StableRNG(FIXEDSEED) Random.default_rng()

"""
    generateParallelRngs(rng::AbstractRNG, n::Integer;reSeed=false)

For multi-threaded models, return n independent random number generators (one per thread) to be used in threaded computations.

Note that each ring is a _copy_ of the original random ring. This means that code that _use_ these RNGs will not change the original RNG state.

Use it with `rngs = generateParallelRngs(rng,Threads.nthreads())` to have a separate rng per thread.
By default the function doesn't re-seed the RNG, as you may want to have a loop index based re-seeding strategy rather than a threadid-based one (to guarantee the same result independently of the number of threads).
If you prefer, you can instead re-seed the RNG here (using the parameter `reSeed=true`), such that each thread has a different seed. Be aware however that the stream  of number generated will depend from the number of threads at run time.
"""
function generateParallelRngs(rng::AbstractRNG, n::Integer;reSeed=false)
    if reSeed
        seeds = [rand(rng,100:18446744073709551615) for i in 1:n] # some RNGs have issues with too small seed
        rngs  = [deepcopy(rng) for i in 1:n]
        return Random.seed!.(rngs,seeds)
    else
        return [deepcopy(rng) for i in 1:n]
    end
end

"""
    shuffle(data;dims,rng)

Shuffle a vector of n-dimensional arrays across dimension `dims` keeping the same order between the arrays

# Parameters
- `data`: The vector of arrays to shuffle
- `dims`: The dimension over to apply the shuffle [def: `1`]
- `rng`:  An `AbstractRNG` to apply for the shuffle

# Notes
- All the arrays must have the same size for the dimension to shuffle

# Example
```
julia> a = [1 2 30; 10 20 30]; b = [100 200 300];
julia> (aShuffled, bShuffled) = shuffle([a,b],dims=2)
2-element Vector{Matrix{Int64}}:
 [1 30 2; 10 30 20]
 [100 300 200]
 ```
"""
function shuffle(data::AbstractArray{T,1};dims=1,rng=Random.GLOBAL_RNG)  where T <: AbstractArray
    Ns = [size(m,dims) for m in data]
    length(Set(Ns)) == 1 || @error "In `shuffle(arrays)` all individual arrays need to have the same size on the dimension specified"
    N    = Ns[1]
    ridx = Random.shuffle(rng, 1:N)
    out = similar(data)
    for (i,a) in enumerate(data)
       aidx = convert(Vector{Union{UnitRange{Int64},Vector{Int64}}},[1:i for i in size(a)])
       #aidx = [collect(1:i) for i in size(a)]
       aidx[dims] = ridx
       out[i] = a[aidx...]
    end
    return out
end
shuffle(rng::AbstractRNG,data::AbstractArray{T,1};dims=1) where T <: AbstractArray = shuffle(data;dims=dims,rng=rng)


text:
function calculate_tgrad(sys::AbstractODESystem;
                         simplify=false)
  isempty(sys.tgrad[]) || return sys.tgrad[]  # use cached tgrad, if possible

  # We need to remove explicit time dependence on the state because when we
  # have `u(t) * t` we want to have the tgrad to be `u(t)` instead of `u'(t) *
  # t + u(t)`.
  rhs = [detime_dvs(eq.rhs) for eq ∈ equations(sys)]
  iv = sys.iv
  xs = states(sys)
  rule = Dict(map((x, xt) -> xt=>x, detime_dvs.(xs), xs))
  rhs = substitute.(rhs, Ref(rule))
  tgrad = [expand_derivatives(ModelingToolkit.Differential(iv)(r), simplify) for r in rhs]
  reverse_rule = Dict(map((x, xt) -> x=>xt, detime_dvs.(xs), xs))
  tgrad = Num.(substitute.(tgrad, Ref(reverse_rule)))
  sys.tgrad[] = tgrad
  return tgrad
end

function calculate_jacobian(sys::AbstractODESystem;
                            sparse=false, simplify=false)
    isempty(sys.jac[]) || return sys.jac[]  # use cached Jacobian, if possible
    rhs = [eq.rhs for eq ∈ equations(sys)]

    iv = sys.iv
    dvs = states(sys)

    if sparse
        jac = sparsejacobian(rhs, dvs, simplify=simplify)
    else
        jac = jacobian(rhs, dvs, simplify=simplify)
    end

    sys.jac[] = jac  # cache Jacobian
    return jac
end

struct ODEToExpr
    sys::AbstractODESystem
    states::Vector
end
ODEToExpr(@nospecialize(sys)) = ODEToExpr(sys,states(sys))
(f::ODEToExpr)(O::Num) = f(value(O))
function (f::ODEToExpr)(O::Term)
    if isa(operation(O), Sym)
        any(isequal(O), f.states) && return tosymbol(O)
        # dependent variables
        return build_expr(:call, Any[operation(O).name; f.(arguments(O))])
    end
    return build_expr(:call, Any[operation(O); f.(arguments(O))])
end
(f::ODEToExpr)(x) = toexpr(x)

function generate_tgrad(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys);
                        simplify=false, kwargs...)
    tgrad = calculate_tgrad(sys,simplify=simplify)
    return build_function(tgrad, dvs, ps, sys.iv;
                          conv = ODEToExpr(sys), kwargs...)
end

function generate_jacobian(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys);
                           simplify=false, sparse = false, kwargs...)
    jac = calculate_jacobian(sys;simplify=simplify,sparse=sparse)
    sub = Dict(value.(dvs) .=> makesym.(value.(dvs)))
    jac = map(d->substitute(d, sub), jac)
    return build_function(jac, dvs, ps, sys.iv;
                          conv = ODEToExpr(sys), kwargs...)
end

function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys); kwargs...)
    # optimization
    dvs′ = makesym.(value.(dvs), states=dvs)
    ps′ = makesym.(value.(ps), states=dvs)

    sub = Dict(dvs .=> dvs′)
    # substitute x(t) by just x
    rhss = [substitute(deq.rhs, sub) for deq ∈ equations(sys)]
    return build_function(rhss, dvs′, ps′, sys.iv;
                          conv = ODEToExpr(sys),kwargs...)
end

function calculate_massmatrix(sys::AbstractODESystem; simplify=false)
    eqs = equations(sys)
    dvs = states(sys)
    M = zeros(length(eqs),length(eqs))
    for (i,eq) in enumerate(eqs)
        if eq.lhs isa Term && operation(eq.lhs) isa Differential
            j = findfirst(x->isequal(tosymbol(x),tosymbol(var_from_nested_derivative(eq.lhs)[1])),dvs)
            M[i,j] = 1
        else
            eq.lhs == 0 || error("Only semi-explicit constant mass matrices are currently supported. Faulty equation: $eq.")
        end
    end
    M = simplify ? ModelingToolkit.simplify.(M) : M
    # M should only contain concrete numbers
    M == I ? I : M
end

jacobian_sparsity(sys::AbstractODESystem) =
    jacobian_sparsity([eq.rhs for eq ∈ equations(sys)],
                      [dv for dv in states(sys)])

function DiffEqBase.ODEFunction(sys::AbstractODESystem, args...; kwargs...)
    ODEFunction{true}(sys, args...; kwargs...)
end

"""
```julia
function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys);
                                     version = nothing, tgrad=false,
                                     jac = false,
                                     sparse = false,
                                     kwargs...) where {iip}
```

Create an `ODEFunction` from the [`ODESystem`](@ref). The arguments `dvs` and `ps`
are used to set the order of the dependent variable and parameter vectors,
respectively.
"""
function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys), u0 = nothing;
                                     version = nothing, tgrad=false,
                                     jac = false,
                                     eval_expression = true,
                                     sparse = false, simplify=false,
                                     kwargs...) where {iip}

    f_gen = generate_function(sys, dvs, ps; expression=Val{eval_expression}, kwargs...)
    f_oop,f_iip = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in f_gen) : f_gen
    f(u,p,t) = f_oop(u,p,t)
    f(du,u,p,t) = f_iip(du,u,p,t)

    if tgrad
        tgrad_gen = generate_tgrad(sys, dvs, ps;
                                   simplify=simplify,
                                   expression=Val{eval_expression}, kwargs...)
        tgrad_oop,tgrad_iip = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in tgrad_gen) : tgrad_gen
        _tgrad(u,p,t) = tgrad_oop(u,p,t)
        _tgrad(J,u,p,t) = tgrad_iip(J,u,p,t)
    else
        _tgrad = nothing
    end

    if jac
        jac_gen = generate_jacobian(sys, dvs, ps;
                                    simplify=simplify, sparse = sparse,
                                    expression=Val{eval_expression}, kwargs...)
        jac_oop,jac_iip = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in jac_gen) : jac_gen
        _jac(u,p,t) = jac_oop(u,p,t)
        _jac(J,u,p,t) = jac_iip(J,u,p,t)
    else
        _jac = nothing
    end

    M = calculate_massmatrix(sys)

    _M = (u0 === nothing || M == I) ? M : ArrayInterface.restructure(u0 .* u0',M)

    ODEFunction{iip}(f,
                     jac = _jac === nothing ? nothing : _jac,
                     tgrad = _tgrad === nothing ? nothing : _tgrad,
                     mass_matrix = _M,
                     jac_prototype = sparse ? similar(sys.jac[],Float64) : nothing,
                     syms = Symbol.(states(sys)))
end

"""
```julia
function DiffEqBase.ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys);
                                     version = nothing, tgrad=false,
                                     jac = false,
                                     sparse = false,
                                     kwargs...) where {iip}
```

Create a Julia expression for an `ODEFunction` from the [`ODESystem`](@ref).
The arguments `dvs` and `ps` are used to set the order of the dependent
variable and parameter vectors, respectively.
"""
struct ODEFunctionExpr{iip} end

function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
                                     ps = parameters(sys), u0 = nothing;
                                     version = nothing, tgrad=false,
                                     jac = false,
                                     linenumbers = false,
                                     sparse = false, simplify=false,
                                     kwargs...) where {iip}

    idx = iip ? 2 : 1
    f = generate_function(sys, dvs, ps; expression=Val{true}, kwargs...)[idx]
    if tgrad
        _tgrad = generate_tgrad(sys, dvs, ps;
                                simplify=simplify,
                                expression=Val{true}, kwargs...)[idx]
    else
        _tgrad = :nothing
    end

    if jac
        _jac = generate_jacobian(sys, dvs, ps;
                                 sparse=sparse, simplify=simplify,
                                 expression=Val{true}, kwargs...)[idx]
    else
        _jac = :nothing
    end

    M = calculate_massmatrix(sys)

    _M = (u0 === nothing || M == I) ? M : ArrayInterface.restructure(u0 .* u0',M)

    jp_expr = sparse ? :(similar($(sys.jac[]),Float64)) : :nothing

    ex = quote
        f = $f
        tgrad = $_tgrad
        jac = $_jac
        M = $_M
        ODEFunction{$iip}(f,
                         jac = jac,
                         tgrad = tgrad,
                         mass_matrix = M,
                         jac_prototype = $jp_expr,
                         syms = $(Symbol.(states(sys))))
    end
    !linenumbers ? striplines(ex) : ex
end

function ODEFunctionExpr(sys::AbstractODESystem, args...; kwargs...)
    ODEFunctionExpr{true}(sys, args...; kwargs...)
end

function DiffEqBase.ODEProblem(sys::AbstractODESystem, args...; kwargs...)
    ODEProblem{true}(sys, args...; kwargs...)
end

"""
```julia
function DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    simplify=false,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...) where iip
```

Generates an ODEProblem from an ODESystem and allows for automatically
symbolically calculating numerical enhancements.
"""
function DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    simplify=false,
                                    linenumbers = true, parallel=SerialForm(),
                                    eval_expression = true,
                                    kwargs...) where iip
    dvs = states(sys)
    ps = parameters(sys)
    u0map′ = lower_mapnames(u0map,sys.iv)
    u0 = varmap_to_vars(u0map′,dvs)

    if !(parammap isa DiffEqBase.NullParameters)
        parammap′ = lower_mapnames(parammap)
        p = varmap_to_vars(parammap′,ps)
    else
        p = ps
    end

    f = ODEFunction{iip}(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,checkbounds=checkbounds,
                        linenumbers=linenumbers,parallel=parallel,simplify=simplify,
                        sparse=sparse,eval_expression=eval_expression,kwargs...)
    ODEProblem{iip}(f,u0,tspan,p;kwargs...)
end

"""
```julia
function DiffEqBase.ODEProblemExpr{iip}(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    linenumbers = true, parallel=SerialForm(),
                                    skipzeros=true, fillzeros=true,
                                    simplify=false,
                                    kwargs...) where iip
```

Generates a Julia expression for constructing an ODEProblem from an
ODESystem and allows for automatically symbolically calculating
numerical enhancements.
"""
struct ODEProblemExpr{iip} end

function ODEProblemExpr{iip}(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    simplify=false,
                                    linenumbers = false, parallel=SerialForm(),
                                    kwargs...) where iip

    dvs = states(sys)
    ps = parameters(sys)
    u0map′ = lower_mapnames(u0map,sys.iv)
    u0 = varmap_to_vars(u0map′,dvs)

    if !(parammap isa DiffEqBase.NullParameters)
        parammap′ = lower_mapnames(parammap)
        p = varmap_to_vars(parammap′,ps)
    else
        p = ps
    end

    f = ODEFunctionExpr{iip}(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,checkbounds=checkbounds,
                        linenumbers=linenumbers,parallel=parallel,
                        simplify=simplify,
                        sparse=sparse,kwargs...)
    ex = quote
        f = $f
        u0 = $u0
        tspan = $tspan
        p = $p
        ODEProblem(f,u0,tspan,p;$(kwargs...))
    end
    !linenumbers ? striplines(ex) : ex
end

function ODEProblemExpr(sys::AbstractODESystem, args...; kwargs...)
    ODEProblemExpr{true}(sys, args...; kwargs...)
end

### Enables Steady State Problems ###
function DiffEqBase.SteadyStateProblem(sys::AbstractODESystem, args...; kwargs...)
    SteadyStateProblem{true}(sys, args...; kwargs...)
end

"""
```julia
function DiffEqBase.SteadyStateProblem(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...) where iip
```
Generates an SteadyStateProblem from an ODESystem and allows for automatically
symbolically calculating numerical enhancements.
"""
function DiffEqBase.SteadyStateProblem{iip}(sys::AbstractODESystem,u0map,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...) where iip
    dvs = states(sys)
    ps = parameters(sys)
    u0 = varmap_to_vars(u0map,dvs)
    p = varmap_to_vars(parammap,ps)
    f = ODEFunction(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,checkbounds=checkbounds,
                        linenumbers=linenumbers,parallel=parallel,
                        sparse=sparse,kwargs...)
    SteadyStateProblem(f,u0,p;kwargs...)
end

"""
```julia
function DiffEqBase.SteadyStateProblemExpr(sys::AbstractODESystem,u0map,tspan,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    skipzeros=true, fillzeros=true,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...) where iip
```
Generates a Julia expression for building a SteadyStateProblem from
an ODESystem and allows for automatically symbolically calculating
numerical enhancements.
"""
struct SteadyStateProblemExpr{iip} end

function SteadyStateProblemExpr{iip}(sys::AbstractODESystem,u0map,
                                    parammap=DiffEqBase.NullParameters();
                                    version = nothing, tgrad=false,
                                    jac = false,
                                    checkbounds = false, sparse = false,
                                    linenumbers = true, parallel=SerialForm(),
                                    kwargs...) where iip
    dvs = states(sys)
    ps = parameters(sys)
    u0 = varmap_to_vars(u0map,dvs)
    p = varmap_to_vars(parammap,ps)
    f = ODEFunctionExpr(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,checkbounds=checkbounds,
                        linenumbers=linenumbers,parallel=parallel,
                        sparse=sparse,kwargs...)
    ex = quote
        f = $f
        u0 = $u0
        p = $p
        SteadyStateProblem(f,u0,p;$(kwargs...))
    end
    !linenumbers ? striplines(ex) : ex
end

function SteadyStateProblemExpr(sys::AbstractODESystem, args...; kwargs...)
    SteadyStateProblemExpr{true}(sys, args...; kwargs...)
end

isdifferential(expr) = istree(expr) && operation(expr) isa Differential
isdiffeq(eq) = isdifferential(eq.lhs)


text:
include("status.jl")
include("card.jl")

struct Player
    name::String
    card::Card
    statusEffects::Vector{Any}
    alive::Bool
end

Player(name) = Player(name, Card("", ""), [ ], true)

function Player_AppendStatus!(player::Player, status::Status)
    push!(player.statusEffects, status)
end

max_stars_repo_name: Flobii/discord-bot

max_issues_repo_name: Flobii/discord-bot

max_forks_repo_name: Flobii/discord-bot


text:
"""
  Rudimentary adaptive function sampling
  Adapted from https://stackoverflow.com/a/14084715
"""
function sample_function(func, nodes, values=func.(nodes), mask=Colon();
    depth=0, min_nodes=16, max_level=16, tol=0.01)

  if depth > max_level
    return nodes, values
  end

  # Calculate the function in between the nodes
  x_c = 0.5 * (nodes[1:end-1][mask] + nodes[2:end][mask])
  y_c = func.(x_c)

  x₂ = vcat(nodes, x_c)
  y₂ = vcat(values, y_c)
  is = sortperm(x₂)

  x₂ = x₂[is]
  y₂ = y₂[is]

  # Determine where to refine
  if length(x₂) < min_nodes
    mask = trues(length(x₂) - 1)
  else
    # Rescale and compute the length of each line segment in the path
    dx = diff(x₂) / ptp(x₂)
    dy = diff(y₂) / abs(ptp(y₂))
    ds = sqrt.(abs2.(dx) + abs2.(dy))

    # Compute the angle between consecutive line segments θ = acos(a⋅b / (|a||b|))
    dx ./= ds
    dy ./= ds

    if eltype(dy) <: Complex
      @views dy_re = reinterpret(Float64, dy)[1:2:end]
      @views dy_im = reinterpret(Float64, dy)[2:2:end]
      dθ = acos.(clamp.(dx[2:end] .* dx[1:end-1] + dy_re[2:end] .* dy_re[1:end-1] + dy_im[2:end] .* dy_im[1:end-1], -1, 1))
    else
      dθ = acos.(clamp.(dx[2:end] .* dx[1:end-1] + dy[2:end] .* dy[1:end-1], -1, 1))
    end

    # Determine where to subdivide (total length of the path is computed ≈ to accuracy `tol`)
    mask = 0.5 * dθ .* (ds[2:end] .+ ds[1:end-1]) .> tol * sum(ds)
    mask = [mask; false]
    mask[2:end] .|= mask[1:end-1]
  end

  # Refine!
  if any(mask)
    return sample_function(func, x₂, y₂, mask; tol=tol, depth=depth+1, min_nodes=min_nodes, max_level=max_level)
  else
    return x₂, y₂
  end
end


text:
# module Bukdu

immutable Scope
    path::String
    host::String
    private::Dict{Symbol,String}
    assigns::Dict{Symbol,String}
end

module RouterScope

import ..Bukdu: ApplicationController
import ..Bukdu: Route, RouterRoute
import ..Bukdu: Scope
import ..Bukdu: Keyword, Naming

stack = Vector{Scope}()

function push_scope!(options::Dict)
    path = Keyword.get(options, :path, "")::String
    path = validate_path(path)
    host = Keyword.get(options, :host, "")::String
    private = Dict{Symbol,String}(Keyword.get(options, :private, Dict()))
    assigns = Dict{Symbol,String}(Keyword.get(options, :assigns, Dict()))
    scope = Scope(path, host, private, assigns)
    push!(RouterScope.stack, scope)
end

function pop_scope!()
    pop!(RouterScope.stack)
end

function find_host(stack::Vector{Scope})
    for scope in stack
        !isempty(scope.host) && return scope.host
    end
    return ""
end

function route{AC<:ApplicationController}(kind::Symbol, verb::Function,
             path::String, ::Type{AC}, action::Function, options::Dict)::Route
    path    = validate_path(path)
    stack = get_stack()
    private = Dict{Symbol,String}(reduce(merge, vcat(extract(stack, :private), Keyword.get(options, :private, Dict()))))
    assigns = Dict{Symbol,String}(reduce(merge, vcat(extract(stack, :assigns), Keyword.get(options, :assigns, Dict()))))
    RouterRoute.build(kind, verb, join_path(stack, path), find_host(stack), AC, action, private, assigns)
end

function validate_path(path)::String
    path
end

function join_path(stack::Vector{Scope}, path::String)
    string(join(extract(stack, :path)), path)
end

function extract(stack::Vector{Scope}, attr::Symbol)
    map(scope->getfield(scope, attr), stack)
end

function get_stack()
    getfield(RouterScope, :stack)
end

end # module Bukdu.RouterScope


text:
using Test
using StaticArrays
using Tri_Tet_Intersections
using LinearAlgebra
using NumericalTricks

include("test_exports.jl")
include("test_static_clip.jl")
include("test_plane_tet_intersection.jl")
include("test_poly_eight.jl")
include("test_utility.jl")
include("test_quadrature.jl")


text:
using PkgBenchmark
using TimeZones
import Dates: DateFormat
import TimeZones.TZData: parse_components

@benchgroup "parse" begin
    @bench "components" parse_components($("2010 Europe/Warsaw"), $(DateFormat("yyyy ZZZ")))

    # https://github.com/JuliaTime/TimeZones.jl/issues/25
    function parse_dates(n)
        df = DateFormat("yyyymmddHH:MM:SS ZZZ")
        arr = Array{ZonedDateTime}(n)
        for (i, s) in enumerate(Iterators.repeated("2016060701:02:03 America/Toronto", n))
            arr[i] = ZonedDateTime(s,df)
        end
        return arr
    end

    @bench "ZonedDateTime" parse_dates(1000)
end


text:
using TimeZones: UTCOffset, value, isdst

@test value(UTCOffset(0, 0)) == 0
@test value(UTCOffset(3600, 0)) == 3600
@test value(UTCOffset(0, 3600)) == 3600
@test value(UTCOffset(-7200, 3600)) == -3600

@test string(UTCOffset(0, 0)) == "+00:00"
@test string(UTCOffset(3600, 0)) == "+01:00"
@test string(UTCOffset(0, 3600)) == "+01:00"
@test string(UTCOffset(-7200, 3600)) == "-01:00"
@test string(UTCOffset(-3661)) == "-01:01:01"
@test string(UTCOffset( 1800)) == "+00:30"
@test string(UTCOffset(-1800)) == "-00:30"

@test !isdst(UTCOffset(0))
@test !isdst(UTCOffset(0, 0))
@test isdst(UTCOffset(0, 3600))
@test isdst(UTCOffset(0, 7200))

# Arithmetic
dt = DateTime(2015, 1, 1)
@test dt + UTCOffset(7200, -3600) == dt + Hour(1)
@test dt - UTCOffset(7200, -3600) == dt - Hour(1)

# Comparisons
let a = UTCOffset(0, 0), b = UTCOffset(0, 0)
    @test a == b
    @test isequal(a, b)
end

let a = UTCOffset(0, 3600), b = UTCOffset(3600, 0)
    @test a != b
    @test !isequal(a, b)
end

let a = UTCOffset(7200, 3600), b = UTCOffset(3600, 7200)
    @test a != b

    # Treated as equal since some time zone representations adjust the standard
    # offset when dealing with midsummer time.
    @test isequal(a, b)
end

@test sprint(show_compact, UTCOffset(0, 0)) == "UTC+0/+0"
@test sprint(show_compact, UTCOffset(3600, 7200)) == "UTC+1/+2"

# Added: https://github.com/JuliaLang/julia/pull/30817
# Reverted for v1.2.0-rc1: https://github.com/JuliaLang/julia/pull/31727
# Reverted for v1.3.0-rc2: https://github.com/JuliaLang/julia/pull/32973
# Active discussion: https://github.com/JuliaLang/julia/pull/33178
if v"1.2.0-DEV.223" <= VERSION < v"1.2.0-pre.39" || v"1.3-DEV" <= VERSION < v"1.3.0-rc1.33" || v"1.4-DEV" <= VERSION
    @test sprint(show, UTCOffset(0, 0)) == "UTCOffset(Second(0), Second(0))"
    @test sprint(show, UTCOffset(3600, 7200)) == "UTCOffset(Second(3600), Second(7200))"
else
    @test sprint(show, UTCOffset(0, 0)) == "UTCOffset(0 seconds, 0 seconds)"
    @test sprint(show, UTCOffset(3600, 7200)) == "UTCOffset(3600 seconds, 7200 seconds)"
end

@test sprint(show, MIME("text/plain"), UTCOffset(0, 0)) == "UTC+0/+0"
@test sprint(show, MIME("text/plain"), UTCOffset(3600, 7200)) == "UTC+1/+2"


text:
## Main function to parse a template This works in several steps: each
## character parsed into a token, the tokens are squashed, then
## nested, then rendered.
function parse(template, tags = ("{{", "}}"))
    tokens = make_tokens(template, tags)
    out = nestTokens(tokens)
    MustacheTokens(out)
end

"""
    mt"string"

Macro to parse tokens from a string.

"""
macro mt_str(s)
    parse(s)
end

macro mt_mstr(s)
    parse(s)
end

## Dict for storing parsed templates
TEMPLATES = Dict{AbstractString, MustacheTokens}()

"""
    Mustache.load(filepath, args...)

Load file specified through  `filepath` and return the compiled tokens.
Tokens are memoized for efficiency,

Additional arguments are passed to `Mustache.parse` (for adjusting the tags).
"""
function load(filepath, args...)

    isfile(filepath) || throw(ArgumentError("File $filepath not found"))

    key = string(mtime(filepath)) * filepath * string(hash(args))
    haskey(TEMPLATES,key) && return  TEMPLATES[key]

    open(filepath) do s
        global tpl = parse(read(s, String), args...)
    end

    TEMPLATES[key] = tpl
    tpl

end

# old names. To be deprecated?
const template_from_file = load


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

function QQplot(rs, data)
    # Example usage:
    # x = 0.5-0.5*log(rand(1,1e3));
    # w = 0.001;
    # rs = sskernel(x, w);
    # QQplot(rs, x)
    #
    # Input argument
    # rs: Estimated intensity function
    # data: Spike time

    R=cumsum(rs*dt);

    x_step = data/dt;
    x_round = Int.(ceil(transpose(x_step))+1);

    R_x = R[x_round][:];
    tau = diff(R_x);

    u_hat = 1-exp(-tau);
    u_hat_sorted = sort(u_hat);

    n = length(u_hat);
    u = ([1:n]-0.5)/n;

    cz =linspace(0,1,100);
    confb95 = 1.96*sqrt(cz.*(1-cz)/n);
    confb99 = 2.575*sqrt(cz.*(1-cz)/n);
    p_x = [0,1]; p_y = [0,1];
    plot(p_x, p_y, color="black", linewidth=2.0, linestyle=:dash, label="")
    plot!(cz, cz-confb95, color="black", linewidth=2.0, linestyle=:solid, label="")
    plot!(cz, cz+confb95, color="black", linewidth=2.0, linestyle=:solid, label="")

    plot!(u, u_hat_sorted, linewidth=0, marker="", label="")
end

max_stars_repo_name: ganow/julia-point-process

max_issues_repo_name: ganow/julia-point-process

max_forks_repo_name: ganow/julia-point-process


text:
using ReviseLoop
using Test

@testset "ReviseLoop.jl" begin
    # Write your own tests here.
end


text:
using DSGE, ModelConstructors, Dates, Optim, Test, Nullables, HDF5
path = dirname(@__FILE__)
m = AnSchorfheide()
m <= Setting(:date_zlb_end, get_setting(m, :date_zlb_start))
m <= Setting(:hessian_path, "")

# Dates, indices, number of periods for each regime
@testset "Test field access functions for AbstractDSGEModel objects" begin
    @test DSGE.get_parameters(m) == m.parameters
    @test DSGE.get_rng(m) == m.rng
    @test DSGE.get_settings(m) == m.settings
    @test DSGE.get_observables(m) == m.observables
    @test DSGE.get_pseudo_observables(m) == m.pseudo_observables
    @test DSGE.get_exogenous_shocks(m) == m.exogenous_shocks
end

@testset "Test setting access functions for AbstractDSGEModel objects" begin
    n_anticipated_shocks(m) == get_setting(m, :n_anticipated_shocks)
    n_anticipated_shocks_padding(m) == get_setting(m, :n_anticipated_shocks_padding)

    @test DSGE.date_presample_start(m) == get_setting(m, :date_presample_start)
    @test DSGE.date_mainsample_start(m) == get_setting(m, :date_mainsample_start)
    @test DSGE.date_zlb_start(m) == get_setting(m, :date_zlb_start)
    @test DSGE.date_zlb_end(m) == get_setting(m, :date_zlb_end)

    @test DSGE.date_presample_end(m) == Dates.lastdayofquarter(get_setting(m, :date_mainsample_start) - Dates.Month(3))
    @test DSGE.date_prezlb_end(m) == Dates.lastdayofquarter(get_setting(m, :date_zlb_start) - Dates.Month(3))
    @test DSGE.date_mainsample_end(m) == Dates.lastdayofquarter(get_setting(m, :date_forecast_start) - Dates.Month(3))
    @test DSGE.date_conditional_end(m) == get_setting(m, :date_conditional_end)

    @test DSGE.index_presample_start(m) == 1
    @test DSGE.index_mainsample_start(m) == subtract_quarters(DSGE.date_mainsample_start(m), DSGE.date_presample_start(m)) + 1
    @test DSGE.index_zlb_start(m) == subtract_quarters(DSGE.date_zlb_start(m), DSGE.date_presample_start(m)) + 1
    @test DSGE.index_forecast_start(m) == subtract_quarters(DSGE.date_forecast_start(m), DSGE.date_presample_start(m)) + 1

    @test DSGE.index_shockdec_start(m) == subtract_quarters(DSGE.date_shockdec_start(m), DSGE.date_mainsample_start(m)) + 1

    @test DSGE.index_shockdec_end(m) == subtract_quarters(DSGE.date_shockdec_end(m), DSGE.date_mainsample_start(m)) + 1

    @test n_presample_periods(m)   == subtract_quarters(DSGE.date_mainsample_start(m), DSGE.date_presample_start(m))
    @test n_prezlb_periods(m)      == subtract_quarters(DSGE.date_zlb_start(m), DSGE.date_mainsample_start(m))
    @test n_zlb_periods(m)         == subtract_quarters(DSGE.date_forecast_start(m), DSGE.date_zlb_start(m))
    @test n_mainsample_periods(m)  == subtract_quarters(DSGE.date_forecast_start(m), DSGE.date_mainsample_start(m))
    @test n_conditional_periods(m) == subtract_quarters(DSGE.date_conditional_end(m), DSGE.date_mainsample_end(m))

    @test inds_presample_periods(m)  == collect(DSGE.index_presample_start(m):(DSGE.index_mainsample_start(m)-1))
    @test inds_prezlb_periods(m)     == collect(DSGE.index_mainsample_start(m):(DSGE.index_zlb_start(m)-1))
    @test inds_zlb_periods(m)        == collect(DSGE.index_zlb_start(m):(DSGE.index_forecast_start(m)-1))
    @test inds_mainsample_periods(m) == collect(DSGE.index_mainsample_start(m):(DSGE.index_forecast_start(m)-1))

    @test cond_vintage(m)    == get_setting(m, :cond_vintage)
    @test cond_id(m)         == get_setting(m, :cond_id)
    @test cond_full_names(m) == get_setting(m, :cond_full_names)
    @test cond_semi_names(m) == get_setting(m, :cond_semi_names)
    @test use_population_forecast(m) == get_setting(m, :use_population_forecast)
    @test DSGE.hpfilter_population(m)     == get_setting(m, :hpfilter_population)

    # Interface for general computation settings
    @test DSGE.use_parallel_workers(m)    == get_setting(m, :use_parallel_workers)

    # Interface for estimation settings
    @test DSGE.reoptimize(m)          == get_setting(m, :reoptimize)
    @test DSGE.calculate_hessian(m) == get_setting(m, :calculate_hessian)
    @test DSGE.hessian_path(m)      == get_setting(m, :hessian_path)
    @test DSGE.n_hessian_test_params(m) == get_setting(m, :n_hessian_test_params)

    # Interface for Metropolis-Hastings settings
    @test DSGE.n_mh_blocks(m)       ==  get_setting(m, :n_mh_blocks)
    @test DSGE.n_mh_param_blocks(m) ==  get_setting(m, :n_mh_param_blocks)
    @test DSGE.n_mh_simulations(m)  ==  get_setting(m, :n_mh_simulations)
    @test DSGE.n_mh_burn(m)         ==  get_setting(m, :n_mh_burn)
    @test DSGE.mh_thin(m)           ==  get_setting(m, :mh_thin)

    # Interface for forecast settings
    @test DSGE.date_forecast_start(m)   == get_setting(m, :date_forecast_start)
    @test DSGE.forecast_block_size(m)   == get_setting(m, :forecast_block_size)
    @test DSGE.forecast_start_block(m)  == get_setting(m, :forecast_start_block)
    @test DSGE.forecast_input_file_overrides(m) == get_setting(m, :forecast_input_file_overrides)
    @test isnull(DSGE.forecast_uncertainty_override(m)) && isnull(get_setting(m, :forecast_uncertainty_override))
    @test DSGE.forecast_smoother(m)     == get_setting(m, :forecast_smoother)
    @test DSGE.forecast_tdist_df_val(m) == get_setting(m, :forecast_tdist_df_val)
    @test DSGE.forecast_tdist_shocks(m) == get_setting(m, :forecast_tdist_shocks)
    @test DSGE.forecast_zlb_value(m)    == get_setting(m, :forecast_zlb_value)
    @test DSGE.impulse_response_horizons(m) == get_setting(m, :impulse_response_horizons)
    @test DSGE.n_shockdec_periods(m)    == DSGE.index_shockdec_end(m) - DSGE.index_shockdec_start(m) + 1

    # Interface for alternative policy settings
    @test DSGE.alternative_policy(m) == get_setting(m, :alternative_policy)
end

@testset "Test other auxiliary setting functions for AbstractDSGEModel objects" begin
    m <= Setting(:population_mnemonic, :test)
    @test get(DSGE.parse_population_mnemonic(m)[1]) == :test
    m <= Setting(:population_mnemonic, Nullable())
    println(DSGE.parse_population_mnemonic(m))
    @test DSGE.parse_population_mnemonic(m)[1].hasvalue == false
    m1 = Model1002("ss10")
    @test @test_matrix_approx_eq inds_states_no_ant(m) collect(1:8)
    @test @test_matrix_approx_eq inds_states_no_ant(m1) vcat(collect(1:62), collect(69:84))
    @test @test_matrix_approx_eq inds_shocks_no_ant(m) collect(1:3)
    @test @test_matrix_approx_eq inds_shocks_no_ant(m1) collect(1:18)
    @test @test_matrix_approx_eq inds_obs_no_ant(m) collect(1:3)
    @test @test_matrix_approx_eq inds_obs_no_ant(m1) collect(1:13)

    m <= Setting(:date_forecast_start, DSGE.quartertodate("2019-Q3"))
    m <= Setting(:forecast_horizons, 10)
    exp_date_forecast_end = DSGE.lastdayofquarter(DSGE.quartertodate("2019-Q3") + Dates.Month(3 * 9))
    @test forecast_horizons(m) == 10
    m <= Setting(:date_conditional_end, DSGE.date_mainsample_end(m) + Dates.Month(3))
    @test forecast_horizons(m; cond_type = :full) == 9
    @test forecast_horizons(m; cond_type = :semi) == 9
    @test DSGE.date_forecast_end(m) == exp_date_forecast_end
    @test DSGE.date_shockdec_start(m) == date_mainsample_start(m)
    m <= Setting(:shockdec_startdate, nothing)
    @test DSGE.date_shockdec_start(m) == DSGE.date_mainsample_start(m)
    m <= Setting(:shockdec_startdate, Nullable{Date}(DSGE.quartertodate("2019-Q3")))
    @test DSGE.date_shockdec_start(m) == DSGE.quartertodate("2019-Q3")
    m <= Setting(:shockdec_enddate, Nullable{Date}())
    @test DSGE.date_shockdec_end(m) == DSGE.date_forecast_end(m)
    m <= Setting(:shockdec_enddate, Nullable{Date}(DSGE.quartertodate("2019-Q3")))
    @test DSGE.date_shockdec_end(m) == DSGE.quartertodate("2019-Q3")

    @test @test_matrix_approx_eq load_parameters_from_file(m, "$path/reference/load_params.h5") h5read("$path/reference/load_params.h5", "params")
    @test_throws ErrorException load_parameters_from_file(m, "$path/reference/smc.h5")
    @test_throws ErrorException load_parameters_from_file(m, "$path/reference/smc.jld2")
    @test_throws AssertionError load_parameters_from_file(m1, "$path/reference/load_params.h5")
    @test_throws AssertionError load_parameters_from_file(m, "$path/reference/load_params_float32.h5")

    m.parameters[1].value = 1.
    specify_mode!(m, "reference/load_params.h5")
    out_params = map(x -> x.value, m.parameters)
    @test @test_matrix_approx_eq out_params load_parameters_from_file(m, "$path/reference/load_params.h5")

    specify_hessian!(m, "reference/hessian.h5")
    @test get_setting(m, :calculate_hessian) == false
    @test get_setting(m, :hessian_path) == joinpath(dirname("$path"), "test", "reference", "hessian.h5")

    para = deepcopy(map(x -> x.value, m.parameters))
    m.parameters[1].value = 1.
    DSGE.update!(m, para)
    @test m.parameters[1].value != 1.
    @test @test_matrix_approx_eq map(x->x.value,m.parameters) para
end

nothing


text:
#
# Querying data using SQL.
#

type SQLTable
    sname::String
    name::String
    nattrs::Int
    attrnames::Vector{String}
    attrtypes::Vector{Type}
    attrisnulls::Vector{Bool}
    attruniques::Vector{Bool}
    sortattrs::Vector{Int}
end

SQLCollectionQuery(uri::String, tab::SQLTable, decors::Decorations=NO_DECORATIONS) =
    let dom = Domain(Symbol(uri, "#", tab.sname, "#", tab.name)) |> setdecorations(decors)
        Query(
            SQLEntitySig(uri, tab),
            Input(Void),
            Output(dom) |> setoptional() |> setplural())
    end

SQLFieldQuery(uri::String, tab::SQLTable, attr::Int, decors::Decorations=NO_DECORATIONS) =
    let idom = Domain(Symbol(uri, "#", tab.sname, "#", tab.name)),
        odom = Domain(tab.attrtypes[attr]) |> setdecorations(decors),
        optional = tab.attrisnulls[attr]
        Query(
            SQLFieldSig(uri, tab, attr),
            Input(idom),
            Output(odom) |> setoptional(optional))
    end

SQLLinkQuery(uri::String, itab::SQLTable, iattr::Int, otab::SQLTable, oattr::Int, decors::Decorations=NO_DECORATIONS) =
    let idom = Domain(Symbol(uri, "#", itab.sname, "#", itab.name)),
        odom = Domain(Symbol(uri, "#", otab.sname, "#", otab.name)) |> setdecorations(decors),
        optional = itab.attrisnulls[iattr],
        plural = !otab.attruniques[oattr]
        Query(
            SQLLinkSig(uri, itab, iattr, otab, oattr),
            Input(idom),
            Output(odom) |> setoptional(optional) |> setplural(plural))
    end

immutable SQLEntitySig <: AbstractPrimitive
    uri::String
    tab::SQLTable
end

immutable SQLFieldSig <: AbstractPrimitive
    uri::String
    tab::SQLTable
    attr::Int
end

immutable SQLLinkSig <: AbstractPrimitive
    uri::String
    itab::SQLTable
    iattr::Int
    otab::SQLTable
    oattr::Int
end

describe(io::IO, sig::SQLEntitySig) = print(io, sig.tab.name)

describe(io::IO, sig::SQLFieldSig) = print(io, sig.tab.attrnames[sig.attr])

describe(io::IO, sig::SQLLinkSig) = print(io, sig.itab.attrnames[sig.iattr])

function ev(sig::SQLEntitySig, ::Input, oty::Output, iflow::InputFlow)
    dv = sqltable(iflow.ctx, sig.uri, sig.tab)
    col = collection_impl(OneTo(length(dv)), length(iflow))
    return OutputFlow(oty, col)
end

function ev(sig::SQLFieldSig, ::Input, oty::Output, iflow::InputFlow)
    dv = sqltable(iflow.ctx, sig.uri, sig.tab)
    col = dv.cols[sig.attr]
    idxs = values(iflow)
    return OutputFlow(oty, col[idxs])
end

function ev(sig::SQLLinkSig, ::Input, oty::Output, iflow::InputFlow)
    idv = sqltable(iflow.ctx, sig.uri, sig.itab)
    odv = sqltable(iflow.ctx, sig.uri, sig.otab)
    icol = idv.cols[sig.iattr][values(iflow)]
    ocol = odv.cols[sig.oattr]
    unique = sig.otab.attruniques[sig.oattr]
    col =
        if unique
            plain_join_impl(icol, ocol)
        else
            join_impl(icol, ocol)
        end
    return OutputFlow(oty, col)
end

function plain_join_impl(icol::Column{false,false}, ocol::Column{false,false})
    ovals = values(ocol)
    perm = sortperm(ovals)
    kptrs = (1:length(ovals))[perm]
    kvals = ovals[perm]
    vals = Vector{Int}(length(icol))
    cr = cursor(icol)
    while !done(icol, cr)
        next!(icol, cr)
        ptr = searchsortedfirst(kvals, cr[1])
        val = kptrs[ptr]
        vals[cr.pos] = val
    end
    return PlainColumn(vals)
end

function plain_join_impl{IOPT,OOPT}(icol::Column{IOPT,false}, ocol::Column{OOPT,false})
    ovals = values(ocol)
    optrs = Vector{Int}(length(ovals))
    ocr = cursor(ocol)
    while !done(ocol, ocr)
        next!(ocol, ocr)
        if !isempty(ocr)
            optrs[ocr.l] = ocr.pos
        end
    end
    perm = sortperm(ovals)
    kptrs = optrs[perm]
    kvals = ovals[perm]
    vals = Vector{Int}(length(icol.vals))
    cr = cursor(icol)
    while !done(icol, cr)
        next!(icol, cr)
        if !isempty(cr)
            ptr = searchsortedfirst(kvals, cr[1])
            val = kptrs[ptr]
            vals[cr.l] = val
        end
    end
    return Column{IOPT,false}(offsets(icol), vals)
end

function join_impl{IOPT,OOPT}(icol::Column{IOPT,false}, ocol::Column{OOPT,false})
    ovals = values(ocol)
    optrs = Vector{Int}(length(ovals))
    ocr = cursor(ocol)
    while !done(ocol, ocr)
        next!(ocol, ocr)
        if !isempty(ocr)
            optrs[ocr.l] = ocr.pos
        end
    end
    perm = sortperm(ovals)
    kptrs = optrs[perm]
    kvals = ovals[perm]
    size = 0
    cr = cursor(icol)
    while !done(icol, cr)
        next!(icol, cr)
        if !isempty(cr)
            ptr = searchsorted(kvals, cr[1])
            size += length(ptr)
        end
    end
    offs = Vector{Int}(length(icol)+1)
    vals = Vector{Int}(size)
    cr = cursor(icol)
    offs[1] = 1
    n = 1
    while !done(icol, cr)
        next!(icol, cr)
        if !isempty(cr)
            ptr = searchsorted(kvals, cr[1])
            copy!(vals, n, kptrs, ptr.start, length(ptr))
            n += length(ptr)
        end
        offs[cr.pos+1] = n
    end
    return Column{IOPT,true}(offs, vals)
end

function sqltable(ctx::Dict{Symbol,Any}, uri::String, tab::SQLTable)
    if !(:sql in keys(ctx))
        ctx[:sql] = Dict{Any,Any}()
    end
    sqlctx = ctx[:sql]
    if !((uri, :conn) in keys(sqlctx))
        sqlctx[(uri, :conn)] = connect(uri)
    end
    conn = sqlctx[(uri, :conn)]
    if !((uri, :tab, tab.sname, tab.name) in keys(sqlctx))
        sqlbuf = IOBuffer()
        print(sqlbuf, "SELECT ")
        comma = false
        for attrname in tab.attrnames
            if comma
                print(sqlbuf, ", ")
            end
            comma = true
            print(sqlbuf, '"', attrname, '"')
        end
        print(sqlbuf, " FROM ", '"', tab.sname, '"', ".", '"', tab.name, '"')
        if !isempty(tab.sortattrs)
            print(sqlbuf, " ORDER BY ")
            comma = false
            for sortattr in tab.sortattrs
                if comma
                    print(sqlbuf, ", ")
                end
                comma = true
                print(sqlbuf, sortattr)
            end
        end
        sql = String(sqlbuf)
        dv = execute(conn, sql)
        cols = Column[]
        for (k, col) in enumerate(dv.cols)
            if !tab.attrisnulls[k]
                col = PlainColumn(col.vals)
            end
            push!(cols, col)
        end
        dv = DataVector(dv.len, cols)
        sqlctx[(uri, :tab, tab.sname, tab.name)] = dv
    end
    return sqlctx[(uri, :tab, tab.sname, tab.name)]::DataVector
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
# This file is part of Kpax3. License is MIT.

function loglikelihood(C::Matrix{UInt8},
                       cl::Vector{Int},
                       k::Int,
                       support::MCMCSupport)
  loglik = 0.0

  for b in 1:support.m, l in 1:k
    loglik += support.lp[C[cl[l], b], cl[l], b]
  end

  loglik
end

function loglikelihood(C::Matrix{UInt8},
                       cl::Vector{Int},
                       k::Int,
                       v::Vector{Int},
                       n1s::Matrix{Float64},
                       priorC::AminoAcidPriorCol)
  m = size(C, 2)
  loglik = 0.0

  # If array A has dimension (d_{1}, ..., d_{l}, ..., d_{L}), to access
  # element A[i_{1}, ..., i_{l}, ..., i_{L}] it is possible to use the
  # following linear index
  # lidx = i_{1} + d_{1} * (i_{2} - 1) + ... +
  #      + (d_{1} * ... * d_{l-1}) * (i_{l} - 1) + ... +
  #      + (d_{1} * ... * d_{L-1}) * (i_{L} - 1)
  #
  # A[i_{1}, ..., i_{l}, ..., i_{L}] == A[lidx]
  lidx = 0
  for b in 1:m
    for l in 1:k
      lidx = C[cl[l], b] + 4 * (b - 1)
      loglik += logmarglik(n1s[cl[l], b], v[cl[l]], priorC.A[lidx],
                           priorC.B[lidx])
    end
  end

  loglik
end

function logmarglikelihood(cl::Vector{Int},
                           k::Int,
                           lp::Array{Float64},
                           priorC::AminoAcidPriorCol)
  lml = 0.0

  g = 0
  tmp1 = zeros(Float64, 3)
  tmp2 = zeros(Float64, 2)
  for b in 1:size(lp, 3)
    tmp1[1] = priorC.logγ[1]
    tmp1[2] = priorC.logγ[2]
    tmp1[3] = priorC.logγ[3]

    for l in 1:k
      g = cl[l]

      tmp1[1] += lp[1, g, b]
      tmp1[2] += lp[2, g, b]

      tmp2[1] = priorC.logω[k][1] + lp[3, g, b]
      tmp2[2] = priorC.logω[k][2] + lp[4, g, b]

      tmp1[3] += (tmp2[1] > tmp2[2]) ?
                 (tmp2[1] + log1p(exp(tmp2[2] - tmp2[1]))) :
                 (tmp2[2] + log1p(exp(tmp2[1] - tmp2[2])))
    end

    lml += if (tmp1[1] >= tmp1[2]) && (tmp1[1] >= tmp1[3])
             tmp1[1] + log1p(exp(tmp1[2] - tmp1[1]) + exp(tmp1[3] - tmp1[1]))
           elseif (tmp1[2] >= tmp1[1]) && (tmp1[2] >= tmp1[3])
             tmp1[2] + log1p(exp(tmp1[1] - tmp1[2]) + exp(tmp1[3] - tmp1[2]))
           else
             tmp1[3] + log1p(exp(tmp1[1] - tmp1[3]) + exp(tmp1[2] - tmp1[3]))
           end
  end

  lml
end


text:
module TestMercer

using Test
using CovarianceFunctions: iscov
using LinearAlgebra

@testset "promotion" begin
    using CovarianceFunctions: Dot, Poly, NN, ExponentialDot
    k_arr = [Dot(), Poly(3), NN()]
    k_strings = ["Dot", "Poly", "NN"]
    for (k, str) in zip(k_arr, k_strings)
        @test k(1, 2.) ≈ k(1., 2.)
    end
    d = 3
    x, y = randn(d), randn(d)
    @test (Dot()^3)(x, y) ≈ Poly(3)(x, y)
    E = ExponentialDot()
    @test E(x, y) ≈ exp(dot(x, y))
    @test E(x, y) ≈ E(y, x)

end

@testset "mercer" begin
    using CovarianceFunctions: FiniteBasis, gramian, Gramian, LazyMatrixProduct

    basis = [sin, cos, identity]
    k = FiniteBasis(basis)
    @test k isa FiniteBasis
    x = randn(16)
    @test gramian(k, x) isa LazyMatrixProduct
    U = [sin.(x) cos.(x) x]
    @test Matrix(gramian(k, x)) ≈ U*U'

    x = randn(2)
    @test gramian(k, x) isa Gramian # if # of functions is larger than data
    U = [sin.(x) cos.(x) x]
    @test Matrix(gramian(k, x)) ≈ U*U'

    # FiniteBasis([]) # throws
end

# @testset "Mercer Kernels" begin
#     k = [Dot, Poly, NN]
#     @test
# end
#
# using CovarianceFunctions: Product, Sum, Power, VerticalRescaling
# @testset "CovarianceFunctions Algebra" begin
#     k = [Dot, Poly, NN]
#     @test
# end
#
# using RescaledMetric, Periodic, Sym, Conditional, SoR
# @testset "CovarianceFunctions Modification" begin
#     k = [Dot, Poly, NN]
#     @test
# end

# @testset "Matrix Kernels" begin
#     k = [Dot, Poly, NN]
#     @test
# end

#
# # projection onto symmetric matrices: A -> (A' + A)/2
# # testing if the rescaled CovarianceFunctions is numerically stable
# function test(k::K, b::T = T(1), n::Int = 1024) where {T, K<:Rescaled{T}}
#     ε = eps(T) # numerical precision of type
#     x = 2b * rand(n) .- b # randomly sample in range
#     ax = k.a.(x)
#     cond = maximum(ax) / minimum(ax)
#     if cond > sqrt(1/ε)
#         println("Warning: Currently specified Rescaled CovarianceFunctions leads to ill-conditioned GP inference. To avoid this, we added 1e-6 to your rescaling function.")
#         f(x) = k.a(x) + T(1e-6)
#         return K(k.k, f)
#     end
#     return k
# end

end # module TestMercer


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

make_value(::Type{T}, i::Integer) where {T<:Integer} = 3*i%T
make_value(::Type{T},i::Integer) where {T<:AbstractFloat} = T(3*i)

Vec{N,T} = NTuple{N,Base.VecElement{T}}

# Crash report for #15244 motivated this test.
@generated function thrice_iota(::Type{Vec{N,T}}) where {N,T}
    :(tuple($([:(Base.VecElement(make_value($T,$i))) for i in 1:N]...)))
end

function call_iota(n::Integer,t::DataType)
    x = thrice_iota(Vec{n,t})
    @test x[1].value === make_value(t,1)
    @test x[n].value === make_value(t,n)
end

# Try various tuple lengths and element types
for i=1:20
    for t in [Bool, Int8, Int16, Int32, Int64, Float32, Float64]
        call_iota(i,t)
    end
end

# Try various large tuple lengths and element types #20961
for i in (34, 36, 48, 64, 72, 80, 96)
    for t in [Bool, Int8, Int16, Int32, Int64, Float32, Float64]
        call_iota(i,t)
    end
end

# Another crash report for #15244 motivated this test.
struct Bunch{N,T}
    elts::NTuple{N,Base.VecElement{T}}
end

unpeel(x) = x.elts[1].value

@test unpeel(Bunch{2,Float64}((Base.VecElement(5.0),
                               Base.VecElement(4.0)))) === 5.0

rewrap(x) = VecElement(x.elts[1].value+0)
b = Bunch((VecElement(1.0), VecElement(2.0)))

@test rewrap(b)===VecElement(1.0)

struct Herd{N,T}
    elts::NTuple{N,Base.VecElement{T}}
    Herd{N,T}(elts::NTuple{N,T}) where {N,T} = new(ntuple(i->Base.VecElement{T}(elts[i]), N))
end

function check(x::Herd{N,T}) where {N,T}
    for i=1:N
        @test x.elts[i].value === N*N+i-1
    end
end

check(Herd{1,Int}((1,)))
check(Herd{2,Int}((4,5)))
check(Herd{4,Int}((16,17,18,19)))

struct Gr{N, T}
    u::T
    v::Bunch{N,T}
    w::T
end

a = Vector{Gr{2,Float64}}(2)
a[2] = Gr(1.0, Bunch((VecElement(2.0), VecElement(3.0))), 4.0)
a[1] = Gr(5.0, Bunch((VecElement(6.0), VecElement(7.0))), 8.0)
@test a[2] == Gr(1.0, Bunch((VecElement(2.0), VecElement(3.0))), 4.0)

@test isa(VecElement((1,2)), VecElement{Tuple{Int,Int}})

# The following test mimic SIMD.jl
const _llvmtypes = Dict{DataType, String}(
    Float64 => "double",
    Float32 => "float",
    Int32 => "i32",
    Int64 => "i64"
)

@generated function vecadd(x::Vec{N, T}, y::Vec{N, T}) where {N, T}
    llvmT = _llvmtypes[T]
    func = T <: AbstractFloat ? "fadd" : "add"
    exp = """
    %3 = $(func) <$(N) x $(llvmT)> %0, %1
    ret <$(N) x $(llvmT)> %3
    """
    return quote
        Base.@_inline_meta
        Base.llvmcall($exp, Vec{$N, $T}, Tuple{Vec{$N, $T}, Vec{$N, $T}}, x, y)
    end
end

function f20961(x::Vector{Vec{N, T}}, y::Vector{Vec{N, T}}) where{N, T}
    @inbounds begin
        a = x[1]
        b = y[1]
        return vecadd(a, b)
    end
end

# Test various SIMD Vectors with known good sizes
for T in (Float64, Float32, Int64, Int32)
    for N in 1:36
        # For some vectortypes Julia emits llvm arrays instead of vectors
        if N %  7 == 0 || N % 11 == 0 || N % 13 == 0 || N % 15 == 0 ||
           N % 19 == 0 || N % 23 == 0 || N % 25 == 0 || N % 27 == 0 ||
           N % 29 == 0 || N % 31 == 0
            continue
        end
        a = ntuple(i->VecElement(T(i)), N)
        result = ntuple(i-> VecElement(T(i+i)), N)
        b = vecadd(a, a)
        @test b == result
        b = f20961([a], [a])
        @test b == result
    end
end

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


text:
include("sampling.jl")
include("id_generator.jl")
include("span_processor.jl")
include("trace_provider.jl")
include("span.jl")


text:
using Documenter, LeafAreaIndex

makedocs(
    modules = [LeafAreaIndex],
    clean    = true,
    format   = :html,
    sitename = "Leaf Area Index Documentation",
    authors  = "Ken Bastiaensen",
    pages    = [
        "LeafAreaIndex.jl" => "index.md",
        "Quick Intro"      => "quick_intro.md",
        "Gap Fraction"     => "gapfraction.md",
        "LAI Inversion"    => "LAI.md",
        "Clumping"         => "clumping.md",
        "Calibration"      => "calibration.md",
    ],
    linkcheck = false,#!("skiplinks" in ARGS),
    # Use clean URLs, unless built as a "local" build
    html_prettyurls = false#!("local" in ARGS)
)
deploydocs(
    #root   = Pkg.dir("LeafAreaIndex","docs"),
    repo   = "github.com/ETC-UA/LeafAreaIndex.jl.git",
    target = "build",
    julia  = "0.6",
    deps   = nothing,
    make   = nothing
)


text:
function LinearAlgebra.pinv(A::Array{T,2}) where T<:ForwardDiff.Dual
    m = size(A,1)
    n = size(A,2)

    Av = map(ForwardDiff.value, A)
    Ap = map(ForwardDiff.partials, A)
    Apt = Ap[reshape(1:m*n,(m,n))'] # taking the transpose

    Pv = pinv(Av)
    Pp = Pv*Pv'*Apt*(I - Av*Pv) + (I - Pv*Av)*Apt*Pv'*Pv - Pv*Ap*Pv

    P = map(T,Pv,Pp)

    return P
end


text:
# ---
# title: 100. Same Tree
# id: problem100
# author: Qling
# date: 2021-05-23
# difficulty: Easy
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/same-tree/description/>
# hidden: true
# ---
#
# Given two binary trees, write a function to check if they are the same or not.
#
# Two binary trees are considered the same if they are structurally identical
# and the nodes have the same value.
#
# **Example 1:**
#
#
#
#     Input:     1         1
#               / \       / \
#              2   3     2   3
#
#             [1,2,3],   [1,2,3]
#
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input:     1         1
#               /           \
#              2             2
#
#             [1,2],     [1,null,2]
#
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input:     1         1
#               / \       / \
#              2   1     1   2
#
#             [1,2,1],   [1,1,2]
#
#     Output: false
#
#
#
## @lc code=start
using LeetCode

function is_same_tree(p::Union{TreeNode, Nothing}, q::Union{TreeNode, Nothing})::Bool
    if isnothing(p) && isnothing(q)
        return true
    elseif isnothing(p) || isnothing(q)
        return false
    elseif p.val != q.val
        return false
    else
        return is_same_tree(p.left, p.left) && is_same_tree(p.right, p.right)
    end
end

## @lc code=end


text:
struct LinearKF
    xk::Vector
	xp::Vector
	P::AbstractMatrix
	K::AbstractMatrix
	S::AbstractMatrix
	y::Vector
	filter::SquareRootFilter
end

"""
	sol = linear_kalman(
		x, x₀, P, F, H, Q, R;
		B=zeros(length(x), length(x)),
		u=zeros(length(x)),
		w=zeros(length(x)),
		v=zeros(length(x)),
		sqrt_filter=CholeskySqrt(),
		)

Performs the simple Kalman filter algorithm for the estimation of a state x.

x:   Vector with one measurement (observation) of m = length(x) variables
x₀:  A priori estimation of the variables
P:   Initial predicted process covariance matrix
F:	 Transition matrix (linear state model)
H:	 Measurement matrix (linear measurement model)
Q:   Process noise covariance matrix
R:   Measurement (observable) covariance matrix

B:	 Control-input matrix (linear state model)
u:	 Control vector (linear state model)
w:	 Process noise vector, w ~ N(0,Q)
v:	 Measurement (observable) noise vector, v ~ N(0,R)

sqrt_filter: Square root filtering technique for the covariance matrix.
			 Accepted values are: CholeskyModSqrt() (default), SVDSqrt(), CholeskySqrt()

sol: LinearKalman type
    xk:     Kalman filter state estimation
    xp:     Predicted state estimate
    P:      Process covariance matrix a posteriori
    K:      Kalman gain
	S:      Innovation (or pre-fit residual) covariance
	y:      Innovation or measurement pre-fit residual
	filter: Type of SquareRootFilter selected

"""
function linear_kalman(
	x::Vector,
    x₀::Vector,
    P::AbstractMatrix,
	F::AbstractMatrix,
	H::AbstractMatrix,
	Q::AbstractMatrix,
	R::AbstractMatrix;
	B::AbstractMatrix=zeros(length(x), length(x)),
	u::Union{Vector, Real}=zeros(length(x)),
	w::Vector=zeros(length(x)),
	v::Vector=zeros(length(x)),
	sqrt_filter::SquareRootFilter=CholeskyModSqrt(),
	)

	m = length(x)

	## Predict process
	# State estimate propagation (Predicted state, a priori)
	x̂ = _linear_state_model(x₀, F, B, u, w)
	# Error covariance propagation (Predicted process covariance matrix, a priori)
	P_ = _process_covariance_priori(P, F, Q)
	# Innovation (or pre-fit residual) covariance
	PHt = P_ * H'
	S = _innovation_covariance(P_, PHt, H, R)
	# Kalman gain
	K = _kalman_gain(P_, PHt, S)

	## Update process
	# New observation (innovation process)
    y = _linear_observation_model(x, H, v)
	# Innovation or measurement pre-fit residual
	y = muladd(-H, x̂, y)
	# State estimate update (Current state, a posteriori)
	xₖ = muladd(K, y, x̂)
	# Error covariance update (Process covariance matrix, a posteriori)
	P_ = _process_covariance_posteriori(P_, H, K, m)

	return LinearKF(vec(xₖ), vec(x̂), P_, K, S, vec(y), sqrt_filter)
end

function _linear_state_model(
	x::Vector, F::AbstractMatrix, B::AbstractMatrix, u::Union{Vector, Real}, w::Vector,
	)
	return vec(F*x .+ B .* u .+ w)
end

function _linear_observation_model(x::Vector, H::AbstractMatrix, v::Vector)
	return vec(muladd(H, x, v))
end

function _process_covariance_priori(P::AbstractMatrix, F::AbstractMatrix, Q::AbstractMatrix)
	return _return_type(P, F*P*F' .+ Q)
end

function _process_covariance_posteriori(
	P::AbstractMatrix, H::AbstractMatrix, K::AbstractMatrix, m::Int64,
	)
	return _return_type(P, (eye(m) .- K*H) *P)
end

function _kalman_gain(P::AbstractMatrix, PHt::AbstractMatrix, S::AbstractMatrix)
	F = qr(S)
	return _return_type(P, F.R \ (F.Q' * PHt))
end

function _innovation_covariance(
	P::AbstractMatrix, PHt::AbstractMatrix, H::AbstractMatrix, R::AbstractMatrix,
	)
	return _return_type(P, H*PHt .+ R)
end


text:
# This file contains the implementation of various operators.
# Tests for them is at test/runtests.jl.

using Base
using Statistics
# TODO: we need kwarg support for many of these

# Generic
get_tuple(x) = (x...,)
get_tuple() = nothing
convert_type(x) = Base.convert(Array{Float32, 1}, x)

ops[:Concat] = function (params, ip1, ip2)
  s = vcall(:ndims, ip1)

  return vcall(:cat, ip1, ip2, Symbol("dims =1"))
end

ops[:Gemm] = function (params, A, B, C)
  if !haskey(params, :transA)
    params[:transA] = 0
  end
  if !haskey(params, :transB)
    params[:transB] = 0
  end
  if !haskey(params, :alpha)
    params[:alpha] = 1
  end
  if !haskey(params, :beta)
    params[:beta] = 1
  end
  if !haskey(params, :broadcast)
    params[:broadcast] = 0
  end
  if (params[:transA] != 1)
    A =  vcall(:permutedims, A, vcall(:reverse, vcall(:range, 1, vcall(:ndims, A))))
  end
  if (params[:transB] != 1)
    B = vcall(:permutedims, B, vcall(:reverse, vcall(:range, 1, vcall(:ndims, B))))
  end
  ip1 = vcall(:*, params[:alpha], A, B)
  ip2 = vcall(:*, params[:beta], C)
  if params[:broadcast] == 0
    ip1 = vcall(:permutedims, ip1, vcall(:reverse, vcall(:range, 1, vcall(:ndims, ip1))))
    res = vcall(:broadcast, :+, ip1, ip2)
    return res
  end
  res = vcall(:broadcast, :+, ip1, ip2)
  return vcall(:permutedims, res, vcall(:reverse, vcall(:range, 1, vcall(:ndims, res))))
end

# Image

ops[:Conv] = function (params, x, w, b...)
  if !haskey(params, Symbol("pads"))
    params[:pads] = [0,0,0,0]
  end
  if !haskey(params, Symbol("strides"))
    params[:strides] = (1,1)
  end
  if !haskey(params, Symbol("dilations"))
    params[:dilations] = (1,1)
  end
  if (haskey(params, Symbol("auto_pad")))
    if (String(params[:auto_pad]) == "SAME_UPPER" || String(params[:auto_pad] == "SAME_LOWER"))
      temp = Base.convert(Array{Int64,1}, (params[:kernel_shape] .- 1)./2) # Only for strides = [1,1]
      params[:pads] = vcat(temp, temp)                                    # To Do: Add support for other stride values.
    elseif String(params[:auto_pad]) == "VALID"
      params[:pads] = [0,0,0,0]
    end
  end
  #if haskey(params, :group)
  #  s = vcall(:Int, vcall(:/, vcall(:size, x, 3), params[:group]))
  #  x = vcall(:reshape, x, vcall(:size, x, 1), vcall(:size, x, 2), s, params[:group], vcall(:size, x, 4))
  #  temp_x = vcall(:getindex, x, :,:,:,1,:)
  #  temp = vcall(vcall(:Conv, Float32[0], :relu,
  #      Symbol("stride=$((params[:strides]...,))"), Symbol("pad=$(pads(params[:pads]))"),
  #        Symbol("dilation=$((params[:dilations]...,))")), temp_x)
  #  if isempty(b)
  #    for i=2:params[:group]
  #      temp = vcall(:cat, 3, temp, vcall(vcall(:Conv, Float32[0], :relu,
  #        Symbol("stride=$((params[:strides]...,))"), Symbol("pad=$(pads(params[:pads]))"),
  #          Symbol("dilation=$((params[:dilations]...,))")), temp))
  #    end
  #
  #  else
  #    for i=2:params[:group]
  #      temp = vcall(:cat, 3, temp, vcall(vcall(:Conv, b[1], :relu,
  #        Symbol("stride=$((params[:strides]...,))"), Symbol("pad=$(pads(params[:pads]))"),
  #          Symbol("dilation=$((params[:dilations]...,))")), temp))
  #    end
  #  end
  #  return temp
  #end
  if isempty(b)
    return vcall(vcall(:CrossCor, w, Float32[0], :relu, Symbol("stride=$((params[:strides]...,))"),
     Symbol("pad=$((params[:pads]...,))"), Symbol("dilation=$((params[:dilations]...,))")), x)
                                 # temp change (Until type fix)
  end
  vcall(vcall(:CrossCor, w, b[1], Symbol("stride=$((params[:strides]...,))"),
    Symbol("pad=$((params[:pads]...,))"),  Symbol("dilation=$((params[:dilations]...,))")), x)
end

ops[:MaxPool] = function (params, x)
  if !(haskey(params, :strides))
    params[:strides] = [1,1]
  end
  if !(haskey(params, :pads))
    params[:pads] = [0,0,0,0]
  end
  strides = params[:strides] == params[:kernel_shape] ? [] : [params[:strides]]
  if length(params[:kernel_shape]) == 1
    push!(params[:kernel_shape], 1)
    n_size = vcall(:Tuple, vcall(:push!, vcall(:collect, vcall(:size, x)), 1))
    new_x = vcall(:reshape, x, n_size)
    return vcall(:dropdims, vcall(:maxpool, new_x, (params[:kernel_shape]...,), Symbol("pad=$(params[:pads])"),
        Symbol("stride=$((params[:strides]...,))")), Symbol("dims=4"))
  end

  vcall(vcall(:MaxPool, (params[:kernel_shape]...,), Symbol("pad=$((params[:pads]...,))"),Symbol("stride=$((params[:strides]...,))")), x)
end

ops[:GlobalAveragePool] = function (params, x)
  vcall(:mean, x, Symbol("dims = (1,2)"))
end

ops[:GlobalMaxPool] = function (params, x)
  vcall(:getindex, vcall(:findmax, x, Symbol("dims=(1,2)")), 1)
end

ops[:AveragePool] = function (params, x)
  length(params[:kernel_shape]) <= 2 || error("Only averagepool2d currently supported")
  if !haskey(params, :strides)
    params[:strides] = [1,1]
  end
  strides = params[:strides] == params[:kernel_shape] ? [] : [params[:strides]]
  if !haskey(params, :pads)
    params[:pads] = [0,0,0,0]
  end
  if length(params[:kernel_shape]) == 1
    push!(params[:kernel_shape], 1)
    n_size = vcall(:Tuple, vcall(:push!, vcall(:collect, vcall(:size, x)), 1))
    new_x = vcall(:reshape, x, n_size)
    return vcall(:dropdims, vcall(:meanpool, new_x, (params[:kernel_shape]...,), Symbol("pad=$((params[:pads]...,))"),
        Symbol("stride=$((params[:strides]...,))")), Symbol("dims=4"))
  end
  if params[:pads] == [0,0,0,0]
    return vcall(vcall(:MeanPool, (params[:kernel_shape]...,), Symbol("pad=$((params[:pads]...,))"),
                                                    Symbol("stride=$((params[:strides]...,))")), x)
  else
    params[:strides_temp] = [1,1]
    params[:kernel_shape_temp] = [1,1]
    params[:pads_temp] = [0,0,0,0]
    temp = vcall(vcall(:MeanPool, (params[:kernel_shape_temp]...,), Symbol("pad=$((params[:pads]...,))"),
                                                    Symbol("stride=$((params[:strides_temp]...,))")), x)
    return vcall(vcall(:MeanPool, (params[:kernel_shape]...,), Symbol("pad=$((params[:pads_temp]...,))"),
                                                    Symbol("stride=$((params[:strides]...,))")), temp)
  end
end

ops[:BatchNormalization] = function (params, x, scale, b, mean, var)
  if !haskey(params ,Symbol("momentum"))
    params[:momentum] = 0.9
  end
  if !haskey(params, Symbol("epsilon"))
    params[:epsilon] = 1e-5
  end
  t = typeof(params[:momentum])
  q = vcall(:broadcast, :+, params[:epsilon], var)
  p = vcall(:broadcast, sqrt ,q)
  r = vcall(:broadcast, Float32, p)
  return vcall(vcall(:BatchNorm,identity, b, scale, vcall(:broadcast, :Float32, mean), r, t(params[:epsilon]), params[:momentum], false), x)
end

function slice(a, s, e)
  return a[s:e]
end

ops[:LSTM] = function(params, ip...)
  if length(ip) == 3
    len = params[:hidden_size]
    arg1 = vcall(reshape, ip[2], (4*len,2))
    arg2 = vcall(reshape, ip[3], (4*len,3))
    ip_ = vcall(reshape, ip[1], vcall(slice ,vcall(:size, ip[1]), 1, 2))

    a = vcall(LSTM, arg1, arg2, zeros(len*4), zeros(len), zeros(len))

    return vcall(a, ip_)
  elseif length(ip) == 4
    len = params[:hidden_size]
    arg1 = vcall(reshape, ip[2], (4*len,3))
    arg2 = vcall(reshape, ip[3], (4*len,4))
    arg3 = ip[4][1:4*len]
    b1 = vcall(:broadcast, Float32, vcall(reinterpret, Float32, vcall(zeros, 2)))
    a = vcall(LSTM, arg1, arg2, arg3, b1, b1)

    ip_ = vcall(reshape, ip[1], vcall(slice ,vcall(:size, ip[1]), 1, 2))
    return vcall(a, ip_)
  end
end

# Regularise

ops[:Dropout] = function (params, x)
  return vcall(:identity, x)        # Inference mode: Dropout just bypasses input.
end

# Activation

iscallp(f, v) = DataFlow.iscall(v) && f(v[1])
islayer(v, name) = iscallp(l -> iscallp(x -> x == constant(name), l), v)

ops[:Identity] = function(params, x)
  vcall(:identity, x)
end

ops[:Flatten] = function(params, x)
  if !haskey(params, :axis)
    params[:axis] = 1
  end
  l = vcall(:length, x)
  rev = vcall(:reverse, vcall(:size, x))
  if (params[:axis] == 0)
    return vcall(:reshape, x, l, 1)
  else
    s = vcall(:prod, vcall(:getindex, rev, 1:params[:axis]))
    return vcall(:reshape, x, vcall(:div, l, s), s)
  end
end

ops[:Relu] = function (params, x)
  vcall(broadcast, :relu, x)
  #end
end

ops[:LeakyRelu] = function(params, x)
  if !haskey(params, :alpha)
    params[:alpha] = 0.01
  end
  vcall(:broadcast, :leakyrelu, x, params[:alpha])
end

ops[:PRelu] = function(params, x, slope)
  ip1 = vcall(:broadcast, :clamp, x, 0, Inf)
  ip2 = vcall(:.*, vcall(:broadcast, :clamp, x, -Inf, 0), slope)
  return vcall(:broadcast, Float32, vcall(:+, ip1, ip2))
end

ops[:ArgMax] = function(params, x)
  return vcall(Flux.argmax, x)
end

ops[:Abs] = function (params, x)
  vcall(:broadcast, abs, x)
end

ops[:Clip] = function (params, x)
  if !haskey(params, :min)
    params[:min] = vcall(:getindex, vcall(:findmin, x), 1)
  end
  if !haskey(params, :max)
    params[:max] = vcall(:getindex, vcall(:findmax, x), 1)
  end
  vcall(:broadcast, clamp, x, params[:min], params[:max])
end

ops[:Equal] = function(params, x, y)
  return vcall(:broadcast, :Int, vcall(:broadcast, :isequal, x, y))
end

ops[:Greater] = function(params, x, y)
  return vcall(:broadcast, :Int, vcall(:broadcast, :isless, y, x))
end

ops[:Sigmoid] = function (params, x)
  vcall(:sigmoid, x)
end

ops[:Softmax] = function (params, x)
  vcall(:softmax, vcall(:vec, x))
end

ops[:Floor] = function (params, x)
  vcall(:broadcast, :floor, x)
end

ops[:Exp] = function(params, x)
  vcall(:broadcast, :exp, x)
end

ops[:Log] = function(params, x)
  vcall(:broadcast, :log, x)
end

ops[:Neg] = function(params, x)
  vcall(:*, -1,  x)
end

ops[:Sum] = function (params, x, y...)
  if isempty(y)
    return vcall(:.+, x, 0)
  end
  vcall(:+, x, y[1])
end

ops[:Cast] = function(params, x)
  if (params[:to] == 1)
    return vcall(:broadcast, :Float32, x)
  elseif params[:to] == 10
    return vcall(:broadcast, :Float16, x)
  elseif params[:to] == 11
    return vcall(:broadcast, :Float64, x)
  end
end

ops[:Constant] = function (params)
  constant(Symbol("weights[\"$(params.name)\"]"))
end

ops[:Ceil] = function (params ,x)
  vcall(:broadcast, :ceil, x)
end

ops[:Unsqueeze] = function(params, x)
  l1 = length(params[:axes])
  l2 = vcall(:+, l1, vcall(:ndims, x))
  temp = x
  for ele in params[:axes]
    temp = vcall(Flux.unsqueeze, temp, vcall(:-, vcall(:+, vcall(:ndims, temp), 1), ele))
  end
  return temp
end

ops[:Reshape] = function(params, tensor1, shape...)
  if haskey(params, :shape)
    return vcall(:reshape, tensor1, vcall(:broadcast, Int64, vcall(:Tuple, params[:shape])))
  end
  vcall(:reshape, tensor1, vcall(:broadcast, Int64, vcall(:Tuple, vcall(:reverse, shape[1]))))
end

ops[:Transpose] = function(params ,tensor)
  temp_tens = vcall(:permutedims, tensor, vcall(:reverse, vcall(:range, 1, vcall(:ndims, tensor))))
  order = vcall(:.+, params[:perm], 1)
  l = vcall(:permutedims, temp_tens, order)
  return vcall(:permutedims, l, vcall(:reverse, vcall(:range, 1, vcall(:ndims, l))))
end

ops[:LRN] = function(params, x)
  if !haskey(params, :bias)
    params[:bias] = 1
  end
  if !haskey(params, :alpha)
    params[:alpha] = 1e-4
  end
  if !haskey(params, :beta)
    params[:beta] = 0.75
  end
  #return vcall(vcall(:LRNorm, params[:bias], params[:size], params[:alpha], params[:beta]), x)
                               # currently, just bypassing the output
  return vcall(:.+, 0, x)
end

#To-Do : add broadcast here (Urgent)
#         Add axis condition here
ops[:Add] = function(params, A, B)
  s1 = vcall(:size, A)
  s2 = vcall(:size, B)
  if (s1==s2)
    return vcall(:Add, params[:axis], A, B)
  else
    return vcall(:.+, A, B)
  end
end

ops[:Sub] = function(params, A , B)
  s1 = vcall(:size, A)
  s2 = vcall(:size, B)
  if (s1==s2)
    return vcall(:-, A, B)
  else
    return vcall(:.-, A, B)
  end
end

ops[:Div] = function(params, A , B)
  if (haskey(params, :broadcast) && params[:broadcast] == 1)
    if !haskey(params, :axis)
      return vcall(:./, A, B)
    end
    return vcall( :Div, params[:axis], A, B)              # To-Do define Div function
  else
    return vcall(:./  , A, B)   # In case of no broadcast, Perform normal div operation.
  end
end

ops[:Mul] = function (params, A, B)
  if (haskey(params, :broadcast) && params[:broadcast] == 1)
    if !haskey(params, :axis)
      return vcall(:.*, A, B)
    end
    return vcall( :Mul, params[:axis], A, B)              # To-Do define Mul function
  else
    return vcall(:.*, A, B)   # In case of no broadcast, Perform normal Mul operation.
  end
end

ops[:Pow] = function (params, A, B)
  if (haskey(params, :broadcast) && params[:broadcast] == 1)
    if !haskey(params, :axis)
      return vcall(:.^, A, B)
    end
    return vcall( :Pow, params[:axis], A, B)              # To-Do define Pow function
  else
    return vcall(:.^, A, B)   # In case of no broadcast, Perform normal Power operation.
  end
end

ops[:MatMul] = function(params, A, B)
  #tempa = vcall(:permutedims, A, vcall(:reverse, vcall(:range, 1, vcall(:ndims, A))))
  #tempb = vcall(:permutedims, B, vcall(:reverse, vcall(:range, 1, vcall(:ndims, B))))
  vcall(:*, B, A)
end

ops[:size] = function(params, A)
  vcall(:prod, vcall(:size, A))
end

ops[:Sqrt] = function(params, A)
  vcall(:broadcast, :sqrt, A)
end

ops[:Reciprocal] = function(params, A)
  vcall(:./ , 1, A)
end

ops[:Xor] = function (params, A, B)
  ip1 = vcall(:broadcast, &, vcall(:Array, vcall(:broadcast, Bool, A)), vcall(:Array,
              vcall(:broadcast, !, vcall(:broadcast, Bool, B))))
  ip2 = vcall(:broadcast, &, vcall(:Array, vcall(:broadcast, Bool, B)), vcall(:Array,
              vcall(:broadcast, !, vcall(:broadcast, Bool, A))))
  return  vcall(:broadcast, :Int, vcall(:broadcast, |, ip1, ip2))
end

ops[:And] = function(params, A, B)
  if (haskey(params, :broadcast) && params[:broadcast] == 1)
    if !haskey(params, :axis)
      return vcall(:.*, vcall(:broadcast, :Bool, A), vcall(:broadcast, :Bool, B))
    end
    return vcall( :And, params[:axis], A, B)              # To-Do define And function
  else
    return vcall(:.*, vcall(:broadcast, :Bool, A), vcall(:broadcast, :Bool, B))   # In case of no broadcast,
                                                                                #Perform normal And operation.
  end
end

ops[:Or] = function(params, A, B)
  if (haskey(params, :broadcast) && params[:broadcast] == 1)
    if !haskey(params, :axis)
      return vcall(:.+, vcall(:broadcast, :Bool, A), vcall(:broadcast, :Bool, B))
    end
    return vcall( :Or, params[:axis], A, B)              # To-Do define Or function
  else
    return vcall(:.+, vcall(:broadcast, :Bool, A), vcall(:broadcast, :Bool, B))   # In case of no broadcast,
                                                                                #Perform normal Or operation.
  end
end

ops[:Expand] = function(params, A, B)
  shape_new = vcall(:reverse, B)
  return vcall(:repeat , A, Symbol("inner=$(vcall(:reverse, B))"))
end
# Preprocessing

ops[:ImageScaler] = function(params, A)
  if !haskey(params, :scale)
    params[:scale] = 1
  end
  vcall(:.*, A, params[:scale])
end

#Trigonometric ops

ops[:Cos] = function(params, A)
  vcall(:broadcast, :cos, A)
end

ops[:Sin] = function(params, A)
  vcall(:broadcast, :sin, A)
end

ops[:Tan] = function(params, A)
  vcall(:broadcast, :tan, A)
end

ops[:Acos] = function(params, A)
  vcall(:broadcast, :acos, A)
end

ops[:Asin] = function(params, A)
  vcall(:broadcast, :asin, A)
end

ops[:Atan] = function(params, A)
  vcall(:broadcast, :atan, A)
end


text:
module Sample

include("mcmc.jl")

greet() = println("Hello World!")
end # module

max_stars_repo_name: bischtob/SEC

max_issues_repo_name: bischtob/SEC

max_forks_repo_name: bischtob/SEC


text:
### A Pluto.jl notebook ###
# v0.19.0

using Markdown
using InteractiveUtils

# ╔═╡ 915f4df0-bc3a-11ec-2560-c9772e143679
begin
	using Pkg#, Revise
	Pkg.activate(Base.current_project())
	Pkg.instantiate()
	using Catlab.CategoricalAlgebra
	using Catlab.Present, Catlab.Theories
	using AlgebraicPetri
	using AlgebraicPetri: Graph
    using ModelExploration
	using Plots
	using Test
	using Logging
	include("../src/ModelSelection.jl")
	Logging.disable_logging(Logging.Warn)
end;

# ╔═╡ 6b70e364-152e-49c4-ae1b-381cbbcb6fc3
md"""Define typing system for models"""

# ╔═╡ 4dcde238-98a4-4982-9409-cc486a70b1e8
begin
	infectious_type_ = LabelledPetriNet([:Pop],
	  :interaction=>((:Pop, :Pop)=>(:Pop, :Pop)),
	  :t_infection=>(:Pop=>:Pop),
	  :t_strata=>(:Pop=>:Pop))

	infectious_type = map(infectious_type_, Name=name->nothing)

	Graph(infectious_type_)
end

# ╔═╡ e2e3c56d-3350-4ca8-b073-6e6b027ca984
IO_help(i::Int) = let d = Dict([j=>j%2+1 for j in 1:2*i]); (I=d,O=d) end;

# ╔═╡ 179b3378-66ff-43b9-920a-100c48b65a40
make_slice(p::LabelledPetriNet, n::NamedTuple) =
    Slice{ACSetTransformation}(
        homomorphism(p, infectious_type; initial=n,
                       type_components=(Name=x->nothing,)))

# ╔═╡ 712853b5-a679-4c5e-a47b-34af8ee8ea0b
md"""Define petri nets for disease dynamics dimension"""

# ╔═╡ 45511aa5-4a45-4d30-94c7-bc47c7fbfe58
begin
	SIR = LabelledPetriNet([:S, :I, :R],
	  :inf => ((:S, :I)=>(:I, :I)),
	  :rec => (:I=>:R),
	  :id => (:S => :S),
	  :id => (:I => :I),
	  :id => (:R => :R)
	)
	SIR_type = make_slice(SIR, merge((T=[1,2,3,3,3],), IO_help(1)))
	Graph(SIR)
end

# ╔═╡ eedc21a6-5bcd-47e4-8e1b-814c232ab479
begin
	SVIIvR = LabelledPetriNet([:S, :I, :R, :Iv, :V],
	  :inf => ((:S, :I)=>(:I, :I)),
	  :inf => ((:V,:I)=>(:I,:I)),
	  :inf => (:S, :Iv) => (:Iv,:Iv),
	  :inf => ((:V, :Iv) => (:Iv,:Iv)),
	  :rec => (:I=>:R),
	  :rec => (:Iv=>:R),
	  :vax => (:S=>:V),
	  :id => (:S => :S),
	  :id => (:I => :I),
	  :id => (:V => :V),
	  :id => (:Iv => :Iv),
	  :id => (:R => :R)
	)

	SVIIvR_type = make_slice(SVIIvR, merge((T=[1,1,1,1,2,2,2,3,3,3,3,3],),IO_help(4)))
	Graph(SVIIvR)
end

# ╔═╡ 78d94905-418c-4cfc-91c5-9dac9bcd4049
md"""Define petri nets for stratification dimension"""

# ╔═╡ 6f2b8b20-a866-49cd-be09-f659269cbcc3
begin
	quarantine = LabelledPetriNet([:Q, :not_Q],
	    :interaction => ((:not_Q, :not_Q) => (:not_Q, :not_Q)),
	    :enter_quarantine => (:not_Q => :Q),
	    :exit_quarantine => (:Q => :not_Q),
	    :id => (:Q => :Q),
	    :id => (:not_Q => :not_Q),
	)
	quarantine_type = make_slice(quarantine, (T=[1, 3, 3, 2, 2], I=Dict(1=>1,2=>2), O=Dict(1=>1,2=>2)))
	Graph(quarantine)
end

# ╔═╡ 075fa8b1-b45c-43b4-8c24-151cbd3ac7f9
begin
	age_stratification = LabelledPetriNet([:Child, :Adult],
	    :interaction => ((:Child, :Child) => (:Child, :Child)),
	    :interaction => ((:Adult, :Adult) => (:Adult, :Adult)),
	    :interaction => ((:Child, :Adult) => (:Child, :Adult)),
	    :id => (:Child => :Child),
	    :id => (:Adult => :Adult),
	)
	age_s_type = make_slice(age_stratification, merge((T=[1,1,1,2,2],),IO_help(3)))
	Graph(age_stratification)
end

# ╔═╡ 63596b64-b7aa-42ed-a980-99a684c98d4b
begin
	@present ThTwo(FreeCategory) begin
	    (X1,X2)::Ob
	end
	@present ThOne(FreeCategory) begin
	  X::Ob
	end

	const ACSetCat{S} = TypeCat{S, ACSetTransformation}
	const ACSetCatSlice{S} = TypeCat{S, ACSetTransformation}
	const Petri = ACSetCat{LabelledPetriNet}
	const PetriHom = SliceCat{ACSetTransformation}

	One, Two = FinCat.([ThOne,ThTwo])

	to_diag(x) = Diagram(FinDomFunctor(x, nothing, Two, PetriHom()))
	to_slicehom(x) = SliceDiagHom(Literal(to_diag(x)))
end;

# ╔═╡ cc7a7deb-af14-4d15-a444-7a021d765168
md"""Construct the model space as the product of two dimensions"""

# ╔═╡ cea17847-544d-4b9d-89a2-2f2a5337df8b
begin
	diag_disease = to_slicehom(Dict(:X1=>SIR_type, :X2=>SVIIvR_type));
	diag_strata = to_slicehom(Dict(:X1=>quarantine_type, :X2=>age_s_type));
end;

# ╔═╡ c5578104-3cbc-4f3d-b468-dfb702759f89
begin
	# This object contains the recipe to build the composite model space
	pb = PullbackSpace(to_model_hom(diag_disease), to_model_hom(diag_strata));
	# This function evaluates the recipe
	upb = unfold(pb);
	# This discard the projection maps, just keeps the product models
	aupb = apex(upb);
end;

# ╔═╡ 5a4b68ce-0546-43ac-911a-00e638b0b0df
md"""For example, combining SIR with quarantine gives us two copies of SIR (one quarantined, one not) with enter/exit quarantine transitions for each"""

# ╔═╡ 131476ed-ed71-4d5f-9fc0-708f1490c19d
Graph(ob_map(aupb, Symbol("(X1, X1)")))

# ╔═╡ 3dcc4486-0961-44ee-ad90-0c286972fa12
md"""Simulate a trajectory based on the combination of the vaccinated model + quarantine. Arbitrary but realistic parameters, 250 time steps."""

# ╔═╡ e15eb36c-1f27-4394-98e9-77a96b75282b
begin
	true_model = ob_map(aupb, Symbol("(X2, X1)"))
	true_rxn = MakeReactionSystem(true_model)

	p_real = vcat(repeat([1e-4], 14), repeat([0.01], 6))
	p_real = vcat([1.1e-4,1.3e-4,1.2e-4,1.4e-4,1.5e-4,1.2e-4,
				   1.e-4,1.3e-4,1.5e-4,1.2e-4,1.7e-4,1.6e-4,
		           1.2e-4,1.1e-4,], [0.01,0.02,0.01,0.01,0.01,0.02])
	u0 = [0.0,0.0,0.0,0.0,0.0,999.0,1.0,0.0,0.0,0.0]
	tspan = (0.0,250.0)
end;

# ╔═╡ 1f47a2fb-8311-45c0-a178-93ddfe7cbcbf
begin
	sample_data, sample_times, prob_real, sol_real = generate_data(true_model, p_real, u0, tspan, 50)

	plt = plot(sol_real, lw=2, label=reshape(map(string, true_model[:, :sname]), 1, ns(true_model)))
	plot!(sample_times, sample_data, seriestype=:scatter, label="")
end

# ╔═╡ b44ff4bc-92c3-4312-aad3-c1d0878cdf46
md"""Equip models with initial guesses ![](https://i.imgur.com/VKwUP9w.png)
"""

# ╔═╡ 8a0d7fb4-49f5-47be-a890-ea7900ddc58d
begin

disease_initial = (S=999., I=1., R=0.,V=0.,Iv=0.)
strata_initial = (Q=0., not_Q=1.,Child=0.5,Adult=0.5)

models = map(ob_generators(dom(diagram(aupb)))) do o
	model = ob_map(aupb, o)
	init_pop =  [disease_initial[n1] * strata_initial[n2]
				 for (n1, n2) in model[:sname]]
	model => init_pop
end
end;

# ╔═╡ f3b0d954-2eef-4f21-9ba5-fd32a2fbdfb8
function explore(models, tspan, sample_data, sample_times)
  losses = zeros(length(models))
  sols = repeat(Any[nothing], length(models))
  Threads.@threads for i in 1:length(models)
    model, u0 = models[i]
    sol, loss = full_train(model, u0, tspan, sample_data, sample_times)
    losses[i] = loss
	sols[i] = sol
  end
  return losses, sols
end;

# ╔═╡ 5e12700b-d7d9-44c3-9728-71d3c08b18c1
md"""Optimizing all models can take minutes"""

# ╔═╡ 3466eca0-f9d6-4ecf-9eeb-5b9eddc2d061
losses, sols = explore(models, tspan, sample_data, sample_times);

# ╔═╡ 6e81b7d5-f225-42ee-a60d-9076c274a33b
function plot_res(i::Int)
	plot(sample_times, sample_data, seriestype=:scatter, label="")
	        plot!(sols[i], lw=2, label=reshape(map(string, models[i][1][:, :sname]), 1, ns(models[i][1])))
end;

# ╔═╡ da5ab695-e5c3-49b7-b849-afb03a81c598
md"""Model (1,1)"""

# ╔═╡ 4cd3c2f4-5ecb-45f5-b6e0-cbeab132da3a
losses[1]

# ╔═╡ 90ee4cd1-2957-49cf-8f29-fb0c26db3219
plot_res(1)

# ╔═╡ 5d3ae1b2-969a-407c-af4f-452212e4c6f8
md"""Model (1,2)"""

# ╔═╡ 9858e6fc-60a6-417a-8092-280ecd9c139c
losses[3]

# ╔═╡ a9d58251-e0f3-4a77-8b0e-afd4f24c31c4
plot_res(3)

# ╔═╡ b836c19d-b5aa-45d1-8487-586042cbcda2
md"""Model (2,2)"""

# ╔═╡ 65cf7a22-c2ca-4a44-9671-51aab7a0f644
losses[4]

# ╔═╡ 8e5b3b93-d023-4e79-b174-1356b560fd20
plot_res(4)

# ╔═╡ a47b8757-8723-4dfe-a81c-6adf6981f356
md"""Model (2,1): The correct model"""

# ╔═╡ 5999a061-ca85-4002-9b72-5e26bfc97939
losses[2]

# ╔═╡ 3322866c-61bd-40c0-9a59-50610fe37a0d
plot_res(2)

# ╔═╡ Cell order:
# ╠═915f4df0-bc3a-11ec-2560-c9772e143679
# ╟─6b70e364-152e-49c4-ae1b-381cbbcb6fc3
# ╠═4dcde238-98a4-4982-9409-cc486a70b1e8
# ╟─e2e3c56d-3350-4ca8-b073-6e6b027ca984
# ╟─179b3378-66ff-43b9-920a-100c48b65a40
# ╟─712853b5-a679-4c5e-a47b-34af8ee8ea0b
# ╠═45511aa5-4a45-4d30-94c7-bc47c7fbfe58
# ╠═eedc21a6-5bcd-47e4-8e1b-814c232ab479
# ╟─78d94905-418c-4cfc-91c5-9dac9bcd4049
# ╠═6f2b8b20-a866-49cd-be09-f659269cbcc3
# ╠═075fa8b1-b45c-43b4-8c24-151cbd3ac7f9
# ╟─63596b64-b7aa-42ed-a980-99a684c98d4b
# ╟─cc7a7deb-af14-4d15-a444-7a021d765168
# ╠═cea17847-544d-4b9d-89a2-2f2a5337df8b
# ╠═c5578104-3cbc-4f3d-b468-dfb702759f89
# ╟─5a4b68ce-0546-43ac-911a-00e638b0b0df
# ╠═131476ed-ed71-4d5f-9fc0-708f1490c19d
# ╟─3dcc4486-0961-44ee-ad90-0c286972fa12
# ╠═e15eb36c-1f27-4394-98e9-77a96b75282b
# ╠═1f47a2fb-8311-45c0-a178-93ddfe7cbcbf
# ╟─b44ff4bc-92c3-4312-aad3-c1d0878cdf46
# ╠═8a0d7fb4-49f5-47be-a890-ea7900ddc58d
# ╟─f3b0d954-2eef-4f21-9ba5-fd32a2fbdfb8
# ╟─5e12700b-d7d9-44c3-9728-71d3c08b18c1
# ╠═3466eca0-f9d6-4ecf-9eeb-5b9eddc2d061
# ╠═6e81b7d5-f225-42ee-a60d-9076c274a33b
# ╟─da5ab695-e5c3-49b7-b849-afb03a81c598
# ╟─4cd3c2f4-5ecb-45f5-b6e0-cbeab132da3a
# ╠═90ee4cd1-2957-49cf-8f29-fb0c26db3219
# ╟─5d3ae1b2-969a-407c-af4f-452212e4c6f8
# ╟─9858e6fc-60a6-417a-8092-280ecd9c139c
# ╠═a9d58251-e0f3-4a77-8b0e-afd4f24c31c4
# ╟─b836c19d-b5aa-45d1-8487-586042cbcda2
# ╟─65cf7a22-c2ca-4a44-9671-51aab7a0f644
# ╠═8e5b3b93-d023-4e79-b174-1356b560fd20
# ╟─a47b8757-8723-4dfe-a81c-6adf6981f356
# ╟─5999a061-ca85-4002-9b72-5e26bfc97939
# ╠═3322866c-61bd-40c0-9a59-50610fe37a0d