text:
function bounds_sanity_check(pomdp::P, b::WPFBelief{S}, L::Float64, U::Float64) where {S,P<:POMDP{S}}
    if L > U
        @warn(@sprintf("L (%-4.1f) > U (%-4.1f)", L, U))
        @info("Try e.g. `IndependentBounds(l, u, consistency_fix_thresh=1e-5)`.", maxlog=1)
    end
    if isinf(L) || isnan(L)
        @warn(@sprintf("L = %-4.1f. Infinite bounds are not supported.", L))
    end
    if isinf(U) || isnan(U)
        @warn(@sprintf("U = %-4.1f. Infinite bounds are not supported.", U))
    end
    return nothing
end

"""
    Dependent Bounds

Specify lower and upper bounds that are not independent.
###
It can take as input a function, `f(pomdp, belief)`, that returns both lower and upper bounds.
It can take as input any object for which a `bounds` function is implemented

"""

init_bounds(bds, pomdp, sol, rng) = bds

# Used to initialize lower and upper bounds with a single function
function bounds(f::Function, pomdp::P, b::WPFBelief{S,A,O}, max_depth::Int, bounds_warning::Bool) where {S,A,O,P<:POMDP{S,A,O}}
    l, u = f(pomdp, b)
    if bounds_warning
        bounds_sanity_check(pomdp, b, l, u)
    end
    return l, u
end

# Used to initialize both the lower and upper bound with an object for which a `bounds` function is implemented
function bounds!(L::Vector{Float64}, U::Vector{Float64}, bd::B, pomdp::P, b::WPFBelief{S,A,O}, W::Vector{Vector{Float64}}, obs::Vector{O}, max_depth::Int, bounds_warning::Bool) where {S,A,O,P<:POMDP{S,A,O},B}
    @inbounds for i in eachindex(W)
        switch_to_sibling!(b, obs[i], W[i])
        L[i], U[i] = bounds(bd, pomdp, b, max_depth, bounds_warning)
    end
    return L, U
end

"""
    IndependentBounds(lower, upper, check_terminal=false, consistency_fix_thresh=0.0)

Specify lower and upper bounds that are independent of each other (the most common case).
A lower or upper bound can be a Number, a Function, `f(pomdp, belief)`, that returns a bound, an object for which a `bound` function
is implemented. Specifically, for FOValue, POValue, FORollout, SemiPORollout, and PORollout, a `bound` function is already implemented.
You can also implement a `bound!` function to initialize sibling beliefs simultaneously, if it could provide a further performace gain.

# Keyword Arguments
- `check_terminal::Bool=false`: if true, then if all the states in the belief are terminal, the upper and lower bounds will be overridden and set to 0. (deprecated option)
- `consistency_fix_thresh::Float64=0.0`: if `upper < lower` and `upper >= lower-consistency_fix_thresh`, then `upper` will be bumped up to `lower`.
"""

mutable struct IndependentBounds{L, U}
    lower::L
    upper::U
    check_terminal::Bool
    consistency_fix_thresh::Float64
end

function IndependentBounds(l, u;
                           check_terminal=false,
                           consistency_fix_thresh=0.0)
    return IndependentBounds(l, u, check_terminal, consistency_fix_thresh)
end

function init_bounds(bds::IndependentBounds, pomdp::POMDP, sol::AdaOPSSolver, rng::AbstractRNG)
    if bds.check_terminal
        @warn("`check_terminal` is deprecated. It is assured that at least one state in the belief is not terminal.")
    end
    return IndependentBounds(convert_estimator(bds.lower, sol, pomdp),
                             convert_estimator(bds.upper, sol, pomdp),
                             bds.check_terminal,
                             bds.consistency_fix_thresh,
                            )
end

function bounds(bds::IndependentBounds, pomdp::P, b::WPFBelief{S,A,O}, max_depth::Int, bounds_warning::Bool) where {S,A,O,P<:POMDP{S,A,O}}
    l = bound(bds.lower, pomdp, b, max_depth)
    u = bound(bds.upper, pomdp, b, max_depth)
    if u < l && u >= l-bds.consistency_fix_thresh
        u = l
    end
    if bounds_warning
        bounds_sanity_check(pomdp, b, l, u)
    end
    return l, u
end

function bounds!(L::Vector{Float64}, U::Vector{Float64}, bds::IndependentBounds{LB,UB}, pomdp::P, b::WPFBelief{S,A,O}, W::Vector{Vector{Float64}}, obs::Vector{O}, max_depth::Int, bounds_warning::Bool) where {LB,UB,S,A,O,P<:POMDP{S,A,O}}
    bound!(L, bds.lower, pomdp, b, W, obs, max_depth)
    bound!(U, bds.upper, pomdp, b, W, obs, max_depth)
    @inbounds for i in eachindex(W)
        if U[i] < L[i] && U[i] >= L[i]-bds.consistency_fix_thresh
            U[i] = L[i]
        end
    end
    if bounds_warning
        @inbounds for i in eachindex(W)
            switch_to_sibling!(b, obs[i], W[i])
            bounds_sanity_check(pomdp, b, L[i], U[i])
        end
    end
    return L, U
end

# Used when the lower or upper bound is a fixed number
bound(n::N, pomdp, b, max_depth) where N<:Real = convert(Float64, n)
bound!(V::Vector{Float64}, n::Float64, pomdp::P, b::WPFBelief{S,A,O}, W::Vector{Vector{Float64}}, obs::Vector{O}, max_depth::Int) where {S,A,O,P<:POMDP{S,A,O}} = fill!(V, n)
bound!(V::Vector{Float64}, n::Int, pomdp::P, b::WPFBelief{S,A,O}, W::Vector{Vector{Float64}}, obs::Vector{O}, max_depth::Int) where {S,A,O,P<:POMDP{S,A,O}} = fill!(V, convert(Float64, n))

# Used when lower or upper bound is a function
bound(f::Function, pomdp::P, b::WPFBelief{S,A,O}, max_depth::Int) where {S,A,O,P<:POMDP{S,A,O}} = f(pomdp, b)

# Used when the lower or upper bound is an object for which a `bound` function is implemented
function bound!(V::Vector{Float64}, bd::B, pomdp::P, b::WPFBelief{S,A,O}, W::Vector{Vector{Float64}}, obs::Vector{O}, max_depth::Int) where {S,A,O,P<:POMDP{S,A,O},B}
    @inbounds for i in eachindex(W)
        switch_to_sibling!(b, obs[i], W[i])
        V[i] = bound(bd, pomdp, b, max_depth)
    end
    return V
end

struct SolvedFOValue{P<:Policy}
    policy::P
    values::Vector{Float64}
end

struct SolvedFORollout{P<:Policy, RNG<:AbstractRNG}
    policy::P
    values::Vector{Float64}
    rng::RNG
end

struct SolvedPORollout{P<:Policy, U<:Updater, RNG<:AbstractRNG}
    policy::P
    values::Vector{Float64}
    updater::U
    rng::RNG
end

struct POValue
    solver::Union{Solver, Policy}
end

struct SolvedPOValue{P<:Policy}
    policy::P
end

struct SemiPORollout
    solver::Union{Solver, Policy}
end

mutable struct SolvedSemiPORollout{S, O, P<:Policy, RNG<:AbstractRNG}
    policy::P
    obs_ind_dict::Vector{Dict{O, Int}}
    probs::Vector{Vector{Float64}}
    states_list::Vector{Vector{Vector{S}}}
    weights_list::Vector{Vector{Vector{Float64}}}
    states::Vector{Vector{S}}
    weights::Vector{Vector{Float64}}
    rng::RNG
end

function bound(bd::SolvedFORollout, pomdp::POMDP, b::WPFBelief{S}, max_depth::Int) where S
    resize!(bd.values, n_particles(b))
    broadcast!((s)->estimate_value(bd, pomdp, s, b, max_depth-b.depth), bd.values, particles(b))
    return dot(weights(b), bd.values)/weight_sum(b)
end

function bound!(V::Vector{Float64}, bd::SolvedFORollout, pomdp::POMDP, b::WPFBelief{S}, W::Vector{Vector{Float64}}, obs::Vector, max_depth::Int) where S
    resize!(bd.values, n_particles(b))
    broadcast!((s)->estimate_value(bd, pomdp, s, b, max_depth-b.depth), bd.values, particles(b))
    @inbounds for i in eachindex(W)
        V[i] = dot(bd.values, W[i]) / sum(W[i])
    end
    return V
end

function bound(bd::SolvedPORollout{P}, pomdp::M, b::WPFBelief{S,A,O}, max_depth::Int) where {P,S,A,O,M<:POMDP{S,A,O}}
    resize!(bd.values, n_particles(b))
    broadcast!((s)->estimate_value(bd, pomdp, s, b, max_depth-b.depth), bd.values, particles(b))
    return dot(weights(b), bd.values)/weight_sum(b)
end

function bound(bd::SolvedSemiPORollout{S,O,P}, pomdp::M, b::WPFBelief{S,A,O}, max_depth::Int) where {P,S,A,O,M<:POMDP{S,A,O}}
    return estimate_value(bd, pomdp, b, max_depth-b.depth, 0, 0)
end

function bound(bd::SolvedFOValue{P}, pomdp::M, b::WPFBelief{S}, max_depth::Int) where {P,S,M<:POMDP{S}}
    resize!(bd.values, n_particles(b))
    broadcast!((s)->value(bd.policy, s), bd.values, particles(b))
    return dot(weights(b), bd.values)/weight_sum(b)
end

function bound!(V::Vector{Float64}, bd::SolvedFOValue{P}, pomdp::M, b::WPFBelief{S}, W::Vector{Vector{Float64}}, obs::Vector, max_depth::Int) where {P,S,M<:POMDP{S}}
    resize!(bd.values, n_particles(b))
    broadcast!((s)->value(bd.policy, s), bd.values, particles(b))
    @inbounds for i in eachindex(W)
        V[i] = dot(bd.values, W[i]) / sum(W[i])
    end
    return V
end

bound(bd::SolvedPOValue, pomdp::POMDP, b::WPFBelief, max_depth::Int) = value(bd.policy, b)

# Convert an unsolved estimator to solved estimator
convert_estimator(ev, sol, mdp) = ev

function convert_estimator(est::FOValue, sol::AdaOPSSolver, pomdp::POMDP)
    policy = MCTS.convert_to_policy(est.solver, UnderlyingMDP(pomdp))
    SolvedFOValue(policy, sizehint!(Float64[], sol.m_max))
end

function convert_estimator(est::FORollout, sol::AdaOPSSolver, pomdp::POMDP)
    policy = MCTS.convert_to_policy(est.solver, UnderlyingMDP(pomdp))
    SolvedFORollout(policy, sizehint!(Float64[], sol.m_max), sol.rng)
end

function convert_estimator(est::POValue, sol::AdaOPSSolver, pomdp::POMDP)
    policy = MCTS.convert_to_policy(est.solver, pomdp)
    SolvedPOValue(policy)
end

function convert_estimator(est::PORollout, sol::AdaOPSSolver, pomdp::POMDP)
    policy = MCTS.convert_to_policy(est.solver, pomdp)
    SolvedPORollout(policy, sizehint!(Float64[], sol.m_max), est.updater, sol.rng)
end

function convert_estimator(est::RolloutEstimator, sol::AdaOPSSolver, pomdp::POMDP)
    policy = MCTS.convert_to_policy(est.solver, pomdp)
    SolvedPORollout(policy, sizehint!(Float64[], sol.m_max), updater(policy), sol.rng)
end

function convert_estimator(est::SemiPORollout, sol::AdaOPSSolver, pomdp::POMDP{S,A,O}) where {S,A,O}
    policy = MCTS.convert_to_policy(est.solver, pomdp)
    # this can be further optimized by memory preallocation - To Do
    obs_ind_dict = Dict{O, Int}[]
    probs = Vector{Float64}[]
    states_list = Vector{Vector{S}}[]
    weights_list = Vector{Vector{Float64}}[]
    states = Vector{S}[]
    weights = Vector{Float64}[]
    for i in 1:sol.max_depth
        push!(obs_ind_dict, Dict{O, Int}())
        push!(probs, Float64[])
        push!(states_list, Vector{S}[])
        push!(weights_list, Vector{Float64}[])
    end
    for i in 1:max(sol.max_depth, 2*sol.m_max)
        push!(states, S[])
        push!(weights, Float64[])
    end
    SolvedSemiPORollout(policy, obs_ind_dict, probs, states_list, weights_list, states, weights, sol.rng)
end

# Estimate the value of state with estimator
function estimate_value(estimator::Union{SolvedPORollout,SolvedFORollout}, pomdp::POMDP, start_state, b::WPFBelief, steps::Int)
    rollout(estimator, pomdp, start_state, b, steps)
end

POMDPLinter.@POMDP_require estimate_value(estimator::Union{SolvedPORollout,SolvedFORollout}, pomdp::POMDP, start_state, b::WPFBelief, steps::Int) begin
    @subreq rollout(estimator, pomdp, start_state, b, steps)
end

# Perform a rollout simulation to estimate the value.
function rollout(est::SolvedPORollout, pomdp::POMDP, start_state, b::WPFBelief, steps::Int)
    sim = RolloutSimulator(est.rng, steps)
    return simulate(sim, pomdp, est.policy, est.updater, b, start_state)
end

POMDPLinter.@POMDP_require rollout(est::SolvedPORollout, pomdp::POMDP, start_state, b::WPFBelief, steps::Int) begin
    sim = RolloutSimulator(est.rng, steps)
    @subreq simulate(sim, pomdp, est.policy, est.updater, b, start_state)
end

function rollout(est::SolvedFORollout, pomdp::POMDP, start_state, b::WPFBelief, steps::Int)
    sim = RolloutSimulator(est.rng, steps)
    return simulate(sim, pomdp, est.policy, start_state)
end

POMDPLinter.@POMDP_require rollout(est::SolvedFORollout, pomdp::POMDP, start_state, b::WPFBelief, steps::Int) begin
    sim = RolloutSimulator(est.rng, steps)
    @subreq simulate(sim, pomdp, est.policy, start_state)
end

function estimate_value(est::SolvedSemiPORollout{S,O,P}, pomdp::M, b::AbstractParticleBelief{S}, steps::Int, inner_ind::Int, leaf_ind::Int) where {P,S,A,O,M<:POMDP{S,A,O}}
    if steps <= 0 || weight_sum(b) == 0.0
        return 0.0
    end
    inner_ind += 1
    obs_ind_dict = est.obs_ind_dict[inner_ind]
    probs = est.probs[inner_ind]
    states_list = est.states_list[inner_ind]
    weights_list = est.weights_list[inner_ind]
    empty!(obs_ind_dict)
    empty!(probs)
    empty!(states_list)
    empty!(weights_list)

    a = action(est.policy, b)

    r_sum = 0.0
    for (k, s) in enumerate(particles(b))
        if !isterminal(pomdp, s)
            sp, o, r = @gen(:sp, :o, :r)(pomdp, s, a, est.rng)

            if !haskey(obs_ind_dict, o)
                push!(probs, 0.0)
                obs_ind_dict[o] = length(probs)
                leaf_ind += 1
                empty!(est.states[leaf_ind])
                empty!(est.weights[leaf_ind])
                push!(states_list, est.states[leaf_ind])
                push!(weights_list, est.weights[leaf_ind])
            end
            obs_ind = obs_ind_dict[o]
            w = weight(b, k)
            push!(states_list[obs_ind], sp)
            push!(weights_list[obs_ind], w * obs_weight(pomdp, s, a, sp, o))
            probs[obs_ind] += w

            r_sum += r * w
        end
    end

    U = 0.0
    for (o, obs_ind) in obs_ind_dict
        if length(states_list[obs_ind]) == 1
            U += probs[obs_ind] * rollout(pomdp, est.policy, ParticleCollection(states_list[obs_ind]), steps-1, est.rng)
        else
            U += probs[obs_ind] * estimate_value(est, pomdp, WeightedParticleBelief(states_list[obs_ind], weights_list[obs_ind]), steps-1, inner_ind, leaf_ind)
        end
    end
    return (r_sum + discount(pomdp)*U)/weight_sum(b)
end

function rollout(pomdp::POMDP{S}, policy::Policy, b::ParticleCollection{S}, steps::Integer, rng::AbstractRNG) where S
    @assert n_particles(b) == 1
    disc = 1.0
    r_total = 0.0
    mem = particles(b)
    s = mem[1]
    while !isterminal(pomdp, s) && steps > 0
        a = action(policy, b)
        s, r = @gen(:sp, :r)(pomdp, s, a, rng)
        r_total += disc*r
        mem[1] = s
        if b._probs !== nothing
            empty!(b._probs)
            b._probs[s] = 1.0
        end
        disc *= discount(pomdp)
        steps -= 1
    end
    return r_total
end


text:
"""

"""
module single_element_examples

using LinearAlgebra
using Arpack
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json, plot_triads

function standard_single_dsg3()
    E = 200e3*phun("MPa")
    nu = 0.3;
    rho= 8000*phun("KG/M^3");
    thickness = 0.1*phun("m");
    L = 1.0*phun("m");

    # Mesh
    tolerance = L/1000
    fens, fes = T3block(L,L,1,1);
    fes = subset(fes, [1])
    fens.xyz = xyz3(fens)

    connected = findunconnnodes(fens, fes);
    fens, new_numbering = compactnodes(fens, connected);
    fes = renumberconn!(fes, new_numbering);
    @show fens.xyz
    @show fes

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    formul = FEMMShellT3FFModule

    # @show Bb[:, [4, 5, 10, 11, 16, 17]]
    # @show Bs[:, [4, 5, 10, 11, 16, 17]]

    # Bb[:, [4, 5, 10, 11, 16, 17]] = [
    # 0.0 -1.0 0.0 1.0 0.0 0.0;
    # 1.0 0.0 0.0 0.0 -1.0 0.0;
    # 1.0 -1.0 -1.0 0.0 0.0 1.0]
    # Bs[:, [4, 5, 10, 11, 16, 17]] = [
    # -1/6 1/3 1/6 1/2 0.0 1/6;
    # -1/3 1/6 -1/6 0.0 -1/2 -1/6]

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    associate = formul.associategeometry!
        stiffness = formul.stiffness
        mass = formul.mass

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # No EBC's
    l1 = collect(1:count(fens))
    for i in [1, 2, 3, 6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);
    @show dchi.dofnums

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    K = Matrix(K)
    dec = eigen(K)
    @show dec.values
    v = dec.vectors

    # mfemm = FEMMDeforLinear(DeforModelRed3D, IntegDomain(fes, TriRule(1), thickness), mater)
    # M = mass(mfemm, geom0, dchi);

    # # Solve1
    # OmegaShift = 0.1*2*pi
    # neigvs = 10
    # d, v, nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
    # d[:] = d .- OmegaShift;
    # fs = real(sqrt.(complex(d)))/(2*pi)
    # @show fs

    for j in 1:6
        @show dec.values[j]
        @show round.(v[:, j], digits=4)
    end

    # # Visualization
    for j in 1:length(dec.values)
        U = v[:, j]
        scattersysvec!(dchi, (L/4)/maximum(abs.(U)).*U)
        Rfield = deepcopy(Rfield0)
        update_rotation_field!(Rfield, dchi)
        plots = cat(plot_space_box([[-L -L -L]; [+L +L +L]]),
            plot_nodes(fens),
            plot_triads(fens; triad_length = 0.2, x = geom0.values, u = dchi.values[:, 1:3], R = Rfield.values),
            plot_triads(fens; triad_length = 0.4, x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values),
            plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values)
            ;
            dims = 1)
        pl = render(plots, title = "$j: Eigenvalue $(dec.values[j])")
    end
end

function distorted_single_dsg3()
    E = 200e3*phun("MPa")
    nu = 0.3;
    rho= 8000*phun("KG/M^3");
    thickness = 0.1*phun("m");
    L = 1.0*phun("m");

    # Mesh
    tolerance = L/1000
    fens, fes = T3block(L,L,1,1);
    fes = subset(fes, [1])
    fens.xyz = [
    0.0 -10.0 0.0;
    13.0 -3.0 0.0;
    -1.0 7.0 0.0]
    fens.xyz = xyz3(fens)

    connected = findunconnnodes(fens, fes);
    fens, new_numbering = compactnodes(fens, connected);
    fes = renumberconn!(fes, new_numbering);
    @show fens.xyz
    @show fes

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    formul = FEMMShellCSDSG3Module
    formul = FEMMShellT3DSGModule

    # @show Bb[:, [4, 5, 10, 11, 16, 17]]
    # @show Bs[:, [4, 5, 10, 11, 16, 17]]

    # Bb[:, [4, 5, 10, 11, 16, 17]] = [
    # 0.0 -1.0 0.0 1.0 0.0 0.0;
    # 1.0 0.0 0.0 0.0 -1.0 0.0;
    # 1.0 -1.0 -1.0 0.0 0.0 1.0]
    # Bs[:, [4, 5, 10, 11, 16, 17]] = [
    # -1/6 1/3 1/6 1/2 0.0 1/6;
    # -1/3 1/6 -1/6 0.0 -1/2 -1/6]

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    associate = formul.associategeometry!
        stiffness = formul.stiffness
        mass = formul.mass

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # No EBC's
    l1 = collect(1:count(fens))
    for i in [1, 2, 3, 6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);
    @show dchi.dofnums

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    K = Matrix(K)
    dec = eigen(K)
    @show dec.values
    v = dec.vectors

    # mfemm = FEMMDeforLinear(DeforModelRed3D, IntegDomain(fes, TriRule(1), thickness), mater)
    # M = mass(mfemm, geom0, dchi);

    # # Solve1
    # OmegaShift = 0.1*2*pi
    # neigvs = 10
    # d, v, nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
    # d[:] = d .- OmegaShift;
    # fs = real(sqrt.(complex(d)))/(2*pi)
    # @show fs

    for j in 1:6
        @show dec.values[j]
        @show round.(v[:, j], digits=4)
    end

    box = reshape(boundingbox(fens.xyz), 1, 6)
    box = vcat(box[1, [1, 3, 5]]' .+ [-5, -5, -5], box[1, [2, 4, 6]]' .- [-5, -5, -5])
    # # Visualization
    for j in 1:length(dec.values)
        U = v[:, j]
        scattersysvec!(dchi, (L/4)/maximum(abs.(U)).*U)
        Rfield = deepcopy(Rfield0)
        update_rotation_field!(Rfield, dchi)
        plots = cat(plot_space_box(box),
            plot_nodes(fens),
            plot_triads(fens; triad_length = 2.2, x = geom0.values, u = dchi.values[:, 1:3], R = Rfield.values),
            plot_triads(fens; triad_length = 3.4, x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values),
            plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values)
            ;
            dims = 1)
        pl = render(plots, title = "$j: Eigenvalue $(dec.values[j])")
    end
end

function allrun()
    println("#####################################################")
    println("# standard_single_dsg3 ")
    standard_single_dsg3()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing

ext: jl


text:
#Doing things another way, even if not better, not over thinking this, which is not what a programming mindset is supposed to do, but jumping in.

#produce a key which counts the total number of significant counts for each pair across the windows were found
function dictTotalThresholdsAdjListWindowCount(winAggDict,startYr,endYr,windowSize)
    #a new dictionary which is non-windows
    #winDictTotal = Dict()

    #create the structure of the adjacency list where we will add the counts found to for each pair
    #produce an initialized country pairing adjacency list for the threshold passings, the window interval is not used, as we need the full country pairings, not the window subset
    thresholdSignificantAdjListTOTAL = initCountryPairAdjList(startYr,endYr)

    #loop through the keys of the Dict to add the significant pairs to the
    yr = startYr
    while( (yr+windowSize) <= endYr )

        threshSigAdjListTMP = winAggDict["$(yr)-$(yr+windowSize)"]["thresholdSigAdjList"]
        for rowIndWin in 1:size(threshSigAdjListTMP,1)#look through every pairing in this window
            for rowIndTotal in 1:size(thresholdSignificantAdjListTOTAL,1)
                if((thresholdSignificantAdjListTOTAL[rowIndTotal,1] == threshSigAdjListTMP[rowIndWin,1]) && (thresholdSignificantAdjListTOTAL[rowIndTotal,2] == threshSigAdjListTMP[rowIndWin,2]))
                    thresholdSignificantAdjListTOTAL[rowIndTotal,3] += threshSigAdjListTMP[rowIndWin,3]
                end

            end

        end

    	yr = yr + windowSize
    end

    #add the accumulated significance total to the dictionary
    winAggDict["thresholdSignificantAdjListTOTAL"] = thresholdSignificantAdjListTOTAL

    #also add the full total country pairing name list for the window set
    winAggDict["countriesNamesTotal"] = subsetCountryNamesArray(startYr,endYr)

    return winAggDict
end

#use the windowConf dictionary which has the threshold values for significance of the average values to determine for each window which country pairings were significant in a new window key 'thresholdSignificantAdjList'
function windowsDictThresholdsAdjList(windowConf, startYr = 1980, endYr = 1990, windowSize = 5)
    #from the csv data build the dictionary for the basic information on the country pair lists, aggregate scores and average of the aggregate scores
    winAggDict = windowsDictScoreAdjList(startYr, endYr, windowSize)
    #println(keys(winAggDict["1980-1985"]))#["1980-1985"]["countries"])
    #println(winAggDict["1980-1985"]["avgScoreAggregateAdjList"])
    #windowConf is a dictionary with the window keys as winDicts here, and each window has a particular threshold
    #extract the threshold, extract the average score matrix, produce an init country pairing adjacency matrix and then put it as a threshold passings component to the dictionary
    yr = startYr
    while( (yr+windowSize) <= endYr )
        #get the threshold value (float)
        thresholdTmp = windowConf["$(yr)-$(yr+windowSize)"]
        #extract the average of the aggregate (adjacency list of country pairs)
        avgAggAdjList = winAggDict["$(yr)-$(yr+windowSize)"]["avgScoreAggregateAdjList"]
        #produce an initialized country pairing adjacency list for the threshold passings
        thresholdSignificantAdjList = initCountryPairAdjList(yr,yr+windowSize)#(startYr,endYr)#!!!!!XXXX
        #now set the pairs of country rows for each significant pair to 1 from 0
        #HERE we do the THRESHOLD COMPARISON
        if(windowConf["tailSide"]=="upper" || windowConf["tailSide"]=="right")
            #println("upper")
            for rowInd in 1:size(avgAggAdjList,1)
                if(avgAggAdjList[rowInd,3] >= thresholdTmp)
                    thresholdSignificantAdjList[rowInd,3] = 1
                end
            end
        else
            #println("lower")
            for rowInd in 1:size(avgAggAdjList,1)
                if(avgAggAdjList[rowInd,3] <= thresholdTmp)
                    thresholdSignificantAdjList[rowInd,3] = 1
                end
            end
        end
        #=
        for rowInd in 1:size(avgAggAdjList,1)
            if(avgAggAdjList[rowInd,3] >= thresholdTmp)
                thresholdSignificantAdjList[rowInd,3] = 1
            end
        end
        =#
        winAggDict["$(yr)-$(yr+windowSize)"]["thresholdSigAdjList"] = thresholdSignificantAdjList
	yr = yr + windowSize
    end
    return winAggDict

end

#Given the total time interval, and windowSize, produce a dictionary for each window
#Returned is a dictionary for the time windows
#"countries" key is the unique list of country names for that window
#"scoreAggregateAdjList" key is the aggregate for the country pairs in each window
function windowsDictScoreAdjList(startYr = 1980, endYr = 1990, windowSize = 5)

    winDict = Dict()
    yr = startYr
    while( (yr+windowSize) <= endYr )
	winDict["$(yr)-$(yr+windowSize)"] = Dict()
	winDict["$(yr)-$(yr+windowSize)"]["countries"] = subsetCountryNamesArray(yr,yr+windowSize)#(startYr,endYr)
        aggAdjList = aggregateAdjList(yr,yr+windowSize)#(startYr,endYr)
	winDict["$(yr)-$(yr+windowSize)"]["scoreAggregateAdjList"] = aggAdjList
        #compute the average over the window for the scores
        avgAggAdjList = initCountryPairAdjList(yr,yr+windowSize)#(startYr,endYr)
        avgAggAdjList[:,3] = aggAdjList[:,3].*(1/(windowSize+1))
        winDict["$(yr)-$(yr+windowSize)"]["avgScoreAggregateAdjList"] = avgAggAdjList

	yr = yr + windowSize
    end
    return winDict
end

#create and return an adjacency country pair list with initially zero for the scores
function initCountryPairAdjList(startYr = 1980, endYr = 1990)

    adjList = yearsScoreAdjList(startYr, endYr)
    cntryNames = subsetCountryNamesArray(startYr,endYr)
    countryNum = length(cntryNames)
    aggregateAdjList = Array{Any}(undef,countryNum^2 - countryNum, 3)
    #INIT: fill the names and initial aggregate scores
    tmpRow = 1
    for ii in 1:length(cntryNames) #over each first cycle of countries
        for jj in 1:length(cntryNames) #over second cycle of countries
            if(ii != jj)
                aggregateAdjList[tmpRow,:] = [cntryNames[ii] cntryNames[jj] 0]
                tmpRow = tmpRow + 1
            end
        end
    end
    return aggregateAdjList
end

#produce the aggregate score for the 'to-and-from' from the adjacency list so that we have the unique pairings and total scores
function aggregateAdjList(startYr = 1980, endYr = 1990)

    adjList = yearsScoreAdjList(startYr, endYr)
    cntryNames = subsetCountryNamesArray(startYr,endYr)
    countryNum = length(cntryNames)
    aggregateAdjList = Array{Any}(undef,countryNum^2 - countryNum, 3)
    #INIT: fill the names and initial aggregate scores
    tmpRow = 1
    for ii in 1:length(cntryNames) #over each first cycle of countries
        for jj in 1:length(cntryNames) #over second cycle of countries
            if(ii != jj)
                aggregateAdjList[tmpRow,:] = [cntryNames[ii] cntryNames[jj] 0]
                tmpRow = tmpRow + 1
            end
        end
    end
    #FILL: the scores for the AdjList
    adjListTotal = yearsScoreAdjList(startYr, endYr)
    for ii in 1:size(aggregateAdjList,1)#look at every country pair
        for rowInd in 1:size(adjListTotal,1)
            if((aggregateAdjList[ii,1]==adjListTotal[rowInd,1]) && (aggregateAdjList[ii,2]==adjListTotal[rowInd,2]))
                aggregateAdjList[ii,3] = aggregateAdjList[ii,3] + adjListTotal[rowInd,4]
            end
        end
    end
    return aggregateAdjList
end

#look at the directory and load every year in the range of the .csv files, then produce a final matrix(adjacency list) which is every year in the countryFrom|countryTo|year|score for this 'subset'
function yearsScoreAdjList(startYr = 1980, endYr = 1990)
    dirFiles = readdir("./dataTables/")
    adjMatScore = []
    for dF in dirFiles
        yrTmp = parse(Int,((split(dF,"."))[1]))
        if(yrTmp < startYr || yrTmp > endYr)
            continue
        end
        fileTmp = open(string("./dataTables/",dF))
        linesTmp = readlines(fileTmp)#read each file lines
        origColNames = split(linesTmp[1],r",|\n",keepempty=false)#get line1 into components
        for ii=2:length(linesTmp)#get the names on the rows
            rowTmp = split(linesTmp[ii],r",|\n",keepempty=false)
            for jj=2:length(origColNames)
                if(ii!=jj)
                    #print("$(rowTmp[1])->$(origColNames[jj])=$(rowTmp[jj]),")
                    if(isempty(adjMatScore) == true)
                        adjMatScore = Array{Any}(undef,1,4)
                        adjMatScore[1,:] = [rowTmp[1],origColNames[jj],yrTmp,parse(Int,rowTmp[jj])]
                    else
                        #println([rowTmp[1],origColNames[jj],yrTmp,parse(Int,rowTmp[jj])])
                        adjMatScore = vcat(adjMatScore , [rowTmp[1] origColNames[jj] yrTmp parse(Int,rowTmp[jj])])
                    end
                end
            end
        end

    end
    adjMatScoreNoShows = missingCountryScoreInsert(adjMatScore)
    return adjMatScoreNoShows
    #return adjMatScore
end

#in this function we will pass the adjMatScore structure (mat)
#in each year, find the average score for all the countries per year
#find countries which are not in each year
#add countries and the average score for each year
function missingCountryScoreInsert(adjMatScore)
    adjMatScoreTmp = adjMatScore#what is returned
    countriesUnique = unique(adjMatScore[:,1])
    yearsUnique = unique(adjMatScore[:,3])
    for yy in yearsUnique

        rowsYear = (adjMatScore[:,3] .== yy)
        yrCountries = unique(adjMatScore[rowsYear,1])
        #println(yrCountries)
        ind = 0
        for cc in countriesUnique
            if( !(cc in yrCountries) )
                ind+=1
            end
        end
        sumTmp = sum(adjMatScore[rowsYear,4])
        muTmp = round(sumTmp/(length(adjMatScore[rowsYear,4])+ind),digits=4)
        #println(muTmp)
        muTmp = convert(Int,round(muTmp))
        #println(muTmp)
        for cc in countriesUnique
            if( !(cc in yrCountries) )
                #println("$(cc)--$(yy)")

                for cTmp in countriesUnique
                    #println([cc cTmp yy muTmp])
                    #println([cTmp cc yy muTmp])
                    adjMatScoreTmp = vcat(adjMatScoreTmp,[cc cTmp yy muTmp])
                    adjMatScoreTmp = vcat(adjMatScoreTmp,[cTmp cc yy muTmp])
                end
            end
        end

    end
    return adjMatScoreTmp

end

#look at the directory and load every .csv file, then produce a final matrix which is the countryFrom|countryTo|year|score
function totalScoreAdjList()
    dirFiles = readdir("./dataTables/")
    adjMatScore = []
    ind = 1
    for dF in dirFiles
        yrTmp = parse(Int,((split(dF,"."))[1]))

        fileTmp = open(string("./dataTables/",dF))
        linesTmp = readlines(fileTmp)#read each file lines
        origColNames = split(linesTmp[1],r",|\n",keepempty=false)#get line1 into components
        for ii=2:length(linesTmp)#get the names on the rows
            rowTmp = split(linesTmp[ii],r",|\n",keepempty=false)
            for jj=2:length(origColNames)
                if(ii!=jj)
                    #print("$(rowTmp[1])->$(origColNames[jj])=$(rowTmp[jj]),")
                    if(isempty(adjMatScore) == true)
                        adjMatScore = Array{Any}(1,4)
                        adjMatScore[1,:] = [rowTmp[1],origColNames[jj],yrTmp,parse(Int,rowTmp[jj])]
                    else
                        #println([rowTmp[1],origColNames[jj],yrTmp,parse(Int,rowTmp[jj])])
                        adjMatScore = vcat(adjMatScore , [rowTmp[1] origColNames[jj] yrTmp parse(Int,rowTmp[jj])])
                    end
                end
            end
        end

    end
    return adjMatScore
end

#return an array of unique country names for the time period
function subsetCountryNamesArray(stYr,endYr)
    #get the dictionary of the year and names
    countryNamesYrDict = subsetHeaderNames(stYr,endYr)
    totalNames = []
    for yearTmp in keys(countryNamesYrDict)
        totalNamesTmp = countryNamesYrDict[yearTmp]
        append!(totalNames,totalNamesTmp)
    end
    return sort(unique(totalNames))

end

#search the directory for the countries in the specific years listed; return dictionary of the year to names of countries
function subsetHeaderNames(stYr,endYr)
    dirFiles = readdir("./dataTables/")
    namesDict = Dict()

    for dF in dirFiles
        yrTmp = parse(Int,((split(dF,"."))[1]))
        if(stYr <= yrTmp <= endYr)
            fileTmp = open(string("./dataTables/",dF))
            linesTmp = readlines(fileTmp)#read each file lines

            origColNames = split(linesTmp[1],r",|\n",keepempty=false)#get line1 into components
            splice!(origColNames,1)#get rid of topleft corner
            rowNAMES = []#temp store
            for ii=2:length(linesTmp)#get the names on the rows
                rowNAMEStmp = split(linesTmp[ii],r",|\n",keepempty=false)
                append!(rowNAMES, [rowNAMEStmp[1]])
            end
            totalNamesYr = sort(unique(append!(rowNAMES,origColNames)))#for yrTmp names
            namesDict[yrTmp] = totalNamesYr#
        end
    end
    return namesDict
end

#search the directory for the complete set of countries having appeared
function totalHeaderNames()
    dirFiles = readdir("./dataTables/")
    namesDict = Dict()

    for dF in dirFiles
        yrTmp = parse(Int,((split(dF,"."))[1]))
        fileTmp = open(string("./dataTables/",dF))
      	linesTmp = readlines(fileTmp)#read each file lines
        origColNames = split(linesTmp[1],r",|\n",keepempty=false)#get line1 into components
        splice!(origColNames,1)#get rid of topleft corner
        rowNAMES = []#temp store
        for ii=2:length(linesTmp)#get the names on the rows
            rowNAMEStmp = split(linesTmp[ii],r",|\n",keepempty=false)
            append!(rowNAMES, [rowNAMEStmp[1]])
        end
        totalNamesYr = sort(unique(append!(rowNAMES,origColNames)))#for yrTmp names
        namesDict[yrTmp] = totalNamesYr#

    end
    return namesDict
end

max_stars_repo_name: mantzaris/eurovisionBias

max_issues_repo_name: mantzaris/eurovisionBias

max_forks_repo_name: mantzaris/eurovisionBias


text:
# Interface with DynamicPPL

function log_posterior_func(mdl, paramnames)
    vi = DynamicPPL.VarInfo(mdl)
    vns = DynamicPPL.VarName.(paramnames)
    dist_vec = [DynamicPPL.getdist(vi, vn) for vn in vns]
    function log_posterior(paramvals)
        for i in eachindex(paramvals)
            DynamicPPL.setindex!(vi, paramvals[i], vns[i])
        end
        return DynamicPPL.logjoint(mdl, vi)
    end
    return log_posterior, dist_vec
end

function extract_samples(chn, mdl)
    vi = DynamicPPL.VarInfo(mdl)
    pnames = DynamicPPL.syms(vi)
    chn_params = get_params(chn)
    nd = length(pnames)
    ns = size(chn,1) * size(chn,3)
    samples = Matrix(undef, nd, ns)
    for i in 1:nd
        samples[i,:] = extract(chn_params[pnames[i]])
    end
    return samples, pnames
end

function extract(V::AbstractArray) # Univariate case
    return V[:]
end

function extract(V::Tuple) # Multivariate case
    S = hcat([vcat(vi...) for vi in V]...)
    ns, nd = size(S)
    return [S[i,:] for i = 1:ns]
end


text:
"""
`Basics.analyze()`  ... analyzes various data structures with regard to requested details.

  + `("level decomposition: % of NR configurations", level::Level)`
    ... to analyze and list the non-relativistic configurations with a weight larger than 5%. A list of NR configurations with their
        corresponding weights are printed, but nothing is returned.
"""
function Basics.analyze(sa::String, level::Level)
    !(sa == "level decomposition: % of NR configurations")   &&   error("Unsupported keystring = $sa")
    !(level.hasBasis)                                        &&   error("Levels without a basis cannot be analyzed.")

    confList   = generate("configuration list: NR, from basis", level.basis)
    percentage = zeros( length(confList) )

    # Now generate for each CSF in this basis a Configuration, compare with confList and add contribution to the corresponding percentage
    for  k = 1:length(level.basis.csfs)
        csf = level.basis.csfs[k]
        shellList = Dict{Shell,Int64}[]
        for  i = 1:length(csf.occupation)
            n = csf.subshell[i].n;    l = Basics.subshell_l( csf.subshell[i] );    sh = Shell(n, l)
            if     haskey(shellList, sh)    shellList(sh) = shellList(sh) + csf.occupation[i]
            else   shellList = merge( shellList, Dict(sh => csf.occupation[i]) )
            end
        end
        confNew = Configuration( shellList, NoElectrons )
        # Compare and add the weight to the right configuration
        for  i = 1:length(confList)
            if   confList[i] == confNew    percentage[i] = percentage[i] + level.mc[k]^2;    break    end
        end
    end

    # Check that all weights sum up properly
    wa = sum( percentage );   abs(1. - wa) > 1.0e-3  &&    error("Total percentage $wa must add to 1.")
    # Sort the percentage and print both, configurations and percentage
    wb = sortperm( percentage )
    for i in wb
        println("Configuration: $(confList[i])  ... with  $(percentage[i]) %")
    end

    return( nothing )
end

"""
  + `("level decomposition: % of jj-coupled CSF", level::Level, N::Int64)`
    ... to anaylze and list (up to) N relativistic CSF, together with their weight |c_n|^2 in the expansion of the given level.
        A list of CSF and their corresponding weights are printed, but nothing is returned otherwise.  **Not yet implemented !**
"""
function Basics.analyze(sa::String, level::Level, N::Int64)
    !(sa == "level decomposition: % of jj-coupled CSF")   &&   error("Unsupported keystring = $sa")
    error("Not yet implemented !")

    return( nothing )
end

"""
`Basics.analyzeConvergence()`
    ... analyzes the convergence for two (subsequently generated) orbitals or related quantities; a value::Float64 is
        returned if not indicated otherwise.

  + `(blocka::JAC.Eigen, blockb::JAC.Eigen)`
    ... to analyze the convergence for two blocks of eigenvalues und eigenvectors; the value > 0 indicates the relative maximum
        difference abs( a-b/a+b ) between two eigenvalues. The value is set arbitrarely set to 1.0e+2 if the number of eigenvalues
        differ in the two blocks.
"""
function Basics.analyzeConvergence(blocka::Basics.Eigen, blockb::Basics.Eigen)
    va = blocka.values;   vb = blockb.values
    if       length(va) !=  length(vb)    wa = 1.0e+2
    else     wa = 0.;       for  i = 1:length(va)   wa = max(wa, abs( (va[i]-vb[i]) / (va[i]+vb[i]) ))   end
    end
    return( wa )
end

"""
  + `(a::JAC.Radial.Orbital, b::JAC.Radial.Orbital)`
    ... to analyze the convergence for two radial orbitals; the value > 0 indicates the relative maximum difference abs( a-b/a+b )
        between the energies or (large or small) components of the orbitals. The function assumes that the orbitals are represented
        on the same grid and that the components are zero if only one of the orbitals is defined beyond some maximum grid point.
"""
function Basics.analyzeConvergence(a::JAC.Radial.Orbital, b::JAC.Radial.Orbital)
    if       a.subshell !=  b.subshell                                             error("stop a")
    elseif   a.useStandardGrid !=  b.useStandardGrid   ||  !(a.useStandardGrid)    error("stop b")
    else     wa = 0.;   wa = max(wa, abs( (a.energy-b.energy) / (a.energy+b.energy) ));   nx = min(length(a.P), length(b.P))
             #=== for  i = 1:nx
                 wp = abs(a.P[i])+abs(b.P[i]);    if  wp > 1.0e-10   wa = max(wa, abs( (a.P[i]-b.P[i]) / wp ))   end
                 ##x if  abs( (a.P[i]-b.P[i]) / wp ) > 1.01    Defaults.warn(AddWarning, "analyzeConv: P, i = $i Pa=$(a.P[i]) Pb=$(b.P[i])")   end
                 wp = abs(a.Q[i])+abs(b.Q[i]);    if  wp > 1.0e-10   wa = max(wa, abs( (a.Q[i]-b.Q[i]) / wp ))   end
                 ##x if  abs( (a.Q[i]-b.Q[i]) / wp ) > 1.01    Defaults.warn(AddWarning, "analyzeConv: Q, i = $i Qa=$(a.Q[i]) Qb=$(b.Q[i]) ")   end
             end ===#
    end
    return( wa )
end


text:
MAPPING_BUSNUMBER2INDEX = Dict{Int64, Int64}()

"""
Takes a dictionary parsed by PowerModels and returns a PowerSystems
dictionary.  Currently Supports MATPOWER and PSSE data files parsed by
PowerModels
"""
function pm2ps_dict(data::Dict{String,Any}; kwargs...)
    if length(data["bus"]) < 1
        throw(DataFormatError("There are no buses in this file."))
    end
    ps_dict = Dict{String,Any}()
    ps_dict["name"] = data["name"]
    ps_dict["baseMVA"] = data["baseMVA"]
    ps_dict["source_type"] = data["source_type"]
    @info "Reading bus data"
    Buses = read_bus(data)
    if !isa(Buses,Nothing)
         ps_dict["bus"] = Buses
    else
        @error "No bus data found" # TODO : need for a model without a bus
    end
    @info "Reading load data"
    Loads= read_loads(data,ps_dict["bus"])
    LoadZones= read_loadzones(data,ps_dict["bus"])
    @info "Reading generator data"
    Generators= read_gen(data, ps_dict["bus"]; kwargs...)
    @info "Reading branch data"
    Branches= read_branch(data,ps_dict["bus"])
    Shunts = read_shunt(data,ps_dict["bus"])
    DCLines= read_dcline(data,ps_dict["bus"])

    ps_dict["load"] = Loads
    if !isa(LoadZones,Nothing)
        ps_dict["loadzone"] = LoadZones
    else
        @info "There are no Load Zones data in this file"
    end
    if !isa(Generators,Nothing)
        ps_dict["gen"] = Generators
    else
        @error "There are no Generators in this file"
    end
    if !isa(Branches,Nothing)
        ps_dict["branch"] = Branches
    else
        @info "There is no Branch data in this file"
    end
    if !isa(Shunts,Nothing)
        ps_dict["shunt"] = Shunts
    else
        @info "There is no shunt data in this file"
    end
    if !isa(DCLines,Nothing)
        ps_dict["dcline"] = DCLines
    else
        @info "There is no DClines data in this file"
    end

    return ps_dict
end

"""
Creates a PowerSystems.Bus from a PowerSystems bus dictionary
"""
function make_bus(bus_dict::Dict{String,Any})
    bus = Bus(bus_dict["number"],
                     bus_dict["name"],
                     bus_dict["bustype"],
                     bus_dict["angle"],
                     bus_dict["voltage"],
                     bus_dict["voltagelimits"],
                     bus_dict["basevoltage"]
                     )
     return bus
 end

"""
Finds the bus dictionary where a Generator/Load is located or the from & to bus
for a line/transformer
"""
function find_bus(Buses::Dict{Int64,Any},device_dict::Dict{String,Any})
    if haskey(device_dict, "t_bus")
        if haskey(device_dict, "f_bus")
            t_bus = Buses[MAPPING_BUSNUMBER2INDEX[Int(device_dict["t_bus"])]]
            f_bus = Buses[MAPPING_BUSNUMBER2INDEX[Int(device_dict["f_bus"])]]
            value =(f_bus,t_bus)
        end
    elseif haskey(device_dict, "gen_bus")
        bus = Buses[MAPPING_BUSNUMBER2INDEX[Int(device_dict["gen_bus"])]]
        value =bus
    elseif haskey(device_dict, "load_bus")
        bus = Buses[MAPPING_BUSNUMBER2INDEX[Int(device_dict["load_bus"])]]
        value =bus
    elseif haskey(device_dict,"shunt_bus")
        bus = Buses[MAPPING_BUSNUMBER2INDEX[Int(device_dict["shunt_bus"])]]
        value =bus
    else
        throw(DataFormatError("Provided Dict missing key/s:  gen_bus or f_bus/t_bus or load_bus"))
    end
    return value
end

function make_bus(bus_name, d, bus_types)
    bus = Dict{String,Any}("name" => bus_name ,
                            "number" => MAPPING_BUSNUMBER2INDEX[d["bus_i"]],
                            "bustype" => bus_types[d["bus_type"]],
                            "angle" => 0, # NOTE: angle 0, tuple(min, max)
                            "voltage" => d["vm"],
                            "voltagelimits" => (min=d["vmin"], max=d["vmax"]),
                            "basevoltage" => d["base_kv"]
                            )
    return bus
end

function read_bus(data)
    Buses = Dict{Int64,Any}()
    bus_types = ["PV", "PQ", "SF","isolated"]
    data = sort(collect(data["bus"]), by = x->parse(Int64,x[1]))
    for (i, (d_key, d)) in enumerate(data)
        # d id the data dict for each bus
        # d_key is bus key
        haskey(d,"bus_name") ? bus_name = d["bus_name"] : bus_name = string(d["bus_i"])
        bus_number = Int(d["bus_i"])
        MAPPING_BUSNUMBER2INDEX[bus_number] = i
        Buses[MAPPING_BUSNUMBER2INDEX[bus_number]] =  make_bus(bus_name, d, bus_types)
    end
    return Buses
end

function make_load(d,bus)
    load =Dict{String,Any}("name" => bus["name"],
                            "available" => true,
                            "bus" => make_bus(bus),
                            "maxactivepower" => d["pd"],
                            "maxreactivepower" => d["qd"],
                            "scalingfactor" => TimeSeries.TimeArray(collect(Dates.DateTime(Dates.today()):Dates.Hour(1):Dates.DateTime(Dates.today()+Dates.Day(1))), ones(25))
                            )
    return load
end

function read_loads(data,Buses)
    if haskey(data,"load")
        Loads = Dict{String,Any}() # Using least constrained Load
        for d_key in keys(data["load"])
            d = data["load"][d_key]
            if d["pd"] != 0.0
                # NOTE: access nodes using index i in case numbering of original data not sequential/consistent
                bus = find_bus(Buses,d)
                Loads[string(d["index"])] = make_load(d,bus)
            end
        end
        return Loads
    else
        @error "There are no loads in this file"
    end
end

function make_loadzones(d_key,d,bus_l,activepower, reactivepower)
    loadzone = Dict{String,Any}("number" => d["index"],
                                "name" => d_key ,
                                "buses" => bus_l,
                                "maxactivepower" => sum(activepower),
                                "maxreactivepower" => sum(reactivepower)
                                )
    return loadzone
end

function read_loadzones(data,Buses)
    if haskey(data,"areas")
        LoadZones = Dict{Int64,Any}()
        for (d_key,d) in data["areas"]
            b_array  = [MAPPING_BUSNUMBER2INDEX[b["bus_i"]] for (b_key, b) in data["bus"] if b["area"] == d["index"] ]
            bus_l = [make_bus(Buses[Int(b_key)]) for b_key in b_array]
            activepower  = [ l["pd"] for (l_key, l) in data["load"] if MAPPING_BUSNUMBER2INDEX[l["load_bus"]] in b_array ] #TODO: Fast Implementations
            reactivepower  = [ l["qd"] for (l_key, l) in data["load"] if MAPPING_BUSNUMBER2INDEX[l["load_bus"]] in b_array]
            LoadZones[d["index"]] = make_loadzones(d_key,d,bus_l,activepower, reactivepower)
        end
        return LoadZones
    else
        return nothing
    end
end

function make_hydro_gen(gen_name, d, bus)
    ramp_agc = get(d, "ramp_agc", get(d, "ramp_10", get(d, "ramp_30", d["pmax"])))
    hydro =  Dict{String,Any}("name" => gen_name,
                            "available" => d["gen_status"], # change from staus to available
                            "bus" => make_bus(bus),
                            "tech" => Dict{String,Any}( "installedcapacity" => float(d["pmax"]),
                                                        "activepower" => d["pg"],
                                                        "activepowerlimits" => (min=d["pmin"], max=d["pmax"]),
                                                        "reactivepower" => d["qg"],
                                                        "reactivepowerlimits" => (min=d["qmin"], max=d["qmax"]),
                                                        "ramplimits" => (up=ramp_agc/d["mbase"],down=ramp_agc/d["mbase"]),
                                                        "timelimits" => nothing),
                            "econ" => Dict{String,Any}("curtailcost" => 0.0,
                                                        "interruptioncost" => nothing),
                            "scalingfactor" => TimeSeries.TimeArray(collect(Dates.DateTime(Dates.today()):Dates.Hour(1):Dates.DateTime(Dates.today()+Dates.Day(1))), ones(25))
                            )
    return hydro
end

function make_ren_gen(gen_name, d, bus)
    gen_re = Dict{String,Any}("name" => gen_name,
                    "available" => d["gen_status"], # change from staus to available
                    "bus" => make_bus(bus),
                    "tech" => Dict{String,Any}("installedcapacity" => float(d["pmax"]),
                                                "reactivepowerlimits" => (min=d["pmin"], max=d["pmax"]),
                                                "powerfactor" => 1),
                    "econ" => Dict{String,Any}("curtailcost" => 0.0,
                                                "interruptioncost" => nothing),
                    "scalingfactor" => TimeSeries.TimeArray(collect(Dates.DateTime(Dates.today()):Dates.Hour(1):Dates.DateTime(Dates.today()+Dates.Day(1))), ones(25))
                    )
    return gen_re
end

function make_thermal_gen(gen_name, d, bus)
    if haskey(d,"model")
        model = GeneratorCostModel(d["model"])
        if model == PIECEWISE_LINEAR::GeneratorCostModel
            cost_component = d["cost"]
            power_p = [i for (ix,i) in enumerate(cost_component) if isodd(ix)]
            cost_p =  [i for (ix,i) in enumerate(cost_component) if iseven(ix)]./power_p
            cost = [(p,c) for (p,c) in zip(cost_p,power_p)]
            fixedcost = cost[1][2]
        elseif model == POLYNOMIAL::GeneratorCostModel
            if d["ncost"] == 0
                cost = x-> 0
            elseif d["ncost"] == 1
                cost = x-> d["cost"][1]
            elseif d["ncost"] == 2
                cost = x-> d["cost"][1]*x + d["cost"][2]
            elseif d["ncost"] == 3
                cost = x-> d["cost"][1]*x^2 + d["cost"][2]*x + d["cost"][3]
            elseif d["ncost"] == 4
                cost = x-> d["cost"][1]*x^3 + d["cost"][2]*x^2 + d["cost"][3]*x + d["cost"][4]
            else
                throw(DataFormatError("invalid value for ncost: $(d["ncost"])"))
            end

            # TODO: Reviewers:  Is this correct?
            fixedcost = cost(d["pmin"])
        end
        startupcost = d["startup"]
        shutdncost = d["shutdown"]
    else
        @warn "Generator cost data not included for Generator: $gen_name"
        tmpcost = EconThermal()
        cost = tmpcost.variablecost
        fixedcost = tmpcost.fixedcost
        startupcost = tmpcost.startupcost
        shutdncost = tmpcost.shutdncost
    end

    # TODO GitHub #148: ramp_agc isn't always present. This value may not be correct.
    ramp_agc = get(d, "ramp_agc", get(d, "ramp_10", get(d, "ramp_30", d["pmax"])))
    thermal_gen = Dict{String,Any}("name" => gen_name,
                                    "available" => d["gen_status"],
                                    "bus" => make_bus(bus),
                                    "tech" => Dict{String,Any}("activepower" => d["pg"],
                                                                "activepowerlimits" => (min=d["pmin"], max=d["pmax"]),
                                                                "reactivepower" => d["qg"],
                                                                "reactivepowerlimits" => (min=d["qmin"], max=d["qmax"]),
                                                                "ramplimits" => (up=ramp_agc/d["mbase"], down=ramp_agc/d["mbase"]),
                                                                "timelimits" => nothing),
                                    "econ" => Dict{String,Any}("capacity" => d["pmax"],
                                                                "variablecost" => cost,
                                                                "fixedcost" => fixedcost,
                                                                "startupcost" => startupcost,
                                                                "shutdncost" => shutdncost,
                                                                "annualcapacityfactor" => nothing)
                                    )
    return thermal_gen
end

"""
Transfer generators to ps_dict according to their classification
"""
function read_gen(data, Buses; kwargs...)

    if :genmap_file in keys(kwargs)
        genmap_file = kwargs[:genmap_file]
    else # use default generator mapping config file
        genmap_file = joinpath(dirname(dirname(pathof(PowerSystems))),
                               "src/parsers/generator_mapping.yaml")
    end
    genmap_dict = open(genmap_file) do file
        YAML.load(file)
    end

    generators = Dict{String,Any}()
    generators["Thermal"] = Dict{String,Any}()
    generators["Hydro"] = Dict{String,Any}()
    generators["Renewable"] = Dict{String,Any}()
    generators["Renewable"]["PV"]= Dict{String,Any}()
    generators["Renewable"]["RTPV"]= Dict{String,Any}()
    generators["Renewable"]["WIND"]= Dict{String,Any}()
    generators["Storage"] = Dict{String,Any}() # not currently used? JJS 3/13/19

    if !haskey(data, "gen")
        return nothing
    end

    fuel = []
    gen_name =[]
    type_gen =[]
    for (d_key,d) in data["gen"]

        fuel = uppercase(get(d, "fuel", "generic"))
        type_gen = uppercase(get(d, "type", "generic"))
        if haskey(d, "name")
            gen_name = d["name"]
        elseif haskey(d, "source_id")
            gen_name = strip(string(d["source_id"][1])*"-"*d["source_id"][2])
        else
            gen_name = d_key
        end

        bus = find_bus(Buses, d)

        assigned = false
        for (rkey, rval) in generators["Renewable"]
            fuelkeys = genmap_dict[rkey]["fuel"]
            typekeys = genmap_dict[rkey]["type"]
            if fuel in fuelkeys && type_gen in typekeys
                generators["Renewable"][rkey][gen_name] = make_ren_gen(gen_name, d, bus)
                assigned = true
                break
            end
        end
        if !assigned
            fuelkeys = genmap_dict["Hydro"]["fuel"]
            typekeys = genmap_dict["Hydro"]["type"]
            if fuel in fuelkeys && type_gen in typekeys
                generators["Hydro"][gen_name] = make_hydro_gen(gen_name, d, bus)
                assigned = true
            end
        end
        if !assigned
            # default to Thermal type if not already assigned
            generators["Thermal"][gen_name] = make_thermal_gen(gen_name, d, bus)
        end

    end # for (d_key,d) in data["gen"]

    return generators
end

function make_transformer(b_name, d, bus_f, bus_t)
    trans = Dict{String,Any}("name" => b_name,
                            "available" => Bool(d["br_status"]),
                            "connectionpoints" => (from=make_bus(bus_f),to=make_bus(bus_t)),
                            "r" => d["br_r"],
                            "x" => d["br_x"],
                            "primaryshunt" => d["b_fr"] ,  #TODO: which b ??
                            "tap" => d["tap"],
                            "rate" => d["rate_a"],
                            "α" => d["shift"]
                            )
    return trans
end

function make_lines(b_name, d, bus_f, bus_t)
    line = Dict{String,Any}("name" => b_name,
                            "available" => Bool(d["br_status"]),
                            "connectionpoints" => (from=make_bus(bus_f),to=make_bus(bus_t)),
                            "r" => d["br_r"],
                            "x" => d["br_x"],
                            "b" => (from=d["b_fr"],to=d["b_to"]),
                            "rate" =>  d["rate_a"],
                            "anglelimits" => (min=rad2deg(d["angmin"]),max =rad2deg(d["angmax"]))
                            )
    return line
end

function read_branch(data,Buses)
    Branches = Dict{String,Any}()
    Branches["Transformers"] = Dict{String,Any}()
    Branches["Lines"] = Dict{String,Any}()
    if haskey(data,"branch")
        b_name = []
        for (d_key,d) in data["branch"]
            haskey(d,"name") ? b_name = d["name"] : b_name = d_key
            (bus_f,bus_t) = find_bus(Buses,d)
            if d["transformer"]  #TODO : 3W Transformer
                Branches["Transformers"][b_name] = make_transformer(b_name, d, bus_f, bus_t)
            else
                Branches["Lines"][b_name] = make_lines(b_name, d, bus_f, bus_t)
            end
        end
        return Branches
    else
        return nothing
    end
end

function make_dcline(l_name, d, bus_f, bus_t)
    dcline = Dict{String,Any}("name" => l_name,
                            "available" =>d["br_status"] ,
                            "connectionpoints" => (from = make_bus(bus_f),to = make_bus(bus_t) ) ,
                            "activepowerlimits_from" => (min= d["pminf"] , max = d["pmaxf"]) ,
                            "activepowerlimits_to" => (min= d["pmint"] , max =d["pmaxt"] ) ,
                            "reactivepowerlimits_from" =>  (min= d["qminf"], max =d["qmaxf"] ),
                            "reactivepowerlimits_to" =>  (min=d["qmint"] , max =d["qmaxt"] ),
                            "loss" =>  (l0=d["loss0"] , l1 =d["loss1"] )
                            )
    return dcline
end

function read_dcline(data,Buses)
    DCLines = Dict{String,Any}()
    if haskey(data,"dcline")
        for (d_key,d) in data["dcline"]
            haskey(d,"name") ? l_name =d["name"] : l_name = d_key
            (bus_f,bus_t) = find_bus(Buses,d)
            DCLines[l_name] = make_dcline(l_name, d, bus_f, bus_t)
        end
        return DCLines
    else
        return nothing
    end
end

function make_shunt(s_name, d, bus)
    shunt = Dict{String,Any}("name" => s_name,
                            "available" => d["status"],
                            "bus" => make_bus(bus),
                            "Y" => (-d["gs"] + d["bs"]im)
                            )
    return shunt
end

function read_shunt(data,Buses)
    Shunts = Dict{String,Any}()
    if haskey(data,"shunt")
        s_name =[]
        for (d_key,d) in data["shunt"]
            haskey(d,"name") ? s_name =d["name"] : s_name = d_key
            bus = find_bus(Buses,d)
            Shunts[s_name] = make_shunt(s_name, d, bus)
        end
        return Shunts
    else
        return nothing
    end
end


text:
using Documenter
using PooksoftAssetModelingKit

makedocs(sitename="PooksoftAssetModelingKit.jl",

    pages = [
        "index.md",
        "Quick Start" => "quick.md",
        "Data" => [
            "data.md"
        ],

        "Stochastic models" => [
            "monte-carlo-gbm.md",
            "sde.md"
        ],

        "Lattice models" => [
            "binomial.md",
            "ternary.md"
        ]
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
# deploydocs(
#     repo = "github.com/Pooksoft/PooksoftAssetModelingKit.jl",
#     devurl = "stable",
# )


text:
is_constraint(::typeof(flow_constraints)) = true

# Branch flow constraints
function _flow_constraints(polar::PolarForm, cons, vmag, vang)
    nlines = PS.get(polar.network, PS.NumberOfLines())
    ev = branch_flow_kernel!(polar.device)(
        cons, vmag, vang,
        polar.topology.yff_re, polar.topology.yft_re, polar.topology.ytf_re, polar.topology.ytt_re,
        polar.topology.yff_im, polar.topology.yft_im, polar.topology.ytf_im, polar.topology.ytt_im,
        polar.topology.f_buses, polar.topology.t_buses, nlines,
        ndrange=nlines,
        dependencies=Event(polar.device)
    )
    wait(ev)
    return
end

function flow_constraints(polar::PolarForm, cons, buffer::PolarNetworkState)
    _flow_constraints(polar, cons, buffer.vmag, buffer.vang)
end

# Specialized function for AD with ForwardDiff
function flow_constraints(polar::PolarForm, cons, vm, va, pbus, qbus)
    _flow_constraints(polar, cons, vm, va)
end

function flow_constraints_grad!(polar::PolarForm, cons_grad, buffer, weights)
    nlines = PS.get(polar.network, PS.NumberOfLines())
    nbus = PS.get(polar.network, PS.NumberOfBuses())
    PT = polar.topology
    fill!(cons_grad, 0)
    adj_vmag = @view cons_grad[1:nbus]
    adj_vang = @view cons_grad[nbus+1:2*nbus]
    ∂edge_vm_fr = similar(cons_grad, nlines)
    ∂edge_vm_to = similar(cons_grad, nlines)
    ∂edge_va_fr = similar(cons_grad, nlines)
    ∂edge_va_to = similar(cons_grad, nlines)
    adj_branch_flow!(weights, buffer.vmag, adj_vmag,
            buffer.vang, adj_vang,
            ∂edge_vm_fr, ∂edge_vm_to,
            ∂edge_va_fr, ∂edge_va_to,
            PT.yff_re, PT.yft_re, PT.ytf_re, PT.ytt_re,
            PT.yff_im, PT.yft_im, PT.ytf_im, PT.ytt_im,
            PT.f_buses, PT.t_buses, nlines, polar.device
    )
    return cons_grad
end

function size_constraint(polar::PolarForm{T, IT, VT, MT}, ::typeof(flow_constraints)) where {T, IT, VT, MT}
    return 2 * PS.get(polar.network, PS.NumberOfLines())
end

function bounds(polar::PolarForm{T, IT, VT, MT}, ::typeof(flow_constraints)) where {T, IT, VT, MT}
    f_min, f_max = PS.bounds(polar.network, PS.Lines(), PS.ActivePower())
    return convert(VT, [f_min; f_min]), convert(VT, [f_max; f_max])
end

function adjoint!(
    polar::PolarForm,
    pbm::AutoDiff.TapeMemory{F, S, I},
    cons, ∂cons,
    vm, ∂vm,
    va, ∂va,
    pinj, ∂pinj,
) where {F<:typeof(flow_constraints), S, I}
    nlines = PS.get(polar.network, PS.NumberOfLines())
    nbus = PS.get(polar.network, PS.NumberOfBuses())
    top = polar.topology

    fill!(pbm.intermediate.∂edge_vm_fr , 0.0)
    fill!(pbm.intermediate.∂edge_vm_to , 0.0)
    fill!(pbm.intermediate.∂edge_va_fr , 0.0)
    fill!(pbm.intermediate.∂edge_va_to , 0.0)

    adj_branch_flow!(
        ∂cons,
        vm, ∂vm,
        va, ∂va,
        pbm.intermediate.∂edge_vm_fr,
        pbm.intermediate.∂edge_vm_to,
        pbm.intermediate.∂edge_va_fr,
        pbm.intermediate.∂edge_va_to,
        top.yff_re, top.yft_re, top.ytf_re, top.ytt_re,
        top.yff_im, top.yft_im, top.ytf_im, top.ytt_im,
        top.f_buses, top.t_buses, nlines, polar.device
    )
end

function matpower_jacobian(polar::PolarForm, X::Union{State,Control}, ::typeof(flow_constraints), V)
    nbus = get(polar, PS.NumberOfBuses())
    nlines = get(polar, PS.NumberOfLines())
    pf = polar.network
    ref = pf.ref ; nref = length(ref)
    pv  = pf.pv  ; npv  = length(pv)
    pq  = pf.pq  ; npq  = length(pq)
    lines = pf.lines

    dSl_dVm, dSl_dVa = PS.matpower_lineflow_power_jacobian(V, lines)

    if isa(X, State)
        j11 = dSl_dVa[:, [pv; pq]]
        j12 = dSl_dVm[:, pq]
        return [j11 j12]
    elseif isa(X, Control)
        j11 = dSl_dVm[:, [ref; pv]]
        j12 = spzeros(2 * nlines, npv)
        return [j11 j12]
    end
end


text:
# This file is a part of Julia, but is derived from
# https://github.com/google/double-conversion which has the following license
#
# Copyright 2006-2014, the V8 project authors. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
#     * Redistributions of source code must retain the above copyright
#       notice, this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above
#       copyright notice, this list of conditions and the following
#       disclaimer in the documentation and/or other materials provided
#       with the distribution.
#     * Neither the name of Google Inc. nor the names of its
#       contributors may be used to endorse or promote products derived
#       from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

function roundweed(buffer,len,rest,tk,unit,kappa)
    unit >= tk && return false, kappa
    tk - unit <= unit && return false, kappa
    tk - rest > rest && (tk - 2 * rest >= 2 * unit) && return true, kappa
    if rest > unit && (tk - (rest - unit) <= (rest - unit))
        buffer[len-1] += 1
        for i = (len-1):-1:2
            buffer[i] != 0x30 + 10 && break
            buffer[i] = 0x30
            buffer[i-1] += 1
        end
        if buffer[1] == 0x30 + 10
            buffer[1] = 0x31
            kappa += 1
        end
        return true, kappa
    end
    return false, kappa
end

function digitgen(w,buffer,requested_digits=1000)
    unit::UInt64 = 1
    one = Float(unit << -w.e, w.e)
    integrals = w.s >> -one.e
    fractionals = w.s & (one.s-1)
    divisor, kappa = bigpowten(integrals, 64 + one.e)
    len = 1
    rest = 0
    while kappa > 0
        digit = div(integrals,divisor)
        buffer[len] = 0x30 + digit
        len += 1
        requested_digits -= 1
        integrals %= divisor
        kappa -= 1
        if requested_digits == 0
            rest = (UInt64(integrals) << -one.e) + fractionals
            r, kappa = roundweed(buffer, len, rest, UInt64(divisor) << -one.e,
                    unit,kappa)
            return r, kappa, len
        end
        divisor = div(divisor,10)
    end
    while requested_digits > 0 && fractionals > unit
        fractionals *= 10
        unit *= 10
        digit = fractionals >> -one.e
        buffer[len] = 0x30 + digit
        len += 1
        requested_digits -= 1
        fractionals &= one.s - 1
        kappa -= 1
    end
    requested_digits != 0 && return false, kappa, len
    r, kappa = roundweed(buffer,len,fractionals,one.s,
                         unit,kappa)
    return r, kappa, len
end

function fastprecision(v, requested_digits, buffer = Vector{UInt8}(undef, 100))
    f = normalize(Float64(v))
    ten_mk_min_exp = kMinExp - (f.e + FloatSignificandSize)
    ten_mk_max_exp = kMaxExp - (f.e + FloatSignificandSize)
    cp = binexp_cache(ten_mk_min_exp,ten_mk_max_exp)
    scaled_w = f * cp
    r, kappa, len = digitgen(scaled_w,buffer,requested_digits)
    decimal_exponent = -cp.de + kappa
    return r, len, decimal_exponent+len-1
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
function mus(z::Bool, s::State, c::Constants, d::Data)
  cs_delta_z = cumsum(s.delta[z])
  return z ? cs_delta_z : -cs_delta_z
end

function mus(z::Bool, l::Integer, s::State, c::Constants, d::Data)
  ds_zl = sum(s.delta[z][1:l])
  return z ? ds_zl : -ds_zl
end

function mus(i::Integer, n::Integer, j::Integer,
             s::State, c::Constants, d::Data)
  k = s.lam[i][n]

  # mean of the noisy class
  mus_inj = 0.0

  # look up mean
  if k > 0
    z = s.Z[j, k]
    l = s.gam[i][n, j]
    mus_inj = mus(z, l, s, c, d)
  end

  return mus_inj
end

function update_delta!(s::State, c::Constants, d::Data)
  for z in 0:1
    for l in 1:c.L[z]
      update_delta!(Bool(z), l, s, c, d)
    end
  end
end

function update_delta!(z::Bool, l::Int, s::State, c::Constants, d::Data)
  (psi_z, tau_z, lower, upper) = params(c.delta_prior[z])

  cardinality = zeros(Int, d.I)
  g_sum = zeros(Float64, d.I)

  for i in 1:d.I
    for j in 1:d.J
      for n in 1:d.N[i]
        k = s.lam[i][n]
        if k > 0
          z_jk = s.Z[j, k]
          gam_inj = s.gam[i][n, j]
          if z_jk == z && gam_inj >= l
            cardinality[i] += 1

            dz = s.delta[z][1:gam_inj][1:end .!= l]
            g_sum[i] += (-1)^(1 - z) * s.y_imputed[i][n, j] - sum(dz)
          end
        end
      end
    end
  end

  denom = 1.0 + tau_z^2 * sum(cardinality ./ s.sig2)
  new_m = (psi_z + tau_z^2 * sum(g_sum ./ s.sig2)) / denom
  new_s = sqrt(tau_z^2 / denom)

  s.delta[z][l] = rand(TruncatedNormal(new_m, new_s, lower, upper))
end

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
using NLSolvers
using LinearAlgebra
using SparseDiffTools
using SparseArrays

# Start with a stupid example
n = 10
A = sprand(10, 10, 0.1)
A = A*A' + I
F(Fx, x) = mul!(Fx, A, x)

x0, = rand(10)
xp = copy(x0)
Fx = copy(xp)

function theta(x)
   if x[1] > 0
       return atan(x[2] / x[1]) / (2.0 * pi)
   else
       return (pi + atan(x[2] / x[1])) / (2.0 * pi)
   end
end

function F_powell!(Fx, x)
        if ( x[1]^2 + x[2]^2 == 0 )
            dtdx1 = 0;
            dtdx2 = 0;
        else
            dtdx1 = - x[2] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
            dtdx2 =   x[1] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
        end
        Fx[1] = -2000.0*(x[3]-10.0*theta(x))*dtdx1 +
            200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 );
        Fx[2] = -2000.0*(x[3]-10.0*theta(x))*dtdx2 +
            200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 );
        Fx[3] =  200.0*(x[3]-10.0*theta(x)) + 2.0*x[3];
    Fx
end

krylov_problem = ResidualKrylovProblem(F_powell!; seed=rand(3))
nlsolve!(krylov_problem, rand(3), ResidualKrylov(FixedForceTerm(0.4)))

function f_2by2!(F, x)
    F[1] = (x[1]+3)*(x[2]^3-7)+18
    F[2] = sin(x[2]*exp(x[1])-1)
    F
end

function g_2by2!(J, x)
    J[1, 1] = x[2]^3-7
    J[1, 2] = 3*x[2]^2*(x[1]+3)
    u = exp(x[1])*cos(x[2]*exp(x[1])-1)
    J[2, 1] = x[2]*u
    J[2, 2] = u
    J
end

nlsolve!(f_2by2!, rand(2), Krylov(FixedForceTerm(0.4)))

function f_sparse!(F, x)
    F[1] = (x[1]+3)*(x[2]^3-7)+18
    F[2] = sin(x[2]*exp(x[1])-1)
end

function g_sparse!(J, x)
    J[1, 1] = x[2]^3-7
    J[1, 2] = 3*x[2]^2*(x[1]+3)
    u = exp(x[1])*cos(x[2]*exp(x[1])-1)
    J[2, 1] = x[2]*u
    J[2, 2] = u
end

function helical_valley()
    tpi = 8*atan(1)
    c7 = 2.5e-1
    c8 = 5e-1

    function f!(fvec, x)
        if x[1] > 0
            temp1 = atan(x[2]/x[1])/tpi
        elseif x[1] < 0
            temp1 = atan(x[2]/x[1])/tpi + c8
        else
            temp1 = c7*sign(x[2])
        end
        temp2 = sqrt(x[1]^2+x[2]^2)
        fvec[1] = 10(x[3] - 10*temp1)
        fvec[2] = 10(temp2 - 1)
        fvec[3] = x[3]
    end
    function fj!(F, J, x)
        if x[1] > 0
            temp1 = atan(x[2]/x[1])/tpi
        elseif x[1] < 0
            temp1 = atan(x[2]/x[1])/tpi + c8
        else
            temp1 = c7*sign(x[2])
        end
        temp2 = sqrt(x[1]^2+x[2]^2)
        if !isa(F, Nothing)
            F[1] = 10*(x[3] - 10*temp1)
            F[2] = 10*(temp2 - 1)
            F[3] = x[3]
        end
        if !isa(J, Nothing)
            temp = x[1]^2 + x[2]^2
            temp1 = tpi*temp
            temp2 = sqrt(temp)
            J[1,1] = 100*x[2]/temp1
            J[1,2] = -100*x[1]/temp1
            J[1,3] = 10
            J[2,1] = 10*x[1]/temp2
            J[2,2] = 10*x[2]/temp2
            J[2,3] = 0
            J[3,1] = 0
            J[3,2] = 0
            J[3,3] = 1
            if !isa(F, nothing)
                return F, J
            else
                return J
            end
        end

    end
    f!, fj!, rand(3), [-1.; 0; 0]), [-1.; 0; 0], "Helical Valley")
end


text:
using HDF5
using Base.Test
using DSGE

path = dirname(@__FILE__)

file = "$path/../reference/solve.h5"
TTT_expected = h5read(file, "TTT")
CCC_expected = h5read(file, "CCC")
RRR_expected = h5read(file, "RRR")

m = AnSchorfheide()
TTT, RRR, CCC = solve(m)

@testset "Check state-space system matches reference" begin
    @test @test_matrix_approx_eq TTT_expected TTT
    @test @test_matrix_approx_eq RRR_expected RRR
    @test @test_matrix_approx_eq CCC_expected CCC
end


text:
## Extend functions for matrix input (col-wise) ##

const vec2vec_inplace = [
                         :lpc2c!,
                         :mc2b!,
                         :mgc2b!,
                         :b2mc!,
                         :b2c!,
                         :gnorm!,
                         :ignorm!,
                         :freqt!,
                         :frqtr!,
                         :gc2gc!,
                         :mgc2mgc!
                         ]

for f in vec2vec_inplace
    @eval begin
        function ($f)(x::AbstractMatrix, args...; kargs...)
            for i = 1:size(x, 2)
                @inbounds $f(view(x, :, i), args...; kargs...)
            end
            x
        end
    end
end

const vec2vec = [
                 :_mcep,
                 :_mgcep,
                 :lpc2c,
                 :mc2b,
                 :mgc2b,
                 :b2mc,
                 :b2c,
                 :gnorm,
                 :ignorm,
                 :c2ir,
                 :freqt,
                 :gc2gc,
                 :mgc2mgc
                 ]

for f in vec2vec
    @eval begin
        function ($f)(x::AbstractMatrix, args...; kargs...)
            r = $f(x[:, 1], args...; kargs...)
            ret = Array(eltype(x), size(r, 1), size(x, 2))
            for i = 1:length(r)
                @inbounds ret[i, 1] = r[i]
            end
            for i = 2:size(x, 2)
                @inbounds ret[:, i] = $f(view(x, :, i), args...; kargs...)
            end
            ret
        end
    end
end


text:
using ReinforcementLearningEnvironmentClassicControl, ReinforcementLearningBase,
Test

test_envinterface(CartPole())
test_envinterface(Pendulum())
test_envinterface(MountainCar())