text:
# --- Initialize variables --- #
#residual parameters. This structure holds the priors hyperparamters for
#the gamma prior on the model residual variances. It also holds the current
#estimate of the residual precision

mutable struct Residuals
    as::Float64
    bs::Float64
    Y::Array{Float64,2}
    p::Int64
    ps::Array{Float64,1}
end

#Factors. This struct holds all information about the latent factors
#including the current number of factors; priors hyperparameters for the
#factor Loadings; as well as current values of the factor loadings; their
#precision; the factor scores; & the genetic heritability of each factor()

mutable struct LatentFactors
    r::Int64
    n::Int64
    p::Int64
    k::Int64
    df::Int64
    ad1::Float64
    bd1::Float64
    ad2::Float64
    bd2::Float64
    h2_divisions::Int64        #discretizations of heritability
    sp::Int64
    nrun::Int64

    psijh::Array{Float64,2}    #individual loadings precisions
    delta::Array{Float64,1}    # components of tauh
    tauh::Array{Float64,1}     #extra shrinkage of each loading column
    Plam::Array{Float64,2}     #total precision of each loading
    Lambda::Array{Float64,2}   #factor loadings
    h2::Array{Float64,1}       #factor heritability

    num::Int64
    no_f::Array{Float64,1}
    nofout::Array{Float64,1}

    scores::Array{Float64, 2}

    function LatentFactors(r, n, p, k, df, ad1, bd1, ad2, bd2, h2_divisions, sp, nrun)
        psijh = rand(Gamma(df/2,2/df), p, k);
        delta = [rand(Gamma(ad1+10,1/bd1), 1); rand(Gamma(ad2,1/bd2), k-1)];
        tauh = cumprod(delta);
        Plam = broadcast(*, psijh, tauh')
        Lambda = randn(p,k) .* sqrt.(1. ./ Plam)
        h2 = rand(k)
        num = 0
        no_f=zeros(nrun)
        nofout = k*ones(nrun+1)
        scores = zeros(k, n)
        new(r, n, p, k, df, ad1, bd1, ad2, bd2, h2_divisions, sp, nrun, psijh, delta, tauh, Plam, Lambda, h2, num, no_f, nofout, scores)
    end
end

#genetic_effects. This structure holds information about latent genetic
#effects. U is latent genetic effects on factor traits. d is genetic
#effects on residuals of the factor traits. Plus prior hyperparameters for
#genetic effect precisions

mutable struct GeneticEffects
    n::Int64
    as::Float64
    bs::Float64
    ps::Array{Float64, 1}
    U::Array{Float64, 2}
    d::Array{Float64, 2}
    function GeneticEffects(n, as, bs, r, p, k, h2)
        ps = rand(Gamma(as,1/bs),p)
        U  = broadcast(*, randn(k,r), sqrt.(h2))
        d  = broadcast(*, randn(p,r), 1. ./ sqrt.(ps))
        new(n, as, bs, ps, U, d)
    end
end

#interaction_effects. Similar to genetic_effects structure except for
#additional random effects that do not contribute to variation in the
#factor traits

mutable struct InteractionEffects
    as::Float64
    bs::Float64
    ps::Array{Float64, 1}
    mean::Array{Float64, 2}
    n::Int64
    W::Array{Float64, 2}
    W_out::Array{Float64, 2}
    function InteractionEffects(as, bs, p, r2)
        ps = rand(Gamma(as, 1/bs), p)
        mean = zeros(p,r2)
        n = r2
        W = broadcast(*, randn(p,r2), 1. ./ sqrt.(ps))
        W_out = zeros(p,r2)
        new(as, bs, ps, mean, n, W, W_out)
    end
end

#fixed_effects hold B
mutable struct FixedEffects
    b::Int64                # number of fixed effects, including intercept
    cov::Array{Float64, 2}  #inverse covariance of fixed effects
    mean::Array{Float64, 2} #mean of fixed effects
    B::Array{Float64, 2}    #current estimate of fixed effects
    function FixedEffects(b, p)
        cov = zeros(b,b)
        mean = zeros(p,b)
        B = randn(p,b)
        new(b, cov, mean, B)
    end
end

#Posterior holds Posterior samples & Posterior means

mutable struct PosteriorSample
    Lambda::Array{Array{Float64,2},1}
    no_f::Array{Float64, 1}
    ps::Array{Float64, 2}
    resid_ps::Array{Float64, 2}
    B::Array{Float64, 2}
    U::Array{Array{Float64,2},1}
    d::Array{Float64, 2}
    W::Array{Float64, 2}
    delta::Array{Array{Float64,1},1}
    G_h2::Array{Array{Float64,1},1}
    function PosteriorSample(n, p, k, sp, r, r2, b)
        Lambda = Array{Array{Float64,2},1}(undef,sp)
        no_f = zeros(sp)
        ps = zeros(p, sp)
        resid_ps = zeros(p ,sp)
        B = zeros(p, b)
        U = Array{Array{Float64,2},1}(undef,sp)
        d = zeros(p, r)
        W = zeros(p, r2)
        delta = Array{Array{Float64,1},1}(undef,sp)
        G_h2 = Array{Array{Float64,1},1}(undef,sp)
        new(Lambda, no_f, ps, resid_ps, B, U, d, W, delta, G_h2)
    end
end

struct InputData
    n::Int64
    p::Int64
    b::Int64
    r1::Int64
    r2::Int64
    Y_full::Array{Float64, 2}
    Y::Array{Float64, 2}
    Mean_Y::Array{Float64, 1}
    Var_Y::Array{Float64, 1}
    X::Array{Float64, 2}
    A::Array{Float64, 2}
    Z_1::Array{Float64, 2}
    Z_2::Array{Float64, 2}
    function InputData(Y, X, A, Z_1)
        n = size(Y, 1)
        p = size(Y, 2)
        b = size(X, 1)
        r1 = size(Z_1, 1)
        r2 = 0
        Y_full = copy(Y);
        Mean_Y = zeros(p);
        Var_Y = zeros(p);
        for i in 1:p
            Mean_Y[i] = mean(Y[:, i]);
            Var_Y[i]  =  var(Y[:, i]);
        end
        Y = broadcast(-, Y, Mean_Y');
        Y = broadcast(*, Y, 1. ./ sqrt.(Var_Y)');
        Z_2 = zeros(0,n);
        new(n, p, b, r1, r2, Y_full, Y, Mean_Y, Var_Y, X, A, Z_1, Z_2)
    end
    function InputData(Y, X, A, Z_1, Z_2)
        n = size(Y, 1)
        p = size(Y, 2)
        b = size(X, 1)
        r1 = size(Z_1, 1)
        r2 = size(Z_2, 1)
        Y_full = copy(Y);
        Mean_Y = zeros(p);
        Var_Y = zeros(p);
        for i in 1:p
            Mean_Y[i] = mean(Y[:, i]);
            Var_Y[i]  =  var(Y[:, i]);
        end
        Y = broadcast(-, Y, Mean_Y');
        Y = broadcast(*, Y, 1. ./ sqrt.(Var_Y)');
        new(n, p, b, r1, r2, Y_full, Y, Mean_Y, Var_Y, X, A, Z_1, Z_2)
    end
end

struct Priors
    burn::Int64
    sp::Int64
    thin::Int64

    b0::Float64
    b1::Float64

    epsilon::Float64
    h2_divisions::Int64

    k_init::Int64

    as::Float64
    bs::Float64

    df::Float64

    ad1::Float64
    bd1::Float64

    ad2::Float64
    bd2::Float64

    k_min::Float64
    prop::Float64

    nrun::Int64
    function Priors(burn, sp, thin, b0, b1, epsilon, h2_divisions, k_init,
        as, bs, df, ad1, bd1, ad2, bd2, k_min, prop)
        nrun = burn+sp*thin
        new(burn, sp, thin, b0, b1, epsilon, h2_divisions, k_init,
            as, bs, df, ad1, bd1, ad2, bd2, k_min, prop, nrun)
    end
end


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

JLLWrappers.@generate_main_file_header("Chuffed")
JLLWrappers.@generate_main_file("Chuffed", UUID("77125aae-c893-5498-99e3-e30470bfa328"))
end  # module Chuffed_jll

ext: jl


text:
struct CircularVector{T}
    v::Vector{T}
    l::Int
    CircularVector{T}(l::Integer) where T = new(zeros(T, l), l)
end
function min_abs_diff(v::CircularVector{T}) where T
    val = Inf
    for i in 1:Base.length(v)
        val = min(val, abs(v[i] - v[i-1]))
    end
    return val
end
function max_abs_diff(v::CircularVector{T}) where T
    val = 0.0
    for i in 1:Base.length(v)
        val = max(val, abs(v[i] - v[i-1]))
    end
    return val
end

function Base.getindex(V::CircularVector{T}, i::Int) where T
    return V.v[mod1(i, V.l)]
end
function Base.setindex!(V::CircularVector{T}, val::T, i::Int) where T
    V.v[mod1(i, V.l)] = val
end
function Base.length(V::CircularVector{T}) where T
    return V.l
end

mutable struct AffineSets
    n::Int  # Size of primal variables
    p::Int  # Number of linear equalities
    m::Int  # Number of linear inequalities
    extra::Int  # Number of adition linear equalities (for disjoint cones)
    A::SparseArrays.SparseMatrixCSC{Float64,Int64}
    G::SparseArrays.SparseMatrixCSC{Float64,Int64}
    b::Vector{Float64}
    h::Vector{Float64}
    c::Vector{Float64}
end

mutable struct SDPSet
    vec_i::Vector{Int}
    mat_i::Vector{Int}
    tri_len::Int
    sq_side::Int
end

mutable struct SOCSet
    idx::Vector{Int}
    len::Int
end

mutable struct ConicSets
    sdpcone::Vector{SDPSet}
    socone::Vector{SOCSet}
end

Base.@kwdef mutable struct Result
    status::Int = 0
    status_string::String = "Problem not solved"
    primal::Vector{Float64} = Float64[]
    dual_cone::Vector{Float64} = Float64[]
    dual_eq::Vector{Float64} = Float64[]
    dual_in::Vector{Float64} = Float64[]
    slack_eq::Vector{Float64} = Float64[]
    slack_in::Vector{Float64} = Float64[]
    primal_residual::Float64 = NaN
    dual_residual::Float64 = NaN
    objval::Float64 = NaN
    dual_objval::Float64 = NaN
    gap::Float64 = NaN
    time::Float64 = NaN
    iter::Int = -1
    final_rank::Int = -1
    primal_feasible_user_tol::Bool = false
    dual_feasible_user_tol::Bool = false
    certificate_found::Bool = false
    result_count::Int = 0
end

mutable struct PrimalDual
    x::Vector{Float64}
    x_old::Vector{Float64}
    y::Vector{Float64}
    y_old::Vector{Float64}

    PrimalDual(aff::AffineSets) = new(
        zeros(aff.n), zeros(aff.n), zeros(aff.m+aff.p), zeros(aff.m+aff.p)
    )
end

mutable struct WarmStart
    x::Vector{Float64}
    y_eq::Vector{Float64}
    y_in::Vector{Float64}
end

mutable struct Residuals
    dual_gap::CircularVector{Float64}
    prim_obj::CircularVector{Float64}
    dual_obj::CircularVector{Float64}
    equa_feasibility::Float64
    ineq_feasibility::Float64
    feasibility::CircularVector{Float64}
    primal_residual::CircularVector{Float64}
    dual_residual::CircularVector{Float64}
    comb_residual::CircularVector{Float64}

    Residuals(window::Int) = new(
        CircularVector{Float64}(2*window),
        CircularVector{Float64}(2*window),
        CircularVector{Float64}(2*window),
        .0,
        .0,
        CircularVector{Float64}(2*window),
        CircularVector{Float64}(2*window),
        CircularVector{Float64}(2*window),
        CircularVector{Float64}(2*window)
    )
end

const ViewVector = SubArray#{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int}}, true}
const ViewScalar = SubArray#{Float64, 1, Vector{Float64}, Tuple{Int}, true}

mutable struct AuxiliaryData
    m::Vector{LinearAlgebra.Symmetric{Float64,Matrix{Float64}}}

    Mty::Vector{Float64}
    Mty_old::Vector{Float64}

    Mx::Vector{Float64}
    Mx_old::Vector{Float64}

    y_half::Vector{Float64}
    y_temp::Vector{Float64}

    soc_v::Vector{ViewVector}
    soc_s::Vector{ViewScalar}

    function AuxiliaryData(aff::AffineSets, cones::ConicSets)
        new(
            [LinearAlgebra.Symmetric(zeros(sdp.sq_side, sdp.sq_side), :U) for sdp in cones.sdpcone],
            zeros(aff.n), zeros(aff.n),
            zeros(aff.p+aff.m), zeros(aff.p+aff.m),
            zeros(aff.p+aff.m), zeros(aff.p+aff.m),
            ViewVector[], ViewScalar[]
        )
    end
end

mutable struct Matrices
    M::SparseArrays.SparseMatrixCSC{Float64,Int64}
    Mt::SparseArrays.SparseMatrixCSC{Float64,Int64}
    c::Vector{Float64}
end

mutable struct Params
    current_rank::Vector{Int}
    target_rank::Vector{Int}
    rank_update::Int
    update_cont::Int
    min_eig::Vector{Float64}
    iter::Int
    stop_reason::Int
    stop_reason_string::String
    iteration::Int
    primal_step::Float64
    primal_step_old::Float64
    dual_step::Float64
    theta::Float64
    beta::Float64
    adapt_level::Float64
    window::Int
    time0::Float64
    norm_c::Float64
    norm_b::Float64
    norm_h::Float64
    sqrt2::Float64

    dual_feasibility::Float64
    dual_feasibility_check::Bool

    certificate_search::Bool
    certificate_search_min_iter::Int
    certificate_found::Bool

    # solution backup

    Params() = new()
end


text:
function initial_state_szegedy(sqrtstochastic::SparseMatrixCSC{<:Number})
   vec(sqrtstochastic)/sqrt(size(sqrtstochastic, 1))
end

function initial_state(qws::QWSearch{<:AbstractSzegedy})
   initial_state_szegedy(qws.model.sqrtstochastic)
end

function evolve(qwd::QWDynamics{<:AbstractSzegedy},
                state::AbstractVector{T}) where T<:Number
   _evolve_szegedy((qwd.parameters[:operators])..., state)
end

function _evolve_szegedy(operator1::AbstractMatrix{T},
                         operator2::AbstractMatrix{T},
                         state::AbstractVector{T}) where T<:Number
   operator2*(operator1*state)
end

function measure(::QWDynamics{<:AbstractSzegedy},
                 state::AbstractVector{<:Number})
   dim = floor(Int, sqrt(length(state)))
   vec(mapslices(sum, reshape((abs.(state)).^2, (dim, dim)), dims=[1]))
end

function measure(::QWDynamics{<:AbstractSzegedy},
                 state::AbstractVector{<:Number},
                 vertices::Vector{Int})
   dim = floor(Int, sqrt(length(state)))
   vec(mapslices(sum, abs.(reshape(state, (dim, dim))[:,vertices]).^2, dims=[1]))
end


text:
include("solve_master.jl")
include("solve_slave.jl")

function has_negative_reduced_cost(slave_solutions, tour_pool)
    """Check if solutions have a negative reduced cost taking acount an error
    precision."""
    if slave_solutions == []
        return(true)
    end
    for i in slave_solutions
        if tour_pool[i]["reduced_cost"]<-10^-3
            return(true)
        end
    end
    return(false)
end

function add_tours(node_tours, slave_solutions, tour_pool)
    """Add slave solutions tours with negative reduced cost to tour pool and
    the index of those solutions to the node tour index array."""
    nb_tours_added = 0
    for i in slave_solutions
        if tour_pool[i]["reduced_cost"]<-10^-3
            push!(node_tours, i)
            nb_tours_added += 1
        end
    end
    return(nb_tours_added, tour_pool)
end

function has_integer_solution(coefficients)
    """Check if all coefficents are integer taking account a precision."""
    for coef in coefficients
        if abs(round(coef)-coef)>10^-3
            return(false)
        end
    end
    return(true)
end

function update_upper_bound(global_upper_bound, coefficients, tours, tour_pool, global_best_node, node)
    """Compute the total cost of a node solution with cost of the tours taken
    into account if all coefficients are integer. Update the actual upper bound
    if needed."""
    if has_integer_solution(coefficients)
        new_upper_bound = sum([coefficients[i]*tour_pool[tour]["cost"] for (i,tour) in enumerate(tours)])
        if new_upper_bound < global_upper_bound
            return(true, new_upper_bound, node)
        end
    end
    return(false, global_upper_bound, global_best_node)
end

function not_feasible(alpha)
    return(isnan(alpha[1]))
end

function solve_relax(node, data, graph, global_upper_bound, tour_pool, branch_pool, global_best_node)
    """For a node, while the relaxation has solutions with a negative reduced
    cost, add those tour to the node tours (and the tour_pool). If solutions
    are integer, check for an upper bound update."""

    slave_solutions = []
    nb_master_solved = 0
    nb_slave_solved = 0
    nb_tours_added = 0

    # Solve node master-slave problem
    if verbose; println("Reduced cost :"); end
    while has_negative_reduced_cost(slave_solutions, tour_pool)
        (nb_tours, tour_pool) = add_tours(node["tours"], slave_solutions, tour_pool)
        nb_tours_added += nb_tours
        dual_values = solve_master(node, data, tour_pool, branch_pool, true)
        if not_feasible(node["alpha"])
            break
        end
        nb_master_solved += 1
        slave_solutions = []
        for depot in data[1, :depots]
            (slave_sol, tour_pool) = solve_slave(graph, data, tour_pool, branch_pool, depot, dual_values, node["branch"])
            push!(slave_solutions, slave_sol)
            nb_slave_solved += 1
            if verbose; overprint(string("Reduced cost : ", last(tour_pool)["reduced_cost"])); end
        end
    end

    # Update node lower bound
    node["lower_bound"] = sum([node["alpha"][i]*tour_pool[tour]["cost"] for (i,tour) in enumerate(node["tours"])])

    # Update upper bound in case of integer solution
    (bound_updated, global_upper_bound, global_best_node) = update_upper_bound(global_upper_bound, node["alpha"], node["tours"], tour_pool, global_best_node, node)

    if verbose; println("Master solved : ", nb_master_solved, " | Slave solved : ",
        nb_slave_solved, " | Tours added : ", nb_tours_added,
        " | Lower bound : ", node["lower_bound"],
        " | Upper bound updated : ", bound_updated,
        " | Number of tours : ", length(node["tours"]),
        " | Upper bound : ", global_upper_bound); end
    return(global_upper_bound, tour_pool, branch_pool, global_best_node)
end

max_stars_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_issues_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_forks_repo_name: TheoGuyard/Groomer-Optimization-Problem


text:
using CUDA.CUSPARSE, SparseArrays

m,n = 2,3
p = 0.5

for elty in [Int32, Int64, Float32, Float64]
    @testset "$typ($elty)" for typ in [CuSparseMatrixCSR, CuSparseMatrixCSC]
        x = sprand(elty, m, n, p)
        dx = typ(x)

        # zero-preserving
        y = x .* elty(1)
        dy = dx .* elty(1)
        @test dy isa typ{elty}
        @test y == SparseMatrixCSC(dy)

        # not zero-preserving
        y = x .+ elty(1)
        dy = dx .+ elty(1)
        @test dy isa CuArray{elty}
        @test y == Array(dy)

        # involving something dense
        y = x .* ones(elty, m, n)
        dy = dx .* CUDA.ones(elty, m, n)
        @test dy isa CuArray{elty}
        @test y == Array(dy)

        # multiple inputs
        y = sprand(elty, m, n, p)
        dy = typ(y)
        z = x .* y .* elty(2)
        dz = dx .* dy .* elty(2)
        @test dz isa typ{elty}
        @test z == SparseMatrixCSC(dz)
    end
end

@testset "bug: type conversions" begin
    x = CuSparseMatrixCSR(sparse([1, 2], [2, 1], [5.0, 5.0]))
    y = Int.(x)
    @test eltype(y) == Int
end


text:
function getTableauHSPARKSymmetricProjection(name, q::CoefficientsRK{T}, p::CoefficientsRK{T}, d=[]; R∞=1) where {T}

    @assert q.s == p.s

    o = min(q.o, p.o)

    a_q = q.a
    a_p = p.a

    α_q = zeros(T, q.s, 2)
    α_q[:,1] .= 0.5

    α_p = zeros(T, p.s, 2)
    α_p[:,1] .= 0.5

    a_q̃ = Array(transpose(hcat(zero(q.b), q.b)))
    a_p̃ = Array(transpose(hcat(zero(p.b), p.b)))

    α_q̃ = [[0.0  0.0]
           [0.5  R∞*0.5]]
    α_p̃ = [[0.0  0.0]
           [0.5  R∞*0.5]]

    b_q = q.b
    b_p = p.b
    β_q = [0.5, R∞*0.5]
    β_p = [0.5, R∞*0.5]

    c_q = q.c
    c_p = p.c
    c_λ = [ 0.0, 1.0]
    d_λ = [ 0.5, 0.5]

    ω_λ  = zeros(T, 1, 3)
    ω_λ .= [0.5 R∞*0.5 0.0]

    δ_λ  = zeros(T, 1, 2)
    δ_λ .= [-1.0 +1.0]

    if length(d) == 0
        return TableauHSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ)
    else
        @assert length(d) == q.s == p.s

        return TableauHSPARKprimary(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            ω_λ, δ_λ, d)
    end

end

"Tableau for Gauss-Lobatto IIIA-IIIB method with two stages and symmetric projection."
function getTableauHSPARKLobIIIAIIIB2pSymmetric()
    d = [+1.0, -1.0]
    getTableauHSPARKSymmetricProjection(:HSPARKLobIIIAIIIB2, getCoefficientsLobIIIA2(), getCoefficientsLobIIIB2(), d; R∞=-1)
end

"Tableau for Gauss-Lobatto IIIA-IIIB method with three stages and symmetric projection."
function getTableauHSPARKLobIIIAIIIB3pSymmetric()
    d = [+0.5, -1.0, +0.5]
    getTableauHSPARKSymmetricProjection(:HSPARKLobIIIAIIIB3, getCoefficientsLobIIIA3(), getCoefficientsLobIIIB3(), d; R∞=+1)
end

"Tableau for Gauss-Legendre method with s stages and symplectic projection."
function getTableauHSPARKGLRKpSymmetric(s)
    glrk = getCoefficientsGLRK(s)
    getTableauHSPARKSymmetricProjection(Symbol("HSPARKGLRK", s), glrk, glrk; R∞=(-1)^s)
end


text:
#!/usr/bin/julia

function count_trees(maze::AbstractArray{Bool, 2}, dx::Integer, dy::Integer)
    x_lim, y_lim = size(maze)
    x = dx
    y = dy
    num_trees = 0
    while x < x_lim
        if maze[begin + x, begin + y]
            num_trees += 1
        end
        x += dx
        y += dy
        y %= y_lim
    end
    return num_trees
end

function main()
    lines = open("input/03.txt", "r") do io
        map(strip, readlines(io))
    end
    maze = hcat(Array(l) for l in lines) .== '#'
    maze = transpose(maze)

    check_list = [
        [1, 1],
        [1, 3],
        [1, 5],
        [1, 7],
        [2, 1],
    ]

    @time count_trees(maze, 1, 3)
    @timev count_trees(maze, 1, 3)
    @info "Found answer" count_trees(maze, 1, 3)
    @timev prod(count_trees(maze, l...) for l in check_list)
    @info "Found answer" prod(count_trees(maze, l...) for l in check_list)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
@testset "Subpopulation report" begin

fName = "report/subpop"

subpop1 = Subpopulation( "Subpop 1" )
subpop2 = Subpopulation( "Subpop 2" )
subpop3 = Subpopulation( "Subpop 3" )

setSubpopulationSourceNode!.( [ subpop1, subpop2, subpop3 ], [ "", "Master",
    "Branch B" ] )

addSubpopulationCondition!( subpop1, MPcondition( "tenure", >=, 60 ) )
addSubpopulationCondition!( subpop2, MPcondition( "started as", ==,
    "A junior" ) )
addSubpopulationCondition!( subpop3, MPcondition( "isCareer", ==, "yes" ) )

excelSubpopReport( mpSim, 12, subpop1, subpop2, subpop3, filename = fName )
excelSubpopAgeReport( mpSim, 12, 12, :age, subpop1, subpop2, subpop3,
    filename = fName, overwrite = false )
excelSubpopAgeReport( mpSim, 12, 12, :tenure, subpop1, subpop2, subpop3,
    filename = fName, overwrite = false )
excelSubpopAgeReport( mpSim, 12, 12, :timeInNode, subpop1, subpop2, subpop3,
    filename = fName, overwrite = false )

end  # @testset "Subpopulation report"


text:
module Triangulations

using Test
using Gridap
using Gridap.Helpers
using StaticArrays

export Triangulation
export CellRefFEs
export CartesianTriangulation
export test_triangulation
export ncells
export @law
import Gridap: CellQuadrature
import Gridap: CellPoints
import Gridap: CellBasis
import Gridap: CellGeomap
import Gridap: CellField
import Gridap: evaluate, gradient, return_size
import Gridap: symmetric_gradient
import Base: div
import Gridap: trace
import Gridap: curl
import Gridap: reindex
import Base: IndexStyle
import Base: size
import Base: getindex
import Base: +, -, *

"""
Minimal interface for a mesh used for numerical integration
I.e., cellwise nodal coordinates and notion about interpolation of these
cell coordinates within the cells.
Z is the dimension of the parametric space and D the dimension of
the physical space
"""
abstract type Triangulation{Z,D} end

function CellPoints(::Triangulation{Z,D})::CellPoints{D,Float64} where {Z,D}
 @abstractmethod
end

function CellRefFEs(
  trian::Triangulation{Z})::CellValue{LagrangianRefFE{Z,Float64}} where Z
  @abstractmethod
end

function CellBasis(trian::Triangulation{Z})::CellBasis{Z,Float64} where Z
 @abstractmethod
end

function CellGeomap(self::Triangulation)
  coords = CellPoints(self)
  basis = CellBasis(self)
  lincomb(basis,coords)
end

function ncells(self::Triangulation)
  coords = CellPoints(self)
  length(coords)
end

# Testers

function test_triangulation(trian::Triangulation{Z,D}) where {Z,D}
  basis = CellBasis(trian)
  @test isa(basis,CellBasis{Z,Float64})
  coords = CellPoints(trian)
  @test isa(coords,CellPoints{D,Float64})
  @test ncells(trian) == length(coords)
  @test ncells(trian) == length(basis)
  phi = CellGeomap(trian)
  @test isa(phi,CellGeomap{Z,D,Float64})
  jac = gradient(phi)
  @test isa(jac,CellField{Z,MultiValue{Tuple{Z,D},Float64,2,Z*D}})
  reffes = CellRefFEs(trian)
  @test isa(reffes,CellValue{LagrangianRefFE{Z,Float64}})
end

# Pretty printing

import Base: show

function show(io::IO,self::Triangulation{Z,D}) where {Z,D}
  print(io,"$(nameof(typeof(self))) object")
end

function show(io::IO,::MIME"text/plain",trian::Triangulation{Z,D}) where {Z,D}
  show(io,trian)
  print(io,":")
  print(io,"\n physdim: $D")
  print(io,"\n refdim: $Z")
  print(io,"\n ncells: $(ncells(trian))")
end

# Factories

function CellField(trian::Triangulation,fun::Function)
  phi = CellGeomap(trian)
  cf = compose(fun,phi)
  _attach_triangulation(cf,trian)
end

function CellField(
  trian::Triangulation{D,Z}, fun::Function, u::CellField{Z}, v...) where {D,Z}
  phi = CellGeomap(trian)
  cf = compose(fun,phi,u,v...)
  _attach_triangulation(cf,trian)
end

function CellBasis(
  trian::Triangulation{D,Z},
  fun::Function,
  b::CellBasis{Z},
  u...) where {D,Z}
  phi = CellGeomap(trian)
  _phi = _setup_cell_field(phi)
  _u = [_setup_cell_field(ui) for ui in u]
  compose(fun,_phi,b,_u...)
end

function CellBasis(
  trian::Triangulation{D,Z},
  fun::Function,
  b::CellField{Z},
  u...) where {D,Z}
  CellField(trian,fun,b,u...)
end

_setup_cell_field(f::CellField) = cellnewaxis(f,dim=1)

_setup_cell_field(f) = f

macro law(fundef)
  s = "The @law macro is only allowed in function definitions"
  @assert isa(fundef,Expr) s
  @assert fundef.head in (:(=), :function) s
  funname = fundef.args[1].args[1]
  nargs = length(fundef.args[1].args)-1
  if nargs == 0
    x =  Symbol[]
  else
    x = fundef.args[1].args[3:end]
  end
  fundef2 = quote
    function $(funname)($(x...))
      trian = Triangulation($(x[1]))
      CellBasis(trian,$(funname),$(x...))
    end
  end
  quote
    $(esc(fundef))
    $(esc(fundef2))
  end
end

"""
Factory function to create CellQuadrature objects in a convenient way
"""
function CellQuadrature(trian::Triangulation;degree::Int=-1,order::Int=-1)
  if order != -1
    s = "`order` key-word argument in CellQuadrature constructor has been deprecated. Use `degree` instead"
    @warn s
    _deg = order
  elseif degree != -1
    _deg = degree
  else
    error("Key-word argument `degree` not assigned in CellQuadrature")
  end
  _quadrature(CellRefFEs(trian),_deg)
end

_quadrature(reffes,order) = @notimplemented

function _quadrature(reffes::ConstantCellValue,order)
  reffe = reffes.value
  poly = polytope(reffe)
  t = poly.extrusion.array.data
  q = Quadrature(t,order=order)
  ncells = length(reffes)
  ConstantCellValue(q,ncells)
  ConstantCellQuadrature(q,ncells)
end

# Concrete implementations

struct CartesianTriangulation{D,C,B,R} <: Triangulation{D,D}
  coords::C
  basis::B
  reffes::R
end

function CartesianTriangulation(partition,domain=nothing)
  _p = _setup_partition(partition)
  _d = _setup_domain(length(_p),domain)
  D = length(_p)
  coords = CartesianCellPoints(_p,_d)
  ncells = prod(_p)
  basis, reffes = _setup_basis(D,ncells)
  C = typeof(coords)
  B = typeof(basis)
  R = typeof(reffes)
  CartesianTriangulation{D,C,B,R}(coords,basis,reffes)
end

CellBasis(trian::CartesianTriangulation) = trian.basis

CellPoints(trian::CartesianTriangulation) = trian.coords

CellRefFEs(trian::CartesianTriangulation) = trian.reffes

_setup_partition(partition::Vector{Int}) = tuple(partition...)

_setup_partition(partition::NTuple{D,Int} where D) = partition

function _setup_domain(D,::Nothing)
  d =  [ i*(-1)^j for i in ones(D) for j in 1:2 ]
  _setup_domain(D,d)
end

_setup_domain(D,domain::Vector{Float64}) = tuple(domain...)

_setup_domain(D,domain::NTuple{D,Int} where D) = domain

function _setup_basis(D,ncells)
  # code = @SVector fill(HEX_AXIS,D)
  code = tuple(fill(HEX_AXIS,D)...)
  polytope = Polytope(code)
  order = 1
  orders = fill(order,D)
  reffe = LagrangianRefFE(Float64,polytope,orders)
  basis = shfbasis(reffe)
  b = ConstantCellMap(basis, ncells)
  r = ConstantCellValue(reffe,ncells)
  (b,r)
end

struct CartesianCellPoints{D} <: IndexCellArray{Point{D,Float64},1,Vector{Point{D,Float64}},D}
  p0::Point{D,Float64}
  dp::Point{D,Float64}
  partition::NTuple{D,Int}
  v::Vector{Point{D,Float64}}
end

function CartesianCellPoints(
  partition::NTuple{D,Int}, domain::NTuple{N,Float64}) where {D,N}

  dim_to_limits = tuple([(domain[2*i-1],domain[2*i]) for i in 1:D ]...)
  p0 = zero(MVector{D,Float64})
  dp = zero(MVector{D,Float64})
  for d in 1:D
    p0[d] = dim_to_limits[d][1]
    dp[d] = (dim_to_limits[d][2] - dim_to_limits[d][1]) / partition[d]
  end
  v = zeros(Point{D,Float64},2^D)
  CartesianCellPoints(Point(p0),Point(dp),partition,v)
end

IndexStyle(::Type{CartesianCellPoints{D}}) where D = IndexCartesian()

size(x::CartesianCellPoints) = x.partition

function getindex(self::CartesianCellPoints{D}, cell::Vararg{<:Integer,D}) where D
  lsize = @SVector fill(2,D)
  cis = CartesianIndices(lsize.data)
  p = zero(MVector{D,Float64})
  p0 = self.p0
  dp = self.dp
  for (k,ci) in enumerate(cis)
    for d in 1:D
      p[d] = p0[d] + ( (cell[d]-1) + (ci[d]-1) )*dp[d]
    end
    self.v[k] = p
  end
  self.v
end

# Helpers

_attach_triangulation(cf,trian) = cf

function _attach_triangulation(cf::IndexCellField,trian)
  IndexCellFieldWithTriangulation(cf,trian)
end

"""
Type used to represent a CellField with a Triangulation as metadata.
"""
struct IndexCellFieldWithTriangulation{
  Z,C<:IndexCellField,F<:Triangulation,R} <: IndexCellValue{R,1}
  cellfield::C
  trian::F
end

function IndexCellFieldWithTriangulation(
  cellfield::IndexCellField,
  trian::Triangulation{D,Z}) where {D,Z}

  C = typeof(cellfield)
  F = typeof(trian)
  R = eltype(cellfield)
  IndexCellFieldWithTriangulation{Z,C,F,R}(cellfield,trian)
end

function evaluate(f::IndexCellFieldWithTriangulation{Z},q::CellPoints{Z}) where Z
  evaluate(f.cellfield,q)
end

for op in (:+,:-,:(gradient),:(symmetric_gradient),:(div),:(trace),:(curl))
  @eval begin
    function ($op)(a::IndexCellFieldWithTriangulation)
      g = $op(a.cellfield)
      IndexCellFieldWithTriangulation(g,a.trian)
    end
  end
end

for op in (:+, :-, :*)
  @eval begin

    function ($op)(a::IndexCellFieldWithTriangulation,b::Function)
      trian = Triangulation(a)
      cf = CellField(trian,b)
      $op(a,cf)
    end

    function ($op)(a::Function,b::IndexCellFieldWithTriangulation)
      trian = Triangulation(b)
      cf = CellField(trian,a)
      $op(cf,b)
    end

  end
end

return_size(f::IndexCellFieldWithTriangulation,s::Tuple{Int}) = return_size(f.cellfield,s)

getindex(f::IndexCellFieldWithTriangulation,i::Integer) = f.cellfield[i]

size(f::IndexCellFieldWithTriangulation) = (length(f.cellfield),)

Triangulation(f::IndexCellFieldWithTriangulation) = f.trian

function reindex(
  values::IndexCellFieldWithTriangulation, indices::CellValue{<:IndexLike})
  reindex(values.cellfield,indices)
end

function reindex(
  values::IndexCellFieldWithTriangulation, indices::IndexCellValue{<:IndexLike})
  reindex(values.cellfield,indices)
end

end # module


text:
using StochasticDiffEq, DiffEqDevTools, Test
using Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems
importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear

f_linear_iip(du,u,p,t) = @.(du = 1.01*u)
σ_linear_iip(du,u,p,t) = @.(du = 0.87*u)

linear_analytic(u0,p,t,W) = @.(u0*exp(0.63155t+0.87W))

prob_sde_linear_iip = SDEProblem(SDEFunction(f_linear_iip,σ_linear_iip,
                             analytic=linear_analytic),σ_linear_iip,[1/2],(0.0,1.0))

Random.seed!(100)
dts = 1 .//2 .^(7:-1:3) #14->7 good plot

println("SROCKC2")
@time sim = test_convergence(dts,prob_sde_linear_iip,SROCKC2(),save_everystep=false,trajectories=Int(5e6),
                        weak_timeseries_errors=false)
@show sim.𝒪est[:weak_final]
@test abs(sim.𝒪est[:weak_final]-2) < 0.35
#@test abs(sim.𝒪est[:weak_l2]-2) < 0.3
#@test abs(sim.𝒪est[:weak_l∞]-2) < 0.3

prob = prob_sde_linear

println("SROCKC2")
@time sim = test_convergence(dts,prob,SROCKC2(),save_everystep=false,trajectories=Int(1e7))
@show sim.𝒪est[:weak_final]
@test abs(sim.𝒪est[:weak_final]-2) < 0.35
#@test abs(sim.𝒪est[:weak_l2]-2) < 0.3
#@test abs(sim.𝒪est[:weak_l∞]-2) < 0.3

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq


text:
export kronfact

## Represents a list of operators with degrees of freedom reduce
immutable ReducedDiscreteOperators{BT,MT}
    bcP::Array{BT,2}  # permute columns of bcs
    bcQ::Array{BT,2}  # bc normalizer
    bcs::Array{BT,2} # normalized bcs

    ops::Vector{SparseMatrixCSC{MT,Int64}}
    opcols::Vector{SparseMatrixCSC{MT,Int64}}

    domainspace::Space
    rangespace::Space
end

numbcs(A::ReducedDiscreteOperators)=size(A.bcs,1)
Base.size(A::ReducedDiscreteOperators,k...)=size(A.ops[1],k...)

function ReducedDiscreteOperators(Bx,Ls,nx)
    Ls=promotespaces(Ls)
    B,LLs=pdetoarray(Bx,Ls,nx)
    K=length(Bx)
    opcols=SparseMatrixCSC{Float64,Int64}[sparse(L[:,1:K]) for L in LLs]
    R,Q,LLs,P=regularize_bcs(B,LLs)
    ops=SparseMatrixCSC{Float64,Int64}[sparse(cont_reduce_dofs(R,full(L))[:,K+1:end]) for L in LLs]
    ReducedDiscreteOperators(P,Q,R,ops,opcols,domainspace(Ls[1]),rangespace(Ls[1]))
end

function Base.kron(Ax::ReducedDiscreteOperators,Ay::ReducedDiscreteOperators)
   @assert length(Ax.ops)==length(Ay.ops)

    M=spzeros(Float64,size(Ax,1)*size(Ay,1),size(Ax,2)*size(Ay,2))
    for k=1:length(Ax.ops)
        M+=kron(Ax.ops[k],Ay.ops[k])
    end

    M
end

regularize_bcs(S::ReducedDiscreteOperators,Gy)=length(Gy)==0?Gy:S.bcQ*Gy

function cont_reduce_dofs!{T<:Fun}(Ax::ReducedDiscreteOperators,Ay::Vector,G::Vector{T},F)
    G=regularize_bcs(Ax,G)
    cont_reduce_dofs!(Ax.opcols,Ay,G,F)
end
function cont_reduce_dofs!{T<:Fun}(Ax::ReducedDiscreteOperators,Ay::ReducedDiscreteOperators,G::Vector{T},F)
    G=regularize_bcs(Ax,G)
    cont_reduce_dofs!(Ax.opcols,Ay.ops,coefficients(G)[numbcs(Ay)+1:end,:],F)
end

immutable PDEOperatorKron{T}
    op::Operator{T} #we need to remember the op to reduce DOFs
    opsx::ReducedDiscreteOperators
    opsy::ReducedDiscreteOperators
    indsBx::Vector{Int}
    indsBy::Vector{Int}
    kron::SparseMatrixCSC{T,Int64}
end

Base.eltype{T}(::PDEOperatorKron{T})=T

for op in (:domainspace,:rangespace)
    @eval $op(P::PDEOperatorKron,k...)=$op(P.op,k...)
end

domain(P::PDEOperatorKron)=domain(domainspace(P))

function PDEOperatorKron(A,nx::Integer,ny::Integer)
    L=A[end]
    @assert iskronsumop(L)
    ops=sumops(L)

    indsBx,Bx=findfunctionals(A,1)
    Ax=ReducedDiscreteOperators(Bx,map(op->dekron(op,1),ops),nx)
    indsBy,By=findfunctionals(A,2)
    Ay=ReducedDiscreteOperators(By,map(op->dekron(op,2),ops),ny)
    PDEOperatorKron(L,Ax,Ay,indsBx,indsBy,kron(Ax,Ay))
end

#TODO: could be confused with size(A.kron)...
Base.size(A::PDEOperatorKron,k)=k==1?size(A.opsx,1):size(A.opsy,1)
Base.size(A::PDEOperatorKron)=size(A,1),size(A,2)
bcinds(A::PDEOperatorKron,k)=k==1?A.indsBx:A.indsBy
numbcs(A::PDEOperatorKron,k)=numbcs(k==1?A.opsx:A.opsy)

kronfact(A::Vector,nx::Integer,ny::Integer)=PDEOperatorKron(A,nx,ny)
kronfact(A::Vector,n::Integer)=kronfact(A,n,n)
kronfact(A::Vector,S::BivariateDomain,n::Integer)=kronfact(A,n)

function pdesolve(K::PDEOperatorKron,G)
    fx,fy,F=pde_standardize_rhs(K,G)

    F=cont_reduce_dofs!(K.opsx,map(op->dekron(op,2),sumops(K.op)),fx,F.').'
    F=cont_reduce_dofs!(K.opsy,K.opsx,fy,F)

    X22=reshape(K.kron\vec(pad(coefficients(F,rangespace(K)),size(K)...)),size(K)...)
    Kx,Ky=numbcs(K,1),numbcs(K,2) # doesn't include boundary rows...could be confusing
    nx,ny=size(K,1)+Kx,size(K,2)+Ky

    Gx=pad(coefficients(regularize_bcs(K.opsx,fx)).',:,ny)
    Gy=pad(coefficients(regularize_bcs(K.opsy,fy)).',:,nx)
    Rx,Ry=K.opsx.bcs,K.opsy.bcs

    Px,Py=K.opsx.bcP,K.opsy.bcP

    # following is copyed from constrained_lyap
    X12=Gx[:,Ky+1:end]-Rx[:,Kx+1:end]*X22
    X21=Gy[:,Kx+1:end].'-X22*Ry[:,Ky+1:end].'
    X11 = Gx[:,1:Ky] - Rx[:,Kx+1:end]*X21
    X11a= Gy[:,1:Kx].' - X12*Ry[:,Ky+1:end].'

    tol = 1e-13
    if !isempty(X11)
        bcerr = norm(X11 - X11a)

        if bcerr>tol
           warn("Boundary conditions differ by " * string(bcerr))
        end
    end

    X = [X11 X12; X21 X22]

    ProductFun(Px*X*Py.',domainspace(K))
end

\(A::PDEOperatorKron,G)=pdesolve(A,G)


text:
### Statistics

function countmap(x::RLEVector{T1,T2}) where {T1,T2}
    tally = Dict{T1,T2}()
    for (v, w) in zip(values(x), widths(x))
        tally[v] = get(tally, v, 0) + w
    end
    tally
end

function mode(x::RLEVector)
    tally = countmap(x)
    which_max = argmax(collect(values(tally)))
    max_val = collect(keys(tally))[which_max]
    max_val
end


text:
using HTTP

import Base.@kwdef, Base.WeakRef
import JSON3.StructTypes

export Config, State, UnknownVariable, ID, Var, JusCmd, ROOT, EMPTYID, Namespace, Connection, VarCommand
export parent, cancel, arg
export PASS, FAIL, SUBSTITUTE

const MODULEDIR = pathof(parentmodule(eval))
const PKGDIR = MODULEDIR !== nothing ? dirname(dirname(MODULEDIR)) : pwd()
const PASS = :pass
const FAIL = :fail
const SUBSTITUTE = :substitute

struct UnknownVariable <:Exception
    name
end

@kwdef mutable struct FakeWS
    isopen = true
    input::Channel = Channel(10)
    output::Channel = Channel(10)
end

other(ws::FakeWS) = FakeWS(input = ws.output, output = ws.input)

Base.isopen(ws::FakeWS) = ws.isopen
Base.eof(ws::FakeWS) = !isopen(ws)

#Base.write(ws::FakeWS, data) = put!(ws.output, data)
Base.write(ws::FakeWS, data) = (@debug("FAKE WRITING: $(repr(data))"); put!(ws.output, data))

Base.flush(ws::FakeWS) = nothing

#Base.readavailable(com::FakeWS) = take!(com.input)
Base.readavailable(com::FakeWS) = (str = take!(com.input); @debug("FAKE READING: $(repr(str))"); str)

function Base.close(com::FakeWS)
    com.isopen = false
    Base.close(com.input, HTTP.WebSockets.WebSocketError(1005, "Closed"))
    Base.close(com.output, HTTP.WebSockets.WebSocketError(1005, "Closed"))
end

@kwdef mutable struct State
    servers::Dict{String, @NamedTuple{host::String,port::UInt16,pid::Int32}} = Dict()
end

StructTypes.StructType(::Type{State}) = StructTypes.Mutable()

@kwdef struct ID
    namespace::String
    number::UInt = 0
end

Base.show(io::IO, id::ID) = print(io, "$(id.namespace)/$(id.number)")

const ROOT = ID("ROOT", UInt(0))
const EMPTYID = ID("", UInt(0))

"""
    Var

A variable:
    id: the variable's unique ID, used in the protocol
    name: the variable's human-readable name, used in UIs and diagnostics
"""
@kwdef mutable struct Var
    parent::ID = EMPTYID
    id::ID = ROOT # root is just the default value and will most likely be changed
    name::Union{Symbol, Integer} = ""
    value = nothing
    metadata::Dict{Symbol, AbstractString} = Dict()
    namedchildren::Dict{Symbol, Var} = Dict()
    indexedchildren::Vector{Var} = []
    properties::Dict{Symbol, Any} = Dict() # misc properties of any type
    active = true # controls refreshing
    internal_value = nothing
    readable::Bool = true
    writeable::Bool = true
    action::Bool = false
    path::Vector{Union{Number, Symbol, Function}} = []
    json_value = nothing
    refresh_exception = nothing
    error_count = 0
    value_conversion = identity
end

@kwdef mutable struct Namespace
    name::String
    secret::String
    curid::Int = 0
end

@kwdef mutable struct Connection
    ws #::HTTP.WebSockets.WebSocket
    namespace::String
    observing::Set{ID} = Set()
    stop::Channel = Channel(1)
    apps::Dict = Dict()
    oid2data::Dict{Int, WeakRef} = Dict()
    data2oid::WeakKeyDict{Any, Int} = WeakKeyDict()
    nextOid::Int = 0
    refresh_queued::Bool = false
    pending_result::NamedTuple = (;)
    vars::Set{Var} = Set{Var}()
    requests::Set = Set()
    sequence::Int = 1
end

"""
    Config

Singleton for this program's state.
    namespace:      this program's unique namespace, assigned by the server
    nextId:         the next available id in this namespace
    vars:           all known variables
    changes:        pending changes to the state to be broadcast to observers
    templates_dir:  directory for looking up viewdef generation templates
    output_dir:     optional path of directory in which to place generated viewdefs
"""
@kwdef mutable struct Config
    namespace::Namespace = Namespace("ROOT", "", 0)
    vars::Dict{ID, Var} = Dict{ID, Var}()
    host::IPv4 = ip"0.0.0.0"
    hostname::AbstractString = "localhost"
    port::UInt16 = UInt16(8181)
    server::Bool = false
    diag::Bool = false
    proxy::Bool = false
    verbose::Bool = false
    cmd::AbstractString = ""
    args::Vector{<: AbstractString} = String[]
    client::AbstractString = ""
    nextmsg::Int = 0
    namespaces::Dict{AbstractString, Namespace} = Dict(namespace.name => namespace)
    serverfunc::Function = serve
    connections::Dict{Any, Connection} = Dict{Any, Connection}()
    changes::Dict{ID, Dict{Symbol, Any}} = Dict{ID, Dict{Symbol, Any}}()
    init_connection::Function = con-> ()
    output_dir::AbstractString = ""
    templates_dir::AbstractString = joinpath(PKGDIR, "templates")
    templates::Dict{Tuple{Symbol, Symbol}, AbstractString} =
        Dict{Tuple{Symbol, Symbol}, AbstractString}()
end

@kwdef mutable struct JusCmd{NAME}
    config::Config
    ws #::HTTP.WebSockets.WebSocket
    namespace::AbstractString
    args::Vector
    cancel::Bool = false
    JusCmd(config, ws, namespace, args::Vector) =
        new{Symbol(lowercase(args[1]))}(config, ws, namespace, args[2:end], false)
end
JusCmd(cfg::Config, con::Connection, name::AbstractString) = JusCmd(cfg, con, [name])
JusCmd(cfg::Config, con::Connection, args::Vector) = JusCmd(cfg, con.ws, con.namespace, args)

connection(cmd::JusCmd) = get(cmd.config.connections, cmd.ws, nothing)

ID(cmd::JusCmd{T}) where T = ID(cmd.namespace, UInt(0))

"""
    VarCommand{Cmd, Path}

Handlers should return a command, either the given one or a new one.
Handlers can alter var during processing.
Each command triggers a refresh.

COMMANDS

Set: set a variable
Get: retrieve a value for a variable
Observe: observe variables
"""
@kwdef mutable struct VarCommand{Cmd, Arg}
    var::Var
    config::Config
    connection::Union{Connection, Nothing}
    cancel::Bool = false
    creating::Bool = false
    arg = nothing
    data = nothing
    path = ()
end

has_path(var::Var) = var.parent != EMPTYID && !isempty(var.path)
internal_value(var::Var) = has_path(var) ? var.internal_value : var.value
parent_value(cfg::Config, var::Var) = var.parent == EMPTYID ? nothing : internal_value(cfg[var.parent])

function VarCommand(cmd::Symbol, path::Union{Tuple{}, Tuple{Vararg{Symbol}}}; args...)
    VarCommand{cmd, path}(; args...)
end
function VarCommand(cmd::VarCommand{Cmd, Arg}; args...) where {Cmd, Arg}
    #VarCommand{Cmd, Arg}(; cmd.var, cmd.config, cmd.connection, cmd.cancel, cmd.creating, cmd.arg, cmd.data, args...)
    VarCommand{Cmd, Arg}(cmd; args...)
end
function VarCommand{Cmd, Arg}(cmd::VarCommand; args...) where {Cmd, Arg}
    VarCommand{Cmd, Arg}(; cmd.var, cmd.config, cmd.connection, cmd.cancel, cmd.creating, cmd.arg, cmd.data, args...)
end
function VarCommand(jcmd::JusCmd, cmd::Symbol, path::Tuple{Vararg{Symbol}}, var; args...)
    VarCommand(cmd, path; var, jcmd.config, connection = connection(jcmd), args...)
end

cancel(cmd::VarCommand) = VarCommand(cmd; cancel = true)

arg(cmd::VarCommand, arg) = VarCommand(cmd; arg)

function parent_var(cmd::VarCommand)
    cmd.var.parent == EMPTYID && throw("No parent for variable $(cmd.var.id)")
    cmd.config[cmd.var.parent]
end

parent(cmd::VarCommand) = VarCommand(cmd; var = parent_var(cmd))

Base.show(io::IO, cmd::VarCommand{Cmd, Path}) where {Cmd, Path} = print(io, "VarCommand{$(repr(Cmd)), $(Path)}($(cmd.creating ? "creating" : "not creating"))")

Base.show(io::IO, var::Var) = print(io, "Var($(var.parent), $(var.id), $(var.name), $(var.path))")

struct NoCause <: Exception end

"""
    CmdException

Error while executing a command

- type: Symbol for the type of exception:
  - path: error using a path
  - not_writeable: variable is not writeable
  - not_readable: variable is not readable
  - refresh: error while refreshing
  - program: error in program
- cmd: the command that caused the problem
- msg: description of the problem
- cause: cause of the problem (if any)
"""
struct CmdException{T} <: Exception
    type::Symbol
    cmd::VarCommand
    msg::AbstractString
    cause::Exception
    CmdException(type, cmd, msg, cause = NoCause()) = new{type}(type, cmd, msg, cause)
end


text:
markpropertyset_spec =
[
	(:x, VegaValueRef, nothing),
	(:x2, VegaValueRef, nothing),
	(:xc, VegaValueRef, nothing),
	(:width, VegaValueRef, nothing),
	(:y, VegaValueRef, nothing),
	(:y2, VegaValueRef, nothing),
	(:yc, VegaValueRef, nothing),
	(:height, VegaValueRef, nothing),
	(:opacity, VegaValueRef, nothing),
	(:fill, VegaValueRef, nothing),
	(:fillOpacity, VegaValueRef, nothing),
	(:stroke, VegaValueRef, nothing),
	(:strokeWidth, VegaValueRef, nothing),
	(:strokeOpacity, VegaValueRef, nothing),
	(:strokeDash, VegaValueRef, nothing),
	(:strokeDashOffset, VegaValueRef, nothing),
	(:cursor, VegaValueRef, nothing),
	(:clip, VegaValueRef, nothing),
	(:size, VegaValueRef, nothing),
	(:shape, VegaValueRef, nothing),
	(:path, VegaValueRef, nothing),
	(:innerRadius, VegaValueRef, nothing),
	(:outerRadius, VegaValueRef, nothing),
	(:startAngle, VegaValueRef, nothing),
	(:endAngle, VegaValueRef, nothing),
	(:interpolate, VegaValueRef, nothing),
	(:tension, VegaValueRef, nothing),
	(:orient, VegaValueRef, nothing),
	(:url, VegaValueRef, nothing),
	(:align, VegaValueRef, nothing),
	(:baseline, VegaValueRef, nothing),
	(:text, VegaValueRef, nothing),
	(:dx, VegaValueRef, nothing),
	(:dy, VegaValueRef, nothing),
	(:radius, VegaValueRef, nothing),
	(:theta, VegaValueRef, nothing),
	(:angle, VegaValueRef, nothing),
	(:font, VegaValueRef, nothing),
	(:fontSize, VegaValueRef, nothing),
	(:fontWeight, VegaValueRef, nothing),
	(:fontStyle, VegaValueRef, nothing),
	# (:title, VegaValueRef, nothing),
	# (:labels, VegaValueRef, nothing),
	# (:legend, VegaValueRef, nothing),
	# (:symbols, VegaValueRef, nothing),
	# (:gradient, VegaValueRef, nothing)

]

primitivefactory(:VegaMarkPropertySet, markpropertyset_spec)

max_stars_repo_name: jpthor/vega


text:
using Pkg
# change pwd
cd(joinpath(dirname(@__FILE__),".."))
# activate environment
VERSION < v"1.6" ? (using GitCommand) : (using Git)
Pkg.activate(pwd())
# fetch package updates
Pkg.update()
# fetch git updates
if VERSION < v"1.6"
    git() do git
        # stash local changes
        GitCommand.run(`$git stash`)
        GitCommand.run(`$git config pull.rebase false`)
        # if on CDN
        #run(`$git config --global http.proxy http://CDNusername:CDNpassword@dmzproxy005.idcn.mil.intra:8080`)
        GitCommand.run(`$git pull`)
    end
else
    Git.run(`$(git()) stash`)
    Git.run(`$(git()) config pull.rebase false`)
    Git.run(`$(git()) pull`)
end


text:
using Pythia
using Test
using CSV
using DataFrames

# Using Examples from Forecasting: Principles and Practice (2nd ed.) by Rob J Hyndman and George Athanasopoulos

### Exponential Smoothing Tests ###

## Simple Exponential Smoothing Tests

# Import Data
oildata = CSV.read("fpp2_datasets/oil.csv", DataFrame)
y = oildata[!, 3]
y = y[32:49]

# Initialize Model and Make Prediction
mdl = SES(y, h = 5)
fitted_mdl = fit(mdl, v = 1)
yhat = predict(fitted_mdl)
y_actual =

print(yhat) # Not sure what to compare this to but I know its correct though.


text:
using Lifting
using Test

@test nothing === makeDays()
@testset "Struct" begin
    @test isdefined(Lifting, :makeDays)
    SampleExercise = Exercise(;
        name = "Swiss Bar Bench",
        modality = ["Block", "Semi-Supinated"],
        equipment = "Swiss Bar",
        size = "1 inch",
        muscles = ["Pecs", "Triceps", "Front Delts"],
        trainingMax = 100,
        roundBase = 2.5,
        roundMode = floor,
    )
    @test SampleExercise == SampleExercise[1]
    @test length(SampleExercise) == 1
    for i in SampleExercise
        @test SampleExercise == i
    end

    SampleExercise2 =
        Exercise(; name = "Squat", trainingMax = 130, roundBase = 2.5, roundMode = floor)
    @test SampleExercise.name == "Swiss Bar Bench"
    @test SampleExercise.modality == ["Block", "Semi-Supinated"]
    @test SampleExercise.equipment == "Swiss Bar"
    @test SampleExercise.size == "1 inch"
    @test SampleExercise.muscles == ["Pecs", "Triceps", "Front Delts"]
    @test SampleExercise.trainingMax == 100
    @test SampleExercise.roundBase == 2.5
    @test SampleExercise.roundMode == floor

    SampleScheme = SetScheme(;
        type = ["Long Rest", "Longer Rest", "Longest Rest"],
        sets = [1, 2, 1],
        reps = [12, 14, 10],
        intensity = [7, 8.5, 9.5],
        roundMode = [floor, floor, ceil],
        rpeMode = true,
    )
    SampleScheme2 = SetScheme(;
        type = ["Long Rest", "Longer Rest", "Longest Rest"],
        sets = [1, 2, 1],
        reps = [12, 14, 10],
        intensity = [7, 8.5, 9.5],
        roundMode = [floor, floor, ceil],
        rpeMode = true,
    )
    @test SampleScheme.type == ["Long Rest", "Longer Rest", "Longest Rest"]
    @test SampleScheme.sets == [1, 2, 1]
    @test SampleScheme.reps == [12, 14, 10]
    @test SampleScheme.intensity ==
          [0.6046188561906853, 0.6000914761292184, 0.6936559685445272]
    @test SampleScheme.roundMode == [floor, floor, ceil]
    @test SampleScheme.rpeMode == true
    @test SampleScheme.addWeight == [0, 0, 0]
    @test SampleScheme.wght == [0, 0, 0]

    SampleScheme = SetScheme(;
        type = ["Long Rest", "Longer Rest", "Longest Rest"],
        sets = [1, 2, 1],
        reps = [12, 14, 10],
        intensity = [0.7, 0.65, 0.75],
        addWeight = [2.5, 5, 0],
    )
    @test SampleScheme.intensity == [0.7, 0.65, 0.75]
    @test SampleScheme.roundMode == [floor, floor, floor]
    @test SampleScheme.rpeMode == false
    @test SampleScheme.wght == [0, 0, 0]
    for i in SampleScheme
        @test SampleScheme == i
    end

    SampleProgression = Progression(;
        type = LinearProgression(),
        name = "Progression Name",
        setScheme = SampleScheme,
    )
    @test SampleProgression == SampleProgression[1]
    @test length(SampleProgression) == 1
    for i in SampleProgression
        @test SampleProgression == i
    end

    SampleProgression2 = Progression(;
        type = LinearProgression(),
        name = "Progression Name",
        sessions = 1,
        period = 2,
        setScheme = [deepcopy(SampleScheme), SampleScheme2],
    )

    @test SampleProgression.type == LinearProgression()
    @test SampleProgression.name == "Progression Name"
    @test SampleProgression.setScheme == SampleScheme

    calcWeights(SampleExercise, SampleProgression)
    calcWeights(SampleExercise2, SampleProgression2)
    @test SampleProgression.setScheme.wght == 100 * [0.7, 0.65, 0.75] + [2.5, 5, 0]

    struct SampleProgramme <: AbstractProgramme end
    exerProg = Dict()
    push!(
        exerProg,
        "SampleExercise" => (exercise = SampleExercise, progression = SampleProgression),
    )
    push!(
        exerProg,
        "SampleExercise2" => (exercise = SampleExercise2, progression = SampleProgression2),
    )

    function makeDays(::SampleProgramme, exerProg)
        week = [[], [], [], [], [], [], []]
        SampleExercise = exerProg["SampleExercise"].exercise
        SampleProgression = exerProg["SampleExercise"].progression
        push!(week[1], SampleExercise, SampleProgression)   # Day 1
        push!(week[1], SampleExercise, SampleProgression)   # Day 1
        insert!(week[1], 1, SampleExercise2, SampleProgression2)   # Day 1
        push!(week[2], "Rest")  # Day 2
        push!(week[3], SampleExercise, SampleProgression)   # Day 3
        push!(week[4], "Rest")  # Day 4
        push!(week[5], SampleExercise, SampleProgression)   # Day 5
        push!(week[5], "Rest")  # Day 6
        push!(week[7], "Rest")  # Day 7
        return week
    end

    week = makeDays(SampleProgramme(), exerProg)
    sampleProgramme = Programme(SampleProgramme(), "SampleProgramme", exerProg, week)
    @test sampleProgramme.type == SampleProgramme()
    @test sampleProgramme.name == "SampleProgramme"
    @test sampleProgramme.exerProg == exerProg
    @test sampleProgramme.days == week
    test = []
    push!(test, SampleExercise2, SampleProgression2)
    @test week[1][1] == test[1]
    @test sampleProgramme[1] == sampleProgramme.days[1]
    @test sampleProgramme[1:2:3] == [sampleProgramme.days[1], sampleProgramme.days[3]]
end

max_stars_repo_name: dcelisgarza/Lifting


text:
using Weave

jmd  = "/mnt/d/julia_dev/FastGroupBy/README.jmd"
out_path = "/mnt/d/julia_dev/FastGroupBy/"
weave(jmd, out_path = out_path, doctype = "pandoc")


text:
using Test
using LinearAlgebra
using DFTK

include("testcases.jl")

@testset "Test compute_fft_size on Silicon" begin
    model = Model(silicon.lattice; n_electrons=1)
    @test compute_fft_size(model,  3, supersampling=2)   == (15, 15, 15)
    @test compute_fft_size(model,  4, supersampling=2)   == (15, 15, 15)
    @test compute_fft_size(model,  5, supersampling=2)   == (18, 18, 18)
    @test compute_fft_size(model, 15, supersampling=2)   == (27, 27, 27)
    @test compute_fft_size(model, 25, supersampling=2)   == (36, 36, 36)
    @test compute_fft_size(model, 30, supersampling=2)   == (40, 40, 40)
    @test compute_fft_size(model, 30, supersampling=1.8) == (36, 36, 36)
end

@testset "Test compute_fft_size on skewed lattice" begin
    lattice = Diagonal([1, 1e-12, 1e-12])
    model   = Model(lattice; n_electrons=1)
    @test compute_fft_size(model, 15, supersampling=2)  == ( 5, 1, 1)
    @test compute_fft_size(model, 300, supersampling=2) == (18, 1, 1)
end

@testset "Test compute_fft_size with :precise" begin
    model = Model(silicon.lattice, silicon.atoms, silicon.positions; terms=[Kinetic()])

    function fft_precise(model, kcoords, Ecut; supersampling=2)
        fft_size_fast = compute_fft_size(model, Ecut; supersampling, algorithm=:fast)
        fft_size      = compute_fft_size(model, Ecut, kcoords;
                                         supersampling, algorithm=:precise)
        @test all(fft_size .<= fft_size_fast)
        fft_size
    end

    @test fft_precise(model, silicon.kcoords, 20, supersampling=2)   == (30, 30, 30)
    @test fft_precise(model, silicon.kcoords, 20, supersampling=1.6) == (24, 24, 24)
end


text:
type hpipm_solver
    # "attribtues"
    py_hpipm_solver
    py_hpipm_data
    py_hpipm_dims

    # "methods"
    solve
    print_sol

    function hpipm_solver(ju_qp_dims, ju_qp_data)
        this = new()

        py_qp_data = hp.hpipm_data()
        this.py_hpipm_data = ju2py_obj(ju_qp_data, py_qp_data)

        py_qp_dims = hp.hpipm_dims()
        this.py_hpipm_dims = ju2py_obj(ju_qp_dims, py_qp_dims)

        py_qp_solver = hp.hpipm_solver(this.py_hpipm_dims, this.py_hpipm_data)
        this.py_hpipm_solver = py_qp_solver

        this.solve = function()
            (this.py_hpipm_solver)[:solve]()
        end

        this.print_sol = function()
            this.py_hpipm_solver[:print_sol]()
        end

        return this
    end
end


text:
using BinaryBuilder

# Collection of sources required to build Nettle
name = "libsnark"
version = v"2019.06.18"
sources = [
    GitSource("https://github.com/scipr-lab/libsnark.git", "477c9dfd07b280e42369f82f89c08416319e24ae"),
    GitSource("https://github.com/scipr-lab/libfqfft.git", "e0183b2cef7d4c5deb21a6eaf3fe3b586d738fe0"),
    GitSource("https://github.com/herumi/ate-pairing.git", "e69890125746cdaf25b5b51227d96678f76479fe"),
    GitSource("https://github.com/google/googletest.git", "3a4cf1a02ef4adc28fccb7eef2b573b14cd59009"),
    GitSource("https://github.com/scipr-lab/libff.git", "f2067162520f91438b44e71a2cab2362f1c3cab4"),
    #GitSource("https://github.com/herumi/xbyak.git", "f0a8f7faa27121f28186c2a7f4222a9fc66c283d"),
    GitSource("https://github.com/mbbarbosa/libsnark-supercop.git", "b04a0ea2c7d7422d74a512ce848e762196f48149"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/

# Link "submodules" in
for submodule in ate-pairing libff libfqfft libsnark-supercop; do
    rm -rf libsnark/depends/${submodule}
    ln -sf $(pwd)/${submodule} libsnark/depends/${submodule}
done
rm -rf libsnark/depends/gtest
ln -sf $(pwd)/googletest libsnark/depends/gtest

# ....recursively....
for submodule in ate-pairing libff; do
    rm -rf libfqfft/depends/${submodule}
    ln -sf $(pwd)/${submodule} libfqfft/depends/${submodule}
done
rm -rf libfqfft/depends/gtest
ln -sf $(pwd)/googletest libfqfft/depends/gtest

# We need the boost headers, so install them
apk update && apk add boost-dev

# For some reason, CMake isn't automatically including `${prefix}/include`....
export CFLAGS=-I${prefix}/include

mkdir build && cd build
cmake ../libsnark -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DMULTICORE=ON \
      -DCURVE=ALT_BN128 \
      -DWITH_PROCPS=OFF
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# Disable FreeBSD for now, because hogweed needs alloca()?
#filter!(!Sys.isfreebsd, platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libsnark", :libsnark),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll", v"6.1.2"),
    Dependency("OpenSSL_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function _transmission_pathway_density(pop::Population,
                                       tn::TNDistribution,
                                       i::I, j::I) where {I <: Integer}
  x = pop.risks[[i; j], :x]
  y = pop.risks[[i; j], :y]
  dens = tn.internal[j, i]
  @debug "Transmission pathway density for individuals $i and $j" xy1 = (x[1], y[1]) xy2 = (x[2], y[2]) density = dens
  return x, y, dens
end

@recipe function f(tn::TNDistribution,
                   pop::Population;
                   show_individuals=true::Bool) where T <: EpidemicModel
  xguide --> ""
  yguide --> ""
  legend --> :none
  xlims --> extrema(pop.risks[!, :x]) .+ (sum(extrema(pop.risks[!, :x]).*(-1,1)) .* (-0.05, 0.05))
  ylims --> extrema(pop.risks[!, :y]) .+ (sum(extrema(pop.risks[!, :y]).*(-1,1)) .* (-0.05, 0.05))
  aspect_ratio := :equal
  axis --> nothing
  framestyle --> :none
  linecolor --> :black
  n = pop.individuals
  for i = 1:n, j =1:n
    if tn.internal[j, i] > 0.0
      @series begin
        x, y, density = _transmission_pathway_density(pop, tn, i, j)
        seriestype --> :path
        linealpha --> density
        x, y
      end
    end
  end
  if show_individuals
    @series begin
      pop
    end
  end
end


text:
# functions to load the data
# from John Hopkins CSSE repository
# https://github.com/CSSEGISandData/COVID-19
#
# the data should already be downloaded from the repository and be on
# your local hard drive. Make sure to adapt the path variables below

using CSV, DataFrames, Dates

function smooth(x)
# make log-transform and smooth data
    lx = log10.(x)
    b = (1/2)*ones(2)     # flat (rectangular) window -> moving average
    lx = filtfilt(b,lx)
    10 .^ lx
end

# read-in covid-19 time series
# for infected, the funcion returns the prevalence (active cases)
# this is calculated by subtracting dead and recovered from cumulative cases
function covid_ts(country, starting_date)

  df_confirmed = CSV.read(csse_path*"time_series_covid19_confirmed_global.csv")
  df_death     = CSV.read(csse_path*"time_series_covid19_deaths_global.csv")
  df_recovered = CSV.read(csse_path*"time_series_covid19_recovered_global.csv")
  first_date = Date(2020,1,22)              # first entry in database: Jan 22, 2020
  countries = df_confirmed[:,2]
  countriesR = df_recovered[:,2]  # recovered still have a different file format

  inds = findall(x->x==country, countries)
  indsR = findall(x->x==country, countriesR)
  maxdays = size(df_confirmed)[2] - 4          # number of reported days
  days = Dates.value(starting_date - first_date) # days we are interested
  I = zeros(Int,maxdays-days)
  D = zeros(Int,maxdays-days)
  R = zeros(Int,maxdays-days)
  for i in inds
      I .+= collect(df_confirmed[i, 5 + days:end])  # sum over all Provinces/States in a country
      D .+= collect(df_death[i,     5 + days:end])
  end

  for i in indsR
      R .+= collect(df_recovered[i, 5 + days:end])
  end
  I = I .- R .- D
  1:length(I),I,R,D
end

covid_ts(country) = covid_ts(country, Date(2020,1,22))

# calculate arrival times: first records of confirmed cases > th for all countries
function arrival_times(th=0, pandemie_start=Date(2019,12,1))
  # th: threshold for arrival
  # pandemie_start: hypothetical starting date of the pandemie
  df_confirmed = CSV.read(csse_path*"time_series_covid19_confirmed_global.csv")
  df_death     = CSV.read(csse_path*"time_series_covid19_deaths_global.csv")
  df_recovered = CSV.read(csse_path*"time_series_covid19_recovered_global.csv")

  first_date = Date(2020,1,22)          # first entry in database: Jan 22, 2020
  countries = df_confirmed[:,2]
  maxdays = size(df_confirmed)[2] - 4          # number of reported days

  println("n countries: ", length(unique(countries)))

  firsts = zeros(Int,length(countries))
  for i = 1:length(countries)
      R = collect(df_confirmed[i,5:end])
      x = findfirst(x->x>th, R)
      firsts[i] = x == nothing ? maxdays : x
  end

  # collect first records over Provinces/States of a country
  # probably can be done with two clever lines of code
  # here the hard way..
  p = sortperm(countries)
  countries = countries[p]
  firsts = firsts[p]

  sep_countries = [countries[1]]
  sep_firsts    = [firsts[1]]
  c = countries[1]
  i = firsts[1]
  for j=2:length(countries)
      c1 = countries[j]
      i1 = firsts[j]
      if c1==c # same countrie
          if i1 < i
             sep_firsts[end] = i1
             i = i1
          end
      else
         push!(sep_countries,c1)
         push!(sep_firsts,i1)
         c = c1
         i = i1
      end
  end

  arrivals = first_date .+ Day.(sep_firsts .- 1)
  time_diff = Dates.value.(arrivals .- pandemie_start)  # time difference in days to start of pandemie

  # return vector of countries, dates of arrival, and time_difference
  sep_countries, arrivals, time_diff
end

function read_all_data()
  df_confirmed = CSV.read(csse_path*"time_series_covid19_confirmed_global.csv")
  df_death     = CSV.read(csse_path*"time_series_covid19_deaths_global.csv")
  #df_recovered = CSV.read(csse_path*"time_series_covid19_recovered_global.csv")
  first_date = Date(2020,1,22)              # first entry in database: Jan 22, 2020
  countries = df_confirmed[:,2]

  unique_countries = unique(countries)
  ncountries = length(unique_countries)
  maxdays = size(df_confirmed)[2] - 4          # number of reported days
  I = zeros(Int,maxdays,ncountries)
  D = zeros(Int,maxdays,ncountries)
  #R = zeros(Int,maxdays,ncountries)

  j = 0
  for c in unique_countries
      j += 1
      inds = findall(x->x==c, countries)
      for i in inds
          I[:,j] .+= collect(df_confirmed[i, 5:end])  # sum over all Provinces/States in a country
          D[:,j] .+= collect(df_death[i,     5:end])
          #R[:,j] .+= collect(df_recovered[i, 5:end])
      end
  end
  1:length(I),I,D,unique_countries
  #1:length(I),I,D,R,unique_countries

end

# this is a path to the Johns Hopkins COVID19 repository
# it is assumed that you have a clone for this on your computer
csse_path = "../../COVID-19/csse_covid_19_data/csse_covid_19_time_series/"

#end # module


text:
function setref!(x::EXPR, binding::Binding)
    if !hasmeta(x)
        x.meta = Meta()
    end
    x.meta.ref = binding
    push!(binding.refs, x)
end

function setref!(x::EXPR, binding)
    if !hasmeta(x)
        x.meta = Meta()
    end
    x.meta.ref = binding
end

# Main function to be called. Given the `state` tries to determine what `x`
# refers to. If it remains unresolved and is in a delayed evaluation scope
# (i.e. a function) it gets pushed to list (.urefs) to be resolved after we've
# run over the entire top-level scope.
function resolve_ref(x, state)
    if !(parentof(x) isa EXPR && headof(parentof(x)) === :quotenode)
        resolve_ref(x, state.scope, state)
    end
end

# The first method that is tried. Searches the current scope for local bindings
# that match `x`. Steps:
# 1. Check whether we've already checked this scope (inifinite loops are
# possible when traversing nested modules.)
# 2. Check what sort of EXPR we're dealing with, separate name from EXPR that
# binds.
# 3. Look in the scope's variable list for a binding matching the name.
# 4. If 3. is unsuccessful, check whether the scope imports any modules then check them.
# 5. If no match is found within this scope check the parent scope.
# The return value is a boolean that is false if x should point to something but
# can't be resolved.

function resolve_ref(x::EXPR, scope::Scope, state::State)::Bool
    # if the current scope is a soft scope we should check the parent scope first
    # before trying to resolve the ref locally
    # if is_soft_scope(scope) && parentof(scope) isa Scope
    #     resolve_ref(x, parentof(scope), state) && return true
    # end

    hasref(x) && return true
    resolved = false

    if is_getfield(x)
        return resolve_getfield(x, scope, state)
    elseif iskwarg(x)
        # Note to self: this seems wronge - Binding should be attached to entire Kw EXPR.
        if isidentifier(x.args[1])
            setref!(x.args[1], Binding(noname, nothing, nothing, []))
        elseif isdeclaration(x.args[1]) && isidentifier(x.args[1].args[1])
            setref!(x.args[1].args[1], Binding(noname, nothing, nothing, []))
        end
        return true
    elseif is_special_macro_term(x) || new_within_struct(x)
        setref!(x, Binding(noname, nothing, nothing, []))
        return true
    end
    mn = nameof_expr_to_resolve(x)
    mn === nothing && return true

    if scopehasbinding(scope, mn)
        setref!(x, scope.names[mn])
        resolved = true
    elseif scope.modules isa Dict && length(scope.modules) > 0
        for m in values(scope.modules)
            resolved = resolve_ref_from_module(x, m, state)
            resolved && return true
        end
    end
    if !resolved && !CSTParser.defines_module(scope.expr) && parentof(scope) isa Scope
        return resolve_ref(x, parentof(scope), state)
    end
    return resolved
end

# Searches a module store for a binding/variable that matches the reference `x1`.
function resolve_ref_from_module(x1::EXPR, m::SymbolServer.ModuleStore, state::State)::Bool
    hasref(x1) && return true

    if CSTParser.ismacroname(x1)
        x = x1
        if valof(x) == "@." && m.name == VarRef(nothing, :Base)
            # @. gets converted to @__dot__, probably during lowering.
            setref!(x, m[:Broadcast][Symbol("@__dot__")])
            return true
        end

        mn = Symbol(valof(x))
        if isexportedby(mn, m)
            setref!(x, maybe_lookup(m[mn], state))
            return true
        end
    elseif isidentifier(x1)
        x = x1
        if Symbol(valof(x)) == m.name.name
            setref!(x, m)
            return true
        elseif isexportedby(x, m)
            setref!(x, maybe_lookup(m[Symbol(valof(x))], state))
            return true
        end
    end
    return false
end

function resolve_ref_from_module(x::EXPR, scope::Scope, state::State)::Bool
    hasref(x) && return true
    resolved = false

    mn = nameof_expr_to_resolve(x)
    mn === nothing && return true

    if scope_exports(scope, mn, state)
        setref!(x, scope.names[mn])
        resolved = true
    end
    return resolved
end

"""
    scope_exports(scope::Scope, name::String)

Does the scope export a variable called `name`?
"""
function scope_exports(scope::Scope, name::String, state)
    if scopehasbinding(scope, name) && (b = scope.names[name]) isa Binding
        initial_pass_on_exports(scope.expr, name, state)
        for ref in b.refs
            if ref isa EXPR && parentof(ref) isa EXPR && headof(parentof(ref)) === :export
                return true
            end
        end
    end
    return false
end

"""
    initial_pass_on_exports(x::EXPR, server)

Export statements need to be (pseudo) evaluated each time we consider
whether a variable is made available by an import statement.
"""

function initial_pass_on_exports(x::EXPR, name, state)
    for a in x.args[3] # module block expressions
        if headof(a) === :export
            for i = 1:length(a.args)
                if isidentifier(a.args[i]) && valof(a.args[i]) == name && !hasref(a.args[i])
                    Delayed(scopeof(x), state.env, state.server)(a.args[i])
                end
            end
        end
    end
end

# Fallback method
function resolve_ref(x::EXPR, m, state::State)::Bool
    return hasref(x)::Bool
end

rhs_of_getfield(x::EXPR) = CSTParser.is_getfield_w_quotenode(x) ? x.args[2].args[1] : x
lhs_of_getfield(x::EXPR) = rhs_of_getfield(x.args[1])

"""
    resolve_getfield(x::EXPR, parent::Union{EXPR,Scope,ModuleStore,Binding}, state::State)::Bool

Given an expression of the form `parent.x` try to resolve `x`. The method
called with `parent::EXPR` resolves the reference for `parent`, other methods
then check whether the Binding/Scope/ModuleStore to which `parent` points has
a field matching `x`.
"""
function resolve_getfield(x::EXPR, scope::Scope, state::State)::Bool
    hasref(x) && return true
    resolved = resolve_ref(x.args[1], scope, state)
    if isidentifier(x.args[1])
        lhs = x.args[1]
    elseif CSTParser.is_getfield_w_quotenode(x.args[1])
        lhs = lhs_of_getfield(x)
    else
        return resolved
    end
    if resolved && (rhs = rhs_of_getfield(x)) !== nothing
        resolved = resolve_getfield(rhs, refof(lhs), state)
    end
    return resolved
end

function resolve_getfield(x::EXPR, parent_type::EXPR, state::State)::Bool
    hasref(x) && return true
    resolved = false
    if isidentifier(x)
        if CSTParser.defines_module(parent_type) && scopeof(parent_type) isa Scope
            resolved = resolve_ref(x, scopeof(parent_type), state)
        elseif CSTParser.defines_struct(parent_type)
            if scopehasbinding(scopeof(parent_type), valof(x))
                setref!(x, scopeof(parent_type).names[valof(x)])
                resolved = true
            end
        end
    end
    return resolved
end

function resolve_getfield(x::EXPR, b::Binding, state::State)::Bool
    hasref(x) && return true
    resolved = false
    if b.val isa Binding
        resolved = resolve_getfield(x, b.val, state)
    elseif b.val isa SymbolServer.ModuleStore || (b.val isa EXPR && CSTParser.defines_module(b.val))
        resolved = resolve_getfield(x, b.val, state)
    elseif b.type isa Binding
        resolved = resolve_getfield(x, b.type.val, state)
    elseif b.type isa SymbolServer.DataTypeStore
        resolved = resolve_getfield(x, b.type, state)
    end
    return resolved
end

function resolve_getfield(x::EXPR, parent_type, state::State)::Bool
    hasref(x)
end

function is_overloaded(val::SymbolServer.SymStore, scope::Scope)
    vr = val.name isa SymbolServer.FakeTypeName ? val.name.name : val.name
    haskey(scope.overloaded, vr)
end

function resolve_getfield(x::EXPR, m::SymbolServer.ModuleStore, state::State)::Bool
    hasref(x) && return true
    resolved = false
    if CSTParser.ismacroname(x) && (val = maybe_lookup(SymbolServer.maybe_getfield(Symbol(valofid(x)), m, getsymbols(state)), state)) !== nothing
        setref!(x, val)
        resolved = true
    elseif isidentifier(x) && (val = maybe_lookup(SymbolServer.maybe_getfield(Symbol(valofid(x)), m, getsymbols(state)), state)) !== nothing
        # Check whether variable is overloaded in top-level scope
        tls = retrieve_toplevel_scope(state.scope)
        # if tls.overloaded !== nothing && (vr = val.name isa SymbolServer.FakeTypeName ? val.name.name : val.name; haskey(tls.overloaded, vr))
        #     @info 1
        #     setref!(x, tls.overloaded[vr])
        #     return true
        # end
        vr = val.name isa SymbolServer.FakeTypeName ? val.name.name : val.name
        if haskey(tls.names, valof(x)) && tls.names[valof(x)] isa Binding && tls.names[valof(x)].val isa SymbolServer.FunctionStore
            setref!(x, tls.names[valof(x)])
            return true
        elseif tls.overloaded !== nothing && haskey(tls.overloaded, vr)
            setref!(x, tls.overloaded[vr])
            return true
        end
        setref!(x, val)
        resolved = true
    end
    return resolved
end

function resolve_getfield(x::EXPR, parent::SymbolServer.DataTypeStore, state::State)::Bool
    hasref(x) && return true
    resolved = false
    if isidentifier(x) && Symbol(valof(x)) in parent.fieldnames
        fi = findfirst(f -> Symbol(valof(x)) == f, parent.fieldnames)
        ft = parent.types[fi]
        val = SymbolServer._lookup(ft, getsymbols(state), true)
        # TODO: Need to handle the case where we get back a FakeUnion, etc.
        setref!(x, Binding(noname, nothing, val, []))
        resolved = true
    end
    return resolved
end

resolvable_macroname(x::EXPR) = isidentifier(x) && CSTParser.ismacroname(x) && refof(x) === nothing

nameof_expr_to_resolve(x) = isidentifier(x) ? valofid(x) : nothing

"""
    valofid(x)

Returns the string value of an expression for which `isidentifier` is true,
i.e. handles NONSTDIDENTIFIERs.
"""
valofid(x::EXPR) = headof(x) === :IDENTIFIER ? valof(x) : valof(x.args[2])

"""
new_within_struct(x::EXPR)

Checks whether x is a reference to `new` within a datatype constructor.
"""
new_within_struct(x::EXPR) = isidentifier(x) && valofid(x) == "new" && is_in_fexpr(x, CSTParser.defines_struct)
is_special_macro_term(x::EXPR) = isidentifier(x) && (valofid(x) == "__source__" || valofid(x) == "__module__") && is_in_fexpr(x, CSTParser.defines_macro)


text:
abstract type PairwisePotential end

struct PseudoHS{T<:AbstractFloat} <: PairwisePotential
    λ::T
    σ::T
    a::T
    b::T
    temp::T
end

function _pseudoHS(λ::T, σ::T, temp::T) where {T<:AbstractFloat}
    b = λ / (λ - 1.0)
    b = convert(T, b)
    a = λ * b^(λ - 1.0)
    a = convert(T, a)
    return PseudoHS{T}(λ, σ, a, b, temp)
end

PseudoHS(λ::T, σ::T, temp::T) where {T<:AbstractFloat} = _pseudoHS(λ, σ, temp)
function PseudoHS(λ::T, temp::T) where {T<:AbstractFloat}
    return _pseudoHS(λ, convert(T, 1.0), temp)
end

@inline function _potential_func(
    x,
    λ::T,
    a::T,
    temp::T,
) where {T<:AbstractFloat}
    energy = (a / temp) * (x^-λ - x^-(λ - 1.0))
    return energy + (1.0 / temp)
end

@inline function _force_func(x, λ::T, a::T, temp::T) where {T<:AbstractFloat}
    force = (λ * x^-(λ + 1.0)) - ((λ - 1.0) * x^-λ)
    return force * (a / temp)
end  # function _force_func

function _apply!(x, λ, a, b, temp)
    if x < b
        energy = _potential_func(x, λ, a, temp)
        energy = oftype(b, energy)
        force = _force_func(x, λ, a, temp)
        force = oftype(b, force)
    else
        energy = zero(b)
        force = zero(b)
    end
    return energy, force
end

function apply!(p::PseudoHS, x::AbstractFloat)
    _apply!(x, p.λ, p.a, p.b, p.temp)
end
# TODO: Work on a cleaner interface for custom potentials


text:
 pv[1] =   x[1]^2*x[2]
 pv[2] =   x[1]*x[2]*x[5]
 pv[3] =   x[1]^3*x[2]
 pv[4] =   x[1]^2*x[2]^2
 pv[5] =   x[1]^2*x[2]*x[5]
 pv[6] =   x[1]^2*x[2]*x[3]
 pv[7] =   x[1]^2*x[5]*x[3]
 pv[8] =   x[1]^4*x[2]
 pv[9] =   x[1]^3*x[2]^2
 pv[10] =   x[1]^3*x[2]*x[5]
 pv[11] =   x[1]^2*x[2]^2*x[5]
 pv[12] =   x[1]^3*x[2]*x[3]
 pv[13] =   x[1]^2*x[2]^2*x[3]
 pv[14] =   x[1]^3*x[5]*x[3]
 pv[15] =   x[1]^2*x[2]*x[5]*x[3]
 pv[16] =   x[1]^5*x[2]
 pv[17] =   x[1]^4*x[2]^2
 pv[18] =   x[1]^3*x[2]^3
 pv[19] =   x[1]^4*x[2]*x[5]
 pv[20] =   x[1]^3*x[2]^2*x[5]
 pv[21] =   x[1]^2*x[2]^2*x[5]^2
 pv[22] =   x[1]^4*x[2]*x[3]
 pv[23] =   x[1]^3*x[2]^2*x[3]
 pv[24] =   x[1]^4*x[5]*x[3]
 pv[25] =   x[1]^3*x[2]*x[5]*x[3]
 pv[26] =   x[1]^2*x[2]^2*x[5]*x[3]
 pv[27] =   x[1]^3*x[5]^2*x[3]
 pv[28] =   x[1]^7
 pv[29] =   x[1]^6*x[2]
 pv[30] =   x[1]^5*x[2]^2
 pv[31] =   x[1]^4*x[2]^3
 pv[32] =   x[1]^5*x[2]*x[5]
 pv[33] =   x[1]^4*x[2]^2*x[5]
 pv[34] =   x[1]^3*x[2]^3*x[5]
 pv[35] =   x[1]^3*x[2]^2*x[5]^2
 pv[36] =   x[1]^5*x[2]*x[3]
 pv[37] =   x[1]^4*x[2]^2*x[3]
 pv[38] =   x[1]^3*x[2]^3*x[3]
 pv[39] =   x[1]^5*x[5]*x[3]
 pv[40] =   x[1]^4*x[2]*x[5]*x[3]
 pv[41] =   x[1]^8
 pv[42] =   x[1]^7*x[2]
 pv[43] =   x[1]^6*x[2]^2
 pv[44] =   x[1]^5*x[2]^3
 pv[45] =   x[1]^9
 pv[46] =   x[1]^8*x[2]

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
using Mimi
using Test

@testset "Mimi-SA" begin

    @info("test_empirical.jl")
    include("test_empirical.jl")

    @info("test_defmcs.jl")
    include("test_defmcs.jl")

    @info("test_defmcs_sobol.jl")
    include("test_defmcs_sobol.jl")

    @info("test_reshaping.jl")
    include("test_reshaping.jl")

    @info("test_payload.jl")
    include("test_payload.jl")

    @info("test_marginalmodel.jl")
    include("test_marginalmodel.jl")
end


text:
# Buffer
# ======

mutable struct Buffer
    data::Vector{UInt8}
    size::Int
end

Buffer(bufsize) = Buffer(Vector{UInt8}(undef, bufsize), 0)
Base.empty!(buffer::Buffer) = buffer.size = 0
Base.getindex(buffer::Buffer, i::Integer) = getindex(buffer.data, i)
Base.setindex!(buffer::Buffer, v, i::Integer) = setindex!(buffer.data, v, i)


text:
import Base: size, getindex, setindex!

#z indices are flipped for even-numbered timepoints
struct BidiImageArray{T} <: AbstractArray{T,4}
    A::AbstractArray{T,4}
    z_size::Int
end

BidiImageArray(A::AbstractArray{T,4}) where {T} = BidiImageArray(A, size(A,3))

size(B::BidiImageArray) = size(B.A)

Base.IndexStyle(::Type{<:BidiImageArray}) = IndexCartesian()

function getindex(B::BidiImageArray, I::Vararg{Int, 4})
    t_ind = I[4]
    z_ind = ifelse(isodd(t_ind), I[3], B.z_size-I[3]+1)
    return getindex(B.A, I[1], I[2], z_ind, t_ind)
end

function setindex!(B::BidiImageArray, v, I::Vararg{Int, 4})
    t_ind = I[4]
    z_ind = ifelse(isodd(t_ind), I[3], B.z_size-I[3]+1)
    setindex!(B.A, v, I[1], I[2], z_ind, t_ind)
end


text:
mutable struct Lattice
    basis::Array{Float64,2} #[12*3]; # fractional coords of each basis atom within unit cell (0 <= coord < 1)
    primitive::Vector{Float64} #[9]; # lattice <-> box transform matrices
    primitinv::Vector{Float64} #[9];
    rotaterow::Vector{Float64} #[9];
    rotatecol::Vector{Float64} #[9];
    spacing::Vector{Float64} #[3]; # lattice scale factors in 3 dims
    origin::Vector{Float64} #[3]; # lattice origin
    sublo::Vector{Float64} #[3];  # sub-box bounds in lattice space on this proc
    subhi::Vector{Float64} #[3];  # sub-box bounds in lattice space on this proc
    a1::Vector{Float64} #[3]; # edge vectors of unit cell
    a2::Vector{Float64} #[3]; # edge vectors of unit cell
    a3::Vector{Float64} #[3]; # edge vectors of unit cell
    scale::Float64

    atomtype::Vector{Int32}   # type of basis atoms
    orientx::Vector{Int32} #[3]; # lattice orientation vecs
    orienty::Vector{Int32} #[3]; # orientx = what lattice dir lies
    orientz::Vector{Int32} #[3]; #           along x dim in box

    style::Int32  # NONE,SC,FCC,etc
    spaceflag::Int32
    nbasis::Int32   # # of basis atoms in unit cell
    ilo::Int32
    ihi::Int32
    jlo::Int32
    jhi::Int32
    klo::Int32
    khi::Int32     # lattice bounds for sub-box on this proc
    natom::Int32                              # # of atoms in the sub-box on this proc

    Lattice() = new();
end

function printout(lat::Lattice)
    display("style, spaceflag, nbasis, natom, ilo, ihi, jlo, jhi, klo, khi, scale")
    display([lat.style, lat.spaceflag, lat.nbasis, lat.natom, lat.ilo, lat.ihi, lat.jlo, lat.jhi, lat.klo, lat.khi, lat.scale])
    display("origin: "); display(lat.origin);
    display("spacing: "); display(lat.spacing);
    display("orientx: "); display(lat.orientx);
    display("orienty: "); display(lat.orienty);
    display("orientz: "); display(lat.orientz);
    display("a1: "); display(lat.a1);
    display("a2: "); display(lat.a2);
    display("a3: "); display(lat.a3);
    display("sublo: "); display(lat.sublo);
    display("subhi: "); display(lat.subhi);
    display("type: "); display(lat.atomtype);
    display("basis: "); display(lat.atombasis);
    display("primitive: "); display(lat.primitive);
    display("primitinv: "); display(lat.primitinv);
    display("rotaterow: "); display(lat.rotaterow);
    display("rotatecol: "); display(lat.rotatecol);
end

function initlattice(stylein, scale, orientx=[1; 0; 0], orienty=[0; 1; 0], orientz=[0; 0; 1],
        spacing=[1.0; 1.0; 1.0], type=nothing, a1=[1.0; 0.0; 0.0], a2=[0.0; 1.0; 0.0], a3=[0.0; 0.0; 1.0],
        basis=nothing)

    lat = Lattice();

    lat.scale = scale;
    lat.spaceflag = 0;

    if (lowercase(stylein) == "none")
        lat.style = 0;
    elseif (lowercase(stylein) == "sc")
        lat.style = 1;
    elseif (lowercase(stylein) == "bcc")
        lat.style = 2;
    elseif (lowercase(stylein) == "fcc")
        lat.style = 3;
    elseif (lowercase(stylein) == "hcp")
        lat.style = 4;
    elseif (lowercase(stylein) == "diamond")
        lat.style = 5;
    elseif (lowercase(stylein) == "sq")
        lat.style = 6;
    elseif (lowercase(stylein) == "sq2")
        lat.style = 7;
    elseif (lowercase(stylein) == "hex")
        lat.style = 8;
    elseif (lowercase(stylein) == "custom")
        lat.style = 9;
    else
        error("Invalid lattice style");
    end

    lat.origin = [0.0; 0.0; 0.0];
    lat.orientx = orientx;
    lat.orienty = orienty;
    lat.orientz = orientz;
    lat.spacing = spacing;

    if (lowercase(stylein) == "custom")
        lat.a1 = a1;
        lat.a2 = a2;
        lat.a3 = a3;
        lat.basis = basis;
    else
        lat.a1 = [1.0; 0.0; 0.0];
        lat.a2 = [0.0; 1.0; 0.0];
        lat.a3 = [0.0; 0.0; 1.0];

        if (lowercase(stylein) === "hex")
            lat.a2[2] = sqrt(3.0);
        end
        if (lowercase(stylein) === "hcp")
            lat.a2[2] = sqrt(3.0);
            lat.a3[3] = sqrt(8.0/3.0);
        end

        if (lowercase(stylein) === "sc")
            lat.basis = [0 0 0];
        elseif (lowercase(stylein) === "bcc")
            lat.basis = [0 0 0; 0.5 0.5 0.5];
        elseif (lowercase(stylein) === "fcc")
            lat.basis = [0 0 0; 0.5 0.5 0.0; 0.5 0.0 0.5; 0.0 0.5 0.5];
        elseif (lowercase(stylein) === "hcp")
            lat.basis = [0 0 0; 0.5 0.5 0.0; 0.5 5.0/6.0 0.5; 0.0 1.0/3.0 0.5];
        elseif (lowercase(stylein) === "diamond")
            lat.basis = [0 0 0; 0.5 0.5 0.0; 0.5 0.0 0.5; 0.0 0.5 0.5; 0.25 0.25 0.25; 0.25 0.75 0.75; 0.75 0.25 0.75; 0.75 0.75 0.25];
        elseif (lowercase(stylein) === "sq")
            lat.basis = [0 0 0];
        elseif (lowercase(stylein) === "sq2")
            lat.basis = [0 0 0; 0.5 0.5 0.0];
        elseif (lowercase(stylein) === "hex")
            lat.basis = [0 0 0; 0.5 0.5 0.0];
        end
    end

    nbasis = size(lat.basis)[1];
    lat.basis = transpose(lat.basis)

    if type===nothing
        lat.atomtype = Int32.(ones(nbasis));
    else
        lat.atomtype = atomtype;
    end

    return lat
end

function latticebounds(lat::Lattice)
    lat.ilo = Int32(lat.sublo[1]) - 1;
    lat.jlo = Int32(lat.sublo[2]) - 1;
    lat.klo = Int32(lat.sublo[3]) - 1;
    lat.ihi = Int32(lat.subhi[1]) + 1;
    lat.jhi = Int32(lat.subhi[2]) + 1;
    lat.khi = Int32(lat.subhi[3]) + 1;

    if (lat.sublo[1] < 0.0)
        lat.ilo = lat.ilo - 1
    end
    if (lat.sublo[2] < 0.0)
        lat.jlo = lat.jlo - 1
    end
    if (lat.sublo[3] < 0.0)
        lat.klo = lat.klo - 1
    end

    lat.natom = 0;
    for k = lat.klo:lat.khi
        for j = lat.jlo:lat.jhi
            for i = lat.ilo:lat.ihi
                for m = 0:(lat.nbasis-1)
                    lat.natom = lat.natom + 1
                end
            end
        end
    end
end

#SetupLattice(lat::Lattice, unitstyle::Int32, dim::Int32)

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using Element
using Base.Test

using Element: lingrid, integrate, findsteps, eig, loggrid
using Element: SphericalHarmonicTable, radschrod1, cart2sph, realsphharm, enumeratelm, randp, radpoisson

# findsteps

f = (x,c) -> round(Int,x^2+c)
intervals = findsteps(f,0.0,2.0,args=1.0)
@test intervals == [(0.0,1.0),(1.0,1.5),(1.5,1.75),(1.75,2.0)]

# eig

x = lingrid(0,1,1000)
dx = x[2]-x[1]
v = 0*x
aE = E -> 2*(E-v)
qmin, qmax = 0,10

evs,efs,q = eig(aE,dx,qmin,qmax,Emin=0.0,Emax=20.0)
n = length(evs)
efs ./= sqrt(integrate(x, efs.^2))'  # normalization

exact_evs = [j^2*pi^2/2 for j in 1:n]
exact_efs = hcat([sqrt(2)*sin(j*pi*x) for j in 1:length(evs)]...)

@test q == collect(qmin:qmax)
@test_approx_eq integrate(x,efs.^2) fill(1.0,n)
@test_approx_eq_eps integrate(x, efs.^2 - exact_efs.^2) fill(0.0,n) 1e-12
@test_approx_eq_eps evs exact_evs 1e-5

# radschrod1

r = loggrid(-22,4,1000)
v = -1./r

n,l = 1,0
ev,ef = radschrod1(r,v,n,l)
@test_approx_eq_eps ev -1/(2*n^2) 1e-8
@test_approx_eq_eps integrate(r, (ef.^2 - (2r.*exp(-1*r)).^2).^2) 0.0 1e-12

n,l = 2,0
ev,ef = radschrod1(r,v,n,l)
@test_approx_eq_eps ev -1/(2*n^2) 1e-8
@test_approx_eq_eps integrate(r, (ef.^2 - (-1*r/sqrt(2).*(1-0.5*r).*exp(-0.5*r)).^2).^2) 0.0 1e-12

n,l = 2,1
ev,ef = radschrod1(r,v,n,l)
@test_approx_eq_eps ev -1/(2*n^2) 1e-8
@test_approx_eq_eps integrate(r, (ef.^2 - (r/3^0.5/2^(3./2).*r.*exp(-0.5*r)).^2).^2) 0.0 1e-12

n,l = 3,1
ev,ef = radschrod1(r,v,n,l)
@test_approx_eq_eps ev -1/(2*n^2) 1e-8
@test_approx_eq_eps integrate(r, (ef.^2 - (r*8./27/sqrt(6).*(1-r/6).*r.*exp(-1.*r/3)).^2).^2) 0.0 1e-12

# radpoisson
r = loggrid(-13,5,500)
rho = exp(-2r)

l = 0
vh_exact = 1/4*exp(-2r).*(1 + 2r) + (1 + exp(-2r).*(-1 - 2r.*(1 + r)))./(4r)  # l = 0
vh = radpoisson(r,rho,l)

@test_approx_eq_eps vh vh_exact 1e-5
@test_approx_eq_eps integrate(r, rho.*vh_exact/2) integrate(r, rho.*vh/2) 1e-6

l = 1
vh_exact = (3*exp(-2r).*(-1 + exp(2r) - 2r - 2r.^2))./(8r.^2)  # l = 1
vh = radpoisson(r,rho,l)

@test_approx_eq_eps vh vh_exact 1e-5
@test_approx_eq_eps integrate(r, rho.*vh_exact/2) integrate(r, rho.*vh/2) 1e-6

# spherical harmonics table

p = randp(3)[]
p = p/norm(p)
r,θ,ϕ = cart2sph(p)

lmax = 8
Y = SphericalHarmonicTable(lmax)
Y(p)  # precompute spherical harmonics

a = [realsphharm(l,m,θ,ϕ) for (l,m) in enumeratelm(lmax)]
b = [Y[l,m] for (l,m) in enumeratelm(lmax)]

@test_approx_eq a b

max_stars_repo_name: johncsnyder/Element


text:
dts = 2.0 .^ collect(-12:2:2) # collect(-9:1:3)
t_end = 1e3 # ms
in_current = 5. # 8

V_REST = -65. # * b2.mV
V_RESET = -65. # * b2.mV
FIRING_THRESHOLD = -50. # * b2.mV
MEMBRANE_RESISTANCE = 10. # * b2.Mohm
MEMBRANE_TIME_SCALE = 10. # * b2.ms
ABSOLUTE_REFRACTORY_PERIOD = 0. # * b2.ms

#doEuler_julia = false

###############################################################################
## EventBasedIntegrator
###############################################################################

using Pkg; Pkg.activate("./../../core/"); Pkg.instantiate()
push!(LOAD_PATH, string(pwd(),"./../../core/src/lifintegrator/core/"))
using EventBasedIntegrator, LinearAlgebra, ProgressMeter
using TimerOutputs, CPUTime, BenchmarkTools
using HDF5

include("./eventbasedlifintegrator_defs/eventbasedlifintegrator_defs.jl")
function single_neuron_julia(dt; doEuler = false)
    network_julia, spike_monitor_julia, profiling_time_julia =
        simulate_LIF_neuron(10 * in_current; # should actually be voltage (R = 1 Ω)
                            simulation_time = t_end,
                            doEuler = doEuler,
                            dt = dt,
                            v_rest=V_REST,
                            v_reset=V_RESET,
                            firing_threshold=FIRING_THRESHOLD,
                            membrane_resistance=MEMBRANE_RESISTANCE,
                            membrane_time_scale=MEMBRANE_TIME_SCALE,
                            abs_refractory_period=ABSOLUTE_REFRACTORY_PERIOD)

    spike_times_julia = [spike[2] for spike in spike_monitor_julia.spikes]
    return profiling_time_julia, spike_times_julia
end

# Event-based

fid = h5open("./data/spike_trains_single_neuron_julia.h5", "w")
write(fid, "dts", dts)
profiling_times_julia = zeros(length(dts))
for i in 1:length(dts)
    t, spts = single_neuron_julia(dts[i]; doEuler = false)
    write(fid, string("spike_train_",dts[i]), spts)
    profiling_times_julia[i] = t
end
write(fid, "profiling_times", profiling_times_julia)
close(fid)

# Euler-forward

fid = h5open("./data/spike_trains_single_neuron_julia_euler.h5", "w")
write(fid, "dts", dts)
profiling_times_julia = zeros(length(dts))
for i in 1:length(dts)
    t, spts = single_neuron_julia(dts[i]; doEuler = true)
    write(fid, string("spike_train_",dts[i]), spts)
    profiling_times_julia[i] = t
end
write(fid, "profiling_times", profiling_times_julia)
close(fid)

###############################################################################
## brian2 using C++ (cpp_standalone) and exact integration
###############################################################################

using PyCall

py"""
import sys
import brian2 as b2

sys.path.append('/Users/Bernd/Documents/PhD/Drafts/draft_mnistbenchmarks/BioPlausibleShallowDeepLearning/scripts/LIF_integrator_benchmarks/brian2_defs/')
from brian2_defs import simulate_LIF_neuron
"""

function single_neuron_brian2(dt; first_init = false)
    py"""
    import sys
    import brian2 as b2

    if $first_init:
        b2.set_device('cpp_standalone', build_on_run=False)
    else:
        #brian2.__init__.clear_cache()
        b2.device.reinit()
        b2.device.activate()
        b2.set_device('cpp_standalone', build_on_run=False)

    spike_monitor_brian2, profiling_summary = simulate_LIF_neuron($in_current * b2.namp,
                            simulation_time = $t_end * b2.ms,
                            dt=$dt,
                            v_rest=$V_REST * b2.mV,
                            v_reset=$V_RESET * b2.mV,
                            firing_threshold=$FIRING_THRESHOLD * b2.mV,
                            membrane_resistance=$MEMBRANE_RESISTANCE * b2.Mohm,
                            membrane_time_scale=$MEMBRANE_TIME_SCALE * b2.ms,
                            abs_refractory_period=$ABSOLUTE_REFRACTORY_PERIOD * b2.ms)
    """
    spike_monitor_brian2, profiling_summary = py"spike_monitor_brian2", py"profiling_summary"
    profiling_times = []
    for x in profiling_summary push!(profiling_times, x[2][1]) end
    profiling_time_brian2 = sum(profiling_times)

    spike_times_brian2 = 1e3 .* spike_monitor_brian2.spike_trains()[0]
    return profiling_time_brian2, spike_times_brian2
end

fid = h5open("./data/spike_trains_single_neuron_brian2.h5", "w")
write(fid, "dts", dts)
profiling_times_brian2 = zeros(length(dts))
for i in 1:length(dts)
    t, spts = single_neuron_brian2(dts[i]; first_init = (i == 1))
    write(fid, string("spike_train_",dts[i]), spts)
    profiling_times_brian2[i] = t
end
write(fid, "profiling_times", profiling_times_brian2)
close(fid)

###############################################################################

# n_spikes_julia =length(spike_times_julia)
# n_spikes_brian2 = length(spike_times_brian2)
# n_spikes = [n_spikes_julia, n_spikes_brian2]
# less_spikes = findmin(n_spikes)
#
# print(string("Time elapsed brian2 (profiling)[ms]: ", 1e3 .* profiling_time_brian2,"\n"))
# print(string("Time elapsed julia (median_of_profiling) [ms]: ", 1e3 .* profiling_time_julia,"\n"))
# print("\n")
# print(string("#spikes brian2: ", n_spikes_brian2,"\n"))
# print(string("#spikes julia: ", n_spikes_julia,"\n"))
# print(string("relative spike count discrepancy: ", abs(n_spikes_julia - n_spikes_brian2) * 100 ./ n_spikes_julia, " %"))
# print("\n")
# print(string("relative abs. error in spike times (Σ |errors| / #spikes * 1/ΔT = Σ |errors| / sim. time) between spike times: \n",
#     sum(abs.(spike_times_julia[1:less_spikes[1]] .- spike_times_brian2[1:less_spikes[1]])) * 100 ./ t_end, " %"))
# print("\n")


text:
userhome = ENV["HOME"]
usershell = ENV["SHELL"]
envregex = r"((^export JULIA)|(^JULIA))_(.*)=(.*)"
shellregex = r"(zsh$|bash$)"

matchshell = match(shellregex, usershell)
rcfile = "."*matchshell[1]*"rc"
shellrcfile = string(userhome)*"/"*string(rcfile)

#=
reading the shells rc file in $HOME and configure the environment with all variables set in the rc file
=#
open(shellrcfile) do f
    for i in eachline(f)
        matchregex = match(envregex, i)
        if matchregex !== nothing
            if matchregex[4] == "DEPOT_PATH"
                pushfirst!(DEPOT_PATH, string(matchregex[5]))
                ENV["JULIA_DEPOT_PATH"] = string(matchregex[5])
            end
        end
    end
end

#=
set DEPOT_PATH to users first slurm account $PC2FS directory if the DEPOT_PATH is not set explicity in the shells rc file
=#
get_groups_cmd = "/usr/bin/groups"
groups_regex = r"(HPC-[A-Z]{3}-.*|hpc-[a-z]{3}-.*|pc2-mitarbeiter|HPC-SNF)"
if length(DEPOT_PATH) <= 1
    groups = read(`$get_groups_cmd`, String)
    groups = replace(string(groups), " " => "\n")
    fgroups = []
    for i in eachline(IOBuffer(groups))
        matchregex = match(groups_regex, i)
        if matchregex !== nothing
            push!(fgroups, matchregex[1])
        end
    end

    if length(fgroups) >= 1
        username = ENV["USER"]
        user_depot_path = "/scratch/"*string(fgroups[1])*"/"*string(username)*"/.julia/"
        pushfirst!(DEPOT_PATH, string(user_depot_path))
        ENV["JULIA_DEPOT_PATH"] = string(user_depot_path)
    end
end


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

mutable struct PrivateEndpointConnectionListResult <: SwaggerModel
    value::Any # spec type: Union{ Nothing, Vector{PrivateEndpointConnection} } # spec name: value

    function PrivateEndpointConnectionListResult(;value=nothing)
        o = new()
        validate_property(PrivateEndpointConnectionListResult, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type PrivateEndpointConnectionListResult

const _property_map_PrivateEndpointConnectionListResult = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"))
const _property_types_PrivateEndpointConnectionListResult = Dict{Symbol,String}(Symbol("value")=>"Vector{PrivateEndpointConnection}")
Base.propertynames(::Type{ PrivateEndpointConnectionListResult }) = collect(keys(_property_map_PrivateEndpointConnectionListResult))
Swagger.property_type(::Type{ PrivateEndpointConnectionListResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PrivateEndpointConnectionListResult[name]))}
Swagger.field_name(::Type{ PrivateEndpointConnectionListResult }, property_name::Symbol) =  _property_map_PrivateEndpointConnectionListResult[property_name]

function check_required(o::PrivateEndpointConnectionListResult)
    true
end

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


text:
using Base: @kwdef

export BotAction
export NoAction
export Reply
export CreateTrialPost

abstract type BotAction end

struct NoAction <: BotAction end

@kwdef struct Reply <: BotAction
    content::AbstractString
end

@kwdef struct CreateTrialPost <: BotAction
    trialname::AbstractString
    trialtime::AbstractString
    trialdate::AbstractString
end