text:
function CompleteObservationFunction(PreObservationFunction::Function)
    numargs = MaximalNumberOfArguments(PreObservationFunction)
    if numargs == 1
        return (u::AbstractArray{<:Number}, t::Real, θ::AbstractVector{<:Number}) -> PreObservationFunction(u)
    elseif numargs == 2
        return (u::AbstractArray{<:Number}, t::Real, θ::AbstractVector{<:Number}) -> PreObservationFunction(u, t)
    elseif numargs == 3
        return (u::AbstractArray{<:Number}, t::Real, θ::AbstractVector{<:Number}) -> PreObservationFunction(u, t, θ)
    else
        throw("Given ObservationFunction should accept either arguments (u) or (u,t) or (u,t,θ). Got function which accepts $numargs arguments.")
    end
end

"""
    InformNames(DS::AbstractDataSet, sys::ODESystem, observables::AbstractVector{<:Int})
Copy the state names saved in `ODESystem` to `DS`.
"""
function InformNames(DS::AbstractDataSet, sys::ModelingToolkit.AbstractSystem, observables::Union{BoolArray,AbstractVector{<:Int}})
    newxnames = xnames(DS) == CreateSymbolNames(xdim(DS),"x") ? [string(ModelingToolkit.get_iv(sys))] : xnames(DS)
    newynames = ynames(DS) == CreateSymbolNames(ydim(DS),"y") ? string.(ModelingToolkit.get_states(sys)[observables]) : ynames(DS)
    InformNames(DS, newxnames, newynames)
end

# No ObservationFunction, therefore try to use sys to infer state names of ODEsys
# Extend for other DEFunctions in the future
function DataModel(DS::AbstractDataSet, sys::Union{ModelingToolkit.AbstractSystem,SciMLBase.AbstractDiffEqFunction}, u0::Union{AbstractArray{<:Number},Function},
                        observables::Union{AbstractVector{<:Int},BoolArray,Function}=collect(1:length(u0)), args...; tol::Real=1e-7, Domain::Union{HyperCube,Nothing}=nothing, kwargs...)
    newDS = (typeof(observables) <: AbstractVector{<:Int} && sys isa ModelingToolkit.AbstractSystem) ? InformNames(DS, sys, observables) : DS
    DataModel(newDS, GetModel(sys, u0, observables; tol=tol, Domain=Domain, kwargs...), args...)
end

# """
# Given `func` is converted to `ODEFunction`. Will probably deprecate this in the future.
# """
# function GetModel(func::Function, u0::Union{AbstractArray{<:Number},Function}, observables::Union{Function,AbstractVector{<:Int},BoolArray}=collect(1:length(u0)); tol::Real=1e-7,
#                     meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), Domain::Union{HyperCube,Nothing}=nothing, inplace::Bool=true)
#     GetModel(ODEFunction{inplace}(func), u0, observables; tol=tol, Domain=Domain, meth=meth, inplace=inplace)
# end

function GetModel(sys::ModelingToolkit.AbstractSystem, u0::Union{AbstractArray{<:Number},Function}, observables::Union{AbstractVector{<:Int},BoolArray,Function}=collect(1:length(u0));
                Domain::Union{HyperCube,Nothing}=nothing, inplace::Bool=true, kwargs...)
    # Is there some optimization that can be applied here? Modellingtoolkitize(sys) or something?
    # sys = Sys isa Catalyst.ReactionSystem ? convert(ODESystem, Sys) : Sys
    Model = if sys isa ModelingToolkit.AbstractODESystem
        GetModel(ODEFunction{inplace}(sys), u0, observables; Domain=Domain, inplace=inplace, kwargs...)
    else
        throw("Not programmed for $(typeof(sys)) yet, please convert to a ModelingToolkit.AbstractODESystem first.")
    end
    if Model isa ModelMap       Model = Model.Map    end
    pnames = ModelingToolkit.get_ps(sys) .|> string
    ylen = if observables isa Function      # ObservationFunction
        # Might still fail if states u are a Matrix.
        argnum = MaximalNumberOfArguments(observables)
        F = if argnum==1  z->observables(z)   elseif argnum==2  z->observables(z,0.1)
            elseif argnum==3 z->observables(z,0.1,GetStartP(length(pnames)))    else throw("Error") end
        num = GetArgLength(F)
        length(F(ones(num)))
    else
        observables isa BoolArray ? sum(observables) : length(observables)
    end
    plen = if Domain isa HyperCube
        length(Domain)
    elseif u0 isa AbstractArray     # Vector / Matrix
        # initial conditions given as array means the parameters are only the ps in sys
        length(pnames)
    else        # SplitterFunction
        # May well fail depending on how splitter function is implemented
        GetArgLength(u0)
    end
    xyp = (1, ylen, plen)
    Domain = isa(Domain, Bool) ? FullDomain(xyp[3], 1e5) : Domain

    pnames = plen - length(pnames) > 0 ? vcat(CreateSymbolNames(plen - length(pnames), "u"), pnames) : pnames
    # new(Map, InDomain, Domain, xyp, pnames, StaticOutput, inplace, CustomEmbedding)
    ModelMap(Model, nothing, Domain, xyp, pnames, Val(false), Val(false), Val(true))
end

# Although specialized methods for constant initial condition and specification of observed components in terms of arrays is faster than using functions, also use general method here.
function GetModel(func::SciMLBase.AbstractDiffEqFunction{T}, Initial, Observables; kwargs...) where T
    SplitterFunction = if Initial isa Function
        Initial
    else
        u0 = PromoteStatic(Initial, T)
        θ -> (u0, θ)
    end
    GetModel(func, SplitterFunction, (Observables isa Function ? Observables : u -> u[Observables]); kwargs...)
end

function GetModel(func::SciMLBase.AbstractDiffEqFunction, SplitterFunction::Function, ObservationFunction::Function; kwargs...)
    throw("If you see this error, it is most likely because no specialized method has been implemented for $(typeof(func)) yet.")
end

# Vanilla version with constant array of initial conditions and vector of observables.
function GetModel(func::AbstractODEFunction{T}, u0::AbstractArray{<:Number}, observables::Union{AbstractVector{<:Int},BoolArray}=collect(1:length(u0)); tol::Real=1e-7,
                    meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), Domain::Union{HyperCube,Nothing}=nothing, inplace::Bool=true) where T
    @assert T == inplace
    u0 = PromoteStatic(u0, inplace)

    function GetSol(θ::AbstractVector{<:Number}, u0::AbstractArray{<:Number}; tol::Real=tol, max_t::Number=10., meth::OrdinaryDiffEqAlgorithm=meth, kwargs...)
        odeprob = ODEProblem(func, u0, (0., max_t), θ)
        solve(odeprob, meth; reltol=tol, abstol=tol, kwargs...)
    end
    function ODEmodel(t::Number, θ::AbstractVector{<:Number}; observables::Union{AbstractVector{<:Int},BoolArray}=observables, u0::AbstractArray{<:Number}=u0,
                                                        tol::Real=tol, max_t::Number=t, meth::OrdinaryDiffEqAlgorithm=meth, FullSol::Bool=false, kwargs...)
        FullSol && return GetSol(θ, u0; tol=tol, max_t=max_t, meth=meth, kwargs...)
        sol = GetSol(θ, u0; tol=tol, max_t=max_t, meth=meth, save_everystep=false, save_start=false, save_end=true, kwargs...)
        sol.u[end][observables]
    end
    function ODEmodel(ts::AbstractVector{<:Number}, θ::AbstractVector{<:Number}; observables::Union{AbstractVector{<:Int},BoolArray}=observables, u0::AbstractArray{<:Number}=u0,
                                                tol::Real=tol, max_t::Number=maximum(ts), meth::OrdinaryDiffEqAlgorithm=meth, FullSol::Bool=false, kwargs...)
        FullSol && return GetSol(θ, u0; tol=tol, max_t=max_t, meth=meth, tstops=ts, kwargs...)
        sol = GetSol(θ, u0; tol=tol, max_t=max_t, meth=meth, saveat=ts, kwargs...)
        length(sol.u) != length(ts) && throw("ODE integration failed, maybe try using a lower tolerance value. θ=$θ.")
        [sol.u[i][observables] for i in 1:length(ts)] |> Reduction
    end
    MakeCustom(ODEmodel, Domain)
end

"""
    GetModel(func::ODEFunction, u0::AbstractArray, ObservationFunction::Function; tol::Real=1e-7, meth::OrdinaryDiffEqAlgorithm=Tsit5(), Domain::Union{HyperCube,Nothing}=nothing, inplace::Bool=true)
Returns a `ModelMap` which evolves the given system of ODEs from the initial configuration `u0` and afterwards applies the `ObservationFunction` to produce its predictions.

`ObservationFunction` should either be of the form `F(u) -> Vector` or `F(u,t) -> Vector` or `F(u,t,θ) -> Vector`.
Internally, the `ObservationFunction` is automatically wrapped as `F(u,t,θ)` if it is not already defined to accept three arguments.

A `Domain` can be supplied to constrain the parameters of the model to particular ranges which can be helpful in the fitting process.
"""
function GetModel(func::AbstractODEFunction{T}, u0::AbstractArray{<:Number}, PreObservationFunction::Function; tol::Real=1e-7,
                    meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), Domain::Union{HyperCube,Nothing}=nothing, inplace::Bool=true) where T
    @assert T == inplace
    u0 = PromoteStatic(u0, inplace)
    ObservationFunction = CompleteObservationFunction(PreObservationFunction)

    function GetSol(θ::AbstractVector{<:Number}, u0::AbstractArray{<:Number}; tol::Real=tol, max_t::Number=10., meth::OrdinaryDiffEqAlgorithm=meth, kwargs...)
        odeprob = ODEProblem(func, u0, (0., max_t), θ)
        solve(odeprob, meth; reltol=tol, abstol=tol, kwargs...)
    end
    function ODEmodel(t::Number, θ::AbstractVector{<:Number}; ObservationFunction::Function=ObservationFunction, u0::AbstractArray{<:Number}=u0,
                                                tol::Real=tol, max_t::Number=t, meth::OrdinaryDiffEqAlgorithm=meth, FullSol::Bool=false, kwargs...)
        FullSol && return GetSol(θ, u0; tol=tol, max_t=max_t, meth=meth, kwargs...)
        sol = GetSol(θ, u0; tol=tol, max_t=max_t, meth=meth, save_everystep=false, save_start=false, save_end=true, kwargs...)
        ObservationFunction(sol.u[end], t, θ)
    end
    function ODEmodel(ts::AbstractVector{<:Number}, θ::AbstractVector{<:Number}; ObservationFunction::Function=ObservationFunction, u0::AbstractArray{<:Number}=u0,
                                            tol::Real=tol, max_t::Number=maximum(ts), meth::OrdinaryDiffEqAlgorithm=meth, FullSol::Bool=false, kwargs...)
        FullSol && return GetSol(θ, u0; tol=tol, max_t=max_t, meth=meth, tstops=ts, kwargs...)
        sol = GetSol(θ, u0; tol=tol, max_t=max_t, meth=meth, saveat=ts, kwargs...)
        length(sol.u) != length(ts) && throw("ODE integration failed, maybe try using a lower tolerance value. θ=$θ.")
        [ObservationFunction(sol.u[i], sol.t[i], θ) for i in 1:length(ts)] |> Reduction
    end
    MakeCustom(ODEmodel, Domain)
end

"""
    GetModel(func::ODEFunction, SplitterFunction::Function, observables::Union{AbstractVector{<:Int},BoolArray}=collect(1:length(u0)); tol::Real=1e-7, meth::OrdinaryDiffEqAlgorithm=Tsit5(), Domain::Union{HyperCube,Nothing}=nothing, inplace::Bool=true)
Returns a `ModelMap` which evolves the given system of ODEs and returns `u[observables]` to produce its predictions.
Here, the initial conditions for the ODEs are produced from the parameters `θ` using the `SplitterFunction` which for instance allows one to estimate them from data.

`SplitterFunction` should be of the form `F(θ) -> (u0, p)`, i.e. the output is a `Tuple` whose first entry is the initial condition for the ODE model and the second entry constitutes the parameters which go on to enter the `ODEFunction`.
Typically, a fair bit of performance can be gained from ensuring that `SplitterFunction` outputs the initial condition `u0` as type `MVector` or `MArray`, if it has less than ~100 components.

A `Domain` can be supplied to constrain the parameters of the model to particular ranges which can be helpful in the fitting process.
"""
function GetModel(func::AbstractODEFunction{T}, SplitterFunction::Function, observables::Union{AbstractVector{<:Int},BoolArray}=collect(1:length(u0)); tol::Real=1e-7,
                    meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), Domain::Union{HyperCube,Nothing}=nothing, inplace::Bool=true) where T
    @assert T == inplace

    function GetSol(θ::AbstractVector{<:Number}, SplitterFunction::Function; tol::Real=tol, max_t::Number=10., meth::OrdinaryDiffEqAlgorithm=meth, kwargs...)
        u0, p = SplitterFunction(θ);        odeprob = ODEProblem(func, u0, (0., max_t), p)
        solve(odeprob, meth; reltol=tol, abstol=tol, kwargs...)
    end
    function ODEmodel(t::Number, θ::AbstractVector{<:Number}; observables::Union{AbstractVector{<:Int},BoolArray}=observables, SplitterFunction::Function=SplitterFunction,
                                tol::Real=tol, max_t::Number=t, meth::OrdinaryDiffEqAlgorithm=meth, FullSol::Bool=false, kwargs...)
        FullSol && return GetSol(θ, SplitterFunction; tol=tol, max_t=max_t, meth=meth, kwargs...)
        sol = GetSol(θ, SplitterFunction; tol=tol, max_t=max_t, meth=meth, save_everystep=false, save_start=false, save_end=true, kwargs...)
        sol.u[end][observables]
    end
    function ODEmodel(ts::AbstractVector{<:Number}, θ::AbstractVector{<:Number}; observables::Union{AbstractVector{<:Int},BoolArray}=observables, SplitterFunction::Function=SplitterFunction,
                                tol::Real=tol, max_t::Number=maximum(ts), meth::OrdinaryDiffEqAlgorithm=meth, FullSol::Bool=false, kwargs...)
        FullSol && return GetSol(θ, SplitterFunction; tol=tol, max_t=max_t, meth=meth, tstops=ts, kwargs...)
        sol = GetSol(θ, SplitterFunction; tol=tol, max_t=max_t, meth=meth, saveat=ts, kwargs...)
        length(sol.u) != length(ts) && throw("ODE integration failed, maybe try using a lower tolerance value. θ=$θ.")
        [sol.u[i][observables] for i in 1:length(ts)] |> Reduction
    end
    MakeCustom(ODEmodel, Domain)
end

"""
    GetModel(func::AbstractODEFunction{T}, SplitterFunction::Function, PreObservationFunction::Function; tol::Real=1e-7, meth::OrdinaryDiffEqAlgorithm=Tsit5(), Domain::Union{HyperCube,Nothing}=nothing, inplace::Bool=true)
Returns a `ModelMap` which evolves the given system of ODEs and afterwards applies the `ObservationFunction` to produce its predictions.
Here, the initial conditions for the ODEs are produced from the parameters `θ` using the `SplitterFunction` which for instance allows one to estimate them from data.

`SplitterFunction` should be of the form `F(θ) -> (u0, p)`, i.e. the output is a `Tuple` whose first entry is the initial condition for the ODE model and the second entry constitutes the parameters which go on to enter the `ODEFunction`.
Typically, a fair bit of additional performance can be gained from ensuring that `SplitterFunction` outputs the initial condition `u0` as type `MVector` or `MArray`, if it has less than ~100 components.

`ObservationFunction` should either be of the form `F(u) -> Vector` or `F(u,t) -> Vector` or `F(u,t,θ) -> Vector`.
Internally, the `ObservationFunction` is automatically wrapped as `F(u,t,θ)` if it is not already defined to accept three arguments.

!!! note
    The vector `θ` passed to `ObservationFunction` is the same `θ` that is passed to `SplitterFunction`, i.e. before splitting.
    This is because `ObservationFunction` might also depend on the initial conditions in general.

A `Domain` can be supplied to constrain the parameters of the model to particular ranges which can be helpful in the fitting process.
"""
function GetModel(func::AbstractODEFunction{T}, SplitterFunction::Function, PreObservationFunction::Function; tol::Real=1e-7,
                    meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), Domain::Union{HyperCube,Nothing}=nothing, inplace::Bool=true) where T
    @assert T == inplace
    ObservationFunction = CompleteObservationFunction(PreObservationFunction)

    function GetSol(θ::AbstractVector{<:Number}, SplitterFunction::Function; tol::Real=tol, max_t::Number=10., meth::OrdinaryDiffEqAlgorithm=meth, kwargs...)
        u0, p = SplitterFunction(θ);        odeprob = ODEProblem(func, u0, (0., max_t), p)
        solve(odeprob, meth; reltol=tol, abstol=tol, kwargs...)
    end
    function ODEmodel(t::Number, θ::AbstractVector{<:Number}; ObservationFunction::Function=ObservationFunction, SplitterFunction::Function=SplitterFunction,
                                                    tol::Real=tol, max_t::Number=t, meth::OrdinaryDiffEqAlgorithm=meth, FullSol::Bool=false, kwargs...)
        FullSol && return GetSol(θ, SplitterFunction; tol=tol, max_t=max_t, meth=meth, kwargs...)
        sol = GetSol(θ, SplitterFunction; tol=tol, max_t=max_t, meth=meth, save_everystep=false, save_start=false, save_end=true, kwargs...)
        ObservationFunction(sol.u[end], sol.t[end], θ)
    end
    function ODEmodel(ts::AbstractVector{<:Number}, θ::AbstractVector{<:Number}; ObservationFunction::Function=ObservationFunction, SplitterFunction::Function=SplitterFunction,
                                            tol::Real=tol, max_t::Number=maximum(ts), meth::OrdinaryDiffEqAlgorithm=meth, FullSol::Bool=false, kwargs...)
        FullSol && return GetSol(θ, SplitterFunction; tol=tol, max_t=max_t, meth=meth, tstops=ts, kwargs...)
        sol = GetSol(θ, SplitterFunction; tol=tol, max_t=max_t, meth=meth, saveat=ts, kwargs...)
        length(sol.u) != length(ts) && throw("ODE integration failed, maybe try using a lower tolerance value. θ=$θ.")
        [ObservationFunction(sol.u[i], sol.t[i], θ) for i in 1:length(ts)] |> Reduction
    end
    MakeCustom(ODEmodel, Domain)
end

"""
    ModifyODEmodel(DM::AbstractDataModel, NewObservationFunc::Function) -> ModelMap
Constructs a new `ModelMap` with new observation function `f(u,t,θ)` from a given ODE-based `DataModel`.
"""
ModifyODEmodel(DM::AbstractDataModel, NewObservationFunc::Function) = ModifyODEmodel(DM, Predictor(DM), NewObservationFunc)
function ModifyODEmodel(DM::AbstractDataModel, Model::ModelMap, NewObservationFunc::Function)
    # Model.Map isa DEModel && return ModifyDEModel(Model, NewObservationFunc)
    Eval = try
        Model(WoundX(DM)[1], MLE(DM); FullSol=true).u[1]
    catch;
        throw("It appears as though the given model is not an ODEmodel.")
    end
    F = CompleteObservationFunction(NewObservationFunc)
    out = F(Eval, WoundX(DM)[1], MLE(DM))
    function NewODEmodel(x::AbstractVector{<:Number}, θ::AbstractVector{<:Number}; FullSol=false, kwargs...)
        FullSol && return Model.Map(x, θ; FullSol=true, kwargs...)
        sol = Model.Map(x, θ; FullSol=true, saveat=x, kwargs...)
        length(sol.u) != length(x) && throw("ODE integration failed, maybe try using a lower tolerance value. θ=$θ.")
        [F(sol.u[i], sol.t[i], θ) for i in 1:length(x)] |> Reduction
    end
    function NewODEmodel(x::Number, θ::AbstractVector{<:Number}; FullSol=false, kwargs...)
        FullSol && return Model.Map(x, θ; FullSol=true, kwargs...)
        sol = Model.Map(x, θ; FullSol=true, save_everystep=false, save_start=false, save_end=true, kwargs...)
        F(sol.u[end], sol.t[end], θ)
    end
    ModelMap(NewODEmodel, Model.InDomain, Model.Domain, (Model.xyp[1], length(out), Model.xyp[3]), Model.pnames, Val(out isa SVector), Model.inplace, Model.CustomEmbedding)
end


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

# usage: perfcomp.jl <baseline file> [<suite>]
# This will run the specified suite (default "all") and compare it to stored
# results in the baseline file. Only test names present in both will be
# compared.
# The file format is the output of running `make` in this directory.

function readperf(f)
    Dict(rstrip(l[1:19])=>[parse(Float64,l[20:27]),parse(Float64,l[29:36]),parse(Float64,l[38:45]),parse(Float64,l[47:54])] for l in eachline(f))
end

function main()
    baseline = readperf(open(ARGS[1]))
    torun = length(ARGS) > 1 ? ARGS[2] : "all"
    e = haskey(ENV,"J") ? "JULIA_EXECUTABLE=$(ENV["J"])" : ""
    io,p = open(`make $e -s $torun`, "r")
    newp = readperf(io)

    names = sort(intersect(keys(baseline),keys(newp)))

    means0 = [ baseline[n][3] for n in names ]
    means1 = [ newp[n][3] for n in names ]

    change = (means0 - means1) ./ means0

    println("test name            old       new     % speedup  % st. dev")
    println("-----------------------------------------------------------")
    for i = 1:length(names)
        print(rpad(names[i],19))
        print(lpad(string(means0[i]),8),"  ")
        print(lpad(string(means1[i]),8),"  ")
        dev = baseline[names[i]][4]/means0[i]
        @printf "%7.2f%%  %7.2f%%\n" change[i]*100 dev*100
    end
    println()

    minname = names[argmin(change)]
    maxname = names[argmax(change)]

    minstd = baseline[minname][4]/baseline[minname][3]
    maxstd = baseline[maxname][4]/baseline[maxname][3]

    @printf "min %7.2f%% (std %7.2f%%) %s\n" minimum(change)*100 minstd*100 minname
    @printf "max %7.2f%% (std %7.2f%%) %s\n" maximum(change)*100 maxstd*100 maxname
    @printf "avg %7.2f%%\n" mean(change)*100
end

main()

ext: jl

max_stars_repo_name: JuliaTagBot/Benchmarks

max_issues_repo_name: JuliaArchive/Benchmarks

max_forks_repo_name: JuliaArchive/Benchmarks


text:
using HuberLR
using Base.Test

let
    # true parameters
    w = [1., 2, 3., 4.]
    b = 5.

    # generate sample points
    srand(20150218)
    x = randn(length(w), 150)
    y = mapslices(x, 1) do x
        w ⋅ x + b + randn() + (rand() > .9 ? abs(randn()) * 5 : 0.)
    end |> vec
    #writedlm("data.tsv", vcat(x, y')')

    # JIT compile
    lr = LinearRegression(.5)
    fit!(lr, x[:,1:5], y[1:5])

    # robust regression
    lr = LinearRegression(1.)
    @time fit!(lr, x, y)
    println(lr)
    @show norm(predict(lr, x) .- y)
    Δw = norm(lr.w .- [b, w])

    # quadratic loss (not robust)
    lr = LinearRegression(Inf)
    @time fit!(lr, x, y)
    println(lr)
    @show norm(predict(lr, x) .- y)
    Δw′ = norm(lr.w .- [b, w])

    #@show Δw Δw′
    @test Δw < Δw′ < 1.
end


text:
# # Maze Solver
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../maze.mp4" type="video/mp4">
# </video>
# ```
# Consider a scenario where a walker agent is stuck in a maze. Finding the shortest path through an
# arbitrary maze or map is simulated using a [`Pathfinding.Pathfinder`](@ref) and its `walkable` map property.

# ## Setup
using Agents, Agents.Pathfinding
using FileIO # To load images you also need ImageMagick available to your project

# The `Walker` agent needs no special property, just the `id` and `position` from [`@agent`](@ref).
@agent Walker GridAgent{2} begin end

# The maze is stored as a simple .bmp image, where each pixel corresponds to a position on the grid.
# White pixels correspond to walkable regions of the maze.
function initalize_model(map_url)
    ## Load the maze from the image file. White values can be identified by a
    ## non-zero red component
    maze = BitArray(map(x -> x.r > 0, load(download(map_url))))
    ## Create a pathfinder by specifying the `walkable` parameter for the pathfinder.
    ## Since we are interested in the most direct path to the end, the default
    ## `DirectDistance` is appropriate.
    ## `diagonal_movement` is set to false to prevent cutting corners by going along
    ## diagonals.
    pathfinder = Pathfinder(walkable=maze, diagonal_movement=false)
    ## The size of the space is the size of the maze
    space = GridSpace(size(maze); pathfinder, periodic = false)
    model = ABM(Walker, space)
    ## Place a walker at the start of the maze
    walker = Walker(1, (1, 4))
    add_agent_pos!(walker, model)
    ## The walker's movement target is the end of the maze
    set_target!(walker, (41, 32), model)

    return model
end

# # Dynamics
# Stepping the agent is a trivial matter of calling [`move_along_route!`](@ref) to move it along it's path to
# the target.
agent_step!(agent, model) = move_along_route!(agent, model)

# ## Visualization
# Visualizing the `Walker` move through the maze is handled through [`InteractiveDynamics.abm_plot`](@ref).
using InteractiveDynamics
using CairoMakie
CairoMakie.activate!(type = "png") # hide

## Our sample walkmap
map_url =
    "https://raw.githubusercontent.com/JuliaDynamics/" *
    "JuliaDynamics/master/videos/agents/maze.bmp"
model = initalize_model(map_url)

# The `heatarray` keyword argument allows plotting the maze as a heatmap behind the agent.
abm_video(
    "maze.mp4",
    model,
    agent_step!;
    resolution=(700,700),
    frames=310,
    framerate=30,
    ac=:red,
    as=11,
    heatarray = model -> walkmap(model),
    add_colorbar = false,
)
nothing # hide

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../maze.mp4" type="video/mp4">
# </video>
# ```


text:
using AnalyticalEOR
using Test

@testset "AnalyticalEOR.jl" begin

    # Test relative permeability functions
    @test water_rel_perm(0.5, 0.2, 0.2, 0.3, 3.0) ≈ 0.0374999999
    @test water_rel_perm([0.5, 0.5], 0.2, 0.2, 0.3, 3.0) ≈ Vector([0.0375, 0.0375])

    @test oil_rel_perm(0.5, 0.2, 0.2, 0.3, 3.0) ≈ 0.0375
    @test oil_rel_perm([0.5, 0.5], 0.2, 0.2, 0.3, 3.0) ≈ Vector([0.0375, 0.0375])

    @test fractional_flow([0.5, 0.6], 0.2, 0.2, 0.3, 1.0, 3.0, 3.0, 1.0, 1.0) ≈ Vector([ 0.2307692307692305,  0.7058823529411761])
end

@testset "Derivative Tests" begin

    @test krw_derivative([0.5, 0.5], 0.2, 0.2, 0.3, 3.0) ≈ [0.375, 0.375]
    @test kro_derivative([0.5, 0.5], 0.2, 0.2, 0.3, 3.0) ≈ [-0.375, -0.375]

    # test derivative of fractional flow
end


text:
# FOR BROADCASTING VARIABLES

Doable = Union{Symbol, Expr}

function yieldfirst(todo::Doable) #required while bug 14445 still around
    return :(yield(); $todo)
end

function sow(pid::Int64, name::Doable, value; mod=Main, callstyle=remotecall_wait)
    todo = Expr(:(=), name, value)
    #callstyle(Main.eval, pid, mod, yieldfirst(todo))
    callstyle(Main.eval, pid, mod, todo)
end

function sow(pids::Array{Int64, 1}, name::Doable, value; mod=Main, callstyle=remotecall_wait)
    refs = Dict{Int64, Any}()
    @sync for p in pids
        @async begin
            rref = sow(p, name, value; mod=mod, callstyle=callstyle)
            refs[p] = rref
        end
    end
    refs
end

function sow(name::Doable, value; mod=Main, callstyle=remotecall_wait)
    pids = workers()
    sow(pids, name, value; mod=mod, callstyle=callstyle)
end

# HARVESTING VALUES OR REFERENCES FROM REMOTES

function reap(pids::Array{Int64, 1}, doable::Doable; callstyle=remotecall_fetch, returntype=Any)
    results = Dict{Int64, returntype}()
    #todo = yieldfirst(doable)
    @sync for pid in pids
        @async results[pid] = callstyle(Main.eval, pid, doable)
    end
    results
end

function reap(doable::Doable; callstyle=remotecall_fetch, returntype=Any)
    pids = workers()
    reap(pids, doable; callstyle=callstyle, returntype=returntype)
end

function reap(pid::Int64, doable::Doable; callstyle=remotecall_fetch, returntype=Any)
    pids = [pid];
    reap(pids, doable; callstyle=callstyle, returntype=returntype)
end

function reaprefs(pids, doable::Doable; returntype=Any)
    reap(pids, doable; callstyle=remotecall_fetch, returntype=returntype)
end


text:
function bam2pileup2consensus(bam_filename, out_stream, ref_filename,
                              min_ac, min_af, max_dp)
    check_samtools_version()
    cmd = `samtools mpileup -f $ref_filename $bam_filename`
    open(cmd) do pileup_stream
        open(ref_filename) do ref_stream
            if endswith(ref_filename, ".gz")
                reader = FASTA.Reader(GzipDecompressorStream(ref_stream))
            else
                reader = FASTA.Reader(ref_stream)
            end

            consensus_fasta = pileup_consensus_fasta(
                pileup_stream, reader,
                min_ac, min_af, max_dp
            )

            w = FASTA.Writer(out_stream)
            for record in consensus_fasta
                write(w, record)
            end
        end
    end
end

function check_samtools_version()
    open(`samtools --version`) do cmd
        version = match(r"samtools (\d+)\.(\d+)", readline(cmd))
        major = parse(Float64, version.captures[1])
        minor = parse(Float64, version.captures[2])
        min_major = 1
        min_minor = 4
        recommended_version = "1.9"
        if major < min_major || (major == min_major && minor < min_minor)
            throw(ErrorException(string(
                "samtools version $version too low.",
                " samtools>=$min_version required;",
                " samtools>=$recommended_version recommended.")))
        elseif "$major.$minor" != recommended_version
            @warn "samtools=$major.$minor. samtools=$recommended_version recommended."
        end
    end
end

function pileup_consensus_fasta(pileup_stream, ref_fasta,
                                ac_lower, af_lower, dp_upper)
    consensus_seq_list = initialize_consensus_seq(ref_fasta)
    consensus_seq_dict = Dict{String, Array{Char}}(
        name => seq for (name, seq) in consensus_seq_list)

    chrom, pos = "", 0

    allele_array = Char[]
    allele_counter = Dict{Char, Int}(a => 0 for a in "ACGTN-")
    for line in eachline(pileup_stream)
        chrom, pos, ref, dp, read_bases, base_quals = parse_pileup_line(
            line, chrom, pos)

        consensus_seq_dict[chrom][pos] = consensus_allele!(
            allele_counter, allele_array,
            read_bases, ref,
            ac_lower, af_lower, dp_upper
        )
    end

    return [
        FASTA.Record(name, LongCharSeq(join(seq)))
        for (name, seq) in consensus_seq_list
    ]
end

function parse_pileup_line(line, prev_chrom, prev_pos)
    fields = split(line)

    chrom = fields[1]
    pos = parse(Int, fields[2])
    @assert chrom != prev_chrom || pos > prev_pos

    ref_idx = 3
    @assert length(fields[ref_idx]) == 1
    ref = fields[ref_idx][1]

    dp = parse(Int, fields[4])
    read_bases = fields[5]
    base_quals = fields[6]

    return (chrom, pos, ref, dp, read_bases, base_quals)
end

function initialize_consensus_seq(ref_fasta)
    consensus_seq = []
    for record in ref_fasta
        name = FASTA.identifier(record)
        seq = sequence(LongCharSeq, record)
        push!(consensus_seq, (name, ['N' for i in 1:length(seq)]))
    end
    return consensus_seq
end

function consensus_allele!(allele_counter, allele_array,
                           read_bases, ref, ac_lower, af_lower, dp_upper)
    fill_alleles_array!(allele_array, read_bases, ref)
    dp = count_alleles!(allele_counter, allele_array)
    major_allele, max_ac = max_allele_and_count(allele_counter)
    max_af = max_ac / dp

    if max_ac >= ac_lower && max_af >= af_lower && dp <= dp_upper
        return major_allele
    else
        return 'N'
    end
end

function ref_and_bases(pileup_line)
    return ref, read_bases
end

function count_alleles!(allele_counter, allele_array)
    for k in keys(allele_counter)
        allele_counter[k] = 0
    end
    dp = 0
    for a in allele_array
        if !haskey(allele_counter, a)
            a = 'N'
        end
        allele_counter[a] += 1
        dp += 1
    end
    return dp
end

function max_allele_and_count(allele_counter)
    major_allele = 'N'
    max_ac = 0
    for (a, ac) in allele_counter
        if ac > max_ac
            major_allele = a
            max_ac = ac
        end
    end
    return major_allele, max_ac
end

function fill_alleles_array!(allele_array, read_bases, ref)
    resize!(allele_array, 0)
    idx = 1
    while idx <= length(read_bases)
        idx = process_next_symbol!(allele_array, read_bases, idx, ref)
    end
end

function process_next_symbol!(allele_array, read_bases, idx, ref)
    next_char = read_bases[idx]
    if next_char == '+' || next_char == '-'
        return indel_bases_end(read_bases, idx) + 1
    elseif next_char == '^'
        return idx + 2
    elseif next_char == '$' || next_char == '<' || next_char == '>'
        return idx + 1
    else
        push!(allele_array, get_allele(next_char, ref))
        return idx + 1
    end
end

function indel_bases_end(read_bases, indel_start)
    digits_str = match(r"\d+", read_bases, indel_start).match
    digits_num = parse(Int, digits_str)
    return indel_start + length(digits_str) + digits_num
end

function get_allele(symbol, ref)
    if symbol == '.' || symbol == ','
        return ref
    elseif symbol == '*'
        return '-'
    else
        return uppercase(symbol)
    end
end

export pileup_consensus_fasta
export bam2pileup2consensus


text:
using  BenchmarkTools, StaticArrays
import NBodyIPs

test(f, x) = (s = 0.0; for n = 1:10; s += f(x)[1][1]; end; s)

@info "Runtimes of single invariants computation"
@info "Devide runtimes by 10"
for N in 2:4
   x = @SVector rand((N*(N-1))÷2)
   print("$N-body BL-invariants: ");
   @btime test(NBodyIPs.BLInvariants.invariants_ed, $x)
   print("$N-body BA-invariants: ");
   @btime test(NBodyIPs.BAInvariants.invariants_ed, $x)
end

x = @SVector rand(10)
@btime NBodyIPs.BLInvariants.invariants($x)
@btime NBodyIPs.BLInvariants.invariants_d($x)
@btime NBodyIPs.BLInvariants.invariants_ed($x)

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1ProjectedVolumeSource <: SwaggerModel
    defaultMode::Any # spec type: Union{ Nothing, Int32 } # spec name: defaultMode
    sources::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1VolumeProjection} } # spec name: sources

    function IoK8sApiCoreV1ProjectedVolumeSource(;defaultMode=nothing, sources=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ProjectedVolumeSource, Symbol("defaultMode"), defaultMode)
        setfield!(o, Symbol("defaultMode"), defaultMode)
        validate_property(IoK8sApiCoreV1ProjectedVolumeSource, Symbol("sources"), sources)
        setfield!(o, Symbol("sources"), sources)
        o
    end
end # type IoK8sApiCoreV1ProjectedVolumeSource

const _property_map_IoK8sApiCoreV1ProjectedVolumeSource = Dict{Symbol,Symbol}(Symbol("defaultMode")=>Symbol("defaultMode"), Symbol("sources")=>Symbol("sources"))
const _property_types_IoK8sApiCoreV1ProjectedVolumeSource = Dict{Symbol,String}(Symbol("defaultMode")=>"Int32", Symbol("sources")=>"Vector{IoK8sApiCoreV1VolumeProjection}")
Base.propertynames(::Type{ IoK8sApiCoreV1ProjectedVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1ProjectedVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1ProjectedVolumeSource }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1ProjectedVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ProjectedVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ProjectedVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1ProjectedVolumeSource)
    (getproperty(o, Symbol("sources")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ProjectedVolumeSource }, name::Symbol, val)
end


text:
using Enzyme
using LLVM
using Test

import Enzyme: typetree, TypeTree, API

const ctx = LLVM.Context()
const dl = string(LLVM.DataLayout(LLVM.JITTargetMachine()))

tt(T) = string(typetree(T, ctx, dl))

struct Atom
  x::Float32
  y::Float32
  z::Float32
  type::Int32
end

struct Composite
    k::Atom
    y::Atom
end

@testset "TypeTree" begin
    @test tt(Float16) == "{[-1]:Float@half}"
    @test tt(Float32) == "{[-1]:Float@float}"
    @test tt(Float64) == "{[-1]:Float@double}"
    @test tt(Symbol) == "{}"
    @test tt(String) ==  "{}"
    @test tt(AbstractChannel) == "{}"
    @test tt(Base.ImmutableDict{Symbol, Any}) == "{[0]:Pointer, [8]:Pointer, [16]:Pointer}"
    @test tt(Atom) == "{[0]:Float@float, [4]:Float@float, [8]:Float@float, [12]:Integer, [13]:Integer, [14]:Integer, [15]:Integer}"
    @test tt(Composite) == "{[0]:Float@float, [4]:Float@float, [8]:Float@float, [12]:Integer, [13]:Integer, [14]:Integer, [15]:Integer, [16]:Float@float, [20]:Float@float, [24]:Float@float, [28]:Integer, [29]:Integer, [30]:Integer, [31]:Integer}"
end


text:
module RevealExpr

using REPL
using REPL.LineEdit

using ReplMaker

function __init__()

  function parse_multi(src::AbstractString)::Union{Missing,Vector{Any}}
    x = Meta.parse("@reveal " * src)
    if Meta.isexpr(x, :incomplete)
      return missing
    end
    @assert Meta.isexpr(x, :macrocall)
    xs = x.args[2:end]
    if length(xs) >= 1 && xs[1] isa LineNumberNode
      xs = xs[2:end]
    end
    return xs
  end

  function dump_expr_by_src(s)
    xs = parse_multi(s)
    for x in xs
      dump(x, maxdepth=30) # TODO: un-hardcode `maxdepth`?
    end
    # admit interpolation
    return Expr(:quote, Expr(:block, xs...))
  end

  function is_complete_julia(s)
    input = String(take!(copy(LineEdit.buffer(s))))
    return parse_multi(input) !== missing
  end

  initrepl(dump_expr_by_src,
    prompt_text="Expr> ",
    prompt_color=:blue,
    start_key=')',
    mode_name="ExprReveal_mode",
    valid_input_checker=is_complete_julia,
  )
end

end # module

max_stars_repo_name: complyue/RevealExpr

max_issues_repo_name: complyue/RevealExpr

max_forks_repo_name: complyue/RevealExpr


text:
#!/bin/julia

#############################################
# Garamond script for TCP client operations #
#############################################
module garw

using Logging
using Sockets
using ArgParse
using HTTP

# Function that parses Garamond's web-socket client arguments
function get_web_socket_client_arguments(args::Vector{String})
	s = ArgParseSettings()
	@add_arg_table! s begin
        "--log-level"
            help = "logging level"
            default = "warn"
        "--web-socket-port", "-w"
            help = "WEB socket data communication port"
            arg_type = UInt16
            default = UInt16(0)
        "--web-socket-ip"
            help = "WEB socket data communication IP"
            default = "127.0.0.1"
        "--http-port", "-p"
            help = "HTTP port for the http server"
            arg_type = Int
            default = 8888
        "--web-page"
            help = "Search web page to serve"
            arg_type = String
        "--return-fields"
            help = "List of fields to return (ignores wrong names)"
            nargs = '*'
	end
	return parse_args(args,s)
end

stringvec(vv) = "[" * join(map(v->"\"" * string(v) * "\"", vv), ",") * "]"

# Function that returns the default webpage
function _default_webpage(ws_ip::String, ws_port::UInt16; fields=[])
    return "
	<!doctype html>

	    <title>~ Garamond Search ~ by 0x0a Research </title>

			<form id=\"search_box\">
				<input id=\"search_box_txt\" type=\"text\" class=\"search\" placeholder=\"Search...\" required>
				<input id=\"search_box_button\" type=\"button\" class=\"button\" value=\"Search\">

 powered by 0x0a Research

			</form>

    "
end

########################
# Main module function #
########################
function julia_main()::Cint
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function real_main()
    # Parse command line arguments
    args = get_web_socket_client_arguments(ARGS)

    # Logging
    log_levels = Dict("debug" => Logging.Debug,
                      "info" => Logging.Info,
                      "warning" => Logging.Warn,
                      "error" => Logging.Error)
    logger = ConsoleLogger(stdout,
                get(log_levels, lowercase(args["log-level"]), Logging.Info))
    global_logger(logger)

    # Start client
    @info "~ GARAMOND ~ (web-socket client)"
    ws_ip = args["web-socket-ip"]
    ws_port = args["web-socket-port"]
    return_fields = args["return-fields"]

    if ws_port > 0
        # Get web page
        webpage_file = args["web-page"]
        if webpage_file != nothing && isfile(webpage_file)
            webpage = read(webpage_file, String)
        else
            webpage = _default_webpage(ws_ip, ws_port; fields=return_fields)
        end

        # Start HTTP server
        http_port = args["http-port"]
        _handler(request::HTTP.Request) = begin
            try
                return HTTP.Response(200, webpage)
            catch e
                return HTTP.Response(404, "Error: $e")
            end
        end
        http_ip = Sockets.localhost
        @info "Serving page on $http_ip:$http_port"
        HTTP.serve(_handler, http_ip, http_port, readtimeout=0)
    else
        @warn "Wrong web-socket port value $ws_port (default is 0). Exiting..."
    end
    return 0
end

################################
# Start main Garamond function #
################################
if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end

end  # module


text:
const DynamicResult{T,E} = Union{Ok{T},Err{E}}

function _ConcreteResult end

struct ConcreteResult{T,E} <: AbstractResult{T,E}
    value::DynamicResult{T,E}

    global _ConcreteResult(::Type{T}, ::Type{E}, value) where {T,E} = new{T,E}(value)
end

# Is this mixture of invariance and covariance good?
const Result{T,E} = Union{ConcreteResult{T,E},DynamicResult{<:T,<:E}}

Try.unwrap(result::ConcreteResult) = Try.unwrap(result.value)
Try.unwrap_err(result::ConcreteResult) = Try.unwrap_err(result.value)
Try.isok(result::ConcreteResult) = result.value isa Ok
Try.iserr(result::ConcreteResult) = result.value isa Err

_concrete(result::Ok) = _ConcreteResult(Try.oktype(result), Union{}, result)
_concrete(result::Err) = _ConcreteResult(Union{}, Try.errtype(result), result)

Base.convert(::Type{ConcreteResult{T,E}}, result::Ok) where {T,E} =
    _ConcreteResult(T, E, convert(Ok{T}, result))
Base.convert(::Type{ConcreteResult{T}}, result::Ok) where {T} =
    _ConcreteResult(T, Union{}, convert(Ok{T}, result))

Base.convert(::Type{ConcreteResult{T,E}}, result::Err) where {T,E} =
    _ConcreteResult(T, E, convert(Err{E}, result))
Base.convert(::Type{ConcreteResult{<:Any,E}}, result::Err) where {E} =
    _ConcreteResult(Union{}, E, convert(Err{E}, result))

function Base.convert(
    ::Type{ConcreteResult{T,E}},
    result::ConcreteResult{T′,E′},
) where {T,E,T′<:T,E′<:E}
    value = result.value
    if value isa Ok
        return _ConcreteResult(T, E, Ok{T}(value.value))
    else
        return _ConcreteResult(T, E, Err{E}(value.value))
    end
end


text:
"""
    allfxq(iset, x, qmu2, n, ichk)

Get all flavour-pdf values for a given `x` and `mu^2`.

# Arguments
- `iset::Integer`: pdf set id (1-24).
- `x::Float64`: input value of `x`.
- `qmu2::Float64`: input value of `qmu2`.
- `n::Integer`: number of additional pdfs to be returned.
- `ichk::Integer`: flag to steer error checking. See QCDNUM
docs. `ichk = -1` makes code faster, at the risk of not
checking certain things.

# Returns
- `pdf::Array{Float64,2}`: pdf values over flavours.
"""
function allfxq(iset::Integer, x::Float64, qmu2::Float64,
                n::Integer, ichk::Integer)

    iset = Ref{Int32}(iset)
    x = Ref{Float64}(x)
    qmu2 = Ref{Float64}(qmu2)
    n = Ref{Int32}(n)
    ichk = Ref{Int32}(ichk)

    if n[] == 0
        pdf = Array{Float64}(undef, 13)
    else
        pdf = Array{Float64}(undef, 13, n+1)
    end

    @ccall allfxq_(iset::Ref{Int32}, x::Ref{Float64}, qmu2::Ref{Float64},
                   pdf::Ref{Float64}, n::Ref{Int32}, ichk::Ref{Int32})::Nothing

    pdf
end

"""
    sumfxq(iset, c, isel, x, qmu2, ichk)

Return the gluon or a weighted sum of quark
densities, depending on the selection flag, `isel`.

# `isel` values:
0: Gluon density |xg>
1: Linear combination `c` summed over _active_ flavours
2-8: Specific singlet/non-singlet quark component
9: Intrinsic heavy flavours
12+i: Additional pdf |xf_i> in iset

See QCDNUM manual for more information.

# Arguments
- `iset::Integer`: pdf set id (1-24).
- `c::Array{Float64}`: Coefficients of quarks/anti-quarks
- `isel::Integer`: Selection flag
- `x::Float64`: input value of `x`.
- `qmu2::Float64`: input value of `qmu2`.
- `n::Integer`: number of additional pdfs to be returned.
- `ichk::Integer`: flag to steer error checking. See QCDNUM
docs. `ichk = -1` makes code faster, at the risk of not
checking certain things.

# Returns
- `pdf::Float64`: pdf value.
"""
function sumfxq(iset::Integer, c::Array{Float64}, isel::Integer, x::Float64,
                qmu2::Float64, ichk::Integer)

    iset = Ref{Int32}(iset)
    isel = Ref{Int32}(isel)
    x = Ref{Float64}(x)
    qmu2 = Ref{Float64}(qmu2)
    ichk = Ref{Int32}(ichk)

    pdf = @ccall sumfxq_(iset::Ref{Int32}, c::Ref{Float64}, isel::Ref{Int32},
                         x::Ref{Float64}, qmu2::Ref{Float64},
                         ichk::Ref{Int32})::Float64

    pdf[]
end

"""
    bvalij(iset, id, ix, iq, ichk)

Get the value of a basis pdf at a given point (`ix`, `iq`).

# Arguments
- `iset::Integer`: pdf set id (1-24)
- `id::Integer`: basis pdf identifier from 0 to 12+n,
where n is the number of additional pdfs in iset.
- `ix::Integer`: x index.
- `iq::Integer`: qq index.
- `ichk::Integer`: flag to steer error checking.
See allfxq().

# Returns
- `pdf::Float64`: pdf values.
"""
function bvalij(iset::Integer, id::Integer, ix::Integer, iq::Integer, ichk::Integer)

    iset = Ref{Int32}(iset)
    id = Ref{Int32}(id)
    ix = Ref{Int32}(ix)
    iq = Ref{Int32}(iq)
    ichk = Ref{Int32}(ichk)

    pdf = @ccall bvalij_(iset::Ref{Int32}, id::Ref{Int32}, ix::Ref{Int32},
                         iq::Ref{Int32}, ichk::Ref{Int32})::Float64

    pdf[]
end

"""
    bvalxq(iset, id, x, qq, ichk)

Get the value of a basis pdf at a given point (`x`, `qq`).

# Arguments
- `iset::Integer`: pdf set id (1-24)
- `id::Integer`: basis pdf identifier from 0 to 12+n,
where n is the number of additional pdfs in iset.
- `x::Float64`: x value.
- `qq::Float64`: qq value.
- `ichk::Integer`: flag to steer error checking.
See allfxq().

# Returns
- `pdf::Float64`: pdf values.
"""
function bvalxq(iset::Integer, id::Integer, x::Float64, qq::Float64, ichk::Integer)

    iset = Ref{Int32}(iset)
    id = Ref{Int32}(id)
    x = Ref{Float64}(x)
    q = Ref{Float64}(qq)
    ichk = Ref{Int32}(ichk)

    pdf = @ccall bvalxq_(iset::Ref{Int32}, id::Ref{Int32}, x::Ref{Float64},
                         qq::Ref{Float64}, ichk::Ref{Int32})::Float64

    pdf[]
end

"""
    fvalij(iset, id, ix, iq, ichk)

Get the value of a flavour momentum density
at a given point (`ix`, `iq`).

# Arguments
- `iset::Integer`: pdf set id (1-24)
- `id::Integer`: basis pdf identifier from -6 to 6+n,
where n is the number of additional pdfs in iset.
- `ix::Integer`: x index.
- `iq::Intger`: qq index.
- `ichk::Integer`: flag to steer error checking.
See allfxq().

# Returns
- `pdf::Float64`: pdf values.
"""
function fvalij(iset::Integer, id::Integer, ix::Integer, iq::Integer, ichk::Integer)

    iset = Ref{Int32}(iset)
    id = Ref{Int32}(id)
    ix = Ref{Int32}(ix)
    iq = Ref{Int32}(iq)
    ichk = Ref{Int32}(ichk)

    pdf = @ccall fvalij_(iset::Ref{Int32}, id::Ref{Int32}, ix::Ref{Int32},
                         iq::Ref{Int32}, ichk::Ref{Int32})::Float64

    pdf[]
end

"""
    fvalxq(iset, id, x, qq, ichk)

Get the value of a flavour momentum density
at a given point (`x`, `qq`).

# Arguments
- `iset::Integer`: pdf set id (1-24)
- `id::Integer`: basis pdf identifier from -6 to 6+n,
where n is the number of additional pdfs in iset.
- `x::Float64`: x value.
- `qq::Float64`: qq value.
- `ichk::Integer`: flag to steer error checking.
See allfxq().

# Returns
- `pdf::Float64`: pdf values.
"""
function fvalxq(iset::Integer, id::Integer, x::Float64, qq::Float64, ichk::Integer)

    iset = Ref{Int32}(iset)
    id = Ref{Int32}(id)
    x = Ref{Float64}(x)
    q = Ref{Float64}(qq)
    ichk = Ref{Int32}(ichk)

    pdf = @ccall fvalxq_(iset::Ref{Int32}, id::Ref{Int32}, x::Ref{Float64},
                         qq::Ref{Float64}, ichk::Ref{Int32})::Float64

    pdf[]
end


text:
# These are functions to convert SimpleGraph's to formats in Julia's
# Graphs module.

# Convert a graph or digraph into a "simple_graph".
#
# We return a triple: the simple_graph H with vertex set 1:n, a
# dictionary d mapping the vertices of G to the vertices of H, and the
# inverse dictionary dinv from V(H) to V(G).

import Graphs.simple_graph, Graphs.add_edge!
export convert_simple

"""
`convert_simple(G)` converts the `SimpleGraph` or `SimpleDigraph` to a
`simple_graph` type from Julia's `Graphs` module. This returns a
three-tuple consisting of:

* a `simple_graph` from which `G` was converted.
* a `Dict` mapping vertices of `G` to the vertices of the output graph.
* another `Dict` mapping vertics of the output graph to vertices of `G`.

Here is an example:
```
julia> G = StringGraph()
SimpleGraphs.SimpleGraph{String} (0 vertices)

julia> add!(G,"alpha", "beta")
true

julia> add!(G,"beta", "gamma")
true

julia> (g, d1, d2) = convert_simple(G);

julia> g
Undirected Graph (3 vertices, 2 edges)

julia> d1
Dict{String,Int64} with 3 entries:
  "alpha" => 1
  "gamma" => 3
  "beta"  => 2

julia> d2
Dict{Int64,String} with 3 entries:
  2 => "beta"
  3 => "gamma"
  1 => "alpha"
```
"""

function convert_simple(G::AbstractSimpleGraph)
    T = vertex_type(G)
    n = NV(G)
    has_dir = isa(G,SimpleDigraph)

    d = vertex2idx(G)
    dinv = Dict{Int,T}()
    for k in keys(d)
        v = d[k]
        dinv[v] = k
    end

    H = simple_graph(n,is_directed=has_dir)

    EE = elist(G)
    for e in EE
        u = d[e[1]]
        v = d[e[2]]
        add_edge!(H,u,v)
    end
    return (H,d,dinv)
end


text:
using JLD2, Distributions, Serialization, Dates, Statistics, DataFrames
import GEAS:sim_QTL, breeding_program, top_QTL, sim_gamma_QTL, sim_pt_QTL

function current_qtl(base, nqtl, op)
    if op == 1                  # gamma
        return sim_gamma_QTL(base, nqtl, nqtl)
    elseif op == 2              # mv normal
        d = MvNormal(zeros(2), [1 0; 0 1])
        return sim_pt_QTL(base, nqtl, d)
    elseif op == 3              # Laplacian
        return sim_QTL(base, nqtl, nqtl)
    end
end

"""
    function the_simulation(jld, nqtl, nrpt; rst = "")
---
## Description
- 3 comparisons
  - Normal genomic selection, no known QTL
  - Known 5 QTL as fixed effects
  - Edit the biggest of known 5 QTL in generation 2-6

## parameters
1. jld: the base
2. nqtl: e.g., 500
4. nrpt: number of repeats
5. rst: result path, in `2021-06-30/`
"""
function the_simulation(jld, nqtl, nrpt, desc, op; rst = "")
    @load jld base
    if length(rst) > 1
        isdir(rst) || mkpath(rst)
    end
    slog = joinpath(rst, "simu.log") # simulation log
    write(slog, desc)

    w = ndigits(nrpt)
    for r in 1:nrpt
        #qtl = sim_QTL(base, nqtl, nqtl)
        qtl = current_qtl(base, nqtl, op)

        suffix = lpad(string(r), w, '0')
        serialize(joinpath(rst, "qtl-$suffix.ser"), qtl)
        par = begin
            Parameters = Dict(
                :nSire=> 100,
                :nDam => 200,   # male:female = 1:2
                :nSib => 40,
                :nC7e => 30,    # number of sib to be challenged
                :nG8n => 11,    # number of generations
                :nQTL => [nqtl, nqtl],
                :h²   => [.5, .5],
                :p8e  => .5,    # death rate after challenge
                :e19e => 1.,    # edit_successsful_rate
                :w4t  => 1,     # weight on the binary trait EBV
                :nk3n => 10,    # number of known QTL on binary trait
                :dd   => 0.01, # A small value to be added to diagonals
                :log  => slog
            )
            (; Parameters...)       # named tuple.  contents as above
        end

        # genomic selection with no known QTL
        prd, snp = breeding_program(base, par, qtl)
        serialize(joinpath(rst, "prd-a-$suffix.ser"), prd)
        serialize(joinpath(rst, "snp-a-$suffix.ser"), snp)

        # genomic selection with known QTL as fixed effects
        prd, snp = breeding_program(base, par, qtl, fixed=true)
        serialize(joinpath(rst, "prd-b-$suffix.ser"), prd)
        serialize(joinpath(rst, "snp-b-$suffix.ser"), snp)

        # genomic selection, edit biggest known QTL each generation
        prd, snp = breeding_program(base, par, qtl, edit=true)
        serialize(joinpath(rst, "prd-c-$suffix.ser"), prd)
        serialize(joinpath(rst, "snp-c-$suffix.ser"), snp)
    end
end

function simulation_01()
    rst = joinpath("dat", string(DateTime(now())))
    @info "Results are in $rst"

    desc = "Uses a Γ(0.4) distritubion.
The Gamma distribution was modified to have half positive and half negative.
Total 100 QTL.
Number of known QTL 10.
Rank QTL fist of MAF(0.1, 0.9), then on abs(effects).

The top QTL of each round of simulation are as below:
"
    the_simulation("dat/run/base.jld2", 100, 10, desc, 1, rst=rst)
end

function simulation_02()
    rst = joinpath("dat", string(DateTime(now())))
    @info "Results are in $rst"

    desc = "Uses a pleiotropic MvNormal distritubion.
Total 100 QTL.
Number of known QTL 10.
Rank QTL fist of MAF(0.1, 0.9), then on abs(effects).

The top QTL of each round of simulation are as below:
"
    the_simulation("dat/run/base.jld2", 100, 10, desc, 2, rst=rst)
end

function simulation_03()
    rst = joinpath("dat", string(DateTime(now())))
    @info "Results are in $rst"

    desc = "Uses a independent Laplacian distribution.
Total 100 QTL.
Number of known QTL 10.
Rank QTL fist of MAF(0.1, 0.9), then on abs(effects).

The top QTL of each round of simulation are as below:
"
    the_simulation("dat/run/base.jld2", 100, 10, desc, 3, rst=rst)
end

function simulation_04()
    rst = joinpath("dat", string(DateTime(now())))
    @info "Results are in $rst"

    desc = "Uses a independent Laplacian distribution.
Total 500 QTL.
Number of known QTL 10.
Rank QTL fist of MAF(0.1, 0.9), then on abs(effects).

The top QTL of each round of simulation are as below:
"
    the_simulation("dat/run/base.jld2", 500, 10, desc, 3, rst=rst)
end

function qfrq(snp)
    fq = zeros(11, 10)
    start = 1
    i = 1
    for stop = 4000:4000:44000
        f = mean(snp[:, start:stop], dims=2)
        fq[i, :] = f'
        i += 1
        start = stop + 1
    end
    fq
end

"""
Generation mean and variance of `prd-abc` on `c`, e.g., `:tbv`.
"""
function mv(c, prd...)
    m = zeros(11, 3)
    v = zeros(11, 3)
    for i in 1:3
        gp = groupby(prd[i], :g8n)
        t = combine(gp, c => mean => :m, c => var => :v)
        m[:, i] = t.m
        v[:, i] = t.v
    end
    return m, v
end

function summarize_202107(dir, rpt)
    ir = lpad(string(rpt), 2, '0')
    @info "Read top QTL of the first repeat"
    log = readlines(joinpath(dir, "simu.log"))
    t = 1
    for i in 1:length(log)
        t += 1
        length(log[i]) > 0 && log[i][end] == ':' && break
    end
    t += 3(rpt - 1)
    tq = abs.(parse.(Int, split(log[t])))

    @info "Read SNP of the first repeat, ordinary GS"
    snp = deserialize(joinpath(dir, "snp-a-$ir.ser"))
    fq = qfrq(snp[tq, :])
    plot(fq, leg=false, title = "ogs")
    savefig(joinpath(dir, "afq.pdf"))

    @info "Read SNP of the first repeat, top QTL as fixed effects"
    snp = deserialize(joinpath(dir, "snp-b-$ir.ser"))
    fq = qfrq(snp[tq, :])
    plot(fq, leg=false, title = "fxq")
    savefig(joinpath(dir, "bfq.pdf"))

    @info "Read SNP of the first repeat, top QTL edited in turn"
    snp = deserialize(joinpath(dir, "snp-c-$ir.ser"))
    fq = qfrq(snp[tq, :])
    plot(fq, leg=false, title = "edt")
    savefig(joinpath(dir, "cfq.pdf"))

    @info "Read the production trait"
    lb = ["ogs" "fxq" "edt"]
    pa = deserialize(joinpath(dir, "prd-a-$ir.ser"))
    pb = deserialize(joinpath(dir, "prd-b-$ir.ser"))
    pc = deserialize(joinpath(dir, "prd-c-$ir.ser"))
    m, v = mv(:t2c, pa, pb, pc)
    plot(m, label = lb, xlabel="generation", ylabel="Bin mean")
    savefig(joinpath(dir, "bin-mean.pdf"))
    plot(v, label = lb, xlabel="generation", ylabel="Bin var")
    savefig(joinpath(dir, "bin-var.pdf"))
    m, v = mv(:tbv, pa, pb, pc)
    plot(m, label = lb, xlabel="generation", ylabel="Prd mean")
    savefig(joinpath(dir, "prd-mean.pdf"))
    plot(v, label = lb, xlabel="generation", ylabel="Prd var")
    savefig(joinpath(dir, "prd-var.pdf"))
end

function sum_mmv(dir, nr)
    mm = zeros(11, 3)
    vm = zeros(11, 3)
    for i in 1:nr
        r = lpad(string(i), 2, '0')
        pa = deserialize(joinpath(dir, "prd-a-$r.ser"))
        pb = deserialize(joinpath(dir, "prd-b-$r.ser"))
        pc = deserialize(joinpath(dir, "prd-c-$r.ser"))
        m, v = mv(:t2c, pa, pb, pc)
        mm += m
        vm += v
    end
    mm ./= nr
    vm ./= nr
    lb = ["ogs" "fxq" "edt"]
    plot(mm, label = lb, xlabel="generation", ylabel="Average bin mean")
    savefig(joinpath(dir, "ave-bin-mean.pdf"))
    plot(vm, label = lb, xlabel="generation", ylabel="Average bin var")
    savefig(joinpath(dir, "ave-bin-var.pdf"))
end

#=
!!!!! Obsolated !!!!!

"""
This is to use Laplacian distribution to have a few large QTL.
Only 100 QTL were simulated to check the editing effects.
Before, 500 QTL scenario show little advantage of the editing.
And, emphasize known QTL effects give negative selection effects.
"""
function simulation_2021_07_03()
    the_simulation("dat/run/base.jld2", 100, 10, "test", rst="dat/2021-07-03")
end

function simulation_2021_07_04()
    rst = joinpath("dat", string(DateTime(now())))
    @info "Results are in $rst"

    desc = "The previous simulation just recorded the top QTL positions in the last repeat.
This simulation changed `w` to `a` for QTL log.
This also uses a Gamma distritubion.
The Gamma distribution was modified to have half positive and half negative.
Number of known QTL changed from 5 to 10.

The top QTL of each round of simulation are as below:
"
    the_simulation("dat/run/base.jld2", 100, 10, desc, rst=rst)
end

function simulation_2021_07_05()
    rst = joinpath("dat", string(DateTime(now())))
    @info "Results are in $rst"

    desc = "The previous simulation just recorded the top QTL positions in the last repeat.
This simulation changed `desc` writing position out of the loop.
This also uses a Gamma distritubion.
The Gamma distribution was modified to have half positive and half negative.
Number of known QTL changed from 5 to 10.
QTL ranking was changed from var to effect only.
Only MAF of [0.1, 0.9] are to be considered.

The top QTL of each round of simulation are as below:
"
    the_simulation("dat/run/base.jld2", 100, 10, desc, rst=rst)
end
=#


text:
using CxxWrap

current_dir = @__DIR__
cxxhome = CxxWrap.prefix().path
juliahome = dirname(Base.Sys.BINDIR)

# Need to compile with the same compiler that nGraph is.
cxx = "clang++"

cxxflags = [
    "-g",
    "-O3",
    "-Wall",
    "-fPIC",
    "-std=c++17",
    "-DPCM_SILENT",
    "-DJULIA_ENABLE_THREADING",
    "-Dexcept_EXPORTS",
    # Surpress some warnings from Cxx
    "-Wno-unused-variable",
    "-Wno-unused-lambda-capture",
]

includes = [
    "-I$(joinpath(cxxhome, "include"))",
    "-I$(joinpath(juliahome, "include", "julia"))",
    "-I$(joinpath(@__DIR__, "usr", "include"))",
]

loadflags = [
    # Linking flags for Julia
    "-L$(joinpath(juliahome, "lib"))",
    "-Wl,--export-dynamic",
    "-Wl,-rpath,$(joinpath(juliahome, "lib"))",
    "-ljulia",
    # Linking Flags for CxxWrap
    "-L$(joinpath(cxxhome, "lib"))",
    "-Wl,-rpath,$(joinpath(cxxhome, "lib"))",
    "-lcxxwrap_julia",
]

src = joinpath(current_dir, "sync.cpp")
so = joinpath(current_dir, "libsync.so")

cmd = `$cxx $cxxflags $includes -shared $src -lpthread -o $so $loadflags`
@show cmd
run(cmd)


text:
import Random

const RNG_SEED = 1234567
const rng = Random.MersenneTwister(RNG_SEED)

function rand_wrapper(args...)
  rand(rng, args...)
end

function gen_rand_clause(; max_variables::Unsigned, max_clause_len::Unsigned)::Clause
  clause_len = rand_wrapper(1:max_clause_len)
  clauses = [rand_wrapper((-1, 1)) * signed(rand_wrapper(1:max_variables)) for _ in 1:clause_len]
  Clause(clauses)
end

function gen_rand_formula(;
                          max_clauses::Unsigned,
                          max_variables::Unsigned,
                          max_clause_len::Unsigned)::Formula
  clause_cnt = rand_wrapper(1:max_clauses)
  Formula([gen_rand_clause(;
                           max_variables=max_variables,
                           max_clause_len=max_clause_len) for _ in 1:clause_cnt])
end


text:
using Gomah

@testset "train" begin
    train()
end

@testset "predict" begin
    predict()
    # after training remove generated things
    rm("result", recursive=true)
end


text:
# Implementation of projection operation under missing data. Probably most of this needs a
# custom gradient implementation.
function project(
    S::Diagonal{T},
    U::AbstractMatrix{T},
    Y::ColVecs{Union{Missing, T}},
    σ²::T,
    D::Diagonal{T},
) where {T<:Real}

    # Compute patterns and assignments of data to patterns.
    patterns, rows, idxs, perm = compute_patterns(Y.X)

    # Construct the projection matrix for each pattern.
    Us = map(row -> U[row, :], rows)

    # Project each pattern in the data.
    sqrtS = cholesky(S).U
    Yproj_blocks = map((U, r, idx) -> sqrtS \ (U \ Y.X[r, idx]), Us, rows, idxs)

    # Construct the single projected data matrix, in the correct order.
    Yproj = hcat(Yproj_blocks...)[:, perm]

    # lens = map(length, idxs)
    lens = map_length(idxs)
    almost_ΣT_blocks = map(
        (U, len) -> repeat(diag(inv(cholesky(Symmetric(U'U + 1e-9I)))), 1, len), Us, lens,
    )

    # Assemble blocks into a single matrix, with a different scalar observation noise
    # for each observation.
    almost_ΣT = hcat(almost_ΣT_blocks...)[:, perm]
    ΣT = σ² .* diag(inv(S)) .* almost_ΣT .+ diag(D)

    return Yproj, ΣT
end

function regulariser(
    S::Diagonal{T},
    U::AbstractMatrix{T},
    σ²::T,
    y::ColVecs{Union{Missing, T}},
) where {T<:Real}

    # Compute patterns and assignments of data to patterns.
    patterns, rows, idxs, _ = compute_patterns(y.X)

    # Construct the projection matrix for each pattern.
    Us = map(row -> U[row, :], rows)

    # Pre-compute one term.
    logdet_S = logdet(cholesky(S))

    # Compute the regularisation term for each block.
    regularisers = map(
        (Uo, row, idx) -> begin
            n = length(idx)
            p, m = size(Uo)

            chol_UotUo = cholesky(Symmetric(Uo'Uo + 1e-9I))
            Yo = eltype(Uo).(y.X[row, idx])

            return -(n * (logdet_S + logdet(chol_UotUo) + (p - m) * log(2π * σ²)) +
                (sum(abs2, Yo) - sum(abs2, chol_UotUo.U' \ Uo'Yo)) / σ²) / 2
        end,
        Us, rows, idxs,
    )
    return sum(regularisers)
end

# Helper function for `project` that handles various bits of non-differentiable stuff
# that can be safely @nograd-ed.
function compute_patterns(Y::AbstractMatrix{Union{Missing, T}} where {T<:Real})

    # Compute unique missing-ness patterns.
    missingness = eachcol(ismissing.(Y))
    patterns = unique(collect.(missingness))

    # For each pattern, compute the rows of `Y` that are not missing.
    available_rows = [filter(n -> !p[n], 1:size(Y, 1)) for p in patterns]

    # Add the location each column of `Y` to mapping from block-to-columns.
    # idxs = Dict(pattern => Int[] for pattern in patterns)
    idxs = [Int[] for pattern in patterns]
    for (n, data_pattern) in enumerate(missingness)
        for (m, pattern) in enumerate(patterns)
            if data_pattern == pattern
                push!(idxs[m], n)
            end
        end
    end

    # Compute the permutation of the data required to restore the original order.
    perm = sortperm(vcat(idxs...))

    return patterns, available_rows, idxs, perm
end

ChainRulesCore.@non_differentiable compute_patterns(::Any)

map_length(x) = map(length, x)

ChainRulesCore.@non_differentiable map_length(::Any)


text:
# Utils

function summaporbroadcast(f, dists::AbstractArray, x::AbstractArray)
    # Broadcasting here breaks Tracker for some reason
    return sum(map(f, dists, x))
end
function summaporbroadcast(f, dists::AbstractVector, x::AbstractMatrix)
    return map(x -> summaporbroadcast(f, dists, x), eachcol(x))
end
@init @require LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" begin
    function summaporbroadcast(f, dists::LazyArrays.BroadcastArray, x::AbstractArray)
        return sum(copy(f.(dists, x)))
    end
    function summaporbroadcast(f, dists::LazyArrays.BroadcastVector, x::AbstractMatrix)
        return vec(sum(copy(f.(dists, x)), dims = 1))
    end
    lazyarray(f, x...) = LazyArrays.LazyArray(Base.broadcasted(f, x...))
    export lazyarray
end

# Univariate

const VectorOfUnivariate = Distributions.Product

function arraydist(dists::AbstractVector{<:UnivariateDistribution})
    return Product(dists)
end

function Distributions.logpdf(dist::VectorOfUnivariate, x::AbstractVector{<:Real})
    return summaporbroadcast(logpdf, dist.v, x)
end
function Distributions.logpdf(dist::VectorOfUnivariate, x::AbstractMatrix{<:Real})
    # eachcol breaks Zygote, so we need an adjoint
    return summaporbroadcast(logpdf, dist.v, x)
end
ZygoteRules.@adjoint function Distributions.logpdf(
    dist::VectorOfUnivariate,
    x::AbstractMatrix{<:Real}
)
    # Any other more efficient implementation breaks Zygote
    f(dist, x) = [sum(logpdf.(dist.v, view(x, :, i))) for i in 1:size(x, 2)]
    return ZygoteRules.pullback(f, dist, x)
end

struct MatrixOfUnivariate{
    S <: ValueSupport,
    Tdist <: UnivariateDistribution{S},
    Tdists <: AbstractMatrix{Tdist},
} <: MatrixDistribution{S}
    dists::Tdists
end
Base.size(dist::MatrixOfUnivariate) = size(dist.dists)
function arraydist(dists::AbstractMatrix{<:UnivariateDistribution})
    return MatrixOfUnivariate(dists)
end
function Distributions.logpdf(dist::MatrixOfUnivariate, x::AbstractMatrix{<:Real})
    return summaporbroadcast(logpdf, dist.dists, x)
end
function Distributions.logpdf(dist::MatrixOfUnivariate, x::AbstractArray{<:AbstractMatrix{<:Real}})
    return map(x -> logpdf(dist, x), x)
end
function Distributions.logpdf(dist::MatrixOfUnivariate, x::AbstractArray{<:Matrix{<:Real}})
    return map(x -> logpdf(dist, x), x)
end
function Distributions.rand(rng::Random.AbstractRNG, dist::MatrixOfUnivariate)
    return rand.(Ref(rng), dist.dists)
end

# Multivariate

struct VectorOfMultivariate{
    S <: ValueSupport,
    Tdist <: MultivariateDistribution{S},
    Tdists <: AbstractVector{Tdist},
} <: MatrixDistribution{S}
    dists::Tdists
end
Base.size(dist::VectorOfMultivariate) = (length(dist.dists[1]), length(dist))
Base.length(dist::VectorOfMultivariate) = length(dist.dists)
function arraydist(dists::AbstractVector{<:MultivariateDistribution})
    return VectorOfMultivariate(dists)
end
function Distributions.logpdf(dist::VectorOfMultivariate, x::AbstractMatrix{<:Real})
    # eachcol breaks Zygote, so we define an adjoint
    return sum(map(logpdf, dist.dists, eachcol(x)))
end
function Distributions.logpdf(dist::VectorOfMultivariate, x::AbstractArray{<:AbstractMatrix{<:Real}})
    return map(x -> logpdf(dist, x), x)
end
function Distributions.logpdf(dist::VectorOfMultivariate, x::AbstractArray{<:Matrix{<:Real}})
    return map(x -> logpdf(dist, x), x)
end
ZygoteRules.@adjoint function Distributions.logpdf(
    dist::VectorOfMultivariate,
    x::AbstractMatrix{<:Real}
)
    return ZygoteRules.pullback(dist, x) do dist, x
        sum(map(i -> logpdf(dist.dists[i], view(x, :, i)), 1:size(x, 2)))
    end
end
function Distributions.rand(rng::Random.AbstractRNG, dist::VectorOfMultivariate)
    init = reshape(rand(rng, dist.dists[1]), :, 1)
    return mapreduce(i -> rand(rng, dist.dists[i]), hcat, 2:length(dist); init = init)
end


text:
function is_cjk_character(ch::Union{Char, Int})::Bool
    #= check wheter a character is cjk charachter
    Args:
      ch: Union{Char, Int}, a character or an integer
    ref: https://en.wikipedia.org/wiki/CJK_Unified_Ideographs_(Unicode_block)
    =#
    ch = Int(ch)
    #return (0x4E00 <= ch <= 0x9FFF) || (0x3400 <= ch <= 0x4DBF) || (0x20000 <= ch <= 0x2A6DF) || (0x2A700 <= ch <= 0x2B73F) || (0x2B740 <= ch <= 0x2B81F) || (0x2B820 <= ch <= 0x2CEAF) || (0xF900 <= ch <= 0xFAFF) || (0x2F800 <= ch <= 0x2FA1F)
    return (19968 <= ch <= 40959) || (13312 <= ch <= 19903) || (131072 <= ch <= 173791) || (173824 <= ch <= 177983) || (177984 <= ch <= 178207) || (178208 <= ch <= 183983) || (63744 <= ch <= 64255) || (194560 <= ch <= 195103)
end

function truncate_text!(texts::Union{Vector{String}, Vector{Vector{String}}}, max_length::Int, offsets::Union{Int, Vector{Int}}, truncation::String)::Union{Vector{String}, Vector{Vector{String}}}
    #= truncate text
    Args:
      texts: Union{Vector{String}, Vector{Vector{String}}}, tokenized text or list of tokenized text
      max_length: Int, max tokens length
      offsets: Union{Int, Vector{Int}}, offsets of texts to delete
      truncation: String, `post` or `pre`
    =#

    @assert truncation in ["pre", "post"] "please specify direction from [`pre`, `post`]"

    if isa(texts, String)
        texts = [texts]
    end

    if isa(offsets, Int)
        offsets = [offsets for _ in 1:length(texts)]
    end

    while true
        lengths = [length(s) for s in texts]
        if sum(lengths) > max_length
            i = argmax(lengths)
            base_index = truncation == "post" ? lastindex(texts[i]) : firstindex(texts[i])
            deleteat!(texts[i], base_index + offsets[i])
        else
            return texts
        end
    end
end


text:
@testset "PolytopeComplements" begin

    @testset "Basic" begin

        HS = HalfSpace([-1.0, -1.0], 0.0)
        PC = PolytopeComplement(HS)

        @test [5.0, 5.0]   ∈ HS && [5.0, 5.0]   ∉ PC
        @test [-1.0, -1.0] ∉ HS && [-1.0, -1.0] ∈ PC

        @test complement(PC) == HS
        @test complement(HS) == PC

        # Hyperrectangle contained in HS
        hr = Hyperrectangle(low = [1.0, 1.0], high = [2.0, 2.0])

        @test is_intersection_empty(hr, PC) == true
        @test is_intersection_empty(hr, HS) == false
        @test (hr ⊆ PC) == false
        @test (hr ⊆ HS) == true
        # Hyperrectangle overlapping with HS
        hr = Hyperrectangle(low = [-1.0, -1.0], high = [1.0, 1.0])
        @test is_intersection_empty(hr, PC) == false
        @test is_intersection_empty(hr, HS) == false
        @test !(hr ⊆ PC)
        @test !(hr ⊆ HS)

        # Test some other sets
        @test @no_error complement(Hyperrectangle(ones(2), ones(2)))
        @test @no_error complement(Ball2(ones(2), 1.0))
        @test @no_error complement(Ball1(ones(3), 1.0))
        @test @no_error complement(Zonotope(ones(4), ones(4, 2)))
        @test @no_error complement(convert(HPolytope, hr))
    end

    @testset "Solvers with PCs" begin

        small_nnet = read_nnet("$(@__DIR__)/../examples/networks/small_nnet_id.nnet")
        in_hyper  = Hyperrectangle(low = [-0.9], high = [0.9])

        # Output sets that are the PolytopeComplements of the complements of the output sets used in the regular tests.
        problem_holds    = Problem(small_nnet, in_hyper, PolytopeComplement(HPolytope([HalfSpace([-1.0], 10.0)])))
        problem_violated = Problem(small_nnet, in_hyper, PolytopeComplement(HPolytope([HalfSpace([1.0], -20.0)])))

        for solver in [NSVerify(), MIPVerify(), ILP(), Reluplex(), Planet()]
            holds    = solve(solver, problem_holds)
            violated = solve(solver, problem_violated)

            @testset "$(typeof(solver))" begin
                @test holds.status    ∈ (:holds, :unknown)
                @test violated.status ∈ (:violated, :unknown)
            end
        end
    end
end


text:
using Test
using Solitaire

@testset "individual moves" begin
    Card = Solitaire.Card
    A = Solitaire.A
    B = Solitaire.B

    @testset "move_A" begin
        @test Solitaire.move_A(Card[A,1,2]) == Card[1,A,2]
        @test Solitaire.move_A(Card[1,A,2]) == Card[1,2,A]
        @test Solitaire.move_A(Card[1,2,A]) == Card[1,A,2]
    end

    @testset "move_B" begin
        @test Solitaire.move_B(Card[B,1,2]) == Card[1,2,B]
        @test Solitaire.move_B(Card[1,B,2]) == Card[1,B,2]
        @test Solitaire.move_B(Card[1,2,B]) == Card[1,2,B]
        @test Solitaire.move_B(Card[1,2,B,3]) == Card[1,B,2,3]
    end

    @testset "cut_AB" begin
        @test Solitaire.cut_AB(Card[A,B]) == Card[A,B]
        @test Solitaire.cut_AB(Card[B,A]) == Card[B,A]
        @test Solitaire.cut_AB(Card[1,A,B]) == Card[A,B,1]
        @test Solitaire.cut_AB(Card[A,1,B]) == Card[A,1,B]
        @test Solitaire.cut_AB(Card[A,B,1]) == Card[1,A,B]
        @test Solitaire.cut_AB(Card[1,2,A,3,4,B,5,6]) == Card[5,6,A,3,4,B,1,2]
        @test Solitaire.cut_AB(Card[1,2,B,3,4,A,5,6]) == Card[5,6,B,3,4,A,1,2]
    end

    @testset "cut_bottom" begin
        @test Solitaire.cut_bottom(Card[4,3,2,1]) == Card[3,2,4,1]
        @test Solitaire.cut_bottom(Card[1,4,3,2]) == Card[3,1,4,2]
        @test Solitaire.cut_bottom(Card[4,1,2,3]) == Card[4,1,2,3]
    end

    @testset "value" begin
        @test Solitaire.value(Card[1, 2, 3]) == 2
        @test Solitaire.value(Card[2, 1, 3]) == 3
        @test Solitaire.value(Card[1, 30, 2]) == 30  # no mod26 yet
        @test Solitaire.value(Card[1, A, 32]) == 53  # no skipping yet
    end

    @testset "cut_given" begin
        @test Solitaire.cut_given(Card[1, 2, 3], 1) == Card[2, 1, 3]
        @test Solitaire.cut_given(Card[1, 2, 3], 2) == Card[1, 2, 3]
        @test Solitaire.cut_given(Card[1, 2, 3, 4], 2) == Card[3, 1, 2, 4]
    end
end

@testset "raw generator output" begin
    # reference output
    ref = [4, 49, 10, 53, 24, 8, 51, 44, 6, 4, 33, 20, 39, 19, 34, 42]
    iter = Solitaire.DeckIterator(Solitaire.sorted_deck())
    for t in zip(ref, iter)
        @test t[1] == t[2]
    end
end

@testset "raw generator output" begin
    # reference output
    ref = [4, 49, 10, 24, 8, 51, 44, 6, 4, 33, 20, 39, 19, 34, 42]
    iter = Solitaire.skip_jokers(Solitaire.DeckIterator(Solitaire.sorted_deck()))
    for t in zip(ref, iter)
        @test t[1] == t[2]
    end
end

@testset "encryption with passphrase" begin
    @test Solitaire.encrypt("AAAAAAAAAAAAAAA", "") == "EXKYIZSGEHUNTIQ"
    @test Solitaire.encrypt("AAAAAAAAAAAAAAA", "f") == "XYIUQBMHKKJBEGY"
    @test Solitaire.encrypt("AAAAAAAAAAAAAAA", "fo") == "TUJYMBERLGXNDIW"
    @test Solitaire.encrypt("AAAAAAAAAAAAAAA", "foo") == "ITHZUJIWGRFARMW"
    @test Solitaire.encrypt("AAAAAAAAAAAAAAA", "a") == "XODALGSCULIQNSC"
    @test Solitaire.encrypt("AAAAAAAAAAAAAAA", "aa") == "OHGWMXXCAIMCIQP"
    @test Solitaire.encrypt("AAAAAAAAAAAAAAA", "aaa") == "DCSQYHBQZNGDRUT"
    @test Solitaire.encrypt("AAAAAAAAAAAAAAA", "b") == "XQEEMOITLZVDSQS"
    @test Solitaire.encrypt("AAAAAAAAAAAAAAA", "bc") == "QNGRKQIHCLGWSCE"
    @test Solitaire.encrypt("AAAAAAAAAAAAAAA", "bcd") == "FMUBYBMAXHNQXCJ"
    @test Solitaire.encrypt("AAAAAAAAAAAAAAAAAAAAAAAAA", "cryptonomicon") == "SUGSRSXSWQRMXOHIPBFPXARYQ"
    @test Solitaire.encrypt("SOLITAIRE", "cryptonomicon") == "KIRAKSFJA"
end

@testset "decryption with passphrase" begin
    @test Solitaire.decrypt("EXKYIZSGEHUNTIQ", "") == "AAAAAAAAAAAAAAA"
    @test Solitaire.decrypt("XYIUQBMHKKJBEGY", "f") == "AAAAAAAAAAAAAAA"
    @test Solitaire.decrypt("TUJYMBERLGXNDIW", "fo") == "AAAAAAAAAAAAAAA"
    @test Solitaire.decrypt("ITHZUJIWGRFARMW", "foo") == "AAAAAAAAAAAAAAA"
    @test Solitaire.decrypt("XODALGSCULIQNSC", "a") == "AAAAAAAAAAAAAAA"
    @test Solitaire.decrypt("OHGWMXXCAIMCIQP", "aa") == "AAAAAAAAAAAAAAA"
    @test Solitaire.decrypt("DCSQYHBQZNGDRUT", "aaa") == "AAAAAAAAAAAAAAA"
    @test Solitaire.decrypt("XQEEMOITLZVDSQS", "b") == "AAAAAAAAAAAAAAA"
    @test Solitaire.decrypt("QNGRKQIHCLGWSCE", "bc") == "AAAAAAAAAAAAAAA"
    @test Solitaire.decrypt("FMUBYBMAXHNQXCJ", "bcd") == "AAAAAAAAAAAAAAA"
    @test Solitaire.decrypt("SUGSRSXSWQRMXOHIPBFPXARYQ", "cryptonomicon") == "AAAAAAAAAAAAAAAAAAAAAAAAA"
    @test Solitaire.decrypt("KIRAKSFJA", "cryptonomicon") == "SOLITAIRE"
end


text:
#  TODO: Change constants to tuples to prevent accidentally mutating the values?

const GLOBAL_ATTRIBUTES_KEYS = [
    "edition", "centre",
    "centreDescription", "subCentre"
]

const DATA_ATTRIBUTES_KEYS = [
    "paramId",
    "shortName",
    "units",
    "name",
    "cfName",
    "cfVarName",
    "dataType",
    "missingValue",
    "numberOfPoints",
    "totalNumber",
    "numberOfDirections",
    "numberOfFrequencies",
    "typeOfLevel",
    "NV",
    "stepUnits",
    "stepType",
    "gridType",
    "gridDefinitionDescription",
]

const GRID_TYPE_MAP = Dict(
    "regular_ll" => [
        "Nx",
        "iDirectionIncrementInDegrees",
        "iScansNegatively",
        "longitudeOfFirstGridPointInDegrees",
        "longitudeOfLastGridPointInDegrees",
        "Ny",
        "jDirectionIncrementInDegrees",
        "jPointsAreConsecutive",
        "jScansPositively",
        "latitudeOfFirstGridPointInDegrees",
        "latitudeOfLastGridPointInDegrees",
    ],
    "rotated_ll" => [
        "Nx",
        "Ny",
        "angleOfRotationInDegrees",
        "iDirectionIncrementInDegrees",
        "iScansNegatively",
        "jDirectionIncrementInDegrees",
        "jPointsAreConsecutive",
        "jScansPositively",
        "latitudeOfFirstGridPointInDegrees",
        "latitudeOfLastGridPointInDegrees",
        "latitudeOfSouthernPoleInDegrees",
        "longitudeOfFirstGridPointInDegrees",
        "longitudeOfLastGridPointInDegrees",
        "longitudeOfSouthernPoleInDegrees",
    ],
    "reduced_ll" => [
        "Ny",
        "jDirectionIncrementInDegrees",
        "jPointsAreConsecutive",
        "jScansPositively",
        "latitudeOfFirstGridPointInDegrees",
        "latitudeOfLastGridPointInDegrees",
    ],
    "regular_gg" => [
        "Nx",
        "iDirectionIncrementInDegrees",
        "iScansNegatively",
        "longitudeOfFirstGridPointInDegrees",
        "longitudeOfLastGridPointInDegrees",
        "N",
        "Ny",
    ],
    "rotated_gg" => [
        "Nx",
        "Ny",
        "angleOfRotationInDegrees",
        "iDirectionIncrementInDegrees",
        "iScansNegatively",
        "jPointsAreConsecutive",
        "jScansPositively",
        "latitudeOfFirstGridPointInDegrees",
        "latitudeOfLastGridPointInDegrees",
        "latitudeOfSouthernPoleInDegrees",
        "longitudeOfFirstGridPointInDegrees",
        "longitudeOfLastGridPointInDegrees",
        "longitudeOfSouthernPoleInDegrees",
        "N",
    ],
    "lambert" => [
        "LaDInDegrees",
        "LoVInDegrees",
        "iScansNegatively",
        "jPointsAreConsecutive",
        "jScansPositively",
        "latitudeOfFirstGridPointInDegrees",
        "latitudeOfSouthernPoleInDegrees",
        "longitudeOfFirstGridPointInDegrees",
        "longitudeOfSouthernPoleInDegrees",
        "DyInMetres",
        "DxInMetres",
        "Latin2InDegrees",
        "Latin1InDegrees",
        "Ny",
        "Nx",
    ],
    "reduced_gg" => ["N", "pl"],
    "sh" => ["M", "K", "J"],
)
const GRID_TYPE_KEYS = unique(vcat(values(GRID_TYPE_MAP)...))

const ENSEMBLE_KEYS = ["number"]
const VERTICAL_KEYS = ["level"]
const DATA_TIME_KEYS = ["dataDate", "dataTime", "endStep"]
const ALL_REF_TIME_KEYS = [
    "time", "step",
    "valid_time", "verifying_time",
    "forecastMonth"
]
const SPECTRA_KEYS = ["directionNumber", "frequencyNumber"]

const ALL_HEADER_DIMS = vcat(
    ENSEMBLE_KEYS,
    VERTICAL_KEYS,
    DATA_TIME_KEYS,
    ALL_REF_TIME_KEYS,
    SPECTRA_KEYS
)

#  TODO: Include the list of included keys in docs automatically
const ALL_KEYS = sort(unique(vcat(
    GLOBAL_ATTRIBUTES_KEYS, DATA_ATTRIBUTES_KEYS,
    GRID_TYPE_KEYS, ALL_HEADER_DIMS
)))

#  TODO: Include the list of coordinate attributes in docs automatically
const COORD_ATTRS = Dict(
    # geography
    "latitude"             => Dict(
        "units"            => "degrees_north",
        "standard_name"    => "latitude",
        "long_name"        => "latitude"
    ),
    "longitude"            => Dict(
        "units"            => "degrees_east",
        "standard_name"    => "longitude",
        "long_name"        => "longitude"
    ),
    # vertical
    "depthBelowLand"       => Dict(
        "units"            => "m",
        "positive"         => "down",
        "long_name"        => "soil depth",
        "standard_name"    => "depth",
    ),
    "depthBelowLandLayer"  => Dict(
        "units"            => "m",
        "positive"         => "down",
        "long_name"        => "soil depth",
        "standard_name"    => "depth",
    ),
    "hybrid"               => Dict(
        "units"            => "1",
        "positive"         => "down",
        "long_name"        => "hybrid level",
        "standard_name"    => "atmosphere_hybrid_sigma_pressure_coordinate",
    ),
    "heightAboveGround"    => Dict(
        "units"            => "m",
        "positive"         => "up",
        "long_name"        => "height above the surface",
        "standard_name"    => "height",
    ),
    "isobaricInhPa"        => Dict(
        "units"            => "hPa",
        "positive"         => "down",
        "stored_direction" => "decreasing",
        "standard_name"    => "air_pressure",
        "long_name"        => "pressure",
    ),
    "isobaricInPa"         => Dict(
        "units"            => "Pa",
        "positive"         => "down",
        "stored_direction" => "decreasing",
        "standard_name"    => "air_pressure",
        "long_name"        => "pressure",
    ),
    "isobaricLayer"        => Dict(
        "units"            => "Pa",
        "positive"         => "down",
        "standard_name"    => "air_pressure",
        "long_name"        => "pressure",
    ),
    # ensemble
    "number"               => Dict(
        "units"            => "1",
        "standard_name"    => "realization",
        "long_name"        => "ensemble member numerical id",
    ),
    # time
    "step"                 => Dict(
        "units"            => "hours",
        "standard_name"    => "forecast_period",
        "long_name"        => "time since forecast_reference_time",
    ),
    "time"                 => Dict(
        "units"            => "seconds since 1970-01-01T00:00:00",
        "calendar"         => "proleptic_gregorian",
        "standard_name"    => "forecast_reference_time",
        "long_name"        => "initial time of forecast",
    ),
    "valid_time"           => Dict(
        "units"            => "seconds since 1970-01-01T00:00:00",
        "calendar"         => "proleptic_gregorian",
        "standard_name"    => "time",
        "long_name"        => "time",
    ),
    "verifying_time"       => Dict(
        "units"            => "seconds since 1970-01-01T00:00:00",
        "calendar"         => "proleptic_gregorian",
        "standard_name"    => "time",
        "long_name"        => "time",
    ),
)

const GRID_TYPES_DIMENSION_COORDS = ["regular_ll", "regular_gg"]
const GRID_TYPES_2D_NON_DIMENSION_COORDS = [
    "rotated_ll",
    "rotated_gg",
    "lambert",
    "albers",
    "polar_stereographic",
]


text:
using Distributions
using SpecialFunctions
using JLD2
using FileIO
using Distributed
using LinearAlgebra

function start_pattern(l,up) # initiate 1 pattern
    sign.(rand(l).-up)
end

function get_patterns(l,n,up) # initiate n patterns
    xis=zeros(l,n)
    for ii in 1:n
        xis[:,ii]=start_pattern(l,up)
    end
    return xis
end

function updateactivestate(states)  #Get the state that is presented to the network
    return rand(1:states)
end

function build_matrix_xi_x_runin(l,n,x,γ,runin)   # initiate J_{ij} and n patterns  % γ = (1-λ)
        xi= get_patterns(l,n,0.5)
            xi,J = let_run_in_Hop_uniform(xi,γ,x,runin,n)
    return J,xi
end

function let_run_in_Hop_uniform(xi,γ,x,runs,states)  # initiate J_{ij}
    activestate=updateactivestate(states)

    J = start_J(xi[:,activestate],γ)

    for ii in 1:runs
        xi[:]=new_XI_rand(xi[:], x)
        J=new_J(J,xi[:,activestate],γ)
        activestate=updateactivestate(states)

    end
    return xi,J
end

function start_J(xi,γ)   # Start initial J_ij
    l =length(xi)
    J=zeros(l,l)
    for jj in 1:l
        for ii in 1:jj-1
            J[ii,jj] = (1-γ)* xi[ii]*xi[jj]
        end
        for ii in jj+1:l
            J[ii,jj] = (1-γ)* xi[ii]*xi[jj]
        end

    end
    return J
end

function new_XI_rand(old_XI, x) # random mutations of patterns
    l=length(old_XI)
    for ii in 1:l
        if (rand() < x)
            old_XI[ii] *=-1
        end
    end
    return old_XI
end

function new_J(old_J::Array{Float64,2},xi::Array{Float64,1},γ::Float64)  # update interaction with stadard Hebbian learning # γ = (1-λ)
    l=length(xi)

    for jj in 1:l
        for ii in 1:jj-1
            old_J[ii,jj] = γ* old_J[ii,jj] + (1-γ)* xi[ii]*xi[jj]
        end
        for ii in jj+1:l
            old_J[ii,jj] = γ* old_J[ii,jj] + (1-γ)* xi[ii]*xi[jj]
        end
    end
    return old_J
end

function energie(J,xi)  #energy of a pettern xi in network J
    l=length(xi)

    return (transpose(xi)*J*xi)/l/2
end

function localfiels(xi,J)
    return sum(sign.(xi.*(J*xi)))
end

function Support_entropy(xi,evecs,energies)   ## Support of the patterns in the eigenspace of the network given in units of entropy
    n=length(energies)
    l=length(xi[:,1])
    entropy_vec=Vector{Float64}()
    entropy_vec_normed=Vector{Float64}()
    overlaps=((transpose(xi)*evecs)./(sqrt(l))).^2

    for p in 1:n
        local_entropy=Vector{Float64}()
        local_entropy_normed=Vector{Float64}()
        norm=sum(overlaps[:,p])
        for ii in 1:n
            if (overlaps[ii,p]>0.000001)
                append!(local_entropy,-overlaps[ii,p]* log2(overlaps[ii,p]))
                append!(local_entropy_normed,-overlaps[ii,p]/norm* log2(overlaps[ii,p]/norm))
            end
        end
        append!(entropy_vec,sum(local_entropy))
        append!(entropy_vec_normed,sum(local_entropy_normed))
    end

    entropy_patts=Vector{Float64}()

    overlepsp=transpose(overlaps)

    for p in 1:n
        local_entropy=Vector{Float64}()
        for ii in 1:n
            if (overlepsp[ii,p]>0.000001)
                append!(local_entropy,-overlepsp[ii,p]* log2(overlepsp[ii,p]))
            end
        end
        append!(entropy_patts,sum(local_entropy))
    end

    return entropy_vec,entropy_vec[1],entropy_vec[n],entropy_vec_normed,entropy_vec_normed[1],entropy_vec_normed[n],entropy_patts,entropy_patts[argmax(energies)],entropy_patts[argmin(energies)]
end

function Support_partisipation(xi,evecs,energies)  ## Support of the patterns in the eigenspace of the network given in units of partisipation
    n=length(energies)
    l=length(xi[:,1])
    entropy_vec=Vector{Float64}()

    overlaps=((transpose(xi)*evecs)./(sqrt(l))).^2

    for p in 1:n
        norm=sum(overlaps[:,p].^4)
        partis=sum(overlaps[:,p].^2)^2
        append!(entropy_vec,partis/norm)
        end

    entropy_patts=Vector{Float64}()

    overlepsp=((transpose(evecs)*xi)./(sqrt(l))).^2

    for p in 1:n
        norm=sum(overlepsp[:,p].^4)
        partis=sum(overlepsp[:,p].^2)^2
        append!(entropy_patts,partis/norm)
    end

    return mean(entropy_vec),std(entropy_vec),entropy_vec[1],entropy_vec[n],mean(entropy_patts),std(entropy_patts),entropy_patts[argmax(energies)],entropy_patts[argmin(energies)]
end


text:
# Functions derived from "private" functions in Base

# print elapsed time, return expression value
const _mem_units = ["bytes", "KB", "MB", "GB", "TB", "PB"]
const _cnt_units = ["", " k", " M", " G", " T", " P"]
function prettyprint_getunits(value, numunits, factor)
    c1 = factor
    c2 = c1 * c1
    if value <= c1 * 100
        return (value, 1)
    end
    unit = 2
    while value > c2 * 100 && (unit < numunits)
        c1 = c2
        c2 *= factor
        unit += 1
    end
    return div(value + (c1 >>> 1), c1), unit
end

const _sec_units = ["nanoseconds ", "microseconds", "milliseconds", "seconds     "]
function prettyprint_nanoseconds(value::UInt64)
    if value < 1000
        return (1, value, 0)    # nanoseconds
    elseif value < 1000000
        mt = 2
    elseif value < 1000000000
        mt = 3
        # round to nearest # of microseconds
        value = div(value+500,1000)
    elseif value < 1000000000000
        mt = 4
        # round to nearest # of milliseconds
        value = div(value+500000,1000000)
    else
        # round to nearest # of seconds
        return (4, div(value+500000000,1000000000), 0)
    end
    frac::UInt64 = div(value,1000)
    return (mt, frac, value-(frac*1000))
end

function time_print(elapsedtime, bytes, gctime, allocs)
    mt, pptime, fraction = prettyprint_nanoseconds(elapsedtime)
    if fraction != 0
        @printf("%4d.%03d %s", pptime, fraction, _sec_units[mt])
    else
        @printf("%8d %s", pptime, _sec_units[mt])
    end
    if bytes != 0 || allocs != 0
        bytes, mb = prettyprint_getunits(bytes, length(_mem_units), 1024)
        allocs, ma = prettyprint_getunits(allocs, length(_cnt_units), 1000)
        @printf(
            " (%d%s allocation%s: %d %s",
            allocs,
            _cnt_units[ma],
            allocs == 1 ? "" : "s",
            bytes,
            _mem_units[mb]
        )
        if gctime > 0
            @printf(", %.2f%% gc time", 100 * gctime / elapsedtime)
        end
        print(")")
    elseif gctime > 0
        @printf(", %.2f%% gc time", 100 * gctime / elapsedtime)
    end
    println()
end


text:
struct getfield_types <:Rule
  enter::Function
  leave::Function
    simbolos::Dict
  function getfield_types()
    simbolos= Dict()
    nombre=""
    function enter(node)
      if (node.kind == "ObjectTypeDefinition")
         nombre = node.name.value
                 push!(simbolos, nombre => Dict())
      end
      if (node.kind == "FieldDefinition")
        nombrecampo= node.name.value
        push!(simbolos[nombre], nombrecampo => Dict() )
        push!(simbolos[nombre][nombrecampo], "tipo" =>node.tipe.name.value )

        if length(node.arguments)>0
          push!(simbolos[nombre][nombrecampo], "args" => Dict() )

                 for data in node.arguments
                     push!(simbolos[nombre][nombrecampo]["args"], data.name.value => data.tipe.name.value )
                 end
        end
      end
      if (node.kind == "OperationTypeDefinition")
        push!(simbolos, node.operation => node.tipe.name.value )
      end
    end
    function leave(node)
    end
    new(enter,leave,simbolos)
  end
end

struct extract_operation <:Rule
  enter::Function
  leave::Function
  found:: Dict
  function extract_operation(operationName::String)
    active= true
    found = Dict()
    function enter(node::Node)
      if (active==true) && (node.kind == "OperationDefinition")
        if node.name != nothing
           if (node.name.value == operationName)
              push!(found, "operation_type" => node.operation, "extracted_operation" => node)
              active=false
           end
        end
      end
    end
    function leave(node::Node)
    end
    new(enter,leave,found)
  end
end


text:
using DiffEqBase
using OrdinaryDiffEq
using SteadyStateDiffEq
using StochasticDiffEq
using DiffEqJump
using SparseArrays

import OrdinaryDiffEq: ODEProblem
import SteadyStateDiffEq: SteadyStateProblem
import StochasticDiffEq: SDEProblem
import DiffEqJump: JumpProblem

funcindex!(list, key, f, vals...) = list[key] = f(list[key],vals...)
valueat(x::Number, u, t) = x
valueat(f::Function, u, t) = try f(u,t) catch e f(t) end
transitionrate(S, T, k, rate, t) = exp(reduce((x,y)->x+log(S[y] <= 0 ? 0 : S[y]),
                                       keys(first(T[k]));
                                       init=log(valueat(rate[k],S,t))))

""" vectorfield(m::Model)

Convert a petri model into a differential equation function that can
be passed into DifferentialEquation.jl or OrdinaryDiffEq.jl solvers
"""
function vectorfield(m::Model)
    S = m.S
    T = m.Δ
    ϕ = Dict()
    f(du, u, p, t) = begin
        for k in keys(T)
          ϕ[k] = transitionrate(u, T, k, p, t)
        end
        for s in S
          du[s] = 0
        end
        for k in keys(T)
            l,r = T[k]
            for s in keys(l)
              funcindex!(du, s, -, ϕ[k] * l[s])
            end
            for s in keys(r)
              funcindex!(du, s, +, ϕ[k] * r[s])
            end
        end
        return du
    end
    return f
end

""" ODEProblem(m::Model, u0, tspan, β)

Generate an OrdinaryDiffEq ODEProblem
"""
ODEProblem(m::Model, u0, tspan, β) = ODEProblem(vectorfield(m), u0, tspan, β)

""" SteadyStateProblem(m::Model, u0, tspan, β)

Generate an SteadyStateDiffEq SteadyStateProblem
"""
SteadyStateProblem(m::Model, u0, tspan, β) = SteadyStateProblem(ODEProblem(m, u0, tspan, β))

function statecb(s)
     cond = (u,t,integrator) -> u[s]
     aff = (integrator) -> integrator.u[s] = 0.0
     return ContinuousCallback(cond, aff)
end

""" SDEProblem(m::Model, u0, tspan, β)

Generate an StochasticDiffEq SDEProblem and an appropriate CallbackSet
"""
function SDEProblem(m::Model, u0, tspan, β)
    S = m.S
    T = m.Δ
    ϕ = Dict()
    Spos = Dict(S[k]=>k for k in keys(S))
    T_keys = collect(keys(T))
    Tpos = Dict(T_keys[k]=>k for k in keys(T_keys))
    nu = spzeros(Float64, length(S), length(T))
    for k in T_keys
      l,r = T[k]
      for i in keys(l)
        nu[Spos[i],Tpos[k]] -= l[i]
      end
      for i in keys(r)
        nu[Spos[i],Tpos[k]] += r[i]
      end
    end
    noise(du, u, p, t) = begin
        sum_u = sum(u)
        for k in keys(T)
          ins = first(T[k])
          ϕ[k] = transitionrate(u, T, k, p, t)
        end

        for k in keys(T)
          l,r = T[k]
          rate = sqrt(abs(ϕ[k]))
          for i in keys(l)
            du[Spos[i],Tpos[k]] = -rate
          end
          for i in keys(r)
            du[Spos[i],Tpos[k]] = rate
          end
        end
        return du
    end
    return SDEProblem(vectorfield(m),noise,u0,tspan,β,noise_rate_prototype=nu),
           CallbackSet([statecb(s) for s in S]...)
end

jumpTransitionRate(T, k) = (u,p,t) -> transitionrate(u, T, k, p, t)

function jumpTransitionFunction(t)
  return (integrator) -> begin
    l,r = t
    for i in keys(l)
      integrator.u[i] -= l[i]
    end
    for i in keys(r)
      integrator.u[i] += r[i]
    end
  end
end

""" JumpProblem(m::Model, u0, tspan, β)

Generate an DiffEqJump JumpProblem
"""
function JumpProblem(m::Model, u0, tspan, β)
  T = m.Δ
  prob = DiscreteProblem(u0, tspan, β)
  return JumpProblem(prob, Direct(), [ConstantRateJump(jumpTransitionRate(T, k),
                                                       jumpTransitionFunction(T[k])
                                                      ) for k in keys(T)]...)
end


text:
@discord_object struct GuildWidget
    enabled::Bool
    channel_id::Snowflake
end


text:
macro labeled(body)
    # @alt `Base.@show' does essentially the same thing
    ##
    bodystr = string(body)
    quote
        out = $(esc(body))
        res = @> map(split(string(out),"\n")) do line
            "\t$line"
        end join("\n")
        println("$($bodystr) =>$(res)")
        out
    end
end
export @labeled


text:
@testset "Test generated structs" begin
    descriptor_file = joinpath(@__DIR__, "..", "src", "descriptors", "structs.json")
    existing_dir = joinpath(@__DIR__, "..", "src", "generated")
    @test IS.test_generated_structs(descriptor_file, existing_dir)
end