text:
import Tracker
"Flux based Hamiltonian Monte Carlo Sampling"
struct HMCFASTAlg <: SamplingAlgorithm end

"Flux based Hamiltonian Monte Carlo Sampling"
const HMCFAST = HMCFASTAlg()
isapproximate(::HMCFASTAlg) = true
defΩ(::Type{HMCFASTAlg}) = SimpleΩ{ID, Tracker.TrackedArray{Float64, 1, Array{Float64,1}}}
defΩ(::HMCFASTAlg) = SimpleΩ{ID, Tracker.TrackedArray{Float64, 1, Array{Float64,1}}}
defcb(::HMCFASTAlg) = default_cbs()
# defcb = default_cbs(n)

abstract type HMCFASTLoop <: Loop end

"""Hamiltonian monte carlo with leapfrog integration:
https://arxiv.org/pdf/1206.1901.pdf"""
function hmcfast(rng, U, ∇U!, qvals, prop_qvals, pvals, ω, prop_ω, nsteps, stepsize)
  # Initialise proposal as unbounded of current state
  foreach(qvals, prop_qvals) do q, prop_q @. prop_q = (q) end

  # Randomize the momentum
  foreach(p -> @.(p = randn()), pvals)

  # Current Kinetic Energy
  current_K =  sum(map(p->sum(p.^2), pvals)) / 2.0
  ∇qvals = [x.grad for x in values(prop_ω)]
  function ∇step()
    foreach(∇qvals) do ∇q @. ∇q = 0 end  # reset gradients
    ∇U!(prop_ω)  # Gradient step
  end

  # Make a half step for momentum at beginning

  ∇step()
  # Unbound
  foreach(prop_qvals) do prop_q @. prop_q = unbound(prop_q) end
  foreach((p, ∇q, prop_q) -> @.(p = p - stepsize * ∇q * jac(prop_q) / 2.0),
  pvals, ∇qvals, prop_qvals)

  for i = 1:nsteps
    # @show prop_qvals
    # Half step p and q
    # @show prop_qvals
    foreach(pvals, prop_qvals) do p, q @. q = q + stepsize * p end
    # @show prop_qvals
    κ = i !=  nsteps ? 1.0 : 0.5
    # Bound q
    foreach(prop_qvals) do prop_q @. prop_q = bound(prop_q) end

    # Gradient step
    ∇step()
    # Unbound q
    foreach(prop_qvals) do prop_q @. prop_q = unbound(prop_q) end
    foreach(pvals, ∇qvals, prop_qvals) do p, ∇q, q
      @. p = p - κ * stepsize * ∇q * jac(q)
    end
  end

  # Make half a step for momentum at the end
  # any(notunit, q) && return current_q, false
  foreach(prop_qvals) do prop_q @. prop_q = bound(prop_q) end

  # @assert false
  # Evaluate the potential and kinetic energies at start and end
  current_U = U(ω)
  proposed_U = U(prop_ω)
  proposed_K = sum(map(p->sum(p.^2), pvals)) / 2.0

  #@show current_U, proposed_U, current_K, proposed_K
  # Accept or reject
  if log(rand(rng)) < current_U - proposed_U + current_K - proposed_K
    (proposed_U, true)
  else
    (current_U, false)
  end
end

"Sample from `x | y == true` with Hamiltonian Monte Carlo"
function Base.rand(rng::AbstractRNG,
                   ΩT::Type{OT},
                   logdensity::RandVar,
                   n::Integer,
                   alg::HMCFASTAlg;
                   takeevery = 1,
                   nsteps = 10,
                   stepsize = 0.001,
                   ωinit = ΩT(),
                   gradalg = Omega.TrackerGrad,
                   offset = 0) where {OT <: Ω}
  ω = ωinit # Current Ω state of chain
  logdensity(ω)  # Initialize omega
  qvals = [x.data for x in values(ω)]   # Values as a vector

  prop_ω = deepcopy(ω)                          # Ω proposal
  prop_qvals = [x.data for x in values(prop_ω)] # as vector

  p = deepcopy(ω)                     # Momentum, deepcopy but could just zero
  pvals = [x.data for x in values(p)] # as vector

  ωsamples = ΩT[]
  U = -logdensity
  ∇U!(ω) = Omega.back!(U, ω, gradalg)

  accepted = 0
  for i = 1:n*takeevery
    p_, wasaccepted = hmcfast(rng, U, ∇U!, qvals, prop_qvals, pvals, ω,
                          prop_ω, nsteps, stepsize)
    if wasaccepted
      i % takeevery == 0 && push!(ωsamples, deepcopy(prop_ω))
      accepted += 1
      foreach(qvals, prop_qvals) do q, prop_q @. q = prop_q  end
    else
      # QVALS need to reflect
      i % takeevery == 0 && push!(ωsamples, deepcopy(ω))
    end
    lens(HMCFASTLoop, (ω = prop_ω, accepted = accepted, p = Tracker.data(p_), i = i + offset))
  end
  ωsamples
end


text:
using DifferentialEquations
using SimpleDiffEq
using Distributions
using Random
using DataFrames
using StatsPlots
using BenchmarkTools

@inline function rate_to_proportion(r::Float64,t::Float64)
    1-exp(-r*t)
end;

function sir_markov!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ,δt) = p
    N = S+I+R
    ifrac = rate_to_proportion(β*c*I/N,δt)
    rfrac = rate_to_proportion(γ,δt)
    infection=rand(Binomial(S,ifrac))
    recovery=rand(Binomial(I,rfrac))
    @inbounds begin
        du[1] = S-infection
        du[2] = I+infection-recovery
        du[3] = R+recovery
    end
    nothing
end;

δt = 0.1
nsteps = 400
tmax = nsteps*δt
tspan = (0.0,nsteps)
t = 0.0:δt:tmax;

u0 = [990,10,0]; # S,I,R

p = [0.05,10.0,0.25,δt]; # β,c,γ,δt

Random.seed!(1234);

prob_markov = DiscreteProblem(sir_markov!,u0,tspan,p)

sol_markov = solve(prob_markov,solver=FunctionMap());

df_markov = DataFrame(sol_markov')
df_markov[!,:t] = t;

@df df_markov plot(:t,
    [:x1 :x2 :x3],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark solve(prob_markov,solver=FunctionMap)

include(joinpath(@__DIR__,"tutorials","appendix.jl"))
appendix()

ext: jl

max_stars_repo_name: FSUcilab/sir-julia

max_issues_repo_name: FSUcilab/sir-julia

max_forks_repo_name: FSUcilab/sir-julia


text:
# Minimal Generation Gap selection for UNDX
function mgg_alternation!(population::Matrix{Float64}, n_population::Int64, n_children::Int64,
                            n_gene::Int64)::Matrix{Float64}
    idx::Vector{Int64} = randperm(n_population)
    ip::Vector{Int64} = zeros(3)
    ip[1] = idx[1]
    ip[2] = idx[2]

    children::Matrix{Float64} = zeros(n_children,n_gene+1)
    for i in 1:n_children
        ip[3] = rand(idx[3:end])
        children[i,:] = get_new_child(
            population[ip,:], n_gene
        )
    end

    family::Matrix{Float64} = zeros(n_children+2,n_gene+1)
    @inbounds for i in 1:n_gene+1
        @simd for j in 1:n_children
            family[j,i] = children[j,i]
        end
        family[n_children+1,i] = population[ip[1],i]
        family[n_children+2,i] = population[ip[2],i]
    end

    family = sortslices(family, dims=1, by=x->x[end])
    ic2::Int64 = rank_selection(n_children+2)
    @inbounds for i in 1:n_gene+1
        population[ip[1],i] = family[1,i]  # Elite
        population[ip[2],i] = family[ic2,i]  # Rank-based Roulette Selection
    end

    population = sortslices(population, dims=1, by=x->x[end])

    return population
end

function get_new_child(parents::Matrix{Float64}, n_gene::Int64)::Vector{Float64}
    #=
    local child::Vector{Float64}
    MAXITER::Int8 = typemax(Int8)
    in_range::Bool = false
    for _ in 1:MAXITER
        child = UNDX(parents, n_gene)
        # Checking whether children are in search region
        if all(x -> 0.0 <= x <= 1.0, child[1:n_gene])
            in_range = true
            break
        end
    end
    if !(in_range)
        for i in 1:n_gene
            @inbounds child[i] = clamp(child[i], 0.0, 1.0)
        end
    end
    =#
    child::Vector{Float64} = UNDX(parents, n_gene)
    for i in 1:n_gene
        @inbounds child[i] = clamp(child[i], 0.0, 1.0)
    end

    child[end] = objective(
        child[1:n_gene]
    )

    return child
end

# Unimodal Normal Distribution Xover
function UNDX(parents::Matrix{Float64},n_gene::Int64)::Vector{Float64}
    child::Vector{Float64} = zeros(n_gene+1)

    ALPHA::Float64 = 0.5
    BETA::Float64 = 0.35 / sqrt(n_gene)

    p1::Vector{Float64} = [parents[1,i] for i=1:n_gene]
    p2::Vector{Float64} = [parents[2,i] for i=1:n_gene]
    p3::Vector{Float64} = [parents[3,i] for i=1:n_gene]
    d1::Float64 = norm(p2-p1)
    d2::Float64 = norm((p3-p1) - (dot((p3-p1),(p2-p1))/(d1^2)) * (p2-p1))
    e1::Vector{Float64} = p1./d1

    t::Vector{Float64} = randn(n_gene) * BETA * d2
    t .= t - dot(t,e1) * e1
    t .= t + randn() * ALPHA * d1 * e1

    for i in 1:n_gene
        @inbounds child[i] = t[i] + (parents[1,i] + parents[2,i]) / 2.0
    end

    return child
end

function rank_selection(n_family::Int64)::Int64
    ranking::Vector{Int64} = fill(2,n_family)

    for i in 3:n_family
        ranking = append!(ranking,fill(i,n_family-i+2))
    end

    # shuffle!(ranking)
    idx::Int64 = rand(1:length(ranking))

    return ranking[idx]
end


text:
"""
GLOBAL_PAGE_VARS

Dictionary of variables assumed to be set for the entire website. Entries have the format
KEY => PAIR where KEY is a string (e.g.: "author") and PAIR is a pair where the first element is
the default value for the variable and the second is a tuple of accepted possible (super)types
for that value. (e.g.: "THE AUTHOR" => (String, Nothing))

DEVNOTE: marked as constant for perf reasons but can be modified since Dict.
"""
const GLOBAL_PAGE_VARS = PageVars()

"""
$(SIGNATURES)

Convenience function to allocate default values of the global site variables. This is called once,
when JuDoc is started.
"""
@inline function def_GLOBAL_PAGE_VARS!()::Nothing
    empty!(GLOBAL_PAGE_VARS)
    GLOBAL_PAGE_VARS["author"]      = Pair("THE AUTHOR",   (String, Nothing))
    GLOBAL_PAGE_VARS["date_format"] = Pair("U dd, yyyy",   (String,))
    GLOBAL_PAGE_VARS["prepath"]     = Pair("",             (String,))
    return nothing
end

"""
LOCAL_PAGE_VARS

Dictionary of variables copied and then set for each page (through definitions). Entries have the
same format as for `GLOBAL_PAGE_VARS`.

DEVNOTE: marked as constant for perf reasons but can be modified since Dict.
"""
const LOCAL_PAGE_VARS = PageVars()

"""
$(SIGNATURES)

Convenience function to allocate default values of page variables. This is called every time a page
is processed.
"""
@inline function def_LOCAL_PAGE_VARS!()::Nothing
    empty!(LOCAL_PAGE_VARS)
    LOCAL_PAGE_VARS["title"]    = Pair(nothing, (String, Nothing))
    LOCAL_PAGE_VARS["hasmath"]  = Pair(true,    (Bool,))
    LOCAL_PAGE_VARS["hascode"]  = Pair(false,   (Bool,))
    LOCAL_PAGE_VARS["date"]     = Pair(Date(1), (String, Date, Nothing))
    LOCAL_PAGE_VARS["lang"]     = Pair("julia", (String,)) # default lang for indented code
    LOCAL_PAGE_VARS["reflinks"] = Pair(true,    (Bool,))   # whether there are reflinks or not

    # page vars used by judoc, should not be accessed or defined
    LOCAL_PAGE_VARS["jd_ctime"]  = Pair(Date(1), (Date,))   # time of creation
    LOCAL_PAGE_VARS["jd_mtime"]  = Pair(Date(1), (Date,))   # time of last modification
    LOCAL_PAGE_VARS["jd_rpath"]  = Pair("",      (String,)) # local path to file src/[...]/blah.md

    # If there are GLOBAL vars that are defined, they take precedence
    local_keys   = keys(LOCAL_PAGE_VARS)
    for k in keys(GLOBAL_PAGE_VARS)
        k in local_keys || continue
        LOCAL_PAGE_VARS[k] = GLOBAL_PAGE_VARS[k]
    end
    return nothing
end

"""
PAGE_HEADERS

Keep track of seen headers. The key amounts to the ordering (~ordered dict), the value contains
the title, the refstring version of the title, the occurence number and the level (1, ..., 6).
"""
const PAGE_HEADERS = Dict{Int,Tuple{AS,AS,Int,Int}}()

"""
$(SIGNATURES)

Empties `PAGE_HEADERS`.
"""
@inline function def_PAGE_HEADERS!()::Nothing
    empty!(PAGE_HEADERS)
    return nothing
end

"""
PAGE_FNREFS

Keep track of name of seen footnotes; the order is kept as it's a list.
"""
const PAGE_FNREFS = String[]

"""
$(SIGNATURES)

Empties `PAGE_FNREFS`.
"""
@inline function def_PAGE_FNREFS!()::Nothing
    empty!(PAGE_FNREFS)
    return nothing
end

"""
PAGE_LINK_DEFS

Keep track of link def candidates
"""
const PAGE_LINK_DEFS = LittleDict{String,String}()

"""
$(SIGNATURES)

Empties `PAGE_LINK_DEFS`.
"""
@inline function def_PAGE_LINK_DEFS!()::Nothing
    empty!(PAGE_LINK_DEFS)
    return nothing
end

"""
GLOBAL_LXDEFS

List of latex definitions accessible to all pages. This is filled when the config file is read
(via manager/file_utils/process_config).
"""
const GLOBAL_LXDEFS = Dict{String, LxDef}()

"""
EMPTY_SS

Convenience constant for an empty substring, used in LXDEFS.
"""
const EMPTY_SS = SubString("")

"""
$(SIGNATURES)

Convenience function to allocate default values of global latex commands accessible throughout
the site. See [`resolve_lxcom`](@ref).
"""
@inline function def_GLOBAL_LXDEFS!()::Nothing
    empty!(GLOBAL_LXDEFS)
    # hyperreferences
    GLOBAL_LXDEFS["\\eqref"]    = LxDef("\\eqref",    1, EMPTY_SS)
    GLOBAL_LXDEFS["\\cite"]     = LxDef("\\cite",     1, EMPTY_SS)
    GLOBAL_LXDEFS["\\citet"]    = LxDef("\\citet",    1, EMPTY_SS)
    GLOBAL_LXDEFS["\\citep"]    = LxDef("\\citep",    1, EMPTY_SS)
    GLOBAL_LXDEFS["\\label"]    = LxDef("\\label",    1, EMPTY_SS)
    GLOBAL_LXDEFS["\\biblabel"] = LxDef("\\biblabel", 2, EMPTY_SS)
    GLOBAL_LXDEFS["\\toc"]      = LxDef("\\toc",      0, subs("~~~{{toc}}~~~"))
    GLOBAL_LXDEFS["\\tableofcontents"] = LxDef("\\tableofcontents", 0, subs("~~~{{toc}}~~~"))
    # inclusion
    GLOBAL_LXDEFS["\\input"]      = LxDef("\\input",      2, EMPTY_SS)
    GLOBAL_LXDEFS["\\output"]     = LxDef("\\output",     1, EMPTY_SS)
    GLOBAL_LXDEFS["\\textoutput"] = LxDef("\\textoutput", 1, EMPTY_SS)
    GLOBAL_LXDEFS["\\textinput"]  = LxDef("\\textinput", 1, EMPTY_SS)
    GLOBAL_LXDEFS["\\figalt"]     = LxDef("\\figalt",     2, EMPTY_SS)
    GLOBAL_LXDEFS["\\fig"]        = LxDef("\\fig",        1, subs("\\figalt{}{#1}"))
    GLOBAL_LXDEFS["\\file"]       = LxDef("\\file",       2, subs("[#1]()"))
    GLOBAL_LXDEFS["\\tableinput"] = LxDef("\\tableinput", 2, EMPTY_SS)
    # text formatting
    GLOBAL_LXDEFS["\\underline"] = LxDef("\\underline", 1,
                            subs("~~~!#1~~~"))
    GLOBAL_LXDEFS["\\textcss"]   = LxDef("\\underline", 2,
                            subs("~~~!#2~~~"))
    return nothing
end

#= ==========================================
Convenience functions related to the jd_vars
============================================= =#

"""
$(SIGNATURES)

Convenience function taking a `DateTime` object and returning the corresponding formatted string
with the format contained in `GLOBAL_PAGE_VARS["date_format"]`.
"""
jd_date(d::DateTime)::AS = Dates.format(d, GLOBAL_PAGE_VARS["date_format"].first)

"""
$(SIGNATURES)

Checks if a data type `t` is a subtype of a tuple of accepted types `tt`.
"""
check_type(t::DataType, tt::NTuple{N,DataType} where N)::Bool = any(<:(t, tᵢ) for tᵢ ∈ tt)

"""
$(SIGNATURES)

Take a var dictionary `dict` and update the corresponding pair. This should only be used internally
as it does not check the validity of `val`. See [`write_page`](@ref) where it is used to store a
file's creation and last modification time.
"""
set_var!(d::PageVars, k::K, v) where K = (d[k] = Pair(v, d[k].second); nothing)

#= =================================================
set_vars, the key function to assign site variables
==================================================== =#

"""
$(SIGNATURES)

Given a set of definitions `assignments`, update the variables dictionary `jd_vars`. Keys in
`assignments` that do not match keys in `jd_vars` are ignored (a warning message is displayed).
The entries in `assignments` are of the form `KEY => STR` where `KEY` is a string key (e.g.:
"hasmath") and `STR` is an assignment to evaluate (e.g.: "=false").

# Example:

```julia-repl
julia> d = Dict("a"=>(0.5=>(Real,)), "b"=>("hello"=>(String,)));
julia> JuDoc.set_vars!(d, ["a"=>"5.0", "b"=>"\"goodbye\""])
Dict{String,Pair{K,Tuple{DataType}} where K} with 2 entries:
  "b" => "goodbye"=>(String,)
  "a" => 5.0=>(Real,)
```
"""
function set_vars!(jd_vars::PageVars, assignments::Vector{Pair{String,String}})::PageVars
    # if there's no assignment, cut it short
    isempty(assignments) && return jd_vars
    # process each assignment in turn
    for (key, assign) ∈ assignments
        # at this point there may still be a comment in the assignment
        # string e.g. if it came from @def title = "blah"
        # so let's strip


text:
@testset "unit tests" begin

    @testset "tmpdir" begin

        #=
        Leftover files from failed tests or errors can cause testing errors.
        Therefore, we remove any files in the `images` directory not pertinent to testing before executing a test.
        =#

        EXCLUDE_FILES = [".keep"]
        for f in readdir("images")
            !(f in EXCLUDE_FILES) && rm("images/$f")
        end
        @test length(readdir("images")) == length(EXCLUDE_FILES)
    end

    @testset "projection" begin
        x0 = Line(Point(0, 10), O)
        p = Point(10, 20)
        @test projection(p, x0) == Point(0, 20)

        y0 = Line(Point(10, 0), O)
        @test projection(p, y0) == Point(10, 0)
    end

    @testset "translation" begin
        video = Video(500, 500)
        Background(1:100, (args...) -> 1)
        object = Object(1:100, (args...) -> O)
        act!(object, Action(1:100, anim_translate(Point(1, 1), Point(100, 100))))
        Javis.preprocess_frames!(video.objects)

        for f in [1, 50, 100]
            Javis.get_javis_frame(video, [object], f)
            @test get_position(object) == Point(f, f)
        end

        # with easing function
        video = Video(500, 500)
        Background(1:100, (args...) -> 1)
        object = Object(1:100, (args...) -> O)
        act!(object, Action(1:100, sineio(), anim_translate(Point(1, 1), Point(100, 100))))

        action = object.actions[1]

        anim = Animation([0.0, 1.0], [1.0, 100.0], [sineio()])
        Javis.preprocess_frames!(video.objects)

        for f in [1, 50, 100]
            m = (f - 1) / 99
            Javis.get_javis_frame(video, [object], f)
            @test get_position(object) == Point(at(anim, m), at(anim, m))
        end
    end

    @testset "Relative frames" begin
        video = Video(500, 500)
        Background(1:110, (args...) -> 1)
        Object(1:100, (args...) -> 1)
        # dummy object doesn't need a real function
        object = Object(RFrames(10), (args...) -> 1)
        test_file = render(video)
        @test Javis.get_frames(object) == 101:110
        rm(test_file)
    end

    @testset "Global frames" begin
        video = Video(500, 500)
        Background(1:110, (args...) -> 1)
        Object(1:100, (args...) -> 1)
        # dummy object doesn't need a real function
        object = Object(RFrames(10), (args...) -> 1)
        # defined globally but will be computed to local time frame -> 5:10
        act!(object, Action(GFrames(105:110), (args...) -> 1))
        render(video; pathname = "")
        @test Javis.get_frames(object.actions[1]) == 5:10
    end

    @testset "translation from origin" begin
        video = Video(500, 500)
        Background(1:100, (args...) -> 1)
        object = Object(1:100, (args...) -> O)
        act!(object, Action(1:100, anim_translate(Point(99, 99))))

        Javis.preprocess_frames!(video.objects)

        for f in [1, 50, 100]
            Javis.get_javis_frame(video, [object], f)
            @test get_position(object) == Point(f - 1, f - 1)
        end

        video = Video(500, 500)
        Background(1:100, (args...) -> 1)
        object = Object(1:100, (args...) -> O)
        act!(object, Action(1:100, anim_translate(99, 99)))

        Javis.preprocess_frames!(video.objects)

        for f in [1, 50, 100]
            Javis.get_javis_frame(video, [object], f)
            @test get_position(object) == Point(f - 1, f - 1)
        end
    end

    @testset "rotations" begin
        video = Video(500, 500)
        Background(1:100, (args...) -> 1)
        object = Object(1:100, (args...) -> O)
        act!(object, Action(1:100, anim_rotate(2π)))
        anim = Animation([0.0, 1.0], [0.0, 2π], [linear()])
        Javis.preprocess_frames!(video.objects)

        for f in [1, 50, 100]
            Javis.get_javis_frame(video, [object], f)
            @test get_angle(object) ≈ at(anim, (f - 1) / 99)
        end

        video = Video(500, 500)
        Background(1:100, (args...) -> 1)
        object = Object(1:100, (args...) -> O)
        rp = Point(2.0, 5.0)
        act!(object, Action(1:100, anim_rotate_around(2π, rp)))
        Javis.preprocess_frames!(video.objects)

        # compute radius
        r = sqrt(rp.x^2 + rp.y^2)
        shifted_start = -rp
        shifted_X = Point(r, 0)
        # compute start angle of O and rotation point
        start_angle = acos(dotproduct(shifted_X, shifted_start) / r^2)

        for f in [1, 50, 100]
            Javis.get_javis_frame(video, [object], f)
            @test get_angle(object) ≈ at(anim, (f - 1) / 99)
            # rotate clockwise because Luxor is flipped
            p = polar(r, start_angle - get_angle(object))
            p1 = Point(p.x, -p.y)
            @test get_position(object) ≈ p1 + rp
        end
    end

    @testset "scaling" begin
        video = Video(500, 500)
        Background(1:100, (args...) -> 1)
        # dummy object doesn't need a real function
        object = Object(1:50, (args...) -> O)
        act!(object, Action(1:50, anim_scale(1.0, 0.5)))
        Javis.preprocess_frames!(video.objects)

        anim = Animation([0.0, 1.0], [1.0, 0.5], [linear()])
        for f in [1, 50]
            Javis.get_javis_frame(video, [object], f)
            @test get_scale(object).x ≈ at(anim, (f - 1) / 49)
            @test get_scale(object).y ≈ at(anim, (f - 1) / 49)
        end
    end

    @testset "Relative frames" begin
        video = Video(500, 500)
        Background(1:110, (args...) -> 1)
        o1 = Object(1:100, (args...) -> 1)

        o2 = Object(RFrames(10), (args...) -> 1)
        test_file = render(video)

        @test Javis.get_frames(o2) == 101:110
        rm(test_file)
    end

    @testset "translation from origin" begin
        video = Video(500, 500)
        Background(1:100, (args...) -> 1)
        object = Object(1:100, (args...) -> O)
        act!(object, Action(1:100, anim_translate(Point(99, 99))))
        Javis.preprocess_frames!(video.objects)

        for f in [1, 50, 100]
            Javis.get_javis_frame(video, [object], f)
            @test get_position(object) == Point(f - 1, f - 1)
        end

        video = Video(500, 500)
        Background(1:100, (args...) -> 1)
        object = Object(1:100, (args...) -> O)
        act!(object, Action(1:100, anim_translate(99, 99)))
        Javis.preprocess_frames!(video.objects)

        for f in [1, 50, 100]
            Javis.get_javis_frame(video, [object], f)
            @test get_position(object) == Point(f - 1, f - 1)
        end
    end

    @testset "action with a single frame" begin
        video = Video(500, 500)
        Background(1:100, (args...) -> 1)
        # dummy object doesn't need a real function
        object = Object(1:100, (args...) -> O)
        act!(object, Action(1:1, anim_translate(Point(10, 10))))
        Javis.preprocess_frames!(video.objects)

        Javis.get_javis_frame(video, [object], 1)
        @test get_position(object) == Point(10, 10)
    end

    @testset "Frames errors" begin
        # throws because a Video object was not previously defined
        empty!(Javis.CURRENT_VIDEO)
        @test_throws ErrorException Object(1:10, (args...) -> star(O, 20, 5, 0.5, 0, :fill))
        # throws because the frames of the first object need to be defined explicitly
        video = Video(500, 500)
        Object((args...) -> 1)
        @test_throws ArgumentError render(video)

        video = Video(500, 500)
        Background(RFrames(10), (args...) -> 1)
        # throws because the frames of the first object need to be defined explicitly
        @test_throws ArgumentError render(video)

        video = Video(500, 500)
        Background(1:100, (args...) -> 1)
        Object(1:100, (args...) -> 1)
        Object(:some, (args...) -> 1)
        # throws because :some is not supported as Symbol for `frames`
        @test_throws ArgumentError render(video)
    end

    @testset "Frame computation" begin
        demo = Video(500, 500)
        back = Background(1:50, (args...) -> 1)
        obj = Object((args...) -> 1)
        act!(obj, Action(anim_scale(1, 2)))

        objects = [back, obj]

        Javis.preprocess_frames!(demo.objects)
        @test Javis.get_frames(objects[1]) == 1:50
        @test Javis.get_frames(objects[2]) == 1:50
        @test Javis.get_frames(objects[2].actions[1]) == 1:50

        demo = Video(500, 500)
        back = Background(1:50, (args...) -> 1)
        obj = Object(RFrames(-19:0), (args...) -> 1)
        act!(obj, Action(1:10, anim_scale(1, 2)))
        act!(obj, Action(RFrames(10), anim_scale(1, 2)))

        objects = [back, obj]
        Javis.preprocess_frames!(demo.objects)
        @test Javis.get_frames(objects[1]) == 1:50
        @test Javis.get_frames(objects[2]) == 31:50
        @test Javis.get_frames(objects[2].actions[1]) == 1:10
        @test Javis.get_frames(objects[2].actions[2]) == 11:20
    end

    @testset "anim_" begin
        s = anim_scale(1, 2)
        @test s.from == Javis.Scale(1, 1)
        @test s.to == Javis.Scale(2, 2)

        s = anim_scale(Javis.Scale(1, 1), 2)
        @test s.from == Javis.Scale(1, 1)
        @test s.to == Javis.Scale(2, 2)

        s = anim_scale(Javis.Scale(1, 1), Javis.Scale(2, 2))
        @test s.from == Javis.Scale(1, 1)
        @test s.to == Javis.Scale(2, 2)

        s = anim_scale((1, 1), Javis.Scale(2, 1))
        @test s.from == Javis.Scale(1, 1)
        @test s.to == Javis.Scale(2, 1)

        s = anim_scale(Javis.Scale(2, 1))
        @test s.from == :current_scale
        @test s.to == Javis.Scale(2, 1)
    end

    @testset "Test warning if background not defined" begin
        video = Video(100, 100)
        Background(1:10, (args...) -> 1)
        Object(1:11, (args...) -> 2)
        @test_logs (:warn,) render(video; pathname = "")
    end

    @testset "Test warning if action outside frame range" begin
        video = Video(100, 100)
        Background(1:20, (args...) -> 1)
        act!(Object(1:11, (args...) -> 2), Action(1:20, anim_translate(0, 10)))
        @test_logs (:warn,) render(video; pathname = "")
    end
end


text:
"""
Module with general utilities and missing standard library features that could be useful in any Julia project
"""
module Utils

export to_long_mi,
    generate_all, generate_data_all

###################
# Misc.
####################

function to_long_mi(m::Matrix{Float64}, min_int, max_int)::Matrix{Int64}
    δmi = maximum(m) - minimum(m)
    δint = max_int - min_int
    return @. round(Int64, m * δint / δmi + min_int)
end

###################
# One-Hot Encoding
####################
"""
One-hot encode data (2-D Array) based on categories (1-D Array)
Each row of the return value is a concatenation of one-hot encoding of elements of the same row in data
Assumption: both input arrays have elements of same type
"""
function one_hot_encode(X::Array{T, 2}, categories::Array{T,1}) where {T<:Any}
    X_dash = zeros(Bool, size(X)[1], length(categories)*size(X)[2])
    for i = 1:size(X)[1], j = 1:size(X)[2]
            X_dash[i, (j-1)*length(categories) + findfirst(==(X[i,j]), categories)] = 1
    end
    X_dash
end

###################
# Testing Utils
####################

"""
Given some missing values generates all possible fillings
"""
function generate_all(row::Array{Int8})
    miss_count = count(row .== -1)
    lits = length(row)
    result = Bool.(zeros(1 << miss_count, lits))

    if miss_count == 0
        result[1, :] = copy(row)
    else
        for mask = 0: (1<<miss_count) - 1
            cur = copy(row)
            cur[row .== -1] = transpose(parse.(Bool, split(bitstring(mask)[end-miss_count+1:end], "")))
            result[mask+1,:] = cur
        end
    end
    result
end

"""
Generates all possible binary configurations of size N
"""
function generate_data_all(N::Int)
    data_all = transpose(parse.(Bool, split(bitstring(0)[end-N+1:end], "")));
    for mask = 1: (1<<N) - 1
        data_all = vcat(data_all,
            transpose(parse.(Bool, split(bitstring(mask)[end-N+1:end], "")))
        );
    end
    data_all
end

end #module


text:
module NNlib
using Requires, TimerOutputs

# Include APIs
include("dim_helpers.jl")
include("activation.jl")
include("softmax.jl")
include("gemm.jl")
include("conv.jl")
include("pooling.jl")

## Include implementations
include("impl/padding_edges.jl")

# Direct implementations of convolutional and depthwise-convolutional algorithms
include("impl/conv_direct.jl")
include("impl/depthwiseconv_direct.jl")
# im2col implementations of convolutional and depthwise-convolutional algorithms
include("impl/conv_im2col.jl")
include("impl/depthwiseconv_im2col.jl")

# Direct implementations of pooling
include("impl/pooling_direct.jl")

to = TimerOutput()

end # module NNlib


text:
import Reactive
const re = Reactive

"""
Overwrites the `fpswhen_connect` function from Reactive defined in time.jl.
This allows to use a singnal object as rate instead of a single value.
"""
function re.fpswhen_connect(rate::re.Signal, switch, output, name)
    prev_time = time()
    outputref = WeakRef(output)
    switchref = WeakRef(switch)
    timer = Timer(identity, 0) # dummy timer to initialise
    function fpswhen_runner()
        # this function will run if switch gets a new value (i.e. is "active")
        # and if output is pushed to (assumed to be by the timer)
        if switch.value
            start_time = time()
            timer = re.setup_next_tick(outputref, switchref, start_time-prev_time, 1.0/re.value(rate))
            prev_time = start_time
        else
            re.close(timer)
        end
        switch.value
    end
    # the fpswhen_aux will start and stop the timer if switch's value updates, it'll
    # also setup the next tick once the first tick is pushed to output
    fpswhen_aux = re.Signal(switch.value, (switch, output); name="fpswhen runner switch: $(switch.name), output: $(output.name)")
    re.preserve(fpswhen_aux)
    re.add_action!(fpswhen_runner, fpswhen_aux)

    fpswhen_runner() # init
    # ensure timer stops when output node is garbage collected
    finalizer(_->re.close(timer), output)
end


text:
export learn_bdd, generate_from_bdd

"Returns an increasing sorted right-linear vtree."
function bdd_vtree!(V::AbstractVector{<:Integer})::Vtree
    passdown(x::Integer)::Vtree = PlainVtreeLeafNode(x)
    passdown(X::AbstractVector{<:Integer})::Vtree = Vtree(passdown(popfirst!(X)),
                                                          passdown(length(X) == 1 ? X[1] : X))
    return passdown(sort!(V))
end

"Returns a fully factorized circuit reusing leaf nodes from `L`."
function factorize_reuse(X::Vector{Int}, L::Dict{Int32, StructProbLiteralNode})::StructProbCircuit
    U = Vtree([PlainVtreeLeafNode(x) for x ∈ X], :random)
    function passdown(V::Vtree)::StructProbCircuit
        if isleaf(V)
            v = convert(Int32, V.var)
            return StructSumNode([get_lit(-v, V, L), get_lit(v, V, L)], V)
        end
        left = passdown(V.left)
        right = passdown(V.right)
        return StructSumNode([StructMulNode(left, right, V)], V)
    end
    return passdown(U)
end

"Returns a structured probabilistic circuit compiled from a binary decision diagram."
function generate_from_bdd(ϕ::Bdd, n::Integer)::StructProbCircuit
    Sc = scope(ϕ)
    X = setdiff!(collect(1:n), Sc)
    U = bdd_vtree!(Sc)
    L = Dict{Int32, StructProbLiteralNode}()
    visited = Dict{UInt64, StructSumNode}()
    function passdown(V::Vtree, α::Bdd)::StructProbCircuit
        sh = hash((shallowhash(α), V))
        if haskey(visited, sh) return visited[sh] end
        if is_⊤(α)
            if isleaf(V)
                v = convert(Int32, V.var)
                return StructSumNode([get_lit(-v, V, L), get_lit(v, V, L)], V)
            end
            v = convert(Int32, V.left.var)
            sub = passdown(V.right, α)
            return StructSumNode([StructMulNode(get_lit(-v, V.left, L), sub, V),
                                  StructMulNode(get_lit(v, V.left, L), sub, V)], V)
        end
        v = convert(Int32, α.index)
        if isleaf(V) && is_lit(α) return get_lit(to_lit(α), V, L) end
        C = StructProbCircuit[]
        if V.left.var != v
            sub = passdown(V.right, α)
            v = convert(Int32, V.left.var)
            push!(C, StructMulNode(get_lit(-v, V.left, L), sub, V))
            push!(C, StructMulNode(get_lit(v, V.left, L), sub, V))
        else
            if !is_⊥(α.low) push!(C, StructMulNode(get_lit(-v, V.left, L), passdown(V.right, α.low), V)) end
            if !is_⊥(α.high) push!(C, StructMulNode(get_lit(v, V.left, L), passdown(V.right, α.high), V)) end
        end
        s = StructSumNode(C, V)
        visited[sh] = s
        return s
    end
    f = passdown(U, ϕ)
    if isempty(X) return f end
    p = factorize_reuse(X, L)
    if isempty(Sc) return p end
    V = Vtree(U, p.vtree)
    pc = StructSumNode([StructMulNode(f, p, V)], V)
    return pc
end

"Learns a structured probabilistic circuit consistent with a binary decision diagram `ϕ`."
function learn_bdd(ϕ::Bdd, D::DataFrame; pseudocount::Real)::StructProbCircuit
    pc = generate_from_bdd(ϕ, ncol(D))
    estimate_parameters(pc, D; pseudocount)
    return pc
end


text:
export MyLinOp

"""
`MyLinOp(domainType::Type, dim_in::Tuple, [domainType::Type,] dim_out::Tuple, Fwd!::Function, Adj!::Function)`

Construct a user defined `LinearOperator` by specifing its linear mapping `Fwd!` and its adjoint `Adj!`. The functions `Fwd!` and `Adj` must be in-place functions consistent with the given dimensions `dim_in` and `dim_out` and the domain and codomain types.

```julia
julia> n,m = 5,4;

julia> A = randn(n,m);

julia> op = MyLinOp(Float64, (m,),(n,), (y,x) -> mul!(y,A,x), (y,x) -> mul!(y,A',x))
A  ℝ^4 -> ℝ^5

julia> op = MyLinOp(Float64, (m,), Float64, (n,), (y,x) -> mul!(y,A,x), (y,x) -> mul!(y,A',x))
A  ℝ^4 -> ℝ^5

```

"""
struct MyLinOp{N,M,C,D} <: LinearOperator
	dim_out::NTuple{N,Int}
	dim_in::NTuple{M,Int}
	Fwd!::Function
	Adj!::Function
end

# Constructors

MyLinOp(domainType::Type, dim_in::NTuple{N,Int}, dim_out::NTuple{M,Int},
	Fwd!::Function, Adj!::Function ) where {N,M} =
MyLinOp{N,M, domainType, domainType}(dim_out, dim_in, Fwd!, Adj! )

MyLinOp(domainType::Type, dim_in::NTuple{N,Int}, codomainType::Type, dim_out::NTuple{M,Int},
	Fwd!::Function, Adj!::Function ) where {N,M} =
MyLinOp{N,M, domainType, codomainType}(dim_out, dim_in, Fwd!, Adj! )

# Mappings

mul!(y::Array{C,N}, L::MyLinOp{N,M,C,D}, b::Array{D,M}) where {N,M,C,D} = L.Fwd!(y,b)
mul!(y::Array{C,N}, L::AdjointOperator{MyLinOp{N,M,C,D}}, b::Array{D,M}) where {N,M,C,D} = L.A.Adj!(y,b)

# Properties

size(L::MyLinOp) = (L.dim_out, L.dim_in)

codomainType(L::MyLinOp{N,M,C,D}) where {N,M,C,D} = C
  domainType(L::MyLinOp{N,M,C,D}) where {N,M,C,D} = D

fun_name(L::MyLinOp)  = "A"


text:
using Debugger

const boundary = [0 100; 0 100; 0 100]
const initTemp = 300.0
const totalTime = 1000
const outputInterval = 1
const cellCutoff = 10.0
const boundaryCondition = "open"

include("Main.jl")

function CustomEvents(universe::Universe)
    if universe.thermo.iterTime == 0
        universe.cache["outputTime"] = 0
        println("Start")
        event = IntroduceRandomDefects(1000,1)
        DefineUniverseEvent(universe, event, 1.0)
        #Introduce!(universe, [IntroduceRandomDefects(1000,1)])
    end

    #if universe.thermo.time > universe.cache["outputTime"]
    if universe.thermo.iterTime % 1000 == 0
        println(universe.thermo.iterTime, " ",
                universe.thermo.time, " ",
                length(universe)," ",
                universe.fre, " ")
        Dump(universe, "run.dump")
        universe.cache["outputTime"] += outputInterval
    end

end

function StopEvents(universe::Universe)
    if universe.thermo.time  >=  totalTime
    #if universe.thermo.iterTime > 100000
        #println(universe.thermo.iterTime)
        Dump(universe, "run.dump")
        universe.stop = true
    end
end

#global debug= Vector()
Random.seed!(20191201)
using BenchmarkTools
Evolve()
#pprof(;webhost="115.25.142.8")

max_stars_repo_name: permissionx/mc-ke

max_issues_repo_name: permissionx/mc-ke

max_forks_repo_name: permissionx/mc-ke


text:
## Graph subtyping

struct GridGraph <: AbstractGraph{Int}
    grid::Matrix{Char}
end

is_passable(c::Char) = (c == '.') || (c == 'G') || (c == 'S')

function node_index(g::GridGraph, i::Integer, j::Integer)
    n, m = size(g)
    if (1 <= i <= n) && (1 <= j <= m)
        v = (i - 1) * m + (j - 1) + 1  # enumerate row by row
        return v
    else
        return 0
    end
end

function node_coord(g::GridGraph, v::Integer)
    n, m = size(g)
    if 1 <= v <= n * m
        i = (v - 1) ÷ m + 1
        j = (v - 1) % m + 1
        return i, j
    else
        return (0, 0)
    end
end

Base.eltype(::GridGraph) = Int
Graphs.edgetype(::GridGraph) = Graphs.SimpleEdge{Int}

Graphs.is_directed(::GridGraph) = true
Graphs.is_directed(::Type{<:GridGraph}) = true

Base.size(g::GridGraph, args...) = size(g.grid, args...)

Graphs.nv(g::GridGraph) = prod(size(g))
Graphs.vertices(g::GridGraph) = 1:nv(g)
Graphs.has_vertex(g::GridGraph, v::Integer) = 1 <= v <= nv(g)

function Graphs.ne(g::GridGraph)
    m = 0
    for s in vertices(g)
        for d in outneighbors(g, s)
            if s < d
                m += 1
            end
        end
    end
    return m
end

function Graphs.has_edge(g::GridGraph, s::Integer, d::Integer)
    if has_vertex(g, s) && has_vertex(g, d) && (d != s)
        is, js = node_coord(g, s)
        id, jd = node_coord(g, d)
        if !is_passable(g.grid[is, js]) || !is_passable(g.grid[id, jd])
            return false
        elseif (abs(is - id) > 1) || (abs(js - jd) > 1)
            return false
        elseif (is == id) || (js == jd)
            return true
        else
            return is_passable(g.grid[is, jd]) || is_passable(g.grid[id, js])
        end
    else
        return false
    end
end

function Graphs.outneighbors(g::GridGraph, s::Integer)
    n, m = size(g)
    i, j = node_coord(g, s)
    possible_neighbors = (
        (i - 1, j - 1),
        (i - 1, j),
        (i - 1, j + 1),
        (i, j - 1),
        (i, j + 1),
        (i + 1, j - 1),
        (i + 1, j),
        (i + 1, j + 1),
    )  # listed in ascending index order!
    neighbors = (
        node_index(g, id, jd) for (id, jd) in possible_neighbors if
        (1 <= id <= n) && (1 <= jd <= m) && has_edge(g, s, node_index(g, id, jd))
    )
    return neighbors
end

Graphs.inneighbors(g::GridGraph, d::Integer) = outneighbors(g, d)

function Graphs.edges(g::GridGraph)
    return (Graphs.SimpleEdge(s, d) for s in vertices(g) for d in outneighbors(g, s))
end

function Graphs.weights(g::GridGraph)
    E = edges(g)
    I = [src(ed) for ed in E]
    J = [dst(ed) for ed in E]
    V = Float64[]
    for ed in E
        s, d = src(ed), dst(ed)
        is, js = node_coord(g, s)
        id, jd = node_coord(g, d)
        if (is == id) || (js == jd)
            push!(V, 1.)
        else
            push!(V, √2)
        end
    end
    w = sparse(I, J, V, nv(g), nv(g))
    return w
end

## Shortest path functions

function shortest_path_grid(g::GridGraph, (is, js), (id, jd))
    # Compute the shortest path
    w = weights(g)
    s = node_index(g, is, js)
    d = node_index(g, id, jd)
    heuristic(v) = sum(abs, (id, jd) .- node_coord(g, v))
    path = a_star(g, s, d, w, heuristic)
    return path
end

Graphs.reverse(g::GridGraph) = g


text:
"""
A rewriter is any function which takes an expression and returns an expression
or `nothing`. If `nothing` is returned that means there was no changes applicable
to the input expression.

The `SymbolicUtils.Rewriters` module contains some types which create and transform
rewriters.

- `Empty()` is a rewriter which always returns `nothing`
- `Chain(itr)` chain an iterator of rewriters into a single rewriter which applies
   each chained rewriter in the given order.
   If a rewriter returns `nothing` this is treated as a no-change.
- `RestartedChain(itr)` like `Chain(itr)` but restarts from the first rewriter once on the
   first successful application of one of the chained rewriters.
- `IfElse(cond, rw1, rw2)` runs the `cond` function on the input, applies `rw1` if cond
   returns true, `rw2` if it retuns false
- `If(cond, rw)` is the same as `IfElse(cond, rw, Empty())`
- `Prewalk(rw; threaded=false, thread_cutoff=100)` returns a rewriter which does a pre-order
   traversal of a given expression and applies the rewriter `rw`. `threaded=true` will
   use multi threading for traversal. `thread_cutoff` is the minimum number of nodes
   in a subtree which should be walked in a threaded spawn.
- `Postwalk(rw; threaded=false, thread_cutoff=100)` similarly does post-order traversal.
- `Fixpoint(rw)` returns a rewriter which applies `rw` repeatedly until there are no changes to be made.
- `PassThrough(rw)` returns a rewriter which if `rw(x)` returns `nothing` will instead
   return `x` otherwise will return `rw(x)`.

"""
module Rewriters
using SymbolicUtils: @timer, is_operation, istree, operation, similarterm, arguments, node_count

export Empty, IfElse, If, Chain, RestartedChain, Fixpoint, Postwalk, Prewalk, PassThrough

# Cache of printed rules to speed up @timer
const repr_cache = IdDict()
cached_repr(x) = Base.get!(()->repr(x), repr_cache, x)

struct Empty end

(rw::Empty)(x) = nothing

struct IfElse{F, A, B}
    cond::F
    yes::A
    no::B
end

function (rw::IfElse)(x)
    rw.cond(x) ?  rw.yes(x) : rw.no(x)
end

If(f, x) = IfElse(f, x, Empty())

struct Chain
    rws
end

function (rw::Chain)(x)
    for f in rw.rws
        y = @timer cached_repr(f) f(x)
        if y !== nothing
            x = y
        end
    end
    return x
end

struct RestartedChain{Cs}
    rws::Cs
end

function (rw::RestartedChain)(x)
    for f in rw.rws
        y = @timer cached_repr(f) f(x)
        if y !== nothing
            return Chain(rw.rws)(y)
        end
    end
    return x
end

@generated function (rw::RestartedChain{<:NTuple{N,Any}})(x) where N
    quote
        Base.@nexprs $N i->begin
            let f = rw.rws[i]
                y = @timer cached_repr(repr(f)) f(x)
                if y !== nothing
                    return Chain(rw.rws)(y)
                end
            end
        end
        return x
    end
end
struct Fixpoint{C}
    rw::C
end

function (rw::Fixpoint)(x)
    f = rw.rw
    y = @timer cached_repr(f) f(x)
    while x !== y && !isequal(x, y)
        isnothing(y) && return x
        x = y
        y = @timer cached_repr(f) f(x)
    end
    return x
end

struct Walk{ord, C, threaded}
    rw::C
    thread_cutoff::Int
end

using .Threads

function Postwalk(rw; threaded::Bool=false, thread_cutoff=100)
    Walk{:post, typeof(rw), threaded}(rw, thread_cutoff)
end

function Prewalk(rw; threaded::Bool=false, thread_cutoff=100)
    Walk{:pre, typeof(rw), threaded}(rw, thread_cutoff)
end

struct PassThrough{C}
    rw::C
end
(p::PassThrough)(x) = (y=p.rw(x); isnothing(y) ? x : y)

passthrough(x, default) = isnothing(x) ? default : x
function (p::Walk{ord, C, false})(x) where {ord, C}
    @assert ord === :pre || ord === :post
    if istree(x)
        if ord === :pre
            x = p.rw(x)
        end
        if istree(x)
            x = similarterm(x, operation(x), map(PassThrough(p), arguments(x)))
        end
        return ord === :post ? p.rw(x) : x
    else
        return p.rw(x)
    end
end

function (p::Walk{ord, C, true})(x) where {ord, C}
    @assert ord === :pre || ord === :post
    if istree(x)
        if ord === :pre
            x = p.rw(x)
        end
        if istree(x)
            _args = map(arguments(x)) do arg
                if node_count(arg) > p.thread_cutoff
                    Threads.@spawn p(arg)
                else
                    p(arg)
                end
            end
            args = map((t,a) -> passthrough(t isa Task ? fetch(t) : t, a), _args, arguments(x))
            t = similarterm(x, operation(x), args)
        end
        return ord === :post ? p.rw(t) : t
    else
        return p.rw(x)
    end
end

end # end module


text:
mutable struct VarInfo5
    name_user::String
    unit::String
    yeardict::Dict{String, Tuple{String, String, Vector{Float64}}}
    iscontinuous::Bool
    labeldict::Dict{String, String}
    VarInfo5() = new()
    VarInfo5(x...) = new(x...)
end
JSON3.StructType(::Type{VarInfo5}) = JSON3.Mutable()

mutable struct VarInput
    name_user::String
    varID::String
    unit::String #family, head, or spouse
    VarInput() = new()
    VarInput(x...) = new(x...)
end
JSON3.StructType(::Type{VarInput}) = JSON3.Mutable()


text:
"""
    is_ordered(T::Type)

Return true iff leaftype `T` is ordered (either totally or partially).
Method: Check if method `<(::T,::T) is available at runtime
"""
is_ordered(::Type{T}) where {T<:Any} = checkm2(T, <) || checkm(T, isless)

"""
    is_partially_ordered(T::Type)

Return true iff leaftype `T` is partially ordered.
Method: Check if method `<(::T,::T) is available and `isless(::T,::T)` is not.
"""
is_partially_ordered(::Type{T}) where {T<:Any} = checkm2(T, <) && ! checkm(T, isless)

"""
    is_totally_ordered(T::Type)

Return true iff `T` is totally ordered.
Method: Check if methods `<(::T,::T) ans `isless(::T,::T)` are available.
"""
is_totally_ordered(::Type{T}) where {T<:Any} = checkm(T, isless)

"""
  is_iterable(::Type)

Return true iff leaf type `T` has methods `start(T)`, `next(T, Any)`, and `done(T, Any)`.
"""
is_iterable(::Type{T}) where {T<:Any} = hasmethod(start, (T,)) && hasmethod(done, (T, Any)) && hasmethod(next, (T, Any))

### Implementation details
# check method for (T,T) is callable
checkm(t::Type, f::Function) = hasmethod(f, (t, t))

# check if method is callable for (T,T) but it is not the default implementation fo (Any,Any)
function checkm2(t::Type, f::Function)
  try
    m1 = which(f, (t, t))
    m2 = which(f, (Any, Any))
    m1 != m2
  catch
    false
  end
end
"""
    hasmethod(f::Function, (type,...)) -> Bool

Check if method of the function `f` and given argument types is callable.
"""
function hasmethod(f::Function, t = Tuple)
  try
    which(f, t)
    true
  catch
    false
  end
end

"""
  partition(f, a::AbstractArray) -> (filter(f, a), filter( !f, a))

Partition the elements of array according to predicate function `f`.
The function is only evaluated once per array element.
"""
function partition(f, As::AbstractArray)
   sel = map(f, As)::AbstractArray{Bool}
   As[sel], As[.!(sel)]
 end


text:
function alloccursin(needles, haystack)
    found = filter(x->occursin(x, haystack), needles)
    length(found) == length(needles)
end

function readpassports(file)
    s = read(file, String)
    split(s, "\n\n", keepempty=false)
end

function containsrequired(passports)
    fields = ["byr", "iyr", "eyr", "hcl", "ecl", "pid"]
    filter(x->alloccursin(fields, x), passports)
end

day04part1(passports) = length(containsrequired(passports))

function passportfield(field, passport)
    rx = Regex("$field:([a-zA-Z0-9#]*)")
    m = match(rx, passport)
    m.captures[1]
end

function validate_string_as_number(s, min, max)
    n = tryparse(Int, s)
    if n == nothing
        return false
    end
    if n < min || n > max
        return false
    end
    true
end

function validateintfield(field, passport, min, max)
    v = passportfield(field, passport)
    validate_string_as_number(v, min, max)
end

validatebyr(passport) = validateintfield("byr", passport, 1920, 2002)

validateiyr(passport) = validateintfield("iyr", passport, 2010, 2020)

validateeyr(passport) = validateintfield("eyr", passport, 2020, 2030)

function validatehgt(passport)
    v = passportfield("hgt", passport)
    trimmed = v[1:end-2]
    if endswith(v, "cm")
        return validate_string_as_number(trimmed, 150, 193)
    end
    if endswith(v, "in")
        return validate_string_as_number(trimmed, 59, 76)
    end
    false
end

function validatehcl(passport)
    v = passportfield("hcl", passport)
    occursin(r"^#[0-9a-f]{6}$", v)
end

function validateecl(passport)
    v = passportfield("ecl", passport)
    occursin(r"^amb|blu|brn|gry|grn|hzl|oth$", v)
end

function validatepid(passport)
    v = passportfield("pid", passport)
    if length(v) != 9
        return false
    end
    tryparse(Int, v) != nothing
end

validatepassport(p) =
    validatebyr(p) &&
    validateiyr(p) &&
    validateeyr(p) &&
    validatehgt(p) &&
    validatehcl(p) &&
    validateecl(p) &&
    validatepid(p)

filtervalid(passports) = filter(x->validatepassport(x), passports)

day04part2(passports) =
    passports |>
    containsrequired |>
    filtervalid |>
    length

function day04(file)
    passports = readpassports(file)
    p1 = day04part1(passports)
    p2 = day04part2(passports)
    (p1, p2)
end

#=
s1 = """
ecl:gry pid:860033327 eyr:2020 hcl:#fffffd
byr:2937 iyr:2017 cid:147 hgt:183cm
"""
alloccursin(["byr", "iyr", "eyr", "hcl", "ecl", "pid"], s1)

file = "day04/input.txt"
passports = readpassports(file)

day04part1(passports)

parse(Int, "hello")

passportfield("ecl", s1)

validatebyr(s1)

day04part2(passports)

=#


text:
using Test
using DataDeps
using ParameterEstimocean

always_accept = get(ENV, "DATADEPS_ALWAYS_ACCEPT", false)
ENV["DATADEPS_ALWAYS_ACCEPT"] = true

function run_datadep_test(dep_str)
    @testset "DataDep $dep_str" begin
        @datadep_str dep_str
        @test true
    end
    return nothing
end

for dep in ["two_day_suite_1m",
            "two_day_suite_2m",
            "two_day_suite_4m",
            "four_day_suite_1m",
            "four_day_suite_2m",
            "four_day_suite_4m",
            "six_day_suite_1m",
            "six_day_suite_2m",
            "six_day_suite_4m"]

    run_datadep_test(dep)
end

ENV["DATADEPS_ALWAYS_ACCEPT"] = always_accept


text:
using HTTP, JSON

const VegaLiteWebsocketFrontEndTemplate = """
<!DOCTYPE html>

  <meta charset="utf-8">

"""

const VegaLiteFitnessOverTimeSpecTemplate = """
{
    "\$schema": "https://vega.github.io/schema/vega-lite/v4.json",
    "description": "Fitness value over time",
    "width": %%WIDTH%%,
    "height": %%HEIGHT%%,
    "padding": {"left": 20, "top": 10, "right": 10, "bottom": 20},
    "data": {
      "name": "table"
    },
    "mark": "line",
    "encoding": {
        "x": {
            "field": "Time",
            "type": "quantitative"
        },
        "y": {
            "field": "Fitness",
            "type": "quantitative",
            "scale": {"type": "log"}
        }
    }
}
"""

rand_websocket_port() = 9000 + rand(0:42)

function vegalite_websocket_frontend(vegalitespec::String = VegaLiteFitnessOverTimeSpecTemplate;
    width::Int = 800, height::Int = 600, port::Int = rand_websocket_port())
    res = replace(VegaLiteWebsocketFrontEndTemplate, "%%VEGALITESPEC%%" => vegalitespec)
    res = replace(res, "%%WIDTH%%" => string(width))
    res = replace(res, "%%HEIGHT%%" => string(height))
    replace(res, "%%SOCKETPORT%%" => string(port))
end

function serve_html_file(file::String)
    HTTP.serve() do request::HTTP.Request
        try
            return HTTP.Response(file)
        catch e
            return HTTP.Response(404, "Error: $e")
        end
    end
end

# Serve a VegaLite front-end on html and then push updates to data
# over a websocket so that the VegaLite graph is updated.
mutable struct VegaLiteGraphFitnessGraph
    verbose::Bool
    httpport::Int
    websocketport::Int
    mindelay::Float64
    data::Vector{Any}
    last_sent_index::Int
    starttime::Float64
    function VegaLiteGraphFitnessGraph(verbose::Bool = true,
        httpport::Int = 8081, websocketport::Int = 9000, mindelay::Float64 = 1.0)
        @assert mindelay > 0.0
        new(verbose, httpport, websocketport, mindelay, Any[], 0, 0.0)
    end
end

timestamp(t = time()) = Libc.strftime("%Y-%m-%d %H:%M.%S", t)
printmsg(vlg::VegaLiteGraphFitnessGraph, msg) = vlg.verbose ? println(timestamp(), ": ", msg) : nothing

function add_data!(vlg::VegaLiteGraphFitnessGraph, newentry::Dict)
    if length(vlg.data) < 1
        vlg.starttime = time()
    end
    if !haskey(newentry, "Time")
        newentry["Time"] = time() - vlg.starttime
    end
    printmsg(vlg, "Adding data $newentry")
    push!(vlg.data, newentry)
end

function serve(vlg::VegaLiteGraphFitnessGraph)
    @async websocket_push_data_loop(vlg.websocketport, vlg.mindelay) do ws
        send_new_data_on_socket(vlg, ws)
    end
    frontend_html = vegalite_websocket_frontend(; port = vlg.websocketport)
    @async serve_html_file(frontend_html)
    println("Serving frontend on http://127.0.0.1:$(vlg.httpport)")
    sleep(3.0) # To give people time to copy the url...
end

hasnewdata(vlg::VegaLiteGraphFitnessGraph) = length(vlg.data) > vlg.last_sent_index

function send_new_data_on_socket(vlg::VegaLiteGraphFitnessGraph, ws)
    if hasnewdata(vlg)
        len = length(vlg.data)
        newdata = vlg.data[(vlg.last_sent_index+1):len]
        printmsg(vlg, "Sending data $newdata")
        write(ws, JSON.json(newdata))
        vlg.last_sent_index = len
    end
end

function websocket_push_data_loop(fn::Function, wsport::Int, mindelay = 1.0)
    HTTP.WebSockets.listen("127.0.0.1", UInt16(wsport)) do ws
        while true
            fn(ws) # Call with websocket
            sleep(mindelay + rand())
        end
    end
end


text:
"""
Nestedness of a single axis (called internally by `η`)
"""
function η_axis(N::AbstractBipartiteNetwork)
  S = richness(N; dims=1)
  n = vec(sum(N.A; dims=2))
  num = 0.0
  den = 0.0
  @simd for j in 2:S
    @simd for i in 1:(j-1)
      @inbounds num += sum(N[i,:].*N[j,:])
      @inbounds den += min(n[i], n[j])
    end
  end
  return sum(num ./ den)
end

"""
    η(N::T, dims::Union{Nothing,Integer}=nothing) where {T <: Union{BipartiteNetwork, BipartiteProbaNetwork}}

Returns the nestedness of a margin of the network, using η. The second argument
can be `1` (for nestedness of rows/top level) or `2` (for nestedness of
columns/bottom level). Leaving it at `nothing` will measure the nestedness of
the entire network.

#### References

Bastolla, U., Fortuna, M.A., Pascual-García, A., Ferrera, A., Luque, B.,
Bascompte, J., 2009. The architecture of mutualistic networks minimizes
competition and increases biodiversity. Nature 458, 1018–1020.
https://doi.org/10.1038/nature07950

Poisot, T., Cirtwill, A.R., Cazelles, K., Gravel, D., Fortin, M.-J., Stouffer,
D.B., 2016. The structure of probabilistic networks. Methods in Ecology and
Evolution 7, 303–312. https://doi.org/10.1111/2041-210X.12468
"""
function η(N::T; dims::Union{Nothing,Integer}=nothing) where {T <: Union{BipartiteNetwork, BipartiteProbabilisticNetwork}}
  dims == 1 && return η_axis(N)
  dims == 2 && return η_axis(permutedims(N))
  if dims === nothing
    return (η(N; dims=1) + η(N; dims=2))/2.0
  end
  throw(ArgumentError("dims can only be 1 (nestedness of rows) or 2 (nestedness of columns), you used $(dims)"))
end

"""
WNODF of a single axis
"""
function nodf_axis(N::BipartiteQuantitativeNetwork)

  # Get the row order
  row_order = sortperm(vec(sum(N.A; dims=2)), rev=true)

  # Extract the ordered matrix as floating point values, so that all other
  # measures will work for both the quanti and deterministic networks
  A = map(Float64, N.A[row_order,:])

  # Initialize the value
  WNODFr = 0.0
  for i in 1:(size(A, 1)-1)
    for j in (i+1):size(A, 1)
      if (sum(A[i,:]) .> sum(A[j,:]))
        Nj = sum(A[j,:] .> 0.0)
        kij = 0.0
        # The following is only applied to interactions that are non-0
        # in j, which is not really clear in the original paper
        for l in eachindex(A[i,:])
          if (A[j,l] > 0.0)
            if A[j,l] < A[i,l]
              kij += 1.0
            end
          end
        end
        WNODFr += kij/Nj
      end
    end
  end

  # Return the value
  return WNODFr

end

"""
NODF of a single axis
"""
function nodf_axis(N::BipartiteNetwork)

  # Get the row order
  row_order = sortperm(vec(sum(N.A; dims=2)), rev=true)

  # Extract the ordered matrix as floating point values, so that all other
  # measures will work for both the quanti and deterministic networks
  A = N.A[row_order,:]
  d = float(vec(sum(A; dims=2)))

  # Initialize the value
  NODFr = 0.0
  for i in 1:(size(A, 1)-1)
    for j in (i+1):size(A, 1)
      DFpaired = d[j] < d[i] ? 1.0 : 0.0
      Npaired = sum(A[i,:] .& A[j,:]) / d[j]
      NODFr += (DFpaired * Npaired)
    end
  end

  # Return the value
  return NODFr
end

"""
    nodf(N::T; dims::Union{Nothing,Integer}=nothing) where {T <: Union{BipartiteNetwork,BipartiteQuantitativeNetwork}}

Returns `nodf` for a margin of the network. The `i` argument can be 1 for
top-level, 2 for bottom level, and the function will throw an `ArgumentError` if
an invalid value is used. For quantitative networks, *WNODF* is used.

#### References

Almeida-Neto, M., Guimarães, P.R., Loyola, R.D., Ulrich, W., 2008. A consistent
metric for nestedness analysis in ecological systems: reconciling concept and
measurement. Oikos 117, 1227–1239.
https://doi.org/10.1111/j.0030-1299.2008.16644.x

Almeida-Neto, M., Ulrich, W., 2011. A straightforward computational approach for
measuring nestedness using quantitative matrices. Environmental Modelling &
Software 26, 173–178. https://doi.org/10.1016/j.envsoft.2010.08.003
"""
function nodf(N::T; dims::Union{Nothing,Integer}=nothing) where {T <: Union{BipartiteNetwork,BipartiteQuantitativeNetwork}}
  if dims == 1
    val = nodf_axis(N)
    correction = (richness(N; dims=1) * (richness(N; dims=1) - 1))
    return (val+val)/correction
  end
  if dims == 2
    val = nodf_axis(permutedims(N))
    correction = (richness(N; dims=2) * (richness(N; dims=2) - 1))
    return (val+val)/correction
  end
  if dims === nothing
    return (nodf(N; dims=1)+nodf(N; dims=2))/2.0
  end
  throw(ArgumentError("dims can only be 1 (nestedness of rows) or 2 (nestedness of columns), you used $(dims)"))
end


text:
using DiffPrivacy, Test

@info "Testing one parallel database"
x = [1,2,3,4,5,6,7,8,9]
y = DiffPrivacy.get_one_parallel_database(x, 2)
@test y == [1 3 4 5 6 7 8 9]

@info "Testing multiple parallel databases"
x = [1,2,3,4]
y = DiffPrivacy.get_multiple_parallel_databases(x)
@test y == [[2 3 4],[1 3 4],[1 2 4],[1 2 3]]

@info "Testing sum sensitivity"
x = [1,2,3,4]
y = DiffPrivacy.sum_sensitivity(x)
@test y == 4

@info "Testing mean sensitivity"
x = [3,6,9,12]
y = DiffPrivacy.mean_sensitivity(x)
@test y == 1.5

@info "All tests completed"


text:
# abstract type AbstractInstrument end  # in abstract_types.jl
struct HPFType    <: AbstractInstrument end
struct NEIDType   <: AbstractInstrument end
struct EXPRESType <: AbstractInstrument end
struct HARPSType  <: AbstractInstrument end
#struct HARPSNType <: AbstractInstrument end

export AbstractInstrument, HPFType, NEIDType, HARPSType
#export HARPSNType

#instrument_symbol_list = Symbol[]
instrument_symbol_list = Symbol[:HPF, :NEID, :EXPRES, :HARPSN, :HARPS ]
export instrument_symbol_list


text:
# Dev code for free space solver for modified stokes
push!(LOAD_PATH, string(pwd(),"/src"))

using ModifiedStokesSolver
using Base.Test

@testset "PeriodicSolver" begin
    # Setup problem
    srand(1)
    L = 4.0
    lambda = 5.6

    sigma = L/20
    x0 = 0.12376938253
    y0 = 0.21938752311

    xi = 1/sqrt(2*sigma^2)
    # Reference solution derived by applying modified biharmonic to Gaussian, plus a constant field
    # ufuncr(r) = exp(-xi^2*r^2)
    ffuncr(r) = -4*xi^2*exp(-xi^2*r^2)*( lambda^2*(1-r^2*xi^2) + 4*xi^2*(2 - 4*r^2*xi^2 + r^4*xi^4) ) - lambda^2
    s1r(r) = -2*xi^2*exp(-xi^2*r^2)*(2*xi^2*r^2-1) + 1
    s2r(r) = 4*xi^4*exp(-xi^2*r^2)
    pr(r) = 2*xi^2*exp(-xi^2*r^2)*(lambda^2 + 8*xi^2 - 4*xi^4*r^2)

    r1(x,y) = x-x0
    r2(x,y) = y-y0
    r(x,y) = sqrt( (x-x0)^2 + (y-y0)^2 )

    c1 = rand()
    c2 = rand()
    ffunc1(x, y) = c1*ffuncr(r(x,y))
    ffunc2(x, y) = c2*ffuncr(r(x,y))

    ufunc1(x, y) = -1*(c1*s1r(r(x,y)) + s2r(r(x,y))*r1(x,y)*(c1*r1(x,y) + c2*r2(x,y)))
    ufunc2(x, y) = -1*(c2*s1r(r(x,y)) + s2r(r(x,y))*r2(x,y)*(c1*r1(x,y) + c2*r2(x,y)))
    pfunc(x, y) = -1*(pr(r(x,y))*(c1*r1(x,y) + c2*r2(x,y)))

    function solve(N)
        # Grid domain
        x = linspace(-L/2, L/2, N+1)
        x = x[1:end-1] # Note that intervals are closed/open!
        y = x
        X,Y = ndgrid(x, y)
        # Nonuniform test points
        M = 1000
        xt = L*(1/2-rand(M, 1))
        yt = L*(1/2-rand(M, 1))

        # Get rhs
        F1 = ffunc1.(X,Y)
        F2 = ffunc2.(X,Y)

        # Get reference solution
        U1ref = ufunc1.(X,Y)
        U2ref = ufunc2.(X,Y)

        # Compute solution with FFT method
        U1, U2, ut1, ut2, U1x, U1y, U2x, U2y = per_modstokes(F1, F2, L, lambda, xt, yt; zerok0=false, ifgrad=true)

        # Compute error on grid
        E1 = U1 - U1ref
        E2 = U2 - U2ref
        Unorm = norm(sqrt.( vec(U1).^2 + vec(U2).^2 ), Inf)
        relerr1 = norm(E1[:], Inf) / Unorm
        relerr2 = norm(E2[:], Inf) / Unorm
        relerr = max(relerr1, relerr2)

        # Compute error at nonuniform test points
        enu1 = ut1-ufunc1.(xt,yt)
        enu2 = ut2-ufunc2.(xt,yt)
        relerrnu1 = norm(enu1, Inf) / norm(ut1, Inf)
        relerrnu2 = norm(enu2, Inf) / norm(ut2, Inf)
        relerrnu = max(relerrnu1, relerrnu2)

        # Compute gradient with spectral differentiation and compare
        U1xd, U1yd = per_gradient(U1, L)
        U2xd, U2yd = per_gradient(U2, L)
        @test maximum(abs.( U1x-U1xd )) / maximum(abs.( U1x )) < 21e-14
        @test maximum(abs.( U1y-U1yd )) / maximum(abs.( U1y )) < 2e-14
        @test maximum(abs.( U2x-U2xd )) / maximum(abs.( U2x )) < 2e-14
        @test maximum(abs.( U2y-U2yd )) / maximum(abs.( U2y )) < 2e-14

        # Check Divergence
        Udiv = U1x+U2y
        diverr_rel = maximum(abs.(Udiv))/Unorm
        @test diverr_rel < 5e-14

        # Compute numerical derivatives to check that we actually satisfy the PDE
        # using finite differences
        h = X[2]-X[1]
        int = 2:N-1
        U1x = (U1[int+1,int]-U1[int-1,int])/(2*h)
        U2y = (U2[int,int+1]-U2[int,int-1])/(2*h)

        P = pfunc.(X,Y)
        Pnorm = norm( vec(P), Inf)

        Udiv = U1x+U2y
        diverr_rel = maximum(abs.(Udiv))/Unorm

        LP = (-4*P[int,int] + P[int+1,int] + P[int-1,int] + P[int,int+1] + P[int,int-1])/h^2
        Px = (P[int+1,int]-P[int-1,int])/(2*h)
        Py = (P[int,int+1]-P[int,int-1])/(2*h)
        LU1 = (-4*U1[int,int] + U1[int+1,int] + U1[int-1,int] + U1[int,int+1] + U1[int,int-1])/h^2
        LU2 = (-4*U2[int,int] + U2[int+1,int] + U2[int-1,int] + U2[int,int+1] + U2[int,int-1])/h^2

        pde1 = lambda^2*U1[int,int] - LU1 + Px - F1[int,int]
        pde2 = lambda^2*U2[int,int] - LU2 + Py - F2[int,int]

        pdeerr1_rel = maximum(abs.(pde1))/Unorm
        pdeerr2_rel = maximum(abs.(pde2))/Unorm
        pdeerr_rel = max(pdeerr1_rel, pdeerr2_rel)

        return relerr, relerrnu, diverr_rel, pdeerr_rel
    end

    # Solve above problem for a variety of N
    iters = 3;
    diverr_list = zeros(iters)
    pdeerr_list = zeros(iters)
    N_list = zeros(iters)
    for i=0:iters-1
        N = 64*2^i
        relerr, relerrnu, diverr_rel, pdeerr_rel = solve(N)
        @test relerr < 1e-13
        @test relerrnu < 1e-13
        N_list[i+1] = N
        diverr_list[i+1] = diverr_rel
        pdeerr_list[i+1] = pdeerr_rel
    end
    # Estimate convergence constant of FD test (should be 2)
    A = [ones(iters) -log.(N_list)]
    div_conv = (A\log.(diverr_list))[2]
    pde_conv = (A\log.(pdeerr_list))[2]
    @test abs(div_conv - 2) < 0.05
    @test abs(pde_conv - 2) < 0.05
end


text:
using Hungarian
using Munkres
using Test
using LinearAlgebra
using DelimitedFiles

@testset "simple examples" begin
    A = [ 0.891171  0.0320582   0.564188  0.8999    0.620615;
          0.166402  0.861136    0.201398  0.911772  0.0796335;
          0.77272   0.782759    0.905982  0.800239  0.297333;
          0.561423  0.170607    0.615941  0.960503  0.981906;
          0.748248  0.00799335  0.554215  0.745299  0.42637]

    assign, cost = @inferred hungarian(A)
    @test assign == [2, 3, 5, 1, 4]

    B = [ 24     1     8;
           5     7    14;
           6    13    20;
          12    19    21;
          18    25     2]

    assign, cost = hungarian(B)
    @test assign == [2, 1, 0, 0, 3]
    @test cost == 8

    assign, cost = hungarian(transpose(B))
    @test assign == [2, 1, 5]
    @test cost == 8

    assign, cost = hungarian(ones(5,5) - I)
    @test assign == [1, 2, 3, 4, 5]
    @test cost == 0
end

@testset "forbidden edges" begin
    # result checked against Python package munkres: https://github.com/bmc/munkres/blob/master/munkres.py
    # Python code:
    #   m = Munkres()
    #   matrix = [[DISALLOWED, 1, 1], [1, 0, 1], [1, 1, 0]]
    #   m.compute(matrix)
    # Result: [(0, 1), (1, 0), (2, 2)]
    A = Union{Int, Missing}[missing 1 1; 1 0 1; 1 1 0]
    assign, cost = hungarian(A)
    @test assign == [2, 1, 3]
    @test cost == 2
end

@testset "issue #6" begin
    # load the cost matrix produced from the `solve_LSAP` function contained in the
    # [clue package](https://cran.r-project.org/web/packages/clue/index.html) in R.
    costMatrix = readdlm(joinpath(@__DIR__, "R-clue-cost.txt"))
    assignmentH, costH = hungarian(copy(transpose(costMatrix)))
    assignmentR = vec(readdlm(joinpath(@__DIR__, "R-clue-solution.txt"), Int64))
    costR = sum(costMatrix[ii, assignmentR[ii]] for ii in 1:size(costMatrix,1))
    @test costH ≤ costR
end

@testset "UInt16" begin
    M=UInt16[28092 44837 19882 39481 59139;
             26039 46258 38932 51057 9;
             11527 59487 61993 29072 8734;
             10691 16977 12796 16370 14266;
             5199  42319 34194 41332 16472]
    assign,cost=hungarian(M)
    @test assign == [3, 5, 4, 2, 1]
    @test cost   == 71139
end

@testset "UInt8" begin
    M=UInt8[67  228 135 197 244;
            112 44  84  206 31;
            225 103 231 225 227;
            170 37  135 9   130;
            110 22  133 77  96]
    assign,cost=hungarian(M)
    @test assign == [1, 5, 2, 4, 3]
    @test cost   == 343
end

@testset "Munkres.jl" begin
    A = rand(300, 300)
    assignH, costH = hungarian(A)
    assignM = munkres(A)

    for i = 1:100
        A = rand(50, 50)
        assignH, costH = hungarian(A)
        assignM = munkres(A)

        costM = 0
        for i in zip(1:size(A, 1), assignM)
            if i[2] != 0
                costM += A[i...]
            end
        end

        @test assignH == assignM
        @test costH == costM
    end
end


text:
####################################################################################################
"""
Plot the branch of solutions during the continuation
"""
function plotBranchCont(contres::ContResult, sol::BorderedArray, contparms, plotuserfunction)
	colorbif = Dict(:fold => :black, :hopf => :red, :bp => :blue, :nd => :magenta, :none => :yellow, :ns => :orange, :pd => :green)
	branch = contres.branch

	if contparms.computeEigenValues == false
		l =  Plots.@layout [a{0.5w} [b; c]]
	else
		l =  Plots.@layout [a{0.5w} [b; c]; e{0.2h}]
	end
	Plots.plot(layout = l)

	# plot the branch of solutions
	plotBranch!(contres; filterbifpoints = true, putbifptlegend = false, xlabel="p",  ylabel="||x||", label="", subplot=1)
	plot!(branch[1, :],	 xlabel="it", ylabel="p", label="", subplot=2)

	# add the bifurcation points along the branch
	# if length(contres.bifpoint)>1
	# 	scatter!(map(x -> x.idx, contres.bifpoint[2:end]),
	# 			 map(x -> x.param, contres.bifpoint[2:end]),
	# 			label="", color = map(x -> colorbif[x.type], contres.bifpoint[2:end]), markersize=3, markerstrokewidth=0, subplot=2)
	# end

	if contparms.computeEigenValues
		ii = length(contres.eig)
		scatter!(real.(contres.eig[ii][1]), imag.(contres.eig[ii][1]), subplot=4, label="", markersize = 1, color=:black)
	end

	plotuserfunction(sol.u, subplot = 3)

	display(title!(""))
end

"""
	plotBranch(contres::ContResult, plot_fold = true; kwargs...)

Plot the branch of solutions from a `ContResult`. You can also pass parameters like `plotBranch(br, marker = :dot)`.
For the continuation diagram, the legend is as follows `(:fold => :black, :hopf => :red, :bp => :blue, :nd => :magenta, :none => :yellow)`
"""
function plotBranch(contres, plot_fold = true; kwargs...)
	# we do not specify the type of contres, not convenient when using JLD2
	plot()
	plotBranch!(contres, plot_fold; kwargs...)
end

"""
	plotBranch(brs::Vector, plot_fold = true; putbifptlegend = true, filterbifpoints = false, kwargs...)

Plot all the branches contained in `brs` in a single figure. Convenient when many bifurcation diagram have been computed.
"""
function plotBranch(brs::Vector, plot_fold = true; putbifptlegend = true, filterbifpoints = false, kwargs...)
	# we do not specify the type of contres, not convenient when using JLD2
	colorbif = Dict(:fold => :black, :hopf => :red, :bp => :blue, :nd => :magenta, :none => :yellow, :ns => :orange, :pd => :green)
	plotBranch(brs[1], plot_fold; putbifptlegend = false, filterbifpoints = filterbifpoints, kwargs...)
	bp = Set(unique([pt.type for pt in brs[1].bifpoint]))
	for ii=2:length(brs)
		plotBranch!(brs[ii], plot_fold; putbifptlegend = false, filterbifpoints = filterbifpoints, kwargs...) |> display
		for pt in brs[ii].bifpoint
			push!(bp, pt.type)
		end
	end
	# add legend for bifurcation points
	if putbifptlegend && length(bp) > 0
		for pt in bp
			scatter!([], [], color = colorbif[pt], label = "$pt")
		end
	end
end

function filterBifurcations(bifpt)
	# this function filters Fold points and Branch points which are located at the same/previous/next point
	res = [(type = :none, idx = 1, param = 1., printsol = 1.)]
	ii = 1
	while ii <= length(bifpt) - 1
		if (abs(bifpt[ii].idx - bifpt[ii+1].idx) <= 1) && bifpt[ii].type ∈ [:fold, :bp]
			if (bifpt[ii].type == :fold && bifpt[ii].type == :bp) ||
				(bifpt[ii].type == :bp && bifpt[ii].type == :fold)
				push!(res, (type = :fold, idx = bifpt[ii].idx, param = bifpt[ii].param, printsol = bifpt[ii].printsol) )
			else
				push!(res, (type = bifpt[ii].type, idx = bifpt[ii].idx, param = bifpt[ii].param, printsol = bifpt[ii].printsol) )
				push!(res, (type = bifpt[ii+1].type, idx = bifpt[ii+1].idx, param = bifpt[ii+1].param, printsol = bifpt[ii+1].printsol) )
			end
			ii += 2
		else
			push!(res, (type = bifpt[ii].type, idx = bifpt[ii].idx, param = bifpt[ii].param, printsol = bifpt[ii].printsol) )
			ii += 1
		end
	end
	0<ii<=length(bifpt) &&	push!(res, (type = bifpt[ii].type, idx = bifpt[ii].idx, param = bifpt[ii].param, printsol = bifpt[ii].printsol) )

	# for p in (res[2:end])
	# 	println(p)
	# end
	res[2:end]
end

"""
	plotBranch!(contres, plot_fold = true; putbifptlegend = true, filterbifpoints::Bool = false, kwargs...)

Append to the current plot, the plot of the branch of solutions from a `ContResult`. You can also pass parameters like `plotBranch!(br, marker = :dot)`. Options to filter the bifurcation points (which are mostly guesses) are provided. For example, `filterbifpoints = true` merges the nearby Fold and Branch points.
"""
function plotBranch!(contres, plot_fold = true; putbifptlegend = true, filterbifpoints::Bool = false, kwargs...)
	# we do not specify the type of contres, not convenient when using JLD2
	colorbif = Dict(:fold => :black, :hopf => :red, :bp => :blue, :nd => :magenta, :none => :yellow, :ns => :orange, :pd => :green)
	branch = contres.branch
	if length(contres.stability) > 2
		# plot!(branch[1, :], branch[2, :], linestyle = map(x -> isodd(x) ? :solid : :dash, contres.stability); kwargs...)
		# plot!(branch[1, :], branch[2, :], color = map(x -> isodd(x) ? :green : :red, contres.stability); kwargs...)
		plot!(branch[1, :], branch[2, :], linewidth = map(x -> isodd(x) ? 3.0 : 1.0, contres.stability); kwargs...)
	else
		plot!(branch[1, :], branch[2, :]; kwargs...)
	end
	scatter!([branch[1, end]], [branch[2, end]], marker=:cross, subplot=1, label = "")
	# add the bifurcation points along the branch
	bifpoints = vcat(contres.bifpoint, filter(x->x.type != :none, contres.foldpoint))
	if length(bifpoints) >= 1
		id = 1
		bifpoints[1].type == :none ? id = 2 : id = 1
		if plot_fold
			bifpt = bifpoints[id:end]
		else
			bifpt = filter(x->x.type != :fold, bifpoints[id:end])
		end
		if filterbifpoints == true
			bifpt = filterBifurcations(bifpt)
		end
		scatter!(map(x -> x.param, bifpt), map(x -> x.printsol, bifpt), label="", color = map(x->colorbif[x.type], bifpt), markersize=2, markerstrokewidth=1 ; kwargs...)
	end
	# add legend for bifurcation points
	if putbifptlegend && length(bifpoints) >= 1
		bp = unique([pt.type for pt in bifpt if pt.type != :none])
		(length(bp) == 0) && return
		for pt in bp
			scatter!([], [], color=colorbif[pt], label="$pt")
		end
	end
end


text:
using CSyntax

@static if Sys.isapple()
    const VERSION_MAJOR = 4
    const VERSION_MINOR = 1
end

include(joinpath(@__DIR__, "glutils.jl"))

# print errors in shader compilation
function shader_info_log(shader::GLuint)
    max_length = GLint(0)
    @c glGetShaderiv(shader, GL_INFO_LOG_LENGTH, &max_length)
    actual_length = GLsizei(0)
    log = Vector{GLchar}(undef, max_length)
    @c glGetShaderInfoLog(shader, max_length, &actual_length, log)
    @info "shader info log for GL index $shader: $(String(log))"
end

# print errors in shader linking
function programme_info_log(program::GLuint)
    max_length = GLint(0)
    @c glGetShaderiv(program, GL_INFO_LOG_LENGTH, &max_length)
    actual_length = GLsizei(0)
    log = Vector{GLchar}(undef, max_length)
    @c glGetShaderInfoLog(program, max_length, &actual_length, log)
    @info "program info log for GL index $program: $(String(log))"
end

# validate shader program
function is_valid(program::GLuint)
    params = GLint(-1)
    glValidateProgram(program)
    @c glGetProgramiv(program, GL_VALIDATE_STATUS, &params)
    @info "program $program GL_VALIDATE_STATUS = $params"
    params == GL_TRUE && return true
    programme_info_log(program)
    return false
end

# print verbose infos
function print_all(shader_prog::GLuint)
    params = GLint(-1)

    @debug "-------------------------"
    @debug "Shader programme $shader_prog info:"
    @c glGetProgramiv(shader_prog, GL_LINK_STATUS, &params)
    @debug "GL_LINK_STATUS = $params"

    @c glGetProgramiv(shader_prog, GL_ATTACHED_SHADERS, &params)
    @debug "GL_ATTACHED_SHADERS = $params"

    @c glGetProgramiv(shader_prog, GL_ACTIVE_ATTRIBUTES, &params)
    @debug "GL_ACTIVE_ATTRIBUTES = $params"

    max_length = GLint(0)
    @c glGetProgramiv(shader_prog, GL_ACTIVE_ATTRIBUTE_MAX_LENGTH, &max_length)
    name = Vector{GLchar}(undef, max_length)
    for i in 0:params-1
        actual_length = GLsizei(0)
        size = GLint(0)
        type = GLenum(0)
        @c glGetActiveAttrib(shader_prog, i, max_length, &actual_length, &size, &type, name)
        if size > 1
            for j = 0:size-1
                longname = @sprintf "%s[%i]" name j
                location = glGetAttribLocation(shader_prog, longname)
                @debug "  $i): type -> $(GLENUM(type).name), name -> $(String(longname)), location -> $location."
            end
        else
            location = glGetAttribLocation(shader_prog, name)
            @debug "  $i): type -> $(GLENUM(type).name), name -> $(String(name)), location -> $location."
        end
    end

    @c glGetProgramiv(shader_prog, GL_ACTIVE_UNIFORMS, &params)
    @debug "GL_ACTIVE_UNIFORMS = $(paramsRef[])"
    for i in 0:params-1
        max_length = GLint(0)
        @c glGetProgramiv(shader_prog, GL_ACTIVE_UNIFORM_MAX_LENGTH, &max_length)
        name = Vector{GLchar}(undef, max_length)
        actual_length = GLsizei(0)
        size = GLint(0)
        type = GLenum(0)
        @c glGetActiveUniform(shader_prog, i, max_length, &actual_length, &size, &type, name)
        if size > 1
            for j = 0:size-1
                longname = @sprintf "%s[%i]" name j
                location = glGetUniformLocation(shader_prog, longname)
                @debug "  $i): type -> $(GLENUM(type).name), name -> $(String(longname)), location -> $location."
            end
        else
            location = glGetUniformLocation(shader_prog, name)
            @debug "  $i): type -> $(GLENUM(type).name), name -> $(String(name)), location -> $location."
        end
    end

    programme_info_log(shader_prog)
end

# init window
width, height = fb_width, fb_height = 640, 480
window = startgl(width, height)

glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LESS)

# load shaders from file
const vert_source = read(joinpath(@__DIR__, "shader.vert"), String)
const frag_source = read(joinpath(@__DIR__, "shader.frag"), String)

# compile shaders and check for shader compile errors
vert_shader = glCreateShader(GL_VERTEX_SHADER)
glShaderSource(vert_shader, 1, Ptr{GLchar}[pointer(vert_source)], C_NULL)
glCompileShader(vert_shader)
# get shader compile status
result = GLint(-1)
@c glGetShaderiv(vert_shader, GL_COMPILE_STATUS, &result)
if result != GL_TRUE
    shader_info_log(vert_shader)
    @error "GL vertex shader(index $vert_shader) did not compile."
end

frag_shader = glCreateShader(GL_FRAGMENT_SHADER)
glShaderSource(frag_shader, 1, Ptr{GLchar}[pointer(frag_source)], C_NULL)
glCompileShader(frag_shader)
# checkout shader compile status
result = GLint(-1)
@c glGetShaderiv(frag_shader, GL_COMPILE_STATUS, &result)
if result != GL_TRUE
    shader_info_log(frag_shader)
    @error "GL fragment shader(index $frag_shader) did not compile."
end

# create and link shader program
shader_prog = glCreateProgram()
glAttachShader(shader_prog, vert_shader)
glAttachShader(shader_prog, frag_shader)
glLinkProgram(shader_prog)
# checkout programe linking status
result = GLint(-1)
@c glGetProgramiv(shader_prog, GL_LINK_STATUS, &result)
if result != GL_TRUE
    programme_info_log(shader_prog)
    @error "Could not link shader programme GL index: $shader_prog"
end

# vertex data
points = GLfloat[ 0.0,  0.5, 0.0,
                  0.5, -0.5, 0.0,
                 -0.5, -0.5, 0.0]

# create buffers located in the memory of graphic card
vbo = GLuint(0)
@c glGenBuffers(1, &vbo)
glBindBuffer(GL_ARRAY_BUFFER, vbo)
glBufferData(GL_ARRAY_BUFFER, sizeof(points), points, GL_STATIC_DRAW)

# create VAO
vao = GLuint(0)
@c glGenVertexArrays(1, &vao)
glBindVertexArray(vao)
glBindBuffer(GL_ARRAY_BUFFER, vao)
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)

# specify color
color_location = glGetUniformLocation(shader_prog, "inputColour")
@assert color_location > -1
glUseProgram(shader_prog)
glUniform4f(color_location, 1.0, 0.0, 0.0, 1.0)

# verbose infos
print_all(shader_prog)
@assert is_valid(shader_prog)

# set background color to gray
glClearColor(0.2, 0.2, 0.2, 1.0)

# render
while !GLFW.WindowShouldClose(window)
    updatefps(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, width, height)
    # drawing
    glUseProgram(shader_prog)
    glBindVertexArray(vao)
    glDrawArrays(GL_TRIANGLES, 0, 3)
    # check and call events
    GLFW.PollEvents()
    # swap the buffers
    GLFW.SwapBuffers(window)
end

GLFW.DestroyWindow(window)

max_stars_repo_path: OpenGL 4 Tutorials/02_shaders/shader.jl

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_path: OpenGL 4 Tutorials/02_shaders/shader.jl

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_path: OpenGL 4 Tutorials/02_shaders/shader.jl

max_forks_repo_name: Gnimuc/Videre


text:
module TestManual

using Test
using Documenter
using Trixi2Vtk

DocMeta.setdocmeta!(Trixi2Vtk, :DocTestSetup, :(using Trixi2Vtk); recursive=true)
doctest(Trixi2Vtk, manual=false)

end # module


text:
for N in [Float64, Float32]
    # random ball
    rand(Ballp)

    # 1D Ball3
    b = Ballp(N(3), N[0], N(1))
    # dimension
    @test dim(b) == 1
    # support vector
    d = N[1]
    @test σ(d, b) == N[1]
    d = N[-1]
    @test σ(d, b) == N[-1]

    # 2D Ball3
    b = Ballp(N(3), N[0, 0], N(2))
    # dimension
    @test dim(b) == 2
    # support vector
    d = N[1, 0]
    @test σ(d, b) == N[2, 0]
    d = N[-1, 0]
    @test σ(d, b) == N[-2, 0]
    d = N[0, 1]
    @test σ(d, b) == N[0, 2]
    d = N[0, -1]
    @test σ(d, b) == N[0, -2]
    d = N[0, 0]
    @test σ(d, b) == N[0, 0]

    # constructors that fall back to specialized set types
    @test Ballp(N(1), N[3], N(4)) isa Ball1
    @test Ballp(N(2), N[3], N(4)) isa Ball2
    @test Ballp(N(Inf), N[3], N(4)) isa BallInf

    # center
    @test center(b) == b.center
    @test an_element(b) isa AbstractVector{N}

    # boundedness
    @test isbounded(b)

    # is_polyhedral
    @test !is_polyhedral(b)

    # isempty
    @test !isempty(b)

    # isuniversal
    answer, w = isuniversal(b, true)
    @test !isuniversal(b) && !answer && w ∉ b

    # membership & an_element
    @test an_element(b) ∈ b

    # translation
    @test translate(b, N[1, 2]) == Ballp(N(3), N[1, 2], N(2))
    bb = Ballp(N(3), N[0, 0], N(1))
    @test translate!(bb, N[1, 1]) == Ballp(N(3), N[1, 1], N(1)) == bb

    # projection
    b4 = Ballp(N(3), N[4, 3, 2, 1], N(2))
    @test project(b4, [2, 4]) == Ballp(N(3), N[3, 1], N(2))
end


text:
###
### Predicates & comparisons
###
### Generalised operations on biological sequence types.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

function Base.cmp(seq1::BioSequence, seq2::BioSequence)
    m = lastindex(seq1)
    n = lastindex(seq2)
    for i in 1:min(m, n)
        c = cmp(inbounds_getindex(seq1, i),
                inbounds_getindex(seq2, i))
        if c != 0
            return c
        end
    end
    return cmp(m, n)
end

function Base.:(==)(seq1::BioSequence, seq2::BioSequence)
    return eltype(seq1)    == eltype(seq2) &&
           length(seq1)    == length(seq2) &&
           cmp(seq1, seq2) == 0
end

Base.isless(seq1::BioSequence, seq2::BioSequence) = cmp(seq1, seq2) < 0

Base.isempty(seq::BioSequence) = length(seq) == 0

"""
    isrepetitive(seq::BioSequence, n::Integer=length(seq))

Return `true` if and only if `seq` contains a repetitive subsequence of length `≥ n`.
"""
function isrepetitive(seq::BioSequence, n::Integer = length(seq))
    if n < 0
        error("repetition must be non-negative")
    elseif isempty(seq)
        return n == 0
    end

    rep = 1
    if rep ≥ n
        return true
    end
    last = first(seq)
    for i in 2:lastindex(seq)
        x = seq[i]
        if x == last
            rep += 1
            if rep ≥ n
                return true
            end
        else
            rep = 1
        end
        last = x
    end

    return false
end

"""
    ispalindromic(seq::BioSequence)

Return `true` if `seq` is a palindromic sequence; otherwise return `false`.
"""
function ispalindromic(seq::BioSequence{<:NucleicAcidAlphabet})
    for i in 1:cld(length(seq), 2)
        if seq[i] != complement(seq[end - i + 1])
            return false
        end
    end

    return true
end

"""
    hasambiguity(seq::BioSequence)

Returns `true` if `seq` has an ambiguous symbol; otherwise return `false`.
"""
function hasambiguity(seq::BioSequence)
    for x in seq
        if isambiguous(x)
            return true
        end
    end
    return false
end
# 2 Bit specialization:
@inline hasambiguity(seq::BioSequence{<:NucleicAcidAlphabet{2}}) = false

"""
    iscanonical(seq::NucleotideSeq)

Returns `true` if `seq` is canonical.

For any sequence, there is a reverse complement, which is the same sequence, but
on the complimentary strand of DNA:

```
------->
ATCGATCG
CGATCGAT
<-------
```

!!! note
    Using the [`reverse_complement`](@ref) of a DNA sequence will give give this
    reverse complement.

Of the two sequences, the *canonical* of the two sequences is the lesser of the
two i.e. `canonical_seq < other_seq`.
"""
function iscanonical(seq::NucleotideSeq)
    i = 1
    j = lastindex(seq)
    @inbounds while i <= j
        f = seq[i]
        r = complement(seq[j])
        f < r && return true
        r < f && return false
        i += 1
        j -= 1
    end
    return true
end


text:
"""
    Timeseries(dataSource::DataSource, intervals::Vector{<:Interval}, granularity::Granularity; <keyword arguments>)

Returns a Druid.TimeseriesResult which is compatible with the Tables.jl interface.

A method with all arguments as keyword arguments is also provided.

# Arguments
- filter::Filter = nothing
- aggregations::Vector{<:Aggregator} = nothing
- postAggregations::Vector{<:PostAggregator} = nothing
- virtualColumns::Vector{<:VirtualColumn} = nothing
- descending::Bool = nothing
- limit::Integer = nothing
- context::Dict = nothing
"""
mutable struct Timeseries <: Query
    queryType::String
    dataSource::DataSource
    intervals::Vector{<:Interval}
    granularity::Granularity
    filter
    aggregations
    postAggregations
    virtualColumns
    descending
    limit
    context
    function Timeseries(
        dataSource, intervals, granularity;
        filter=nothing, aggregations=nothing, postAggregations=nothing,
        virtualColumns=nothing, descending=nothing, limit=nothing, context=nothing
    )
        nothing_or_type(filter, Filter)
        nothing_or_type(aggregations, Vector{<:Aggregator})
        nothing_or_type(postAggregations, Vector{<:PostAggregator})
        nothing_or_type(virtualColumns, Vector{<:VirtualColumn})
        nothing_or_type(descending, Bool)
        limit === nothing || (isa(limit, Integer) && limit >= 0) || error("limit must be a non-negative integer")
        nothing_or_type(context, Dict)
        new(
            "timeseries", dataSource, intervals, granularity,
            filter, aggregations, postAggregations,
            virtualColumns, descending, limit, context
        )
    end
end
Timeseries(
    ; dataSource, intervals, granularity,
    filter=nothing, aggregations=nothing, postAggregations=nothing,
    virtualColumns=nothing, descending=nothing, limit=nothing, context=nothing
) = Timeseries(
    dataSource, intervals, granularity;
    filter, aggregations, postAggregations,
    virtualColumns, descending, limit, context
)

function execute(client::Client, query::Timeseries; pretty=false)
    res = convert(Vector{Dict}, JSON.parse(execute_native_query(client, query; pretty)))
    if length(res) != 0
        names = [:timestamp, Symbol.(keys(res[1]["result"]))...]
    else
        names = Symbol[]
    end
    TimeseriesResult(names, res)
end

struct TimeseriesResult <: QueryResult
    names::Vector{Symbol}
    inner_array::Vector{Dict}
end

names(tr::TimeseriesResult) = getfield(tr, :names)

Base.summary(io::IO, tr::TimeseriesResult) = print(io, string(length(tr)) * "-element " * string(typeof(tr)))
Base.getindex(tr::TimeseriesResult, i::Int) = (i <= length(tr) || throw(BoundsError(tr, i))) && TimeseriesRow(i, tr)

Tables.rowaccess(::TimeseriesResult) = true
Tables.rows(tr::TimeseriesResult) = tr

Base.eltype(::TimeseriesResult) = TimeseriesRow
Base.length(tr::TimeseriesResult) = length(getfield(tr, :inner_array))
Base.iterate(tr::TimeseriesResult, state=1) = state > length(tr) ? nothing : (TimeseriesRow(state, tr), state + 1)

function Tables.schema(tr::TimeseriesResult)
    if length(tr) != 0
        types = [typeof(cl) for cl in tr[1]]
    else
        types = DataType[]
    end
    Tables.Schema(names(tr), types)
end

struct TimeseriesRow <: Tables.AbstractRow
    row::Int
    source::TimeseriesResult
end

function Tables.getcolumn(tr::TimeseriesRow, name::Symbol)
    if name == :timestamp
        getfield(getfield(tr, :source), :inner_array)[getfield(tr, :row)][string(name)]
    else
        getfield(getfield(tr, :source), :inner_array)[getfield(tr, :row)]["result"][string(name)]
    end
end

Tables.getcolumn(tr::TimeseriesRow, i::Int) = Tables.getcolumn(tr, names(getfield(tr, :source))[i])
Tables.columnnames(tr::TimeseriesRow) = names(getfield(tr, :source))


text:
# Windows GDI version
function getpixelcolors(x, y)
    hdc = ccall((:GetDC, "user32.dll"), UInt32, (UInt32,), 0)
    pxl = ccall((:GetPixel, "gdi32.dll"), UInt32, (UInt32, Cint, Cint), hdc, Cint(x), Cint(y))
    return pxl & 0xff, (pxl >> 8) & 0xff, (pxl >> 16) & 0xff
end

const x = 120
const y = 100
cols = getpixelcolors(x, y)
println("At screen point (x=$x, y=$y) the color RGB components are red: $(cols[1]), green: $(cols[2]), and blue: $(cols[3])")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# project.jl

using BenchmarkTools: @btime
using SPECTrecon: SPECTplan
using SPECTrecon: project!
using MATLAB

function call_SPECTproj_matlab(mpath, image, mumap, psfs, dy)

    mat"""
    addpath($mpath)
    SPECTproj_matlab($image, $mumap, $psfs, $dy);
    """
end

function project_time()
    T = Float32
    nx = 64
    ny = 64
    nz = 40
    nview = 60

    mumap = rand(T, nx, ny, nz)

    nx_psf = 19
    nz_psf = 19
    psfs = rand(T, nx_psf, nz_psf, ny, nview)
    psfs = psfs .+ mapslices(reverse, psfs, dims = [1, 2])
    psfs = psfs .+ mapslices(transpose, psfs, dims = [1, 2])
    psfs = psfs ./ mapslices(sum, psfs, dims = [1, 2])

    dy = T(4.7952)

    for interpmeth in (:one, :two)
        for mode in (:fast, :mem)
            plan = SPECTplan(mumap, psfs, dy; interpmeth, mode)
            xtrue = rand(T, nx, ny, nz)
            views = zeros(T, nx, nz, nview)
            println(string(interpmeth)*", "*string(mode))
            @btime project!($views, $xtrue, $plan)
        end
    end

    mpath = pwd()
    xtrue = rand(T, nx, ny, nz)
    println("project-matlab")
    println("Warning: Check if MIRT is installed")
    call_SPECTproj_matlab(mpath, xtrue, mumap, psfs, dy) # 216.518 ms, about 0.01 GiB
    nothing
end

# run all functions, time may vary on different machines, will allocate ~4MB memory.
project_time()

#= one, fast
    378.348 ms (25983 allocations: 1.37 MiB)
one, mem
    793.918 ms (31257 allocations: 1.76 MiB)
two, fast
    289.804 ms (25982 allocations: 1.37 MiB)
two, mem
    618.369 ms (31268 allocations: 1.76 MiB)
MIRT
    230.125 ms
=#