text:
module NonstationaryModelLearning

# greet() = print("Hello World!")

using Random
const ENV_RNG = MersenneTwister(0)

abstract type BayesFilter end

include("env_specs.jl")
include("environments.jl")
include(joinpath("multinomialtask","multinomialtask.jl"))
include(joinpath("gaussiantask","gaussiantask.jl"))

include("approximations.jl")

include(joinpath("multinomialtask","tlearners", "tparticlefilter.jl"))
include(joinpath("multinomialtask","tlearners", "tvarsmile.jl"))
include(joinpath("multinomialtask","tlearners", "tsmileOriginal.jl"))
include(joinpath("multinomialtask","tlearners", "tleakyintegrator.jl"))
 include(joinpath("multinomialtask","tlearners", "tmpn.jl"))
include(joinpath("multinomialtask","tlearners", "tsor.jl"))
include(joinpath("multinomialtask","tlearners", "tbayesfilter.jl"))

include(joinpath("gaussiantask","glearners", "gvarsmile.jl"))
include(joinpath("gaussiantask","glearners", "gsmileOriginal.jl"))
include(joinpath("gaussiantask","glearners", "gleakyintegrator.jl"))
include(joinpath("gaussiantask","glearners", "gmpn.jl"))
include(joinpath("gaussiantask","glearners", "gsor.jl"))
include(joinpath("gaussiantask","glearners", "gbayesfilter.jl"))
include(joinpath("gaussiantask","glearners", "gparticlefilter.jl"))
include(joinpath("gaussiantask","glearners", "gnassar12original.jl"))
include(joinpath("gaussiantask","glearners", "gnassar10original.jl"))
include(joinpath("gaussiantask","glearners", "gnassar12.jl"))
include(joinpath("gaussiantask","glearners", "gnassar10.jl"))

include(joinpath("gaussiantask", "calcgaussianterms.jl"))

include("tools.jl")
include("runtools.jl")
include("callbacks.jl")
include("learn.jl")
include("getparameterset.jl")
include("errormeasures.jl")
include(joinpath("multinomialtask","gridsearch_multinomial.jl"))
include(joinpath("multinomialtask","run_gridsearch_multinomial.jl"))
include(joinpath("multinomialtask","compare_multinomial.jl"))
include(joinpath("gaussiantask","gridsearch_gaussian.jl"))
include(joinpath("gaussiantask","run_gridsearch_gaussian.jl"))
include(joinpath("gaussiantask","compare_gaussian.jl"))

include("run_comparisons.jl")
include("examples.jl")

include(joinpath("processandplot","gridplot.jl"))
include(joinpath("processandplot","run_gridplot.jl"))
include(joinpath("processandplot","df_load.jl"))
include(joinpath("processandplot","df_analyze.jl"))
include(joinpath("processandplot","df_analyze_heatmap.jl"))
include(joinpath("processandplot","process_comparisons.jl"))
include(joinpath("processandplot","process_robustness.jl"))

include("grid_makeparameterDict.jl")

include(joinpath("processandplot","df_edit.jl"))

include(joinpath("processandplot","jitterplotdata.jl"))

include(joinpath("experiment","experiment.jl"))
include(joinpath("experiment","experiment_first.jl"))
include(joinpath("experiment","experiment_second.jl"))
include(joinpath("experiment","run_experiments.jl"))
include(joinpath("experiment","experiment_analysis_first.jl"))
include(joinpath("experiment","experiment_analysis_second.jl"))

include(joinpath("processandplot","reportsignificance.jl"))

end # module


text:
# This file can be run by the following command run in the folder CGP.jl/test/
# test>  julia -p 4 -L ../src/CGP.jl runtests.jl
println("test/runtests.jl")
include("test_Entropy.jl")
include("test_degen_complexity.jl")
include("test_inf_alleles.jl")
include("test_mut_evolve.jl")
include("test_pop_evolve.jl")
include("test_robust_evolve.jl")
include("test_run_mut_evolution.jl")
include("test_generate_random_functions.jl")

ext: jl


text:
module saft

using Unitful
using NLsolve

include("general.jl")

∑(x)=sum(x)
const q01, q11, q21 = [3.7039,-0.3226,0.6907]u"Å"
const q02, q12, q22 = [0.06233,-0.02236,-0.01563]u"mol/g"
const q03, q13, q23 = [150.03,80.68,38.96]u"K"
const MCH4 = 16.043u"g/mol"
const R = 8.31446261815324u"MPa*cm^3/(K*mol)"
const Nₐ = 6.02214076e23u"1/(mol)"
const kᵦ = R/Nₐ

const a0 = [0.9105631445,0.6361281449,2.6861347891,-26.547362491,97.759208784,-159.59154087,91.297774084]
const a1 = [-0.3084016918,0.1860531159,-2.5030047259,21.419793629,-65.255885330,83.318680481,-33.746922930]
const a2 = [-0.0906148351,0.4527842806,0.5962700728,-1.7241829131,-4.1302112531,13.776631870,-8.6728470368]
const b0 = [0.7240946941,2.2382791861,-4.0025849485,-21.003576815,26.855641363,206.55133841,-355.60235612]
const b1 = [-0.5755498075,0.6995095521,3.8925673390,-17.215471648,192.67226447,-161.82646165,-165.20769346]
const b2 = [0.0976883116,-0.2557574982,-9.1558561530,20.642075974,-38.804430052,93.626774077,-29.666905585]

"""
    mutable struct Mix

gfhgfh
"""
mutable struct Mix
    n::Int
    componentes::Vector{String}
    MW::typeof([1.0]u"g/mol")
    x::typeof([1.0])
    y::typeof([1.0])
    m::typeof([1.0])
    σ::typeof([1.0]u"Å")
    ϵ::typeof([1.0]u"K")
    σᵢⱼ::typeof([1.0 1.0;1.0 1.0]u"Å")
    ϵᵢⱼ::typeof([1.0 1.0;1.0 1.0]u"K")
    k::Array{Float64,2}
end

function Mix(MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),k::Array{Float64,2})
    n = length(MW)
    y = x
    componentes = ["Componente $i" for i=1:n]
    Mix(n,componentes,MW,x,y,k)
end

function Mix(componentes::Vector{String},MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),k::Array{Float64,2})
    n = length(MW)
    y = x
    Mix(n,componentes,MW,x,y,k)
end

function Mix(MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),y::typeof([1.0]),k::Array{Float64,2})
    n = length(MW)
    componentes = ["Componente $i" for i=1:n]
    Mix(n,componentes,MW,x,y,k)
end

function Mix(n::Int,MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),k::Array{Float64,2})
    y = x
    componentes = ["Componente $i" for i=1:n]
    Mix(n,componentes,MW,x,y,k)
end

function Mix(n::Int,componentes::Vector{String},MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),k::Array{Float64,2})
    y = x
    Mix(n,componentes,MW,x,y,k)
end

function Mix(componentes::Vector{String},MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),y::typeof([1.0]),k::Array{Float64,2})
    n = length(MW)
    Mix(n,componentes,MW,x,y,k)
end

function Mix(n::Int,MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),y::typeof([1.0]),k::Array{Float64,2})
    componentes = ["Componente $i" for i=1:n]
    Mix(n,componentes,MW,x,y,k)
end

function Mix(n::Int,componentes::Vector{String},MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),y::typeof([1.0]),k::Array{Float64,2})
    σ(MW) = q01+(MW-MCH4)*q11/MW+(MW-MCH4)*(MW-2*MCH4)*q21/MW^2
    m(MW) =(q02 +(MW - MCH4)*q12/MW+(MW - MCH4)*(MW - 2*MCH4)*q22/MW^2)*MW
    ϵ(MW) = q03 +(MW - MCH4)*q13/MW+(MW - MCH4)*(MW - 2*MCH4)*q23/MW^2

    m = [m(MW[i]) for i=1:n]
    σ = [σ(MW[i]) for i=1:n]
    ϵ = [ϵ(MW[i]) for i=1:n]

    σᵢⱼ = Array{Float64}(undef,n,n)*u"Å"
    ϵᵢⱼ = Array{Float64}(undef,n,n)*u"K"
    for i=1:n
        for j=1:n
            σᵢⱼ[i,j] = 0.5*(σ[i]+σ[j])
            ϵᵢⱼ[i,j] = sqrt(ϵ[i]*ϵ[j])*(1-k[i,j])
        end
    end

    Mix(n,componentes,MW,x,y,m,σ,ϵ,σᵢⱼ,ϵᵢⱼ,k)
end

function Mix(MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),m::typeof([1.0]),σ::typeof([1.0]u"Å"),ϵ::typeof([1.0]u"K"),k::Array{Float64,2})
    n = length(MW)
    y = x
    componentes = ["Componente $i" for i=1:n]
    Mix(n,componentes,MW,x,y,m,σ,ϵ,k)
end

function Mix(componentes::Vector{String},MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),m::typeof([1.0]),σ::typeof([1.0]u"Å"),ϵ::typeof([1.0]u"K"),k::Array{Float64,2})
    n = length(MW)
    y = x
    Mix(n,componentes,MW,x,y,m,σ,ϵ,k)
end

function Mix(MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),y::typeof([1.0]),m::typeof([1.0]),σ::typeof([1.0]u"Å"),ϵ::typeof([1.0]u"K"),k::Array{Float64,2})
    n = length(MW)
    componentes = ["Componente $i" for i=1:n]
    Mix(n,componentes,MW,x,y,m,σ,ϵ,k)
end

function Mix(n::Int,MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),m::typeof([1.0]),σ::typeof([1.0]u"Å"),ϵ::typeof([1.0]u"K"),k::Array{Float64,2})
    y = x
    componentes = ["Componente $i" for i=1:n]
    Mix(n,componentes,MW,x,y,m,σ,ϵ,k)
end

function Mix(n::Int,componentes::Vector{String},MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),m::typeof([1.0]),σ::typeof([1.0]u"Å"),ϵ::typeof([1.0]u"K"),k::Array{Float64,2})
    y = x
    Mix(n,componentes,MW,x,y,m,σ,ϵ,k)
end

function Mix(componentes::Vector{String},MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),y::typeof([1.0]),m::typeof([1.0]),σ::typeof([1.0]u"Å"),ϵ::typeof([1.0]u"K"),k::Array{Float64,2})
    n = length(MW)
    Mix(n,componentes,MW,x,y,m,σ,ϵ,k)
end

function Mix(n::Int,MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),y::typeof([1.0]),m::typeof([1.0]),σ::typeof([1.0]u"Å"),ϵ::typeof([1.0]u"K"),k::Array{Float64,2})
    componentes = ["Componente $i" for i=1:n]
    Mix(n,componentes,MW,x,y,m,σ,ϵ,k)
end

function Mix(n::Int,componentes::Vector{String},MW::typeof([1.0]u"g/mol"),x::typeof([1.0]),y::typeof([1.0]),m::typeof([1.0]),σ::typeof([1.0]u"Å"),ϵ::typeof([1.0]u"K"),k::Array{Float64,2})
    σᵢⱼ = Array{Float64}(undef,n,n)*u"Å"
    ϵᵢⱼ = Array{Float64}(undef,n,n)*u"K"
    for i=1:n
        for j=1:n
            σᵢⱼ[i,j] = 0.5*(σ[i]+σ[j])
            ϵᵢⱼ[i,j] = sqrt(ϵ[i]*ϵ[j])*(1-k[i,j])
        end
    end

    Mix(n,componentes,MW,x,y,m,σ,ϵ,σᵢⱼ,ϵᵢⱼ,k)
end

const mix = Mix

function zdisp(T,mix,ρ,phase)
    n = mix.n
    m,σ,ϵ,k = mix.m,mix.σ,mix.ϵ,mix.k
    σᵢⱼ,ϵᵢⱼ = mix.σᵢⱼ,mix.ϵᵢⱼ
    phase == "gas" ? x = mix.y : x = mix.x
    T = T |> u"K"

    d = σ.*(1 .-0.12*exp.(-3*ϵ/T))
    m̄ = ∑(x.*m)

    a = a0+ (m̄-1)*a1/m̄ + (m̄-1)*(m̄-2)*a2/m̄^2
    b = b0+ (m̄-1)*b1/m̄ + (m̄-1)*(m̄-2)*b2/m̄^2

    η = ∑(x.*m.*d.^3)
    η = η*π/6*ρ

    ∂ηI₁ = ∑(a[j]*j*η^(j-1) for j=1:7)
    ∂ηI₂ = ∑(b[j]*j*η^(j-1) for j=1:7)

    term1 = m̄*(8*η-2*η^2)/(1-η)^4;
    term2 = (1-m̄)*(20*η-27*η^2+12*η^3-2*η^4)/((1-η)*(2-η))^2;
    C₁ = (1+term1 + term2)^-1

    term1 = m̄*(-4*η^2+20*η+8)/(1-η)^5;
    term2 = (1-m̄)*(2*η^3+12*η^2-48*η+40)/((1-η)*(2-η))^3;
    C₂ = -C₁^2*(term1 + term2)

    m²ϵσ³ = 0u"Å^3"
    m²ϵ²σ³ = 0u"Å^3"

    for j=1:n
        for i=1:n
            m²ϵσ³ += x[i]*x[j]*m[i]*m[j]*(ϵᵢⱼ[i,j]/T)*σᵢⱼ[i,j]^3
            m²ϵ²σ³ += x[i]*x[j]*m[i]*m[j]*(ϵᵢⱼ[i,j]/T)^2*σᵢⱼ[i,j]^3
        end
    end

    I₂ = ∑(b[i]*η^(i-1) for i=1:7)

    z_disp = -2*π*ρ*∂ηI₁*m²ϵσ³ -π*ρ*m̄*(C₁*∂ηI₂+C₂*η*I₂)*m²ϵ²σ³

end

function z_hc(T,mix,ρ,phase)
    n = mix.n
    m,σ,ϵ,k = mix.m,mix.σ,mix.ϵ,mix.k
    phase == "gas" ? x = mix.y : x = mix.x
    T = T |> u"K"

    d = σ.*(1 .-0.12*exp.(-3*ϵ/T))
    m̄ = ∑(x.*m)

    ξ(nn) = (π*ρ/6)*∑(x[i]*m[i]*d[i]^nn for i=1:n)
    #ρ=ρ
    #η=η
    ξ = [ξ(i) for i=0:3]
    dᵢⱼ(i,j) = d[i]*d[j]/(d[i]+d[j])

    gʰˢ = Array{Float64}(undef,n,n)
    for j=1:n
        for i=1:n
            term1 = 1/(1-ξ[4]);
            term2 = dᵢⱼ(i,j)*3*ξ[3]/(1-ξ[4])^2;
            term3 = (dᵢⱼ(i,j))^2*2*ξ[3]^2/(1-ξ[4])^3;
            gʰˢ[i,j] = term1 + term2 + term3
        end
    end

    term1 = ξ[4]/(1-ξ[4]);
    term2 = 3*ξ[2]*ξ[3]/(ξ[1]*(1-ξ[4])^2);
    term3 = (3*ξ[3]^3-ξ[4]*ξ[3]^3)/(ξ[1]*(1-ξ[4])^3);
    Zʰˢ = term1 + term2 + term3

    ∂g = Array{Float64}(undef,n,n)
    for j=1:n
        for i=1:n
            term1 = ξ[4]/(1-ξ[4])^2;
            term2 = dᵢⱼ(i,j)*(3*ξ[3]/(1-ξ[4])^2+6*ξ[3]*ξ[4]/(1-ξ[4])^3);
            term3 = (dᵢⱼ(i,j))^2*(4*ξ[3]^2/(1-ξ[4])^3 + 6*ξ[3]^2*ξ[4]/(1-ξ[4])^4);
            ∂g[i,j] = term1 + term2 + term3
        end
    end

    z_hc = m̄*Zʰˢ - ∑(x[i]*(m[i]-1)*∂g[i,i]*gʰˢ[i,i]^-1 for i=1:n)

end

function obj_saft(η,T,P,mix,phase)
    n = mix.n
    m,σ,ϵ,k = mix.m,mix.σ,mix.ϵ,mix.k
    phase == "gas" ? x = mix.y : x = mix.x
    T = T |> u"K"
    P = P |> u"Pa"

    d = σ.*(1 .-0.12*exp.(-3*ϵ/T))
    m̄ = ∑(x.*m)

    η_new = ∑(x.*m.*d.^3)
    ρ = 6/π*η*η_new^-1
    Zʰᶜ = z_hc(T,mix,ρ,phase)
    Zᵈⁱˢᵖ = zdisp(T,mix,ρ,phase)

    Zcalc = 1 + Zʰᶜ + Zᵈⁱˢᵖ
    Pcalc = Zcalc*T*ρ*1.380649e-23u"J/K"
    Pcalc = Pcalc |> u"Pa"

    ustrip(abs(P-Pcalc))
end

function compresivilidad(η,T,mix,phase)
    n = mix.n
    m,σ,ϵ,k = mix.m,mix.σ,mix.ϵ,mix.k
    phase == "gas" ? x = mix.y : x = mix.x
    T = T |> u"K"

    d = σ.*(1 .-0.12*exp.(-3*ϵ/T))
    m̄ = ∑(x.*m)

    η_new = ∑(x.*m.*d.^3)
    ρ = 6/π*η*η_new^-1
    Zʰᶜ = z_hc(T,mix,ρ,phase)
    Zᵈⁱˢᵖ = zdisp(T,mix,ρ,phase)

    Zcalc = 1 + Zʰᶜ + Zᵈⁱˢᵖ
end

function compr(T,P,mix,phase)
    n = mix.n
    m,σ,ϵ,k = mix.m,mix.σ,mix.ϵ,mix.k
    x,M,y = mix.x,mix.MW,mix.y
    T = T |> u"K"

    if phase == "gas"
        ρ = P/(0.8*(R/Nₐ)*T) |> u"Å^-3"
        ηini = ∑(y.*m.*σ.^3)
        ηini = ηini*π/6*ρ
    elseif phase == "liq"
        MW = ∑(M)/n
        ρ = (800000u"g/m^3")*Nₐ/MW |> u"Å^-3"
        ηini = ∑(x.*m.*σ.^3)
        ηini = ηini*π/6*ρ
    else
        error("No especifico el tipo de fase")
    end
    obj_saft1(η) = obj_saft(abs(η),T,P,mix,phase)
    ηsol = nlsolve(n_ary(obj_saft1),[ηini],iterations=100)
    η = ηsol.zero
    #println("P, $P, x, $x, η, $η, $phase")
    z = compresivilidad(η...,T,mix,phase)
    return z,η
end

function fugF(T,P,mix,phase)
    T = T |> u"K"

    Z,η = compr(T,P,mix,phase)
    ρ = P/(Z*kᵦ*T) |> u"Å^-3"
    μʰᶜ = mu_hc(T,mix,ρ,phase)
    μᵈⁱˢᵖ = mu_disp(T,mix,ρ,phase)
    μ = exp.(μʰᶜ + μᵈⁱˢᵖ .-log(Z))
    return μ
end

function mu_hc(T,mix,ρ,phase)
    n = mix.n
    m,σ,ϵ,k = mix.m,mix.σ,mix.ϵ,mix.k
    phase == "gas" ? x = mix.y : x = mix.x
    T = T |> u"K"

    d = σ.*(1 .-0.12*exp.(-3*ϵ/T))
    m̄ = ∑(x.*m)

    ξ(nn) = (π*ρ/6)*∑(x[i]*m[i]*d[i]^nn for i=1:n)
    #ρ=ρ
    #η=η
    ξ = [ξ(i) for i=0:3]
    dᵢⱼ(i,j) = d[i]*d[j]/(d[i]+d[j])
    dᵢⱼ = d.*d./(d+d)

    gʰˢ = Vector{Float64}(undef,n)
    for i=1:n
        term1 = 1/(1-ξ[4]);
        term2 = dᵢⱼ[i]*3*ξ[3]/(1-ξ[4])^2;
        term3 = (dᵢⱼ[i])^2*2*ξ[3]^2/(1-ξ[4])^3;
        gʰˢ[i] = term1 + term2 + term3
    end

    Zʰᶜ = z_hc(T,mix,ρ,phase)

    term1 = 3*ξ[2]*ξ[3]/(1-ξ[4]);
    term2 = ξ[3]^3/(ξ[4]*(1-ξ[4])^2);
    term3 = (ξ[3]^3/ξ[4]^2-ξ[1])*log(1-ξ[4]);
    aʰˢ = (1/ξ[1])*(term1 + term2 + term3)

    Aʰᶜ = m̄*aʰˢ-∑(x[i]*(m[i]-1)*log(gʰˢ[i]) for i=1:n)

    ξₖ = [(π*ρ/6)*m[i]*d[i]^(j-1) for j=1:4,i=1:n]

    ∂aʰˢ∂x = Vector{Float64}(undef,n)
    for i=1:n
        term1 = -ξₖ[1,i]/ξ[1]*aʰˢ;
        term2 = 3*(ξₖ[2,i]*ξ[3] + ξ[2]*ξₖ[3,i])/(1-ξ[4]);
        term3 = 3*ξ[2]*ξ[3]*ξₖ[4,i]/(1-ξ[4])^2;
        term4 = 3*ξ[3]^2*ξₖ[3,i]/(ξ[4]*(1-ξ[4])^2);
        term5 = ξ[3]^3*ξₖ[4,i]*(3*ξ[4]-1)/(ξ[4]^2*(1-ξ[4])^3);
        term6 = ((3*ξ[3]^2*ξₖ[3,i]*ξ[4]-2*ξ[3]^3*ξₖ[4,i])/ξ[4]^3-ξₖ[1,i])*log(1-ξ[4]);
        term7 = (ξ[1]-ξ[3]^3/ξ[4]^2)*ξₖ[4,i]/(1-ξ[4]);

        ∂aʰˢ∂x[i] = term1 + 1/ξ[1]*(term2 + term3 + term4 + term5 + term6 +term7)
    end

    ∂g∂x = Array{Float64}(undef,n,n)
    for k=1:n
        for i=1:n
            term1 = ξₖ[4,k]/(1-ξ[4])^2;
            term2 = dᵢⱼ[i]*(3*ξₖ[3,k]/(1-ξ[4])^2 + 6*ξ[3]*ξₖ[4,k]/(1-ξ[4])^3);
            term3 = (dᵢⱼ[i])^2*(4*ξ[3]*ξₖ[3,k]/(1-ξ[4])^3 + 6*ξ[3]^2*ξₖ[4,k]/(1-ξ[4])^4);
            ∂g∂x[i,k] = term1 + term2 + term3
        end
    end

    ∂aʰᶜ∂x = m*aʰˢ + m̄*∂aʰˢ∂x -(m.-1).*log.(gʰˢ) -∑(x[i]*(m[i]-1)*∂g∂x[i,:]*gʰˢ[i]^(-1) for i=1:n)

    μ = Aʰᶜ + Zʰᶜ .+ ∂aʰᶜ∂x .- ∑(x[i]*∂aʰᶜ∂x[i] for i=1:n)
end

function HelmholtzDisp(T,mix,ρ,phase)
    n = mix.n
    m,σ,ϵ,k = mix.m,mix.σ,mix.ϵ,mix.k
    σᵢⱼ,ϵᵢⱼ = mix.σᵢⱼ,mix.ϵᵢⱼ
    phase == "gas" ? x = mix.y : x = mix.x
    T = T |> u"K"

    d = σ.*(1 .-0.12*exp.(-3*ϵ/T))
    m̄ = ∑(x.*m)

    a = a0+ (m̄-1)*a1/m̄ + (m̄-1)*(m̄-2)*a2/m̄^2
    b = b0+ (m̄-1)*b1/m̄ + (m̄-1)*(m̄-2)*b2/m̄^2

    η = ∑(x.*m.*d.^3)
    η = η*π/6*ρ

    term1 = m̄*(8*η-2*η^2)/(1-η)^4;
    term2 = (1-m̄)*(20*η-27*η^2+12*η^3-2*η^4)/((1-η)*(2-η))^2;
    C₁ = (1+term1 + term2)^-1

    m²ϵσ³ = 0u"Å^3"
    m²ϵ²σ³ = 0u"Å^3"

    for j=1:n
        for i=1:n
            m²ϵσ³ += x[i]*x[j]*m[i]*m[j]*(ϵᵢⱼ[i,j]/T)*σᵢⱼ[i,j]^3
            m²ϵ²σ³ += x[i]*x[j]*m[i]*m[j]*(ϵᵢⱼ[i,j]/T)^2*σᵢⱼ[i,j]^3
        end
    end

    I₁ = ∑(a[i]*η^(i-1) for i=1:7)
    I₂ = ∑(b[i]*η^(i-1) for i=1:7)

    A_disp = -2*π*ρ*I₁*m²ϵσ³-π*ρ*m̄*C₁*I₂*m²ϵ²σ³

end

function mu_disp(T,mix,ρ,phase)
    n = mix.n
    m,σ,ϵ,k = mix.m,mix.σ,mix.ϵ,mix.k
    σᵢⱼ,ϵᵢⱼ = mix.σᵢⱼ,mix.ϵᵢⱼ
    phase == "gas" ? x = mix.y : x = mix.x
    T = T |> u"K"

    d = σ.*(1 .-0.12*exp.(-3*ϵ/T))
    m̄ = ∑(x.*m)

    a = a0+ (m̄-1)*a1/m̄ + (m̄-1)*(m̄-2)*a2/m̄^2
    b = b0+ (m̄-1)*b1/m̄ + (m̄-1)*(m̄-2)*b2/m̄^2

    η = ∑(x.*m.*d.^3)
    η = η*π/6*ρ

    Zᵈⁱˢᵖ = zdisp(T,mix,ρ,phase)
    Aᵈⁱˢᵖ = HelmholtzDisp(T,mix,ρ,phase)

    ξₖ = [(π*ρ/6)*m[i]*d[i]^(j-1) for j=1:4,i=1:n]

    m²ϵσ³ = 0u"Å^3"
    m²ϵ²σ³ = 0u"Å^3"

    for j=1:n
        for i=1:n
            m²ϵσ³ += x[i]*x[j]*m[i]*m[j]*(ϵᵢⱼ[i,j]/T)*σᵢⱼ[i,j]^3
            m²ϵ²σ³ += x[i]*x[j]*m[i]*m[j]*(ϵᵢⱼ[i,j]/T)^2*σᵢⱼ[i,j]^3
        end
    end

    m²ϵσ³ₓₖ = 2*m.*∑(x[j]*m[j]*(ϵᵢⱼ[:,j]/T).*σᵢⱼ[:,j].^3 for j=1:n)
    m²ϵ²σ³ₓₖ = 2*m.*∑(x[j]*m[j]*((ϵᵢⱼ[:,j]/T).^2).*σᵢⱼ[:,j].^3 for j=1:n)

    I₁ = ∑(a[i]*η^(i-1) for i=1:7)
    I₂ = ∑(b[i]*η^(i-1) for i=1:7)

    term1 = m̄*(8*η-2*η^2)/(1-η)^4;
    term2 = (1-m̄)*(20*η-27*η^2+12*η^3-2*η^4)/((1-η)*(2-η))^2;
    C₁ = (1+term1 + term2)^-1

    term1 = m̄*(-4*η^2+20*η+8)/(1-η)^5;
    term2 = (1-m̄)*(2*η^3+12*η^2-48*η+40)/((1-η)*(2-η))^3;
    C₂ = -C₁^2*(term1 + term2)

    C₁ₓₖ = C₂*ξₖ[4,:] - (C₁^2)*(m*(8*η-2*η^2)/(1-η)^4 -m*(20*η-27*η^2+12*η^3-2*η^4)/((1-η)*(2-η))^2)

    aₖ = a1*m'/m̄^2 +a2*m'*(3-4/m̄)/m̄^2
    bₖ = b1*m'/m̄^2 +b2*m'*(3-4/m̄)/m̄^2

    I₁ₓₖ = Vector{Float64}(undef,n)
    I₂ₓₖ = Vector{Float64}(undef,n)
    for j=1:n
        sum1 = ∑((a[i]*(i-1)*ξₖ[4,j]*η^(i-2))+aₖ[i,j]*η^(i-1) for i=1:7)
        sum2 = ∑((b[i]*(i-1)*ξₖ[4,j]*η^(i-2))+bₖ[i,j]*η^(i-1) for i=1:7)
        I₁ₓₖ[j] = sum1
        I₂ₓₖ[j] = sum2
    end

    ∂aᵈⁱˢᵖ∂x = -2*π*ρ*(I₁ₓₖ*m²ϵσ³+I₁*m²ϵσ³ₓₖ) -π*ρ*((m*C₁*I₂+m̄*C₁ₓₖ*I₂+m̄*C₁*I₂ₓₖ)*m²ϵ²σ³ +m̄*C₁*I₂*m²ϵ²σ³ₓₖ)

    μ = Aᵈⁱˢᵖ + Zᵈⁱˢᵖ .+ ∂aᵈⁱˢᵖ∂x .- ∑(x[i]*∂aᵈⁱˢᵖ∂x[i] for i=1:n)

end

function obj_buble(T,P,mix)
    iter = 0
    error = 100
    K = Vector{Float64}(undef, mix.n)
    while error > 1e-5 && iter < 100
        iter +=1
        fliq = fugF(T,P,mix,"liq")
        fgas = fugF(T,P,mix,"gas")
        K = fliq./fgas
        y_calc = mix.x.*K
        y_sum = ∑(y_calc)
        #println("mix.y= $(mix.y), y_calc= $y_calc")
        error = ∑(abs.(y_calc-mix.y))
        mix.y = y_calc/y_sum
        P = y_sum*P
    end

    return P,mix.y

end

"""
    pxy(T,P,mix::mix,xx;uni=u"Torr")

**Los campos de entrada son:**
- `T :: Float` Es la temperatura del sistema
- `P :: m-element Array{Float64,1}` presión inicial
- `mix` Datos de entrada de la mezcla
- `xx :: StepRange ó Vector` Son los puntos en el liquido donde se buscara el equilibrio con el vapor
- `uni` Son las unidades en las que se desea el resultado, por defecto esta en Torr

**Salida:**

Regresa las fracciones de liquido `x`,
las fracciones del vapor `y` y la presión `P` en Torr (`m` es la longitud de `xx`).

- `x :: n×m Array{Float64,2}` liquido
- `y :: n×m Array{Float64,2}` vapor
- `P :: m-element Array{Float64,1}` presión
"""
function pxy(T,P,mix::mix,xx;uni=u"Torr")
    T = T |> u"K"
    P = P |> u"Torr"
    xx = general.nc(xx...)
    m = size(xx,2)
    y = Array{Float64}(undef,mix.n,m)
    p = Vector{Float64}(undef,m)*u"Torr"

    mix.x = xx[:,1]
    mix.y = mix.x
    p[1],y[:,1] = obj_buble(T,P,mix)
    for i=2:m
        mix.x = xx[:,i]
        mix.y = y[:,i-1]
        p[i],y[:,i] = obj_buble(T,p[i-1],mix)
    end

    p = p .|> uni
    return xx,y,p

end

end

max_stars_repo_name: EmilioAlvizo/TermoQuimica

max_issues_repo_name: EmilioAlvizo/TermoQuimica

max_forks_repo_name: EmilioAlvizo/TermoQuimica


text:
using JLD
simdata = load("CTSV.jld", "simdata")
writedlm("simdata", simdata)


text:
#************************ LOCALS TO CONFIGURE!!!! **************************
email= ENV["NCBI_EMAIL"] #Enviroment variable that need to be setup
umls_user = ENV["UMLS_USER"]
umls_pswd = ENV["UMLS_PSSWD"]
search_term="(obesity[MeSH Major Topic]) AND (\"2010\"[Date - Publication] : \"2012\"[Date - Publication])"
max_articles = 10
overwrite_db=true
verbose = false

#************************ MYSQL **********************************************
host="localhost" #If want to hide - use enviroment variables instead
mysql_usr="root"
mysql_pswd=""
dbname="pubmed_processes_test"
#*****************************************************************************

db = nothing

@testset "Search and Save" begin
    println("-----------------------------------------")
    println("       Testing Search and Save")
    db_config = Dict(:host=>host,
                     :dbname=>dbname,
                     :username=>mysql_usr,
                     :pswd=>mysql_pswd,
                     :overwrite=>overwrite_db)

    println()
    @time db = pubmed_search_and_save(email, search_term, max_articles,
    save_efetch_mysql, db_config, verbose)

    #query the article table and make sure the count is correct
    all_pmids = BioMedQuery.Entrez.DB.all_pmids(db)
    @test length(all_pmids) == max_articles

end

@testset "MESH2UMLS" begin
    println("-----------------------------------------")
    println("       Testing MESH2UMLS")
    credentials = Credentials(umls_user, umls_pswd)
    append = false

    @time begin
        map_mesh_to_umls_async!(db, credentials; append_results=append)
    end

    all_pairs_query = db_query(db, "SELECT mesh FROM mesh2umls;")
    all_pairs = all_pairs_query[1]
    @test length(all_pairs) > 0
    @test isa(all_pairs, DataArrays.DataArray{AbstractString,1})

end

@testset "Occurrences" begin
    println("-----------------------------------------")
    println("       Testing Occurrences")
    umls_concept = "Disease or Syndrome"
    @time begin
        labels2ind, occur = umls_semantic_occurrences(db, umls_concept)
    end

    @test length(keys(labels2ind)) > 0
    @test length(find(x->x=="obesity", collect(keys(labels2ind)))) ==1
end

db_query(db, "DROP DATABASE IF EXISTS $dbname;")

println("------------End Test Processes MySQL-----------")
println("-----------------------------------------------")


text:
# This file defines the hidden dynamics system problem
# with (1) background noise and sparse large
# deviations, (2) "broken sensor" noise, and
# (3) an added "bump". L2, robust (huber), and
# robust (trimming) fits are
# run for several background noise levels
# and the quality of the resulting eigenvalues is
# returned

using Optim

################################################################
################################################################
################################################################
# Define some helper routines

function hidden_dynamics_system_data(;nt=2^9,
                                     nx=300,k1=1.0,
                                     omega1=1.0,
                                     gamma1=1.0,
                                     k2=0.4,
                                     omega2=3.7,
                                     gamma2=-0.2)
################################################
# Generate data for hidden dynamics system
#
# z(t,x) = exp(gamma1*t)*sin(k1*x-omega1*t) +
#      exp(gamma2*t)*sin(k2*x-omega2*t) +
#
# as gamma1 is positive and gamma2 is negative,
# the mode for gamma2 becomes hidden, especially
# with significant noise.
#
################################################

    evals = [gamma1 + im*omega1; gamma1 - im*omega1;
             gamma2 + im*omega2; gamma2 - im*omega2]

    x = collect(range(0,stop=15.0,length=nx))
    t = complex(collect(range(0,stop=2*pi,length=nt)))

    z = complex(exp.(gamma1*t) .* sin.(k1*transpose(x) .- omega1*t) +
         exp.(gamma2*t) .* sin.(k2*transpose(x) .- omega2*t))

    return z, evals, t, x
end

function make_corrupt_data_sparse_spikes(xclean,sigma,mu,p)
    #------------------------------------------
    # Generate Corrupted Data - Spikes
    #------------------------------------------
    #
    # sigma - background noise
    # mu - spike size
    # p - spike frequency

    m, n = size(xclean)

    noise = sigma*randn(m,n) .+ mu*randn(m,n).*( rand(m,n) .< p )
    xdat = xclean .+ noise

    return xdat
end

function make_corrupt_data_bump(xclean,x,t,sigma,A,w)
    #------------------------------------------
    # Generate Corrupted Data - bump
    #------------------------------------------
    #
    # sigma - background noise
    # A - bump amplitude
    # w - bump width

    m, n = size(xclean)

    i1 = rand(1:m)
    j1 = rand(1:n)
    x1 = x[j1]
    t1 = t[i1]

    dx = x[2]-x[1]
    dt = t[2]-t[1]

    noise = (sigma*randn(m,n) .+
             A*exp.( real(-(t .- t1).^2/(dt*w)^2 .- (transpose(x) .- x1).^2/(dx*w)^2 )))
    xdat = xclean .+ noise

    return xdat
end

function make_corrupt_data_broken_sensors(xclean,sigma,mu,p)
    #---------------------------------------------------------
    # Generate Corrupted Data - "broken sensors"
    #---------------------------------------------------------
    #
    # sigma - background noise
    # mu - spike size
    # p - spike frequency

    m, n = size(xclean)

    ncol = Integer(round(n*p))
    icols = sample(1:n,ncol,replace=false)

    noise = sigma*randn(m,n)
    noise[:,icols] .+= mu*randn(m,length(icols))
    xdat = xclean .+ noise

    return xdat
end

function dofits(xdat,t,eigs,kappa,optionsl2,optionshuber,optionstrim,
                nkeep)
    # this wraps up the calls to the DMD solver
    # Here we use the robust but generally slower
    # proximal gradient routine for both L2 and Huber
    # fits. On return, we provide the error for these
    # fits and the exact DMD error.

    k = length(eigs)

    #-------------------------------------------------------------
    # Initial Guess
    #-------------------------------------------------------------
    m,n = size(xdat);

    a0, B0 = dmdexactestimate(m, n, k, xdat, t);

    ainit, Binit = dmdexactestimate(m,n,k,xdat,t,dmdtype="trap")

    #----------------------------------------------------
    # l2 Trial
    #----------------------------------------------------
    #println("running l2 experiment...");
    # loss functions
    lossf = (z) -> l2_func(z);
    lossg = (z) -> l2_grad!(z);

    l2params = DMDParams(k, xdat, t, lossf, lossg,
                         inner_directl2=true);
    copyto!(l2params.a,ainit)

    # apply solver
    obj_his1, err_his1 = solveDMD_withPG(l2params, optionsl2);
    res = solveDMD_withOptimBFGS(l2params)

    #----------------------------------------------------
    # huber Trial
    #----------------------------------------------------
    #println("running huber experiment...");
    # loss functions
    lossf = (z) -> huber_func(z,kappa);
    lossg = (z) -> huber_grad!(z,kappa);

    huberparams = DMDParams(k, xdat, t, lossf, lossg);
    copyto!(huberparams.a,ainit)

    # apply solver
    obj_his2, err_his2 = solveDMD_withPG(huberparams, optionshuber);
    res = solveDMD_withOptimBFGS(huberparams)
    #----------------------------------------------------
    # l2 trimming trial
    #----------------------------------------------------
    #println("running l2 trimming experiment...");
    # loss functions
    lossf = (z) -> l2_func(z);
    lossg = (z) -> l2_grad!(z);

    trimparams = DMDParams(k, xdat, t, lossf, lossg,
                           inner_directl2=true,nkeep=nkeep);

    copyto!(trimparams.a,ainit)

    # apply solver
    obj_his3, err_his3 = solveDMD_withPG(trimparams, optionstrim);
    res = solveDMD_withOptimBFGS(trimparams)
    #--------------------------------------------------------------------
    # Evaluate Errors and Get Best Eigenvalue Permutation
    #--------------------------------------------------------------------

    if any(isnan,l2params.a)
        err1 = Inf
        p1 = 1:k
    else
        err1, p1 = besterrperm_wi(l2params.a, eigs);
    end
    if any(isnan,huberparams.a)
        err2 = Inf
        p2 = 1:k
    else
        err2, p2 = besterrperm_wi(huberparams.a, eigs);
    end
    if any(isnan,trimparams.a)
        err3 = Inf
        p3 = 1:k
    else
        err3, p3 = besterrperm_wi(trimparams.a, eigs);
    end
    err4, p4 = besterrperm_wi(a0, eigs);

    return (err1, err2, err3, err4, l2params.a[p1], huberparams.a[p2],
            trimparams.a[p3], a0[p4])

end


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

export BlobRef, BlobFuture, BatchFuture, fetch, fetch!, migrate, migrate!

###################################################################################################
# Blob reference: a string that describes a blob in blob storage

mutable struct BlobRef
    name::Union{String, Tuple}
end

###################################################################################################
# Batch Future: reference to batch resource on worker's disk (worker-side fetch)

mutable struct BatchFuture
    container::Union{String, Nothing}
    blob
end

BatchFuture(container) = BatchFuture(container, nothing)
BatchFuture() = BatchFuture(nothing, nothing)

# Fetch on a batch worker (if executed on user-side -> fetch data from blob)
function fetch(arg::BatchFuture)
    try
        iostream = open(arg.blob.name, "r")
        data = deserialize(iostream)
        close(iostream)
        return data
    catch
        iostream = __clients__[1]["blob_client"].get_blob_to_bytes(arg.container, arg.blob.name)
        return deserialize(IOBuffer(iostream.content))
    end
end

function fetch!(arg::BatchFuture)
    data = nothing
    try
        iostream = open(arg.blob.name, "r")
        data = deserialize(iostream)
        close(iostream)
    catch
        iostream = __clients__[1]["blob_client"].get_blob_to_bytes(arg.container, arg.blob.name)
        data = deserialize(IOBuffer(iostream.content))
    end
    arg.blob = data
    return arg.blob
end

###################################################################################################
# Blob Future: reference to blob in blob storage (user-side fetch)

mutable struct BlobFuture
    container::Union{String, Nothing}
    blob
    client_index::Integer
end

BlobFuture(container, blob) = BlobFuture(container, blob, 1)
BlobFuture(container) = BlobFuture(container, nothing, 1)
BlobFuture() = BlobFuture(nothing, nothing, 1)

function fetch(arg::BlobFuture; destroy_blob=false)

    # Try to read future from disk, otherwise fetch from blob
    i = arg.client_index
    try
        num_files = length(arg.blob.name)
        out_files = []
        for blob in arg.blob.name

            # Fetch blob and add to collection
            if ~isnothing(__clients__[i]["blob_client"])
                iostream = __clients__[i]["blob_client"].get_blob_to_bytes(arg.container, blob)
                push!(out_files, deserialize(IOBuffer(iostream.content)))

                # Delete blob (default is true)
                destroy_blob && __clients__[i]["blob_client"].delete_blob(arg.container, blob)
            end
        end

        if num_files > 1
            return tuple(out_files...)
        elseif length(out_files) == 1
            return out_files[1]
        else
            return nothing
        end
    catch
        throw("Blob does not (yet) exist or BlobFuture does not contain a proper blob reference.")
    end
end

function fetch!(arg::BlobFuture; destroy_blob=false)

    # Try to read future from disk, otherwise fetch from blob
    i = arg.client_index
    try
        num_files = length(arg.blob.name)
        out_files = []
        for blob in arg.blob.name

            # Fetch blob and add to collection
            if ~isnothing(__clients__[i]["blob_client"])
                iostream = __clients__[i]["blob_client"].get_blob_to_bytes(arg.container, blob)
                push!(out_files, deserialize(IOBuffer(iostream.content)))

                # Delete blob (default is true)
                destroy_blob && __clients__[i]["blob_client"].delete_blob(arg.container, blob)
            end
        end

        if num_files > 1
            arg.blob = tuple(out_files...)
        elseif length(out_files) == 1
            arg.blob = out_files[1]
        else
            arg.blob = nothing
        end
        return arg.blob
    catch
        throw("Blob does not (yet) exist or BlobFuture does not contain a proper blob reference.")
    end
end

"""
    `migrate(blob_future::BlobFuture, destination)`

 Migrate a blob future from the storage account of one pool to another. The destination is the index of the pool to whose storage account the future is moved.

 *Arguments*:

 - `blob_future` (BlobFuture): Blob future that is to be moved to a different storage account.

 - `destination` (Integer): Pool number to whose storage account the blob future will be moved.

 *Output*:

 - `new_blob_future::BlobFuture`: New blob future at new location.

 See also:  [`fetch`](@ref)
 """
function migrate(arg::BlobFuture, destination; delete_blob=true)

    # Fetch to master
    blobs = fetch(arg)

    # Upload to new destination
    for (name, blob) in zip(arg.blob.name, blobs)

        # Serialize variable
        iobuff = IOBuffer()
        AzureClusterlessHPC.serialize(iobuff, blob)   # eval in local scope
        binary = iobuff.data

        # Create new resource
        upload_bytes_to_container(__clients__[destination]["blob_client"], arg.container, name, binary; verbose=__verbose__)

        # Remove old blob
        if delete_blob
            __clients__[arg.client_index]["blob_client"].delete_blob(arg.container, name)
        end
    end

    return BlobFuture(arg.container, arg.blob, destination)
end

"""
    `migrate!(blob_future::BlobFuture, destination)`

 Migrate a blob future from the storage account of one pool to another. The destination is the index of the pool to whose storage account the future is moved.

 *Arguments*:

 - `blob_future` (BlobFuture): Blob future that is to be moved to a different storage account.

 - `destination` (Integer): Pool number to whose storage account the blob future will be moved.

 See also:  [`fetch`](@ref)
 """
function migrate!(arg::BlobFuture, destination; delete_blob=true)

    # Fetch to master
    blobs = fetch(arg)

    # Upload to new destination
    for (name, blob) in zip(arg.blob.name, blobs)

        # Serialize variable
        iobuff = IOBuffer()
        AzureClusterlessHPC.serialize(iobuff, blob)   # eval in local scope
        binary = iobuff.data

        # Create new resource
        upload_bytes_to_container(__clients__[destination]["blob_client"], arg.container, name, binary; verbose=__verbose__)

        # Remove old blob
        if delete_blob
            __clients__[arg.client_index]["blob_client"].delete_blob(arg.container, name)
        end
    end
    arg.client_index = destination
    return arg
end


text:
module covid19abm
using Parameters, Distributions, StatsBase, StaticArrays, Random

@enum HEALTH UNDEF=0 SUSC=1 EXP=2 INFMILD=4 INFHOSP=5 INFICU=6
Base.@kwdef mutable struct Human
    idx::Int64 = 0
    health::HEALTH = SUSC
    swap::HEALTH = UNDEF
    age::Int64   = 0    # in years.
    ag::Int64    = 0
    tis::Int64   = 0   # time in state
    exp::Int64   = 0    # max statetime
    iso_time::Int64 = 0
    is_iso::Bool = false
end

@with_kw mutable struct ModelParameters @deftype Float64    ## use @with_kw from Parameters
    st::Int64 = 500 ## sim time in days
    β = 0.0
    I₀::Int64 = 1    ## initial infected
    σ::Distribution = LogNormal(log(5.2), 0.1)  ## distribution for incubation period
    γ::Int64 = 5 ## 5 days infectious period
    q::Int64 = 0.05 ## percentage of immediate self-isolation
end
Base.show(io::IO, ::MIME"text/plain", z::Human) = dump(z)

## constants
const agedist =  Categorical(@SVector [0.190674655, 0.41015843, 0.224428346, 0.17473857])
const agebraks = @SVector [0:19, 20:49, 50:64, 65:99]
const humans = Array{Human}(undef, 1000)
const p = ModelParameters()
export ModelParameters, HEALTH, Human, humans

## initialization functions
function initialize()
    @inbounds for i = 1:length(humans)
        x = humans[i]
        x = Human()              ## create an empty human
        x.idx = i
        x.ag = rand(agedist)
        x.age = rand(agebraks[x.ag])
    end
end
export initialize

function insert_exposed(num)
    ## inserts a number of infected people in the population randomly
    l = findall(x -> x.health == SUSC, humans)
    if length(l) > 0
        h = sample(l, num; replace = false)
        @inbounds for i in h
            move_to_latent(humans[i])
        end
    end
end
export insert_exposed

function move_to_latent(h::Human)
    ## make the i'th human exposed
    h.health = EXP
    h.swap = UNDEF
    h.exp = Int(ceil(rand(p.σ)))
    h.tis = 0
end

function move_to_infected(h::Human)
    ## to do: determine if infection is mild, hosp, or ICU
    h.health = INFMILD
    h.swap = UNDEF
    h.exp = p.γ
    h.tis = 0

    ## determine qurantine levels.
    # determine, f, tau, hosp, icu.
    error("not implemented")

end

## references:
# critical care capacity in Canada https://www.ncbi.nlm.nih.gov/pubmed/25888116


text:
using ProgressMeter

function isfile_drive(parentId, fileName)
    length(readall(`drive list -n -q "'$parentId' in parents and title = '$fileName'"`)) != 0
end

function dir2drive(sourceDir, destName)
    destId = split(readall(`drive list -n --title $destName`))[1]
    @showprogress for fileName in readdir(sourceDir)
        if !isfile_drive(destId, fileName)
            readall(`drive upload -f $sourceDir/$fileName --title $fileName --parent $destId`)
        end
    end
end


text:
using HypothesisTests, Test
using HypothesisTests: default_tail

@testset "Wilcoxon" begin
@testset "Basic exact test" begin
	@test abs(pvalue(ExactSignedRankTest([1:10;], [2:2:20;])) - 0.0020) <= 1e-4
	@test abs(pvalue(ExactSignedRankTest([2:2:20;], [1:10;])) - 0.0020) <= 1e-4
	@test abs(pvalue(ExactSignedRankTest([1:10;], [2:2:16; -1; 1])) - 0.4316) <= 1e-4
	@test abs(pvalue(ExactSignedRankTest([2:2:16; -1; 1], [1:10;])) - 0.4316) <= 1e-4
	@test default_tail(ExactSignedRankTest([1:10;], [2:2:20;])) == :both
	show(IOBuffer(), ExactSignedRankTest([1:10;], [2:2:20;]))
end

@testset "Exact with ties" begin
	@test abs(pvalue(ExactSignedRankTest([1:10;], [1:10;])) - 1) <= 1e-4
	@test abs(pvalue(ExactSignedRankTest([1:10;], [2:11;])) - 0.0020) <= 1e-4
	@test abs(pvalue(ExactSignedRankTest([2:11;], [1:10;])) - 0.0020) <= 1e-4
	@test abs(pvalue(ExactSignedRankTest(1:10, [1:5; ones(5)])) - 0.0625) <= 1e-4
	@test abs(pvalue(ExactSignedRankTest([1:5; ones(5)], [1:10;])) - 0.0625) <= 1e-4
	show(IOBuffer(), ExactSignedRankTest([1:10;], [1:10;]))
end

@testset "Approximate test" begin
	@test abs(pvalue(ApproximateSignedRankTest([1:10;], [2:2:20;])) - 0.005922) <= 1e-6
	@test abs(pvalue(ApproximateSignedRankTest([2:2:20;], [1:10;])) - 0.005922) <= 1e-6
	@test abs(pvalue(ApproximateSignedRankTest([1:10;], [2:2:16; -1; 1])) - 0.4148) <= 1e-4
	@test abs(pvalue(ApproximateSignedRankTest([2:2:16; -1; 1], [1:10;])) - 0.4148) <= 1e-4
	@test default_tail(ApproximateSignedRankTest([1:10;], [2:2:20;])) == :both
	show(IOBuffer(), ApproximateSignedRankTest([1:10;], [2:2:20;]))
end

@testset "Approximate with ties" begin
	@test abs(pvalue(ApproximateSignedRankTest([1:10;], [1:10;])) - 1) <= 1e-4
	@test abs(pvalue(ApproximateSignedRankTest([1:10;], [2:11;])) - 0.001904) <= 1e-6
	@test abs(pvalue(ApproximateSignedRankTest([2:11;], [1:10;])) - 0.001904) <= 1e-6
	@test abs(pvalue(ApproximateSignedRankTest([1:10;], [1:5; ones(5)])) - 0.05906) <= 1e-5
	@test abs(pvalue(ApproximateSignedRankTest([1:5; ones(5)], 1:10)) - 0.05906) <= 1e-5
	show(IOBuffer(), ApproximateSignedRankTest([1:10;], [1:10;]))
end

@testset "Tests for automatic selection" begin
	@test abs(pvalue(SignedRankTest([1:10;], [2:2:20;])) - 0.0020) <= 1e-4
	@test abs(pvalue(SignedRankTest([1:10;], [2:11;])) - 0.0020) <= 1e-4
	@test default_tail(SignedRankTest([1:10;], [2:2:20;])) == :both
	show(IOBuffer(), SignedRankTest([1:10;], [2:2:20;]))
end

@testset "One Sample tests" begin
	# P-value computed using R wilcox.test
	@test abs(pvalue(SignedRankTest([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15] .- 10.1)) - 0.09460449) <= 1e-4
	# P-value computed using R wilcox.test
	@test abs(pvalue(SignedRankTest([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16] .- 10.1)) - 0.1928101) <= 1e-4
end

@testset "One Sample tests with ties" begin
	# P-value computed using R package exactRankTests wilcox.exact
	@test abs(pvalue(SignedRankTest([1,2,3,4,5,6,7,10,10,10,10,10,13,14,15] .- 10.1)) - 0.04052734) <= 1e-4
	# P-value computed using R wilcox.test
	@test abs(pvalue(SignedRankTest([1,2,3,4,5,6,7,10,10,10,10,10,13,14,15,16] .- 10.1)) - 0.1021964) <= 1e-4
end

@testset "Issue 128" begin
	@test pvalue(SignedRankTest([54.5, 54.5, 95.0, 51.5]), tail=:left)  == 1
	@test pvalue(SignedRankTest([54.5, 54.5, 95.0, 51.5]), tail=:right) == 0.0625
end

@testset "Test confidence interval" begin
x = [-7.8, -6.9, -4.7, 3.7, 6.5, 8.7, 9.1, 10.1, 10.8, 13.6, 14.4, 16.6, 20.2, 22.4, 23.5]
	@test isapprox(confint(ExactSignedRankTest(x))[1], 3.3, atol=1e-4)
	@test isapprox(confint(ExactSignedRankTest(x))[2], 15.5, atol=1e-4)
	@test isapprox(confint(ApproximateSignedRankTest(x))[1], 3.3, atol=1e-4)
	@test isapprox(confint(ApproximateSignedRankTest(x))[2], 15.5, atol=1e-4)
	@test isapprox(confint(SignedRankTest(x); tail=:left)[1], 4.45, atol=1e-4)
	@test isapprox(confint(SignedRankTest(x); tail=:right)[2], 14.45, atol=1e-4)
end
end


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 0a774706-1d0e-11eb-1d3a-6d641500d34d
using Pkg, DrWatson

# ╔═╡ 0a7776e0-1d0e-11eb-2a05-795e2470b260
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using PlutoUI
	using Plots
end

# ╔═╡ bc93b1e6-1d0d-11eb-3a1b-cf047fc8d621
md"## Listing 1.1"

# ╔═╡ 0a77d522-1d0e-11eb-1dfa-af3a546454ae
"There is more than one way to say hello:"

# ╔═╡ 0a83a6a4-1d0e-11eb-27cf-6b1ae5c75d6c
md"##### This is an array consisting of three strings."

# ╔═╡ 5cce7a9c-1d0e-11eb-0a59-43dceb80ee34
helloArray = ["Hello","G'day","Shalom"]

# ╔═╡ 0970cc02-1d35-11eb-0636-4d61aba864fe
md"##### Printing intermediate results in Pluto notebooks is a bit tedious currently, e.g.:"

# ╔═╡ ae5f130c-1d32-11eb-0322-4387a35fcae7
map(1:length(helloArray)) do i
	gr = helloArray[i]
    md" $gr world"
end

# ╔═╡ 33ffe1f2-2688-11eb-071d-ef4974b52b73
helloArray[1], "universe"

# ╔═╡ 5eee59a6-1d35-11eb-0e55-e53b9fef1b9c
md"##### Note the space in front of the first interpolation ( `$gr` )."

# ╔═╡ 5a9ab6da-1d33-11eb-27b2-a3ec1b2d8e02
map(1:3) do i
    R = i+100
    T = R^2
    md"hello $R and $T"
end

# ╔═╡ 988eab90-1d33-11eb-1278-0f39acc4dcd0
map(1:3) do i
    R = i+100
    T = R^2
    Text("hello $R and $T")
end

# ╔═╡ 0a8f6b9c-1d0e-11eb-0109-21ef005e63fe
md"##### Below we use a construct called a `comprehension' (or 'list comprehension')"

# ╔═╡ 0a8eef00-1d0e-11eb-11b7-d566161ac295
md"##### These squares are just perfect:"

# ╔═╡ 3d114318-1d10-11eb-2b93-ed04fda820ca
squares = [i^2 for i in 0:10]

# ╔═╡ 5eaa07b2-2689-11eb-23ce-55972c6a2861
md"##### Or, the above printing example:"

# ╔═╡ 8c1e52ee-2688-11eb-3fc7-ad83de5ec8c7
[("$(helloArray[i]) universe") for i in 1:length(helloArray)]

# ╔═╡ 0a9a2898-1d0e-11eb-1583-3bfdf887a967
md"##### You can loop on elements of arrays without having to use indexing. The last line of every code snippet is also evaluated as output (in addition to any figures and printing output shown previously)."

# ╔═╡ a0e7ab92-1d22-11eb-27fd-ef499034dd88
map(squares) do s
	sqrt(s)
end

# ╔═╡ 8b5571fa-1d35-11eb-1fef-1749f6da2d5e
md"##### Another important concept in Julia is `broadcasting` over an Array:"

# ╔═╡ 464824a8-1d0e-11eb-08b0-39b95ed0c30c
sqrt.(squares)

# ╔═╡ aca76098-1d35-11eb-1cf3-1b68437e17f7
md"##### Finally, the most tricky difference between the Julia REPL and other notebook formats is the reactivity of Pluto. It's cool, but also prevents multiple definitions across cells. If this happens, Pluto will suggest to bracket the assigments in a single cell with `begin` and `end`"

# ╔═╡ 01420810-1d36-11eb-1846-b579d2238fcc
begin
	x = 5
	x = x^2
end

# ╔═╡ 9e47c86e-1d38-11eb-0dae-ed14a8c961ad
md"##### Example of Pluto reactivity."

# ╔═╡ 88c09530-1d37-11eb-1b90-1991d07df0c2
@bind y Slider(3:12; default=7, show_value=true)

# ╔═╡ eaf674e0-1d37-11eb-234e-cbd28ef92cda
y^2

# ╔═╡ 569dad26-1d38-11eb-3a0b-69598098bf7c
plot(rand(y^3), lab="y^3 samples", leg=:topleft)

# ╔═╡ 02efc864-1d0e-11eb-027c-2f992de0069e
md"## End of listing 1. 1"

# ╔═╡ Cell order:
# ╟─bc93b1e6-1d0d-11eb-3a1b-cf047fc8d621
# ╠═0a774706-1d0e-11eb-1d3a-6d641500d34d
# ╠═0a7776e0-1d0e-11eb-2a05-795e2470b260
# ╠═0a77d522-1d0e-11eb-1dfa-af3a546454ae
# ╟─0a83a6a4-1d0e-11eb-27cf-6b1ae5c75d6c
# ╠═5cce7a9c-1d0e-11eb-0a59-43dceb80ee34
# ╟─0970cc02-1d35-11eb-0636-4d61aba864fe
# ╠═ae5f130c-1d32-11eb-0322-4387a35fcae7
# ╠═33ffe1f2-2688-11eb-071d-ef4974b52b73
# ╟─5eee59a6-1d35-11eb-0e55-e53b9fef1b9c
# ╠═5a9ab6da-1d33-11eb-27b2-a3ec1b2d8e02
# ╠═988eab90-1d33-11eb-1278-0f39acc4dcd0
# ╟─0a8f6b9c-1d0e-11eb-0109-21ef005e63fe
# ╟─0a8eef00-1d0e-11eb-11b7-d566161ac295
# ╠═3d114318-1d10-11eb-2b93-ed04fda820ca
# ╟─5eaa07b2-2689-11eb-23ce-55972c6a2861
# ╠═8c1e52ee-2688-11eb-3fc7-ad83de5ec8c7
# ╟─0a9a2898-1d0e-11eb-1583-3bfdf887a967
# ╠═a0e7ab92-1d22-11eb-27fd-ef499034dd88
# ╟─8b5571fa-1d35-11eb-1fef-1749f6da2d5e
# ╠═464824a8-1d0e-11eb-08b0-39b95ed0c30c
# ╟─aca76098-1d35-11eb-1cf3-1b68437e17f7
# ╠═01420810-1d36-11eb-1846-b579d2238fcc
# ╟─9e47c86e-1d38-11eb-0dae-ed14a8c961ad
# ╠═88c09530-1d37-11eb-1b90-1991d07df0c2
# ╠═eaf674e0-1d37-11eb-234e-cbd28ef92cda
# ╠═569dad26-1d38-11eb-3a0b-69598098bf7c
# ╟─02efc864-1d0e-11eb-027c-2f992de0069e


text:
# Helper function to plot results for the
# deterministic covariance estimation problem
using Plots
using JLD

#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(1.5)

@load "./out/exp2/cvest_d_60_maxIter_100000.jld"
NUM_GRADS = 100000
d = 40

cvgPlot = plot(xlim=(0:NUM_GRADS),
                 xlabel="Number of forward passes",
                 ylabel="\$\\log_{10}(\\|V V^{T} - \\Sigma\\|_{Fro}^2 + \\|W W^{T}\\|_{Fro}^2)/2\$",
                 legend=:right)

  for (sgrads, errs, η) in zip(sgradsCGDA, errsCGDA, ηList)
    mxg = findlast(x -> ~isnan(x), errs)
    if findfirst(isnan, errs) == nothing
      plot!(cvgPlot, sgrads[1:mxg], log10.(errs[1:mxg]),
            linestyle=:solid, linecolor=cdict[η],
            label="CGD, \\eta = $η")
    end
  end

  for (sgrads, errs, η) in zip(sgradsOGDA, errsOGDA, ηList)
    mxg = findlast(x -> ~isnan(x), errs)
    if findfirst(isnan, errs) == nothing
      plot!(cvgPlot, sgrads[1:mxg], log10.(errs[1:mxg]),
            linestyle=:dot, linecolor=cdict[η],
            label="OGDA, \\eta = $η")
    end
  end

  for (sgrads, errs, η) in zip(sgradsSGA, errsSGA, ηList)
    mxg = findlast(x -> ~isnan(x), errs)
    if findfirst(isnan, errs) == nothing
      plot!(cvgPlot, sgrads[1:mxg], log10.(errs[1:mxg]),
            linestyle=:dash, linecolor=cdict[η],
            label="SGA, \\eta = $η")
    end
  end

  for (sgrads, errs, η) in zip(sgradsConOpt, errsConOpt, ηList)
    mxg = findlast(x -> ~isnan(x), errs)
    if findfirst(isnan, errs) == nothing
      plot!(cvgPlot, sgrads[1:mxg], log10.(errs[1:mxg]),
            linestyle=:dashdot, linecolor=cdict[η],
            label="ConOpt, \\eta = $η")
    end
  end

  savefig(cvgPlot, "./out/cvest_d_$(d)")


text:
a, b = parse.(Int, split(readline()))
println(a*b, " ", 2a+2b)


text:
using irtfun
using Test

@testset "irtfun.jl" begin
    # Write your tests here.
end


text:
########################################################################################
# Structure following https://github.com/SciML/SimpleDiffEq.jl/tree/master/src/rk4
########################################################################################

########################################################################################
# Algorithm
########################################################################################
abstract type AbstractODEFilter <: DiffEqBase.AbstractODEAlgorithm end
mutable struct ODEFilter <: AbstractODEFilter end
mutable struct EKF0 <: AbstractODEFilter end
mutable struct EKF1 <: AbstractODEFilter end

########################################################################################
# Options
########################################################################################
mutable struct DEOptions{MI, absType, relType, QT, F1, tType, F2}
    maxiters::MI
    adaptive::Bool
    abstol::absType
    reltol::relType
    gamma::QT
    qmin::QT
    qmax::QT
    beta1::QT
    beta2::QT
    qoldinit::QT
    internalnorm::F1
    unstable_check::F2
    dtmin::tType
    dtmax::tType
    force_dtmin::Bool
    verbose::Bool
end

########################################################################################
# Integrator
########################################################################################
mutable struct ODEFilterIntegrator{
    IIP, S, T, P, F, QT, O, constantsType, cacheType,
    sigmaestType, stepruleType, proposalsType,
    xType, sigmaType, probType, algType} <: DiffEqBase.AbstractODEIntegrator{algType, IIP, S, T}
    sol
    f::F                               # eom
    u::S                               # current functionvalue
    # x::X                             # current state
    t::T                               # current time
    t_new::T                           # next time
    t0::T                              # initial time, only for reinit
    tmax::T                            # end of the interval
    dt::T                              # step size
    p::P                               # parameter container
    EEst::QT                           # (Scaled) error estimate
    qold::QT                           #

    constants::constantsType
    cache::cacheType

    # Options
    opts::O                            # General (not PN-specific) solver options
    sigma_estimator::sigmaestType
    steprule::stepruleType
    smooth::Bool                       # Smooth the solution or not

    # Save into
    proposals::proposalsType           # List of proposals
    state_estimates::Vector{xType}     # List of state estimates, used to build the solution
    times::Vector{T}
    sigmas::Vector{sigmaType}

    # Misc
    iter::UInt                         # Current iteration count
    success_iter::UInt                 # Number of successful steps
    accept_step::Bool                  # If the current step is accepted
    retcode::Symbol                    # Current return code, used to build the solution
    prob::probType                     # Only used to build the solution
    alg::algType                       # Only used to build the solution
    destats::DiffEqBase.DEStats        # To track stats like the number of f evals
end
DiffEqBase.isinplace(::ODEFilterIntegrator{IIP}) where {IIP} = IIP

########################################################################################
# Caches
########################################################################################
abstract type ProbNumODECache <: DiffEqBase.DECache end
abstract type ProbNumODEConstantCache <: ProbNumODECache end
abstract type ProbNumODEMutableCache <: ProbNumODECache end
struct GaussianODEFilterConstantCache{RType, EType, F1, F2} <: ProbNumODEMutableCache
    d::Int                  # Dimension of the problem
    q::Int                  # Order of the prior
    # dm                    # Dynamics Model
    # mm                    # Measurement Model
    A!
    Q!
    h!
    H!
    R::RType
    # Precond
    # Precond_inv
    E0::EType
    E1::EType
    jac
    Precond::F1
    InvPrecond::F2
end
function GaussianODEFilterConstantCache(prob, q, prior, method)
    d, f = length(prob.u0), prob.f
    @assert prior == :ibm
    Precond, InvPrecond = preconditioner(d, q)
    A!, Q! = ibm(d, q)

    E0 = kron([i==1 ? 1 : 0 for i in 1:q+1]', diagm(0 => ones(d)))
    # E0 = E0 * InvPrecond
    E1 = kron([i==2 ? 1 : 0 for i in 1:q+1]', diagm(0 => ones(d)))
    # E1 = E1 * InvPrecond

    @assert method in (:ekf0, :ekf1) ("Type of measurement model not in [:ekf0, :ekf1]")
    jac = method == :ekf1 ? f.jac : nothing
    h!(h, du, m) = h .= E1*m - du
    H!(H, ddu) = H .= E1 - ddu * E0
    # R = diagm(0 => eps(eltype(prob.u0)) .* ones(d))
    R = zeros(d, d)

    return GaussianODEFilterConstantCache{typeof(R), typeof(E0), typeof(Precond), typeof(InvPrecond)}(
        d, q, A!, Q!, h!, H!, R, E0, E1, jac, Precond, InvPrecond)
end

mutable struct GaussianODEFilterCache{uType, xType, matType, sigmaType} <: ProbNumODEMutableCache
    u::uType
    u_pred::uType
    u_filt::uType
    u_tmp::uType
    # tmp
    x::xType
    x_pred::xType
    x_filt::xType
    x_tmp::xType
    measurement
    Ah::matType
    Qh::matType
    h::uType
    H::matType
    du::uType
    ddu::matType
    K::matType
    σ_sq::sigmaType
    σ_sq_prev::sigmaType
    P_tmp::matType
    err_tmp::uType
end
function GaussianODEFilterCache(d, q, prob, constants, σ0, initialize_derivatives=true)
    @unpack E0, E1 = constants
    u0 = prob.u0

    uType = typeof(u0)
    uElType = eltype(u0)
    matType = Matrix{uElType}

    t0 = prob.tspan[1]
    # Initial states
    m0, P0 = initialize_derivatives ?
        get_initial_states_forwarddiff(prob, q) :
        initialize_without_derivatives(prob, q)
    # P0 += eps(eltype(P0)) * I
    x0 = Gaussian(m0, P0)

    Ah_empty = diagm(0=>ones(uElType, d*(q+1)))
    Qh_empty = zeros(uElType, d*(q+1), d*(q+1))
    h = E1 * x0.μ
    H = uElType.(zeros(d, d*(q+1)))
    du = copy(h)
    ddu = uElType.(zeros(d, d))
    sigma = σ0
    v, S = copy(h), copy(ddu)
    measurement = Gaussian(v, S)
    K = copy(H')

    GaussianODEFilterCache{uType, typeof(x0), matType, typeof(sigma)}(
        copy(u0), copy(u0), copy(u0), copy(u0),
        copy(x0), copy(x0), copy(x0), copy(x0),
        measurement,
        Ah_empty, Qh_empty, h, H, du, ddu, K, sigma, sigma,
        copy(P0),
        copy(u0),
    )

end

function initialize_without_derivatives(prob, order, var=1e-3)
    q = order
    u0 = prob.u0
    d = length(u0)
    p = prob.p
    t0 = prob.tspan[1]

    m0 = zeros(d*(q+1))
    m0[1:d] = u0
    if !isinplace(prob)
        m0[d+1:2d] = prob.f(u0, p, t0)
    else
        prob.f(m0[d+1:2d], u0, p, t0)
    end
    P0 = [zeros(d, d) zeros(d, d*q);
          zeros(d*q, d) diagm(0 => var .* ones(d*q))]
    return m0, P0
end

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
function BoundedUntil(imdp, phi1, phi2, k, imdp_filepath; synthesis_flag=false)

    labels_vector = Array{String}(undef, length(imdp.states))
    for label_key in keys(imdp.labels)
        labels_vector[label_key] = imdp.labels[label_key]
    end

    # Get the Qyes and Qno states
    Qyes = findall(x->x==phi2, labels_vector)
    @info Qyes
    Qno = isnothing(phi1) ? [] : findall(x->!(x==phi1 || x==phi2), labels_vector)
    @info Qno

    # Write them to the file
    write_imdp_to_file_bounded(imdp, Qyes, Qno, imdp_filepath)

    mode1 = synthesis_flag ? "maximize" : "minimize"
    # Do verification
    if !isnothing(phi1)
        result_mat = run_bounded_imdp_verification(imdp_filepath, k, mode1=mode1)
    elseif synthesis_flag
        result_mat = run_imdp_synthesis(imdp_filepath, k; ep=1e-6, mode1="maximize", mode2="pessimistic", save_mats=true)
    else
        result_mat = run_imdp_synthesis(imdp_filepath, k; ep=1e-6, mode1="minimize", mode2="pessimistic", save_mats=true)
    end
        # Done-zoe
    return result_mat
end

function Globally(imdp, phi, k, imdp_filepath; synthesis_flag=false)

    phi1 = nothing
    phi2 = "!$phi"

    result_mat = BoundedUntil(imdp, phi1, phi2, k, imdp_filepath, synthesis_flag=synthesis_flag)
    safety_result_mat = zeros(size(result_mat))
    safety_result_mat[:, 1] = result_mat[:, 1]
    safety_result_mat[:, 2] = result_mat[:, 2]
    safety_result_mat[:, 3] = 1 .- result_mat[:, 4]
    safety_result_mat[:, 4] = 1 .- result_mat[:, 3]
    return safety_result_mat
end

"""
Call the synthesis tool with the given parameters.
"""
function run_imdp_synthesis(imdp_file, k; ep=1e-6, mode1="maximize", mode2="pessimistic", save_mats=true)
    exe_path = "/usr/local/bin/synthesis"  # Assumes that this program is on the user's path
    @assert isfile(imdp_file)
    res = read(`$exe_path $mode1 $mode2 $k 0.000001 $imdp_file`, String)
    filter_res = replace(res, "\n"=>" ")
    res_split = split(filter_res)
    dst_dir = dirname(imdp_file)
    open("$dst_dir/verification-result-$k.txt", "w") do f
        print(f, res)
    end
    print(res)
    res_mat = res_to_numbers(res)
    return res_mat
end

"""
Call the synthesis tool with unbounded time horizon.
"""
function run_unbounded_imdp_synthesis(imdp_file)
    res_mat = run_imdp_synthesis(imdp_file, -1)
    return res_mat
end

"""
Call the synthesis tool to perform verification on a bounded horizon.
"""
function run_bounded_imdp_verification(imdp_file, k; mode1="minimize")
    res_mat = run_imdp_synthesis(imdp_file, k, mode1=mode1, mode2="pessimistic")
    return res_mat
end

"""
Call the MATLAB verification tool.
"""
function run_MATLAB_verification(script_path, exp_dir, k, verification_mode)
    @info "Performing BMDP Verification..."
    # mat"cd($script_path)"
    # mat"addpath(genpath('./'))"
    # mat"generateVerification($exp_dir,$verification_mode,$k)"
end

"""
Call the synthesis tool to perform verification on an unbounded horizon.
"""
function run_unbounded_imdp_verification(imdp_file)
    res_mat = run_imdp_synthesis(imdp_file, -1, mode1="minimize", mode2="pessimistic")
    return res_mat
end

"""
Confirm that the verification of a data-driven system is consistent with the known system.
"""
function verify_experiment(path_to_experiment, path_to_truth)

    @info "Verifying experiment at $path_to_experiment"
    # load verification matrices
    exp_mats = MAT.matread("$path_to_experiment/globally-safe-1-steps.mat")
    # load truth matrices
    true_mats = MAT.matread("$path_to_truth/globally-safe-1-steps.mat")

    if sum(exp_mats["indVmin"] .> true_mats["indVmin"]) != 0
        @error "Error in the following regions: " findall(x->x, exp_mats["indVmin"] .> true_mats["indVmin"])
    end
    @assert sum(exp_mats["indVmax"] .< true_mats["indVmax"]) == 0

end

"""
Creates a matrix from the string output of the BMDP synthesis tool.
"""
function res_to_numbers(res_string)
    filter_res = replace(res_string, "\n"=>" ")
    res_split = split(filter_res)
    num_rows = Int(length(res_split)/4)

    res_mat = zeros(num_rows, 4)
    for i=1:num_rows
        res_mat[i, 1] = parse(Int, res_split[(i-1)*4+1])+1.
        res_mat[i, 2] = parse(Int, res_split[(i-1)*4+2])+1.
        res_mat[i, 3] = parse(Float64, res_split[(i-1)*4+3])
        res_mat[i, 4] = parse(Float64, res_split[(i-1)*4+4])
    end

    return res_mat
end

"""
Saves the output of the synthesis tool in a legacy way.
"""
function save_legacy_mats(res_mat, dst_dir, k)
    policy = res_mat[:, 2]
    indVmax = res_mat[:, 4]
    indVmin = res_mat[:, 3]

    filename = "$dst_dir/globally-safe-$k-steps.mat"
    matwrite(filename, Dict("indVmin" => indVmin, "indVmax" => indVmax, "policy" => policy))
end

"""
Assign a label to the extent based on its center.
"""
function general_label_fcn(state_extent, default_label, unsafe_label, labels_dict; unsafe=false)
    if unsafe
        return unsafe_label
    end
    state_label = default_label
    for label in keys(labels_dict)
        for extent in labels_dict[label]
            flags = [sum(state_extent[dim])/2∈extent[dim] for dim in keys(extent)]
            if sum(flags) == extent.count
                state_label = label
                @debug state_label
                break
            end
        end
    end
    return state_label
end


text:
module TracingTests #end

include("../examples/Tracing.jl")
using .Tracing

@static if VERSION >= v"1.3-"
    import Base.Threads: @spawn
else
    macro spawn(e) esc(:(@async(e))) end
end

using Test

@testset "Tracing example" begin

    edges = tracetasks() do
        @sync begin
            @async begin
                # Embedded trace within the parent trace
                inneredges = @tracetasks 2+2
                @assert length(inneredges) == 1
            end
            Threads.@spawn begin
            end
        end
    end

    @test length(edges) == 3

end

end  # module


text:
# Clifford module

# This module contains code for generating sequences of clifford gates.  The
# exported functions benchmark1Qubit and benchmark2Qubit generate sequences
# equivalent to an idle gate for one qubit, and to either an idle or swap for
# two qubits.  Additional arguments allow for restricting the control sequences
# to only XY control and for interleaved random benchmarking of a specific gate.

#============= Matrices =============#

# Mathematical background for using integers mod 5 instead of complexes:

# 2^n x 2^n Clifford matrices can be normalized to have all entries one of
# {0, +1, -1, +i, -i}, and the first nonzero entry in the first row may always
# be 1.  This is a 5-element set closed under multiplication, so we can just use
# the 5-element field. Explicitly: 0 <=> 0, 1 <=> 1, -1 <=> 4, i <=> 2, -i <=> 3
# respects element multiplication.  Matrix multiplication involves addition as
# well, and while 1 + 2 = 3, 1 + i is not -i.  However, under the normalization
# convention above, the entire matrix is divided by the first nonzero entry, or
# equivalently multiplied by its inverse.  Thus, one need only verify that all
# encodings of {1,2,3,4} with magnitudes sqrt(2)^k for integer k have the same
# angular orientation (that is, they read 1, 2, 4, 3 taken clockwise).  This
# may be easily proven graphically, showing that this isomorphism is sound.

# The clifford matrix is 2^2n elements and all are part of the 5-element field,
# so in principle an enitre 4x4 matrix is one UInt64's worth of information.
# In practice, that kind of bit twiddling is likely slower than machine
# vectorized integer operations, so instead we will just use matrices
# of UInt8s, taking things mod 5 at the end of a calculation.

module Clifford

export benchmark1Qubit
export benchmark2Qubit
export Pulse
export gateNames

import Base: convert
import Base: getindex, size
import Base: show, showarray, summary
import Base: normalize, *, inv, kron

using StaticArrays
using Core.Intrinsics: box, unbox

immutable CMatrix{S,L} <: StaticMatrix{UInt8}
    data::NTuple{L,UInt8}

    function CMatrix(d::NTuple{L,UInt8})
        StaticArrays.check_smatrix_params(Val{S}, Val{S}, UInt8, Val{L})
        new(d)
    end

    function CMatrix(d::NTuple{L})
        StaticArrays.check_smatrix_params(Val{S}, Val{S}, UInt8, Val{L})
        new(StaticArrays.convert_ntuple(UInt8, d))
    end
end

function CMatrix(x::Tuple)
    L = length(x)
    S = sqrt(L)
    if !isinteger(S)
        error("Must provide a square number of elements.")
    end
    CMatrix{Int(S),L}(x)
end

@generated function (::Type{CMatrix{S1}}){S1,L}(x::NTuple{L})
    if S1*S1 != L
        error("Incorrect matrix size: $S1 * $S1 != $L.")
    end
    return quote
        $(Expr(:meta, :inline))
        CMatrix{S1, L}(x)
    end
end

# StaticArrays interfacing
@inline size{S,L}(::Union{CMatrix{S,L}, Type{CMatrix{S,L}}}) = (S,S)
@inline getindex(A::CMatrix, i::Integer) = A.data[i]

# Here's how we are going to display these matrices in a meaningful way.
# Unnormalized matrices result in a delicious slice of pizza.
bitstype 8 CEntry  # has to be a multiple of 8 at the moment
convert(::Type{CEntry}, x::CEntry) = x
convert(::Type{CEntry}, x::Number) = box(CEntry, unbox(UInt8, UInt8(x)))
convert{T<:Number}(::Type{T}, x::CEntry) = T(box(UInt8, unbox(CEntry, x)))
function show(io::IO, c::CEntry)
    v = UInt8(c)
    if v == 0x00
        print(io, " 0")
    elseif v == 0x01
        print(io, " 1")
    elseif v == 0x02
        print(io, " 𝒊")
    elseif v == 0x03
        print(io, "-𝒊")
    elseif v == 0x04
        print(io, "-1")
    else
        print(io, " 🍕")
    end
end

show{S}(io::IO, ::MIME"text/plain", x::CMatrix{S}) =
    showarray(io, SMatrix{S,S,CEntry}(x.data), false)
summary{S}(m::SMatrix{S,S,CEntry}) = "Clifford matrix"

# Our normalization convention is that the first nonzero entry be a 1.  This
# normalization reconciles the 5-element field and the +1, -1, +i, -i
# representations, showing both are isomorphic to the Clifford group.
# NOTE: normalize only produces correct results on matrices already modded by 5.
function normalize(A::CMatrix)
  normConst = A[findfirst(A)]
  inverse = normConst $ ((normConst >>> 0x01) & 0x01) # Inverse mod 5
  mod.(A*inverse, 0x05)
end

# Implement multiplication such that the matrix product is properly normalized.
# Note that this function throws an error if called with no arguments, because
# it cannot infer what dimensions to use.
function *(A::CMatrix, B::CMatrix...)
    mult(X::CMatrix, Y::CMatrix) =
        invoke(*, (StaticMatrix{UInt8}, StaticMatrix{UInt8}), X, Y)
    normalize(reduce((D,C) -> mod.(mult(D,C), 0x05), eye(A), (A,B...)))
end

# We also have to modify matrix inversion.  My present choice is to compute the
# standard inverse and multiply it by the determinant to get a strictly-integral
# matrix.  Then we normalize.  More direct routes could be taken to get the
# minors of each element but I feel it would amount to rewriting too much for
# something not performance-critical.
function inv{S}(A::CMatrix{S})
    B = SMatrix{S,S,Float64}(A) # avoid problem in StaticArrays.jl
    normalize(CMatrix{S,S*S}(map(x -> UInt8(round(mod.(x, 0x05))), inv(B) * det(B))))
end

# Output a clifford matrix corresponding to A on qubit 1 and B on qubit 2.
function kron{S,T}(A::CMatrix{S}, B::CMatrix{T})
    #TODO: implement kron for StaticArrays in a pull request
    r = invoke(kron, (StaticMatrix{UInt8}, StaticMatrix{UInt8}), A, B)
    A = S*T
    B = A*A
    CMatrix{A,B}(mod(r, 0x05))
end

# I'll explicitly write out the x pi/2 and y pi/2 single-qubit matrices, using
# these to generate the rest.  I will refer to them as f and g, respectively.
const f = CMatrix{2,4}([1 2; 2 1])
const g = CMatrix{2,4}([1 4; 1 1])
const h = *(g, f, g, g, g) # Gives the z/2 rotation

# For later use, the 2-qubit matrices cz and swap
const cz = CMatrix{4,16}([1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 4])
const swap = CMatrix{4,16}([1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1])

# Single qubit clifford group is 24 gates, I'll just do this the long way.
# In the Mathematica notebook these are shown to correspond 1:1 with the pulses
# listed below (keepiing in mind that pulses are time-ordered left to right,
# while matrices are applied right to left)
const SQClif = CMatrix[eye(f),
               f,
               g,
               h,
               *(f, f),
               *(g, g),
               *(h, h),
               *(f, f, f),
               *(g, g, g),
               *(h, h, h),
               *(f, g),
               *(f, g, g, g),
               *(f, f, f, g),
               *(f, f, f, g, g, g),
               *(f, h),
               *(f, h, h, h),
               *(f, f, f, h),
               *(f, f, f, h, h, h),
               *(f, f, g),
               *(f, f, g, g, g),
               *(g, g, h),
               *(g, g, h, h, h),
               *(h, h, f),
               *(h, h, f, f, f)]

# 2-qubit clifford group can be constructed as a SQClif gate on each qubit,
# followed by one of the ten operations below, followed by an optional swap.
# We use only the "easy" half of the clifford group, modding out by swaps.
const EntanglingGate = CMatrix[eye(cz),
                       cz,
                       *(kron(eye(f), f), cz),
                       *(kron(eye(f), g), cz),
                       *(kron(f, eye(f)), cz),
                       *(kron(f, f), cz),
                       *(kron(f, g), cz),
                       *(kron(g, eye(f)), cz),
                       *(kron(g, f), cz),
                       *(kron(g, g), cz)]

#============= Pulses =============#

# We need an understandable way to designate the gates, so that control lines
# can quickly convert between this program's output and pulse sequencing.  A
# pulse matrix has each column a different qubit, and each row a different time
# step.  Functions closer to the instruments will convert the entries for a
# specific qubit into microwave control pulses, while it is the responsibility
# of the higher-up code to ensure that the requested pulses are all timed
# appropriately.

# Initially Pulse was simply a typealias to Matrix{Int8}, but I wanted to treat
# Pulses differently in some special cases for dispatch.  The definitions below
# are sufficient to allow Julia to treat Pulse objects identically to Matrices
# for indexing and concatenation, and show Pulses nicely with named gates.  It
# also allows one to specify a sequence with one pulse repeated many times, for
# example [Xpi; Idle*10; Ypi] for a pause between Xpi and Ypi of 10 Idle gates.
import Base: size, getindex, setindex!, linearindexing, similar, showarray, *
type Pulse <: AbstractMatrix{Int8}
    m::Matrix{Int8}
end
size(x::Pulse) = size(x.m)
getindex(x::Pulse, i::Int) = getindex(x.m, i)
getindex(x::Pulse, i::Int...) = getindex(x.m, i...)
linearindexing(::Type{Pulse}) = Base.LinearFast()
getindex(x::Pulse, I...) = getindex(x.m, I...)
setindex!(x::Pulse, v, i::Int) = x.m[i] = v
setindex!(x::Pulse, v, I::Int...) = x.m[I...] = v
setindex!(x::Pulse, v, I...) = x.m[I...] = v
similar{T}(x::Pulse, ::Type{T}, dims::Dims) = Pulse(Array{Int8}(dims))
*(x::Pulse, I::Integer) = vcat(fill(x, I)...)
*(I::Integer, x::Pulse) = vcat(fill(x, I)...)

# This scheme for getting pulses to display nicely in the REPL is hackish and
# unstable.  Hopefully in later releases it won't be necessary.  TODO: better..
showarray(io::IO, x::Pulse, repr::Bool) = prettyPrint(x)
function prettyPrint(x::Pulse)
  println("Length-"*string(size(x,1))*" pulse sequence:")
  [println(reduce(*,"",[gateNames[x[y,z]]*"\t" for z in 1:size(x,2)]))
                                               for y in 1:size(x,1)]
  nothing
end

# This is essentially an enum, but the numbers are chosen to have some later
# significance.  Note the null gate is empty rather than idle.
const nul = Pulse(Array(Int8,(0,1)))
const Idle =  Pulse(fill(Int8(7),(1,1)))
const Xpi2 =  Pulse(fill(Int8(1),(1,1)))
const Xpi =   Pulse(fill(Int8(5),(1,1)))
const X3pi2 = Pulse(fill(Int8(3),(1,1)))
const Ypi2 =  Pulse(fill(Int8(2),(1,1)))
const Ypi =   Pulse(fill(Int8(6),(1,1)))
const Y3pi2 = Pulse(fill(Int8(4),(1,1)))
const Zpi2 =  Pulse(fill(Int8(8),(1,1)))
const Zpi =   Pulse(fill(Int8(9),(1,1)))
const Z3pi2 = Pulse(fill(Int8(10),(1,1)))
const CZ =    Pulse(Int8[11 11])  # CZ requires 2 qubits to implement.

# While the code works with numbers, it is handy to see the gates as strings.
const gateNames = ["Xpi2", "Ypi2", "X3pi2", "Y3pi2", "Xpi", "Ypi", "Idle",
"Zpi2", "Zpi", "Z3pi2", "CZ"]

# There are two potentially desirable ways to perform clifford gates.  One is
# in the fewest number of gates, assuming XY and Z lines are present (WithZ).
# The other is in the fewest number of XY-only gates (NoZ).  In the WithZ case,
# the program will be able to do RB seqeunces with the same effect in fewer
# gates, and the Qubit type will ask for the channel for all 3 lines: XYI, XYQ
# and Z.  In the NoZ case, the gates performed will only exist on the
# XYI and XYQ lines, and the Qubit type will not ask about a Z line.  The latter
# is necessary for doing single qubit RB with only the AWG, since there are 4
# lines and 2 are used for readout.
const SQPulseWithZ = Pulse[nul,
      Xpi2,
      Ypi2,
      Zpi2,
      Xpi,
      Ypi,
      Zpi,
      X3pi2,
      Y3pi2,
      Z3pi2,
      [Ypi2; Xpi2],
      [Y3pi2; Xpi2],
      [Ypi2; X3pi2],
      [Y3pi2; X3pi2],
      [Zpi2; Xpi2],
      [Z3pi2; Xpi2],
      [Zpi2; X3pi2],
      [Z3pi2; X3pi2],
      [Ypi2; Xpi],
      [Y3pi2; Xpi],
      [Zpi2; Ypi],
      [Z3pi2; Ypi],
      [Xpi2; Zpi],
      [X3pi2; Zpi]]

const SQPulseNoZ = Pulse[nul,
      Xpi2,
      Ypi2,
      [Xpi2; Ypi2; X3pi2],
      Xpi,
      Ypi,
      [Xpi; Ypi],
      X3pi2,
      Y3pi2,
      [Xpi2; Y3pi2; X3pi2],
      [Ypi2; Xpi2],
      [Y3pi2; Xpi2],
      [Ypi2; X3pi2],
      [Y3pi2; X3pi2],
      [Xpi2; Ypi2],
      [Xpi2; Y3pi2],
      [X3pi2; Y3pi2],
      [X3pi2; Ypi2],
      [Ypi2; Xpi],
      [Y3pi2; Xpi],
      [Y3pi2; Xpi2; Y3pi2],
      [Ypi2; Xpi2; Ypi2],
      [Ypi; Xpi2],
      [Ypi; X3pi2]]

# With the above defined, we can make a lookup dictionary for the recovery pulse
SQLookupWithZ = Dict{CMatrix, Pulse}(map(=>, SQClif, SQPulseWithZ))
SQLookupNoZ = Dict{CMatrix, Pulse}(map(=>, SQClif, SQPulseNoZ))

# And with only slightly more effort, we can do the same for the 2-qubit case.
# We need to list the pulses for the entangling gates:
const EntanglingPulse = Pulse[[nul nul],
                   CZ,
                   [CZ; [Idle Xpi2]],
                   [CZ; [Idle Ypi2]],
                   [CZ; [Xpi2 Idle]],
                   [CZ; [Xpi2 Xpi2]],
                   [CZ; [Xpi2 Ypi2]],
                   [CZ; [Ypi2 Idle]],
                   [CZ; [Ypi2 Xpi2]],
                   [CZ; [Ypi2 Ypi2]]]

# Coordinate pulses of unknown lengths on two qubits in parallel.  This is in
# complement to "outer" defined above, in that if pulse A has matrix mA, B has
# matrix mB, then pulseOuter(A, B) has matrix given by outer(mA, mB).
function pulseOuter(A::Pulse, B::Pulse)
  if A == nul
    if B != nul
      A = Idle
    else
      return [A B]
    end
  elseif B == nul
    B = Idle
  end

  pad = size(A, 1) - size(B, 1)
  if pad == 2
    [A [B; Idle; Idle]]
  elseif pad == 1
    [A [B; Idle]]
  elseif pad == -1
    [[A; Idle] B]
  elseif pad == -2
    [[A; Idle; Idle] B]
  else
    [A B]
  end
end

# And then the analogous (5760-element) lists and lookup table.  Note that we
# only employ the WithZ case here, assuming qubits that are coupled need to be
# Z-tunable.  With a gmon scheme this might not be accurate.
const TQClif = collect([*(EntanglingGate[k], kron(SQClif[i], SQClif[j]))
          for i=1:24, j=1:24, k=1:10])
const TQPulse= collect([[pulseOuter(SQPulseWithZ[i], SQPulseWithZ[j]); EntanglingPulse[k]]
          for i=1:24, j=1:24, k=1:10])

# The lookup table finds the pulse corresponding to either the matrix given or
# that matrix plus a swap operation.  Since we use this lookup table to reach
# the ground state, both will work.  We will not be turning error states into
# clean ones, just switching which qubit fails to read ground state.
TQLookup = Dict{CMatrix, Pulse}(map(=>, [TQClif map(A->*(swap, A), TQClif)]
                                      , [TQPulse TQPulse]))

# From this we can build a function which chooses uniformly from the group,
# building sequences of nClifs operations ultimately ending in the identity.
# While we treat the 1-qubit and 2-qubit cases separately, it is mostly
# copy-paste.

# These functions take a gate to benchmark, given by its index in the SQClif or
# TQClif lists.  These functions are not meant to be called directly by users,
# who would not know this index.  Calling functions should find the index using
# find() or findfirst() and some more meaningful description of the gate, either
# the matrix in *Clif or the pulse sequence in *Pulse.  For a baseline, the null
# gate should be used, which has empty lists as pulses.  This is index 1.

# This leads to a bit of complication if we want to benchmark the idle gate,
# which explicitly idles for a short time instead of proceeding directly to the
# next clifford.  One possibility is to add to the gateset an Idle pulse which
# has the same matrix but a different pulse and different index.  It would be
# left out of the lookup dictionary.  TODO: add support for benchmarking "idle."

function benchmark1Qubit(nClifs, pulseIndex::Int = 1, ZControl::Bool = true)
  selection = rand(1:24, nClifs-1)
  recovery = inv(*([fill(SQClif[pulseIndex], (1,nClifs-1));
                     SQClif[reverse(selection)']][:]...))
  if ZControl
    [vcat(SQPulseWithZ[[selection fill(pulseIndex, nClifs-1)]'[:]]...); SQLookupWithZ[recovery]]
  else
    [vcat(SQPulseNoZ[[selection fill(pulseIndex, nClifs-1)]'[:]]...); SQLookupNoZ[recovery]]
  end
end

# A deterministic selection is specified for debugging
function altbenchmark1Qubit(nClifs, selection, pulseIndex::Int = 1, ZControl::Bool = true)
    recovery = inv(*([fill(SQClif[pulseIndex], (1,nClifs-1));
                       SQClif[reverse(selection)']][:]...))
    if ZControl
      [vcat(SQPulseWithZ[[selection fill(pulseIndex, nClifs-1)]'[:]]...); SQLookupWithZ[recovery]]
    else
      [vcat(SQPulseNoZ[[selection fill(pulseIndex, nClifs-1)]'[:]]...); SQLookupNoZ[recovery]]
    end
  end

# Convenience method allowing ommitted index but specified ZControl status
function benchmark1Qubit(nclifs, ZControl::Bool, pulseIndex::Int = 1)
  benchmark1Qubit(nclifs, pulseIndex, ZControl)
end

# 2 qubit gates cannot be done without Z control unless gmon-coupled, and this
# code doesn't presently support sending pulses to interqubit couplings.  As
# such, ZControl is assumed true always.
function benchmark2Qubit(nClifs, pulseIndex::Int = 1)
  selection = rand(1:5760, nClifs-1)
  recovery = inv(*([TQClif[reverse(selection)]
                          fill(TQClif[pulseIndex], nClifs-1)]'[:]...))
  [vcat(TQPulse[[selection fill(pulseIndex, nClifs-1)]'[:]]...); TQLookup[recovery]]
end

end # End module Clifford

max_stars_repo_name: PainterQubits/PulseController

max_issues_repo_name: PainterQubits/PulseController

max_forks_repo_name: PainterQubits/PulseController


text:
struct Pot
    tea_bag::Any
    pid::PotID
    require::ResourceInfo
    logger::Any
end

function Pot(
    tea_bag;
    name=string(uuid4()),
    cpu=eps(),
    gpu=0,
    logger=DEFAULT_LOGGER
)
    pid = name isa PotID ? name : PotID(name)
    require = ResourceInfo(cpu, gpu)
    Pot(tea_bag, pid, require, logger)
end

macro pot(tea, kw...)
    tea_bag = esc(:(() -> ($(tea))))
    xs = [esc(x) for x in kw]
    quote
        p = Pot($tea_bag; $(xs...))
        register(p)
        p.pid
    end
end

mutable struct PotState
    pid::PotID
    ch::Channel
    create_time::DateTime
    last_update::DateTime
    n_processed::UInt
end

_self() = get!(task_local_storage(), KEY, PotState(USER, current_task()))
self() = _self().pid

local_scheduler() = SCHEDULER/"local_scheduler_$(myid())"

Base.parent() = parent(self())
Base.parent(p::PotID) = PotID(getfield(p, :path[1:end-1]))

children() = children(self())


text:
function load_invalid_file(cur_package::Module, invalid_file::FileInfo)

  # read in file

  cur_file = invalid_file.name

  opened_file = open(cur_file)

  read_file = readstring(opened_file)

  seekstart(opened_file)

  # find shard in file

  for cur_shard in invalid_file.loaded_shards
    parse(opened_file)
  end

  beg_position = position(opened_file)

  iszero(beg_position) && ( beg_position = 1 )

  parse(opened_file)

  end_position = position(opened_file)

  close(opened_file)

  # build include string

  line_count = length(
    matchall(
      r"\n",
      read_file[1:beg_position-1]
    )
  )

  cur_string = "\n" ^ line_count

  cur_string *= read_file[beg_position:end_position]

  # make call to bad file

  cur_package.include_string(cur_string, cur_file)

end


text:
using Documenter, SpacetimeFields

makedocs(modules=[SpacetimeFields])


text:
defaultTo(x, default) = x ≠ nothing ? x : default

walk(f, leaf) = f(leaf)
function walk(f, ex::Expr)
    new_ex = copy(ex)
    new_ex.args = [walk(f, x) for x ∈ new_ex.args]
    new_ex
end

contains(x, ex) = x == ex
function contains(x, ex::Expr)
    found = false
    walk(y -> found |= y==x, ex)
    found
end

function partition(xs, pred)
    a, b = [], []
    for x ∈ xs
        push!(pred(x) ? a : b, x)
    end
    a, b
end


text:
using Test, LightGraphs, Parsers, LightGraphsIO
using LightGraphs.Experimental

@testset "DOT tests" begin
    g = SimpleGraph(10,10)
    dg = SimpleDiGraph(10,10)

    read_g = load_graph("graph.dot", DOTFormat)
    read_dg = load_graph("digraph.dot", DOTFormat)

    # Single graph read/write
    @test_broken has_isomorph(g, read_g)
    @test_broken has_isomorph(dg, read_dg)

    # Multiple graphs read/write
    #
end


text:
type DynamicFactorModel
    # A dynamic factor model is the same as a factor model but it includes time dependence structure which we model by using lags of the static factor in the forecasting equation
    factor_model::FactorModel
    number_of_factor_lags::Int64  # lag level to use for forecasting equation

    function DynamicFactorModel(factor_model_args::Tuple)
        factor_model = FactorModel(factor_model_args...)
        number_of_factor_lags = estimate_number_of_lags(factor_model)  # if number of factor lags are not given we estimate them
        warn("Automatic recognition of the number of factor lags is not implemented yet")
        # TODO: this is undone! number_of_factor_lags has to be estimated from frequency domain transformation
        new(factor_model)
    end
    function DynamicFactorModel(factor_model_args::Tuple, number_of_factor_lags::Int64)
        factor_model = FactorModel(factor_model_args...)
        new(factor_model, number_of_factor_lags)
    end
end

# transforms x to the space spanned by the factors and optionally only selects active factors
#   type="active" returns only the active factors (which explain enough of the variance)
function get_factors(dfm::DynamicFactorModel, x::Matrix, factors="active")
    # TODO: this function makes no sense like this, there are no break indices in the model anymore
    break_indices = [0, dfm.break_indices, size(dfm.x, 1)]
    rotations = [loadings * inv(loadings'loadings) for loadings in dfm.loadings] # simplifies to loadings if T>N and loadings = inv(loadings')  which makes sense given x = F * loadings'
    [(normalize(x, (mean(dfm.factor_model.x), std(dfm.factor_model.x)))*rotations[i-1])[break_indices[i-1]+1:break_indices[i], factors=="active" ? (1:dfm.number_of_factors) : (1:end)] for i in 2:length(break_indices)]
end

function make_factor_model_design_matrix(y, number_of_lags, factors, number_of_factors, number_of_factor_lags, break_indices)
    w = lagged_matrix(y, Int64[1:number_of_lags])
    return(w[:, 1], hcat(w[:, 2:end], apply(vcat, factors)[end-size(w,1)+1:end, 1:number_of_factors]))
    # TODO: unfinished business: factor lags to get dynamic factor models (how many? see page
end

#p_values(fm::FactorModel) = map(t_stat -> 2*(1-cdf(Distributions.TDist(size(fm.x, 1)-length(fm.coefficients)), abs(t_stat))), fm.t_stats)
#make_stars(p_value::Float64) = p_value<=0.001 ? "***" : (p_value<=0.01 ? "**" : (p_value<=0.05 ? "*" : (p_value<= 0.1 ? "." : " ")))
#add_stars(p_values::Array{Float64, 1}) = [string(p_values[i], map(make_stars, p_values)[i]) for i in 1:length(p_values)]

function Base.show(io::IO, dfm::DynamicFactorModel)
    @printf io "Dynamic Factor Model\n"
    @printf io "Dimensions of X: %s\n" size(dfm.factor_model.x)
    @printf io "Number of static factors used: %s\n" dfm.factor_model.number_of_factors
    @printf io "Number of lags used (dynamic factors): %s\n" dfm.number_of_factor_lags
    @printf io "Factors calculated by: %s\n" dfm.factor_model.factor_type
    @printf io "Factor model residual variance: %s\n" sum(dfm.factor_model.factor_residuals.^2)/apply(*, size(dfm.factor_model.x))
    # TODO: show coefficient values rather than p_values
    # TODO: visual separation between coefficients for factors and lags
end

function predict(dfm::DynamicFactorModel, y::Array{Float64, 1}, h::Int64=1, number_of_lags::Int64=5, number_of_factors::Int64=0)
    if number_of_factors == 0  # number of factors can be given to set them to a different value in the forecasting equation than in the factor equation
        #println("using the same number of factors for forecasting as in the factor equation: ", dfm.factor_model.number_of_factors)
        number_of_factors = dfm.factor_model.number_of_factors
    end
    # makes a h step ahead forecast of y using a linear regression on lags of y and static factors
    start_index_w = number_of_lags > dfm.number_of_factor_lags ? 1 : (dfm.number_of_factor_lags-number_of_lags+1)  # we might have to cut off the first few obs due to lagging
    start_index_factor_lags = number_of_lags > dfm.number_of_factor_lags ? (number_of_lags-dfm.number_of_factor_lags+1) : 1  # and we have to check if more due to factors of due to lags of y
    w = lagged_matrix(y, [0, [h:number_of_lags+h-1]])[start_index_w:end, :]

    # matrix of factors and their lags
    factor_lags = apply(hcat, [lagged_matrix(dfm.factor_model.factors[:, i], Int64[0:dfm.number_of_factor_lags]) for i in 1:number_of_factors])[start_index_factor_lags:end, :]

    y, x = w[2:end, 1], hcat(w[:, 2:end], factor_lags)
    x = hcat(ones(size(x, 1)), x)  # add a constant term to the regression
    new_x, design_matrix = x[end,:], x[1:end-1, :]  # we reserve the last observation for prediction and dont use it for learning

    #println(size(design_matrix, 1), " observations vs ", size(design_matrix, 2), " variables")
    if size(design_matrix, 1) < size(design_matrix, 2)
        warn("more columns than rows in regression. Maybe try to reduce the number of common factors in the factor model.")
    end
    coefficients = inv(design_matrix'design_matrix)*design_matrix'y
    prediction = new_x*coefficients
    residuals = y - design_matrix*coefficients  # returned for e.g. info criteria
    #residuals = y - prediction
    #hat_matrix = design_matrix*inv(design_matrix'design_matrix)*design_matrix'
    #residual_variance = (residuals.^2)./(1.-diag(hat_matrix))  # HC 2
    #coefficient_covariance = inv(design_matrix'design_matrix)*(design_matrix'diagm(residual_variance)design_matrix)*inv(design_matrix'design_matrix)

    return(residuals, prediction[1])
end

include("criteria.jl")  # defines the criteria in Bai and Ng 2002


text:
@testset "NormRel" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  f = x^8 - x^4 + 1
  K, a = NumberField(f, "a", cached = false)
  S = prime_ideals_up_to(maximal_order(K), 1000)
  class_group(maximal_order(K))
  C, mC = Hecke.sunit_group_fac_elem(S)
  Q, mQ = quo(C, 3)
  CC, mCC = Hecke.NormRel._sunit_group_fac_elem_quo_via_brauer(K, S, 3)
  elts = GrpAbFinGenElem[]
  for i in 1:ngens(CC)
    u = mCC(CC[i])
    push!(elts, mQ(mC\u))
  end
  V = sub(Q, elts)[1]
  @test order(V) == order(Q)
  S = prime_ideals_up_to(maximal_order(K), Hecke.factor_base_bound_grh(maximal_order(K)))
  c, U = Hecke.NormRel._sunit_group_fac_elem_quo_via_brauer(K, S, 2)

  # Test a non-normal group

  Qx, x = PolynomialRing(FlintQQ, "x");
  K, a = NumberField(x^4-2*x^2+9);
  OK = maximal_order(K);
  lP = NfOrdIdl[]
  push!(lP, ideal(OK, 43, OK(a^2 - 12)));
  push!(lP, ideal(OK, 47, OK(a^2 - 14*a + 3)));
  push!(lP, ideal(OK, 53, OK(a^2 - 7*a - 3)));
  push!(lP, ideal(OK, 5, OK(a^2 - 1*a + 2)));
  push!(lP, ideal(OK, 2, OK(1//12*a^3 + 1//4*a^2 + 7//12*a + 7//4)));
  push!(lP, ideal(OK, 3, OK(a^2 + 1)));
  push!(lP, ideal(OK, 3, OK(8*a^3 + 4*a^2 + 2*a + 6)));
  push!(lP, ideal(OK, 37, OK(a^2 + 12*a - 3)));
  push!(lP, ideal(OK, 41, OK(a-15)));
  push!(lP, ideal(OK, 5, OK(a^2 + 1*a + 2)))
  U, mU = Hecke.NormRel._sunit_group_fac_elem_quo_via_brauer(K, lP, 8)
  S, mS = Hecke.sunit_group_fac_elem(lP)
  Q, mQ = quo(S, 8)
  V = quo(Q, [mQ(mS\(mU(U[i]))) for i in 1:ngens(U)])
  @test order(V[1]) == 1

  # Test a non-normal number field (with a C2 x C2 subgroup of automorphisms)
  # We do not test the non-quotient S-unit group, since the saturation is
  # killing us (because of assertions to full level)
  f = x^8 - 8*x^6 + 80*x^4 + 512*x^2 + 1024
  K, a = NumberField(f)
  OK = lll(maximal_order(K))
  # S invariant
  lP = prime_ideals_up_to(OK, 50)
  U, mU = Hecke.NormRel._sunit_group_fac_elem_via_brauer(K, lP)
  S, mS = Hecke.sunit_group_fac_elem(lP)
  V = quo(S, [(mS\(mU(U[i]))) for i in 1:ngens(U)])
  @test order(V[1]) == 1
end


text:
include("TestFunctions.jl")

# vim:filetype=julia:foldmethod=syntax

max_stars_repo_name: devendra-ghate/surrogateModels

max_issues_repo_name: devendra-ghate/surrogateModels

max_forks_repo_name: devendra-ghate/surrogateModels


text:
module Microclimate

using Unitful, NCDatasets
using Unitful: W, m, °C, hr, mol, K, s, J, Mg, g, kg, L, kPa, Pa

export AbstractEnvironment,
       AbstractMicroclimEnvironment,
       AbstractMicroclimGrid, MicroclimGrid,
       AbstractMicroclimPoint, MicroclimPoint,
       AbstractMicroclimInstant, MicroclimInstant,
       LinearLayerInterpolator, LinearLayerInterpolators

export radiation, snowdepth, airtemperature, relhumidity, windspeed,
       soiltemperature, soilwaterpotential, soilwatercontent

export mean_soiltemperature, mean_soilwaterpotential, mean_soilwatercontent

export weightedmean, layermax

export load_grid, load_point

radiation_type = typeof(0.1W*m^-2)
snowdepth_type = typeof(0.1m)
airtemperature_type = typeof(0.1K)
relhumidity_type = typeof(0.001)
windspeed_type = typeof(0.1m*s^-1)
soiltemperature_type = typeof(0.1K)
soilwaterpotential_type = typeof(0.0kPa)
soilwatercontent_type = typeof(0.001)

abstract type AbstractEnvironment end
abstract type AbstractMicroclimEnvironment{S,I,A,R} <: AbstractEnvironment end
abstract type AbstractMicroclimGrid{S,I,A,R} <: AbstractMicroclimEnvironment{S,I,A,R} end
abstract type AbstractMicroclimPoint{S,I,A,R} <: AbstractMicroclimEnvironment{S,I,A,R} end

struct MicroclimGrid{S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC} <: AbstractMicroclimGrid{S,I,A,R}
    radiation::RA
    snowdepth::SN
    airtemperature::AT
    relhumidity::RH
    windspeed::WS
    soiltemperature::ST
    soilwaterpotential::WP
    soilwatercontent::WC

    MicroclimGrid{S,I,A,R}(radiation::RA, snowdepth::SN, airtemperature::AT, relhumidity::RH,
                            windspeed::WS, soiltemperature::ST, soilwaterpotential::WP,
                            soilwatercontent::WC) where {S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC} =
        new{S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC}(radiation, snowdepth, airtemperature,
            relhumidity, windspeed, soiltemperature, soilwaterpotential, soilwatercontent)

    MicroclimGrid{S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC}(radiation::RA, snowdepth::SN, airtemperature::AT, relhumidity::RH,
                            windspeed::WS, soiltemperature::ST, soilwaterpotential::WP,
                            soilwatercontent::WC) where {S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC} =
        new{S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC}(radiation, snowdepth, airtemperature,
            relhumidity, windspeed, soiltemperature, soilwaterpotential, soilwatercontent)
end

struct MicroclimPoint{S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC} <: AbstractMicroclimPoint{S,I,A,R}
    radiation::RA
    snowdepth::SN
    airtemperature::AT
    relhumidity::RH
    windspeed::WS
    soiltemperature::ST
    soilwaterpotential::WP
    soilwatercontent::WC

    MicroclimPoint{S,I,A,R}(radiation::RA, snowdepth::SN, airtemperature::AT, relhumidity::RH,
                            windspeed::WS, soiltemperature::ST, soilwaterpotential::WP,
                            soilwatercontent::WC) where {S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC} =
        new{S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC}(radiation, snowdepth, airtemperature,
            relhumidity, windspeed, soiltemperature, soilwaterpotential, soilwatercontent)

    MicroclimPoint{S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC}(radiation::RA, snowdepth::SN, airtemperature::AT, relhumidity::RH,
                            windspeed::WS, soiltemperature::ST, soilwaterpotential::WP,
                            soilwatercontent::WC) where {S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC} =
        new{S,I,A,R,RA,SN,AT,RH,WS,ST,WP,WC}(radiation, snowdepth, airtemperature,
            relhumidity, windspeed, soiltemperature, soilwaterpotential, soilwatercontent)
end

abstract type AbstractMicroclimInstant end

struct MicroclimInstant{M,I,T} <: AbstractMicroclimInstant
    microclimate::M
    interp::I
    t::T
end
MicroclimInstant(microclim, height::Number, t) = MicroclimInstant(microclim, LinearLayerInterpolators(microclim, height), t)

# From external data without units
MicroclimPoint{S,I,A,R}(radiation::Vector{<:AbstractFloat}, snowdepth::Vector{<:AbstractFloat},
                      airtemperature::Array{<:AbstractFloat}, relhumidity::Array{<:AbstractFloat},
                      windspeed::Array{<:AbstractFloat}, soiltemperature::Array{<:AbstractFloat},
                      soilwaterpotential::Array{<:AbstractFloat}, soilwatercontent::Array{<:AbstractFloat}
                     ) where {S,I,A,R} = begin

    rad = to_radiation.(radiation)
    snow = to_snowdepth.(snowdepth)
    airt = to_airtemperature.(airtemperature)
    relh = to_relhumidity.(relhumidity)
    wind = to_windspeed.(windspeed)
    soilt = to_soiltemperature.(soiltemperature)
    soilwp = to_soilwaterpotential.(soilwaterpotential)
    soilwc = to_soilwatercontent.(soilwatercontent)

    args = (rad, snow, airt, relh, wind, soilt, soilwp, soilwc)
    MicroclimPoint{I,A,R,typeof.(args)...}(args...)
end

# From microclim dataset in tuples and a cartesian index
# Requires conversion to float from Int
MicroclimPoint{S,I,A,R}(radiation::Vector, snowdepth::Vector, airtemperature::Tuple,
                        relhumidity::Tuple, windspeed::Tuple, soiltemperature::Tuple,
                        soilwaterpotential::Tuple, soilwatercontent::Tuple, i::CartesianIndex) where {S,I,A,R} = begin
    rad = to_radiation.(catpoint(radiation, i) .* 0.1)
    snow = to_snowdepth.(catpoint(snowdepth, i) .* 0.001)
    airt = to_airtemperature.(catpoint(airtemperature, i) .* 0.1)
    relh = to_relhumidity.(catpoint(relhumidity, i) .* 0.001)
    wind = to_windspeed.(catpoint(windspeed, i) .* 0.1)
    soilt = to_soiltemperature.(catpoint(soiltemperature, i) .* 0.1)
    soilwp = to_soilwaterpotential.(catpoint(soilwaterpotential, i) .* 0.1)
    soilwc = to_soilwatercontent.(catpoint(soilwatercontent, i) .* 0.001)

    args = (rad, snow, airt, relh, wind, soilt, soilwp, soilwc)
    MicroclimPoint{S,I,A,R,typeof.(args)...}(args...)
end

# From MiriclimGrid and a cartesian index
MicroclimPoint(env::AbstractMicroclimGrid{S,I,A,R}, i::CartesianIndex) where {S,I,A,R} =
    MicroclimPoint{S,I,A,R}(radiation(env), snowdepth(env), airtemperature(env),
                   relhumidity(env), windspeed(env), soiltemperature(env),
                   soilwaterpotential(env), soilwatercontent(env), i)

struct MicroclimateLayers{L<:Tuple}
    layers::L
end

include("interpolation.jl")
include("netcdf.jl")

# unit conversions
to_radiation(x) = x * W*m^-2
to_snowdepth(x) = x * m
to_airtemperature(x) = x * °C |> K
to_relhumidity(x) = x
to_windspeed(x) = x * m*s^-1
to_soiltemperature(x) = x * °C |> K
to_soilwaterpotential(x) = x * kPa
to_soilwatercontent(x) = x

# Get layer specifications from type system
get_shade(layers::AbstractMicroclimEnvironment{S,I,A,R}) where {S,I,A,R} = S
get_increments(layers::AbstractMicroclimEnvironment{S,I,A,R}) where {S,I,A,R} = I
get_nextlayer(layers::AbstractMicroclimEnvironment{S,I,A,R}) where {S,I,A,R} = A
get_range(layers::AbstractMicroclimEnvironment{S,I,A,R}) where {S,I,A,R} = R

@inline catpoint(layers::Tuple{}, i::CartesianIndex) = Float64[]
@inline catpoint(layers::Tuple, i::CartesianIndex) =
    cat([[cat((getindex(y, i, :) for y in l)..., dims=1)...] for l in layers]..., dims=2)
@inline catpoint(layers::Vector, i::CartesianIndex) = cat((getindex(y, i, :) for y in layers)..., dims=1)

# Accessors for single layer data
for name in (:radiation, :snowdepth)
    eval(quote
        $name(env::AbstractEnvironment) = env.$name
        $name(env::AbstractMicroclimPoint, interp, i) = lin_interp($name(env), i)
        $name(env::AbstractMicroclimPoint, i) = lin_interp($name(env.microclimate), i)
        $name(env::AbstractMicroclimInstant) = lin_interp($name(env.microclimate), env.t)
    end)
end

# Accessor functions for range and increment data
for (name, mode) in ((:airtemperature, :range), (:relhumidity, :range), (:windspeed, :range),
                     (:soiltemperature, :increment), (:soilwaterpotential, :increment),
                     (:soilwatercontent, :increment))
    eval(quote
        $name(env::AbstractMicroclimEnvironment) = env.$name
        $name(env::AbstractMicroclimPoint, interp::LinearLayerInterpolators, i) =
            $name(env::AbstractMicroclimPoint, interp.$mode, i)
        $name(env::AbstractMicroclimPoint, interp::LinearLayerInterpolator, i) =
            interp_layer($name(env), interp, i)
        $name(env::AbstractMicroclimPoint, height::Number, i) =
            $name(env, $(Symbol(mode, :_interpolator))(height), i)
        $name(env::AbstractMicroclimInstant) =
            $name(env.microclimate, env.interp, env.t)
        $name(env::AbstractMicroclimInstant, i::Number) =
            $name(env.microclimate, env.interp, i)
        $name(env::AbstractMicroclimInstant, interp::AbstractLayerInterpolator) =
            $name(env.microclimate, interp, env.t)
    end)
end

# Accessors for weighted mean of 8 layer data
for name in (:soiltemperature, :soilwatercontent, :soilwaterpotential)
    eval(quote
        $(Symbol(:mean_, name))(env::AbstractMicroclimPoint, x, i) =
            weightedmean(env, $name(env), x, i)
        $(Symbol(:mean_, name))(env::AbstractMicroclimInstant) =
            weightedmean(env, $name(env.microclimate), env.interp, env.t)
        $(Symbol(:mean_, name))(env::AbstractMicroclimInstant, i::Number) =
            weightedmean(env, $name(env.microclimate), env.interp, i)
        $(Symbol(:mean_, name))(env::AbstractMicroclimInstant, interp::AbstractLayerInterpolator) =
            weightedmean(env, $name(env.microclimate), interp, env.t)
    end)
end

end # module


text:
using MinFEM

function solve_poisson(file_name::String)

  mesh = import_mesh(file_name)

  L = asmLaplacian(mesh)
  M = asmMassMatrix(mesh)

  n=3
  m=2
  f(x) = sin(n*x[1]*pi)*sin(m*x[2]*pi)
  s = evaluateMeshFunction(mesh, f)

  boundary = union(mesh.Boundaries[1001].Nodes,
                   mesh.Boundaries[1002].Nodes,
                   mesh.Boundaries[1003].Nodes,
                   mesh.Boundaries[1004].Nodes)

  pde = PDESystem(A=L, b=M*s, bc=zeros(mesh.nnodes), DI=boundary)

  solve(pde)

  v = ((n*pi)^2 + (m*pi)^2)*pde.state - s
  return sqrt(v'*M*v)

end

@test solve_poisson("test_v2.msh") == solve_poisson("test_v4.msh")
@test solve_poisson("test_v2.msh") < 4e-3


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

import Test
Test.@testset "Package LegendHDF5IO" begin

# test code goes here

end # testset