text:
#import Base.isless
#import PyPlot
#using SparseArrays
#Allow sorting between Numbers and Strings
isless(A::Number,B::AbstractString) = false
isless(A::AbstractString,B::Number) = true

StringOrNumArray = Union{AbstractString,Array,Number}
UnionArray = Array{Union{AbstractString,Number}}
SparseOrTranspose = Union{AbstractSparseMatrix,Adjoint{<:Any,<:SparseMatrixCSC},Transpose{<:Any,<:SparseMatrixCSC}}

#Creation of Assoc require StrUnique to split Single-Character-separated String Sequence.
#include("stringarrayhelpers.jl")

#=
Type Assoc (Associative Array)
Support a
=#

struct Assoc
# TODO: "type" being depreciated, should change to struct or mutable struct
# Should be struct- operations on A return a new Assoc, not a changed A.
    row::UnionArray
    col::UnionArray
    val::UnionArray
    A::SparseOrTranspose

    # default uses min, should it be sum?
    Assoc(rowIn::StringOrNumArray,colIn::StringOrNumArray,valIn::StringOrNumArray) = Assoc(rowIn,colIn,valIn,min)
    Assoc(row::Array{Int64}, col::Array{Int64},val::Array{Int64},A::SparseOrTranspose) = new(row,col,val,A)
    #Setting Default Function

    #Assoc(rowIn::UnionArray,colIn::UnionArray,A::Adjoint{<:Any,<:SparseMatrixCSC}) = new(row,col,A)

    function Assoc(rowIn::Array{Union{AbstractString,Number}}, colIn::Array{Union{AbstractString,Number}}, valIn::Array{Union{AbstractString,Number}}, AIn::SparseOrTranspose)
        if (!isempty(valIn) && isassigned(valIn)) && isa(valIn[1],Number)
            valIn = convert(Array{Union{AbstractString,Number}},[1.0])
        end
        return new(rowIn,colIn,valIn,AIn)
    end

    function Assoc(rowIn::StringOrNumArray,colIn::StringOrNumArray,valIn::StringOrNumArray,funcIn::Function)
        if isempty(rowIn) || isempty(colIn) || isempty(valIn)
            # testing needed for isemtpy, for Matlab isemtpy is always possible TODO
            # Seems to work okay with String or NumArray type hard defined, Union type untested.
            # Should keep an eye.
            x = Array{Union{AbstractString,Number}}(undef)
            return Assoc(x,x,x,spzeros(0,0));
        end

        # Convert any scalar numbers to an array
        if isa(rowIn,Number)
            rowIn = Array{Union{AbstractString,Number},1}([rowIn])
        end
        if isa(colIn,Number)
            colIn = Array{Union{AbstractString,Number},1}([colIn])
        end
        if isa(valIn,Number)
            valIn = Array{Union{AbstractString,Number},1}([valIn])
        end

        # row, col, and val are always searchsortedfirst
        # i, j, and v are passed into the sparse constructor
        i = rowIn;
        j = colIn;
        v = valIn;
        row = rowIn;
        col = colIn;
        val = valIn;

        if isa(rowIn,AbstractString)
            row, i_out2in, i = StrUnique(rowIn);
        else
            # Get unique sorted keys
            row = unique(i)
            sort!(row)

            # Find index of row keys from triples in row (int indices for sparse matrix)
            i = convert(AbstractArray{Int64},[searchsortedfirst(row,x) for x in i])
        end

        if isa(colIn,AbstractString)
            col, j_out2in, j = StrUnique(colIn);
        else
            # Get unique sorted keys
            col = unique(j)
            sort!(col)

            # Find index of row keys from triples in row (int indices for sparse matrix)
            j = convert(AbstractArray{Int64},[searchsortedfirst(col,x) for x in j])
        end

        if isa(valIn,AbstractString)
            val, v_out2in, v = StrUnique(valIn);
        else
            val = unique(v)
            sort!(val)
            if (isa(valIn[1],AbstractString))

                # This bit ensures zeros are placed in any location where there are empty strings for value
                if val[1] == ""
                    val = val[2:end]
                    emptyidx = v .== ""
                else
                    emptyidx = []
                end
                v = convert(AbstractArray{Int64},[searchsortedfirst(val,x) for x in v])

                v[emptyidx] .= 0

            else # convert v to Int64 or Float64
                if any(isa.(v,Float64)) # If there are any floats, convert all to float
                    v = convert(AbstractArray{Float64},v)
                else
                    v = convert(AbstractArray{Int64},v)
                end
            end
        end

        # If any of r,c,v are length 1, expands to array with length of others
        NMax = maximum([length(i) length(j) length(v)]);
        if length(i) == 1
            i = convert(AbstractArray{Int64},repeat(i,NMax))
        end
        if length(j) == 1
            j = convert(AbstractArray{Int64},repeat(j,NMax))
        end
        if length(v) == 1
            v = convert(AbstractArray{Int64},repeat(v,NMax))
        end

        # Create the sparse matrix
        # If the values are string, assume that there are duplicates
        # and take the earliest one ( the numbers should be the same)
        if isa(val[1],AbstractString)
            A = sparse(i,j,v,length(row),length(col),min);
        else
            A = sparse(i,j,v,length(row),length(col),(+));
        end

        # If values are Numbers (not strings), actual values are stored in the sparse matrix
        # and "val" is set to [1.0] as a signal
        if isa(val[1],Number)
            val = convert(Array{Union{AbstractString,Number}},[1.0])
        end
        return new(row,col,val,A)
    end
end

#=
Returns an empty Associative Array
=#
function emptyAssoc()
    Assoc([],[],[])
end

#=
size: Return the dimensions of the Associative Array
=#
function size(A::Assoc)
    return size(A.A)
end

#=
nnz: Return the number of nonzeros in an Associative Array
=#
function nnz(A::Assoc)

    if isa(A.A,LinearAlgebra.Adjoint) || isa(A.A, LinearAlgebra.Transpose)
        return nnz(A.A.parent)
    else
        return nnz(A.A)
    end
end

#=
isempty : check if given Assoc is empty.
Note: Assoc can be considered empty even if there are mapping for potential or past values.
=#
function isempty(A::Assoc)
    return isempty(A.A)
end

#=
Adding related operations for Assoc_orig
=#
include("./Assoc/getindex.jl")
include("./Assoc/condense.jl")
include("./Assoc/operations.jl")
include("./Assoc/print.jl")
include("./Assoc/accessor.jl")
include("./Assoc/put.jl")
include("./Assoc/convert.jl")
include("./Assoc/broadcast.jl")
include("./Assoc/io.jl")
include("./Assoc/convertvals.jl")
include("./Assoc/bfs.jl")
#include("./Assoc_orig/find.jl")
#include("./Assoc_orig/no.jl")
#include("./Assoc_orig/isempty.jl")
#include("./Assoc_orig/logical.jl")
#include("./Assoc_orig/and.jl")
#include("./Assoc_orig/transpose.jl")
#include("./Assoc_orig/multiply.jl")
#include("./Assoc_orig/sqIn.jl")
#include("./Assoc_orig/sqOut.jl")
#include("./Assoc_orig/diag.jl")
#include("./Assoc_orig/plus.jl")
#include("./Assoc_orig/deepcondense.jl")
#include("./Assoc_orig/lt.jl")
#include("./Assoc_orig/gt.jl")
#include("./Assoc_orig/jld.jl")
#include("./Assoc_orig/equal.jl")
#include("./Assoc_orig/minus.jl")
#include("./Assoc_orig/emptyAssoc.jl")
#include("./Assoc_orig/size.jl")
#include("./Assoc_orig/printTriple.jl")
#include("./Assoc_orig/nnz.jl")
#include("./Assoc_orig/printFull.jl")
#include("./Assoc_orig/str2num.jl")

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
#overlap funcs {{{
## for binary states {{{
export overlap, overlapm

function overlap{T<:Real}(x::AbstractVector{T},y::AbstractVector{T})
    return dot(x, y)/length(x)
end

function overlap{T<:Real}(x::AbstractVector{Complex{T}},y::AbstractVector{Complex{T}})
    rsum = dot(real(x), real(y))
    isum = dot(imag(x), imag(y))
    return (rsum+isum)/length(x)*0.5
end

function overlap{T<:Real}(x::AbstractVector{Quaternion{T}},y::AbstractVector{Quaternion{T}})
    rsum = dot(real(x), real(y))
    isum = dot(imagi(x), imagi(y))
    jsum = dot(imagj(x), imagj(y))
    ksum = dot(imagk(x), imagk(y))
    return (rsum+isum+jsum+ksum)/length(x)*0.25
end
## }}}

## for multi states {{{
function overlapm{T<:Real}(x::AbstractVector{T},y::AbstractVector{T})
    return sum(x .== y)/length(x)
end

function overlapm{T<:Real}(x::AbstractVector{Complex{T}},y::AbstractVector{Complex{T}})
    rsum = sum(real(x).== real(y))
    isum = sum(imag(x).== imag(y))
    return (rsum+isum)/length(x)*0.5
end

function overlapm{T<:Real}(x::AbstractVector{Quaternion{T}},y::AbstractVector{Quaternion{T}})
    rsum = sum(real(x).== real(y))
    isum = sum(imagi(x).== imagi(y))
    jsum = sum(imagj(x).== imagj(y))
    ksum = sum(imagk(x).== imagk(y))
    return (rsum+isum+jsum+ksum)/length(x)*0.25
end
## }}}

ext: jl


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot(sin, 0, 25, Coord.cartesian(aspect_ratio=1.0))


text:
p=50;setprecision(p); setrounding(BigFloat, RoundUp);
L=BigFloat("10"); n=BigInt(100); iN=convert(Int64,n);
dx=BigFloat("2")*L/(n-1); x=convert(Array{Float64},-L:dx:L)
I=zeros(BigFloat, n, n);for i=1:n I[i,i]=BigFloat("1") end

V=TaylorBigF.calc_Vandermonde(x,n,n)
xi=(x.+dx/BigFloat("2"))[1:end-1]
Vi=TaylorBigF.calc_Vandermonde(xi,n-1,n)
dV=TaylorBigF.calc_Vandermonde_dot(x,n,n)
ddV=TaylorBigF.calc_Vandermonde_dot_dot(x,n,n)
dddV=TaylorBigF.calc_Vandermonde_dot_dot_dot(x,n,n)
ddddV=TaylorBigF.calc_Vandermonde_dot_dot_dot_dot(x,n,n)
intV=TaylorBigF.calc_Vandermonde_int(x,n,n)

f=sin.(x)+sin.(x./BigFloat("2"))
fi=sin.(xi)+sin.(xi./BigFloat("2"))
f+=convert(Array{BigFloat},(rand(n).-1/2)./1)
int_f=zeros(BigFloat,n)
for i=2:n
    int_f[i]=int_f[i-1]+(f[i]+f[i-1])*dx/2
end

d1=f[2:end]-f[1:end-1]
d2=(d1[2:end]+d1[1:end-1])/2
d=[d1[1];d2;d1[end]]

scatter(x,f)
plot!(x,f)

Eq1=V[1,:]
q1=f[1]
Eq2=V[end,:]
q2=f[end]

VV=[intV;Eq1';Eq2']
q_all=[f;q1;q2]
a=VV\q_all

r_root=Array{Float64}(undef,0); for i=1:Int64(n-1) push!(r_root,abs(a[i])^(1/(n-1))) end; r_rootF=convert(Array{Float64},r_root);
plot(r_rootF)
error1=maximum(abs.(f-V*a))
plot(x,f)
scatter!(x,V*a)
Vi=TaylorBigF.calc_Vandermonde(xi,n-1,BigInt(150));
fi_comp=Vi*a
error2=maximum(abs.(fi_comp-fi))
scatter!(xi,Vi*a)


text:
# Test sample functions

using StatsBase
using Base.Test

srand(1234)

function est_p(x, K)
	h = zeros(Int, K)
	for xi in x
		h[xi] += 1
	end
	p = h / length(x)
end

#### randi

const randi = StatsBase.randi
n = 1_000_000

x = [randi(10) for i = 1:n]
@test isa(x, Vector{Int})
@test minimum(x) == 1
@test maximum(x) == 10

x = [randi(3, 12) for i = 1:n]
@test isa(x, Vector{Int})
@test minimum(x) == 3
@test maximum(x) == 12

#### sample with replacement

## individual

x = [sample(1:5) for _ = 1:10^5]
p = fill(0.2, 5)
@test isa(x, Vector{Int})
@test_approx_eq_eps est_p(x, 5) p 0.02

## unordered

x = sample(1:5, 10^5)
@test isa(x, Vector{Int})
@test length(x) == 10^5
@test !issorted(x)  # extremely unlikely to be sorted for n=10^5
@test_approx_eq_eps est_p(x, 5) p 0.02
@test_approx_eq_eps est_p(x[1:50000], 5) p 0.03

## ordered

x = sample(1:5, 10^5; ordered=true)
@test isa(x, Vector{Int})
@test length(x) == 10^5
@test issorted(x)
@test_approx_eq_eps est_p(x, 5) p 0.02

#### sample without replacement

## unordered (using samplepair)

x = Array(Int, 2, n)
for i = 1:n
	x[:,i] = sample(1:5, 2; replace=false)
end

@test all(x[1,:] .!= x[2,:])
@test_approx_eq_eps est_p(x[1,:], 5) p 0.02
@test_approx_eq_eps est_p(x[2,:], 5) p 0.02

## unordered (using Fisher-Yates)

x = Array(Int, 3, n)
for i = 1:n
	x[:,i] = sample(1:5, 3; replace=false)
end

@test all(x[1,:] .!= x[2,:])
@test all(x[1,:] .!= x[3,:])
@test all(x[2,:] .!= x[3,:])
@test_approx_eq_eps est_p(x[1,:], 5) p 0.02
@test_approx_eq_eps est_p(x[2,:], 5) p 0.02
@test_approx_eq_eps est_p(x[3,:], 5) p 0.02

## unordered (using Self-avoid method)

x = Array(Int, 3, 10000)
for i = 1:10000
	x[:,i] = sample(1:1000, 3; replace=false)
end

@test all(x[1,:] .!= x[2,:])
@test all(x[1,:] .!= x[3,:])
@test all(x[2,:] .!= x[3,:])

## ordered

x = Array(Int, 3, n)
for i = 1:n
	x[:,i] = sample(1:5, 3; replace=false, ordered=true)
end

@test all(x[1,:] .< x[2,:])
@test all(x[2,:] .< x[3,:])
@test_approx_eq_eps est_p(vec(x), 5) p 0.02

#### weighted sample with replacement

wv = weights([1.0, 2.0, 4.5, 2.5])
p = [0.10, 0.20, 0.45, 0.25]

## individual

x = Int[sample(wv) for _ = 1:10^5]
@test isa(x, Vector{Int})
@test_approx_eq_eps est_p(x, 4) p 0.02

## unordered

x = sample(1:5, wv, 10^5)
@test isa(x, Vector{Int})
@test length(x) == 10^5
@test !issorted(x)
@test_approx_eq_eps est_p(x, 4) p 0.02

## ordered

x = sample(1:5, wv, 10^5; ordered=true)
@test isa(x, Vector{Int})
@test length(x) == 10^5
@test issorted(x)
@test_approx_eq_eps est_p(x, 4) p 0.02


text:
module MovcolNTests

using MovcolN
using FactCheck

# write your own tests here

facts("Generating Q") do
    kmax  = 4; nder = 3;
    Ql,Qr,x = MovcolN.generateQats(nder,1//2,kmax=kmax)

    @fact x => 1//2

    @fact size(Ql) => (kmax,nder)

    @fact Ql => [[1//2,5//32,1//64] [-15//8,-7//16,-1//32] [0//1,-3//2,-1//4] [30//1,15//1,3//2]]'

    @fact Qr => [[1//2,-5//32,1//64] [15//8,-7//16,1//32] [0//1,3//2,-1//4] [-30//1,15//1,-3//2]]'

end

facts("Computing derivatives") do
    kmax    = 4; nder = 3;
    Ql,Qr,_ = MovcolN.generateQats(nder,0.5,kmax=kmax)

    ul = hcat(Float64[ 0, 0, 2   ])
    ur = hcat(Float64[ 4, 4, 2   ])
    utl= hcat(Float64[ 0, 0, 2   ])
    utr= hcat(Float64[12, 8, 2   ])
    uc = hcat(Float64[ 1, 2, 2, 0])
    utc=      Float64[ 1, 2, 2   ]
    h  = 2.; ht = 2.
    H  = [h^j for j = 0:kmax]
    Ht = [j*ht*h^(j-1) for j = 0:kmax]

    uc1 = Array(Float64,kmax,1)
    utc1 = Array(Float64,kmax,1)

    MovcolN.computeux!(uc1,(Ql,Qr,0.5),H,ul,ur)
    MovcolN.computeutx!(utc1,(Ql,Qr,0.5),H,Ht,0.,utl,utr,uc,ul,ur)

    @fact uc1  => uc
    @fact utc1[1:3] => utc
end

facts("Generating AB matrix") do
    relerr = Array(Float64,0)
    abserr = Array(Float64,0)

    # for Float64 this method converges for up to 8 collocatin points
    nconv = 8
    for n = 2:nconv
        AB  =MovcolN.generateAB(n)
        sl,_=MovcolN.lobatto(n+1)
        sg,_=MovcolN.gauss(n)

        @fact size(AB) => (n,n+1)

        push!(abserr,norm((AB*sin(sl)-cos(sg)),         Inf))
        push!(relerr,norm((AB*sin(sl)-cos(sg))./sin(sg),Inf))
    end

    relconv = exp(diff(log(relerr)))
    absconv = exp(diff(log(abserr)))

    @fact relconv.<0.2 => all  # relative error is at most ~O(0.2^n)
    @fact absconv.<0.1 => all  # absolute error is at most ~O(0.1^n)
end

end


text:
function permittivity(w)
    # E.R. Peck and K. Reeder "Dispersion of Air" JOSA, 62, 958 (1972)
    if w == 0
        sig = 0
    else
        lam = 2 * pi * C0 / w
        sig = 1 / lam
        sig = sig * 1e-6   # 1/m -> 1/um
    end
    dum = 8060.51 + 2480990 / (132.247 - sig^2) + 17455.7 / (39.32957 - sig^2)
    n = 1 + dum * 1e-8
    eps = n^2
    return eps
end

function permeability(w)
    mu = 1
    return mu
end

n2 = 0.5e-23   # [m**2/W] nonlinear index
N0 = 2.5e25   # [1/m^3] neutrals density


text:
using MultipleScattering, Test

@testset "run examples and generate README.md" begin
    examplestoconvert = ["random_particles", "particles_in_circle"]
    folder = "../docs/src/example/"

    for str in examplestoconvert
        include(string(folder,str,"/",str,".jl"))
        @test true
    end

    # auto generate their README.md to guarantee the code in these READMEs works.

    include(folder*"generate-READMEs.jl")

    gernerate_README(examplestoconvert, folder)

    @test true
end


text:
# Control functions

function u_f(u, q)
    # convert torque input to wrench
    Qb = q[4:7]
    Q0 = q[11:14]
    Q2 = q[25:28]

    u0 = u[1]  # joint0
    u2 = u[2]  # joint2

    τb = [0; -u0-u2; 0]
    τ0 = [0; u0;     0]
    τ2 = [0; u2;     0]

    fb = f_applied(τb, l_cb)  # force in body frame
    f0 = f_applied(τ0, l_c0)
    f2 = f_applied(τ2, l_c2)

    #Corresponding wrench "F" for each link
    Fk = [rotate(Qb, fb); τb;        # body
          rotate(Q0, f0); τ0;        # link0
          zeros(3);       zeros(3);  # link1
          rotate(Q2, f2); τ2;        # link2
          zeros(3);       zeros(3)]  # link3

    return Fk
end

function f_applied(τ, r)
    cross(τ, r)/(norm(r)^2)
end

function a_control(a_target, a_pos, a_vel, q)
    # a_target: joint angle target
    # a_pos: joint angles
    kp = 5
    kd = copy(kp)*0.03
    # print("a_pos = ", a_pos.*180/pi, "\n")
    u = kp*(a_pos-a_target) + kd*(a_vel)
    # B = Diagonal([1; 0; 1; 0])  # actuator selection matrix
    Fk = u_f(u, q)  # convert torque input to wrench
    return Fk
end

function max_to_min(q)
    # convert max coords state to min coords
    rb = q[1:3]
    Qb = q[4:7]
    a = a_joint(q)
    return [rb; Qb; a]
end

# Kinematics
function kin_ee_min(q_min)
    # forward kinematics of the end effector from body pose and joint positions
    rb = q_min[1:3]
    Qb = q_min[4:7]
    a = q_min[8:11]
    q2 = a[3]
    q3 = a[4]

    pb = rb + rotate(Qb, l_cb)  # position vector from world frame to *JOINTS* 0 and 2

    Qb2 = Expq([0, q2, 0])  # quaternion from base to link 2
    Q2 = L(Qb)*Qb2  # quaternion from world frame to link 2

    Q23 = Expq([0, q3, 0])  # quaternion from link 2 to link 3
    Q3 = L(Q2)*Q23  # quaternion from world frame to link 2
    ree = pb + rotate(Q2, l2) + rotate(Q3, lee)

    return ree
end

max_stars_repo_name: bbokser/hopper-sim

max_issues_repo_name: bbokser/hopper-sim

max_forks_repo_name: bbokser/hopper-sim


text:
@testset "Conversion from and to SBML model" begin
    sbmlm = load_sbml_model(model_paths["ecoli_core_model.xml"])
    sm = convert(StandardModel, sbmlm)
    sbmlm2 = convert(SBMLModel, sm)

    @test Set(reactions(sbmlm)) == Set(reactions(sbmlm2))
    @test Set(reactions(sbmlm)) == Set(reactions(sm))
    @test Set(metabolites(sbmlm)) == Set(metabolites(sbmlm2))
    @test all([
        sbmlm.sbml.reactions[i].stoichiometry == sbmlm2.sbml.reactions[i].stoichiometry for
        i in reactions(sbmlm2)
    ])
end

@testset "SBMLModel generic interface" begin
    model = load_model(model_paths["e_coli_core.xml"])

    @test reaction_stoichiometry(model, "R_EX_ac_e") == Dict("M_ac_e" => -1)
end


text:
using LinearAlgebraicRepresentation,ViewerGL
Lar = LinearAlgebraicRepresentation
GL = ViewerGL

function randomarrangement2d(V,EV)
	VV = [[k] for k=1:size(V,2)]
	#GL.VIEW( GL.numbering(0.1)((V,[VV, EV])) );

	W = convert(Lar.Points, V')
	cop_EV = Lar.coboundary_0(EV::Lar.Cells)
	cop_EW = convert(Lar.ChainOp, cop_EV)
	V, copEV, copFE = Lar.Arrangement.planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp)

	triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
	FVs = convert(Array{Lar.Cells}, triangulated_faces)
	V = convert(Lar.Points, V')
	#GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,99));

	W, copEV, copFE = Lar.Arrangement.planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp)
	EVs = Lar.FV2EVs(copEV, copFE) # polygonal face fragments
	V = convert(Lar.Points, W')
	#GL.VIEW(GL.GLExplode(V,EVs,1.2,1.2,1.2,1));

	return V,EVs,FVs
end

V,EV = Lar.randomcuboids(100, .4)
V = GL.normalize2(V,flag=true)

W,EVs,FVs = randomarrangement2d(V,EV)
GL.VIEW(GL.GLExplode(W,EVs,1.2,1.2,1.2,1,1));
GL.VIEW(GL.GLExplode(W,FVs,1.2,1.2,1.2,99,1));
GL.VIEW(GL.GLExplode(W,FVs,1,1,1,99,1));


text:
"""
    AbstractSolver
Abstract type to encapsulate the different solver types such as `Solver` or `_SubSolver`.
"""
abstract type AbstractSolver end

# Dummy method to (not) add a TimeStamps to a solver
add_time!(::AbstractSolver, i) = nothing

function solver(ms, id, role; pool = pool(), strats = MetaStrategy(ms))
    mlid = make_id(meta_id(ms), id, Val(role))
    return solver(mlid, ms.model, ms.options, pool, ms.rc_report, ms.rc_sol, ms.rc_stop, strats, Val(role))
end

# Forwards from model field
@forward AbstractSolver.model add!
@forward AbstractSolver.model add_value!
@forward AbstractSolver.model add_var_to_cons!
@forward AbstractSolver.model constraint!
@forward AbstractSolver.model constriction
@forward AbstractSolver.model delete_value!
@forward AbstractSolver.model delete_var_from_cons!
@forward AbstractSolver.model describe
@forward AbstractSolver.model domain_size
@forward AbstractSolver.model draw
@forward AbstractSolver.model get_cons_from_var
@forward AbstractSolver.model get_constraint
@forward AbstractSolver.model get_constraints
@forward AbstractSolver.model get_domain
@forward AbstractSolver.model get_name
@forward AbstractSolver.model get_objective
@forward AbstractSolver.model get_objectives
@forward AbstractSolver.model get_variable
@forward AbstractSolver.model get_variables
@forward AbstractSolver.model get_vars_from_cons
@forward AbstractSolver.model is_sat
@forward AbstractSolver.model is_specialized
@forward AbstractSolver.model length_cons
@forward AbstractSolver.model length_objs
@forward AbstractSolver.model length_var
@forward AbstractSolver.model length_vars
@forward AbstractSolver.model max_domains_size
@forward AbstractSolver.model objective!
@forward AbstractSolver.model sense
@forward AbstractSolver.model sense!
@forward AbstractSolver.model state
@forward AbstractSolver.model update_domain!
@forward AbstractSolver.model variable!
@forward AbstractSolver.model _best_bound
@forward AbstractSolver.model _inc_cons!
@forward AbstractSolver.model _is_empty
@forward AbstractSolver.model _max_cons
@forward AbstractSolver.model _set_domain!

# Forwards from state field
@forward AbstractSolver.state get_error
@forward AbstractSolver.state get_value
@forward AbstractSolver.state get_values
@forward AbstractSolver.state set_error!
@forward AbstractSolver.state set_value!
@forward AbstractSolver.state _best
@forward AbstractSolver.state _best!
@forward AbstractSolver.state _cons_cost
@forward AbstractSolver.state _cons_cost!
@forward AbstractSolver.state _cons_costs
@forward AbstractSolver.state _cons_costs!
@forward AbstractSolver.state _inc_last_improvement!
@forward AbstractSolver.state _last_improvement
@forward AbstractSolver.state _optimizing
@forward AbstractSolver.state _optimizing!
@forward AbstractSolver.state _satisfying!
@forward AbstractSolver.state _set!
@forward AbstractSolver.state _solution
@forward AbstractSolver.state _swap_value!
@forward AbstractSolver.state _reset_last_improvement!
@forward AbstractSolver.state _value
@forward AbstractSolver.state _value!
@forward AbstractSolver.state _values
@forward AbstractSolver.state _values!
@forward AbstractSolver.state _var_cost
@forward AbstractSolver.state _var_cost!
@forward AbstractSolver.state _vars_costs
@forward AbstractSolver.state _vars_costs!

# Forward from options
@forward AbstractSolver.options get_option
@forward AbstractSolver.options set_option!
@forward AbstractSolver.options _verbose

# Forwards from pool (of solutions)
@forward AbstractSolver.pool best_config
@forward AbstractSolver.pool best_value
@forward AbstractSolver.pool best_values
@forward AbstractSolver.pool has_solution

# Forwards from strategies
@forward AbstractSolver.strategies check_restart!
@forward AbstractSolver.strategies decay_tabu!
@forward AbstractSolver.strategies decrease_tabu!
@forward AbstractSolver.strategies delete_tabu!
@forward AbstractSolver.strategies empty_tabu!
@forward AbstractSolver.strategies insert_tabu!
@forward AbstractSolver.strategies length_tabu
@forward AbstractSolver.strategies tabu_list

"""
    specialize!(solver)
Replace the model of `solver` by one with specialized types (variables, constraints, objectives).
"""
specialize!(s) = s.model = specialize(s.model)

"""
    _draw!(s)
Draw a random (re-)starting configuration.
"""
function _draw!(s)
    foreach(x -> _set!(s, x, draw(s, x)), keys(get_variables(s)))
end

"""
    _compute_cost!(s, ind, c)

Compute the cost of constraint `c` with index `ind`.
"""
function _compute_cost!(s, ind, c)
    old_cost = _cons_cost(s, ind)
    new_cost = compute_cost(c, _values(s), s.state.icn_computations)
    _cons_cost!(s, ind, new_cost)
    foreach(x -> _var_cost!(s, x, _var_cost(s, x) + new_cost - old_cost), c.vars)
end

"""
    _compute_costs!(s; cons_lst::Indices{Int} = Indices{Int}())

Compute the cost of constraints `c` in `cons_lst`. If `cons_lst` is empty, compute the cost for all the constraints in `s`.
"""
function _compute_costs!(s; cons_lst=Indices{Int}())
    if isempty(cons_lst)
        foreach(((id, c),) -> _compute_cost!(s, id, c), pairs(get_constraints(s)))
    else
        foreach(
            ((id, c),) -> _compute_cost!(s, id, c),
            pairs(view(get_constraints(s), cons_lst))
        )
    end
    set_error!(s, sum(_cons_costs(s)))
end

"""
    _compute_objective!(s, o::Objective)
    _compute_objective!(s, o = 1)

Compute the objective `o`'s value.
"""
function _compute_objective!(s, o::Objective)
    val = sense(s) * apply(o, _values(s).values)
    set_value!(s, val)
    if is_empty(s.pool) || val < best_value(s)
        s.pool = pool(s.state.configuration)
    end
end
_compute_objective!(s, o=1) = _compute_objective!(s, get_objective(s, o))

"""
    _compute!(s; o::Int = 1, cons_lst = Indices{Int}())

Compute the objective `o`'s value if `s` is satisfied and return the current `error`.

# Arguments:
- `s`: a solver
- `o`: targeted objective
- `cons_lst`: list of targeted constraints, if empty compute for the whole set
"""
function _compute!(s; o::Int=1, cons_lst=Indices{Int}())
    _compute_costs!(s; cons_lst)
    if get_error(s) == 0.0
        _optimizing(s) && _compute_objective!(s, o)
        is_sat(s) && (s.pool = pool(s.state.configuration))
        return true
    end
    return false
end

"""
    _neighbours(s, x, dim = 0)

DOCSTRING

# Arguments:
- `s`: DESCRIPTION
- `x`: DESCRIPTION
- `dim`: DESCRIPTION
"""
function _neighbours(s, x, dim = 0)
    if dim == 0
        is_discrete = typeof(get_variable(s, x).domain) <: ContinuousDomain
        return is_discrete ? map(_ -> draw(s, x), 1:(length_vars(s)*length_cons(s))) : get_domain(s, x)
    else
        neighbours = Set{Int}()
        foreach(
            c -> foreach(y ->
                begin
                    b = _value(s, x) ∈ get_variable(s, y) && _value(s, y) ∈ get_variable(s, x)
                    b && push!(neighbours, y)
                end, get_vars_from_cons(s, c)),
            get_cons_from_var(s, x)
        )
        return delete!(neighbours, x)
    end
end

state!(s) = s.state = state(s)

function _init!(s, ::Val{:global})
    if !is_specialized(s) && get_option(s, "specialize")
        specialize!(s)
        set_option!(s, "specialize", true)
    end
    put!(s.rc_stop, nothing)
    foreach(i -> put!(s.rc_report, nothing), setdiff(workers(), [1]))
end
_init!(s, ::Val{:meta}) = foreach(id -> push!(s.subs, solver(s, id-1, :sub)), 2:nthreads())

function _init!(s, ::Val{:remote})
    for w in setdiff(workers(), [1])
        ls = remotecall(solver, w, s, w, :lead)
        remote_do(set_option!, w, fetch(ls), "print_level", :silent)
        remote_do(set_option!, w, fetch(ls), "threads",remotecall_fetch(Threads.nthreads, w))
        push!(s.remotes, w => ls)
    end
end

function _init!(s, ::Val{:local}; pool = pool())
    get_option(s, "tabu_time") == 0 && set_option!(s, "tabu_time", length_vars(s) ÷ 2) # 10?
    get_option(s, "tabu_local") == 0 && set_option!(s, "tabu_local", get_option(s, "tabu_time") ÷ 2)
    get_option(s, "tabu_delta") == 0 && set_option!(s, "tabu_delta", get_option(s, "tabu_time") - get_option(s, "tabu_local")) # 20-30
    state!(s)
    return has_solution(s)
end

_init!(s, role::Symbol) = _init!(s, Val(role))

"""
    _restart!(s, k = 10)

Restart a solver.
"""
function _restart!(s, k=10)
    _verbose(s, "\n============== RESTART!!!!================\n")
    _draw!(s)
    empty_tabu!(s)
    δ = ((k - 1) * get_option(s, "tabu_delta")) + get_option(s, "tabu_time") / k
    set_option!(s, "tabu_delta", δ)
    _compute!(s) ? _optimizing!(s) : _satisfying!(s)
end

"""
    _check_restart(s)

Check if a restart of `s` is necessary. If `s` has subsolvers, this check is independent for all of them.
"""
function _check_restart(s)
    a = _last_improvement(s) > length_vars(s)
    b = check_restart!(s; tabu_length = length_tabu(s))
    return a || b
end

"""
    _select_worse(s::S) where S <: Union{_State, AbstractSolver}
Within the non-tabu variables, select the one with the worse error .
"""
function _select_worse(s)
    nontabu = setdiff(keys(_vars_costs(s)), keys(tabu_list(s)))
    return _find_rand_argmax(view(_vars_costs(s), nontabu))
end

"""
    _move!(s, x::Int, dim::Int = 0)

Perform an improving move in `x` neighbourhood if possible.

# Arguments:
- `s`: a solver of type S <: AbstractSolver
- `x`: selected variable id
- `dim`: describe the dimension of the considered neighbourhood
"""
function _move!(s, x::Int, dim::Int=0)
    best_values = [begin old_v = _value(s, x) end]; best_swap = [x]
    tabu = true # unless proved otherwise, this variable is now tabu
    best_cost = old_cost = get_error(s)
    copy_to!(s.state.fluct, _cons_costs(s), _vars_costs(s))
    for v in _neighbours(s, x, dim)
        dim == 0 && v == old_v && continue
        dim == 0 ? _value!(s, x, v) : _swap_value!(s, x, v)

        _verbose(s, "Compute costs: selected var(s) x_$x " * (dim == 0 ? "= $v" : "⇆ x_$v"))

        cons_x_v = union(get_cons_from_var(s, x), dim == 0 ? [] : get_cons_from_var(s, v))
        _compute!(s, cons_lst=cons_x_v)

        cost = get_error(s)
        if cost < best_cost
            _verbose(s, "cost = $cost < $best_cost")
            tabu = false
            best_cost = cost
            dim == 0 ? best_values = [v] : best_swap = [v]
        elseif cost == best_cost
            _verbose(s, "cost = best_cost = $cost")
            push!(dim == 0 ? best_values : best_swap, v)
        end

        if cost == 0 && is_sat(s)
            s.pool == pool(s.state.configuration)
            return best_values, best_swap, tabu
        end

        copy_from!(s.state.fluct, _cons_costs(s), _vars_costs(s))
        set_error!(s, old_cost)

        # swap/change back the value of x (and y/)
        dim == 0 ? _value!(s, x, old_v) : _swap_value!(s, x, v)
    end
    return best_values, best_swap, tabu
end

"""
    _step!(s)

Iterate a step of the solver run.
"""
function _step!(s)
    # select worst variables
    x = _select_worse(s)
    _verbose(s, "Selected x = $x")

    # Local move (change the value of the selected variable)
    best_values, best_swap, tabu = _move!(s, x)
    # _compute!(s)

    # If local move is bad (tabu), then try permutation
    if tabu
        _, best_swap, tabu = _move!(s, x, 1)
        _compute!(s)
    else # compute the costs changes from best local move
        _compute!(s; cons_lst=get_cons_from_var(s, x))
    end

    # decay tabu list
    decay_tabu!(s)

    # update tabu list with either worst or selected variable
    insert_tabu!(s, x, tabu ? :tabu : :pick)
    _verbose(s, "Tabu list: $(tabu_list(s))")

    # Inc last improvement if tabu
    tabu ? _inc_last_improvement!(s) : _reset_last_improvement!(s)

    # Select the best move (value or swap)
    if x ∈ best_swap
        _value!(s, x, rand(best_values))
        _verbose(s, "best_values: $best_values")
    else
        _swap_value!(s, x, rand(best_swap))
        _verbose(s, "best_swap : $best_swap")
    end

    # Compute costs and possibly evaluate objective functions
    # return true if a solution for sat is found
    if _compute!(s)
        !is_sat(s) ? _optimizing!(s) : return true
    end

    # Restart if necessary
    _check_restart(s) && _restart!(s)

    return false # no satisfying configuration or optimizing
end

"""
    _check_subs(s)

Check if any subsolver of a main solver `s`, for
- *Satisfaction*, has a solution, then return it, resume the run otherwise
- *Optimization*, has a better solution, then assign it to its internal state
"""
_check_subs(::AbstractSolver) = 0 # Dummy method

"""
    stop_while_loop()
Check the stop conditions of the `solve!` while inner loop.
"""
stop_while_loop(::AbstractSolver) = nothing

"""
    solve_while_loop!(s, )
Search the space of configurations.
"""
function solve_while_loop!(s, stop, sat, iter, st)
    while stop_while_loop(s, stop, iter, st)
        iter += 1
        _verbose(s, "\n\tLoop $(iter) ($(_optimizing(s) ? "optimization" : "satisfaction"))")
        _step!(s) && sat && break
        _verbose(s, "vals: $(length(_values(s)) > 0 ? _values(s) : nothing)")
        best_sub = _check_subs(s)
        if best_sub > 0
            bs = s.subs[best_sub]
            s.pool = deepcopy(bs.pool)
            sat && break
        end
    end
end

"""
    remote_dispatch!(solver)
Starts the `LeadSolver`s attached to the `MainSolver`.
"""
remote_dispatch!(::AbstractSolver) = nothing # dummy method

"""
    solve_for_loop!(solver, stop, sat, iter)
First loop in the solving process that starts `LeadSolver`s from the `MainSolver`, and `_SubSolver`s from each `MetaSolver`.
"""
solve_for_loop!(s, stop, sat, iter, st) = solve_while_loop!(s, stop, sat, iter, st)

function update_pool!(s, pool)
    is_empty(pool) && return nothing
    if is_sat(s) || best_value(s) > best_value(pool)
        s.pool = deepcopy(pool)
    end
end

"""
    remote_stop!!(solver)
Fetch the pool of solutions from `LeadSolvers` and merge it into the `MainSolver`.
"""
remote_stop!(::AbstractSolver) = nothing

"""
    post_process(s::MainSolver)
Launch a serie of tasks to round-up a solving run, for instance, export a run's info.
"""
post_process(::AbstractSolver) = nothing

function solve!(s, stop = Atomic{Bool}(false))
    start_time = time()
    add_time!(s, 1) # only used by MainSolver
    iter = 0 # only used by MainSolver
    sat = is_sat(s)
    _init!(s) && (sat ? (iter = typemax(0)) : _optimizing!(s))
    add_time!(s, 2) # only used by MainSolver
    solve_for_loop!(s, stop, sat, iter, start_time)
    add_time!(s, 5) # only used by MainSolver
    remote_stop!(s)
    add_time!(s, 6) # only used by MainSolver
    post_process(s) # only used by MainSolver
end

"""
    solution(s)
Return the only/best known solution of a satisfaction/optimization model.
"""
solution(s) = has_solution(s) ? best_values(s) : _values(s)


text:
using Memoize

abstract type Descriptor{T<:Number, DN} end

Base.eltype(::Descriptor{T, DN}) where {T, DN} = T

"""
idx : Index of self
specie : Periodic number of self
position : Cartesian position of self
envidcs : Index of atoms within BallTree
envnumbers : Cartesian coordinates of atoms within BallTree
envpositions : Periodic number of each atoms
cutoff
descriptor : cache
"""
struct AtomicDescriptor{T, DN} <: Descriptor{T, DN}
    idx::Int
    number::Int
    position::Array{T, 1}
    ##
    # numbers
    # cell
    # pbc
    # balltree
    #
    ## below I want to delte it later
    envidcs::Array{Int, 1}
    envnumbers::Array{Int, 1}
    envpositions::Array{T, 2}
    cutoff::Number
    descriptor::Dict
end

@inline AtomicDescriptor() = AtomicDescriptor{Float64, 1}(1, 1, [0., 0., 0.], [1], [1], zeros(1, 3), 7, Dict())

Base.convert(::Type{AtomicDescriptor{T, DN}}, not::Nothing) where {T<:Number, DN} = AtomicDescriptor()
Base.convert(::Type{AtomicDescriptor{T}}, not::Nothing) where {T<:Number} = AtomicDescriptor()
Base.convert(::Type{AtomicDescriptor}, not::Nothing) = AtomicDescriptor()

#### <!--- Utility function neeeded for manybody descriptors
function degeneracy(cluster)
    idx = Dict()
    multiple = Dict()
    for (i, a) in enumerate(cluster)
        if get(multiple, a, nothing) == nothing
            multiple[a] = 0
            idx[a] = []
        end
        multiple[a] += 1
        push!(idx[a], i)
    end

    total_mul = 1
    for (k, v) in pairs(multiple)
        total_mul *= factorial(v)
    end

    permutedidx = Dict()
    for s in keys(idx)
        permutedidx[s] = permutations(idx[s])
    end

    allpermutations = []
    orgidx = 1:length(cluster) |> collect
    for combinations in Iterators.product(values(permutedidx)...)
        copied = copy(orgidx)
        for (s, combo) in zip(keys(permutedidx), combinations)
            copied[idx[s]] = combo
        end

        push!(allpermutations, copied)
    end

    # return multiple, total_mul, allpermutations
    return allpermutations
end

@memoize function getchart(numbers::Tuple, n::Int)
    species = Set(numbers) |> collect
    chart = Vector{Tuple}([])
    # multiplicity = Dict()
    permutations = Dict()
    for cluster in with_replacement_combinations(species, n)
        srtdcluster = Tuple(sort(cluster))
        push!(chart, srtdcluster)
        permutations[srtdcluster] = degeneracy(srtdcluster)
    end
    return (chart, permutations)
end

macro ini(chart)
    return :(Dict([(key, Vector{Number}([])) for key in $chart]))
end
#### -->

"""
idx : Index of self
specie : Periodic number of self
position : Cartesian position of self
envidcs : Index of atoms within BallTree
envnumbers : Cartesian coordinates of atoms within BallTree except self
envpositions : Periodic number of each atoms except self
cutoff
descriptor : cache
"""
struct TwoBody{T, DN} <: Descriptor{T, DN}
    idx::Int
    number::Int
    position::Array{T, 1}
    envidcs::Array{Int, 1}
    envnumbers::Array{Int, 1}
    envpositions::Array{T, 2}
    cutoff::Number
    descriptor::Dict
end

@inline TwoBody() = TwoBody{Float64, 1}(1, 1, [0., 0., 0.], [1], [1], zeros(1, 3), 7, Dict())

Base.convert(::Type{TwoBody{T, DN}}, not::Nothing) where {T<:Number, DN} = TwoBody()
Base.convert(::Type{TwoBody{T}}, not::Nothing) where {T<:Number} = TwoBody()
Base.convert(::Type{TwoBody}, not::Nothing) = TwoBody()

Base.copy(d::TwoBody{T, DN}) where {T<:Number, DN} = #=
    =# TwoBody{T, DN}(d.idx, d.number, copy(d.position), copy(d.envidcs),
                      copy(d.envnumbers), copy(d.envpositions), d.cutoff,
                      copy(d.descriptor))

struct TwoBodyFromAtomicDescriptor <: Transformation; twobodycutoff; end
Base.convert(::Type{TwoBody}, atomic::AtomicDescriptor) = TwoBodyFromAtomicDescriptor(twobodycutoff)(atomic)
function (twobody::TwoBodyFromAtomicDescriptor)(atomic)
    atom = atomic
    number = atom.number
    numbers = atom.descriptor["numbers"]
    # chart = atom.descriptor["chart"]
    A = length(atom.descriptor["numbers"])
    idx = atom.idx
    position = atom.position

    balltree = atom.descriptor["balltree"]
    miccell = atom.descriptor["miccell"]

    envidcs_plusself = inrange(balltree, position, twobody.twobodycutoff, true)

    selfidx = searchsortedfirst(envidcs_plusself, idx)
    envidcs = deleteat!(envidcs_plusself, selfidx)
    envnumbers = numbers[(envidcs .% A) .+ 1]
    envpositions = view(miccell, :, envidcs)'

    M = length(envnumbers)

    # Create a dictionary maps sends species to distances
    # print(getchart(Tuple(numbers), 2))
    chart, degeneracy = getchart(Tuple(numbers), 2)
    descriptor = Dict([(
        key, Dict("idx"   =>Vector{Number}([]),
                  "xi"    =>Vector{Number}([]),
                  "yi"    =>Vector{Number}([]),
                  "zi"    =>Vector{Number}([]),
                  "ri"    =>Array{Array{Float64, 1}}([]),
                  "xi1"   =>Vector{Number}([]),
                  "yi1"   =>Vector{Number}([]),
                  "zi1"   =>Vector{Number}([]),
                  "ri1"   =>Array{Array{Float64, 1}}([]),
                  "xi_i1" =>Vector{Number}([]),
                  "yi_i1" =>Vector{Number}([]),
                  "zi_i1" =>Vector{Number}([]),
                  "rii1"  =>Vector{Number}([]),
                  "permutations" => degeneracy[key])) for key in chart])

    s1 = number
    for j=1:M
        s2 = envnumbers[j]

        if s1 <= s2
            cluster = (s1, s2)
            p1, p2 = position, envpositions[j, :]
        else
            cluster = (s2, s1)
            p1, p2 = envpositions[j, :], position
        end

        xi, yi, zi = p1
        xi1, yi1, zi1 = p2
        xi_i1, yi_i1, zi_i1 = xi-xi1, yi-yi1, zi-zi1
        rii1 = sqrt(xi_i1*xi_i1 + yi_i1*yi_i1 + zi_i1*zi_i1)

        push!(descriptor[cluster]["xi"], xi)
        push!(descriptor[cluster]["yi"], yi)
        push!(descriptor[cluster]["zi"], zi)

        push!(descriptor[cluster]["ri"], p1)

        push!(descriptor[cluster]["xi1"], xi1)
        push!(descriptor[cluster]["yi1"], yi1)
        push!(descriptor[cluster]["zi1"], zi1)

        push!(descriptor[cluster]["ri1"], p2)

        push!(descriptor[cluster]["xi_i1"], xi_i1)
        push!(descriptor[cluster]["yi_i1"], yi_i1)
        push!(descriptor[cluster]["zi_i1"], zi_i1)

        push!(descriptor[cluster]["rii1"], rii1)
    end
    TwoBody{eltype(atomic), 1}(idx, number,position, envidcs,
                    envnumbers, envpositions, twobody.twobodycutoff, descriptor)
end

"""
idx : Index of self
specie : Periodic number of self
position : Cartesian position of self
envidcs : Index of atoms within BallTree
envnumbers : Cartesian coordinates of atoms within BallTree exept self
envpositions : Periodic number of each atoms exept self
cutoff
descriptor : cache
"""
struct ThreeBody{T, DN} <: Descriptor{T, DN}
    idx::Int
    number::Int
    position::Array{T, 1}
    envidcs::Array{Int, 1}
    envnumbers::Array{Int, 1}
    envpositions::Array{T, 2}
    cutoff::Number
    descriptor::Dict
end

@inline ThreeBody() = ThreeBody{Float64, 1}(1, 1, [0., 0., 0.], [1], [1], zeros(1, 3), 7, Dict())

Base.convert(::Type{ThreeBody{T, DN}}, not::Nothing) where {T<:Number, DN} = ThreeBody()
Base.convert(::Type{ThreeBody{T}}, not::Nothing) where {T<:Number} = ThreeBody()
Base.convert(::Type{ThreeBody}, not::Nothing) = ThreeBody()

Base.copy(d::ThreeBody{T, DN}) where {T<:Number, DN} = #=
=# ThreeBody{T, DN}(d.idx, d.number, copy(d.position), copy(d.envidcs),
      copy(d.envnumbers), copy(d.envpositions), d.cutoff, copy(d.descriptor))

struct ThreeBodyFromAtomicDescriptor; threebodycutoff; end
Base.convert(::Type{ThreeBody}, atomic::AtomicDescriptor) = ThreeBodyFromAtomicDescriptor(threebodycutoff)(atomic)

function (threebody::ThreeBodyFromAtomicDescriptor)(atomic)
    atom = atomic
    number = atom.number
    numbers = atom.descriptor["numbers"]
    A = length(atom.descriptor["numbers"])
    idx = atom.idx
    position = atom.position

    balltree = atom.descriptor["balltree"]
    miccell = atom.descriptor["miccell"]

    envidcs_plusself = inrange(balltree, position, threebody.threebodycutoff, true)

    selfidx = searchsortedfirst(envidcs_plusself, idx)
    envidcs = deleteat!(envidcs_plusself, selfidx)
    envnumbers = numbers[(envidcs .% A) .+ 1]
    envpositions = view(miccell, :, envidcs)'

    M = length(envnumbers)

    # Create a dictionary maps sends species to distances
    chart, degeneracy = getchart(Tuple(numbers), 3)
    descriptor = Dict([(
        key, Dict("idx"   =>Vector{Number}([]),
                  "xi"    =>Vector{Number}([]),
                  "yi"    =>Vector{Number}([]),
                  "zi"    =>Vector{Number}([]),
                  "ri"    =>Array{Array{Float64, 1}, 1}([]),
                  "xi1"   =>Vector{Number}([]),
                  "yi1"   =>Vector{Number}([]),
                  "zi1"   =>Vector{Number}([]),
                  "ri1"    =>Array{Array{Float64, 1}, 1}([]),
                  "xi2"   =>Vector{Number}([]),
                  "yi2"   =>Vector{Number}([]),
                  "zi2"   =>Vector{Number}([]),
                  "ri2"    =>Array{Array{Float64, 1}, 1}([]),
                  "xi_i1" =>Vector{Number}([]),
                  "yi_i1" =>Vector{Number}([]),
                  "zi_i1" =>Vector{Number}([]),
                  "xi_i2" =>Vector{Number}([]),
                  "yi_i2" =>Vector{Number}([]),
                  "zi_i2" =>Vector{Number}([]),
                  "xi1_i2"=>Vector{Number}([]),
                  "yi1_i2"=>Vector{Number}([]),
                  "zi1_i2"=>Vector{Number}([]),
                  "rii1"  =>Vector{Number}([]),
                  "rii2"  =>Vector{Number}([]),
                  "ri1i2" =>Vector{Number}([]),
                  "permutations" => degeneracy[key])) for key in chart])
    # descriptor["chart"] = chart
    # descriptor["degeneracy"] = degeneracy

    s1 = number
    for j=1:M
        s2 = envnumbers[j]

        if s1 <= s2
            cluster = (s1, s2)
            p1, p2 = position, envpositions[j, :]
        else
            cluster = (s2, s1)
            s2, s1 = s1, s2
            p1, p2 = envpositions[j, :], position
        end

        xi, yi, zi = p1
        xi1, yi1, zi1 = p2
        xi_i1, yi_i1, zi_i1 = xi-xi1, yi-yi1, zi-zi1
        _rii1 = sqrt(xi_i1*xi_i1 + yi_i1*yi_i1 + zi_i1*zi_i1)

        for k=j+1:M
            s3 = envnumbers[k]
            p3 = envpositions[k, :]
            if s2 <= s3
                cluster = (s1, s2, s3)
                xi2, yi2, zi2 = p3

                xi_i2, yi_i2, zi_i2 = xi-xi2, yi-yi2, zi-zi2
                xi1_i2, yi1_i2, zi1_i2 = xi1-xi2, yi1-yi2, zi1-zi2
                # rii1 = sqrt(xi_i1*xi_i1 + yi_i1*yi_i1 + zi_i1*zi_i1)
                rii1 = _rii1
                rii2 = sqrt(xi_i2*xi_i2 + yi_i2*yi_i2 + zi_i2*zi_i2)
                ri1i2 = sqrt(xi1_i2*xi1_i2 + yi1_i2*yi_i2 + zi1_i2*zi_i2)
            elseif s3 <= s1
                cluster = (s3, s1, s2)
                xi, yi, zi = p3
                xi1, yi1, zi1 = p1
                xi2, yi2, zi2 = p2

                xi_i2, yi_i2, zi_i2 = xi-xi2, yi-yi2, zi-zi2
                xi1_i2, yi1_i2, zi1_i2 = xi_i1, yi_i1, zi_i1
                xi_i1, yi_i1, zi_i1 = xi-xi1, yi-yi1, zi-zi1

                rii1 = sqrt(xi_i1*xi_i1 + yi_i1*yi_i1 + zi_i1*zi_i1)
                rii2 = sqrt(xi_i2*xi_i2 + yi_i2*yi_i2 + zi_i2*zi_i2)
                ri1i2 = _rii1

                xi_i2, yi_i2, zi_i2 = xi-xi2, yi-yi2, zi-zi2
            else
                cluster = (s1, s3, s2)
                xi1, yi1, zi1 = p3
                xi2, yi2, zi2 = p2

                xi_i2, yi_i2, zi_i2 = xi_i1, yi_i1, zi_i1
                xi_i1, yi_i1, zi_i1 = xi-xi2, yi-yi2, zi-zi2
                xi1_i2, yi1_i2, zi1_i2 = xi1-xi2, yi1-yi2, zi1-zi2

                rii1 = sqrt(xi_i1*xi_i1 + yi_i1*yi_i1 + zi_i1*zi_i1)
                rii2 = _rii1
                ri1i2 = sqrt(xi1_i2*xi_i2 + yi1_i2*yi1_i2 + zi1_i2*zi1_i2)
            end
            push!(descriptor[cluster]["xi"], xi)
            push!(descriptor[cluster]["yi"], yi)
            push!(descriptor[cluster]["zi"], zi)
            push!(descriptor[cluster]["ri"], p1)

            push!(descriptor[cluster]["xi1"], xi1)
            push!(descriptor[cluster]["yi1"], yi1)
            push!(descriptor[cluster]["zi1"], zi1)
            push!(descriptor[cluster]["ri1"], p2)

            push!(descriptor[cluster]["xi2"], xi2)
            push!(descriptor[cluster]["yi2"], yi2)
            push!(descriptor[cluster]["zi2"], zi2)
            push!(descriptor[cluster]["ri2"], p3)

            push!(descriptor[cluster]["xi_i1"], xi_i1)
            push!(descriptor[cluster]["yi_i1"], yi_i1)
            push!(descriptor[cluster]["zi_i1"], zi_i1)

            push!(descriptor[cluster]["xi_i2"], xi_i2)
            push!(descriptor[cluster]["yi_i2"], yi_i2)
            push!(descriptor[cluster]["zi_i2"], zi_i2)

            push!(descriptor[cluster]["xi1_i2"], xi1_i2)
            push!(descriptor[cluster]["yi1_i2"], yi1_i2)
            push!(descriptor[cluster]["zi1_i2"], zi1_i2)

            push!(descriptor[cluster]["rii1"], rii1)
            push!(descriptor[cluster]["rii2"], rii2)
            push!(descriptor[cluster]["ri1i2"], ri1i2)

        end
    end
    ThreeBody{eltype(atomic), 1}(idx, number, position, envidcs, envnumbers,
          envpositions, threebody.threebodycutoff, descriptor)
end

"""
idx : Index of self
specie : Periodic number of self
position : Cartesian position of self
envidcs : Index of atoms within BallTree
envnumbers : Cartesian coordinates of atoms within BallTree
envpositions : Periodic number of each atoms
cutoff
descriptor : cache
"""
struct TwoThreeBody{T, DN} <: Descriptor{T, DN}
    twobody::TwoBody
    threebody::ThreeBody
    descriptor::Dict
    # idx::Int
    # number::Int
    # position::Array{T, 1}
    # envidcs::Array{Int, 1}
    # envnumbers::Array{Int, 1}
    # envpositions::Array{T, 2}
    # cutoff::Number
    # descriptor::Dict
end

@inline TwoThreeBody() = TwoThreeBody{Float64, 1}(TwoBody(), ThreeBody(), Dict())

Base.convert(::Type{TwoThreeBody{T, DN}}, not::Nothing) where {T<:Number, DN} = TwoThreeBody()
Base.convert(::Type{TwoThreeBody{T}}, not::Nothing) where {T<:Number} = TwoThreeBody()
Base.convert(::Type{TwoThreeBody}, not::Nothing) = TwoThreeBody()

Base.copy(d::TwoThreeBody{T, DN}) where {T<:Number, DN} = #=
=# TwoThreeBody{T, DN}(copy(d.twobody), copy(d.threebody), copy(d.descriptor))

struct TwoThreeBodyFromAtomicDescriptor; twobodycutoff; threebodycutoff; end
Base.convert(::Type{TwoThreeBody}, atomic::AtomicDescriptor;
             twobodycutoff=7, threebodycutoff=3.5) =  #=
    =#TwoThreeBodyFromAtomicDescriptor(twobodycutoff, threebodycutoff)(atomic)

function (twothreebody::TwoThreeBodyFromAtomicDescriptor)(atomic)

    twobody = TwoBodyFromAtomicDescriptor(twothreebody.twobodycutoff)(atomic)
    threebody = ThreeBodyFromAtomicDescriptor(twothreebody.threebodycutoff)(atomic)
    dtype = eltype(atomic)
    descriptor = Dict()
    for (k, v) in pairs(twobody.descriptor)
        descriptor[k] = v
    end
    for (k, v) in pairs(threebody.descriptor)
        descriptor[k] = v
    end
    TwoThreeBody{dtype, 1}(twobody, threebody, descriptor)
end

max_stars_repo_name: schinavro/taps

max_issues_repo_name: schinavro/taps

max_forks_repo_name: schinavro/taps


text:
# FIFO queue

mutable struct Queue{T}
    store::Deque{T}
end

"""
    Queue{T}([blksize::Integer=1024])

Create a `Queue` object containing elements of type `T`.
"""
Queue{T}() where {T} = Queue(Deque{T}())
Queue{T}(blksize::Integer) where {T} = Queue(Deque{T}(blksize))

@deprecate Queue(::Type{T}) where {T} Queue{T}()
@deprecate Queue(::Type{T}, blksize::Integer) where {T} Queue{T}(blksize)

isempty(s::Queue) = isempty(s.store)
length(s::Queue) = length(s.store)

front(s::Queue) = front(s.store)
back(s::Queue) = back(s.store)

"""
    enqueue!(s::Queue, x)

Inserts the value `x` to the end of the queue `s`.
"""
function enqueue!(s::Queue, x)
    push!(s.store, x)
    s
end

"""
    dequeue!(s::Queue)

Removes an element from the front of the queue `s` and returns it.
"""
dequeue!(s::Queue) = popfirst!(s.store)

# Iterators

start(q::Queue) = start(q.store)
next(q::Queue, s::Tuple) = next(q.store, s)
done(q::Queue, s::Tuple) = done(q.store, s)

reverse_iter(q::Queue) = reverse_iter(q.store)


text:
# ~/znb/corona/nb/corona.jl  -- downloaded from ~/znb/corona/nb/corona.ipynb
# 200515 004

#= working (notebook) directory: ~/znb/corona/nb

We start out by getting the data running debian/linux bash shell:
cd ~/znb/corona/
$ clone https://github.com/bumbeishvili/covid19-daily-data

# For later updating:
$ cd ~/znb/corona/covid19-daily-data
$ git pull origin

# These are the files we get:
$ ls ~/znb/corona/covid19-daily-data
misc  README.md  time_series_19-covid-Confirmed.csv  time_series_19-covid-Deaths.csv  time_series_19-covid-Recovered.csv

# Returning to working (notebook) directory
$ cd ~/znb/corona/nb
Now we can look into the time series of confirmed cases, running Julia:
=#

using CSV, Missings, DataFrames, Dates,  Printf, Plots
using SpecialFunctions, Distributions, RandomNumbers

df_cfm_src = CSV.read("../covid19-daily-data/time_series_19-covid-Confirmed.csv");
first(df_cfm_src, 10)

names(df_cfm_src)[end]

last(df_cfm_src)

countries = df_cfm_src[:,2];

# Specify which countries to extract:
clist  = ["Argentina" "Belgium" "Brazil" "Norway" "Spain" "Sweden" "Thailand"]

# Extracting the epidemic data:
cindx = Array{Int64}(undef, length(clist))
cfm   = zeros(size(df_cfm_src)[2]-4, length(clist)); # Data starts in column 4
for i in 1:length(clist)
    indx = findall(x->x==clist[i], df_cfm_src[:,2])
    cfm[:,i] = convert(Array, df_cfm_src[indx,5:end])
end
df_cfm = DataFrame(cfm)
for i in 1:length(clist)
     rename!(df_cfm, names(df_cfm)[i] => Symbol(clist[i]));
end

last(df_cfm, 6)

# Converting the Dates:
colstring = Array{Union{Nothing, String}}(nothing, length(names(df_cfm_src))-4,2)
coldates  = Array{Date,length(names(df_cfm_src))-4}
for i in 1:length(names(df_cfm_src))-4 colstring[i,1] = String(names(df_cfm_src)[i+4])*"20" end #  adding 20 to get 2020
colstring2=colstring
for i=1:length(colstring[:,1]) colstring2[i,1]=replace(colstring[i,1], "/" => "-") end
coldates=Date.(colstring2[:,1], Dates.DateFormat("mm-dd-yyyy"))
df_00 = DataFrame([coldates])
rename!(df_00, :x1 => :Date);

# Joining Dates and Data:
df_cfm = hcat(df_00, df_cfm);
# The Epidemic takes hold for these countries from late March
# To show the last week we use
last(df_cfm, 7)

CSV.write("dfs/df_cfm.csv", df_cfm)

plot(df_cfm[:,1], cfm,
     label = clist,
     legend=:outertopright, xrotation=45,
     title = "Corona 2020\nConfirmed Cases - Selected Countries")

savefig("figs/cfm_01.png")

th = 100;
# Finding first occurrence of more than 'th' confirmed cases
firsts = zeros(Int, size(df_cfm)[2]-1);
for j in 1:size(df_cfm)[2]-1 firsts[j] = findfirst(x->x>th, cfm[:,j]) end
firsts

# Extracting the data that is above the treshold:
mtx = Array{Union{Missing, Float64}}(missing,size(df_cfm)[1]-minimum(firsts)+1, size(df_cfm)[2]-1)
for j in 1:size(mtx)[2]
    for i in 1:size(cfm)[1]-firsts[j]+1   mtx[i,j] = cfm[firsts[j]-1+i,j] end
end
df_cfm_th = DataFrame(mtx)
for i in 1:length(clist)
    rename!(df_cfm_th, names(df_cfm_th)[i] => Symbol(clist[i]));
end
size(df_cfm_th)

first(df_cfm_th,7)

last(df_cfm_th,7)

CSV.write("dfs/df_cfm_th.csv", df_cfm_th)

# PLOTTING:
    updated = @sprintf "%10s" df_cfm.Date[end]

plot(mtx,
     label = clist,
     legend=:outertopright,
     xlabel = "Days since 100th Case\nEnd of data: "*updated, xguidefontsize=8,
     title = "Corona 2020\nConfirmed Cases - Selected Countries")

savefig("figs/cfm_02.png")

# Differences
cfm_th = Matrix(df_cfm_th);
dcfm_th = diff(cfm_th, dims=1)

# 1: Check if zero:
for i = 1:size(dcfm_th)[1]
    for j = 1:size(dcfm_th)[2]-1
        if !ismissing(dcfm_th[i,j])
            if dcfm_th[i,j] == 0.0 && !ismissing(dcfm_th[i+1,j])
                dcfm_th[i,j]   = 0.5*dcfm_th[i+1,j]÷1
                dcfm_th[i+1,j] = dcfm_th[i,j]+2*(dcfm_th[i+1,j]%1)
            end
         end
    end
end
dcfm_th

# 2: Check if negativ: (Norway!)
for i = 1:size(dcfm_th)[1]
    for j = 1:size(dcfm_th)[2]-1
        if !ismissing(dcfm_th[i,j])
            if dcfm_th[i,j] < 0.0 && !ismissing(dcfm_th[i+1,j])
                dcfm_th[i,j]   = 0.5*(dcfm_th[i,j] + dcfm_th[i+1,j])÷1
                dcfm_th[i+1,j] = dcfm_th[i,j]
            end
         end
    end
end

df_dcfm_th = DataFrame(dcfm_th, names(df_cfm)[2:end])
CSV.write("dfs/df_dcfm_th.csv", df_dcfm_th)

# Moving average:
n=5
dcfm_thmv= Array{Union{Missing, Float64}}(missing,size(dcfm_th)[1]-(n-1), size(dcfm_th)[2])
vsm(vs,n) = [sum(@view vs[i:(i+n-1)])/n for i in 1:(length(vs)-(n-1))]
for j in 1:size(dcfm_th)[2]
    dcfm_thmv[:,j] = vsm(dcfm_th[:,j],5).÷1
end
df_dcfm_thmv = DataFrame(dcfm_thmv, names(df_cfm)[2:end])
CSV.write("dfs/df_dcfm_thmv.csv", df_dcfm_thmv)

# PLOTTING:
    updated = @sprintf "%10s" df_cfm.Date[end]

plot(dcfm_thmv,
     label = clist,
     legend=:outertopright,
     xlabel = "Days since 100th Case\nEnd of data: "*updated, xguidefontsize=8,
     title = "Corona 2020\nDaily New Confirmed Cases - Selected Countries")

savefig("figs/cfm_03.png")

function dgamma(mid, std) # gt Discretized Gamma distribution
    midv = mid  #3.0
    stdv = std  #1.5
    tmax = ceil(midv+10*stdv)
    a = midv*midv/(stdv*stdv)
    s = stdv*stdv / midv
    gd=Gamma(a,s)
    tscale = range(0.5, stop=13.5, step=1)
    gtg = zeros(size(tscale)[1])
    for i in 2:size(gtg)[1] gtg[i] = cdf(gd,tscale[i])-cdf(gd,tscale[i-1]) end
    return gtg
end # function dgamma

function dweibull(mid, std) # gt Discretized Weibull distribution
    midv = mid  #8.4
    stdv = std  #3.8
    tmax = ceil(midv+10*stdv)
    #a = midv*midv/(stdv*stdv)
    a = 1.2785 * (midv/stdv) - 0.5004
    # s = stdv*stdv / midv
    s = midv / gamma(1 + ( 1/midv))
    gw=Weibull(a,s)
    tscale = range(0.5, stop=13.5, step=1)
    gtw = zeros(size(tscale)[1])
    for i in 2:size(gtw)[1] gtw[i] = cdf(gw,tscale[i])-cdf(gw,tscale[i-1]) end
    return gtw
end # function dweibull

function R0( gt, ep)
    Tmax = length(ep)
    nsim = 100
    q = [0.025, 0.975]

    gtpad = gt
    if length(gtpad) < Tmax gtpad = vcat(gtpad, zeros(Tmax-length(gtpad))) end

    P = zeros(Tmax, Tmax)
    p = zeros(Tmax, Tmax)
    ms  = Matrix{Float64}[]
    ms0 = zeros(Tmax, nsim)
    for s in 1:Tmax push!(ms, ms0) end

    # Loop on epidemic duration:
    for s in 2:Tmax
        ms[s] = ms0
        if ep[s] > 0
            gttmp=convert(Array, view(gtpad, s:-1:1))
            wcs = ep[1:s].*gttmp[1:s]
            wcs = wcs./sum(wcs)
            pcs = wcs .* (ep[s]./ep[1:s])
            for i in 1:(s-1) if ep[i]==0 pcs[i] = 0.0 end end
            if ep[s] == 1 pcs[s] = 0.0 end
            P[1:s,s] = pcs
            p[1:s,s] = wcs
            nsims=zeros(s, nsim)

            for i in 1:nsim  nsims[1:s,i] = rand(Multinomial(convert(Int64,ep[s]), p[1:s,s])) end
            ms[s][1:s,1:nsim] = ms[s-1][1:s,1:nsim] + nsims
        else
            P[1:s,s] .= 0.0
            p[1:s,s] .= 0.0
            ms[s][1:s,1:nsim] = ms[s-1][1:s,1:nsim]
        end
    end

        #We now have enough data to compute R (from infection network P)
    #along with its 5% and 95% quantiles (from multiple simulations and p)
    # R.WT<-apply(P,1,sum)  # Hvorfor gir dim=2 korrekt resultat?
    rwt = mapslices(sum, P, dims = [2])
    rcr = rwt./(cumsum(gtpad[1:Tmax]))[Tmax:-1:1]
    if isnan(rcr[length(ep)])  rcr[length(ep)] = 0.0 end

    #Simulated incidence at each time unit is the sum of all cases,
    #stored in the last element of multinom.sim list
    tnf   = ms[length(ep)]
    rsm   = tnf./ep            # Belgium s06 ->NaN divisjon på 0 -> Glatte input data nødvendig.
    rsmcr = rsm./(cumsum(gtpad[1:Tmax]))[Tmax:-1:1]
    rsmcr[length(ep),:] = zeros(nsim)

    # Quantiles:
    qsm    = zeros(Tmax,2)
    qsmcr  = zeros(Tmax,2)
    for s in 1:Tmax
        if (ep[s] == 0)
          rwt[s] <- 0
          rsm[s] <- 0
          rcr[s] <- 0
          rsmcr[s] <- 0
        end
        qsm[s,:]   = quantile(rsm[s,:], q)
        qsmcr[s,:] = quantile(rsmcr[s,:], q)
    end

   # Confidence Interval:
    # multinomial simulations at each time step with the expected value of R.
    correct=true
    cint = zeros(Tmax,2) # colnames?
    if correct == true
        R = rcr
        cint = qsmcr
    else
        R =rwt
        cint = qsm
    end

    # Theoretical Predictions:
    pred = vcat(ep, gt)
    for i in 2:length(ep)+length(gt) pred[i]=0.0 end
    for s in 1:Tmax
        pred[s:s+length(gt)-1] = pred[s:s+length(gt)-1] + R[s] .* ep[s] * gt
    end
    pred = pred[1:Tmax]

return R, pred, cint
end # function R0

# MAIN:
df       = df_dcfm_thmv

df_R0    = DataFrame(Array{Union{Missing, Float64}}(missing,size(df)), names(df))
df_Pred  = DataFrame(Array{Union{Missing, Float64}}(missing,size(df)), names(df))
df_Clow  = DataFrame(Array{Union{Missing, Float64}}(missing,size(df)), names(df))
df_Chgh  = DataFrame(Array{Union{Missing, Float64}}(missing,size(df)), names(df))
mtxCint =Array{Tuple{Union{Missing,Float64},Union{Missing,Float64}},2}(undef, size(df))
for i in 1:size(mtxCint)[1] for j in 1:size(mtxCint)[2] mtxCint[i,j] = (missing, missing) end end
df_Cint = DataFrame(mtxCint, names(df))

#gtx = dgamma(8.4, 3.8) # Default Discrete Gamma Distribution
gtx = dweibull(8.4, 3.8) # Default Discrete Weibull Distribution
# COUNTRY LOOP:
for j in 1:size(df)[2]
    cname=names(df)[j]
    epx = dropmissing(df, cname)[:,j]
    println(j," ", cname," ", sum(epx), " ",   length(epx))
    R, pred, cint = R0(gtx,epx)
    for i in 1:length(R)-1      df_R0[i,j] = R[i] end
    for i in 1:length(pred)   df_Pred[i,j] = pred[i] end
    for i in 1:size(cint)[1]-1  df_Clow[i,j] = cint[i,1] end
    for i in 1:size(cint)[1]-1  df_Chgh[i,j] = cint[i,2] end
    for i in 1:size(cint)[1]-1  df_Cint[i,j] = (cint[i,1], cint[i,2]) end
end

CSV.write("dfs/df_R0.csv", df_R0)

R = Matrix(df_R0)

plot(R,
     label = clist,
     legend=:outertopright,
     xlabel = "Days since 100th Case\nEnd of data: "*updated, xguidefontsize=8,
     title = "Corona 2020\nReproduction Number R - Selected Countries")

savefig("figs/cfm_04.png")

R = Matrix(df_R0)

plot([20:size(R)[1]],R[20:end,:],
     label = clist,
     legend=:outertopright,
     xlabel = "Days since 100th Case\nEnd of data: "*updated, xguidefontsize=8,
     title = "Corona 2020\nReproduction Number R - Selected Countries")

savefig("figs/cfm_05.png")

Rlast = [findfirst( x-> ismissing(x), R[:,i]) for i in 1:size(R)[2]]
for i in 1:length(Rlast)
    if Rlast[i] == nothing
       Rlast[i] = size(R,1)
    else
       Rlast[i]= Rlast[i] - 1
    end
end

labels1 = String[]
for i in 1:length(names(df_R0)) push!(labels1, String(names(df_R0)[i])) end
labels = reshape(labels1,1,:)

plot([20:size(R)[1]],R[20:end,:],
     label = clist,
     legend=:outertopright,
     xlabel = "Days since 100th Case\nEnd of data: "*updated, xguidefontsize=8,
     annotations=([(Rlast[i], R[Rlast[i],i],text(labels[i][1:2], :left, 8)) for i in 1:length(Rlast)]),
     title = "Corona 2020\nReproduction Number R - Selected Countries")

savefig("figs/cfm_06.png")

max_stars_repo_name: dyrehaugen/jdt

max_issues_repo_name: dyrehaugen/jdt

max_forks_repo_name: dyrehaugen/jdt


text:
y = 7
x = 4
text = typeof(y)
println(text)
z = 3x + y
println(z)
println(0 < x < 3) # false

# One line comment

#=
Mulitple line comment
=#

text = isa(x, Int64)
println(text)

# casting
text = trunc(Int64, 3.14159) # small case type name as a function
println(text)
text = trunc(Int64, -3.14159) # small case type name as a function
println(text)

# display an object
#display(z)

arr = [100,10,0]
println("typeof(arr): $(typeof(arr))")

arr3 = Float64[]
println("typeof(arr3): $(typeof(arr3))")

function f()
  # type can only be specified in function
  x::Int64 = 50
  println("typeof(x): $(typeof(x))")
  println("x = $x")

  s::String = "A string"
  println("typeof(s): $(typeof(s))")
end

f()


text:
# Jonas August, 2021

#=
This problem was asked by Facebook.

Given a string of round, curly, and square open and closing brackets, return whether the brackets are balanced (well-formed).

For example, given the string "([])[]({})", you should return true.

Given the string "([)]" or "((()", you should return false.
=#

#=

Observations:

- Closing bracket follows matching opening bracket.

- Brackets come in matched pairs.

- Cannot close an outer pair before inner is closed.

Recursive definition of well-formed string:

good_string
= no brackets OR
 good_string, opening bracket, good_string, matched closing bracket, good_string

=#

closer = Dict("(" => r"\)", "{" => r"\}", "[" => r"\]")

function isgood(s::String)
  mnone = match(r"^[^(){}[\]]*$", s)
  if !isnothing(mnone)
    return true
  end
  mleft = match(r"[({[]", s)
  if isnothing(mleft)
      return false
  end
  mright = match(closer[mleft.match], s)
  if isnothing(mright)
      return false
  end
  ileft = mleft.offset
  iright = mright.offset
  if ileft >= iright
      return false
  end
  sleft = s[1:ileft-1]
  smid = s[ileft+1:iright-1]
  sright = s[iright+1:end]
  isgood(sleft) && isgood(smid) && isgood(sright)
end

using Test

well_formed = ["", "ab", "()", "[]", "()[]", "{}", "([])", "([])[]", "([])[]({})"]
@testset "well formed string \"$s\"" for s in well_formed
    @test isgood(s)
end

badly_formed = ["(", "}", ")(", "([)]", "((()"]
@testset "badly formed string \"$s\"" for s in badly_formed
    @test !isgood(s)
end

max_stars_repo_path: Balanced brackets - Problem 27.jl

max_issues_repo_path: Balanced brackets - Problem 27.jl

max_forks_repo_path: Balanced brackets - Problem 27.jl


text:
using SafeBuffers
using Test

@testset "SafeBuffers.jl" begin
    @test_nowarn SafeBuffer(2, Ref{Int})
    @test_nowarn SafeBuffer(2, Vector{Int}, 2)
    @test_nowarn SafeBuffer(2, Vector{Int}, (2,))
    @test_nowarn SafeBuffer(2, Matrix{Int}, 2, 2)
    @test_nowarn SafeBuffer(2, Matrix{Int}, (2,2))

    buf = Matrix{Int}(undef, 2,2)
    @test_nowarn SafeBuffer(2, buf)

    b1 = SafeBuffer(2, Vector{Int}, 2)
    @test 5 == withbuffer(b1) do buf
        buf .= 1:2
        sum(abs2, buf)
    end

    @test_nowarn begin
        t1 = @async withbuffer(b1) do buf
            buf .= 1
            sleep(1)
            yield()
            sum(log, buf)
        end
        t2 = @async withbuffer(b1) do buf
            buf .= -1
            nothing
        end
        t3 = @async withbuffer(b1) do buf
            buf .= -1
            nothing
        end
        wait.([t1, t2, t3])
    end

end


text:
"""
creating asset tree hierarchy

"""

abstract type Asset end

abstract type Property <: Asset end
abstract type Investment <: Asset end
abstract type Cash <: Asset end

abstract type House <: Property end
abstract type Apartment <: Property end

abstract type FixedIncome <: Investment end
abstract type Equity <: Investment end

#composite type
#by default, is immutable
#can be made mutable by adding mutable in front of struct
mutable struct Stock <: Equity
    symbol::String
    name::String
end

function describe(s::Stock)
    return s.symbol * "(" * s.name * ")"
end

struct BasketOfStocks
    stocks::Vector{Stock}
    reason::String
end

abstract type Art end

struct Painting <: Art
    artist::String
    title::String
end

#demonstrating how Union works: enables multiple data types to be allowed in a field
const Thing = Union{Painting,Stock}

struct BasketOfThings
    thing::Vector{Thing}
    reason::String
end

# Display the entire type hierarchy starting from the specified `roottype`
function subtypetree(roottype, level = 1, indent = 4)
    level == 1 && println(roottype)
    for s in subtypes(roottype)
       println(join(fill(" ", level * indent)) * string(s))
       subtypetree(s, level + 1, indent)
    end
end

#using a type parameter enables flexible typing - ex: be able to hold int and fractional shares
struct StockHolding{T <: Real}
    stock::Stock
    quantity::T
end

#using two type parameters: this enforces price and marketvalue to both be of the same data type
struct StockHolding2{T <: Real, P <: AbstractFloat}
    stock::Stock
    quantity::T
    price::P
    marketvalue::P
end

#parametric abstract types
abstract type Holding{P} end

struct StockHolding3{T, P} <: Holding{P}

    stock::Stock
    quantity::T
    price::P
    marketvalue::P
   end

   struct CashHolding{P} <: Holding{P}
    currency::String
    amount::P
    marketvalue::P
   end


text:
using Documenter, MyPkgTemplates

makedocs(;
    modules=[MyPkgTemplates],
    format=:html,
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/garborg/MyPkgTemplates.jl/blob/{commit}{path}#L{line}",
    sitename="MyPkgTemplates.jl",
    authors="Sean Garborg",
    assets=[],
)

deploydocs(;
    repo="github.com/garborg/MyPkgTemplates.jl",
    target="build",
    julia="1.0",
    deps=nothing,
    make=nothing,
)


text:
"""
$(SIGNATURES)

Checks for a `config.md` file and uses it to set the global variables
referenced in `GLOBAL_VARS` it also sets the global latex commands via
`GLOBAL_LXDEFS`. If the configuration file is not found a warning is shown.
The keyword `init` is used internally to distinguish between the first call
where only structural variables are considered (e.g. controlling folder
structure).
"""
function process_config()::Nothing
    FD_ENV[:SOURCE] = "config.md"
    config_path = joinpath(FOLDER_PATH[], "config.md")
    if isfile(config_path)
        convert_md(read(config_path, String); isconfig=true)
    else
        config_warn()
    end
    return nothing
end

"""
$(SIGNATURES)

Checks for a `utils.jl` file and uses it to set global computed variables,
functions and html functions. Whatever is defined in `utils.jl` takes
precedence over things defined internally in Franklin or in the global vars;
in particular users can redefine the behaviour of `hfuns` though that's not
recommended.
"""
function process_utils()::Nothing
    FD_ENV[:SOURCE] = "utils.jl"
    utils = joinpath(FOLDER_PATH[], "utils.jl")

    if !isfile(utils)
        FD_ENV[:UTILS_COUNTER] = 0
        FD_ENV[:UTILS_HASH] = nothing
        return nothing
    end

    # has the hash changed?
    if FD_ENV[:UTILS_HASH] !== nothing
        new_hash    = hash(read(utils, String))
        has_changed = (new_hash != FD_ENV[:UTILS_HASH])
        has_changed || return nothing
        FD_ENV[:UTILS_HASH] = new_hash
    end
    FD_ENV[:UTILS_COUNTER] += 1

    # create new Utils module
    newmodule(utils_name())
    Base.include(utils_module(), utils)

    # keep track of utils names
    ns = String.(names(utils_module(), all=true))
    filter!(n -> n[1] != '#' && n ∉ ("eval", "include", "Utils"), ns)
    empty!(UTILS_NAMES)
    append!(UTILS_NAMES, ns)
    return nothing
end

"""
$(SIGNATURES)

See [`process_file_err`](@ref).
"""
function process_file(case::Symbol, fpair::Pair{String,String}, args...)::Int
    if FD_ENV[:DEBUG_MODE]
        process_file_err(case, fpair, args...)
        return 0
    end

    try
        process_file_err(case, fpair, args...)
    catch err
        rp = fpair.first
        rp = rp[end-min(20, length(rp))+1 : end]
        if !FD_ENV[:QUIET_TEST]
            FD_ENV[:SOURCE] = fpair.second
            print_warning("""
                Encountered an issue processing '$(fpair.second)' in $rp.
                Verify, then re-start the Franklin server.
                The error is displayed below:
                $err
                """)
        end
        FD_ENV[:SUPPRESS_ERR] || throw(err)
        return -1
    end
    return 0
end

"""
$(SIGNATURES)

Considers a source file which, depending on `case` could be a HTML file or a
file in Franklin-Markdown etc, located in a place described by `fpair`,
processes it by converting it and adding appropriate header and footer and
writes it to the appropriate place. It can throw an error which will be
caught in `process_file(args...)`.
"""
function process_file_err(case::Symbol, fpair::Pair{String, String},
                          head::AS="", pgfoot::AS="", foot::AS="",
                          t::Float64=0.0)::Nothing
    # depending on the file extension, either full process (.md), partial
    # process (.html) or no process (everything else)
    inp  = joinpath(fpair...)
    outp = form_output_path(fpair.first, fpair.second, case)
    if case == :md
        FD_ENV[:SOURCE] = get_rpath(inp)
        convert_and_write(fpair..., head, pgfoot, foot, outp)
    elseif case == :html
        FD_ENV[:SOURCE] = get_rpath(inp)
        set_cur_rpath(joinpath(fpair...))
        set_page_env()
        raw_html  = read(inp, String)
        # add the item *before* the conversion so that the conversion
        # can affect the page itself with {{...}}
        cond_add = globvar(:generate_sitemap) && FD_ENV[:FULL_PASS]
        cond_add && add_sitemap_item(html=true)
        proc_html = convert_html(raw_html) |> postprocess_page
        write(outp, proc_html)
    else # case in (:other, :infra)
        # there's a bunch of thing we don't want to copy over
        if startswith(inp, path(:layout))   ||
            startswith(inp, path(:literate)) ||
            endswith(inp, "config.md") ||
            endswith(inp, "search.md")
            # skip
            @goto end_copyblock
        end
        # NOTE: some processing may be added here later on (e.g. parsing of
        # CSS files). Only copy again if necessary (file is not there or
        # has changed)
        if !isfile(outp) || (mtime(outp) < t && !filecmp(inp, outp))
            cp(inp, outp, force=true)
        end
        @label end_copyblock
    end
    FD_ENV[:SOURCE] = ""
    FD_ENV[:FULL_PASS] || FD_ENV[:SILENT_MODE] || rprint("→ page updated [✓]")
    return nothing
end

"""
$(SIGNATURES)

Convenience function to replace the extension of a filename with another.
"""
change_ext(fname::AS, ext=".html")::String = splitext(fname)[1] * ext

"""
    get_rpath(fpath)

Extracts the relative file system path out of the full system path to a file
currently being processed. Does not start with a path separator.
So `[some_fs_path]/blog/page.md` --> `blog/page.md` (keeps the extension).
"""
function get_rpath(fpath::String)
    root = path(:folder)
    return fpath[lastindex(root)+length(PATH_SEP)+1:end]
end

"""
    set_cur_rpath(fpath)

Takes the path to the current file and sets the `fd_rpath` local page variable
as well as the `FD_ENV[:CUR_PATH]` variable (used for conditional blocks
depending on URL for instance).
"""
function set_cur_rpath(fpath::String; isrelative::Bool=false)
    if isrelative
        rpath = fpath
    else
        rpath = get_rpath(fpath)
    end
    FD_ENV[:CUR_PATH] = rpath
    set_var!(LOCAL_VARS, "fd_rpath", rpath)
    return nothing
end


text:
include("2048_fn.jl")
using BenchmarkTools
@benchmark simulate_game()


text:
  line = "2005 1     0 5500   540 120031   1DBBX       0310          102452  60  76        74    71    A               165 14255992114 5 0                   111FF1377AAA11A118AA1      594    1                 4                                                       9863"
  val = imma(rstrip(line))
  @test val[:SST] == 7.1000000000000005
  @test val[:YR] == 2005.0


text:
@testset "Basicfunctions Module Tests" begin

  param_file = "params.yml"
  model = Science_ABM.initialize_model(param_file)
  problem = Science_ABM.choose_problem(model[1], model)

  @testset "choose_problem Tests" begin
    @test isdefined(Science_ABM, :choose_problem) == true
    @test isdefined(problem, :complexity)
    @inferred Science_ABM.choose_problem(model[1], model)
    @test length(fieldnames(typeof(problem))) ==  5
    @test problem.success_probability < 1
    @test problem.success_probability > 0
  end

  @testset "pushproblem! Tests" begin
    @test isdefined(Science_ABM, :pushproblem!) == true
    @test typeof(model[1].publication_success[1]) == Bool
  end
end


text:
# Functions for querying network properties.

######### Accessors: #########

"""
    species(network)
Given an `AbstractReactionNetwork`, return a vector of species symbols.
"""
function species(network)
    network.syms
end

"""
    params(network)
Given an `AbstractReactionNetwork`, return a vector of parameter symbols.
"""
function params(network)
    network.params
end

"""
    speciesmap(network)

Given an `AbstractReactionNetwork`, return a Dictionary mapping from species
symbol to species index.
"""
function speciesmap(network)
    network.syms_to_ints
end

"""
    paramsmap(network)

Given an `AbstractReactionNetwork`, return a Dictionary mapping from parameter
symbol to parameter index.
"""
function paramsmap(network)
    network.params_to_ints
end

"""
    numspecies(network)

Return the number of species within the given `AbstractReactionNetwork`.
"""
function numspecies(network)
    length(speciesmap(network))
end

"""
    numreactions(network)

Return the number of reactions within the given `AbstractReactionNetwork`.
"""
function numreactions(network)
    length(network.reactions)
end

"""
    numparams(network)

Return the number of parameters within the given `AbstractReactionNetwork`.
"""
function numparams(network)
    length(paramsmap(network))
end

######### Generated Functions: #########

"""
    oderhsfun(network)

Given an `AbstractReactionNetwork`, return a function, `f!(du,u,p,t)`, that
evaluates the current value of the ODE model derivative functions, ``du/dt = f(u,t)``,
within `du`.

*Note,* for a network generated with the `@min_reaction_network` macro `addodes!`
must be called first.
"""
function oderhsfun(network)
    isnothing(network.f) && error("Error, call addodes! first.")
    network.f
end

"""
    jacfun(network)

Given an `AbstractReactionNetwork`, return a function, `jac!(J,u,p,t)`, that
evaluates the current Jacobian matrix, `J`, of the ODE model, ``du/dt = f(u,t)``.
The Jacobian matrix has entries

``J_{i,j} = \\partial f_i(u,t) / \\partial u_j``.

*Note,* for a network generated with the `@min_reaction_network` macro `addodes!`
must be called first.
"""
function jacfun(network)
    isnothing(network.jac) && error("Error, call addodes! first.")
    network.jac
end

"""
    paramjacfun(network)

Given an `AbstractReactionNetwork`, return a function, `pjac(pJ,u,p,t)`, that
evaluates the current parameter Jacobian matrix, `pJ`, of the ODE model, ``du/dt = f(u,t)``.
The parameter Jacobian matrix has entries

``pJ_{i,j} = \\partial f_i(u,t) / \\partial p_j``.

*Note,* for a network generated with the `@min_reaction_network` macro `addodes!`
must be called first.
"""
function paramjacfun(network)
    isnothing(network.paramjac) && error("Error, call addodes! first.")
    network.paramjac
end

"""
    odefun(network)

Given an `AbstractReactionNetwork`, return a `DiffEqBase.ODEFunction` encoding
an ODE model for the reaction network.

*Note,* for a network generated with the `@min_reaction_network` macro `addodes!`
must be called first.
"""
function odefun(network)
    isnothing(network.odefun) && error("Error, call addodes! first.")
    network.odefun
end

"""
    noisefun(network)

Given an `AbstractReactionNetwork`, return a function, `g(η,u,p,t)`, that
evaluates the current noise coefficients for each reaction in the Chemical
Langevin Equation representation within `η`.

*Note,* for a network generated with the `@min_reaction_network` macro `addsdes!`
must be called first.
"""
function noisefun(network)
    isnothing(network.g) && error("Error, call addsdes! first.")
    network.g
end

"""
    sdefun(network)

Given an `AbstractReactionNetwork`, return a `DiffEqBase.SDEFunction` encoding
a Chemical Langevin Equation SDE model for the reaction network.

*Note,* for a network generated with the `@min_reaction_network` macro `addsdes!`
must be called first.
"""
function sdefun(network)
    isnothing(network.sdefun) && error("Error, call addsdes! first.")
    network.sdefun
end

"""
    jumps(network)

Given an `AbstractReactionNetwork`, return a tuple of `AbstractJumps` encoding
a stochastical chemical kinetics representation for the reaction network.

*Note,* for a network generated with the `@min_reaction_network` macro `addjumps!`
must be called first.
"""
function jumps(network)
    isnothing(network.jumps) && error("Error, call addjumps! first.")
    network.jumps
end

"""
    regularjumps(network)

Given an `AbstractReactionNetwork`, return a `RegularJump` encoding a
stochastical chemical kinetics representation of the reaction network for use in
``\\tau``-leaping approximations.

*Note,* for a network generated with the `@min_reaction_network` macro `addjumps!`
must be called first.
"""
function regularjumps(network)
    isnothing(network.regular_jumps) && error("Error, call addjumps! first.")
    network.regular_jumps
end

######### Generated Expressions: #########

"""
    odeexprs(network)

Given an `AbstractReactionNetwork`, return a vector of the ODE expressions.

*Note,* for a network generated with the `@min_reaction_network` macro `addodes!`
must be called first.
"""
function odeexprs(network)
    isnothing(network.f_func) && error("Error, call addodes! first.")
    network.f_func
end

"""
    jacobianexprs(network)

Given an `AbstractReactionNetwork`, return a matrix with the ODE Jacobian
expressions.

*Note,* for a network generated with the `@min_reaction_network` macro `addodes!`
must be called first.
"""
function jacobianexprs(network)
    isnothing(network.symjac) && error("Error, call addodes! first.")
    network.symjac
end

"""
    noiseexprs(network)

Given an `AbstractReactionNetwork`, return a vector of the SDE noise expressions
for each reaction.

*Note,* for a network generated with the `@min_reaction_network` macro `addsdes!`
must be called first.
"""
function noiseexprs(network)
    isnothing(network.g_func) && error("Error, call addsdes! first.")
    network.g_func
end

"""
    jumpexprs(network)

Given an `AbstractReactionNetwork`, return a tuple of the jump rates and affects
expressions.

*Note,* for a network generated with the `@min_reaction_network` macro `addjumps!`
must be called first.
"""
function jumpexprs(network)
    (isnothing(network.jump_rate_expr) || isnothing(network.jump_affect_expr)) && error("Error, call addjumps! first.")
    network.jump_rate_expr,network.jump_affect_expr
end

"""
    rateexpr(network, rxidx)

Given an `AbstractReactionNetwork`, return the reaction rate expression for the
reaction with index `rxidx`. Note, for a reaction defined by

`k*X*Y, X+Z --> 2X + Y`

the expression that is returned will be `:(k*X*Y)`, while the *rate law* used in
ODEs and SDEs would be `k*X^2*Y*Z`.
"""
function rateexpr(network, rxidx)
    network.reactions[rxidx].rate_org
end

"""
    oderatelawexpr(network, rxidx)

Given an `AbstractReactionNetwork`, return the reaction rate law expression used
in generated ODEs for the reaction with index `rxidx`. Note, for a reaction
defined by

`k*X*Y, X+Z --> 2X + Y`

the expression that is returned will be `:(k*X^2*Y*Z)`. For a reaction of the
form

`k, 2X+3Y --> Z`

the expression that is returned will be `:(k * (X^2/2) * (Y^3/6))`.
"""
function oderatelawexpr(network, rxidx)
    network.reactions[rxidx].rate_DE
end

"""
    ssaratelawexpr(network, rxidx)

Given an `AbstractReactionNetwork`, return the reaction rate law expression used
in generated stochastic chemical kinetic model SSAs for the reaction with index
`rxidx`. Note, for a reaction defined by

`k*X*Y, X+Z --> 2X + Y`

the expression that is returned will be `:(k*X^2*Y*Z)`. For a reaction of the
form

`k, 2X+3Y --> Z`

the expression that is returned will be `:(k * binomial(X,2) * binomial(Y,3))`.
"""
function ssaratelawexpr(network, rxidx)
    network.reactions[rxidx].rate_SSA
end

######### Reaction Properties: #########

function substratestoich(rs::ReactionStruct, specmap)
    sort!( [specmap[s.reactant] => s.stoichiometry for s in rs.substrates] )
end

"""
    substratestoich(network, rxidx)

Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
vector of pairs, mapping ids of species that serve as substrates in the reaction
to the corresponding stoichiometric coefficient as a substrate.

Allocates a new vector to store the pairs.
"""
function substratestoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
    substratestoich(rn.reactions[rxidx], speciesmap(rn))
end

"""
    substratesymstoich(network, rxidx)

Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
Vector of `ReactantStruct`s, mapping the symbols of species that serve as
substrates in the reaction to the corresponding stoichiometric coefficient as a
substrate.

Non-allocating, returns underlying field within the reaction_network.
"""
function substratesymstoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
    rn.reactions[rxidx].substrates
end

function productstoich(rs::ReactionStruct, specmap)
    sort( [specmap[p.reactant] => p.stoichiometry for p in rs.products] )
end

"""
    productstoich(network, rxidx)

Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
vector of pairs, mapping ids of species that are products in the reaction to the
corresponding stoichiometric coefficient as a product.

Allocates a new vector to store the pairs.
"""
function productstoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
    productstoich(rn.reactions[rxidx], speciesmap(rn))
end

"""
    productsymstoich(network, rxidx)

Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
Vector of `ReactantStruct`s, mapping the symbols of species that are products in
the reaction to the corresponding stoichiometric coefficient as a product.

Non-allocating, returns underlying field within the reaction_network.
"""
function productsymstoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
    rn.reactions[rxidx].products
end

function netstoich(rs::ReactionStruct, specmap)
    sort!([specmap[ns.reactant] => ns.stoichiometry for ns in rs.netstoich])
end

"""
    netstoich(network, rxidx)

Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
vector of pairs, mapping ids of species that change numbers due to the reaction
to the net stoichiometric coefficient of the species (i.e. net change in the
species due to the reaction).

Allocates a new vector to store the pairs.
"""
function netstoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
    netstoich(rn.reactions[rxidx], speciesmap(rn))
end

"""
    netsymstoich(network, rxidx)

Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
Vector of `ReactantStruct`s, mapping the symbols of species that change numbers
due to the reaction to the net stoichiometric coefficient of the species (i.e.
net change in the species due to the reaction).

Non-allocating, returns underlying field within the reaction_network.
"""
function netsymstoich(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
    rn.reactions[rxidx].netstoich
end

"""
    ismassaction(network, rxidx)

Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
boolean indicating whether the given reaction is of mass action form. For
example, the reaction

`2*k, 2X + 3Y --> 5Z + W`

would return true, while reactions with state-dependent rates like

`k*X, X + Y --> Z`

would return false.

Non-allocating, returns underlying field within the reaction_network.
"""
function ismassaction(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
    rn.reactions[rxidx].is_pure_mass_action
end

"""
    dependents(network, rxidx)

Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
vector of symbols of species the *reaction rate law* depends on. i.e. for

`k*W, 2X + 3Y --> 5Z + W`

the returned vector would be `[:W,:X,:Y]`.

Non-allocating, returns underlying field within the reaction_network.
"""
function dependents(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
    rn.reactions[rxidx].dependants
end

"""
    dependants(network, rxidx)

See documentation for [`dependents(network, rxidx)`](@ref).
"""
function dependants(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
    dependents(rn, rxidx)
end

"""
    substrates(network, rxidx)

Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
vector of symbols of species that correspond to substrates in the reaction.
i.e. for

`k*W, X + 3Y --> X + W`

the returned vector would be `[:X,:Y]`.

Allocates a new vector to store the symbols.
"""
function substrates(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
    [sub.reactant for sub in rn.reactions[rxidx].substrates]
end

"""
    products(network, rxidx)

Given an `AbstractReactionNetwork` and a reaction index, `rxidx`, return a
vector of symbols of species that correspond to products in the reaction.
i.e. for

`k*W, X + 3Y --> X + W`

the returned vector would be `[:X,:W]`.

Allocates a new vector to store the symbols.
"""
function products(rn::DiffEqBase.AbstractReactionNetwork, rxidx)
    [prod.reactant for prod in rn.reactions[rxidx].products]
end

######### Network Properties: #########

"""
    rxtospecies_depgraph(network)

Given an `AbstractReactionNetwork`, returns a `Vector{Vector{Int}}` mapping each
reaction index to the indices of species that depend on it.
"""
function rxtospecies_depgraph(network)
    specmap = speciesmap(network)

    # map from a reaction to vector of species that depend on it
    [sort!( [ns.first for ns in netstoich(rx, specmap)] ) for rx in network.reactions]
end

"""
    speciestorx_depgraph(network)

Given an `AbstractReactionNetwork`, returns a `Vector{Vector{Int}}` mapping each
species index to the indices of reactions that depend on it.
"""
function speciestorx_depgraph(network)
    numrxs  = numreactions(network)
    specmap = speciesmap(network)

    # map from a species to reactions that depend on it
    spectorxs = [Vector{Int}() for i=1:numspecies(network)]
    for rx in 1:numrxs
        for specsym in dependents(network, rx)
            push!(spectorxs[specmap[specsym]], rx)
        end
    end
    foreach(stor -> sort!(stor), spectorxs)

    spectorxs
end

"""
    rxtorx_depgraph(network)

Given an `AbstractReactionNetwork`, returns a `Vector{Vector{Int}}` mapping each
reaction index to the indices of reactions that depend on it.
"""
function rxtorx_depgraph(network, sptorxs=speciestorx_depgraph(network))
    numrxs = numreactions(network)
    dep_sets = [SortedSet{Int}() for n = 1:numrxs]

    # dg as vector of sets
    for rx in 1:numrxs
        net_stoich = netstoich(network, rx)

        for (spec,stoich) in net_stoich
            for dependent_rx in sptorxs[spec]
                push!(dep_sets[rx], dependent_rx)
            end
        end

        # every reaction should depend on itself
        push!(dep_sets[rx], rx)
    end

    # dg as vector of vectors
    dep_graph = Vector{Vector{Int}}(undef, numrxs)
    for rx in 1:numrxs
        dep_graph[rx] = [dep for dep in dep_sets[rx]]
    end

    dep_graph
end


text:
using DoubleExponentialFormulas
using LinearAlgebra: norm
using Test

quadde32 = QuadDE(Float32)
quadde64 = QuadDE(Float64)
quaddeBF = QuadDE(BigFloat)

# Test integral interval [-1, 1] with QuadDE
let
    f(x::AbstractFloat) = 2/(1 + x^2)
    expect = π

    I, E = quadde32(f, -1, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, -1, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -1, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)

    # Split into multiple integral intevals
    I, E = quadde32(f, -1, 0, 1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, -1, 0, 1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -1, 0, 1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test integral interval [a, b] with QuadDE
# a, b are arbitrary finite float numbers
let
    f(x::AbstractFloat) = 1/(1 + (x/2)^2)
    expect = π

    I, E = quadde32(f, -2, 2)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, -2, 2)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -2, 2, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)

    # Split into multiple integral intevals
    I, E = quadde32(f, -2, 0, 2)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, -2, 0, 2)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -2, 0, 2, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test integral interval [0, ∞] with QuadDE
let
    f(x::AbstractFloat) = exp(-x)
    expect = 1

    I, E = quadde32(f, 0, Inf)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, Inf, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)

    # Split into multiple integral intevals
    I, E = quadde32(f, 0, 1, Inf)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 0, 1, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, 1, Inf, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test integral interval [a, ∞] with QuadDE
let
    f(x::AbstractFloat) = exp(-(x - 1))
    expect = 1

    I, E = quadde32(f, 1, Inf)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 1, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 1, Inf, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)

    # Split into multiple integral intevals
    I, E = quadde32(f, 1, 2, Inf)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, 1, 2, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 1, 2, Inf, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test integral interval [-∞, 0] with QuadDE
let
    f(x::AbstractFloat) = exp(x)
    expect = 1

    I, E = quadde32(f, -Inf, 0)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, -Inf, 0)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -Inf, 0, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)

    # Split into multiple integral intevals
    I, E = quadde32(f, -Inf, -1, 0)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, -Inf, -1, 0)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -Inf, -1, 0, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test integral interval [-∞, b] with QuadDE
let
    f(x::AbstractFloat) = exp(x + 1)
    expect = 1

    I, E = quadde32(f, -Inf, -1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, -Inf, -1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -Inf, -1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)

    # Split into multiple integral intevals
    I, E = quadde32(f, -Inf, -2, -1)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, -Inf, -2, -1)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -Inf, -2, -1, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test integral interval [-∞, ∞] with QuadDE
let
    f(x::AbstractFloat) = exp(-x^2)
    expect = sqrt(BigFloat(π))

    I, E = quadde32(f, -Inf, Inf)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, -Inf, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -Inf, Inf, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)

    # Split into multiple integral intevals
    I, E = quadde32(f, -Inf, 0, Inf)
    @test I isa Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    I, E = quadde64(f, -Inf, 0, Inf)
    @test I isa Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(typeof(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -Inf, 0, Inf, rtol=rtol)
    @test I isa BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end

# Test non-scalar output with QuadDE
let
    f(x::AbstractFloat) = [1/(1 + x^2), 2/(1 + x^2)]
    expect = [BigFloat(π)/2, BigFloat(π)]

    I, E = quadde32(f, -1, 1)
    @test eltype(I) == Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadde64(f, -1, 1)
    @test eltype(I) == Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -1, 1, rtol=rtol)
    @test eltype(I) == BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)

    I, E = quadde32(f, -1, 0, 1)
    @test eltype(I) == Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadde64(f, -1, 0, 1)
    @test eltype(I) == Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -1, 0, 1, rtol=rtol)
    @test eltype(I) == BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)

    f(x::AbstractFloat) = [exp(-x), 2exp(-x)]
    expect = [1, 2]

    I, E = quadde32(f, 0, Inf)
    @test eltype(I) == Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadde64(f, 0, Inf)
    @test eltype(I) == Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, 0, Inf, rtol=rtol)
    @test eltype(I) == BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)

    f(x::AbstractFloat) = [exp(-x^2), 2*exp(-x^2)]
    expect = [sqrt(BigFloat(π)), 2*sqrt(BigFloat(π))]

    I, E = quadde32(f, -Inf, Inf)
    @test eltype(I) == Float32
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    I, E = quadde64(f, -Inf, Inf)
    @test eltype(I) == Float64
    @test I ≈ expect
    @test E ≤ sqrt(eps(eltype(I)))*norm(I)

    rtol = 1e-30
    I, E = quaddeBF(f, -Inf, Inf, rtol=rtol)
    @test eltype(I) == BigFloat
    @test isapprox(I, expect, rtol=10rtol)
    @test E ≤ rtol*norm(I)
end