text:
function miexperimento(prob::Number; genj::GeneticJulia = GenJ)
    # Generates a population

    genj._population = genPopulation(genj._generator, genj._experimentInfo)
    evaluate!(genj._evaluator, genj._experimentInfo, genj._population, genj._stopCondition)

    while !reached(genj._stopCondition)
        selectedParents = selectParents(genj._selector, genj._population, genj._crossoverOp._nParents, genj._evaluator, genj._experimentInfo._rng)

        offspring = genj._mutationOp._probability > 0.0 ? mutate(genj._mutationOp, selectedParents, genj._experimentInfo) : selectedParents

        evaluate!(genj._evaluator, genj._experimentInfo, offspring, genj._stopCondition)

        random = rand(genj._experimentInfo._rng)
        if random < prob
            genj._population = replacePopulation(genj._replacementOp, genj._population, offspring, genj._evaluator, genj._experimentInfo._rng)
        end

        notifyIteration(genj._stopCondition)
    end
end

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
insts = OrderedDict()

insts["test"] = [
    ((true, false, 3, 4, true),),
    ((true, false, 3, 4, false),),
    ((false, false, 3, 4, true),),
    ((false, false, 3, 4, false),),
    ((true, true, 6, 6, true),),
    ((true, true, 6, 6, false),),
    ((false, true, 6, 6, true),),
    ((false, true, 6, 6, false),),
]

function matrixcompletion_insts(use_nat::Bool, nrow_max::Int)
    return [((true, true, nrow, nrow, use_nat),) for nrow in vcat(20, 20:20:nrow_max)]
end

insts["nat"] = matrixcompletion_insts(true, 220)
insts["ext"] = matrixcompletion_insts(false, 120)

return (MatrixCompletion, insts)

ext: jl


text:
include("../../../170220_RSsite/bin/GetNormalizedCoverage.jl")

pathBamFiles = ARGS[1]
pathRegions = ARGS[2]
odirPrefix = ARGS[3]

mkpath(odirPrefix)

getNormalizedCoverage(pathBamFiles, pathRegions, odirPrefix)


text:
# DIC
deviance(loglikelihood::Real) = -2 * loglikelihood
function dic(loglikelihood::AbstractVector{<:Real})
  D = deviance.(loglikelihood)
  return mean(D) + 0.5 * var(D)
end

"""
log(1 / mean(exp(x)))
"""
logharmmeanexp(x::AbstractVector{<:Real}) = log(length(x)) - logsumexp(-x)

"""
  log_bayes_factor(lls::Vector{AbstractVector{<:Real}})

Compute log bayes factor for all models with respect to first model.
"""
function log_bayes_factor(lls::Vector{<:AbstractVector{<:Real}})
  harmean_of_loglike = map(logharmmeanexp, lls)

  # Ratio of harmonic means.
  return harmean_of_loglike .- harmean_of_loglike[1]
end

"""
    log_bayes_factor(ll1::AbstractVector{<:Real}, ll2::AbstractVector{<:Real})

Returns: log bayes factor in favor of model2 (vs. model1).

# Arguments

`ll1`: loglikelihood evaluations for model 1 for each posterior sample.
`ll2`: loglikelihood evaluations for model 2 for each posterior sample.
"""
function log_bayes_factor(ll1::AbstractVector{<:Real}, ll2::AbstractVector{<:Real})
  return log_bayes_factor([ll1, ll2])[2]
end

"""
Numerically stable computation of `log(mean(exp.(xs)))`.
"""
logmeanexp(xs) = logsumexp(xs) - log(length(xs))

"""
Bhattacharyya distance. This is numerically the same as `-log(BC(F, G, n))`,
but more computationally efficent.
"""
function negativeLogBC(F::Distribution, G::Distribution, n::Integer)
  samps = rand(F, n)
  return -logmeanexp((logpdf.(G, samps) - logpdf.(F, samps)) / 2)
end

"""
Monte Carlo (MC) approximation of Bhattacharyya coefficient between two distributions
`F` and `G`. `n` is the number of MC samples to use for the approximation.

BC(F, G) = ∫ √(f(x) ⋅ √g(x) dx
         = ∫ √g(x) / √f(x) ⋅ f(x) dx

Thus, we sample from `F`, then approximate the integral by evaluating
√g(x) / √f(x).

"""
function BC(F::Distribution, G::Distribution, n::Integer)
  return exp(-negativeLogBC(F, G, n))
end

"""
Monte Carlo (MC) approximation of Hellinger distance between two distributions
`F` and `G`. `n` is the number of MC samples to use for the approximation.

H²(F, G) = 1 - ∫ √(f(x) ⋅ √g(x) dx
         = 1 - ∫ √g(x) / √f(x) ⋅ f(x) dx

Thus, we sample from `F`, then approximate the integral by evaluating
√g(x) / √f(x).
"""
hellinger(F::Distribution, G::Distribution, n::Integer) = sqrt(hellinger2(F, G, n))

"""
Hellinger squared distance. Note that `hellinger(F, G, n) == sqrt(hellinger2(F, G, n))`.
"""
function hellinger2(F::Distribution, G::Distribution, n::Integer)
  samps = rand(F, n)
  log_mc_approx_integral = logmeanexp((logpdf.(G, samps) - logpdf.(F, samps)) / 2)
  return 1 - BC(F, G, n)
end


text:
using Random, Test, WindingNelderMead
using WindingNelderMead: Vertex, Simplex, windingnumber, windingangle
using WindingNelderMead: centre, assessconvergence, position, value
using WindingNelderMead: bestvertex, issortedbyangle, hypervolume

@testset "WindingNelderMead tests" begin

  Random.seed!(0)

  @testset "Vertices" begin
    v1 = Vertex([0.0, 1.0], 2.0 - 3.0im)
    v2 = Vertex([0.0, 1.0], 2.0 - 3.0im)
    @test !(v1 == v2)
    @test isequal(v1, v2)
    @test hash(v1) == hash(v2)
  end

  @testset "Simplex tests" begin

    irrelevant = [0.0, 0.0]
    @testset "Simplex encloses zero" begin
      v1 = Vertex(irrelevant, 1.0 - im)
      v2 = Vertex(irrelevant, 0.0 + im)
      v3 = Vertex(irrelevant, -1.0 - im)
      encloseszero = Simplex([v1, v2, v3])
      @test isapprox(1, abs(windingangle(encloseszero)) / (2π))
      @test windingnumber(encloseszero) == 1
    end

    @testset "Simplex doesn't enclose zero" begin
      v1 = Vertex(irrelevant, 1.0 + im)
      v2 = Vertex(irrelevant, 2.0 + im)
      v3 = Vertex(irrelevant, 1.0 + im * 2)
      doesntenclosezero = Simplex([v1, v2, v3])
      @test abs(windingangle(doesntenclosezero)) / (2π) < 1.0e-3
      @test windingnumber(doesntenclosezero) == 0
      v1 = Vertex(irrelevant, 1.0 + im)
      v2 = Vertex(irrelevant, 0.0 + 2*im)
      v3 = Vertex(irrelevant, -1.0 + im)
      doesntenclosezero = Simplex([v1, v2, v3])
      @test abs(windingangle(doesntenclosezero)) / (2π) < 1.0e-3
      @test windingnumber(doesntenclosezero) == 0
      v1 = Vertex(irrelevant, 1.0 + 0im)
      v2 = Vertex(irrelevant, 0.0 + im)
      v3 = Vertex(irrelevant, 1.0 + 0im)
      doesntenclosezero = Simplex([v1, v2, v3])
      @test abs(windingangle(doesntenclosezero)) / (2π) < 1.0e-3
      @test windingnumber(doesntenclosezero) == 0
    end

    @testset "Simplices with identical vertices is converged" begin
      dim = 2
      T = Float64
      U = ComplexF64
      pos = T[1.0, 1.0]
      val = rand(U)
      v1 = Vertex(pos, val)
      v2 = Vertex(pos, val)
      v3 = Vertex(pos, val)
      s = Simplex([v1, v2, v3])

      defaults = WindingNelderMead.convergenceconfig(dim, T)
      returncode = assessconvergence(s, defaults)
      @test returncode == :XTOL_REACHED
    end

    @testset "Vertices in a simplex eps apart are converged" begin
      dim = 2
      T = Float64
      U = ComplexF64
      v1 = Vertex([one(T), one(T)], rand(U))
      v2 = Vertex([one(T), one(T) + eps(T)], rand(U))
      v3 = Vertex([one(T) + eps(T), one(T)], rand(U))
      s = Simplex([v1, v2, v3])
      defaults = WindingNelderMead.convergenceconfig(dim, T)
      returncode = assessconvergence(s, defaults)
      @test returncode == :XTOL_REACHED
    end

    @testset "Vertices in a chain eps apart are converged" begin
      dim = 2
      T = Float64
      U = ComplexF64
      v1 = Vertex([one(T), one(T)], rand(U))
      v2 = Vertex([one(T), one(T) + eps(T)], rand(U))
      v3 = Vertex([one(T), one(T) + eps(T) + eps(T)], rand(U))
      s = Simplex([v1, v2, v3])
      defaults = WindingNelderMead.convergenceconfig(dim, T)
      returncode = assessconvergence(s, defaults)
      @test returncode == :XTOL_REACHED
      @assert all(isapprox(position(v1)[d], position(v2)[d],
                        atol=defaults[:xtol_abs][d],
                        rtol=defaults[:xtol_rel][d]) for d ∈ 1:dim)
      @assert all(isapprox(position(v2)[d], position(v3)[d],
                        atol=defaults[:xtol_abs][d],
                        rtol=defaults[:xtol_rel][d]) for d ∈ 1:dim)
      @assert !all(isapprox(position(v1)[d], position(v3)[d],
                        atol=defaults[:xtol_abs][d],
                        rtol=defaults[:xtol_rel][d]) for d ∈ 1:dim)
      @test returncode == :XTOL_REACHED
    end

    @testset "Simplex vertices are in order" begin
      for _ ∈ 1:10
        dim = 2
        T = Float64
        U = ComplexF64
        v1 = Vertex(rand(dim), rand(U))
        v2 = Vertex(rand(dim), rand(U))
        v3 = Vertex(rand(dim), rand(U))
        s = Simplex([v1, v2, v3])
        p = s.permabs
        @test abs(value(s[p[1]])) < abs(value(s[p[2]]))
        @test abs(value(s[p[2]])) < abs(value(s[p[3]]))
        @test issortedbyangle(s)
      end
    end

    @testset "Simplex centres" begin
      v1 = Vertex([0.0, 0.0], one(ComplexF64))
      v2 = Vertex([1.0, 0.0], one(ComplexF64))
      v3 = Vertex([0.0, 1.0], one(ComplexF64))
      s = Simplex([v1, v2, v3])
      @test all(centre(s) .== [1/3, 1/3])
    end

    @testset "Simplex extremas" begin
      s = Simplex(x->im, [1.0, 3.0], 1.0)
      exs = WindingNelderMead.extrema(s)
      @test exs[1] == (1.0, 2.0)
      @test exs[2] == (3.0, 4.0)
    end

    @testset "Simplex hypervolumes" begin
      x0 = rand(2)
      a, b = rand(2)
      v1 = Vertex(x0 .+ [0.0, 0.0], one(ComplexF64))
      v2 = Vertex(x0 .+ [a, 0.0], one(ComplexF64))
      v3 = Vertex(x0 .+ [0.0, b], one(ComplexF64))
      s = Simplex([v1, v2, v3])
      @test hypervolume(s) ≈ (a * b) / 2
      x0 = rand(3)
      a, b, c = rand(3)
      v1 = Vertex(x0 .+ [0, 0, 0], one(ComplexF64))
      v2 = Vertex(x0 .+ [a, 0, 0], one(ComplexF64))
      v3 = Vertex(x0 .+ [0, b, 0], one(ComplexF64))
      v4 = Vertex(x0 .+ [0, 0, c], one(ComplexF64))
      s = Simplex([v1, v2, v3, v4])
      @test hypervolume(s) ≈ (a * b * c) / 6
    end

    @testset "atan gives same as angle for sortby" begin
      v1 = Vertex([0.0, 0.0], rand(ComplexF64))
      v2 = Vertex([1.0, 0.0], rand(ComplexF64))
      v3 = Vertex([0.0, 1.0], rand(ComplexF64))
      s = Simplex([v1, v2, v3])
      c = centre(s)
      for v in s
        pv = position(v)
        a = angle(Complex(pv[1] - c[1], pv[2] - c[2]))
        @test a == WindingNelderMead.sortby(s)(v)
      end
    end

    @testset "root (pole) has winding number 1 (-1) in 2D" begin
      objective(x) = (x[1] + im * x[2])
      ps = ([-1.0, -1.0], [1.0, -1.0], [0.0, 1.0])
      s = Simplex([Vertex(p, objective(p)) for p in ps])
      @test windingnumber(s) == 1
      invobjective(x) = 1 / objective(x)
      s = Simplex([Vertex(p, invobjective(p)) for p in ps])
      @test windingnumber(s) == -1
    end

  end

  @testset "End-to-end tests roots" begin
    defaults = WindingNelderMead.convergenceconfig(2, Float64)
    @testset "Single root" begin
      for _ in 1:10
        stopval=10.0^(-rand(3:14))
        root = rand(ComplexF64)
        objective(x) = (x[1] + im * x[2]) - root
        ics = rand(2)
        sizes = rand(2)
        solution = WindingNelderMead.optimise(objective,
          ics, sizes, stopval=stopval, maxiters=10_000)
        (s, n, returncode, its) = solution
        @test n == 1
        if returncode == :STOPVAL_REACHED
          @test abs(value(bestvertex(s))) <= stopval #defaults[:stopval]
        elseif returncode == :XTOL_REACHED && n != 0
          @test isapprox(bestvertex(s).position[1], real(root),
                         rtol=defaults[:xtol_rel][1], atol=defualts[:xtol_abs][1])
          @test isapprox(bestvertex(s).position[2], imag(root),
                         rtol=defaults[:xtol_rel][2], atol=defaults[:xtol_abs][1])
        else
          @show stopval, root, ics, sizes, returncode
          @test false
        end
      end
    end
  end

  @testset "check errors are caught" begin
    @test_throws ArgumentError WindingNelderMead.optimise(x->im, [0.0], [0.0])
    @test_throws ArgumentError WindingNelderMead.optimise(x->im, [0.0], [1.0,
                                                                         2.0])
    @test_throws ArgumentError WindingNelderMead.optimise(x->im, [[0],[1, 2]])
  end

end


text:
using Test
using ForestBiometrics
using StatPlots; gr()

@test typeof(sdi_chart(500,8)) == Plots.Plot{Plots.GRBackend}
@test typeof(sdi_chart(500,8;max_sdi=600)) == Plots.Plot{Plots.GRBackend}


text:
using Printf, Test
using FunctionAccuracyTests

MAX_EXP(n::Val{2}, ::Type{Float32}) = 128.0f0
MAX_EXP(n::Val{2}, ::Type{Float64}) = 1024.0
MAX_EXP(n::Val{:ℯ}, ::Type{Float32}) = 88.72284f0
MAX_EXP(n::Val{:ℯ}, ::Type{Float64}) = 709.7827128933841
MAX_EXP(n::Val{10}, ::Type{Float32}) = 38.53184f0
MAX_EXP(n::Val{10}, ::Type{Float64}) = 308.25471555991675

# min_exp = T(-(exponent_bias(T)+significand_bits(T)) * log(base, big(2)))
MIN_EXP(n::Val{2}, ::Type{Float32}) = -150.0f0
MIN_EXP(n::Val{2}, ::Type{Float64}) = -1075.0
MIN_EXP(n::Val{:ℯ}, ::Type{Float32}) = -103.97208f0
MIN_EXP(n::Val{:ℯ}, ::Type{Float64}) = -745.1332191019412
MIN_EXP(n::Val{10}, ::Type{Float32}) = -45.1545f0
MIN_EXP(n::Val{10}, ::Type{Float64}) = -323.60724533877976

for (func, base) in (exp2=>Val(2), exp=>Val(:ℯ), exp10=>Val(10))
    for T in (Float32, Float64)
        xx = range(MIN_EXP(base,T),  MAX_EXP(base,T), length = 10^6);
        test_acc(func, xx)
    end
end

exp10(Float16(-3.764))
exp10(Float16(-1.76))
asinh(Float16(258.5))
asinh(Float16(-258.5))
@testset "Float16" begin
    @testset "$func" for func in (atan,sinh,cosh,tanh,asinh,
                                  exp,exp2,exp10,expm1,cbrt)
        test_acc(func, -Inf16, Inf16, tol=.501)
    end
    @testset "$func" for func in (sin, cos, tan)
        test_acc(func, -floatmax(Float16),floatmax(Float16), tol=.501)
    end
    @testset "$func" for func in (asin, acos, atanh)
        test_acc(func, -one(Float16),one(Float16), tol=.501)
    end
    test_acc(acosh, one(Float16), Inf16, tol=.501)
    @testset "$func" for func in (log, log2, log10, sqrt)
        test_acc(func, zero(Float16), Inf16, tol=.501)
    end
    test_acc(log1p, -one(Float16), Inf16, tol=.501)
end


text:
module PortfolioAnalytics

using DataFrames
using Dates
using LinearAlgebra
using DataFramesMeta
using Statistics
using LoopVectorization

include("helpers.jl")

include("return_functions.jl")

export plusTwo

plusTwo(x) = return x+2

end # module


text:
struct GebbianPosition <: Position
    x::Float64
    y::Float64
end

struct GebbianEnvironment <: Environment end

get_time(gebbian_model::GebModel) = gebbian_model.time
get_organisms(gebbian_model::GebModel) = [o for o in gebbian_model.organisms]

get_logger(gebbian_model::GebModel) = gebbian_model.logger
function set_logger!(gebbian_model::GebModel, logger::Logger)
    gebbian_model.logger = logger
end

get_id(::GebModel, gebbian_organism::GebOrganism) = gebbian_organism.key
get_genotype(::GebModel, gebbian_organism::GebOrganism) = gebbian_organism.genotype
function get_position(::GebModel, gebbian_organism::GebOrganism)
    GebbianPosition(gebbian_organism.coordinates[1], gebbian_organism.coordinates[2])
end
function get_environment(::GebModel, ::GebOrganism)
    GebbianEnvironment()
end
function get_time_birth(::GebModel, gebbian_organism::GebOrganism)
    gebbian_organism.time_birth
end
function get_time_death(::GebModel, gebbian_organism::GebOrganism)
    gebbian_organism.time_birth + gebbian_organism.age
end
function get_age(::GebModel, gebbian_organism::GebOrganism)
    gebbian_organism.age
end
function clean_snapshot!(model::GebModel)
    model.logger = DoNothingLogger()
end

get_daughters(::GebModel, gebbian_organism::GebOrganism) = gebbian_organism.daughters

function get_abstracted_organism(model::GebModel, model_organism::GebOrganism, parent_ids, parent_genotypes)
    return Organism{GebbianPosition, GebbianEnvironment, String}(
                    get_id(model, model_organism),
                    get_genotype(model, model_organism),
                    parent_ids, parent_genotypes,
                    get_position(model, model_organism),
                    get_environment(model, model_organism),
                    get_time_birth(model, model_organism),
                    get_time_death(model, model_organism),
                    -1
    )
end

function replace_organism!(model::GebModel, key::UInt64, new_genotype::String)
    # remove old organism first

    old_organism = get_organism(model, key)
    kill!(model, old_organism)

    # add new organism

    new_organism = GebOrganism(new_genotype, old_organism.coordinates, old_organism.parent_genotypes[1],
                                old_organism.parent_genotypes[2], model.time)

    add_organism!(model, new_organism)

    return new_organism
end

function run_until!(termination_predicate, model::GebModel)
    while termination_predicate(model) == false
        execute!(model)
    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
module IsingLite

using PyPlot

export neighbors,
       nspins,
       spingrid,
       namefunc,
       heatbath!,
       metropolis!,
       wolff!,
       diagram

include("utils.jl")
include("heatbath.jl")
include("metropolis.jl")
include("wolff.jl")

end # module


text:
function volume(ϵ::Ellipsoid)
    n = size(ϵ.H, 1)
    return (n^(n / 2) * _unit_ball_volume(n)) / sqrt(det(ϵ.H))
end

function _unit_ball_volume(n::Integer)
    return π^(n / 2) / gamma(n / 2 + 1)
end


text:
function calc_Energies( LF::LF3dGrid, Vpot, psi::Array{Float64,2} )
  #
  Ncol = size(psi)[2]
  ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h
  Etot = 0.0
  Ekin = 0.0
  Epot = 0.0
  #
  for ic = 1:Ncol
    nabla2_v = op_nabla2( LF, psi[:,ic] )
    Ekin = Ekin + -0.5*dot( psi[:,ic], nabla2_v ) * ΔV
    Epot = Epot + sum( Vpot .* psi[:,ic].^2 ) * ΔV
  end
  Etot = Ekin + Epot
  return EnergiesT( Etot, Ekin, Epot )
end

function calc_Energies( LF::LF3dGrid, ∇2::SparseMatrixCSC{Float64,Int64},
                        Vpot, psi::Array{Float64,2} )
  #
  Ncol = size(psi)[2]
  ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h
  Etot = 0.0
  Ekin = 0.0
  Epot = 0.0
  #
  for ic = 1:Ncol
    nabla2_v = ∇2 * psi[:,ic]
    Ekin = Ekin + -0.5*dot( psi[:,ic], nabla2_v ) * ΔV
    Epot = Epot + sum( Vpot .* psi[:,ic].^2 ) * ΔV
  end
  Etot = Ekin + Epot
  return EnergiesT( Etot, Ekin, Epot )
end


text:
using MAT

@testset "scaleABCD" begin
    # MATLAB output, generated using output of Julia randn
    test = matread("resources/scaleABCD.mat")

    H = synthesizeNTF(5, 42, 1)
    a, g, b, c = realizeNTF(H)
    ABCD = stuffABCD(a, g, b, c)
    ABCDs, umax, S = scaleABCD(ABCD, 2, 0, 1, NaN, NaN, 10000)
    @test isapprox(test["ABCDs"], ABCDs)
    @test isapprox(test["umax"], umax)
    @test isapprox(test["S"], S)
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
immutable InverseGaussian <: ContinuousUnivariateDistribution
    mu::Float64   # mean
    lambda::Float64 # shape
    function InverseGaussian(mu::Real, lambda::Real)
        mu > zero(mu) || error("mu must be positive")
        lambda > zero(lambda) || error("lambda must be positive")
        new(float64(mu), float64(lambda))
    end
end

InverseGaussian() = InverseGaussian(1.0, 1.0)

insupport(::InverseGaussian, x::Real) = zero(x) <= x < Inf
insupport(::Type{InverseGaussian}, x::Real) = zero(x) <= x < Inf

mean(d::InverseGaussian) = d.mu
mode(d::InverseGaussian) = (r=d.mu/d.lambda; d.mu*(sqrt(1.0+2.25*r*r)-1.5*r))
var(d::InverseGaussian) = d.mu*d.mu*d.mu/d.lambda
skewness(d::InverseGaussian) = 3.0*sqrt(d.mu/d.lambda)
kurtosis(d::InverseGaussian) = 15.0*d.mu/d.lambda

function pdf(d::InverseGaussian, x::Real)
    if x <= 0.0
        return 0.0
    end
    dd = x-d.mu
    sqrt(d.lambda/(twoπ*x*x*x))*exp(-d.lambda*dd*dd/(2.0*d.mu*d.mu*x))
end
function logpdf(d::InverseGaussian, x::Real)
    dd = x-d.mu
    0.5*log(d.lambda/(twoπ*x*x*x))-d.lambda*dd*dd/(2.0*d.mu*d.mu*x)
end

function cdf(d::InverseGaussian, x::Real)
    u = sqrt(d.lambda/x)
    v = x/d.mu
    Φ(u*(v-1.0)) + exp(2.0*d.lambda/d.mu)*Φ(-u*(v+1.0))
end
function ccdf(d::InverseGaussian, x::Real)
    u = sqrt(d.lambda/x)
    v = x/d.mu
    Φc(u*(v-1.0)) - exp(2.0*d.lambda/d.mu)*Φ(-u*(v+1.0))
end
function logcdf(d::InverseGaussian, x::Real)
    u = sqrt(d.lambda/x)
    v = x/d.mu
    a = logΦ(u*(v-1.0))
    b = 2.0*d.lambda/d.mu + logΦ(-u*(v+1.0))
    a + log1pexp(b-a)
end
function logccdf(d::InverseGaussian, x::Real)
    u = sqrt(d.lambda/x)
    v = x/d.mu
    a = logΦc(u*(v-1.0))
    b = 2.0*d.lambda/d.mu + logΦ(-u*(v+1.0))
    a + log1mexp(b-a)
end

function quantile(d::InverseGaussian, p::Real)
    if p <= 0.0 || p >= 1.0
        if p == 1.0
            return inf(Float64)
        elseif p == 0.0
            return 0.0
        else
            return nan(Float64)
        end
    end
    # Whitmore and Yalovsky (1978) approximation
    w = Φinv(p)
    phi = d.lambda/d.mu
    x = d.mu*exp(w/sqrt(phi) - 0.5/phi)
    # one multiplicitive Newton iteration
    xn = x*exp((p-cdf(d,x))/(pdf(d,x)*exp(x)))
    # additive Newton iterations
    while !isapprox(xn,x)
        x = xn
        xn = x + (p-cdf(d,x))/pdf(d,x)
    end
    xn
end

# rand method from:
#   John R. Michael, William R. Schucany and Roy W. Haas (1976)
#   Generating Random Variates Using Transformations with Multiple Roots
#   The American Statistician , Vol. 30, No. 2, pp. 88-90
function rand(d::InverseGaussian)
    z = randn()
    v = z*z
    w = d.mu*v
    x1 = d.mu + d.mu/(2.0*d.lambda)*(w - sqrt(w*(4.0*d.lambda + w)))
    p1 = d.mu / (d.mu + x1)
    u = rand()
    u >= p1 ? d.mu*d.mu/x1 : x1
end


text:
module encoder

using conv_code

export current_out
export encode

function shift_mems!(new_bit::Int64, extended_mems::Array{Int64})
	for k=length(extended_mems):-1:2
		extended_mems[k] = extended_mems[k-1]
	end
	extended_mems[1] = new_bit
end

function current_out(code::ConvolutionalCode, extended_mems::Array{Int64})
	return map(
		indices -> mod(reduce(
			(previous, index) -> previous + extended_mems[index] ,
			indices;
			init=0
		), 2),
		code.outs
	)
end

function encode(code::ConvolutionalCode, info_bits::Array{Int64})
	out_sequence = Array{Array{Int64, 1}, 1}(undef, length(info_bits))
	mems = Array{Int64}(undef, code.quant_mem + 1)
	fill!(mems, 0)

	for index=1:length(info_bits)
		shift_mems!(info_bits[index], mems)
		out_sequence[index] = current_out(code, mems)
	end

	return out_sequence
end

end


text:
"""
    ElasticContactPairMaterial(; coefficientOfRestitution=0.7,
                                 slidingFrictionCoefficient=0.5,
                                 rotationalResistanceCoefficient=0.001,
                                 vsmall=0.01, wsmall=0.01)

Generates an `ElasticContactPairMaterial` object by providing the material properties of two objects that are in contact to each other. Those properties are needed for the [Contact Force Law](@ref).

# Arguments
- `coefficientOfRestitution`: Coefficient of restitution between two objects (=0: inelastic ... =1: fully elastic).
- `slidingFrictionCoefficient`: Kinetic/sliding friction force coefficent between two objects (`>= 0.0`).
- `rotationalResistanceCoefficient`: Rotational resistance torque coefficient between two objects (`>= 0.0`).
  Its effect is that the contact torque is computed to reduce the relative angular velocity between two objects.
  For a ball, `rotationalResistanceCoefficient` is the (standard) rolling resistance coefficient.
- `vsmall` in [m/s]: Used for regularization when computing the unit vector in direction of the relative tangential velocity to avoid a division by zero.
- `wsmall` in [rad/s]: Used for regularization when computing the unit vector in direction of the relative angular velocity to avoid a division by zero.
"""
mutable struct ElasticContactPairMaterial <: Modia3D.AbstractContactPairMaterial
    coefficientOfRestitution::Float64
    slidingFrictionCoefficient::Float64
    rotationalResistanceCoefficient::Float64
    vsmall::Float64
    wsmall::Float64

    function ElasticContactPairMaterial(; coefficientOfRestitution=0.0,
                                          slidingFrictionCoefficient=0.0,
                                          rotationalResistanceCoefficient=0.0,
                                          vsmall=0.01,
                                          wsmall=0.01)
        @assert(coefficientOfRestitution >= 0.0)
        @assert(coefficientOfRestitution <= 1.0)
        @assert(slidingFrictionCoefficient >= 0.0)
        @assert(rotationalResistanceCoefficient >= 0.0)
        new(coefficientOfRestitution, slidingFrictionCoefficient,
            rotationalResistanceCoefficient, vsmall, wsmall)
end; end

struct NoContactPairMaterial <: Modia3D.AbstractContactPairMaterial
end

mutable struct ObserverContactPairMaterial <: Modia3D.AbstractContactPairMaterial
    printAlarm::Bool

    function ObserverContactPairMaterial(; printAlarm=false)
        new(printAlarm)
end; end

mutable struct ImpulseContactPairMaterial <: Modia3D.AbstractContactPairMaterial
end

mutable struct WheelRailContactPairMaterial <: Modia3D.AbstractContactPairMaterial
end

@enum ResponseType ElasticResponse ObserverResponse ImpulseResponse WheelRailResponse NoResponse

function getResponseMaterial(responseType::String)
    if responseType == "ElasticResponse"
        return ElasticContactPairMaterial()
    end
    if responseType == "ObserverResponse"
        return ObserverContactPairMaterial()
    end
    if responseType == "ImpulseResponse"
        return ImpulseContactPairMaterial()
    end
    if responseType == "WheelRailResponse"
        return WheelRailContactPairMaterial()
    end
    if responseType == "NoResponse"
        return NoContactPairMaterial()
    end
end

"""
    material = combineContactPairMaterials(material1, material2)

Return a new `ElasticContactPairMaterial` object from the mean values of the ElasticContactPairMaterial
objects `material1` and `material2`.
"""
function combineContactPairMaterials(mat1::ElasticContactPairMaterial, mat2::ElasticContactPairMaterial)::ElasticContactPairMaterial
    mat = ElasticContactPairMaterial()
    mat.coefficientOfRestitution        = (mat1.coefficientOfRestitution + mat2.coefficientOfRestitution)/2
    mat.slidingFrictionCoefficient      = (mat1.slidingFrictionCoefficient + mat2.slidingFrictionCoefficient)/2
    mat.rotationalResistanceCoefficient = (mat1.rotationalResistanceCoefficient + mat2.rotationalResistanceCoefficient)/2
    mat.vsmall = (mat1.vsmall + mat2.vsmall)/2
    mat.wsmall = (mat1.wsmall + mat2.wsmall)/2
    return mat
end

struct TwoNamesKey
    name1::String
    name2::String
    function TwoNamesKey(name1::String, name2::String)
        if name1 < name2
            return new(name1, name2)
        else
            return new(name2, name1)
end; end; end

Base.:isequal(keyA::TwoNamesKey, keyB::TwoNamesKey) = keyA.name1 == keyB.name1 && keyA.name2 == keyB.name2
Base.:isless(keyA::TwoNamesKey, keyB::TwoNamesKey) = keyA.name1 < keyB.name1 || ( keyA.name1 == keyB.name1 && keyA.name2 < keyB.name2 )

"""
    readContactPairMaterialFromJSON(filename)

Read a JSON file consisting of a dictionary of ContactPairMaterial instances from `filename` and
return a `Dict{String, ContactPairMaterial}` dictionary.
"""
function readContactPairMaterialFromJSON(filename::AbstractString)
    dataJSON = JSON.parsefile(filename)
    palette = Dict{TwoNamesKey, Modia3D.AbstractContactPairMaterial}()
    for (namesPair,dataPair) in dataJSON
        # Split a key "name1,name2" in two parts
        i = findfirst(",", namesPair)
        if isnothing(i) || i.start < 2 || i.start >= length(namesPair)
            println("... Wrong key \"", namesPair,"\" not included in contactPairMaterialPalette.")
        else
            namesKey = TwoNamesKey( namesPair[1:i.start-1], namesPair[i.start+1:end] )
            responseMat = getResponseMaterial(pop!(dataPair, "responseType", NoResponse))
            for (propName,propData) in dataPair
                if isdefined(responseMat, Symbol(propName))
                    setfield!(responseMat, Symbol(propName), propData)
                else
                    error("... from key: ", namesKey, " . There might be spelling mistakes. Check ", propName, ".")
                end
            end
            palette[namesKey] = responseMat
    end; end
    return palette
end

"""
    contactPairMaterialPalette

Dictionary of contact pair material data, see [`ElasticContactPairMaterial`](@ref)
"""
contactPairMaterialPalette = readContactPairMaterialFromJSON( joinpath(Modia3D.path, "palettes", "contactPairMaterials.json") )

"""
    material = getContactPairMaterial(name1, name2)

Return a contact pair material object of the materials with the
names `name1` and `name2` from dictionary `Modia3D.contactPairMaterialPalette`.
"""
function getContactPairMaterial(name1::AbstractString, name2::AbstractString )::Modia3D.AbstractContactPairMaterial
    value = get(contactPairMaterialPalette, TwoNamesKey(name1,name2), NoContactPairMaterial())
    if typeof(value) == NoContactPairMaterial
        # Combination name1,name2 is not present
        value1 = get(contactPairMaterialPalette, TwoNamesKey(name1,name1), NoContactPairMaterial())
        value2 = get(contactPairMaterialPalette, TwoNamesKey(name2,name2), NoContactPairMaterial())
        if typeof(value1) == NoContactPairMaterial || typeof(value2) == NoContactPairMaterial
            error("No contact pair material (\"$name1\",\"$name2\") or\n",
                  "(\"$name1\",\"$name1\") or (\"$name2\",\"$name2\")\n",
                  "in Modia3D.contactPairMaterialPalette dictionary.")
        end
        return combineContactPairMaterials(value1, value2)
    end
    return value
end


text:
module Sim
include("./name2idx/parameters.jl")
include("./name2idx/species.jl")
include("./set_model.jl")
include("./observable.jl")

using .C
using .V

using DelayDiffEq

normalization = Dict{String,Dict{}}()
for observable in observables
    normalization[observable] = Dict(
        "timepoint" => nothing,
        "condition" => ["WT"]
    )
end

const dt = 1.0

t = collect(0.0:1.0:360.0)  # 0, 1, 2, ..., 360 [min.]

const sstime = 1000.0  # time to reach steady state

const conditions = ["WT"]

simulations = Array{Float64,3}(
    undef, length(observables), length(t), length(conditions)
)

function solvedde(
        f::Function,u0::Vector{Float64},history::Vector{Float64},
        tspan::Tuple{Float64,Float64},p::Vector{Float64},tau::Float64)
    h(p, t) = history
    lags = [tau]
    prob = DDEProblem(f, u0, h, tspan, p; constant_lags=lags)
    alg = MethodOfSteps(BS3())
    sol = solve(
        prob,alg,saveat=dt,abstol=1e-8,reltol=1e-8,verbose=false
    )
    return sol
end

function get_steady_state(
        p::Vector{Float64},u0::Vector{Float64},
        sstime::Float64,tau::Float64)::Vector{Float64}
    # get steady state (t<0)
    p[C.term] = 1.0
    history::Vector{Float64} = u0
    tspan::Tuple{Float64,Float64} = (0.0, sstime)
    try
        sol = solvedde(diffeq!, u0, history, tspan, p, tau)
        if sol.retcode === :Success
            return sol[:, end]
        else
            return []
        end
    catch
        return []
    end
end

function get_time_course(
        p::Vector{Float64},u0::Vector{Float64},
        sstime::Float64,tau::Float64)
    p1::Vector{Float64} = copy(p)
    p1[C.term] = 0.0
    u1::Vector{Float64} = get_steady_state(p, u0, sstime, tau)
    if isempty(u1)
        return nothing
    end
    history::Vector{Float64} = u1
    tspan::Tuple{Float64,Float64} = (0.0, t[end])
    try
        sol = solvedde(diffeq!, u1, history, tspan, p1, tau)
        return ifelse(sol.retcode === :Success, sol, nothing)
    catch
        return nothing
    end
end

function simulate!(
        p::Vector{Float64},
        u0::Vector{Float64})::Union{Bool,Nothing}
    for (i, condition) in enumerate(conditions)
        # if condition == "WT"
        #    pass
        # end
        sol = get_time_course(p, u0, sstime, p[C.delayrnae])
        if sol === nothing
            return false
        else
            @inbounds @simd for j in eachindex(t)
                simulations[observables_index("Nuclear_NFkB"),j,i] = (
                    sol.u[j][V.NFKBn]
                )
            end
        end
    end
end
end # module


text:
module YStockData

using Requests
using Base.Dates
using DataFrames

const fields = Dict(
"1_year_target" => "t8",
"200_sma" => "m4",
"50_sma" => "m3",
"52_week_high" => "k",
"52_week_low" => "j",
"52_week_range" => "w",
"after_hours_change" => "c8",
"annualized_gain" => "g3",
"ask_realtime" => "b2",
"ask_size" => "a5",
"average_daily_volume" => "a2",
"bid_realtime" => "b3",
"bid_size" => "b6",
"book_value" => "b4",
"change" => "c1",
"change_200_sma" => "m5",
"change_50_sma" => "m7",
"change_from_52_week_high" => "k4",
"change_from_52_week_low" => "j5",
"change_percent" => "p2",
"change_percent_change" => "c",
"change_percent_realtime" => "k2",
"change_realtime" => "c6",
"commission" => "c3",
"company_name" => "n",
"dividend_pay_date" => "r1",
"dividend_per_share" => "d",
"dividend_yield" => "y",
"ebitda" => "j4",
"eps" => "e",
"eps_estimate_current_year" => "e7",
"eps_estimate_next_quarter" => "e9",
"eps_estimate_next_year" => "e8",
"ex_dividend_date" => "q",
"float_shares" => "f6",
"high_limit" => "l2",
"holdings_gain" => "g4",
"holdings_gain_percent" => "g1",
"holdings_gain_percent_realtime" => "g5",
"holdings_gain_realtime" => "g6",
"holdings_value" => "v1",
"holdings_value_realtime" => "v7",
"last_trade_date" => "d1",
"last_trade_price" => "l1",
"last_trade_realtime_time" => "k1",
"last_trade_size" => "k3",
"last_trade_time" => "t1",
"last_trade_time_plus" => "l",
"low_limit" => "l3",
"market_cap" => "j1",
"market_cap_realtime" => "j3",
"more_info" => "v",
"notes" => "n4",
"order_book_realtime" => "i5",
"pe" => "r",
"pe_realtime" => "r2",
"peg" => "r5",
"percent_change_200_sma" => "m6",
"percent_change_50_sma" => "m8",
"percent_change_from_52_week_high" => "k5",
"percent_change_from_52_week_low" => "j6",
"previous_close" => "p",
"price_book" => "p6",
"price_eps_estimate_current_year" => "r6",
"price_eps_estimate_next_year" => "r7",
"price_paid" => "p1",
"price_sales" => "p5",
"revenue" => "s6",
"shares_outstanding" => "j2",
"shares_owned" => "s1",
"short_ratio" => "s7",
"stock_exchange" => "x",
"ticker_trend" => "t7",
"today_open" => "o",
"todays_high" => "h",
"todays_low" => "g",
"todays_range" => "m",
"todays_range_realtime" => "m2",
"todays_value_change" => "w1",
"todays_value_change_realtime" => "w4",
"trade_date" => "d2",
"trade_links" => "t6",
"volume" => "v"
)

const inverted_fields = Dict{ASCIIString, ASCIIString}()
for (n,v) in fields
  sn = symbol("get_"*n)
  @eval begin
      function $(sn)(sym)
        get(sym, $v)
      end
  end
  inverted_fields[v] = n
end

get(sym::ByteString, fld::ByteString) =
  to_nullable(request_yahoo(sym, fld))
get(sym::ByteString, fld::Vector) =
  to_nullable_dict(fld, split(request_yahoo(sym, reduce(*, fld)), ","))
get(sym::ByteString) = get(sym, collect(values(fields)))

function request_yahoo(sym, fld)
  url = "http://finance.yahoo.com/d/quotes.csv?s=$sym&f=$fld"
  resp = Requests.get(url)
  r = readall(resp)
  return r
end

function get_historical(sym, startd, endd)

    q = Dict(
     "s" => sym,
     "a" => month(startd) - 1,
     "b" => dayofmonth(startd),
     "c" => year(startd),
     "d" => month(endd) - 1,
     "e" => dayofmonth(endd),
     "f" => year(endd),
     "g" => "d",
     "ignore" => ".csv"
    )
    url = "http://real-chart.finance.yahoo.com/table.csv"
    resp = Requests.get(url; query=q)
    return readtable(IOBuffer(readall(resp)))
end

function get_historical(sym::Vector, startd, endd)
  r = DataFrame()
  println("Fetching historical data from Yahoo \nfrom $startd to $endd for $(length(sym)) tickers")
  for (i,x) in enumerate(sym)
    println("$i of $(length(sym)) : $x")
    df=get_historical(x, startd, endd)
    df[:Ticker] = repmat([x], size(df, 1))
    if length(r) == 0
      r = df
    else
      r = vcat(r, df)
    end
  end
  return r
end

function to_nullable_dict{T<:AbstractString}(fld::Vector, x::Vector{T})
  r = Dict{ASCIIString, Nullable}()
  for (i, y)  in enumerate(x)
    r[inverted_fields[fld[i]]] = to_nullable(y)
  end
  return r
end

function to_nullable{T<:AbstractString}(x::T)
  try
    return Nullable(float(x))
  catch
    if x == "N/A"
      return Nullable{Float64}()
    else
      return Nullable(strip(replace(x, '\"', ' ')))
    end
  end
end

end # module


text:
using Flux
using Random

mutable struct HedgeModel{M, F}
    h::M
    f::F
end

function HedgeModel(h::M; init=Flux.glorot_uniform) where {M<:Chain}
    out = size(h[end].W)[1]
    f = [Dense(size(h[1].W)[2], out; init=init)]
    for i in 1:length(h)
        push!(f, Dense(size(h[1.W])[1], out; init=init))
    end
    HedgeModel(model, f)
end

mutable struct HedgeDQNAgent{M, TN, O, LU, AP<:AbstractValuePolicy, Φ, ERP} <: AbstractAgent
    model::M
    target_network::TN
    opt::O
    lu::LU
    ap::AP
    er::ExperienceReplay{ERP}
    γ::Float32
    batch_size::Int64
    tn_counter_init::Int64
    target_network_counter::Int64
    action::Int64
    prev_s::Φ
end

HedgeDQNAgent(model, target_network, opt, lu, ap, size_buffer, γ, batch_size, tn_counter_init, s) =
    HedgeDQNAgent(model,
                  target_network,
                  opt,
                  lu,
                  ap,
                  ExperienceReplay(100000,
                                   (Array{Float32, 1}, Int64, Array{Float32, 1}, Float32, Bool),
                                   (:s, :a, :sp, :r, :t)),
                  γ,
                  batch_size,
                  tn_counter_init,
                  tn_counter_init,
                  0,
                  s)

function JuliaRL.start!(agent::HedgeDQNAgent, env_s_tp1, rng::AbstractRNG; kwargs...)
    # Start an Episode
    agent.action = sample(agent.ap,
                          agent.model(env_s_tp1),
                          rng)
    return agent.action
end

function JuliaRL.step!(agent::HedgeDQNAgent, env_s_tp1, r, terminal, rng::AbstractRNG; kwargs...)

    add!(agent.er, (Float32.(agent.prev_s), agent.action, copy(Float32.(env_s_tp1)), r, terminal))

    if size(agent.er)[1] > 1000
        e = sample(agent.er, agent.batch_size; rng=rng)
        update_params!(agent, e)
    end

    agent.prev_s .= env_s_tp1
    agent.action = sample(agent.ap,
                          agent.model(agent.prev_s),
                          rng)

    return agent.action
end

function update_params!(agent::HedgeDQNAgent, e)

    if agent.tn_counter_init > 0
        update!(agent.model, agent.lu, agent.opt, e.s, e.a, e.sp, e.r, e.t, agent.target_network)

        if agent.target_network_counter == 1
            agent.target_network_counter = agent.tn_counter_init
            agent.target_network = mapleaves(
                Flux.Tracker.data,
                deepcopy(agent.model))
        else
            agent.target_network_counter -= 1
        end
    else
        update!(agent.model, agent.lu, agent.opt, e.s, e.a, e.sp, e.r, e.t)
    end
    return nothing

end


text:
module TopOptLevelSet

using TupleTools, StaticArrays, Reexport, SparseArrays, LinearAlgebra
using ChainRulesCore, ForwardDiff, SparseDiffTools
@reexport using Meshes
const TOLS = TopOptLevelSet

export TOLS

include("wedge.jl")
include("fractions.jl")
include("levelset.jl")

end


text:
function Logistic(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble)
    return Logistic((jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5)
end

function value(obj::Logistic, arg0::DerivativeStructure)
    return jcall(obj, "value", DerivativeStructure, (DerivativeStructure,), arg0)
end

function value(obj::Logistic, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text:
# This file is part of EndemicModel project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

abstract type AbstractEndemicModel end

variables(model::AbstractEndemicModel) = model.variables
variables!(model::AbstractEndemicModel, value) = model.variables = value

derivatives(model::AbstractEndemicModel) = model.derivatives
derivatives!(model::AbstractEndemicModel, value) = model.derivatives = value

parameters(model::AbstractEndemicModel) = model.parameters
parameters!(model::AbstractEndemicModel, value) = model.parameters = value

modeldata(model::AbstractEndemicModel) = model.modeldata
modeldata!(model::AbstractEndemicModel, value::AbstractDataFrame) = model.modeldata = value

modeldata!(model::AbstractEndemicModel; kwargs...) =
    modeldata!(model::AbstractEndemicModel, to_dataframe(model; kwargs...))

default_kwargs(model::AbstractEndemicModel) = model.kwargs
default_kwargs!(model::AbstractEndemicModel, value) = model.kwargs = value

default_kwarg(model::AbstractEndemicModel, key::Symbol) = model.kwargs[key]
default_kwarg(model::AbstractEndemicModel, key::Symbol, default) =
    get(model.kwargs, key, default)
default_kwarg!(model::AbstractEndemicModel, key::Symbol, value) = model.kwargs[key] = value

realdata(model::AbstractEndemicModel) = model.realdata
realdata!(model::AbstractEndemicModel, value) = model.realdata = value

groupnames(model::AbstractEndemicModel) = model.groupnames
groupnames!(model::AbstractEndemicModel, value) = model.groupnames = value

ngroups(model::AbstractEndemicModel) = model.ngroups
ngroups!(model::AbstractEndemicModel, value) = model.ngroups = value

export_data(model::AbstractEndemicModel; kwargs...) =
    export_data(modeldata(model); kwargs...)

export_data(source, model::AbstractEndemicModel; kwargs...) =
    export_data(source, modeldata(model); kwargs...)

const SEIR_VARIABLES = (:S, :E, :I, :R)
const SEIR_DERIVATIVES = (:dS, :dE, :dI, :dR)
const SEIR_PARAMETERS = (:M, :β, :γ, :α, :μ)
const SEIRVariables = NamedTuple{SEIR_VARIABLES,NTuple{4,Vector{Float}}}
const SEIRDerivatives = NamedTuple{SEIR_DERIVATIVES,NTuple{4,Vector{Float}}}
const SEIRParameters = NamedTuple{SEIR_PARAMETERS,NTuple{5,Vector{Float}}}

SEIRVariables(S, E, I, R) = SEIRVariables((S, E, I, R))
SEIRDerivatives(dS, dE, dI, dR) = SEIRDerivatives((dS, dE, dI, dR))
SEIRParameters(M, β, γ, α, μ) = SEIRParameters((M, β, γ, α, μ))

function _SEIR_derivative(vars::SEIRVariables, params::SEIRParameters)
    (S, E, I, R) = vars
    (M, β, γ, α, μ) = params

    _SEIR_derivative(S, E, I, R, M, β, γ, α)
end

"""
Compute the derivative of a SEIR-type differential equation.
In this system, `exposed` people are also considered infections, yet not tested confirmed.
`S`, `E`, `I`, `R` must be vectors of the same length `n`.
`M`, `γ` must be vectors of length `n`.
`β` must be Float64.

Equations:
M = S + E + I + R
dM = 0
dS = - β * (I ./ sum(M)) .* S
dE = β * (I ./ sum(M)) .* S - γ .* E
dI = γ .* E - α .* I
dR = α .* I
"""
function _SEIR_derivative(
    S::AbstractVector{Float},
    E::AbstractVector{Float},
    I::AbstractVector{Float},
    R::AbstractVector{Float},
    M::AbstractVector{Float},
    β::AbstractVector{Float},
    γ::AbstractVector{Float},
    α::AbstractVector{Float};
    backward_limit::Bool = false,
)
    pop = sum(M)
    StoE = β .* sum(I ./ pop) .* S
    EtoI = γ .* E
    ItoR = α .* I
    if backward_limit
        StoE = min.(E, StoE)
        EtoI = min.(I, EtoI)
        ItoR = min.(R, ItoR)
    end

    dS = -StoE
    dE = StoE - EtoI
    dI = EtoI - ItoR
    dR = ItoR

    # @debug "Computed SEIR derivatives" (S, E, I, R) (dS, dE, dI, dR)

    (; dS = dS, dE = dE, dI = dI, dR = dR)
end

initial_M(kind::Symbol; kwargs...) = initial_M(Val(kind); kwargs...)

function initial_M(kind::Val{:nogroup}; estimated_population = 1e6, kwargs...)
    Float[estimated_population]
end

"""
Represent a SEIR differential equation.

S, E, I, R must be vectors of the same length `n`.
M, β, γ, α, μ must be vectors of length `n`.

Variables:
S = susceptible
E = exposed (infected, but not infeccious)
I = active (infeccious, confirmed)
R = recovered / deaths

Parameters:
M = total susceptible population
β = transmission rate
γ = evolution rate, i.e., E -> I rate
α = recovery/death rate, i.e., Ic -> R rate
μ = death rate

Equations:
M = S + E + I + R
dM = 0
dS = - β * (I ./ M) .* S
dE = β * (I ./ M) .* S - γ .* E
dI = γ .* E - α .* I
dR = α .* I
"""
mutable struct SEIRModel <: AbstractEndemicModel
    variables::SEIRVariables
    parameters::SEIRParameters
    derivatives::SEIRDerivatives
    ngroups::Int
    groupnames::Vector{Symbol}
    realdata::Union{Nothing,AbstractDataFrame}
    modeldata::Union{Nothing,AbstractDataFrame}
    kwargs::Dict{Symbol,Any}

    function SEIRModel(
        vars::SEIRVariables,
        params::SEIRParameters;
        ngroups::Int = length(vars[1]),
        groupnames = nothing,
        realdata::Union{Nothing,AbstractDataFrame} = nothing,
        modeldata::Union{Nothing,AbstractDataFrame} = nothing,
        kwargs...,
    )
        groupnames =
            something(groupnames, ngroups == 1 ? [:all] : "group_" .* string.(1:ngroups))
        groupnames = collect(Symbol.(groupnames))
        deriv = _SEIR_derivative(vars, params)
        info = Dict{Symbol,Any}(kwargs)

        model = new(vars, params, deriv, ngroups, groupnames, realdata, modeldata, info)
        if isnothing(modeldata)
            modeldata!(model, to_dataframe(model; kwargs...))
        end
        model
    end
end

function Base.copy(model::SEIRModel)
    SEIRModel(
        variables(model),
        parameters(model);
        ngroups = ngroups(model),
        groupnames = groupnames(model),
        realdata = realdata(model),
        modeldata = modeldata(model),
        default_kwargs(model)...,
    )
end

SEIRModel(model::SEIRModel) = model

function SEIRModel(
    S::AbstractVector{Float},
    E::AbstractVector{Float},
    I::AbstractVector{Float},
    R::AbstractVector{Float},
    M::AbstractVector{Float},
    β::AbstractVector{Float},
    γ::AbstractVector{Float},
    α::AbstractVector{Float},
    μ::AbstractVector{Float};
    kwargs...,
)
    SEIRModel(SEIRVariables(S, E, I, R), SEIRParameters(M, β, γ, α, μ); kwargs...)
end

function model_step(model::SEIRModel, n::Int = 1; initial_date = nothing, kwargs...)
    (S, E, I, R) = variables(model)
    (dS, dE, dI, dR) = derivatives(model)
    params = parameters(model)
    (M, β, γ, α, μ) = params
    initial_date = something(
        initial_date,
        default_kwarg(model, :initial_date, today() - Day(15)) + Day(n),
    )
    data = realdata(model)
    if n > 0
        while n > 0
            for i ∈ 1:4
                (S, E, I, R) = (S + 0.25dS, E + 0.25dE, I + 0.25dI, R + 0.25dR)
                (dS, dE, dI, dR) = _SEIR_derivative(S, E, I, R, M, β, γ, α)
            end
            n -= 1
        end
    else
        z = zeros(Float, length(S))
        if !isnothing(data)
            idx = Dates.days(initial_date - data.date[1]) + 1
            if idx >= 1
                n = 0 # skip loop below
            else
                n, idx = idx - 1, 1
            end
            E = [data.exposed[idx]]
            I = [data.active[idx]]
            R = [data.closed[idx]]
            S = [data.estimated_population[idx] - E[1] - I[1] - R[1]]
        end
        while n < 0
            for i ∈ 1:4
                (S, E, I, R) = (S - 0.25dS, E - 0.25dE, I - 0.25dI, R - 0.25dR)
                (dS, dE, dI, dR) =
                    _SEIR_derivative(S, E, I, R, M, β, γ, α; backward_limit = true)
            end
            n += 1
        end

        # optimize_variables!(model; kwargs...)
        modeldata!(model, to_dataframe(model; kwargs...))
    end
    SEIRModel(
        SEIRVariables(S, E, I, R),
        params;
        initial_date = initial_date,
        realdata = data,
        kwargs...,
    )
end

function pack_vars(S, E, I, R)
    vectors = zip(S, E, I, R)
    result = Float[]
    for vec ∈ vectors
        append!(result, vec)
    end
    result
end

function pack_vars(tup::Union{Tuple,NamedTuple})
    pack_vars(tup...)
end

function pack_params(M, β, γ, α, μ)
    args = filter(!isnothing, (M, β, γ, α, μ))
    vectors = zip(args...)
    result = Float[]
    for v ∈ vectors
        append!(result, v)
    end
    result
end

function pack_params(tup::Tuple)
    pack_params(tup...)
end

function pack_params(tup::NamedTuple; packed_params = keys(tup))
    if packed_params == SEIR_PARAMETERS
        pack_params(tup...)
    else
        M = :M ∈ packed_params ? tup[:M] : nothing
        β = :β ∈ packed_params ? tup[:β] : nothing
        γ = :γ ∈ packed_params ? tup[:γ] : nothing
        α = :α ∈ packed_params ? tup[:α] : nothing
        μ = :μ ∈ packed_params ? tup[:μ] : nothing
        pack_params(M, β, γ, α, μ)
    end
end

function unpack_vars(X::Vector{Float})
    nv = 4
    Xlen = length(X)
    S = view(X, 1:nv:Xlen)
    E = view(X, 2:nv:Xlen)
    I = view(X, 3:nv:Xlen)
    R = view(X, 4:nv:Xlen)
    SEIRVariables(S, E, I, R)
end

function unpack_params(
    P::Vector{Float},
    ngroups::Int;
    packed_params = SEIR_PARAMETERS,
    default_params = (),
)
    vecvars = [:M, :β, :γ, :α, :μ] ∩ packed_params
    Plen = length(vecvars) * ngroups
    index = 1
    np = length(vecvars)
    if :M ∈ vecvars
        M = view(P, index:np:Plen)
        index += 1
    else
        M = default_params[:M]
    end
    if :β ∈ vecvars
        β = view(P, index:np:Plen)
        index += 1
    else
        β = default_params[:M]
    end
    if :γ ∈ vecvars
        γ = view(P, index:np:Plen)
        index += 1
    else
        γ = default_params[:γ]
    end
    if :α ∈ vecvars
        α = view(P, index:np:Plen)
        index += 1
    else
        α = default_params[:α]
    end
    if :μ ∈ vecvars
        μ = view(P, index:np:Plen)
        index += 1
    else
        μ = default_params[:μ]
    end
    SEIRParameters(M, β, γ, α, μ)
end

function SEIR_ODE_fun(X, P, t)
    (S, E, I, R) = unpack_vars(X)
    n = length(S)
    (M, β, γ, α, μ) = unpack_params(P, n)

    (dS, dE, dI, dR) =
        _SEIR_derivative(SEIRVariables(S, E, I, R), SEIRParameters(M, β, γ, α, μ))
    pack_vars(dS, dE, dI, dR)
end

function SEIR_ODEProblem(
    inivars::SEIRVariables,
    params::SEIRParameters;
    maxtime = 360,
    kwargs...,
)
    ODEProblem(SEIR_ODE_fun, pack_vars(inivars), (0.0, Float(maxtime)), pack_params(params))
end

function model_problem(model::SEIRModel; kwargs...)
    SEIR_ODEProblem(variables(model), parameters(model); kwargs...)
end

function model_solution(model::SEIRModel; dt = 1.0, kwargs...)
    solve(model_problem(model; kwargs...), Euler(); dt = dt)
end

function to_dataframe(model::SEIRModel; diff_columns::Bool = true, kwargs...)
    solution = model_solution(model; dt = 0.25, kwargs...)
    (M, β, γ, α, μ) = parameters(model)
    days = solution.t[1:4:end]
    solution = solution[1:4:end]
    initial_date = default_kwarg(model, :initial_date, today() - Day(15))
    date = initial_date .+ Day.(days)
    n = length(days)
    ng = model.ngroups
    exposed = zeros(Int, n)
    active = zeros(Int, n)
    recovered = zeros(Int, n)
    deaths = zeros(Int, n)
    infected = zeros(Int, n)
    if diff_columns
        diff_exposed = zeros(Int, n)
        diff_active = zeros(Int, n)
        diff_recovered = zeros(Int, n)
        diff_deaths = zeros(Int, n)
    end

    _round(x) = isnan(x) ? 0 : round(Int, max(0, min(1e12, x)))

    for i ∈ 1:n
        (S, E, I, R) = unpack_vars(solution[i])
        (dS, dE, dI, dR) = _SEIR_derivative(S, E, I, R, M, β, γ, α)
        exposed[i] = sum(E) |> _round
        active[i] = sum(I) |> _round
        recovered[i] = sum(R .* (1.0 .- μ)) |> _round
        deaths[i] = sum(R .* μ) |> _round
        infected[i] = active[i] + recovered[i] + deaths[i] |> _round
        if diff_columns
            diff_exposed[i] = sum(dE) |> _round
            diff_active[i] = sum(dI) |> _round
            diff_recovered[i] = sum(dR .* (1.0 .- μ)) |> _round
            diff_deaths[i] = sum(dR .* μ) |> _round
        end
    end
    dcols = if diff_columns
        [
            :diff_exposed => diff_exposed,
            :diff_active => diff_active,
            :diff_recovered => diff_recovered,
            :diff_deaths => diff_deaths,
        ]
    else
        []
    end
    DataFrame(
        :date => date,
        :infected => infected,
        :exposed => exposed,
        :active => active,
        :recovered => recovered,
        :deaths => deaths,
        dcols...;
        copycols = false,
    )
end


text:
@testset "Spin Sampling" begin
    function test_sampling_spin(
        sampling_type::Symbol, pg::ParticleGroup{D,V}, df::CosSumGaussianSpin
    ) where {D,V}
        mean = zeros(2)
        sigma = zeros(2)

        n_particles = pg.n_particles

        sampling = ParticleSampler{D,V}(sampling_type, n_particles)

        sample!(pg, sampling, df, mesh)

        for i_part in 1:n_particles
            xi = get_x(pg, i_part)
            vi = get_v(pg, i_part)
            mean[1] += xi[1]
            mean[2] += vi[1]
        end

        mean = mean / n_particles

        for i_part in 1:n_particles
            xi = get_x(pg, i_part)
            vi = get_v(pg, i_part)
            sigma[1] += (xi[1] - mean[1])^2
            sigma[2] += (vi[1] - mean[2])^2
        end

        sigma = sigma / (n_particles - 1)

        return mean, sigma
    end

    σ, μ = 0.02, 0.0
    kx, α = 1.004355, 0.001
    n_particles = 100000
    xmin = 0.0
    xmax = 2π / kx
    Lx = xmax - xmin
    nx = 512

    mesh = OneDGrid(xmax, xmin, nx)

    pg = ParticleGroup{1,1}(n_particles; n_spin=3)

    params = (k=[[kx]], α=[α], σ=[[σ]], μ=[[μ]])

    df = CosSumGaussianSpin(params...)

    mean, sigma = test_sampling_spin(:sobol, pg, df)

    @test true
end


text:
"""
    function kitty()

A kitten of sorts, I think

# Example

```jldoctest
julia> cowsay("Meow", cow=Cowsay.kitty)
 ______
< Meow >
 ------
     \\
      \\
       ("`-'  '-/") .___..--' ' "`-._
         ` *_ *  )    `-.   (      ) .`-.__. `)
         (_Y_.) ' ._   )   `._` ;  `` -. .-'
      _.. `--'_..-_/   /--' _ .' ,4
   ( i l ),-''  ( l i),'  ( ( ! .-'
```
"""
function kitty(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
     $thoughts
      $thoughts
       ("`-'  '-/") .___..--' ' "`-._
         ` *_ *  )    `-.   (      ) .`-.__. `)
         (_Y_.) ' ._   )   `._` ;  `` -. .-'
      _.. `--'_..-_/   /--' _ .' ,4
   ( i l ),-''  ( l i),'  ( ( ! .-'
"""
return the_cow
end


text:
using MusicXML

# Reads musicxml file and then extracts the data, builds all the types and stores them in proper format.
doc = readmusicxml(joinpath("examples", "musescore.musicxml"))

# Example 1:
# Prints Each instrument name and then the pitches

# Extracting each instrument information
scprts = doc.scorepartwise.partlist.scoreparts

# Extracting parts
prts = doc.scorepartwise.parts

# Extracting each part
for prt in prts

    ind = findfirst(x -> prt.id == x.id, scprts) # returns the index of scorepart that matches the id of part

    # printing the instrument name
    println(scprts[ind].name)

    # Extracting each measure of the part
    iMsr=1
    for msr in prt.measures
        # Extracting notes of each measure
        for nt in msr.notes
            if !isnothing(nt.pitch)

                println("Measure no. $iMsr")    # print measure number
                println(nt.pitch)     # print pitch of the note
                println(nt.duration)  # print duration of the note
            elseif !isnothing(nt.unpitched)
                println("Measure no. $iMsr")    # print measure number
                println(nt.unpitched) # print unpitched of the note
                println(nt.duration)  # print duration of the note
            end

        end

        iMsr+=1
    end
end


text:
#!/usr/bin/julia

rho0 = [1 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 0]
mua = [0 1 1 0; 1 0 0 0; 1 0 0 0; 0 0 0 0]
mub = [0 0 0 0; 0 0 0 1; 0 0 0 1; 0 1 1 0]

println("rho0")
display(rho0)
print("mu g-e")
display(mua)
print("mu e-f")
display(mub)

println("GSB")
println("R")
println("step0")
display(rho0)
println("step1")
display(rho0 * mua)
println("step2")
display((rho0 * mua) * mua)
println("step3")
display(mua * ((rho0 * mua) * mua))

max_stars_repo_name: lbmnky/CMDS

max_issues_repo_name: lbmnky/CMDS

max_forks_repo_name: lbmnky/CMDS


text:
"""
    abstract type AbstractCredentials

Abstract representation of credentials.
"""
abstract type AbstractCredentials end

"""
    struct Credentials

Represents account information for a Reddit script application.
"""
struct Credentials <: AbstractCredentials
    id::AbstractString
    secret::AbstractString
    useragent::AbstractString
    username::AbstractString
    password::AbstractString
end

"""
    struct Authorized

Represents a credentials that has been authorized and recieved an
access token.
"""
struct Authorized <: AbstractCredentials
    id::AbstractString
    secret::AbstractString
    useragent::AbstractString
    username::AbstractString
    password::AbstractString
    token::AbstractString
end

"""
    struct Session

Represents an authorized session with Reddit's API.
"""
mutable struct Session
    creds::Union{Authorized, Array{Authorized}}
    Session() = new()
end

"""
    authorize(c::Credentials)

Use Credentials to request an acess token and return Authorized.
"""
function authorize(c::Credentials)
    Authorized(c.id, c.secret, c.useragent, c.username, c.password, token(c))
end

"""
    authorize(id::AbstractString,
              secret::AbstractString,
              useragent::AbstractString,
              username::AbstractString,
              password::AbstractString)

Use reddit application account information to request an acess token and
return Authorized.
"""
function authorize(id::AbstractString,
                   secret::AbstractString,
                   useragent::AbstractString,
                   username::AbstractString,
                   password::AbstractString)
    Authorized(id, secret, useragent, username, password,
               token(id, secret, useragent, username, password))
end

"""
    readconfig(name::AbstractString)

Create new Credentials from the specified client in the default ini file.
"""
function readconfig(name::AbstractString)
    readconfig(name, CONFIG)
end

"""
    readconfig(name::AbstractString, config::AbstractString)

Create new Credentials from the specified client in the specified ini file.
"""
function readconfig(name::AbstractString, config::AbstractString)
    conf = ConfParse(config)
    parse_conf!(conf)
    id = retrieve(conf, name, "client_id")
    secret = retrieve(conf, name, "client_secret")
    useragent = retrieve(conf, name, "user_agent")
    username = retrieve(conf, name, "username")
    password = retrieve(conf, name, "password")
    Credentials(id, secret, useragent, username, password)
end

"""
    default()

Get the default credentials.
"""
function default()
    if isdefined(DEFAULT_SESSION, :creds)
        DEFAULT_SESSION.creds
    else
        println(
        """
        Error: Default credentials not set.
        Use default!(creds::Union{Authorized, Array{Authorized}})
        to set default credentials.
        """)
    end
end

"""
    default!(creds::Union{Authorized, Array{Authorized}})

Set the default credentials.
"""
function default!(creds::Union{Authorized, Array{Authorized}})
    DEFAULT_SESSION.creds = creds
end

"""
    token(c::Credentials)

Get token with Credentials.
"""
function token(c::Credentials)
    token(c.id, c.secret, c.username, c.password)
end

"""
token(id::AbstractString,
      secret::AbstractString,
      username::AbstractString,
      password::AbstractString)

Get token with client_id, client_secret, username, and password.
"""
function token(id::AbstractString,
               secret::AbstractString,
               username::AbstractString,
               password::AbstractString)
    auth = encode("$id:$secret")
    resp = HTTP.request("POST", REDDIT_URL*"/api/v1/access_token",
        ["Authorization" => "Basic $auth"],
        "grant_type=password&username=$username&password=$password")
    body = String(resp.body)
    JSON.parse(body)["access_token"]
end


text:
using Random

""" Energy of the Dicke model """
function Energy(x, parameters)
    P, p, Q, q = x
    λ, δ, ω, ω₀  = parameters

    return 0.5 * ω * (p*p + q*q) + 0.5 * ω₀ * (P*P + Q*Q) + λ * sqrt(1.0 - 0.25 * (P*P + Q*Q)) * ((1.0 + δ) * P * q - (1.0 - δ) * p * Q)
end

""" Filling initial conditions (at a given energy calculates randomly all missing points) """
function InitialCondition!(x0, e, parameters)
    numMissingValues = length(x0) - length(collect(skipmissing(x0)))

    P, p, Q, q = x0
    λ, δ, ω, ω₀  = parameters

    if numMissingValues == 1 && ismissing(q) && p == 0.0
        s2 = 1.0 - 0.25 * (P*P + Q*Q)
        if s2 < 0
            return false
        end

        s = sqrt(s2)
        b = (1.0 + δ) * λ * P * s

        discriminant = b*b + 2.0 * ω * (e - 0.5 * ω₀ * (P*P + Q*Q))
        if discriminant < 0
            return false
        end

        sign = rand() >= 0.5 ? 1.0 : -1.0

        q = (sign * sqrt(discriminant) - b) / ω
        x0[4] = q

        return true
    end

    if numMissingValues == 1 && ismissing(Q) && P == 0.0
        a = λ * (1.0 - δ) * p
        b = 0.5 * ω * (p * p + q * q)

        discriminant1 = a * a - (e - b) * (e - 2.0 * ω₀ - b)

        if discriminant1 < 0
            return false
        end

        c = e * ω₀ + a * a - ω₀ * b
        d = ω₀ * ω₀ + a * a

        discriminant2p = 2.0 * (c + a * sqrt(discriminant1)) / d
        discriminant2m = 2.0 * (c - a * sqrt(discriminant1)) / d

        if discriminant2p < 0 && discriminant2m < 0
            return false
        end

        # Not all of the found initial conditions have the desired energy; let's select the correct ones
        ics = []

        x0[3] = sqrt(discriminant2m)
        if isapprox(e, Energy(x0, parameters))
            append!(ics, x0[3])
        end

        x0[3] = -sqrt(discriminant2m)
        if isapprox(e, Energy(x0, parameters))
            append!(ics, x0[3])
        end

        x0[3] = sqrt(discriminant2p)
        if isapprox(e, Energy(x0, parameters))
            append!(ics, x0[3])
        end

        x0[3] = -sqrt(discriminant2p)
        if isapprox(e, Energy(x0, parameters))
            append!(ics, x0[3])
        end

        if length(ics) == 0
            x0[3] = missing
            return false
        end

        x0[3] = rand(ics)
        return true
    end

    return false
end

""" True if given trajectory is within the kinematically accessible domain """
function CheckDomain(x, parameters, t)
    P, p, Q, q = x[1:4]
    s2 = 1.0 - 0.25 * (P*P + Q*Q)
    return s2 < 0
end

function EquationOfMotion!(dx, x, parameters, t)
    P, p, Q, q = x[1:4]             # Equations of motion for trajectories

    s2 = 1.0 - 0.25 * (P*P + Q*Q)
    if s2 < 0
        s2 = 0                      # With isoutofdomain=CheckDomain the calculation shouldn't enter here, but in enters anyway
        @error("s negative!")
    end

    Φ = reshape(x[5:end], 4, 4)     # Tangent dynamics
    λ, δ, ω, ω₀ = parameters.modelParameters

    s = sqrt(s2)
    b = (1.0 + δ) * P * q - (1.0 - δ) * p * Q

    ls_4 = 0.25 * λ * s
    l_4s = 0.25 * λ/ s

    Q2_4s2 = 0.25 * Q * Q / s2
    P2_4s2 = 0.25 * P * P / s2

    a = l_4s * P * Q

    HPQ = -l_4s * (0.25 * b * P * Q / s2 + (1.0 + δ) * Q * q - (1.0 - δ) * P * p)
    HpQ = λ * s * (1.0 - δ) * (Q2_4s2 - 1.0)
    HQQ = ω₀ - l_4s * b * (Q2_4s2 + 1.0) + 2.0 * l_4s * (1.0 - δ) * p * Q
    HQq = -a * (1.0 + δ)

    HPq = -λ * s * (1.0 + δ) * (P2_4s2 - 1.0)
    Hpq = 0.0
    Hqq = ω

    HPP = ω₀ - l_4s * b * (P2_4s2 + 1.0) - 2.0 * l_4s * (1.0 + δ) * P * q
    HPp = a * (1.0 - δ)

    Hpp = ω

    dx[1] = -ω₀ * Q + λ * (0.25 * b * Q / s + (1.0 - δ) * p * s)
    dx[2] = -ω * q - λ * (1.0 + δ) * P * s
    dx[3] = ω₀ * P + λ * (-0.25 * b * P / s + (1.0 + δ) * q * s)
    dx[4] = ω * p - λ * (1.0 - δ) * Q * s

    j = [-HPQ -HpQ -HQQ -HQq;
         -HPq -Hpq -HQq -Hqq;
         HPP HPp HPQ HPq;
         HPp Hpp HpQ Hpq]

    dx[5:end] = j * Φ
end

max_stars_repo_name: PavelStransky/ClassicalDynamics

max_issues_repo_name: PavelStransky/ClassicalDynamics

max_forks_repo_name: PavelStransky/ClassicalDynamics


text:
using .Variational

export @hilbertspace, @varform, @discretise
export solve

type DiscreteEquation
  equation
  trial_space_dict # dictionary mapping indices into trial space to FE spaces
  test_space_dict  # dictionary mapping indices into test space to FE spaces
end

function discretise(eq, space_mappings::Pair...)
    trial_space_dict = Dict()
    test_space_dict = Dict()
    for sm in space_mappings

        found = false
        sm.first.space == eq.lhs.trial_space && (dict = trial_space_dict; found = true)
        sm.first.space == eq.lhs.test_space  && (dict = test_space_dict;  found = true)
        @assert found "Vector $(sm.first) neither in test nor in trial space"

        @assert !haskey(dict, sm.first.idx) "multiple mappings for $(sm.first)"
        dict[sm.first.idx] = sm.second
    end

    # check that all symbols where mapped
    for p in 1:length(eq.lhs.trial_space) @assert haskey(trial_space_dict,p) end
    for p in 1:length(eq.lhs.test_space)  @assert haskey(test_space_dict, p) end

    DiscreteEquation(eq, trial_space_dict, test_space_dict)
end

# export discretisation, @discretisation
# function discretisation(space_mappings::Pair...)
#     space_dict = Dict()
#     for sm in space_mappings
#         space_dict[sm.first.idx] = sm.second
#     end
#     space_dict
# end
#
# macro discretisation(pairs...)
#     r = :(discretisation())
#     for p in pairs
#         x = p.args[2]
#         X = p.args[3]
#         push!(r.args, :($x=>$X))
#     end
#     return esc(r)
# end

"""
    discr(eq, pairs...)

This macro provides syntactical sugar for the definition of a discretisation
of a varational formulation. Given a variational equation EQ: Find j ∈ X such
that for all k ∈ Y a(k,j) = f(k) can be discretised by stating:

    eq = @discretise EQ j∈X k∈Y
"""
macro discretise(eq, pairs...)
    r = :(BEAST.discretise($eq))
    for p in pairs
        x = p.args[2]
        X = p.args[3]
        push!(r.args, :($x=>$X))
    end
    return esc(r)
end

export sysmatrix, rhs

sysmatrix(eq::DiscreteEquation) = assemble(eq.equation.lhs, eq.test_space_dict, eq.trial_space_dict)
rhs(eq::DiscreteEquation) = assemble(eq.equation.rhs, eq.test_space_dict)

function assemble(lform::LinForm, test_space_dict)

    terms = lform.terms
    T = Complex128

    I = Int[1]
    for p in 1:length(lform.test_space)
        X = test_space_dict[p]
        push!(I, last(I) + numfunctions(X))
    end
    B = zeros(T, last(I)-1)

    for t in terms

        α = t.coeff
        a = t.functional
        m = t.test_id
        X = test_space_dict[m]
        o = t.test_ops

        # act with the various ops on X
        for op in reverse(o)
            Y = X;
            X = op[1](Y, op[2:end]...)
        end

        b = assemble(a, X)
        B[I[m] : I[m+1]-1] = b
    end

    return B
end

function assemble(bilform::BilForm, test_space_dict, trial_space_dict)

  lhterms = bilform.terms
  T = Complex128 # TDOD: Fix this

  # determine the offsets of the different blocks in the sys matrix
  I = Int[1]
  J = Int[1]

  for p in 1:length(bilform.test_space)
    X = test_space_dict[p]
    push!(I, last(I) + numfunctions(X))
  end

  for q in 1:length(bilform.trial_space)
    Y = trial_space_dict[q]
    push!(J, last(J) + numfunctions(Y))
  end

  # allocate the memory for the matrices
  Z = zeros(T, last(I)-1, last(J)-1)

  # For each block, compute the interaction matrix
  for t in lhterms

      α = t.coeff
      a = t.kernel

      m = t.test_id
      x = test_space_dict[m]
      for op in reverse(t.test_ops)
          x = op[1](x, op[2:end]...)
      end

      n = t.trial_id
      y = trial_space_dict[n]
      for op in reverse(t.trial_ops)
          y = op[1](y, op[2:end]...)
      end

      r = I[m] : (I[m+1] - 1)
      c = J[n] : (J[n+1] - 1)

      z = assemble(a, x, y)
      Z[r,c] += α * z
  end

  return Z
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

module GeoClustering

using Meshes
using Tables
using TableOperations
using Distances
using Statistics
using LinearAlgebra

import Meshes: partition

include("slic.jl")

export
  SLIC,
  partition

end


text:
# testing nearal ODE
# https://discourse.julialang.org/t/diffeqflux-with-time-as-additional-input-to-neural-ode/26456/3
# https://discourse.julialang.org/t/defining-and-input-signal-for-a-differential-equation/23703/10
# https://github.com/JuliaDiffEq/DiffEqFlux.jl#universal-differential-equations-for-neural-optimal-control
# https://docs.juliadiffeq.org/stable/tutorials/ode_example/#Defining-Parameterized-Functions-1

using Pkg
Pkg.activate(@__DIR__)
pkg"instantiate"
0
pkg"precompile"

using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--individuals"
            help = "Number of simultaneous training individuals to compute"
            arg_type = Int
            default = 10
        "--maxiters"
            help = "Max ADAM optimization iterations"
            arg_type = Int
            default = 100
        "--stop"
            help = "number of encoder readings to consume for trahectory length"
            arg_type = Int
            default = 1300
        # "--flag1"
        #     help = "an option without argument, i.e. a flag"
        #     action = :store_true
        # "arg1"
        #     help = "a positional argument"
        #     required = true
    end

    return parse_args(s)
end

pargs = parse_commandline()

# using Revise

# using Flux
# using DiffEqFlux, OrdinaryDiffEq, Optim

using Dates, DelimitedFiles
using DiffEqFlux, Flux, Optim, OrdinaryDiffEq
using DifferentialEquations

# Requires Julia 1.3 and up
import Base.Threads.@spawn

## Fix issue Zygote.jl #440 #516
using Zygote

Zygote.@adjoint function vcat(xs::Union{Number, AbstractVector}...)
    vcat(xs...), dy -> begin
        d = 0
        map(xs) do x
            x isa Number && return dy[d+=1]
            l = length(x)
            dx = dy[d+1:d+l]
            d += l
            return dx
        end
    end
end

## Load MAT VicPrk data

using TransformUtils
using Plots
using MAT
using Interpolations

vpdata = matread(ENV["HOME"]*"/data/dataset_victoria_park/original_MATLAB_dataset/aa3_dr.mat")

# utility functions

include(joinpath(@__DIR__, "..", "wheeled", "parametricWheelOdo.jl"))

vpdata["time"] .*= 1e-3
vpdata["time"] = round.(vpdata["time"][:], digits=5)

vpsensors = hcat(vpdata["time"],
                 vpdata["speed"],
                 vpdata["steering"])
#

# see the input data
# plot(vpsensors[:,2])
# plot(vpsensors[:,3])

# # Likely unused
# steer =  LinearInterpolation( (vpdata["time"],) , vpdata["steering"][:])
# speed =  LinearInterpolation( (vpdata["time"],) , vpdata["speed"][:])

## Get baseline solution over all data
odo = allOdoEasy(vpsensors)
# vcat(ode_data, zeros(Float32, size()))

# calculate local speed
deltadist = sum( (odo[1:end-1,1:2] - odo[2:end,1:2]).^2, dims=2) .|> sqrt |> vec
push!(deltadist, deltadist[end])
bodySpeed = deltadist*40

# Callback to show current cost during training
cb = function (θ,l)
  println("thread=$(Threads.threadid()), cost=$l")
  return false
end

# Train a single instancwe
function doTraining(model::FastChain,
                    u0::Vector{<:Real};
                    p::Vector{<:Real}=initial_params(model),
                    maxiters=50,
                    α=0.05/(10*rand()) )
  #
  θlocal = Float32[u0;p]

  function dudt_(u,p,t)
    tidx = round(Int, (t-tspan[1])*40+1) # interpolation approach has issues with train
    tidx = 1 < tidx ? tidx : 1
    sp = vpdata["speed"][tidx, 1] + 1e-6*randn()
    st = vpdata["steering"][tidx, 1]
    input = [u; sp; st]
    model(input, p)
  end

  prob = ODEProblem(dudt_, u0, tspan, p)

  function predict_adjoint(θ)
    Array(  concrete_solve(prob, Tsit5(), θ[1:5], θ[6:end], saveat=ts)  )
  end

  function loss_adjoint(θ)
    pred = predict_adjoint(θ)
    loss = 0.0
    for idx in 1:size(ode_data,2)
      loss = 10*(ode_data[1,idx]-pred[1,idx])^2 +
             10*(ode_data[2,idx]-pred[2,idx])^2 +
             5*(cos(ode_data[3,idx])-pred[3,idx])^2 +
             5*(sin(ode_data[3,idx])-pred[4,idx])^2 +
             2*(bodySpeed[idx]-pred[5,idx])^2 +
             (1-pred[3,idx]^2-pred[4,idx]^2)^2
    end
    loss
  end

  # l = loss_adjoint(θlocal)
  # cb(θlocal,l)

  res = DiffEqFlux.sciml_train(loss_adjoint, θlocal, ADAM(α), cb = cb, maxiters=maxiters)
  # res = DiffEqFlux.sciml_train(loss_adjoint, θ, BFGS(initial_stepnorm=0.01), cb = cb)

  return res
end

function trainGeneration(models::Vector{FastChain}, u0; maxiters::Int=50, MC::Int=10)

  TASKS = Vector{Any}(undef,MC)
  for i in 1:MC
    TASKS[i] = Threads.@spawn doTraining(models[i], u0, maxiters=maxiters)
  end

  RES = Vector{Any}(undef,MC)
  PP = Vector{Vector{Float32}}(undef, MC)
  CO = zeros(MC)
  for i in 1:MC
    RES[i] = fetch(TASKS[i])
    CO[i] = RES[i].minimum
    PP[i] = deepcopy(RES[i].minimizer)
    @show i, RES[i].minimum
  end

  # TASKS = Vector{Any}(undef,MC)
  # RES = Vector{Any}(undef,MC)
  # GC.gc()
  ord = sortperm(CO)
  return PP[ord], CO[ord]
end

function main(vpdata, ode_data, usrmodel::FastChain; MC::Int=10, maxiters::Int=50)
  global ts, tspan

  t0 = tspan[1]
  # t0 = vpdata["time"][1] |> Float32

  # x_i, y_i, cosθ_i, sinθ_i, vx_i
  u0 = Float32[0.; 0.; 1.; 0.; 0.]
  # u0 = zeros(Float32, 4) #Float32(1.1)

  ## Specialized Neural model with time dependent inputs==========================
  models = Vector{FastChain}(undef, MC)
  for thr in 1:MC
    models[thr] = deepcopy(usrmodel)
  end

  PPs, COs = trainGeneration(models, u0; maxiters=maxiters, MC=MC)

  return PPs, COs
end

START=1
STOP=pargs["stop"]
ode_data = odo[START:STOP,:]' .|> Float32
datasize = size(ode_data,2)
tspan = (Float32(vpdata["time"][1]), Float32(vpdata["time"][datasize]))
ts = range(tspan[1],tspan[2],length=datasize)

# 7 inputs are x, y, cθ, sθ, speed, speed, steer
usrmdl = FastChain(FastDense(7,20,tanh),
                   FastDense(20,20,tanh),
                   FastDense(20,5,tanh))
#

# plot(odo[START:STOP,1], odo[START:STOP,2])
# plot(odo[:,1], odo[:,2])
# plot(bodySpeed[START:STOP])

PPs, COs = main(vpdata, ode_data, usrmdl, MC=pargs["individuals"], maxiters=pargs["maxiters"])

# GC.gc()

##==============================================================================
## Plot Store result

dtim = now()
resultsFolder = ENV["HOME"]*"/Documents/networks/$dtim/"
mkpath(resultsFolder)

fid = open(resultsFolder*"model.txt", "w")
for i in 1:length(usrmdl.layers)
  ll = usrmdl.layers[i]
  println(fid, "layer $i")
  println(fid, "in: $(ll.in)")
  println(fid, "out: $(ll.out)")
  println(fid, "act: $(ll.σ)")
end
close(fid)

PPm = hcat(PPs...)
# COs

writedlm(resultsFolder*"loss.txt", COs)
writedlm(resultsFolder*"params.txt", PPm')

heatmap(PPm)
savefig(resultsFolder*"heatmap.svg")

function dudt_test(u,p,t)
  global usrmdl, tspan, vpdata
  tidx = round(Int, (t-tspan[1])*40+1) # interpolation approach has issues with train
  tidx = 1 < tidx ? tidx : 1
  sp = vpdata["speed"][tidx, 1] + 1e-6*randn()
  st = vpdata["steering"][tidx, 1]
  input = [u; sp; st]
  usrmdl(input, p)
end

u0_test = Float32[ode_data[:,1]; 0.]
prob_test = ODEProblem(dudt_test, u0_test, tspan) #, p

function predict_adjoint_test(θ)
  Array(  concrete_solve(prob_test, Tsit5(), θ[1:5], θ[6:end], saveat=ts)  )
end

for i in 1:length(COs)

XX = predict_adjoint_test(PPs[i])

DRx = [odo[START:STOP,1]';XX[1,:]']'
DRy = [odo[START:STOP,2]';XX[2,:]']'
DRv = [bodySpeed[START:STOP]';XX[5,:]']'

DRct = [cos.(odo[START:STOP,3]');XX[3,:]']'
DRst = [sin.(odo[START:STOP,3]');XX[4,:]']'

fn = plot(DRx, DRy, fmt=:svg)

savefig(resultsFolder*"result_xy_$i.svg")

fn = plot(DRv, fmt=:svg)
savefig(resultsFolder*"result_v_$i.svg")

fn = plot(DRct, fmt=:svg)
savefig(resultsFolder*"result_costh_$i.svg")
fn = plot(DRst, fmt=:svg)
savefig(resultsFolder*"result_sinth_$i.svg")

# plot(XX[1:2,:]')
# plot(XX[1,:], XX[2,:])

end

#


text:
# # [Testing LiterateTest.jl](@id tests)

using Test
using LiterateTest

# This file tests LiterateTest.jl using plain Literate.jl-compatible
# source code.

# ## Helper function

decode_output(output) = join((strip(ln, '|') for ln in split(output, "\n")), "\n")

@assert decode_output("""
    |a|
    |b|
    """) === """
    a
    b
    """

# ## Removing `ans = begin` and `end`

input = """
    ans = begin
        1 + 2
    end
    """

output = decode_output("""
    |ans = begin # hide|
    |1 + 2|
    |end # hide|
    """)

@assert LiterateTest.preprocess(input) == output

# ## Removing `@testset`

input = """
    @testset begin
        @test 1 + 2 == 3
    end
    """

output = """
    """

@assert LiterateTest.preprocess(input) == output

# Above two patterns can be combined to write tests like this that are
# not visible in the Documenter.jl output:

ans = begin
    1 + 2  # only this line is visible in the documentation
end
@testset begin
    @test ans == 3
end
nothing  # hide

# ## Removing `@test begin` and `end ...`

input = """
    @test begin
        1 + 2
    end == 3
    """

output = """
    1 + 2
    """

@assert LiterateTest.preprocess(input) == output

# When extracting the result of `@test`, it's useful to use `global`.
# However, `global` is not needed when `@test begin ... end` is
# removed.  Thus, `global` is removed from the code as well:

input = """
    @test begin
        global y = 1 + 2
    end == 3
    y + 3
    """

output = """
    y = 1 + 2
    y + 3
    """

@assert LiterateTest.preprocess(input) == output

# `@test` itself can be indented

input = """
    if VERSION >= v"1.4"
        @test begin
            code
        end == 1
    end
    """

output = """
    if VERSION >= v"1.4"
    code
    end
    """

@assert LiterateTest.preprocess(input) == output

# ## Extracting `$code` from `@evaltest "$code"`

input = """
    @evaltest "1 + 2" begin
        @test ans == 3
        @test ans isa Int
    end
    """

output = """
    1 + 2
    """

@assert LiterateTest.preprocess(input) == output

# Note that `"$code"` will be used as the `@testset` title when it's
# evaluated withtout the transformation.

@evaltest "1 + 2" begin
    @test ans == 3
    @test ans isa Int
end
nothing  # hide

# !!! note
#
#     `@evaltest "1 + 2" begin ...` above is shown in this
#     Documenter.jl output because this file is *not* processed by
#     [`LiterateTest.preprocess`](@ref).

#  When using `$` in the code, `raw""` can be useful:

input = raw"""
    x = 1
    @evaltest raw":(1 + $x)" begin
        @test ans == :(1 + $x)
    end
    """

output = raw"""
    x = 1
    :(1 + $x)
    """

@assert LiterateTest.preprocess(input) == output

# Demo:

x = 1
@evaltest raw":(1 + $x)" begin
    @test ans == :(1 + $x)
end
nothing  # hide

# ## Showing error thrown from `$code` in `@evaltest_throw "$code"`

input = """
    @evaltest_throw \"\"\"error("msg")\"\"\" begin
        @test ans == ErrorException(1)
    end
    """

Text(LiterateTest.preprocess(input))

# It's a bit tricky to test this in Literate.jl:

output = decode_output("""
    |err = try # hide|
    |error("msg")|
    |catch _err; _err; end # hide|
    |print(stdout, "ERROR: ") # hide|
    |showerror(stdout, err) # hide|
    |#-|
    """)

@assert LiterateTest.preprocess(input) == output

# Demo:

@evaltest_throw """error("msg")""" begin
    @test ans == ErrorException("msg")
end
nothing  # hide

# ## Testing an exception from multi-line code with `@testset_error`

input = """
    @testset_error try
        if true
            error(1)
        end
    catch err
        @test err == ErrorException(1)
    end
    """

Text(LiterateTest.preprocess(input))

#-

output = decode_output("""
    |err = try # hide|
    |if true|
    |    error(1)|
    |end|
    |catch _err; _err; end # hide|
    |print(stdout, "ERROR: ") # hide|
    |showerror(stdout, err) # hide|
    |#-|
    """)

Text(output)

#-

@assert LiterateTest.preprocess(input) == output

# Demo:

@testset_error "@testset_error" try
    if true
        error(1)
    end
catch err
    @test err == ErrorException("1")
end
nothing  # hide

# ## De-indenting `code` in `@dedent begin code end`

input = """
    @dedent @time begin
        code
    end
    """

output = decode_output("""
    |@time begin # hide|
    |code|
    |end # hide|
    """)

@assert LiterateTest.preprocess(input) == output

# Demo:

@dedent @time begin
    1 + 1
end

# Expressions after `end` are ignored

input = """
    @dedent @assert begin
        code
    end == 1
    """

output = decode_output("""
    |@assert begin # hide|
    |code|
    |end == 1 # hide|
    """)

@assert LiterateTest.preprocess(input) == output

# Demo:

@dedent @assert begin
    1 + 1
end == 2

# `@dedent` works with assignments

input = """
    @dedent y1 = begin
        1 + 1
    end
    @dedent y2 = begin
        3 - 1
    end
    @test y1 == y2
    """

output = decode_output("""
    |y1 = begin # hide|
    |1 + 1|
    |end # hide|
    |y2 = begin # hide|
    |3 - 1|
    |end # hide|
    """)

@assert LiterateTest.preprocess(input) == output

# Demo:

@dedent y1 = begin
    1 + 1
end
@dedent y2 = begin
    3 - 1
end
@test y1 == y2
nothing  # hide

# `@dedent` itself can be indented

input = """
    if VERSION >= v"1.4"
        @dedent @assert begin
            code
        end == 1
    end
    """

output = decode_output("""
    |if VERSION >= v"1.4"|
    |@assert begin # hide|
    |code|
    |end == 1 # hide|
    |end|
    """)

@assert LiterateTest.preprocess(input) == output

# ## Combine `@testset_error`, `@dedent` and `@eval` to test macro expansion

input = """
    @testset_error try
        @dedent @eval begin
            LiterateTest.@this_macro_does_not_exist
        end
    catch err
        msg = sprint(showerror, err)
        @test occursin("UndefVarError: @this_macro_does_not_exist not defined", msg)
    end
    """

Text(LiterateTest.preprocess(input))

# Note that all lines except `LiterateTest.@this_macro_does_not_exist` ends with `# hide`.

output = decode_output("""
    |err = try # hide|
    |@eval begin # hide|
    |LiterateTest.@this_macro_does_not_exist|
    |end # hide|
    |catch _err; _err; end # hide|
    |print(stdout, "ERROR: ") # hide|
    |showerror(stdout, err) # hide|
    |#-|
    """)

@assert LiterateTest.preprocess(input) == output

# Demo:

@testset_error "@testset_error + @dedent + @eval" try
    @dedent @eval begin
        LiterateTest.@this_macro_does_not_exist
    end
catch err
    msg = sprint(showerror, err)
    @test occursin("UndefVarError: @this_macro_does_not_exist not defined", msg)
end
nothing  # hide

# ## Removing JuliaFormatter.jl on/off toggles

# [The on/off comments for JuliaFormatter.jl](
# https://domluna.github.io/JuliaFormatter.jl/dev/skipping_formatting/
# ) are removed by `LiterateTest.preprocess`:

input = """
    a
    #! format: off
    b
    #! format: on
    c
    """

output = """
    a
    b
    c
    """

@assert LiterateTest.preprocess(input) == output

# There is also [`LiterateTest.preprocess_juliaformatter`](@ref) that
# does only this preprocessing:

@assert LiterateTest.preprocess_juliaformatter(input) == output


text:
module RcVolume
using Parameters
export ExpVol, ExpVolResults

struct ExpVol
end

struct ExpVolResults
end

end


text:
module HartreeFock

using BlockArrays
using JSON
using LinearAlgebra
using SpecialFunctions
using StaticArrays

include("constants.jl")
include("parser.jl")
include("utils.jl")
include("particles.jl")
include("basis.jl")
include("mole.jl")
include("Libcint/Libcint.jl")
include("cintor.jl")
include("integral.jl")
include("intor.jl")
include("diis.jl")
include("scf.jl")
include("electrostatics.jl")
include("population.jl")

using .Libcint

end # module


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0-beta3
#     language: julia
#     name: julia-1.7
# ---

# %%
using BenchmarkTools

f(n) = 1 + sum(abs2, 1:n)
g(n) = 1 + sum(abs2(k) for k in 1:n)

@show VERSION
@show f(10^6) == g(10^6)
print("f(10^6):")
@btime f(10^6)
print("g(10^6):")
@btime g(10^6)

# %%
using BenchmarkTools

f(n) = 1 + sum(abs2, 1:n)
g(n) = 1 + sum(abs2(k) for k in 1:n)

@show VERSION
@show f(10^6) == g(10^6)
print("f(10^6):")
@btime f(10^6)
print("g(10^6):")
@btime g(10^6)

# %%
import Base.Broadcast.preprocess_args
import Base.Broadcast.preprocess
@inline preprocess_args(dest, args::Tuple) = (Base.Broadcast.preprocess(dest, args[1]), Base.Broadcast.preprocess_args(dest, Base.tail(args))...)
@inline preprocess_args(dest, args::Tuple{Any}) = (Base.Broadcast.preprocess(dest, args[1]),)
@inline preprocess_args(dest, args::Tuple{}) = ()

using BenchmarkTools

f(n) = 1 + sum(abs2, 1:n)
g(n) = 1 + sum(abs2(k) for k in 1:n)

@show VERSION
@show f(10^6) == g(10^6)
print("f(10^6):")
@btime f(10^6)
print("g(10^6):")
@btime g(10^6)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using Flux
using Statistics
using Flux: onehotbatch, onecold, crossentropy, throttle
using Base.Iterators: repeated

path1=realpath(dirname(@__FILE__)*"/../..")
include(string(path1,"/src/ANNBN.jl"))
using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random, MLDatasets, Plots, Printf, Debugger
train_x, yy_train_all = MLDatasets.MNIST.traindata()
test_x,  yy_test_all  = MLDatasets.MNIST.testdata()
xx_train=MLDatasets.MNIST.convert2features(MLDatasets.MNIST.traintensor())'
xx_test=MLDatasets.MNIST.convert2features(test_x)'
xx_train=convert(Array{Float64,2},xx_train)
xx_test=convert(Array{Float64,2},xx_test)
yy_train_all=convert(Array{Float64,1},yy_train_all)
yy_test_all=convert(Array{Float64,1},yy_test_all)
i_train=size(xx_train,1)
i_test=size(xx_test,1)
vars=size(xx_train,2)

tol1=0.02
num_to_test=0
yy_train=copy(yy_train_all)
i1=abs.(yy_train.-num_to_test).<0.01;i2=abs.(yy_train.-num_to_test).>0.01;yy_train[i1].=1;;yy_train[i2].=0;
mi=minimum(yy_train);yy_train.-=mi;ma=maximum(yy_train);yy_train./=ma;yy_train.*=1.0-2tol1;yy_train.+=tol1

yy_test=copy(yy_test_all)
i1=abs.(yy_test.-num_to_test).<0.01;i2=abs.(yy_test.-num_to_test).>0.01;yy_test[i1].=1;yy_test[i2].=0;
yy_test.-=mi;yy_test./=ma;yy_test.*=1.0-2tol1;yy_test.+=tol1

inds_all=1:i_train;neurons=1000;items_per_neuron=(Int64(floor(i_train/(neurons))))*ones(Int64,neurons);n_per_part=[0;cumsum(items_per_neuron)];n_per_part[end]=i_train
@time a_all,a_layer1,layer1=ANNBN.train_layer_1_sigmoid_fast(neurons,vars,i_train,n_per_part,inds_all,xx_train,yy_train)
predl1=[layer1 ones(i_train)]*a_layer1
maetr=mean(abs.(yy_train-predl1))
maximum(abs.(yy_train-predl1))
predl1_test,layer1_test=ANNBN.predict_new(a_all,a_layer1,xx_test,i_test,neurons);
maete=mean(abs.(yy_test-predl1_test))
maximum(abs.(yy_test-predl1_test))
i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).<0.01;maetr=100sum(i2)/length(i2)
i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<0.01;
maete=100sum(i2)/length(i2)

neurons_all=Vector{Int64}(undef,0)
push!(neurons_all,length(a_layer1)-1)
weights=Vector{Array{Float32}}(undef,0)
aa=reduce(hcat, a_all)'[:,2:end]
push!(weights,convert(Array{Float32,2},aa))
bb=reduce(hcat, a_all)'[:,1]
push!(weights,convert(Array{Float32,1},bb))
# aa=reshape(a_layer1[1:end-1],1,neurons)
# push!(weights,convert(Array{Float32,2},aa))
# bb=convert(Vector{Float32},[a_layer1[end]])
# push!(weights,bb)

# _todo_Eq_15
nof_layers=1
for i=2:nof_layers
    global neurons_all,layer1,layer1_test,weights,a_layer1
    neurons=neurons_all[end]
    inds_all=1:i_train;items_per_neuron=(Int64(floor(i_train/(neurons))))*ones(Int64,neurons);n_per_part=[0;cumsum(items_per_neuron)];n_per_part[end]=i_train
    @time a_all,a_layer1,layer1=ANNBN.train_layer_1_sigmoid_fast(neurons,neurons,i_train,n_per_part,inds_all,layer1,yy_train)
    predl1=[layer1 ones(i_train)]*a_layer1
    maetr=mean(abs.(yy_train-predl1))
    maximum(abs.(yy_train-predl1))
    predl1_test,layer1_test=ANNBN.predict_new(a_all,a_layer1,layer1_test,i_test,neurons);
    maete=mean(abs.(yy_test-predl1_test))
    maximum(abs.(yy_test-predl1_test))
    i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).<0.01;maetr=100sum(i2)/length(i2)
    i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<0.01;
    maete=100sum(i2)/length(i2)

    println(" >>>>>>>>>>>>>>>>>>>>    ",i," ",maetr," ",maete," ",length(a_layer1)-1)

    push!(neurons_all,length(a_layer1)-1)
    aa=reduce(hcat, a_all)'[:,2:end]
    push!(weights,convert(Array{Float32,2},aa))
    bb=reduce(hcat, a_all)'[:,1]
    push!(weights,convert(Array{Float32,1},bb))

end

aa=reshape(a_layer1[1:end-1],1,neurons_all[end])
push!(weights,convert(Array{Float32,2},aa))
bb=convert(Vector{Float32},[a_layer1[end]])
push!(weights,bb)

weights00=copy(weights)

X=convert(Array{Float32,2},xx_train')
Y=convert(Array{Float32,2},yy_train')
m = Chain(
  Dense(vars, neurons_all[1], ANNBN.sigm1),
  # Dense(neurons_all[1], neurons_all[2], ANNBN.sigm1),
  # Dense(neurons_all[2], neurons_all[3], ANNBN.sigm1),
  # Dense(neurons_all[3], neurons_all[4], ANNBN.sigm1),
  # Dense(neurons_all[4], neurons_all[5], ANNBN.sigm1),
  # Dense(neurons_all[5], neurons_all[6], ANNBN.sigm1),
  # Dense(neurons_all[6], neurons_all[7], ANNBN.sigm1),
  # Dense(neurons_all[7], neurons_all[8], ANNBN.sigm1),
  # Dense(neurons_all[8], neurons_all[9], ANNBN.sigm1),
  # Dense(neurons_all[9], neurons_all[10], ANNBN.sigm1),
  Dense(neurons_all[1], 1))
loss(x, y) = mean((m(x).-y).^2)
# opt = ADAM()
opt = Descent(1e-2)
dataset = repeated((X,Y),10)
# evalcb = () -> @show(loss(X,Y),loss(Xte,Yte))
evalcb = () -> @show(100mean(abs.(round.(Tracker.data(m(X)[1,:])).-round.(yy_train)).<0.1), 100mean(abs.(round.(Tracker.data(m(Xte)[1,:])).-round.(yy_test)).<0.1))
# weights = Tracker.data.(params(m))

Flux.loadparams!(m, weights00)

loss(X, Y)
predl1=Tracker.data(m(X)[1,:])
weights = Tracker.data.(params(m))
maetr=mean(abs.(yy_train-predl1))
maximum(abs.(yy_train-predl1))
Xte=convert(Array{Float32,2},xx_test')
Yte=convert(Array{Float32,2},yy_test')
loss(Xte, Yte)
predl1_test=Tracker.data(m(Xte)[1,:])
maete=mean(abs.(yy_test-predl1_test))
maximum(abs.(yy_test-predl1_test))
i1=predl1.<=0.5;predl1[i1].=tol1;i1=predl1.>0.5;predl1[i1].=1.0-tol1;i2=abs.(predl1.-yy_train).<0.01;maetr=100sum(i2)/length(i2)
i1=predl1_test.<=0.5;predl1_test[i1].=tol1;i1=predl1_test.>0.5;predl1_test[i1].=1.0-tol1;i2=abs.(predl1_test.-yy_test).<0.01;
maete=100sum(i2)/length(i2)

100mean(abs.(round.(predl1).-round.(yy_train)).<0.1)
100mean(abs.(round.(predl1_test).-round.(yy_test)).<0.1)

println(">>>>>>>>>>>>>>>>>>>>>>>>>>START>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>")
@time for ep=1:1 Flux.train!(loss, params(m), dataset, opt, cb = throttle(evalcb, 1)); end


text:
# Displays

struct Inline end
struct Clipboard end

struct Editor end
struct Console end # deprecated
struct PlotPane end

# View data structures

mutable struct Model
  data
end

struct Tree
  head
  children::Vector{Any}
end

struct LazyTree
  head
  children::Function
end

struct SubTree
  label
  child
end

limit(s::AbstractString) = length(s) ≤ 5000 ? s : s[nextind(s, 0):prevind(s, 5000)]*"..."

struct Copyable
  view
  text::String
  Copyable(view, text::String) = new(view, limit(text))
end

Copyable(view, text) = Copyable(view, render(Clipboard(), text))
Copyable(view) = Copyable(view, view)

struct Link
  file::String
  line::Int
  contents::Vector{Any}
  function Link(file::AbstractString, line::Integer, contents...)
    isempty(contents) && (contents = ("$file:$line",))
    new(file, line, [contents...])
  end
end

Link(file::AbstractString, contents...) = Link(file, 0, contents...)

link(a...) = Link(a...)

mutable struct Row
  xs::Vector{Any}
  Row(xs...) = new(collect(xs))
end

mutable struct Table
  xs::Matrix{Any}
end

"""
    defaultrepr(x, lazy = false)

`render` fallback for types without any specialized `show` methods.

If `lazy` is true, then the type's fields will be loaded lazily when expanding the tree.
This is useful when the fields contain big elements that might need to be inspectable.

Can be used by packages to restore Juno's default printing if they have defined
a `show` method that should *not* be used by Juno:
```julia
Juno.render(i::Juno.Inline, x::myType) = Juno.render(i, Juno.defaultrepr(x))
```
"""
defaultrepr(x, lazy = false) = Main.Atom.defaultrepr(x, lazy)


text:
module Automa

import TranscodingStreams: TranscodingStream

mutable struct State{S <: TranscodingStream}
    # Stream
    stream::S
    # Machine state
    state::Int
    # Line number
    linenum::Int
    # Is record filled?
    filled::Bool
end

end # module

max_forks_repo_name: BioJulia/BioGenerics


text:
module Agents

using Distributed
using LightGraphs
using DataFrames
using Random
import Base.Iterators.product
import Base.iterate
import Base.length

include("core/model.jl")
include("core/space.jl")
include("core/agent_space_interaction.jl")
include("simulations/data_collector.jl")
include("simulations/step.jl")
include("simulations/parallel.jl")
include("simulations/sample.jl")
include("CA1D.jl")
include("CA2D.jl")

end # module


text:
julia_type_to_valtype(julia_type)::Ptr{wasm_valtype_t} =
    julia_type_to_valkind(julia_type) |> wasm_valtype_new

# TODO: the other value types
function WasmInt32(i::Int32)
    val = Ref(wasm_val_t(tuple((zero(UInt8) for _ = 1:16)...)))
    ptr = Base.unsafe_convert(Ptr{wasm_val_t}, Base.pointer_from_objref(val))
    ptr.kind = WASM_I32
    ptr.of.i32 = i

    val[]
end
function WasmInt64(i::Int64)
    val = Ref(wasm_val_t(tuple((zero(UInt8) for _ = 1:16)...)))
    ptr = Base.unsafe_convert(Ptr{wasm_val_t}, Base.pointer_from_objref(val))
    ptr.kind = WASM_I64
    ptr.of.i64 = i

    val[]
end
function WasmFloat32(f::Float32)
    val = Ref(wasm_val_t(tuple((zero(UInt8) for _ = 1:16)...)))
    ptr = Base.unsafe_convert(Ptr{wasm_val_t}, Base.pointer_from_objref(val))
    ptr.kind = WASM_F32
    ptr.of.f32 = f

    val[]
end
function WasmFloat64(f::Float64)
    val = Ref(wasm_val_t(tuple((zero(UInt8) for _ = 1:16)...)))
    ptr = Base.unsafe_convert(Ptr{wasm_val_t}, Base.pointer_from_objref(val))
    ptr.kind = WASM_F64
    ptr.of.f64 = f

    val[]
end

function julia_type_to_valkind(julia_type::Type)::wasm_valkind_enum
    if julia_type == Int32
        WASM_I32
    elseif julia_type == Int64
        WASM_I64
    elseif julia_type == Float32
        WASM_F32
    elseif julia_type == Float64
        WASM_F64
    else
        error("No corresponding valkind for type $julia_type")
    end
end

function valkind_to_julia_type(valkind::wasm_valkind_enum)
    if valkind == WASM_I32
        Int32
    elseif valkind == WASM_I64
        Int64
    elseif valkind == WASM_F32
        Float32
    elseif valkind == WASM_F64
        Float64
    else
        error("No corresponding type for kind $valkind")
    end
end

Base.convert(::Type{wasm_val_t}, i::Int32) = WasmInt32(i)
Base.convert(::Type{wasm_val_t}, i::Int64) = WasmInt64(i)
Base.convert(::Type{wasm_val_t}, f::Float32) = WasmFloat32(f)
Base.convert(::Type{wasm_val_t}, f::Float64) = WasmFloat64(f)

Base.convert(::Type{wasm_val_t}, val::wasm_val_t) = val
function Base.convert(julia_type, wasm_val::wasm_val_t)
    valkind = julia_type_to_valkind(julia_type)
    @assert valkind == wasm_val.kind "Cannot convert a value of kind $(wasm_val.kind) to corresponding kind $valkind"
    ctag = Ref(wasm_val.of)
    ptr = Base.unsafe_convert(Ptr{LibWasmtime.__JL_Ctag_4}, ctag)
    jl_val = GC.@preserve ctag unsafe_load(Ptr{julia_type}(ptr))
    jl_val
end

function Base.show(io::IO, wasm_val::wasm_val_t)
    name, maybe_val = if wasm_val.kind == WASM_I32
        "WasmInt32", wasm_val.of.i32 |> string
    elseif wasm_val.kind == WASM_I64
        "WasmInt64", wasm_val.of.i64 |> string
    elseif wasm_val.kind == WASM_F32
        "WasmFloat32", wasm_val.of.f32 |> string
    elseif wasm_val.kind == WASM_F64
        "WasmFloat64", wasm_val.of.f64 |> string
    else
        "WasmAny", ""
    end

    print(io, "$name($maybe_val)")
end


text:
@diffusion_process JansenRit begin
    :dimensions
    process --> 6
    wiener --> 3

    :parameters
    (A, a, B, b, C, νmax, v, r, μx, μy, μz, σx, σy, σz) --> Float64

    :conjugate
    phi(t, x) --> (
        (P.A*P.a*0.8P.C*sigm(P.C*x[1], P) - 2P.a*x[5] - P.a*P.a*x[2],),
        (P.A*P.a,),
        (0,),
        (0,),
        (0,),
        (0,),
        (0,),
        (0,),
        (0,),
        (0,),
        (0,),
        (0,),
        (0,),
    )

    :additional
    constdiff --> true
end

sigm(x, P::JansenRit) = P.νmax / (1 + exp(P.r*(P.v0 - x)))
μx(t, P::JansenRit) = P.μx
μy(t, P::JansenRit) = P.μy
μz(t, P::JansenRit) = P.μz

function b(t, x, P::JansenRit)
    ℝ{6}(
        x[4],
        x[5],
        x[6],
        P.δ*x[1]*x[2] - P.γ*x[2],
        P.A*P.a*(μx(t, P) + sigm(x[2] - x[3], P)) - 2P.a*x[4] - P.a*P.a*x[1],
        P.A*P.a*(μy(t, P) + 0.8P.C*sigm(P.C*x[1], P)) - 2P.a*x[5] - P.a*P.a*x[2],
        P.B*P.b*(μz(t, P) + 0.25P.C*sigm(0.25P.C*x[1], P)) - 2P.b*x[6] - P.b*P.b*x[3])
    )
end

function σ(t, x, P::JansenRit)
    @SMatrix [
        0.0  0.0  0.0;
        0.0  0.0  0.0;
        0.0 0.0  0.0;
        P.σx  0.0  0.0;
        0.0  P.σy 0.0;
        0.0  0.0  P.σz
    ]
end


text:
"""
    TrackingResult

This type stores a vector of `blobs` that were active when tracking terminated, a vector `dead` of blobs that died during tracking and a vector of all measurements obtained during tracking. This type supports `allblobs, trace`
"""
TrackingResult

Base.@kwdef mutable struct TrackingResult
    blobs::Vector{Blob} = Blob[]
    dead::Vector{Blob} = Blob[]
    coordinates::Vector{Trace} = Trace[]
end

Base.isempty(result::TrackingResult) = isempty(result.blobs) && isempty(result.dead)

"""
    allblobs(tr)

Get all blobs in a tracking result
"""
allblobs(tr::TrackingResult) = [tr.dead; tr.blobs]

"""
    trace(tr::TrackingResult; minlife=0)

Get all traces in a tracking result. Optionally filter based on minimum lifetime of a blob (the blob must have seen this many measurements).
"""
function trace(tr::TrackingResult; minlife=0)
    blobs = allblobs(tr)
    filter!(b->lifetime(b)>=minlife, blobs)
    trace.(blobs)
end

function tracem(tr::TrackingResult; minlife=0)
    blobs = allblobs(tr)
    filter!(b->lifetime(b)>=minlife, blobs)
    tracem.(blobs)
end

log_measurement(result::TrackingResult, m::Measurement) = push!(result.coordinates, m.coordinates)
log_measurement(result::TrackingResult, ms::Vector{Measurement}) = push!(result.coordinates, ms[1].coordinates)
Distributions.location(result::TrackingResult) = result.coordinates