text:
# ==============================================================================
# arrondi la solution correspondant au generateur (pas d'historique donc)
# version avec cone inferieur seulement
function roundingSolution!(vg::Vector{tGenerateur}, k::Int64, c1::Array{Int,1}, c2::Array{Int,1}, d::tListDisplay)

    nbvar = length(vg[k].sInt.x)
    nbgen = size(vg,1)

    # --------------------------------------------------------------------------
    # identifie les variables fractionnaires et marquage par valeur sentinelle -1
    for i in 1:nbvar
        if  isapprox(vg[k].sPrj.x[i] , 0.0, atol=1e-3)
            vg[k].sInt.x[i] = 0
        elseif isapprox(vg[k].sPrj.x[i] , 1.0, atol=1e-3)
            vg[k].sInt.x[i] = 1
        else
            vg[k].sInt.x[i] = -1
        end
    end

    # --------------------------------------------------------------------------
    # Selectionne les points pour le cone pointe sur le generateur k (pCour) et ouvert vers Y
    pPrec, pCour, pSuiv = selectionPoints(vg, k)

    # --------------------------------------------------------------------------
    # elabore le point pointeur du cone ouvert vers L
    pC = elaborePointConeOuvertversL(vg, k, pPrec, pSuiv)

    # --------------------------------------------------------------------------
    # Arrondi les valeurs non-entieres d'une solution fractionnaire
    # 1) applique un arrondi qui fixe 1 variable (la premiere venue) par iteration

    vg[k].sInt.y[1] = 0
    vg[k].sInt.y[2] = 0

    # calcule en differentiel par rapport a l'image de la solution fractionnaire
    z1 = vg[k].sPrj.y[1]
    z2 = vg[k].sPrj.y[2]
    pM = tPoint( z1 , z2 )

    print("Depart a arrondir : ")
    @show pPrec, pCour, pSuiv, pM
    if length(vg[k].sPrj.x) ≤ 20 @show vg[k].sPrj.x end
    if length(vg[k].sInt.x) ≤ 20 @show vg[k].sInt.x end

    nbVarNonEntiere = 0
    for i in 1:nbvar
        if vg[k].sInt.x[i] == -1

            # la variable i est non entiere
            nbVarNonEntiere += 1

            # defalque la contribution fractionnaire de la variable aux objectifs => pose x[i]=0
            pM = tPoint( z1 - vg[k].sPrj.x[i] * c1[i] , z2 - vg[k].sPrj.x[i] * c2[i])

            if inSector(pM, pCour, pSuiv, pPrec)
                # le point pM obtenu est hors du cone => valide x[i]=1
                vg[k].sInt.x[i] = 1
                z1 = z1 - (vg[k].sPrj.x[i] * c1[i]) + c1[i]
                z2 = z2 - (vg[k].sPrj.x[i] * c2[i]) + c2[i]
                @printf("  x[%2d]=1 |  : [ %12.5f , %12.5f ] \n",i, z1, z2)
            else
                # le point pM obtenu est dans le cone => valide x[i]=0
                vg[k].sInt.x[i] = 0
                z1 = z1 - (vg[k].sPrj.x[i] * c1[i])
                z2 = z2 - (vg[k].sPrj.x[i] * c2[i])
                @printf("  x[%2d]=0 |  : [ %12.5f , %12.5f ] \n",i, z1, z2)
            end
        end

        # calcule la performance de la solution entiere
        vg[k].sInt.y[1] += vg[k].sInt.x[i] * c1[i]
        vg[k].sInt.y[2] += vg[k].sInt.x[i] * c2[i]
    end

    if length(vg[k].sInt.x) ≤ 20 @show vg[k].sInt.x end
    @printf("→ #round : %4d → [ %5d , %5d ] ", nbVarNonEntiere, vg[k].sInt.y[1], vg[k].sInt.y[2])

    # archive le point obtenu pour les besoins d'affichage
    if generateurVisualise == -1
        # archivage pour tous les generateurs
        push!(d.XInt,vg[k].sInt.y[1])
        push!(d.YInt,vg[k].sInt.y[2])
    elseif generateurVisualise == k
        # archivage seulement pour le generateur k
        push!(d.XInt,vg[k].sInt.y[1])
        push!(d.YInt,vg[k].sInt.y[2])
    end

end

# ==============================================================================
# arrondi la solution correspondant au generateur (pas d'historique donc)
# version avec cone inferieur et superieur
function roundingSolutionnew24!(vg::Vector{tGenerateur}, k::Int64, c1::Array{Int,1}, c2::Array{Int,1}, d::tListDisplay)

    nbvar = length(vg[k].sInt.x)
    nbgen = size(vg,1)

    # --------------------------------------------------------------------------
    # identifie les variables fractionnaires et marquage par valeur sentinelle -1
    for i in 1:nbvar
        if  isapprox(vg[k].sPrj.x[i] , 0.0, atol=1e-3)
            vg[k].sInt.x[i] = 0
        elseif isapprox(vg[k].sPrj.x[i] , 1.0, atol=1e-3)
            vg[k].sInt.x[i] = 1
        else
            vg[k].sInt.x[i] = -1
        end
    end

    # --------------------------------------------------------------------------
    # Selectionne les points pour le cone pointe sur le generateur k (pCour) et ouvert vers Y
    pPrec, pCour, pSuiv = selectionPoints(vg, k)

    # --------------------------------------------------------------------------
    # elabore le point pointeur du cone ouvert vers L
    pC = elaborePointConeOuvertversL(vg, k, pPrec, pSuiv)

    # --------------------------------------------------------------------------
    # Arrondi les valeurs non-entieres d'une solution fractionnaire
    # 1) applique un arrondi qui fixe 1 variable (la premiere venue) par iteration

    vg[k].sInt.y[1] = 0
    vg[k].sInt.y[2] = 0

    # calcule en differentiel par rapport a l'image de la solution fractionnaire
    z1 = vg[k].sPrj.y[1]
    z2 = vg[k].sPrj.y[2]
    pM = tPoint( z1 , z2 )

#    print("Depart a arrondir : ")
#    @show pPrec, pCour, pSuiv, pM
#    if length(vg[k].sPrj.x) ≤ 20 @show vg[k].sPrj.x end
#    if length(vg[k].sInt.x) ≤ 20 @show vg[k].sInt.x end

    verbose ? @printf("  %2dR : [ %8.2f , %8.2f ] ", k, z1, z2) : nothing

    nbVarNonEntiere = 0
    for i in 1:nbvar
        if vg[k].sInt.x[i] == -1

            # la variable i est non entiere
            nbVarNonEntiere += 1

            # defalque la contribution fractionnaire de la variable aux objectifs => pose x[i]=0
            pM0 = tPoint( z1 - vg[k].sPrj.x[i] * c1[i] , z2 - vg[k].sPrj.x[i] * c2[i])
            pM1 = tPoint( z1 - vg[k].sPrj.x[i] * c1[i]+c1[i] , z2 - vg[k].sPrj.x[i] * c2[i]+c2[i])

            if !inCone1VersZ(pCour, pSuiv, pPrec, pM0)
                # le point pM0 obtenu est hors du cone inferieur => valide x[i]=1
                vg[k].sInt.x[i] = 1
                z1 = z1 - (vg[k].sPrj.x[i] * c1[i]) + c1[i]
                z2 = z2 - (vg[k].sPrj.x[i] * c2[i]) + c2[i]
#                @printf("  x[%2d]=1 |  : [ %12.5f , %12.5f ] \n",i, z1, z2)
            elseif inCone2Vers0(pC, pSuiv, pPrec, pM0)
                # le point pM0 obtenu est dans le cone => valide x[i]=0
                vg[k].sInt.x[i] = 0
                z1 = z1 - (vg[k].sPrj.x[i] * c1[i])
                z2 = z2 - (vg[k].sPrj.x[i] * c2[i])
#                @printf("  x[%2d]=0 |  : [ %12.5f , %12.5f ] \n",i, z1, z2)
            elseif inCone2Vers0(pC, pSuiv, pPrec, pM1)
                # le point pM1 obtenu est dans le cone => valide x[i]=1
                vg[k].sInt.x[i] = 1
                z1 = z1 - (vg[k].sPrj.x[i] * c1[i]) + c1[i]
                z2 = z2 - (vg[k].sPrj.x[i] * c2[i]) + c2[i]
#                @printf("  x[%2d]=0 |  : [ %12.5f , %12.5f ] \n",i, z1, z2)
            else
#                println("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
                # choisir entre pM0 et pM1 selon la distance a la droite [O;OO]
                O = tPoint(vg[k].sRel.y[1], vg[k].sRel.y[2]) # generateur
                OO = tPoint(vg[k].sPrj.y[1], vg[k].sPrj.y[2]) # projection
                # les points pM0 et pM1 obtenus sont hors du cone => choix de 0 ou 1 avec la plus courte distance à O et OO
                distM0 = sqrt((O.x-pM0.x)^2+(O.y-pM0.y)^2)+sqrt((OO.x-pM0.x)^2+(OO.y-pM0.y)^2)
                distM1 = sqrt((O.x-pM1.x)^2+(O.y-pM1.y)^2)+sqrt((OO.x-pM1.x)^2+(OO.y-pM1.y)^2)
                if distM0 < distM1
                    vg[k].sInt.x[i] = 0
                    z1 = z1 - (vg[k].sPrj.x[i] * c1[i])
                    z2 = z2 - (vg[k].sPrj.x[i] * c2[i])
                else
                    vg[k].sInt.x[i] = 1
                    z1 = z1 - (vg[k].sPrj.x[i] * c1[i]) + c1[i]
                    z2 = z2 - (vg[k].sPrj.x[i] * c2[i]) + c2[i]
                end
            end
        end

        # calcule la performance de la solution entiere
        vg[k].sInt.y[1] += vg[k].sInt.x[i] * c1[i]
        vg[k].sInt.y[2] += vg[k].sInt.x[i] * c2[i]
    end

    if length(vg[k].sInt.x) ≤ 20 @show vg[k].sInt.x end
    @printf("→ #round : %4d → [ %5d , %5d ] ", nbVarNonEntiere, vg[k].sInt.y[1], vg[k].sInt.y[2])

    # archive le point obtenu pour les besoins d'affichage
    if generateurVisualise == -1
        # archivage pour tous les generateurs
        push!(d.XInt,vg[k].sInt.y[1])
        push!(d.YInt,vg[k].sInt.y[2])
    elseif generateurVisualise == k
        # archivage seulement pour le generateur k
        push!(d.XInt,vg[k].sInt.y[1])
        push!(d.YInt,vg[k].sInt.y[2])
    end

end

# ==============================================================================
# arrondi la solution correspondant au generateur (pas d'historique donc)
# version avec voisinage et selection d'un voisin selon distance L1 avec generateur
function roundingSolutionNew23!(vg::Vector{tGenerateur}, k::Int64, c1::Array{Int,1}, c2::Array{Int,1}, d::tListDisplay)

    nbvar = length(vg[k].sInt.x)
    nbgen = size(vg,1)
    lstIdxFrac =(Int64)[]

    vg[k].sInt.y[1] = 0
    vg[k].sInt.y[2] = 0

    # --------------------------------------------------------------------------
    # identifie les variables fractionnaires et marquage par valeur sentinelle -1
    for i in 1:nbvar
        if  isapprox(vg[k].sPrj.x[i] , 0.0, atol=1e-3)
            vg[k].sInt.x[i] = 0
        elseif isapprox(vg[k].sPrj.x[i] , 1.0, atol=1e-3)
            vg[k].sInt.x[i] = 1
            vg[k].sInt.y[1] += c1[i] # Comptabilise la contribution des var a 1
            vg[k].sInt.y[2] += c2[i]
        else
            vg[k].sInt.x[i] = -1
            push!(lstIdxFrac,i)
        end
    end

    # --------------------------------------------------------------------------
    # Selectionne les points pour le cone pointe sur le generateur k (pCour) et ouvert vers Y
    pPrec, pCour, pSuiv = selectionPoints(vg, k)

    # --------------------------------------------------------------------------
    # elabore le point pointeur du cone ouvert vers L
    pC = elaborePointConeOuvertversL(vg, k, pPrec, pSuiv)

    # --------------------------------------------------------------------------
    # Arrondi les valeurs non-entieres d'une solution fractionnaire
    # 1) applique un arrondi qui fixe 1 variable (la premiere venue) par iteration

    # calcule en differentiel par rapport a l'image de la solution fractionnaire
    z1 = vg[k].sPrj.y[1]
    z2 = vg[k].sPrj.y[2]
    pM = tPoint( z1 , z2 )

#    if length(vg[k].sPrj.x) ≤ 20 @show vg[k].sPrj.x end
#    if length(vg[k].sInt.x) ≤ 20 @show vg[k].sInt.x end

    nbVarNonEntiere = length(lstIdxFrac)
    # traite toutes les variables non entieres

#=
    @printf("\n\n")
    @show c1
    @show c2
    @printf("VARIABLES FIXEES : \n")
    @show vg[k].sPrj.x
    @show vg[k].sInt.x
    @show vg[k].sInt.y
    @show pM
    @show lstIdxFrac
=#
    verbose ? @printf("  %2dR : [ %8.2f , %8.2f ] ", k, z1, z2) : nothing

    # traitement de toutes les variables fractionnaires une a une
    for i= 1:nbVarNonEntiere
#        @printf("\n")

        distMin = Inf; ind = -1; val = -1
        # voisinage a 1 variable sur toutes les variables
        for iFrac in [i for i=1:nbvar if vg[k].sInt.x[i] == -1]

            # Evalue x[iFrac]=0
            # defalque la contribution fractionnaire de la variable aux objectifs => pose x[iFrac]=0
            pM = tPoint( z1 - c1[iFrac] * vg[k].sPrj.x[iFrac] , z2 - c2[iFrac] * vg[k].sPrj.x[iFrac])
            if inCone1VersZ(pCour, pSuiv, pPrec, pM)
                distL1 = abs(vg[k].sRel.y[1]-pM.x) + abs(vg[k].sRel.y[2]-pM.y)
                if distL1 < distMin
                    distMin = distL1; ind = iFrac; val = 0
                end
            end

            # Evalue x[iFrac]=1
            # ajoute la contribution entiere de la variable aux objectifs => pose x[iFrac]=1
            pM = tPoint( z1 - c1[iFrac] * vg[k].sPrj.x[iFrac] + c1[iFrac], z2 - c2[iFrac] * vg[k].sPrj.x[iFrac] + c2[iFrac])
            distL1 = abs(vg[k].sRel.y[1]-pM.x) + abs(vg[k].sRel.y[2]-pM.y)
            if distL1 < distMin
                distMin = distL1; ind = iFrac; val = 1
            end
        end

        # applique le voisin le plus interessant trouve
        vg[k].sInt.x[ind] = val
        vg[k].sInt.y[1] += vg[k].sInt.x[ind] * c1[ind] * val
        vg[k].sInt.y[2] += vg[k].sInt.x[ind] * c2[ind] * val
        z1 = z1 - (vg[k].sPrj.x[ind] * c1[ind]) + c1[ind] * val
        z2 = z2 - (vg[k].sPrj.x[ind] * c2[ind]) + c2[ind] * val
    #    @printf("  x[%2d]=%2d |  : [ %12.5f , %12.5f ] \n",ind, val, z1, z2)

    end

#    if length(vg[k].sInt.x) ≤ 20 @show vg[k].sInt.x end
    @printf("→ #round : %4d → [ %5d , %5d ] ", nbVarNonEntiere, vg[k].sInt.y[1], vg[k].sInt.y[2])

    # archive le point obtenu pour les besoins d'affichage
    if generateurVisualise == -1
        # archivage pour tous les generateurs
        push!(d.XInt,vg[k].sInt.y[1])
        push!(d.YInt,vg[k].sInt.y[2])
    elseif generateurVisualise == k
        # archivage seulement pour le generateur k
        push!(d.XInt,vg[k].sInt.y[1])
        push!(d.YInt,vg[k].sInt.y[2])
    end

end

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)
sess = Session(); init(sess)

# TODO: specify your input parameters

nc = 100
n = 10
xc = rand(nc)
yc = rand(nc)
zc = rand(nc)
x = rand(n)
y = rand(n)
z = rand(n)
c = rand(nc)
d = rand(4)
e = rand(nc)

# problem with eps, xc
for kind in [0,1,2,3]
    # TODO: change your test parameter to `m`
    #       in the case of `multiple=true`, you also need to specify which component you are testings
    # gradient check -- v
    function scalar_function(m)
        rbf = RBF3D(xc, m, zc, c = c, d = d, eps = e, kind = kind)
        return sum(rbf(x,y,z)^2)
    end

    # TODO: change `m_` and `v_` to appropriate values
    m_ = constant(rand(size(xc)...))
    v_ = rand(size(m_)...)
    y_ = scalar_function(m_)
    dy_ = gradients(y_, m_)
    ms_ = Array{Any}(undef, 5)
    ys_ = Array{Any}(undef, 5)
    s_ = Array{Any}(undef, 5)
    w_ = Array{Any}(undef, 5)
    gs_ =  @. 1 / 10^(1:5)

    for i = 1:5
        g_ = gs_[i]
        ms_[i] = m_ + g_*v_
        ys_[i] = scalar_function(ms_[i])
        s_[i] = ys_[i] - y_
        w_[i] = s_[i] - g_*sum(v_.*dy_)
    end

    sess = Session(); init(sess)
    sval_ = run(sess, s_)
    wval_ = run(sess, w_)
    close("all")
    loglog(gs_, abs.(sval_), "*-", label="finite difference")
    loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
    loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
    loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

    plt.gca().invert_xaxis()
    legend()
    xlabel("\$\\gamma\$")
    ylabel("Error")
    savefig("test$kind.png")
end

ext: jl


text:
function EuropeanSpreadOptionMC(K, T, S1, S2, sig1, sig2, div1, div2, rho, r, N, M)

## setting parameters

dt = T/N
nu1dt = (r - div1 - 0.5 * sig1^2) * dt
nu2dt = (r - div2 - 0.5 * sig2^2) * dt
sig1sdt = sig1 * sqrt(dt)
sig2sdt = sig2 * sqrt(dt)
srho = sqrt(1 - rho^2)

## variables to store the sums of the option values

sum_CT = 0
sum_CT2 = 0

## outer loop is the number of simulations bigger is better

for j = 1:M

	St1 = S1
	St2 = S2

	## inner loop is simulating each price path of length N
	for i = 1:N

		error1 = randn()
		error2 = randn()

		z1 = error1
		z2 = rho * error1 + srho * error2

		St1 = St1 * exp(nu1dt + sig1sdt * z1)
		St2 = St2 * exp(nu2dt + sig2sdt * z2)

	end

    	## wraps up and stores all the values of each path
	CT = max(0, St1 - St2 - K)
	sum_CT = sum_CT + CT
	sum_CT2 = sum_CT2 + CT^2

end

## discounting the average call value
call_value = sum_CT / M * exp(-r*T)

SD = (sqrt(sum_CT2 - 1/M * (sum_CT^2)) * exp(-2*r*T)) / (M-1)
SE = SD / sqrt(M)

println("The Call Value is $call_value")
println("The SD is $SD")
println("The SE is $SE")

MC = [call_value, SD, SE]

end


text:
module ClusterUtils

export lookup, filterpaths, dictify, save, load
export describepids, localpids, invtopology, chunkit
export Doable, sow, reap, reaprefs
export @reap, @sow
export getrepresentation, display, broadcast_shared

#export collectmsgs, collectmsgsatpid, collectmsgsatmaster, swapmsgs
#export stitch

include("topology.jl") #finding, decribing and setting vars on other pids
include("reap.jl") #main message passing utilities
include("utils.jl") #misc utilities
include("experimental.jl") #other styles of message passing, not tested as much
include("patches.jl") #patches for issues to do with HPC

end  #module ClusterUtils


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

const VarType      = Pair{Symbol,DataType}
const VarOrVarType = Union{Symbol,VarType}

"""
    SimulationProblem(sdata, sdomain, vars, nreals)
    SimulationProblem(sdomain, vars, nreals)

A spatial simulation problem on a given spatial domain `sdomain`
in which the variables to be simulated are listed in `vars`.

For conditional simulation, the data of the problem is stored in
spatial data `sdata`.

For unconditional simulation, a list of pairs `vars` must be
provided mapping variable names to their types.

In both cases, a number `nreals` of realizations is requested.

## Examples

Create a conditional simulation problem for porosity and permeability
with 100 realizations:

```julia
julia> SimulationProblem(sdata, sdomain, (:porosity,:permeability), 100)
```

Create an unconditional simulation problem for porosity and facies type
with 100 realizations:

```julia
julia> SimulationProblem(sdomain, (:porosity => Float64, :facies => Int), 100)
```
"""
struct SimulationProblem{S,D,N} <: AbstractProblem
  sdata::S
  sdomain::D
  vars::NTuple{N,Variable}
  nreals::Int

  function SimulationProblem{S,D,N}(sdata, sdomain, vars, nreals) where {S,D,N}
    pnames = name.(vars)

    @assert !isempty(pnames) "target variables must be specified"
    @assert nreals > 0 "number of realizations must be positive"

    new(sdata, sdomain, vars, nreals)
  end
end

function SimulationProblem(sdata::S, sdomain::D, vars::NTuple{N,VarOrVarType}, nreals::Int) where {S,D,N}
  @assert coordtype(sdata) == coordtype(sdomain) "data and domain must have the same coordinate type"

  datavars = Dict(name(var) => mactype(var) for var in variables(sdata))

  # pairs with variable names and types
  varstypes = map(vars) do vt
    if vt isa Symbol # for variables without type, find the type in spatial data
      var = vt
      if var ∈ keys(datavars)
        var => datavars[var]
      else
        @error "please specify the type of target variable $var"
      end
    else # for variables with type, make sure the type is valid
      var, T = vt
      if var ∈ keys(datavars)
        U = datavars[var]
        @assert U <: T "type $T for variable $var cannot hold values of type $U in spatial data"
      end
      vt
    end
  end

  vars = Tuple(Variable(var, V) for (var, V) in varstypes)

  SimulationProblem{S,D,N}(sdata, sdomain, vars, nreals)
end

SimulationProblem(sdata::S, sdomain::D, var::VarOrVarType, nreals::Int) where {S,D} =
  SimulationProblem(sdata, sdomain, (var,), nreals)

SimulationProblem(sdomain::D, vars::NTuple{N,VarType}, nreals::Int) where {D,N} =
  SimulationProblem{Nothing,D,N}(nothing, sdomain, Tuple(Variable(n, t) for (n, t) in vars), nreals)

SimulationProblem(sdomain::D, var::VarType, nreals::Int) where {D} =
  SimulationProblem(sdomain, (var,), nreals)

"""
    data(problem)

Return the spatial data of the simulation `problem`.
"""
data(problem::SimulationProblem) = problem.sdata

"""
    domain(problem)

Return the spatial domain of the simulation `problem`.
"""
domain(problem::SimulationProblem) = problem.sdomain

"""
    variables(problem)

Return the target variables of the simulation `problem` and their types.
"""
variables(problem::SimulationProblem) = problem.vars

"""
    hasdata(problem)

Return `true` if simulation `problem` has data.
"""
hasdata(problem::SimulationProblem) = !isnothing(problem.sdata)

"""
    nreals(problem)

Return the number of realizations of the simulation `problem`.
"""
nreals(problem::SimulationProblem) = problem.nreals

# ------------
# IO methods
# ------------
function Base.show(io::IO, problem::SimulationProblem)
  N = ncoords(problem.sdomain)
  kind = hasdata(problem) ? "conditional" : "unconditional"
  print(io, "$(N)D SimulationProblem ($kind)")
end

function Base.show(io::IO, ::MIME"text/plain", problem::SimulationProblem)
  vars = ["$(name(v)) ($(mactype(v)))" for v in problem.vars]
  println(io, problem)
  if problem.sdata ≠ nothing
    println(io, "  data:      ", problem.sdata)
  end
  println(io, "  domain:    ", problem.sdomain)
  println(io, "  variables: ", join(vars, ", ", " and "))
  print(  io, "  N° reals:  ", problem.nreals)
end


text:
## index.jl : representation of a Level and an Index
#
# Representation of levels and multi-dimensional indices.
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

## Index ##
"""
    Index(i::Integer...)

Return a multi-index.

# Examples
```jldoctest; setup = :(using MultilevelEstimators)
julia> index = Index(2, 1)
(2, 1)
```
See also: [`Level`](@ref)
"""
Index = CartesianIndex

show(io::IO, index::Index) = print(io, string("(", join(index.I, ", "), ")"))
show(io::IO, ::Type{Index}) = print(io, "Index")

## Level ##
"""
    Level(l::Integer)

Return a level.

# Examples
```jldoctest; setup = :(using MultilevelEstimators)
julia> level = Level(2)
2
```
See also: [`Index`](@ref)
"""
Level = CartesianIndex{1}

show(io::IO, level::Level) = print(io, level[1])
show(io::IO, ::Type{Level}) = print(io, "Level")

## utilities ##
function diff(index::Index{d}) where d
    D = Dict{Index{d}, Int}()
    Istart = max(zero(index), index - one(index))
    Iend = index
	R = Istart:Iend
	for I in Base.Iterators.take(R, length(R)-1)
        D[I] = isodd(sum(I - index)) ? -1 : 1
    end
    return D
end

sum(I::CartesianIndex) = sum(Tuple(I))

if VERSION < v"1.1.0-"
	(:)(I::CartesianIndex{N}, J::CartesianIndex{N}) where N = CartesianIndices(map((i,j) -> i:j, Tuple(I), Tuple(J)))
end


text:
using Whirl
using Profile
using ProfileView

Profile.init(delay = 1e-2)

stepper = build(sod_shock_tube_builder(number_of_particles = 1000))

function profiling(n)
   for _ = 1:n
      step!(stepper, 1e-4)
   end
end

@profview profiling(1)
@profview profiling(10000)

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
module Levels

include("optics_utils.jl")

include("states.jl")

include("species.jl")
include("species_data.jl")

include("rates.jl")

end # module


text:
"""
    dataprocess(data_file; target_satnum=missing, T=10e-3, t_length=4.,
                output_dir=missing, prns=:all, data_start_t=5e-3,
                fd_range=5000., fd_center=0., Δdays=5, showplot=true,
                site_lla=(40.01, -105.2437, 1655), figsize=missing,
                use_gps_orbit_only=false)

Perform course acquisition and/or tracking over entire data file using
coherent integration `T` with length of data to be processed being defined
by `t_length`.
"""
function dataprocess(data_file; target_satnum=missing, T=1e-3, t_length=4.,
                     output_dir=missing, prns=:all, skip_to=0.1,
                     fd_range=5000., fd_center=0., Δdays=5, showplot=true,
                     site_lla=(40.01, -105.2437, 1655), figsize=missing,
                     use_gps_orbit_only=false)
    # Load data & infer f_s, f_if, and data type from file name
    file_info = data_info_from_name(data_file)
    println("-----------------------------------------------------")
    println("                      File Info")
    println("-----------------------------------------------------")
    println("Timestamp (UTC):            $(file_info.timestamp)")
    println("Sampling frequency (Hz):    $(file_info.f_s)")
    println("Signal type:                $(gnsstypes[file_info.sigtype])")
    println("Data center frequency (Hz): $(file_info.sig_freq)")
    println("IF frequency (Hz):          $(file_info.f_if)")
    println("Data type:                  $(gnsstypes[file_info.data_type])")
    println("-----------------------------------------------------")
    start_data_idx = Int(file_info.f_s * data_start_t)+1
    data = loaddata(file_info.data_type, data_file, file_info.f_s,
                    file_info.f_if, T;
                    skip_to=skip_to)
    # Initialize replica struct
    replica = definesignal(file_info.sigtype, file_info.f_s, T)
    # Calculate Doppler bin spacing for course acquisition
    Δfd = 1/replica.t_length  # Hz
    # Allocate space for correlation result
    corr_result = gencorrresult(fd_range, Δfd, replica.sample_num)
    # Get list of PRNs and corresponding NORAD IDs at given time
    if prns == :all
        if (file_info.sigtype == Val(:l5q)) || (file_info.sigtype == Val(:l5i))
            prns = [1, 24, 27, 30, 6, 9, 3, 26, 8, 10, 32]  # all Block IIR-F satellites
        else
            prns = Array(1:32)
        end
    end
    gps_satnums = getGPSSatnums(file_info.timestamp_JD, prns)
    prns = collect(keys(gps_satnums))
    # Initialize orbit struct if target object TLE given
    if ismissing(target_satnum) && ~use_gps_orbit_only
        use_orbitinfo = false
    else
        use_orbitinfo = true
        if use_gps_orbit_only
            satnum_list = collect(values(gps_satnums))
        else
            satnum_list = [collect(values(gps_satnums)); target_satnum]
        end
        # Get TLEs for each PRN and target satnum
        tles = getTLEs(file_info.timestamp_JD,
                       satnum_list,
                       Δdays=Δdays)
        orbitinfo = Dict{Int,OrbitInfo}()
        for prn in prns
            if use_gps_orbit_only
                orbitinfo[prn] = initorbitinfo(tles[gps_satnums[prn]],
                                               file_info.timestamp,
                                               site_lla)
            else
                orbitinfo[prn] = initorbitinfo(tles[gps_satnums[prn]],
                                               tles[target_satnum],
                                               file_info.timestamp,
                                               site_lla)
            end
        end
    end
    # Allocate space for t, fd_est, n₀_est, and SNR_est arrays
    N = Int(floor(t_length/T))
    t = calctvector(N, 1/T)
    results = Dict{Int,Dict{String,Array{Float64}}}()
    for prn in prns
        results[prn] = Dict("n0_est"=>Array{Float64}(undef, N),
                            "fd_est"=>Array{Float64}(undef, N),
                            "SNR_est"=>Array{Float64}(undef, N),
                            "fd_exp"=>Array{Float64}(undef, N))
    end
    # Perform course acquisition on data at every Nᵗʰ T segment for each
    # `prn` in `prns`
    p = Progress(N*length(prns), 1, "Processing...")
    for n in 1:N
        Δt_JD_low = t[n]/60/60/24  # Days
        Δt_JD_cen = (t[n] + T/2)/60/60/24  # Days
        Δt_JD_high = (t[n] + T)/60/60/24  # Days
        for prn in prns
            # Calculate Doppler at (t + T/2), the middle of the replica signal
            # Assume that the Doppler is constant
            if use_orbitinfo
                if use_gps_orbit_only
                    fd_exp_low = calcdoppler(orbitinfo[prn].orb,
                                             orbitinfo[prn].start_time_julian_date+Δt_JD_low,
                                             orbitinfo[prn].eop,
                                             orbitinfo[prn].site_loc_ecef,
                                             file_info.sig_freq)
                    fd_exp_high = calcdoppler(orbitinfo[prn].orb,
                                              orbitinfo[prn].start_time_julian_date+Δt_JD_high,
                                              orbitinfo[prn].eop,
                                              orbitinfo[prn].site_loc_ecef,
                                              file_info.sig_freq)
                else
                    fd_exp_low = calcdoppler(orbitinfo[prn].orb[1],
                                             orbitinfo[prn].orb[2],
                                             orbitinfo[prn].start_time_julian_date+Δt_JD_low,
                                             orbitinfo[prn].eop,
                                             orbitinfo[prn].site_loc_ecef,
                                             file_info.sig_freq)
                    fd_exp_high = calcdoppler(orbitinfo[prn].orb[1],
                                              orbitinfo[prn].orb[2],
                                              orbitinfo[prn].start_time_julian_date+Δt_JD_high,
                                              orbitinfo[prn].eop,
                                              orbitinfo[prn].site_loc_ecef,
                                              file_info.sig_freq)
                end
                fd_rate = (fd_exp_high - fd_exp_low)/T
            else
                fd_exp_low = fd_center
                fd_rate = 0.
            end
            # Perform course acquisition
            courseacquisition!(corr_result, data, replica, prn;
                               fd_center=fd_exp_low, fd_range=fd_range,
                               fd_rate=fd_rate, Δfd=Δfd, threads=nthreads(),
                               showprogressbar=false)
            # Estimate n₀_est, fd_est, and SNR_est from course acquisition
            # result
            n₀, f_d, snr = course_acq_est(corr_result, fd_exp_low, fd_range, Δfd)
            results[prn]["n0_est"][n] = n₀
            results[prn]["fd_est"][n] = f_d
            results[prn]["fd_exp"][n] = fd_exp_low
            results[prn]["SNR_est"][n] = snr
            next!(p)
        end
        # reloaddata!(data, start_data_idx+n*N+1)
        data = loaddata(file_info.data_type, data_file, file_info.f_s,
                        file_info.f_if, T;
                        skip_to=skip_to+n*T)
    end
    # Save results to HDF5 file
    if ~ismissing(output_dir)
        h5open(output_dir, "w") do file
            for prn in prns
                write(file, "$prn/n0_est", results[prn]["n0_est"])
                write(file, "$prn/fd_est", results[prn]["fd_est"])
                write(file, "$prn/SNR_est", results[prn]["SNR_est"])
                write(file, "$prn/fd_exp", results[prn]["fd_exp"])
            end
        end
    end
    # Plot results for only 1 PRN. More than 1 are not plotted.
    if showplot && (length(prns) == 1)
        if t_length > 60
            t .*= 1/60
            t_label = "Minutes"
        else
            t_label = "Seconds"
        end
        if ~ismissing(figsize)
            fig = figure(figsize=figsize)
        else
            fig = figure()
        end
        ax1 = subplot(3, 1, 1)
        plot(t, results[prns[1]]["n0_est"], "k.")
        xlabel("Time ($t_label)")
        ylabel("n₀ (Samples)")
        ax2 = subplot(3, 1, 2)
        plot(t, results[prns[1]]["fd_est"]./1000, "k.")
        xlabel("Time ($t_label)")
        ylabel("Doppler Frequency (kHz)")
        ax3 = subplot(3, 1, 3)
        plot(t, results[prns[1]]["SNR_est"], "k.")
        xlabel("Time ($t_label)")
        ylabel("Correlation Peak SNR (dB)")
        suptitle("PRN $(prns[1])")
        show()
    end
    return results
end


text:
export CenterFirst, CenterLast, CenterMiddle, CenterFT
export center

"""
    Center

Abstract supertype for all Center types.

See [`CenterFirst`](@ref), [`CenterLast`](@ref),
[`CenterMiddle`](@ref), [`CenterFT`](@ref).
"""
abstract type Center end

"""
    CenterFirst

Type to indicate that the center should be
at the first entry of the array.
This corresponds to the indices `(1,1,1,...)`.
"""
struct CenterFirst <: Center end

"""
    CenterEnd

Type to indicate that the center should be
at the last entry of the array.
This corresponds in an array with size `(5,4,3)` to
the indices `(5,4,3)`.
"""
struct CenterLast <: Center end

"""
    CenterMiddle

Type to indicate that the center should be
at the mathematical center of the array (if interpreted
as a N dimensional volume).
This corresponds in an array with size `(5,4,3)` to
the indices `(3,2.5,2)`.
"""
struct CenterMiddle <: Center end

"""
    CenterFT

Type to indicate that the center should be
at the center defined in the FFT sense.
This corresponds in an array with size `(5,4,3)` to
the indices `(3,3,2)`.
"""
struct CenterFT <: Center end

"""
    center(sz::NTuple{N, T}, ::Type{<:Center})

Return the corresponding center of an array with size `sz`.
Depending on the `Center` type the center is chosen.
See [`CenterFirst`](@ref), [`CenterLast`](@ref),
[`CenterMiddle`](@ref), [`CenterFT`](@ref).

```jldoctest
julia> center((1,2,3,4), CenterFirst)
(1, 1, 1, 1)

julia> center((1,2,3,4), CenterLast)
(1, 2, 3, 4)

julia> center((1,2,3,4), CenterMiddle)
(1.0, 1.5, 2.0, 2.5)

julia> center((1,2,3,4), CenterFT)
(1, 2, 2, 3)
```
"""
center(size::NTuple{N, T},
       ::Type{CenterFirst}) where {T, N} = ntuple(i -> one(T), Val(N))

center(size::NTuple{N, T},
       ::Type{CenterLast}) where {T, N} = ntuple(i -> size[i], Val(N))

center(size::NTuple{N, T},
       ::Type{CenterMiddle}) where {T, N} = ntuple(i -> (1 + size[i]) / 2, Val(N))

center(size::NTuple{N, T},
       ::Type{CenterFT}) where {T, N} = ntuple(i -> size[i] ÷ 2 + 1, Val(N))

"""
    center(sz::NTuple{N, T}, j<:Number)

Return an tuple with the same length as `sz` and
with entries `j`.

 # Examples
```jldoctest
julia> center((2,2), 1)
(1, 1)

julia> center((2,2, 5), 3)
(3, 3, 3)
```
"""
center(size::NTuple{N, T},
       j::Number) where {T, N} = ntuple(i -> j, Val(N))

"""
    center(sz::NTuple{N, T}, ctr::NTuple{N, T})

Return ctr.

 # Examples
```jldoctest
julia> center((1,2), (1,1))
(1, 1)

julia> center((1,3), (1,1))
(1, 1)
```
"""
center(size::NTuple{N, T},
       ctr::NTuple{N, T}) where {T, N} = ctr


text:
module CommunalHelperDreamTunnelEntry

using ..Ahorn, Maple
using Ahorn.CommunalHelper

@mapdef Entity "CommunalHelper/DreamTunnelEntry" DreamTunnelEntry(
    x::Integer,
    y::Integer,
    overrideAllowStaticMovers=false,
    depth::Integer=-13000,
)

const placements = Ahorn.PlacementDict(
    "Dream Tunnel Entry ($dir) (Communal Helper)" => Ahorn.EntityPlacement(
        DreamTunnelEntry,
        "rectangle",
        Dict{String,Any}(
            "orientation" => dir,
        ),
    ) for dir in Maple.spike_directions
)

Ahorn.editingOptions(entity::DreamTunnelEntry) = Dict{String,Any}(
    "orientation" => Maple.spike_directions,
    "depth" => CommunalHelper.depths,
)

Ahorn.minimumSize(entity::DreamTunnelEntry) = 8, 8

const resizeDirections = Dict{String,Tuple{Bool,Bool}}(
    "Up" => (true, false),
    "Down" => (true, false),
    "Left" => (false, true),
    "Right" => (false, true),
)

function Ahorn.resizable(entity::DreamTunnelEntry)
    orientation = get(entity.data, "orientation", "Up")
    return resizeDirections[orientation]
end

function Ahorn.selection(entity::DreamTunnelEntry)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, width, height)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::DreamTunnelEntry)
    orientation = get(entity.data, "orientation", "Up")

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    x, y = Ahorn.position(entity)
    ox, oy = x + width * (orientation == "Right"), y + height * (orientation == "Down")
    cx, cy = x + width * in(orientation, ("Up", "Down", "Right")),
    y + height * in(orientation, ("Down", "Left", "Right"))

    Ahorn.drawRectangle(ctx, x, y, width, height, (0.0, 0.0, 0.0, 0.4))
    Ahorn.drawLines(ctx, ((ox, oy), (cx, cy)), (1.0, 1.0, 1.0, 1.0), thickness=1)
end

end


text:
export PseudoPeriodicTwin

import TimeseriesSurrogates: Surrogate, SurrogateGenerator, surrogenerator
import StatsBase: sample, pweights
import Distances: Euclidean, PreMetric, pairwise
import DelayEmbeddings: Dataset, embed

"""
    PseudoPeriodicTwin(d::Int, τ::Int, δ = 0.2, ρ = 0.1, metric = Euclidean())
    PseudoPeriodicTwin(δ = 0.2, ρ = 0.1, metric = Euclidean())

A pseudoperiodic twin surrogate[^Miralles2015], which is a fusion of the
twin surrogate[^Thiel2006] and the pseudo-periodic surrogate[^Small2001].

## Input parameters

A delay reconstruction of the input timeseries is constructed using embedding
dimension `d` and embedding delay `τ`. The threshold `δ ∈ (0, 1]` determines which
points are "close" (neighbors) or not, and is expressed as a fraction of the
attractor diameter, as determined by the input data. The authors of the original
twin surrogate paper recommend `0.05 ≤ δ ≤ 0.2`[^Thiel2006].

If you have pre-embedded your timeseries, and timeseries is already a `::Dataset`, use the
three-argument constructor (so that no delay reconstruction is performed).
If you want a surrogate for a scalar-valued timeseries, use the five-argument constructor
to also provide the embedding delay `τ` and embedding dimension `d`.

## Null hypothesis

Pseudo-periodic twin surrogates generate signals similar to the original data if the
original signal is (quasi-)periodic. If the original signal is not
(quasi-)periodic, then these surrogates will have different
recurrence plots than the original signal, but preserve the overall
shape of the attractor. Thus, `PseudoPeriodicTwin` surrogates can be used to test
null hypothesis that the observed timeseries (or orbit) is consistent
with a quasi-periodic orbit[^Miralles2015].

## Returns

A `d`-dimensional surrogate orbit (a `Dataset`) is returned. Sample the first
column of this dataset if a scalar-valued surrogate is desired.

[^Small2001]: Small et al., Surrogate test for pseudoperiodic timeseries data, [Physical Review Letters, 87(18)](https://doi.org/10.1103/PhysRevLett.87.188101)
[^Thiel2006]: Thiel, Marco, et al. "Twin surrogates to test for complex synchronisation." EPL (Europhysics Letters) 75.4 (2006): 535.
[^Miralles2015]: Miralles, R., et al. "Characterization of the complexity in short oscillating timeseries: An application to seismic airgun detonations." The Journal of the Acoustical Society of America 138.3 (2015): 1595-1603.
"""
struct PseudoPeriodicTwin{T<:Real, P<:Real, D<:PreMetric} <: Surrogate
    d::Union{Nothing, Int}
    τ::Union{Nothing, Int}
    δ::T
    ρ::P
    metric::D

    function PseudoPeriodicTwin(d::Int, τ::Int, δ::T = 0.2, ρ::P = 0.1, metric::D = Euclidean()) where {T, P, D}
        new{T, P, D}(d, τ, δ, ρ, metric)
    end

    function PseudoPeriodicTwin(δ::T = 0.2, ρ::P = 0.1, metric::D = Euclidean()) where {T, P, D}
        new{T, P, D}(nothing, nothing, δ, ρ, metric)
    end
end

"""
    _prepare_embed(x::AbstractVector, d, τ) → Dataset
    _prepare_embed(x::Dataset, d, τ) → Dataset

Prepate input data for surrogate generation. If input is a vector, embed it using
the provided parameters. If input as a dataset, we assume it already represents an
orbit.
"""
function _prepare_embed end
_prepare_embed(x::AbstractVector, d, τ) = embed(x, d, τ)
_prepare_embed(x::Dataset, d, τ) = x

function surrogenerator(x::Union{AbstractVector, Dataset}, pp::PseudoPeriodicTwin, rng = Random.default_rng())
    d, τ, δ, metric = getfield.(Ref(pp), (:d, :τ, :δ, :metric))
    ρ = getfield.(Ref(pp), (:ρ))

    pts = _prepare_embed(x, d, τ)

    Nx = length(x)
    Npts = length(pts)

    dists = pairwise(metric, Matrix(pts), dims = 1)
    normalisedδ = δ*maximum(dists)

    T = eltype(pts)
    R = zeros(T, Npts, Npts)

    # Recurrence matrix
    for j = 1:Npts
        for i = 1:Npts
            R[i, j] = normalisedδ - dists[i, j] >= 0 ? 1.0 : 0.0
        end
    end

    # Identify twins
    #println("R contains $(count(R .== 1)/length(R)*100)% black dots")
    twins_i = Vector{Int}(undef, 0)
    twins_j = Vector{Int}(undef, 0)
    for j = 1:Npts
        for i = 1:Npts
            if i !== j && all(R[:, i] .≈ R[:, j])
                push!(twins_i, i)
                push!(twins_j, j)
            end
        end
    end
    #println("Found $(length(twins)) twins")

    twins = Dict{Int,Vector{Int}}()

    # For every point that has a twin, store the indices of all of its twins
    for twi in unique(twins_i)
        twins[twi] = twins_j[findall(twins_i .== twi)]
    end
    for twj in unique(twins_j)
        twins[twj] = twins_i[findall(twins_j .== twj)]
    end

    # Sampling weights (exclude the point itself)
    W = [pweights(exp.(-dists[setdiff(1:Npts, i), i] / ρ)) for i = 1:Npts]

    # The surrogate will be a vector of vectors (if pts is a Dataset, then
    # the eltype is SVector).
    PT = eltype(pts.data)
    s = Vector{PT}(undef, Nx)

    init = (pts = pts, Nx = Nx, Npts = Npts, dists = dists, R = R, twins = twins, W = W)
    return SurrogateGenerator(pp, x, s, init, rng)
end

function (sg::SurrogateGenerator{<:PseudoPeriodicTwin})()
    pts, Nx, Npts, twins, W = getfield.(Ref(sg.init), (:pts, :Nx, :Npts, :twins, :W))
    ρ = getfield.(Ref(sg.method), (:ρ))
    s = sg.s

    # Randomly pick a point from the state space as the starting point for the surrogate.
    n = 1
    i = rand(sg.rng, 1:Npts)
    s[n] = pts[i]

    while n < Nx
        # Look for possible twins of the point xᵢ. If any twin exists, jump to one of them
        # with probability 1/nⱼ, where nⱼ are the number of twins for the point xᵢ.
        if haskey(twins, i)
            # sample uniformly (with probability 1/ntargettwins) over possible target twins
            j = twins[i][rand(sg.rng, 1:length(twins[i]))]
            s[n] = pts[j]
            i = j
            n += 1
        end

        # The orbit moves on from the current point xᵢ to a randomly selected point
        # on the attractor. Closer points are more likely to be selected, and points
        # further away are less likely to be selected. The sampling probability
        # of the next point xⱼ decreases exponentially with increasing distance
        # from the current point xᵢ. Probabilities for jumping from xᵢ to any other
        # point have been pre-computed, and are stored in W[i].
        j = sample(sg.rng, 1:Npts, W[i])
        s[n] = pts[j]
        i = j
        n += 1
    end

    s[Nx] = pts[sample(1:Npts, W[i])[1]]

    return Dataset(s)
end


text:
# This file is a part of ROOTFramework.jl, licensed under the MIT License (MIT).

import Compat.Test
Test.@testset "Package ROOTFramework" begin
    include("test_json.jl")
    include("test_ttree.jl")
end


text:
using Pkg
Pkg.activate(pwd())

# # Exception handling

function fact(n)
    isinteger(n) && n >= 0 || error("argument must be non-negative integer")
    return n == 0 ? 1 : n * fact(n - 1)
end

#+

fact(1.4)
fact(-5)

#+

function fact(n)
    isinteger(n) && n >= 0 || throw(DomainError(n, "argument must be non-negative integer"))
    return n == 0 ? 1 : n * fact(n - 1)
end

#+

fact(1.4)
fact(-5)
fact("a")

#+

function fact_new(n::Real)
    isinteger(n) && n >= 0 || throw(DomainError(n, "argument must be non-negative integer"))
    return n == 0 ? 1 : n * fact(n - 1)
end

#+

fact_new("aaa")
methods(fact_new)


text:
# options checking specific to the Euler Equations

"""
  This function checks the options in the options dictionary after default
  values have been supplied and throws exceptions if unsupported/incompatable
  options are specified

  **Inputs**

   * opts: options dictionary
   * comm: MPI communicator

  **Outputs**

    none
"""
function checkOptions(opts, comm=MPI.COMM_WORLD)

  if opts["physics"] != PhysicsName
    error("physics not specified as $PhysicsName, are you lost?")
  end

  if opts["use_edge_res"]
    error("use_edge_res does not work correctly, do not use")
  end

  if opts["volume_integral_type"] != 1
    integral_type = opts["volume_integral_type"]
    error("SimpleODE equation supports volume integral type 1 only, type $integral_type not supported")
  end

  if opts["face_integral_type"] != 1
    integral_type = opts["face_integral_type"]
    error("SimpleODE equation supports face integral type 1 only, type $integral_type not supported")
  end

  if opts["use_edgestab"]
    error("SimpleODE does not support edge stabilization")
  end

  if opts["use_filter"]
    error("SimpleODE does not support filter stabilization")
  end

  if opts["use_dissipation"] || opts["use_dissipation_prec"]
    error("SimpleODE does not support dissipation stabilization")
  end

  if opts["numBC"] != 0
    error("SimpleODE does not support boundary conditions")
  end

  get!(opts, "calc_jac_explicit", false)
  get!(opts, "preallocate_jacobian_coloring", true)

  return nothing
end


text:
using Distributions
using Gadfly
using HypothesisTests
using Random
using RecipesBase
using Statistics
using Test
using ProjectManagement

Random.seed!(1)

@testset "ProjectManagement.jl" begin
    @testset "$f" for f in (
        "project_examples.jl", # must be first
        "project.jl",
        "timing_distributions.jl",
        "timing_graph.jl",
        "paths.jl",
        "graph_viz.jl",
    )
        include(f)
    end
end


text:
function createcasetable(elementdf, wobrange, inclinationrange, computefunction, pdir)

  # The casetable contains a case for every requested wob and inclination

  casetable = Array{Tuple}(undef, length(wobrange) * length(inclinationrange))
  for (i, wob) in enumerate(wobrange)
    for (j, incl) in enumerate(inclinationrange)

      # Set up directory structure to store results

      wobstring = string(Int(wob))
      inclstring = string(Int(incl))
      cd(pdir) do
        if isdir("tp0/$wobstring/$inclstring")
          rm("tp0/$wobstring/$inclstring", recursive=true)
        end
        (!isdir("tp0")) && mkdir("tp0")
        (!isdir("tp0/$wobstring")) && mkdir("tp0/$wobstring")
        (!isdir("tp0/$wobstring/$inclstring")) && mkdir("tp0/$wobstring/$inclstring")
        if isdir("tp/$wobstring/$inclstring")
          rm("tp/$wobstring/$inclstring", recursive=true)
        end
        (!isdir("tp")) && mkdir("tp")
        (!isdir("tp/$wobstring")) && mkdir("tp/$wobstring")
        (!isdir("tp/$wobstring/$inclstring")) && mkdir("tp/$wobstring/$inclstring")
      end

      # And create the casetable to be used by pmap()

      casetable[ (i-1) * length(inclinationrange) + j ] = (elementdf, wob, incl, computefunction, pdir)

    end
  end

  casetable
end


text:
# -- PRIVATE METHODS BELOW ------------------------------------------------------------------------------ #
# Function to build symbol list -
function build_metabolic_symbol_array(reaction_order_array::Array{String,1}, reaction_dictionary::Dict{String, VLMetabolicReaction})

  # Method variables -
  species_symbol_array::Array{VLSpeciesSymbol,1} = []
  reaction_array = Array{VLMetabolicReaction,1}()
  partitioned_symbol_array = VLSpeciesSymbol[]
  final_partitioned_symbol_array = VLSpeciesSymbol[]

  # build the reaction array -
  for reaction_key in reaction_order_array

      # grab -
      reaction_object = reaction_dictionary[reaction_key]

      # cache -
      push!(reaction_array, reaction_object)
  end

  # Helper function to parse the reaction phrases, split out the symbols
  function _classify_species_type(species_symbol)::Symbol

      # check: gene -
      if (occursin("GENE_",species_symbol) == true)
        if (occursin("OPEN",species_symbol) == true)
          return :UNCLASSIFIED
        else
          return :GENE
        end
      end

      # check: mRNA -
      if (occursin("mRNA_",species_symbol) == true)
          return :MRNA
      end

      # check: protein -
      if (occursin("PROTEIN_",species_symbol) == true)
          return :PROTEIN
      end

      # check: METABOLITE
      if (occursin("M_",species_symbol) == true)
          if (occursin("tRNA", species_symbol) == true)
              return :TRNA
          end

          # else -
          return :METABOLITE
      end

      # check: tRNA -
      if (occursin("tRNA",species_symbol) == true)
          return :TRNA
      end

      # default -
      return :UNCLASSIFIED
  end

  function _parse_phrase(reaction_phrase::String)

    # Method variables -
    local_species_array::Array{VLSpeciesSymbol} = []

    # Split around + -
    fragment_array = split(reaction_phrase,"+")
    for fragment in fragment_array

      if (contains(fragment,"*"))

          local_fragment_array = split(fragment,"*");
          species_symbol = VLSpeciesSymbol();
          species_symbol.species_symbol = local_fragment_array[2];

          # classify -
          species_symbol.species_type = _classify_species_type(species_symbol.species_symbol)
      else

          # Build -
          species_symbol = VLSpeciesSymbol()
          species_symbol.species_symbol = fragment

          # classify -
          species_symbol.species_type = _classify_species_type(species_symbol.species_symbol)
      end

      # grab -
      push!(local_species_array,species_symbol)
    end

    # return -
    return local_species_array
  end

  function _isequal(species_model_1::VLSpeciesSymbol,species_model_2::VLSpeciesSymbol)
    if (species_model_1.species_symbol == species_model_2.species_symbol)
      return true
    end
    return false
  end

  function _add_symbol!(species_symbol_array,species_symbol)

    contains_species_already = false
    for cached_species_model in species_symbol_array

      if (_isequal(cached_species_model,species_symbol))
        contains_species_already = true
        break
      end
    end

    if (contains_species_already == false)
      push!(species_symbol_array,species_symbol)
    end
  end

  # iterate through and get the symbols -
  for reaction in reaction_array

    tmp_species_array_left = _parse_phrase(reaction.left_phrase)
    tmp_species_array_right = _parse_phrase(reaction.right_phrase)
    append!(tmp_species_array_left,tmp_species_array_right)

    for species_model in tmp_species_array_left
      _add_symbol!(species_symbol_array,species_model)
    end
  end

  # ok, so the species symbol array is *not* sorted)
  # let's sort the species -
  tmp_array = String[]
  for species_symbol::VLSpeciesSymbol in species_symbol_array
      push!(tmp_array,species_symbol.species_symbol)
  end

  # generate permutation array -
  idxa_sorted = sortperm(tmp_array)
  sorted_symbol_array = species_symbol_array[idxa_sorted]

  # find all the species of type METABOLITE -
  idx_all_e = findall(x->x.species_type == :METABOLITE, sorted_symbol_array)
  for index in idx_all_e
      push!(partitioned_symbol_array,sorted_symbol_array[index])
  end

  # find all the species of type TRNA -
  idx_all_e = findall(x->x.species_type == :TRNA, sorted_symbol_array)
  for index in idx_all_e
      push!(partitioned_symbol_array,sorted_symbol_array[index])
  end

  # find all the species of type GENE -
  idx_all_e = findall(x->x.species_type == :GENE, sorted_symbol_array)
  for index in idx_all_e
      push!(partitioned_symbol_array,sorted_symbol_array[index])
  end

  # find all the species of type MRNA -
  idx_all_e = findall(x->x.species_type == :MRNA, sorted_symbol_array)
  for index in idx_all_e
      push!(partitioned_symbol_array,sorted_symbol_array[index])
  end

  # find all the species of type PROTEIN -
  idx_all_e = findall(x->x.species_type == :PROTEIN, sorted_symbol_array)
  for index in idx_all_e
      push!(partitioned_symbol_array,sorted_symbol_array[index])
  end

  # find all the species of type UNCLASSIFIED -
  idx_all_e = findall(x->x.species_type == :UNCLASSIFIED, sorted_symbol_array)
  for index in idx_all_e
      push!(partitioned_symbol_array,sorted_symbol_array[index])
  end

  # lastly - filter out [] -
  idx_all_e = findall(x->x.species_symbol != "[]", partitioned_symbol_array)
  for index in idx_all_e
      push!(final_partitioned_symbol_array, partitioned_symbol_array[index])
  end

  # return -
  return final_partitioned_symbol_array
end

# Function to build the stoichiometric_matrix from reaction list -
function build_metabolic_stoichiometric_matrix(species_symbol_array::Array{VLSpeciesSymbol,1}, reaction_array::Array{VLMetabolicReaction,1})

  # Method variables -
  number_of_species = length(species_symbol_array)
  number_of_reactions = length(reaction_array);
  stoichiometric_matrix = zeros(number_of_species,number_of_reactions);

  function _parse_reaction_phrase(lexeme,reaction_phrase)

    # Split around + -
    coefficient = 0.0;
    fragment_array = split(reaction_phrase,"+")
    for fragment in fragment_array

      if (contains(fragment,"*"))
          local_fragment_array = split(fragment,"*");
          test_lexeme = local_fragment_array[2];

          if (lexeme == test_lexeme)
            coefficient = parse(Float64,local_fragment_array[1]);
            break
          end

      else

        # Build -
        test_lexeme = fragment;
        if (lexeme == test_lexeme)
          coefficient = 1.0;
          break
        end
      end
    end

    return coefficient;
  end

  function _find_stoichiometric_coefficient(species_model::VLSpeciesSymbol, reaction::VLMetabolicReaction)

    # Method variables -
    stoichiometric_coefficient = 0.0

    # Check the left and right phrase -
    stoichiometric_coefficient += -1.0*(_parse_reaction_phrase(species_model.species_symbol,reaction.left_phrase))
    stoichiometric_coefficient += _parse_reaction_phrase(species_model.species_symbol,reaction.right_phrase)
    return stoichiometric_coefficient;
  end

  # setup counters -
  for (row_index,species_symbol) in enumerate(species_symbol_array)
    for (col_index,reaction) in enumerate(reaction_array)

      # Is this species involved in this reaction?
      stoichiometric_matrix[row_index,col_index] = _find_stoichiometric_coefficient(species_symbol,reaction);

    end
  end

  # return -
  return stoichiometric_matrix
end
# -- PRIVATE METHODS ABOVE ------------------------------------------------------------------------------ #

# -- PUBLIC METHODS BELOW ------------------------------------------------------------------------------- #
function generate_stoichiomteric_matrix_program_component(problem_dictionary::Dict{String,Any})::VLProgramComponent

  # initialize -
  filename = "Network.dat"
  program_component = VLProgramComponent()
  reaction_array = Array{VLMetabolicReaction,1}()

  # grab the array of reactions -
  reaction_order_array = problem_dictionary["metabolic_reaction_order_array"]
  reaction_dictionary = problem_dictionary["metabolic_reaction_dictionary"]
  for (index, reaction_key) in enumerate(reaction_order_array)

    # check do we have this key?
    if (haskey(reaction_dictionary,reaction_key) == true)

      # grab this reaction object -
      reaction_object = reaction_dictionary[reaction_key]

      # cache -
      push!(reaction_array, reaction_object)
    else
      error("$(reaction_key) not found in the problem dictionary")
    end
  end

  # get the species array -
  list_of_species_objects = problem_dictionary["list_of_species"]

  # build the stoichiometric_matrix -
  stm_buffer = build_metabolic_stoichiometric_matrix(list_of_species_objects, reaction_array)

  # package and return -
  program_component.matrix = stm_buffer
  program_component.filename = filename
  program_component.type = :matrix

  return program_component
end
# -- PUBLIC METHODS ABOVE ------------------------------------------------------------------------------- #

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
export PDSearch

function lcs(X::String, Y::String)
    D = zeros(Int, length(X)+1, length(Y)+1)
    for i = 1:length(X)
        x = X[i]
        for j = 1:length(Y)
            y = Y[j]
            if i == 0 || j == 0
                L[i,j] = 0
            elseif (X[i-1] == Y[j-1])
                L[i,j] = L[i-1,j-1] + 1
            else
                L[i,j] = max(L[i-1,j], L[i,j-1])
            end
        end
    end

    # row 0 and column 0 are initialized to 0 already
    for (i, x) in enumerate(a), (j, y) in enumerate(b)
        if x == y
            lengths[i+1, j+1] = lengths[i, j] + 1
        else
            lengths[i+1, j+1] = max(lengths[i+1, j], lengths[i, j+1])
        end
    end

    result = ""
    x, y = length(a) + 1, length(b) + 1
    while x > 1 && y > 1
        if lengths[x, y] == lengths[x-1, y]
            x -= 1
        elseif lengths[x, y] == lengths[x, y-1]
            y -= 1
        else
            @assert a[x-1] == b[y-1]
            result = string(a[x-1], result)
            x -= 1
            y -= 1
        end
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Xorg_xcb_util_wm_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Xorg_xcb_util_wm")
JLLWrappers.@generate_main_file("Xorg_xcb_util_wm", UUID("c22f9ab0-d5fe-5066-847c-f4bb1cd4e361"))
end  # module Xorg_xcb_util_wm_jll


text:
"""
    Tian()

[Tian binomial model](http://simulations.lpsm.paris/binomial_trees/tian).
"""
struct Tian <: Model end

"""
    evaluate(O, Tian(), N = 1000)

Evaluate option `O` using `Tian` binomial model.

# Arguments
`N`: number of paths to simulate
"""
function evaluate(O::Option, m::Tian, N::Int64 = 1000)
    Δt = O.t / N
    R = exp(O.r * Δt)
    v = exp(O.σ^2 * Δt)
    U = 0.5 * R * v * (v + 1 + sqrt(v^2 + 2 * v - 3))
    D = 0.5 * R * v * (v + 1 - sqrt(v^2 + 2 * v - 3))
    p = (R - D) / (U - D)
    q = (U - R) / (U - D)

    if O.call == -1
        Z = [max(0, O.k - O.s * exp((2 * i - N) * O.σ * √Δt)) for i = 0:N]
    elseif O.call == 1
        Z = [max(0, O.s * exp((2 * i - N) * O.σ * √Δt) - O.k) for i = 0:N]
    end

    for n = N-1:-1:0, i = 0:n
        if O.call == -1
            x = O.k - O.s * exp((2 * i - n) * O.σ * √Δt)
        elseif O.call == 1
            x = O.s * exp((2 * i - n) * O.σ * √Δt) - O.k
        end
        y = (q * Z[i+1] + p * Z[i+2]) / R
        Z[i+1] = max(x, y)
    end

    return Z[1]
end


text:
# INSTRUCTIONS: Start a Julia session on Hoffman and run code below chunk-by-chunk

target_data = ["data1.vcf.gz", "data4.vcf.gz"]
memory = 47
missingprop = 0.1
samples = 1000

# filter/mask/unphase data.
# THIS CHUNK HAS TO RUN BEFORE RUNNING ANY CHUNKS BELOW
cd("/u/home/b/biona001/haplotype_comparisons/data")
for data in target_data
    # if unphased target genotype file exist already, move to next step
    if isfile(data * "./tgt_masked_unphased_" * data)
        continue
    end

    open("filter.sh", "w") do io
        println(io, "#!/bin/bash")
        println(io, "#\$ -cwd")
        println(io, "# error = Merged with joblog")
        println(io, "#\$ -o joblog.\$JOB_ID")
        println(io, "#\$ -j y")
        println(io, "#\$ -l arch=intel-X5650,exclusive,h_rt=24:00:00,h_data=$(memory)G")
        println(io, "# Email address to notify")
        println(io, "#\$ -M \$USER@mail")
        println(io, "# Notify when")
        println(io, "#\$ -m a")
        println(io)
        println(io, "# load the job environment:")
        println(io, ". /u/local/Modules/default/init/modules.sh")
        println(io, "module load julia/1.2.0")
        println(io, "module load R/3.5.1")
        println(io, "module load java/1.8.0_111")
        println(io)
        println(io, "# filter/mask data")
        println(io, "julia ./filter_and_mask.jl $data $samples")
    end

    # submit job
    run(`qsub filter.sh`)
    rm("filter.sh", force=true)
    sleep(2)
end

# AFTER above jobs finish, run prephasing using beagle 4.1.
# THIS CHUNK HAS TO FINISH RUNNING BEFORE MINIMAC4 AND BEAGLE 5 CHUNKS CAN RUN
for data in target_data
    open("prephase.sh", "w") do io
        println(io, "#!/bin/bash")
        println(io, "#\$ -cwd")
        println(io, "# error = Merged with joblog")
        println(io, "#\$ -o joblog.\$JOB_ID")
        println(io, "#\$ -j y")
        println(io, "#\$ -l arch=intel-X5650,exclusive,h_rt=24:00:00,h_data=$(memory)G")
        println(io, "# Email address to notify")
        println(io, "#\$ -M \$USER@mail")
        println(io, "# Notify when")
        println(io, "#\$ -m a")
        println(io)
        println(io, "# load the job environment:")
        println(io, ". /u/local/Modules/default/init/modules.sh")
        println(io, "module load julia/1.2.0")
        println(io, "module load R/3.5.1")
        println(io, "module load java/1.8.0_111")
        println(io)
        println(io, "# run prephasing using beagle 4.1")
        println(io, "java -Xss5m -Xmx$(memory)g -jar beagle4.1.jar gt=./tgt_masked_unphased_$(data) ref=./ref_$(data) niterations=0 out=./tgt_masked_phased_$(data)")
    end
    # submit job
    run(`qsub prephase.sh`)
    rm("prephase.sh", force=true)
    sleep(2)
end

# run MendelImpute (fast version and dynamic programming version)
widths = [400; 800; 1600]
for data in target_data, width in widths
    # fast version
    open("mendel$width.sh", "w") do io
        println(io, "#!/bin/bash")
        println(io, "#\$ -cwd")
        println(io, "# error = Merged with joblog")
        println(io, "#\$ -o joblog.\$JOB_ID")
        println(io, "#\$ -j y")
        println(io, "#\$ -l arch=intel-X5650,exclusive,h_rt=24:00:00,h_data=$(memory)G")
        println(io, "# Email address to notify")
        println(io, "#\$ -M \$USER@mail")
        println(io, "# Notify when")
        println(io, "#\$ -m a")
        println(io)
        println(io, "# load the job environment:")
        println(io, ". /u/local/Modules/default/init/modules.sh")
        println(io, "module load julia/1.2.0")
        println(io, "module load R/3.5.1")
        println(io, "module load java/1.8.0_111")
        println(io)
        println(io, "# run MendelImpute (fast)")
        println(io, "julia mendel_fast.jl $data $width")
    end
    # submit job
    run(`qsub mendel$width.sh`)
    rm("mendel$width.sh", force=true)
    sleep(2)

    # dynamic programming version
    open("dp_mendel$width.sh", "w") do io
        println(io, "#!/bin/bash")
        println(io, "#\$ -cwd")
        println(io, "# error = Merged with joblog")
        println(io, "#\$ -o joblog.\$JOB_ID")
        println(io, "#\$ -j y")
        println(io, "#\$ -l arch=intel-X5650,exclusive,h_rt=24:00:00,h_data=$(memory)G")
        println(io, "# Email address to notify")
        println(io, "#\$ -M \$USER@mail")
        println(io, "# Notify when")
        println(io, "#\$ -m a")
        println(io)
        println(io, "# load the job environment:")
        println(io, ". /u/local/Modules/default/init/modules.sh")
        println(io, "module load julia/1.2.0")
        println(io, "module load R/3.5.1")
        println(io, "module load java/1.8.0_111")
        println(io)
        println(io, "# run MendelImpute (dynamic programming)")
        println(io, "julia mendel_dp.jl $data $width")
    end
    # submit job
    run(`qsub dp_mendel$width.sh`)
    rm("dp_mendel$width.sh", force=true)
    sleep(2)
end

# run minimac4
for data in target_data
    open("minimac4.sh", "w") do io
        println(io, "#!/bin/bash")
        println(io, "#\$ -cwd")
        println(io, "# error = Merged with joblog")
        println(io, "#\$ -o joblog.\$JOB_ID")
        println(io, "#\$ -j y")
        println(io, "#\$ -l arch=intel-X5650,exclusive,h_rt=24:00:00,h_data=$(memory)G")
        println(io, "# Email address to notify")
        println(io, "#\$ -M \$USER@mail")
        println(io, "# Notify when")
        println(io, "#\$ -m a")
        println(io)
        println(io, "# load the job environment:")
        println(io, ". /u/local/Modules/default/init/modules.sh")
        println(io, "module load julia/1.2.0")
        println(io, "module load R/3.5.1")
        println(io, "module load java/1.8.0_111")
        println(io)
        println(io, "# Impute using minimac4")
        println(io, "java -Xmx$(memory)g -jar beagle5.0.jar gt=./tgt_masked_unphased_$(data) ref=./ref_$(data) out=./beagle_imputed_$(data)")
    end
    # submit job
    run(`qsub minimac4.sh`)
    rm("minimac4.sh", force=true)
    sleep(2)
end

# run beagle 5.0
for data in target_data
    open("beagle.sh", "w") do io
        println(io, "#!/bin/bash")
        println(io, "#\$ -cwd")
        println(io, "# error = Merged with joblog")
        println(io, "#\$ -o joblog.\$JOB_ID")
        println(io, "#\$ -j y")
        println(io, "#\$ -l arch=intel-X5650,exclusive,h_rt=24:00:00,h_data=$(memory)G")
        println(io, "# Email address to notify")
        println(io, "#\$ -M \$USER@mail")
        println(io, "# Notify when")
        println(io, "#\$ -m a")
        println(io)
        println(io, "# load the job environment:")
        println(io, ". /u/local/Modules/default/init/modules.sh")
        println(io, "module load julia/1.2.0")
        println(io, "module load R/3.5.1")
        println(io, "module load java/1.8.0_111")
        println(io)
        println(io, "# Impute using beagle 5.0")
        println(io, "java -Xmx$(memory)g -jar beagle5.0.jar gt=./tgt_masked_unphased_$(data) ref=./ref_$(data) out=./beagle_imputed_$(data)")
    end
    # submit job
    run(`qsub beagle.sh`)
    rm("beagle.sh", force=true)
    sleep(2)
end

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
mutable struct Model <: Renderable
	mesh::Mesh
	material::Material
	transform::Matrix{Float32}
	boundDirty::Bool
	bound::Shapes.Shape

	Model(mesh::Mesh, material::Material) = new(mesh, material, Matrix{Float32}(I, 4, 4), true, similar(mesh.bound))
end

isvalid(model::Model) = isvalid(model.mesh) && isvalid(model.material)

gettransform(model::Model) = model.transform
function settransform(model::Model, m::Matrix{Float32})
	@assert isvalid(model)
	model.transform[:] = m
	model.boundDirty = true
end

localbound(model::Model) = model.mesh.bound

function getbound(model::Model)
	@assert isvalid(model)
	if model.boundDirty
		Shapes.transform(model.bound, model.transform, model.mesh.bound)
		model.boundDirty = false
	end
	model.bound
end

function setbound(model::Model, bound::Shapes.Shape)
	@assert isvalid(model)
	model.bound = bound
	model.boundDirty = false
end

function render(model::Model, renderer::Renderer)
	@assert isvalid(model)
	model.material.shader.setupMaterial(model)
	apply(model.material, renderer)
	render(model.mesh)
end


text:
function adjusted_sharpe_ratio{T}(Ra::TimeArray{T,1})
function appraisal_ratio{T}(Ra::TimeArray{T,1}, Rb::TimeArray{T,1})
function bernardo_ledoit_ratio{T}(Ra::TimeArray{T,1})
function calmar_ratio{T}(Ra::TimeArray{T,1})
function d_ratio{T}(Ra::TimeArray{T,1})
function kelly_ratio{T}(Ra::TimeArray{T,1})
function martin_ratio{T}(Ra::TimeArray{T,1})
function omega_sharpe_ratio{T}(Ra::TimeArray{T,1})
function pain_ratio{T}(Ra::TimeArray{T,1})
function prospect_ratio{T}(Ra::TimeArray{T,1}, mar=0)
function sharpe_ratio{T}(Ra::TimeArray{T,1}, method::Function=std)
function sharpe_ratio{T}(Ra::TimeArray{T,1}, method::Function=value_at_risk)
function sharpe_ratio{T}(Ra::TimeArray{T,1}, method::Function=expected_shortfall)
function skewness_kurtosis_ratio{T}(Ra::TimeArray{T,1})
function sortino_ratio{T}(Ra::TimeArray{T,1})
function treynor_ratio{T}(Ra::TimeArray{T,1},Rb::TimeArray{T,1})
function upside_potential_ratio{T}(Ra::TimeArray{T,1})
function hurst_index{T}(Ra::TimeArray{T,1})
function pain_index{T}(Ra::TimeArray{T,1})
function smoothing_index{T}(Ra::TimeArray{T,1})
function ulcer_index{T}(Ra::TimeArray{T,1})
function active_premium{T}(Ra::TimeArray{T,1},Rb::TimeArray{T,1})
function downside_deviation{T}(Ra::TimeArray{T,1})
function downside_frequency{T}(Ra::TimeArray{T,1})
function expected_shortfall{T}(Ra::TimeArray{T,1})
function fama_beta{T}(Ra::TimeArray{T,1},Rb::TimeArray{T,1})
function lower_partial_moment{T}(Ra::TimeArray{T,1},Rb::TimeArray{T,1})
function m2_sortina{T}(Ra::TimeArray{T,1},Rb::TimeArray{T,1})
function msquared{T}(Ra::TimeArray{T,1},Rb::TimeArray{T,1})
function msquared_excess{T}(Ra::TimeArray{T,1},Rb::TimeArray{T,1})
function modigliani{T}(Ra::TimeArray{T,1},Rb::TimeArray{T,1})
function max_drawdown{T}(Ra::TimeArray{T,1})
function mean_absolute_deviation{T}(Ra::TimeArray{T,1})
function net_selectivity{T}(Ra::TimeArray{T,1},Rb::TimeArray{T,1})
function omega{T}(Ra::TimeArray{T,1})
function selectivity{T}(Ra::TimeArray{T,1},Rb::TimeArray{T,1})
function semi_deviation{T}(Ra::TimeArray{T,1})
function specific_risk{T}(Ra::TimeArray{T,1},Rb::TimeArray{T,1})
function total_risk{T}(Ra::TimeArray{T,1},Rb::TimeArray{T,1})
function upside_frequency{T}(Ra::TimeArray{T,1})
function upside_risk{T}(Ra::TimeArray{T,1})
function value_at_risk{T}(Ra::TimeArray{T,1})
function volatility_skewness{T}(Ra::TimeArray{T,1})


text:
using PowerSpectra
using Healpix
using CSV
using Test
using LinearAlgebra
using DataFrames
using DelimitedFiles
using NPZ

@testset "Basic Covmat Mode Decoupling" begin
    A = [1.0 0.2 0.3; 0.4 2.0 0.15; 0.3 0.1 1.44]
    B1 = [1.2 0.6 0.1; 0.3 1.4 0.5; 0.44 0.2 1.3]
    B2 = [1.6 0.4 0.1; 0.3 1.4 0.9; 0.45 0.8 1.7]
    Cref = inv(B1) * (A) * (inv(B2)')
    C = decouple_covmat(SpectralArray(deepcopy(A)), SpectralArray(B1), SpectralArray(B2))
    @test all(Cref .≈ parent(C))

    # A₀ = deepcopy(A)
    # A₀[1,1] = 0.0
    # Cref = inv(B1) * A₀ * (inv(B2)')
    # C = decouple_covmat(SpectralArray(deepcopy(A)), SpectralArray(B1), SpectralArray(B2);
    #     lmin1=1, lmin2=1)
    # @test all(Cref .≈ parent(C))
end

##
@testset "Covariance Matrix Diagonal in the Isotropic Noise Limit" begin
    nside = 256
    mask1_T = readMapFromFITS("data/mask1_T.fits", 1, Float64)
    mask2_T = readMapFromFITS("data/mask2_T.fits", 1, Float64)
    mask1_P = readMapFromFITS("data/mask1_P.fits", 1, Float64)
    mask2_P = readMapFromFITS("data/mask2_P.fits", 1, Float64)
    unit_var = HealpixMap{Float64, RingOrder}(ones(nside2npix(nside)))
    flat_mask = HealpixMap{Float64, RingOrder}(ones(nside2npix(nside)) )
    beam1 = SpectralVector(ones(3nside))
    beam2 = SpectralVector(ones(3nside))
    theory = CSV.read("data/theory.csv", DataFrame)
    noise = CSV.read("data/noise.csv", DataFrame)
    identity_spectrum = SpectralVector(ones(3nside));

    cltt = SpectralVector(convert(Vector, theory.cltt))
    clte = SpectralVector(convert(Vector, theory.clte))
    clee = SpectralVector(convert(Vector, theory.clee))
    nlee = SpectralVector(convert(Vector, noise.nlee))
    nltt = SpectralVector(convert(Vector, noise.nltt))

    # this test specifies a map with unit variance. the corresponding white noise level is divided out in r_coeff
    Ω_p = PowerSpectra.pixsize(mask1_T)
    N_white = Ω_p
    r_coeff = Dict{PowerSpectra.SpectrumName, SpectralVector{Float64, Vector{Float64}}}(
        (:TT, "143_hm1", "143_hm1") => sqrt.(nltt ./ N_white),
        (:TT, "143_hm1", "143_hm2") => identity_spectrum,
        (:TT, "143_hm2", "143_hm1") => identity_spectrum,
        (:TT, "143_hm2", "143_hm2") => sqrt.(nltt ./ N_white),

        (:EE, "143_hm1", "143_hm1") => sqrt.(nlee ./ N_white),
        (:EE, "143_hm1", "143_hm2") => identity_spectrum,
        (:EE, "143_hm2", "143_hm1") => identity_spectrum,
        (:EE, "143_hm2", "143_hm2") => sqrt.(nlee ./ N_white))

    spectra = Dict{PowerSpectra.SpectrumName, SpectralVector{Float64, Vector{Float64}}}(
        (:TT, "143_hm1", "143_hm1") => cltt,
        (:TT, "143_hm1", "143_hm2") => cltt,
        (:TT, "143_hm2", "143_hm1") => cltt,
        (:TT, "143_hm2", "143_hm2") => cltt,

        (:EE, "143_hm1", "143_hm1") => clee,
        (:EE, "143_hm1", "143_hm2") => clee,
        (:EE, "143_hm2", "143_hm1") => clee,
        (:EE, "143_hm2", "143_hm2") => clee ,

        (:TE, "143_hm1", "143_hm1") => clte,
        (:TE, "143_hm1", "143_hm2") => clte,
        (:TE, "143_hm2", "143_hm1") => clte,
        (:TE, "143_hm2", "143_hm2") => clte,
    )

    σ² = PolarizedHealpixMap(unit_var, unit_var, unit_var)
    m1 = CovField("143_hm1", mask1_T, mask1_P, σ², beam1, beam1)
    m2 = CovField("143_hm2", mask2_T, mask2_P, σ², beam2, beam2)
    workspace = CovarianceWorkspace(m1, m2, m1, m2)

    C = coupledcov(:TT, :TT, workspace, spectra, r_coeff)
    C_ref = npzread("data/covar_TT_TT.npy")
    @test isapprox(diag(parent(C))[3:end], diag(C_ref)[3:end])

    C = coupledcov(:TT, :TE, workspace, spectra, r_coeff)
    C_ref = npzread("data/covar_TT_TE.npy")
    @test isapprox(diag(parent(C))[3:end], diag(C_ref)[3:end])

    C = coupledcov(:TE, :TE, workspace, spectra, r_coeff)
    C_ref = npzread("data/covar_TE_TE.npy")
    @test isapprox(diag(parent(C))[3:end], diag(C_ref)[3:end])

    C = coupledcov(:TT, :EE, workspace, spectra, r_coeff)
    C_ref = npzread("data/covar_TT_EE.npy")
    @test isapprox(diag(parent(C))[3:end], diag(C_ref)[3:end])

    PowerSpectra.coupledcovTEEE!(C, workspace, spectra, r_coeff; planck=false);
    C_ref = npzread("data/covar_TE_EE.npy")
    @test isapprox(diag(parent(C))[3:end], diag(C_ref)[3:end])

    C = coupledcov(:EE, :EE, workspace, spectra, r_coeff)
    C_ref = npzread("data/covar_EE_EE.npy")
    @test isapprox(diag(parent(C))[3:end], diag(C_ref)[3:end])

    # test that planck approx is kind of close at high ell
    C = coupledcov(:TE, :EE, workspace, spectra, r_coeff)
    C_ref = npzread("data/covar_TE_EE.npy")
    @test isapprox(diag(parent(C))[30:end], diag(C_ref)[30:end], rtol=0.01)
end

##
@testset "Covariance Matrix Decoupling" begin
    mask1_T = readMapFromFITS("data/mask1_T.fits", 1, Float64)
    mask2_T = readMapFromFITS("data/mask2_T.fits", 1, Float64)
    # test decoupling
    𝐌 = mcm(:TT, mask1_T, mask2_T)
    C_coupled = similar(𝐌)
    X=rand(size(parent(𝐌))...)
    C_coupled.parent.parent .= X'*X
    C_decoupled = decouple_covmat(C_coupled, 𝐌, 𝐌)

    @test isapprox(parent(C_coupled),
        parent(𝐌) * parent(C_decoupled) * parent(𝐌)' )
end


text:
function model_step!(model::Agents.ABM)
    model.new_cases_real = 0
    model.new_cases_detected = 0
    update_seir_states!(model)
    start_quarantines!(model)
    end_quarantines!(model)
    update_statecounts!(model)
    if model.behavior_activated
        update_authority_recommendation!(model)
        update_authority_policy!(model)
    end
    model.clock += 1
    return model
end

function update_seir_states!(model::Agents.ABM)
    for a in values(model.agents)
        advance_state_internal!(a, model)
    end
    for n in LightGraphs.vertices(model.space.graph)
        simulate_node_encounters!(model, n)
    end
    return model
end

function start_quarantines!(model::Agents.ABM)
    for a in values(model.agents)
        if a.quarantined & !(a.pos == a.home)
            Agents.move_agent!(a, a.home, model)
        end
    end
    return model
end

function end_quarantines!(model::Agents.ABM)
    for agent in values(model.agents)
        if agent.quarantined & (agent.ticks_quarantined >= model.quarantine_duration)
            agent.quarantined = false
        end
    end
    return model
end

function update_statecounts!(model::Agents.ABM)
    model.Id_cumulative = model.Id_cumulative + model.new_cases_detected
    model.I_cumulative = model.I_cumulative + model.new_cases_real
    model.S_count = sum([model.agents[i].state == :S for i in 1:model.agentcount])
    model.E_count = sum([model.agents[i].state == :E for i in 1:model.agentcount])
    model.Iu_count = sum([model.agents[i].state == :Iu for i in 1:model.agentcount])
    model.Id_count = sum([model.agents[i].state == :Id for i in 1:model.agentcount])
    model.I_count = model.E_count + model.Iu_count + model.Id_count
    model.R_count = sum([model.agents[i].state == :R for i in 1:model.agentcount])
    model.B_count = sum([model.agents[i].behavior for i in 1:model.agentcount])
    return model
end

function update_authority_recommendation!(model::Agents.ABM)
    if model.clock >= model.tick_authority_recommendation
        model.authority_recommendation = 1.0
    else
        model.authority_recommendation = 0.0
    end
    return model
end

function update_authority_policy!(model::Agents.ABM)
    if model.clock >= model.tick_authority_policy
        model.behavior_policy = true
    else
        model.behavior_policy = false
    end
    return model
end


text:
"""
	parse_sfs(;sample_size::Int64,data::String)

Function to parse polymorphism and divergence by subset of genes. The input data is based on supplementary material described at [Uricchio et al. 2019](https://doi.org/10.1038/s41559-019-0890-6). Please be sure the file is tabulated.

| GeneId | Pn | DAF seppareted by commas | Ps | DAF separated by commas | Dn | Ds |
|--------|----|--------------------------|----|-------------------------|----|----|
| XXXX   | 0  | ,0,0,0,0,0,0,0,0         | 0  | ,0,0,0,0,0,0,0,0        | 0  | 0  |
| XXXX   | 0  | ,0,0,0,0,0,0,0,0         | 0  | ,0,0,0,0,0,0,0,0        | 0  | 0  |
| XXXX   | 0  | ,0,0,0,0,0,0,0,0         | 0  | ,0,0,0,0,0,0,0,0        | 0  | 0  |

# Arguments
 - `sample_size::Int64`: data sample size.
 - `data`: File containing polymorphic and divergence data
 - `gene_list`: File containing gene IDs to subset. You can perform multiple subset input a file where eachrow contain the gene IDs to subset.
 - `sfs_columns::Array{Int64,1}`: non-synonymous and synonymous daf columns. Please introduce first the non-synonymous number.
 - `div_columns::Array{Int64,1}`: non-synonymous and synonymous divergence columns. Please introduce first the non-synonymous number.
 - `bins::Array{Int64,1}`: bin to collapse the SFS from `sample_size` to new `bin` size

# Returns
 - `Array{Float64,1}`: α values
 - `Array{Float64,2}`: Site Frequency Spectrum
 - `Array{Float64,1}`: Synonymous and non-synonymous divergence counts
"""
function parse_sfs(;param::parameters,data::S,gene_list::Union{Nothing,S}=nothing,sfs_columns::Vector{Int64}=[3,5],div_columns::Vector{Int64}=[6,7],m_columns::Vector{Int64}=[8,9],bins::Union{Nothing,Int64}=nothing,isolines::Bool=false) where S <: AbstractString

	@unpack n, cutoff = param

	if isolines
		s_size = n
	else
		s_size = (n*2)
	end

	df = CSV.read(data,header=false,delim='\t',DataFrame)

	if(!isnothing(gene_list))
		ids = @view df[:,1];

		gene_matrix = Array(CSV.read(gene_list,header=false,DataFrame))
		m,n = size(gene_matrix)

		if n == 1
			gene_matrix = permutedims(gene_matrix)
		end

		out = SubDataFrame[]
		for c ∈ eachrow(gene_matrix)
			tmp = @view df[filter(!isnothing,indexin(c,ids)),:];
			push!(out,tmp);
		end

		α, sfs, divergence, m = unzip(map(i-> get_pol_div(i,s_size,cutoff,sfs_columns,div_columns,m_columns,bins),out));

	else
		α, sfs, divergence,m = get_pol_div(df,s_size,cutoff,sfs_columns,div_columns,m_columns,bins);
		α = [α]; sfs = [sfs]; divergence = [divergence]; m = [m]
	end

	return α, sfs, divergence, m
end

function get_pol_div(df_subset::Union{DataFrame,SubDataFrame},s_size::Int64,cutoff::Vector{Float64},sfs_columns::Vector{Int64},div_columns::Vector{Int64},m_columns::Vector{Int64},bins::Union{Nothing,Int64})

	g(x) = parse.(Float64,x[2:end-1])

	freq = OrderedDict(round.(collect(1:(s_size-1))/s_size,digits=4) .=> 0)
	# freq = OrderedDict(collect(1:(s_size-1)) .=> 0)

	tmp  = split.(df_subset[:,sfs_columns], ",")

	pn   = sort!(OrderedDict(countmap(round.(reduce(vcat,tmp[:,1] .|> g),digits=4))))
	ps   = sort!(OrderedDict(countmap(round.(reduce(vcat,tmp[:,2] .|> g),digits=4))))

	# Dn, Ds, Pn, Ps, sfs
    Dn     = sum(df_subset[:,div_columns[1]])
    Ds     = sum(df_subset[:,div_columns[2]])
    Pn     = sum(values(pn))
    Ps     = sum(values(ps))
    sfs_pn = reduce(vcat,values(merge(+,freq,pn)))
    sfs_ps = reduce(vcat,values(merge(+,freq,ps)))

	if(!isnothing(bins))
		sfs_pn = reduce_sfs(hcat(collect(1:(s_size-1)),sfs_pn),bins)[:,2]
		sfs_ps = reduce_sfs(hcat(collect(1:(s_size-1)),sfs_ps),bins)[:,2]

		sfs   = hcat(collect(1:(bins-1))./bins,sfs_pn,sfs_ps,1:(bins-1))
	else
		sfs   = hcat(freq.keys,sfs_pn,sfs_ps,(1:s_size-1))
	end

	# Filtering SFS and changing frequency to DAC
	sfs = sfs[sfs[:,1] .>= cutoff[1].&& sfs[:,1] .<= cutoff[2],[4,2,3]]

	scumu = cumulative_sfs(sfs)
	α     = round.(1 .- (Ds/Dn .*  scumu[:,2] ./scumu[:,3]),digits=5)

	m::Matrix{Int64} = try
		[sum(df_subset[:,m_columns[1]]) sum(df_subset[:,m_columns[2]])]
	catch
		m = [0 0]
	end

	return (α,sfs,[Dn Ds],m)
end

"""
	ABCreg(analysis_folder, replicas, S, tol, abcreg)

Performing ABC inference using ABCreg. Please, be sure your analysis_folder contain the files alphas.txt and summaries.txt produced by Analytical.summaryStatsFromRates()

# Arguments
 - `analysis_folder::String` : Folder containing the observed data and summary estatistics. It will be used to output the posterior distributions
 - `P::Int64` : Number of parameters to perform the inference.
 - `S::Int64` : Number of summary stastitics to perform the inference.
 - `tol::Float64` : Tolerance value. It define the number of accepted values at ABC inference
 - `abcreg::String` : Path to ABCreg binary

# Output
Files containing posterior distributions from ABCreg

"""
function ABCreg(;analysis_folder::String,S::Int64,P::Int64=5,tol::Float64,abcreg::String)

	# List alphas and summstat files
	a_file     = filter(x -> occursin("alphas",x), readdir(analysis_folder,join=true));
	sum_file   = filter(x -> occursin("summstat",x), readdir(analysis_folder,join=true));

	# Creating output names
	# out = [analysis_folder .* "/out"]
	out = analysis_folder .* "/out_" .* string.(1:size(a_file,1))

	r(a,s,o,abcreg=abcreg,P=P,S=S,tol=tol) = run(`$abcreg -d $a -p $s -P $P -S $S -t $tol -b $o`)

	progress_pmap(r,a_file,sum_file,out);

end

"""
	Bootstrap data following polyDFE manual
"""
function data_to_poisson(sfs::Vector{Matrix{Float64}},divergence::Vector{Matrix{Int64}},dac::Vector{Int64},bootstrap::Union{Bool,Int64})

	if bootstrap != false
		sfs         = repeat(sfs,bootstrap)
		divergence  = repeat(divergence,bootstrap)
		pr(x)       = hcat(x[:,1],pois_rand.(x[:,2:end]))
		sfs[2:end] .= pr.(sfs[2:end])
	end

    scumu      = cumulative_sfs.(sfs)

    f(x,d=dac) = map(z->sum(x[x[:,1].==z,2:3]),d)
    al(a,b)    = hcat(a[:,1],@. round(1 - (b[2]/b[1] * a[:,2]/a[:,3]),digits=5))

    s          = f.(scumu)
    d          = [[sum(divergence[i][1:2])] for i ∈ eachindex(divergence)]

    α_x        = al.(scumu,divergence)
    α          = map(x -> permutedims(x[in(dac).(x[:,1]),2]),α_x)

	return(α,s,d)
end

function open_sfs_div(x::Array{String,1},y::Array{String,1},dac::Vector{Int64},bootstrap::Union{Bool,Int64})

	sfs = Array.(CSV.read.(x,DataFrame,header=false))
	divergence = Array.(CSV.read.(y,DataFrame,header=false))

	if bootstrap != false
		sfs = repeat(sfs,bootstrap)
		divergence = repeat(divergence,bootstrap)
		pr(x) = hcat(x[:,1],pois_rand.(x[:,2:end]))
		sfs[2:end] .= pr.(sfs[2:end])
	end

	scumu = cumulative_sfs.(sfs)
	f(x,d=dac) = sum(x[:,2:3],dims=2)[d]
	s = f.(scumu)

	d = [[sum(divergence[i][1:2])] for i ∈ eachindex(divergence)]
	al(a,b,c=dac) = @. round(1 - (b[2]/b[1] * a[:,2]/a[:,3])[c],digits=5)
	α = permutedims.(al.(scumu,divergence))
	return(s,d,α)
end

"""
	Function to download and source plotMap function. We do not include at part of the module to avoid external dependecies. Once the function is execute properly you will have a function called *plotMap which used R to
		estimate and plot the Maximum A Posterior following ABCreg example. It uses locfit and ggplot2 libraries.
#Arguments
- `script_path::String`: Path and name to save MAP scritp.
"""
function source_plot_map_r(script_path::String)

	try
		download("https://raw.githubusercontent.com/jmurga/MKtest.jl/main/scripts/plot_map.jl",script_path)
		include(script_path)
	catch
		println("\nPlease be sure you have R and the pacakges data.table, ggplot2 and locfit installed in your system. You can install them using Conda inside julia with the following commands:\n\n\tusing Pkg\n\tPkg.add(\"conda\")\n\tENV[\"R_HOME\"]=\"*\"\n\tPkg.add(\"Conda\")\n\tusing Conda\n\tConda.add(\"r-base\",channel=\"conda-forge\")\n\tConda.add([\"r-locfit\",\"r-ggplot2\",\"r-data.table\",\"r-r.utils\"],channel=\"conda-forge\")\n\tPkg.add(\"RCall\")\n")
	end
end

function grapes(sfs::Vector{Matrix{Float64}},divergence::Vector{Matrix{Int64}},m::Vector{Matrix{Int64}},model::String,folder::String,bins::Int64)

	if(isnothing(CondaPkg.which("grapes")))
		CondaPkg.add("grapes-static",channel="genomedk")
	end

	grapes = CondaPkg.which("grapes")

    sfs = reduce_sfs.(sfs,bins);

    pn  = map(x-> permutedims(x[:,2]),sfs);
    ps  = map(x-> permutedims(x[:,3]),sfs);

    dn  = map(x->x[1],divergence);
    ds  = map(x->x[2],divergence);

    mn  = map(x->x[1],m);
    ms  = map(x->x[2],m);

    idx = string.(collect(1:length(sfs)));

	f(pn,ps,dn,ds,mn,ms,w) = DataFrame(hcat("dofe_"*string(w),bins,mn,pn,ms,ps,mn,dn,ms,ds...),:auto)

    dofe          = f.(pn,ps,dn,ds,mn,ms,idx);
    h             = fill(DataFrame(["" ""; "#unfolded" ""],:auto),length(sfs))
    output_dofe   = folder .* "/dofe_" .* idx .* ".txt"
    output_grapes = folder .* "/dofe_" .* idx .* "." .* model

	w(x,name,a=false) = CSV.write(name,x,delim='\t',header=false,append=a);

	w.(h,output_dofe)
	w.(dofe,output_dofe,fill(true,length(sfs)))

	r(d,o,m=model,gr=grapes) = run(`$gr -in $d -out $o -model $m`)

	@suppress_out begin
		progress_pmap(r,output_dofe,output_grapes);
	end

	df = @suppress begin
		CSV.read.(output_grapes,DataFrame,footerskip=1,skipto=3);
	end

	return(vcat(df...))
end


text:
pushfirst!(LOAD_PATH, joinpath(@__DIR__, "../.."))

using OceanTurbulenceParameterEstimation, LinearAlgebra, CairoMakie
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity
using OceanTurbulenceParameterEstimation.EnsembleKalmanInversions: NaNResampler, FullEnsembleDistribution

include("utils/lesbrary_paths.jl")
include("utils/parameters.jl")
include("utils/visualize_profile_predictions.jl")

examples_path = joinpath(pathof(OceanTurbulenceParameterEstimation), "../../examples")
include(joinpath(examples_path, "intro_to_inverse_problems.jl"))

# Pick a parameter set defined in `./parameters.jl`
parameter_set = CATKEParametersRiDependent
closure = closure_with_parameters(CATKEVerticalDiffusivity(Float64;), parameter_set.settings)

# Pick the secret `true_parameters`
true_parameters = (Cᵟu = 0.5, CᴷRiʷ = 1.0, Cᵂu★ = 2.0, CᵂwΔ = 1.0, Cᴷeʳ = 5.0, Cᵟc = 0.5, Cᴰ = 2.0, Cᴷc⁻ = 0.5, Cᴷe⁻ = 0.2, Cᴷcʳ = 3.0, Cᴸᵇ = 1.0, CᴷRiᶜ = 1.0, Cᴷuʳ = 4.0, Cᴷu⁻ = 0.8, Cᵟe = 0.5)
true_closure = closure_with_parameters(closure, true_parameters)

# Generate and load synthetic observations
kwargs = (tracers = (:b, :e), Δt = 10.0, stop_time = 4hours)

# Note: if an output file of the same name already exist, `generate_synthetic_observations` will return the existing path and skip re-generating the data.
observ_path1 = generate_synthetic_observations("perfect_model_observation1"; Qᵘ = 2e-4, Qᵇ = 4e-8, f₀ = 1e-4, closure = true_closure, kwargs...)
observ_path2 = generate_synthetic_observations("perfect_model_observation2"; Qᵘ = -1e-4, Qᵇ = 0, f₀ = 0, closure = true_closure, kwargs...)
observation1 = SyntheticObservations(observ_path1, field_names = (:b, :e, :u, :v), normalize = ZScore)
observation2 = SyntheticObservations(observ_path2, field_names = (:b, :e, :u, :v), normalize = ZScore)
observations = [observation1, observation2]

# Set up ensemble model
ensemble_simulation, closure = build_ensemble_simulation(observations; Nensemble = 50)

# Build free parameters
build_prior(name) = ConstrainedNormal(0.0, 1.0, bounds(name) .* 0.5...)
free_parameters = FreeParameters(named_tuple_map(names(parameter_set), build_prior))

# Pack everything into Inverse Problem `calibration`
calibration = InverseProblem(observations, ensemble_simulation, free_parameters, output_map = ConcatenatedOutputMap());

# Make sure the forward map evaluated on the true parameters matches the observation map
x = forward_map(calibration, true_parameters)[:, 1:1];
y = observation_map(calibration);
@show x == y

directory = "calibrate_catke_to_perfect_model/"
isdir(directory) || mkpath(directory)

visualize!(calibration, true_parameters;
    field_names = (:u, :v, :b, :e),
    directory = directory,
    filename = "perfect_model_visual_true_params.png"
)

# Calibrate

noise_covariance = Matrix(Diagonal(observation_map_variance_across_time(calibration)[1, :, 1] .+ 1e-5))
resampler = NaNResampler(; abort_fraction=0.5, distribution=FullEnsembleDistribution())

eki = EnsembleKalmanInversion(calibration; noise_covariance = noise_covariance,
                                           resampler = resampler)
params = iterate!(eki; iterations = 5)

###
### Summary Plots
###

using CairoMakie
using LinearAlgebra

visualize!(calibration, params;
    field_names = [:u, :v, :b, :e],
    directory = directory,
    filename = "perfect_model_visual_calibrated.png"
)
@show params

# Last, we visualize the outputs of EKI calibration.

### Parameter convergence plot

# Vector of NamedTuples, ensemble mean at each iteration
ensemble_means = getproperty.(eki.iteration_summaries, :ensemble_mean)

# N_param x N_iter matrix, ensemble covariance at each iteration
θθ_std_arr = sqrt.(hcat(diag.(getproperty.(eki.iteration_summaries, :ensemble_cov))...))

N_param, N_iter = size(θθ_std_arr)
iter_range = 0:(N_iter-1)
pnames = calibration.free_parameters.names

n_cols = 3
n_rows = Int(ceil(N_param / n_cols))
ax_coords = [(i, j) for i = 1:n_rows, j = 1:n_cols]

f = Figure(resolution = (500n_cols, 200n_rows))
for (i, pname) in enumerate(pnames)
    coords = ax_coords[i]
    ax = Axis(f[coords...],
        xlabel = "Iteration",
        xticks = iter_range,
        ylabel = string(pname))

    ax.ylabelsize = 20

    mean_values = [getproperty.(ensemble_means, pname)...]
    lines!(ax, iter_range, mean_values)
    band!(ax, iter_range, mean_values .+ θθ_std_arr[i, :], mean_values .- θθ_std_arr[i, :])
    hlines!(ax, [true_parameters[pname]], color = :red)
end
save(joinpath(directory, "perfect_model_calibration_parameter_convergence.png"), f);

### Pairwise ensemble plots

N_param, N_iter = size(θθ_std_arr)
for pname1 in pnames, pname2 in pnames
    if pname1 != pname2

        f = Figure()
        axtop = Axis(f[1, 1])
        axmain = Axis(f[2, 1],
            xlabel = string(pname1),
            ylabel = string(pname2)
        )
        axright = Axis(f[2, 2])
        scatters = []
        for iteration in [0, 1, 2, N_iter - 1]
            ensemble = eki.iteration_summaries[iteration].parameters
            ensemble = [[particle[pname1], particle[pname2]] for particle in ensemble]
            ensemble = transpose(hcat(ensemble...)) # N_ensemble x 2
            push!(scatters, scatter!(axmain, ensemble))
            density!(axtop, ensemble[:, 1])
            density!(axright, ensemble[:, 2], direction = :y)
        end
        vlines!(axmain, [true_parameters[pname1]], color = :red)
        vlines!(axtop, [true_parameters[pname1]], color = :red)
        hlines!(axmain, [true_parameters[pname2]], color = :red, alpha = 0.6)
        hlines!(axright, [true_parameters[pname2]], color = :red, alpha = 0.6)
        colsize!(f.layout, 1, Fixed(300))
        colsize!(f.layout, 2, Fixed(200))
        rowsize!(f.layout, 1, Fixed(200))
        rowsize!(f.layout, 2, Fixed(300))
        Legend(f[1, 2], scatters,
            ["Initial ensemble", "Iteration 1", "Iteration 2", "Iteration $N_iter"],
            position = :lb)
        hidedecorations!(axtop, grid = false)
        hidedecorations!(axright, grid = false)
        xlims!(axright, 0, 10)
        ylims!(axtop, 0, 10)
        save(joinpath(directory, "eki_$(pname1)_$(pname2).png"), f)
    end
end

# Compare EKI result to true values

weight_distances = [norm(collect(ensemble_means[iter]) .- collect(true_parameters)) for iter = 0:N_iter-1]
output_distances = [mapslices(norm, (forward_map(calibration, [ensemble_means...])[:, 1:N_iter] .- y), dims = 1)...]

f = Figure()
lines(f[1, 1], iter_range, weight_distances, color = :red, linewidth = 2,
    axis = (title = "Parameter distance",
        xlabel = "Iteration",
        ylabel = "|θ̅ₙ - θ⋆|",
        yscale = log10))
lines(f[1, 2], iter_range, output_distances, color = :blue, linewidth = 2,
    axis = (title = "Output distance",
        xlabel = "Iteration",
        ylabel = "|G(θ̅ₙ) - y|",
        yscale = log10))

save(joinpath(directory, "error_convergence_summary.png"), f);


text:
# Representation of SQL entities.

"""
    SQLTable(; schema = nothing, name, columns)
    SQLTable(name; schema = nothing, columns)
    SQLTable(name, columns...; schema = nothing)

The structure of a SQL table or a table-like entity (TEMP TABLE, VIEW, etc) for
use as a reference in assembling SQL queries.

The `SQLTable` constructor expects the table `name`, a vector `columns` of
column names, and, optionally, the name of the table `schema`.  A name can be
provided as a `Symbol` or `String` value.

# Examples

```jldoctest
julia> t = SQLTable(:location,
                    :location_id, :address_1, :address_2, :city, :state, :zip);

julia> show(t.name)
:location

julia> show(t.columns)
[:location_id, :address_1, :address_2, :city, :state, :zip]
```

```jldoctest
julia> t = SQLTable(schema = "public",
                    name = "person",
                    columns = ["person_id", "birth_datetime", "location_id"]);

julia> show(t.schema)
:public

julia> show(t.name)
:person

julia> show(t.columns)
[:person_id, :birth_datetime, :location_id]
```
"""
struct SQLTable
    schema::Union{Symbol, Nothing}
    name::Symbol
    columns::Vector{Symbol}
    column_set::Set{Symbol}

    function SQLTable(;
                      schema::Union{Symbol, AbstractString, Nothing} = nothing,
                      name::Union{Symbol, AbstractString},
                      columns::AbstractVector{<:Union{Symbol, AbstractString}})
        schema = schema !== nothing ? Symbol(schema) : nothing
        name = Symbol(name)
        columns =
            !isa(columns, Vector{Symbol}) ?
                Symbol[Symbol(col) for col in columns] :
                columns
        column_set = Set{Symbol}(columns)
        new(schema, name, columns, column_set)
    end
end

SQLTable(name; schema = nothing, columns) =
    SQLTable(schema = schema, name = name, columns = columns)

SQLTable(name, columns...; schema = nothing) =
    SQLTable(schema = schema, name = name, columns = [columns...])

Base.show(io::IO, tbl::SQLTable) =
    print(io, quoteof(tbl, limit = true))

Base.show(io::IO, ::MIME"text/plain", tbl::SQLTable) =
    pprint(io, tbl)

function PrettyPrinting.quoteof(tbl::SQLTable; limit::Bool = false)
    ex = Expr(:call, nameof(SQLTable))
    push!(ex.args, quoteof(tbl.name))
    if tbl.schema !== nothing
        push!(ex.args, Expr(:kw, :schema, quoteof(tbl.schema)))
    end
    if !limit
        push!(ex.args, Expr(:kw, :columns, tbl.columns))
    else
        push!(ex.args, :…)
    end
    ex
end


text:
using Keldysh, Test

@testset "dos" begin
  dos_integrator = Keldysh.GaussKronrodDOSIntegrator()
  # Flat DOS
  let ν=5.0, D=2.0, dos = Keldysh.flat_dos(ν=ν, D=D)
    α = ν*D
    moments = [dos_integrator(ω -> ω^n, dos) for n = 0:3]
    moments_ref = (1/π)*[D*(1 + coth(α)),
                         0.,
                         D^3*(π^2 + α^2)*(1 + coth(α))/(3α^2),
                         0.]
    @test isapprox(moments, moments_ref, atol=1e-10, rtol=1e-10)
  end

  # Gaussian DOS
  let ν = 2.0, ϵ = 1.0, dos = Keldysh.gaussian_dos(ν=ν, ϵ=ϵ)
    moments = [dos_integrator(ω -> ω^n, dos) for n = 0:6]
    moments_ref = [1.0,
                   ϵ,
                   ϵ^2 + 2ν,
                   ϵ^3 + 6ϵ*ν,
                   ϵ^4 + 12ϵ^2*ν + 12*ν^2,
                   ϵ^5 + 20ϵ^3*ν + 60ϵ*ν^2,
                   ϵ^6 + 30ϵ^4*ν + 180ϵ^2*ν^2 + 120ν^3
                  ]
    @test isapprox(moments, moments_ref, atol=1e-10, rtol=1e-10)
  end

  # Bethe DOS
  let t=2.0, dos = Keldysh.bethe_dos(t=t)
    moments = [dos_integrator(ω -> ω^n, dos) for n = 0:6]
    moments_ref = [1.0, 0., t^2, 0., 2t^4, 0., 5t^6]
    @test isapprox(moments, moments_ref, atol=1e-10, rtol=1e-10)
  end

  # Linear chain DOS
  let t=2.0, dos = Keldysh.chain_dos(t=t)
    moments = [dos_integrator(ω -> ω^n, dos) for n = 0:6]
    moments_ref = [1.0, 0., 2t^2, 0., 6t^4, 0., 20t^6]
    @test isapprox(moments, moments_ref, atol=1e-10, rtol=1e-10)
  end

  # Square lattice DOS
  let t=2.0, dos = Keldysh.square_dos(t=t)
    moments = [dos_integrator(ω -> ω^n, dos) for n = 0:6]
    moments_ref = [1.0, 0., 4t^2, 0., 36t^4, 0., 400t^6]
    @test isapprox(moments, moments_ref, atol=1e-10, rtol=1e-10)
  end
end


text:
using ManifoldML, Documenter

makedocs(
    format = Documenter.HTML(prettyurls = false, assets = ["assets/favicon.ico"]),
    modules = [ManifoldML],
    authors = "Ronny Bergmann.",
    sitename = "ManifoldML.jl",
    pages = [
        "Home" => "index.md",
        "K means" => "kmeans.md",
        "K Nearest Neighbors" => "knn.md",
        "Tangent space models" => "tangent-models.md",
    ],
)

deploydocs(repo = "github.com/JuliaManifolds/ManifoldML.jl.git", push_preview = true)

max_stars_repo_name: JuliaManifolds/ManifoldML

max_issues_repo_name: JuliaManifolds/ManifoldML

max_forks_repo_name: JuliaManifolds/ManifoldML


text:
using Arborist
using Test

@testset "Arborist.jl" begin
    test_files = (
        "ast_tools.jl",
        "reflection.jl",
    )
    @testset "$file" for file in test_files
        include(file)
    end
end