text:
###
### Sampler states
###

struct HMCState{
    TV<:AbstractVarInfo,
    TTraj<:AHMC.AbstractTrajectory,
    THam<:AHMC.Hamiltonian,
    PhType<:AHMC.PhasePoint,
    TAdapt<:AHMC.Adaptation.AbstractAdaptor,
}
    vi::TV
    i::Int
    traj::TTraj
    hamiltonian::THam
    z::PhType
    adaptor::TAdapt
end

##########################
# Hamiltonian Transition #
##########################

struct HMCTransition{T,NT<:NamedTuple,F<:AbstractFloat}
    θ::T
    lp::F
    stat::NT
end

function HMCTransition(vi::AbstractVarInfo, t::AHMC.Transition)
    theta = tonamedtuple(vi)
    lp = getlogp(vi)
    return HMCTransition(theta, lp, t.stat)
end

function metadata(t::HMCTransition)
    return merge((lp = t.lp,), t.stat)
end

DynamicPPL.getlogp(t::HMCTransition) = t.lp

###
### Hamiltonian Monte Carlo samplers.
###

"""
    HMC(ϵ::Float64, n_leapfrog::Int)

Hamiltonian Monte Carlo sampler with static trajectory.

Arguments:

- `ϵ::Float64` : The leapfrog step size to use.
- `n_leapfrog::Int` : The number of leapfrog steps to use.

Usage:

```julia
HMC(0.05, 10)
```

Tips:

- If you are receiving gradient errors when using `HMC`, try reducing the leapfrog step size `ϵ`, e.g.

```julia
# Original step size
sample(gdemo([1.5, 2]), HMC(0.1, 10), 1000)

# Reduced step size
sample(gdemo([1.5, 2]), HMC(0.01, 10), 1000)
```
"""
struct HMC{AD, space, metricT <: AHMC.AbstractMetric} <: StaticHamiltonian{AD}
    ϵ::Float64 # leapfrog step size
    n_leapfrog::Int # leapfrog step number
end

HMC(args...; kwargs...) = HMC{ADBackend()}(args...; kwargs...)
function HMC{AD}(ϵ::Float64, n_leapfrog::Int, ::Type{metricT}, space::Tuple) where {AD, metricT <: AHMC.AbstractMetric}
    return HMC{AD, space, metricT}(ϵ, n_leapfrog)
end
function HMC{AD}(
    ϵ::Float64,
    n_leapfrog::Int,
    ::Tuple{};
    kwargs...
) where AD
    return HMC{AD}(ϵ, n_leapfrog; kwargs...)
end
function HMC{AD}(
    ϵ::Float64,
    n_leapfrog::Int,
    space::Symbol...;
    metricT=AHMC.UnitEuclideanMetric
) where AD
    return HMC{AD}(ϵ, n_leapfrog, metricT, space)
end

DynamicPPL.initialsampler(::Sampler{<:Hamiltonian}) = SampleFromUniform()

# Handle setting `nadapts` and `discard_initial`
function AbstractMCMC.sample(
    rng::AbstractRNG,
    model::AbstractModel,
    sampler::Sampler{<:AdaptiveHamiltonian},
    N::Integer;
    chain_type=MCMCChains.Chains,
    resume_from=nothing,
    progress=PROGRESS[],
    nadapts=sampler.alg.n_adapts,
    discard_adapt=true,
    discard_initial=-1,
    kwargs...
)
    if resume_from === nothing
        # If `nadapts` is `-1`, then the user called a convenience
        # constructor like `NUTS()` or `NUTS(0.65)`,
        # and we should set a default for them.
        if nadapts == -1
            _nadapts = min(1000, N ÷ 2)
        else
            _nadapts = nadapts
        end

        # If `discard_initial` is `-1`, then users did not specify the keyword argument.
        if discard_initial == -1
            _discard_initial = discard_adapt ? _nadapts : 0
        else
            _discard_initial = discard_initial
        end

        return AbstractMCMC.mcmcsample(rng, model, sampler, N;
                                       chain_type=chain_type, progress=progress,
                                       nadapts=_nadapts, discard_initial=_discard_initial,
                                       kwargs...)
    else
        return resume(resume_from, N; chain_type=chain_type, progress=progress,
                      nadapts=0, discard_adapt=false, discard_initial=0, kwargs...)
    end
end

function DynamicPPL.initialstep(
    rng::AbstractRNG,
    model::AbstractModel,
    spl::Sampler{<:Hamiltonian},
    vi::AbstractVarInfo;
    init_params=nothing,
    nadapts=0,
    kwargs...
)
    # Transform the samples to unconstrained space and compute the joint log probability.
    link!(vi, spl)
    model(rng, vi, spl)

    # Extract parameters.
    theta = vi[spl]

    # Create a Hamiltonian.
    metricT = getmetricT(spl.alg)
    metric = metricT(length(theta))
    ∂logπ∂θ = gen_∂logπ∂θ(vi, spl, model)
    logπ = gen_logπ(vi, spl, model)
    hamiltonian = AHMC.Hamiltonian(metric, logπ, ∂logπ∂θ)

    # Compute phase point z.
    z = AHMC.phasepoint(rng, theta, hamiltonian)

    # If no initial parameters are provided, resample until the log probability
    # and its gradient are finite.
    if init_params === nothing
        while !isfinite(z)
            model(rng, vi, SampleFromUniform())
            link!(vi, spl)
            theta = vi[spl]

            hamiltonian = AHMC.Hamiltonian(metric, logπ, ∂logπ∂θ)
            z = AHMC.phasepoint(rng, theta, hamiltonian)
        end
    end

    # Cache current log density.
    log_density_old = getlogp(vi)

    # Find good eps if not provided one
    if iszero(spl.alg.ϵ)
        ϵ = AHMC.find_good_stepsize(hamiltonian, theta)
        @info "Found initial step size" ϵ
    else
        ϵ = spl.alg.ϵ
    end

    # Generate a trajectory.
    traj = gen_traj(spl.alg, ϵ)

    # Create initial transition and state.
    # Already perform one step since otherwise we don't get any statistics.
    t = AHMC.step(rng, hamiltonian, traj, z)

    # Adaptation
    adaptor = AHMCAdaptor(spl.alg, hamiltonian.metric; ϵ=ϵ)
    if spl.alg isa AdaptiveHamiltonian
        hamiltonian, traj, _ =
            AHMC.adapt!(hamiltonian, traj, adaptor,
                        1, nadapts, t.z.θ, t.stat.acceptance_rate)
    end

    # Update `vi` based on acceptance
    if t.stat.is_accept
        vi[spl] = t.z.θ
        setlogp!(vi, t.stat.log_density)
    else
        vi[spl] = theta
        setlogp!(vi, log_density_old)
    end

    transition = HMCTransition(vi, t)
    state = HMCState(vi, 1, traj, hamiltonian, t.z, adaptor)

    return transition, state
end

function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model::Model,
    spl::Sampler{<:Hamiltonian},
    state::HMCState;
    nadapts=0,
    kwargs...
)
    # Get step size
    @debug "current ϵ" getstepsize(spl, state)

    # Compute transition.
    hamiltonian = state.hamiltonian
    z = state.z
    t = AHMC.step(rng, hamiltonian, state.traj, z)

    # Adaptation
    i = state.i + 1
    if spl.alg isa AdaptiveHamiltonian
        hamiltonian, traj, _ =
            AHMC.adapt!(hamiltonian, state.traj, state.adaptor,
                        i, nadapts, t.z.θ, t.stat.acceptance_rate)
    else
        traj = state.traj
    end

    # Update variables
    vi = state.vi
    if t.stat.is_accept
        vi[spl] = t.z.θ
        setlogp!(vi, t.stat.log_density)
    end

    # Compute next transition and state.
    transition = HMCTransition(vi, t)
    newstate = HMCState(vi, i, traj, hamiltonian, t.z, state.adaptor)

    return transition, newstate
end

function get_hamiltonian(model, spl, vi, state, n)
    metric = gen_metric(n, spl, state)
    ℓπ = gen_logπ(vi, spl, model)
    ∂ℓπ∂θ = gen_∂logπ∂θ(vi, spl, model)
    return AHMC.Hamiltonian(metric, ℓπ, ∂ℓπ∂θ)
end

"""
    HMCDA(n_adapts::Int, δ::Float64, λ::Float64; ϵ::Float64=0.0)

Hamiltonian Monte Carlo sampler with Dual Averaging algorithm.

Usage:

```julia
HMCDA(200, 0.65, 0.3)
```

Arguments:

- `n_adapts::Int` : Numbers of samples to use for adaptation.
- `δ::Float64` : Target acceptance rate. 65% is often recommended.
- `λ::Float64` : Target leapfrog length.
- `ϵ::Float64=0.0` : Inital step size; 0 means automatically search by Turing.

For more information, please view the following paper ([arXiv link](https://arxiv.org/abs/1111.4246)):

- Hoffman, Matthew D., and Andrew Gelman. "The No-U-turn sampler: adaptively
  setting path lengths in Hamiltonian Monte Carlo." Journal of Machine Learning
  Research 15, no. 1 (2014): 1593-1623.
"""
struct HMCDA{AD, space, metricT <: AHMC.AbstractMetric} <: AdaptiveHamiltonian{AD}
    n_adapts    ::  Int         # number of samples with adaption for ϵ
    δ           ::  Float64     # target accept rate
    λ           ::  Float64     # target leapfrog length
    ϵ           ::  Float64     # (initial) step size
end
HMCDA(args...; kwargs...) = HMCDA{ADBackend()}(args...; kwargs...)
function HMCDA{AD}(n_adapts::Int, δ::Float64, λ::Float64, ϵ::Float64, ::Type{metricT}, space::Tuple) where {AD, metricT <: AHMC.AbstractMetric}
    return HMCDA{AD, space, metricT}(n_adapts, δ, λ, ϵ)
end

function HMCDA{AD}(
    δ::Float64,
    λ::Float64;
    init_ϵ::Float64=0.0,
    metricT=AHMC.UnitEuclideanMetric
) where AD
    return HMCDA{AD}(-1, δ, λ, init_ϵ, metricT, ())
end

function HMCDA{AD}(
    n_adapts::Int,
    δ::Float64,
    λ::Float64,
    ::Tuple{};
    kwargs...
) where AD
    return HMCDA{AD}(n_adapts, δ, λ; kwargs...)
end

function HMCDA{AD}(
    n_adapts::Int,
    δ::Float64,
    λ::Float64,
    space::Symbol...;
    init_ϵ::Float64=0.0,
    metricT=AHMC.UnitEuclideanMetric
) where AD
    return HMCDA{AD}(n_adapts, δ, λ, init_ϵ, metricT, space)
end

"""
    NUTS(n_adapts::Int, δ::Float64; max_depth::Int=5, Δ_max::Float64=1000.0, ϵ::Float64=0.0)

No-U-Turn Sampler (NUTS) sampler.

Usage:

```julia
NUTS()            # Use default NUTS configuration.
NUTS(1000, 0.65)  # Use 1000 adaption steps, and target accept ratio 0.65.
```

Arguments:

- `n_adapts::Int` : The number of samples to use with adaptation.
- `δ::Float64` : Target acceptance rate for dual averaging.
- `max_depth::Int` : Maximum doubling tree depth.
- `Δ_max::Float64` : Maximum divergence during doubling tree.
- `init_ϵ::Float64` : Inital step size; 0 means automatically searching using a heuristic procedure.

"""
struct NUTS{AD,space,metricT<:AHMC.AbstractMetric} <: AdaptiveHamiltonian{AD}
    n_adapts::Int         # number of samples with adaption for ϵ
    δ::Float64        # target accept rate
    max_depth::Int         # maximum tree depth
    Δ_max::Float64
    ϵ::Float64     # (initial) step size
end

NUTS(args...; kwargs...) = NUTS{ADBackend()}(args...; kwargs...)

function NUTS{AD}(
    n_adapts::Int,
    δ::Float64,
    max_depth::Int,
    Δ_max::Float64,
    ϵ::Float64,
    ::Type{metricT},
    space::Tuple
) where {AD, metricT}
    return NUTS{AD, space, metricT}(n_adapts, δ, max_depth, Δ_max, ϵ)
end

function NUTS{AD}(
    n_adapts::Int,
    δ::Float64,
    ::Tuple{};
    kwargs...
) where AD
    NUTS{AD}(n_adapts, δ; kwargs...)
end

function NUTS{AD}(
    n_adapts::Int,
    δ::Float64,
    space::Symbol...;
    max_depth::Int=10,
    Δ_max::Float64=1000.0,
    init_ϵ::Float64=0.0,
    metricT=AHMC.DiagEuclideanMetric
) where AD
    NUTS{AD}(n_adapts, δ, max_depth, Δ_max, init_ϵ, metricT, space)
end

function NUTS{AD}(
    δ::Float64;
    max_depth::Int=10,
    Δ_max::Float64=1000.0,
    init_ϵ::Float64=0.0,
    metricT=AHMC.DiagEuclideanMetric
) where AD
    NUTS{AD}(-1, δ, max_depth, Δ_max, init_ϵ, metricT, ())
end

function NUTS{AD}(kwargs...) where AD
    NUTS{AD}(-1, 0.65; kwargs...)
end

for alg in (:HMC, :HMCDA, :NUTS)
    @eval getmetricT(::$alg{<:Any, <:Any, metricT}) where {metricT} = metricT
end

#####
##### HMC core functions
#####

getstepsize(sampler::Sampler{<:Hamiltonian}, state) = sampler.alg.ϵ
getstepsize(sampler::Sampler{<:AdaptiveHamiltonian}, state) = AHMC.getϵ(state.adaptor)

"""
    gen_∂logπ∂θ(vi, spl::Sampler, model)

Generate a function that takes a vector of reals `θ` and compute the logpdf and
gradient at `θ` for the model specified by `(vi, spl, model)`.
"""
function gen_∂logπ∂θ(vi, spl::Sampler, model)
    function ∂logπ∂θ(x)
        return gradient_logp(x, vi, model, spl)
    end
    return ∂logπ∂θ
end

"""
    gen_logπ(vi, spl::Sampler, model)

Generate a function that takes `θ` and returns logpdf at `θ` for the model specified by
`(vi, spl, model)`.
"""
function gen_logπ(vi, spl::AbstractSampler, model)
    function logπ(x)::Float64
        x_old, lj_old = vi[spl], getlogp(vi)
        vi[spl] = x
        model(vi, spl)
        lj = getlogp(vi)
        vi[spl] = x_old
        setlogp!(vi, lj_old)
        return lj
    end
    return logπ
end

gen_metric(dim::Int, spl::Sampler{<:Hamiltonian}, state) = AHMC.UnitEuclideanMetric(dim)
function gen_metric(dim::Int, spl::Sampler{<:AdaptiveHamiltonian}, state)
    return AHMC.renew(state.hamiltonian.metric, AHMC.getM⁻¹(state.adaptor.pc))
end

gen_traj(alg::HMC, ϵ) = AHMC.StaticTrajectory(AHMC.Leapfrog(ϵ), alg.n_leapfrog)
gen_traj(alg::HMCDA, ϵ) = AHMC.HMCDA(AHMC.Leapfrog(ϵ), alg.λ)
gen_traj(alg::NUTS, ϵ) = AHMC.NUTS(AHMC.Leapfrog(ϵ), alg.max_depth, alg.Δ_max)

####
#### Compiler interface, i.e. tilde operators.
####
function DynamicPPL.assume(
    rng,
    spl::Sampler{<:Hamiltonian},
    dist::Distribution,
    vn::VarName,
    vi,
)
    DynamicPPL.updategid!(vi, vn, spl)
    r = vi[vn]
    # acclogp!(vi, logpdf_with_trans(dist, r, istrans(vi, vn)))
    # r
    return r, logpdf_with_trans(dist, r, istrans(vi, vn))
end

function DynamicPPL.dot_assume(
    rng,
    spl::Sampler{<:Hamiltonian},
    dist::MultivariateDistribution,
    vns::AbstractArray{<:VarName},
    var::AbstractMatrix,
    vi,
)
    @assert length(dist) == size(var, 1)
    DynamicPPL.updategid!.(Ref(vi), vns, Ref(spl))
    r = vi[vns]
    var .= r
    return var, sum(logpdf_with_trans(dist, r, istrans(vi, vns[1])))
end
function DynamicPPL.dot_assume(
    rng,
    spl::Sampler{<:Hamiltonian},
    dists::Union{Distribution, AbstractArray{<:Distribution}},
    vns::AbstractArray{<:VarName},
    var::AbstractArray,
    vi,
)
    DynamicPPL.updategid!.(Ref(vi), vns, Ref(spl))
    r = reshape(vi[vec(vns)], size(var))
    var .= r
    return var, sum(logpdf_with_trans.(dists, r, istrans(vi, vns[1])))
end

function DynamicPPL.observe(
    spl::Sampler{<:Hamiltonian},
    d::Distribution,
    value,
    vi,
)
    return DynamicPPL.observe(SampleFromPrior(), d, value, vi)
end

function DynamicPPL.dot_observe(
    spl::Sampler{<:Hamiltonian},
    ds::Union{Distribution, AbstractArray{<:Distribution}},
    value::AbstractArray,
    vi,
)
    return DynamicPPL.dot_observe(SampleFromPrior(), ds, value, vi)
end

####
#### Default HMC stepsize and mass matrix adaptor
####

function AHMCAdaptor(alg::AdaptiveHamiltonian, metric::AHMC.AbstractMetric; ϵ=alg.ϵ)
    pc = AHMC.MassMatrixAdaptor(metric)
    da = AHMC.StepSizeAdaptor(alg.δ, ϵ)

    if iszero(alg.n_adapts)
        adaptor = AHMC.Adaptation.NoAdaptation()
    else
        if metric == AHMC.UnitEuclideanMetric
            adaptor = AHMC.NaiveHMCAdaptor(pc, da)  # there is actually no adaptation for mass matrix
        else
            adaptor = AHMC.StanHMCAdaptor(pc, da)
            AHMC.initialize!(adaptor, alg.n_adapts)
        end
    end

    return adaptor
end

AHMCAdaptor(::Hamiltonian, ::AHMC.AbstractMetric; kwargs...) = AHMC.Adaptation.NoAdaptation()

##########################
# HMC State Constructors #
##########################

function HMCState(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Hamiltonian},
    vi::AbstractVarInfo;
    kwargs...
)
    # Link everything if needed.
    !islinked(vi, spl) && link!(vi, spl)

    # Get the initial log pdf and gradient functions.
    ∂logπ∂θ = gen_∂logπ∂θ(vi, spl, model)
    logπ = gen_logπ(vi, spl, model)

    # Get the metric type.
    metricT = getmetricT(spl.alg)

    # Create a Hamiltonian.
    θ_init = Vector{Float64}(spl.state.vi[spl])
    metric = metricT(length(θ_init))
    h = AHMC.Hamiltonian(metric, logπ, ∂logπ∂θ)

    # Find good eps if not provided one
    if iszero(spl.alg.ϵ)
        ϵ = AHMC.find_good_stepsize(h, θ_init)
        @info "Found initial step size" ϵ
    else
        ϵ = spl.alg.ϵ
    end

    # Generate a trajectory.
    traj = gen_traj(spl.alg, ϵ)

    # Generate a phasepoint. Replaced during sample_init!
    h, t = AHMC.sample_init(rng, h, θ_init) # this also ensure AHMC has the same dim as θ.

    # Unlink everything.
    invlink!(vi, spl)

    return HMCState(vi, 0, 0, traj, h, AHMCAdaptor(spl.alg, metric; ϵ=ϵ), t.z)
end


text:
"""
    py_initial_state(mdp::GenerativeMergingMDP)
sample an initial state, returns the observation vector associated to the initial state.
This function is intended to be called by python pyjulia.
Note that it is using the GLOBAL RNG to sample the initial state.
"""
function py_initial_state(mdp::GenerativeMergingMDP, seed::Int64=1)
    s0 = initialstate(mdp, MersenneTwister(seed))
    return scene_to_vec(mdp, s0)
end

"""
    py_generate_s(mdp::GenerativeMergingMDP, o::Vector{Float64}, acc::Float64)
It performs one simulation step.
This funciton is intended to be called by python pyjulia.
"""
function py_generate_s(mdp::GenerativeMergingMDP, o::Vector{Float64}, acc::Float64)
    s = vec_to_scene(mdp, o)
    scene = deepcopy(s.scene)
    ego_acc = LaneFollowingAccel(acc)
    models = init_driver_models(mdp, s.ego_info.acc)
    models[EGO_ID].a = ego_acc
    # mdp.driver_models[EGO_ID].a = ego_acc
    acts = Vector{LaneFollowingAccel}(undef, length(scene))

    # call driver models
    for i=EGO_ID+1:EGO_ID+mdp.n_cars_main
        models[i].other_acc = s.ego_info.acc
    end
    get_actions!(acts, scene, mdp.env.roadway, models)

    # update scene
    tick!(scene, mdp.env.roadway, acts, mdp.dt, true)

    # clamp speed
    for (i, veh) in enumerate(scene)
        # scene[i] = clamp_speed(mdp.env, veh)
        scene[i] = wrap_around(mdp.env, scene[i])
    end
    sp = AugScene(scene, (acc=acc,))
    return scene_to_vec(mdp, sp)
end

function init_driver_models(mdp::GenerativeMergingMDP, acc::Float64)
    models = deepcopy(mdp.driver_models)
    # call driver models
    for i=EGO_ID+1:EGO_ID+mdp.n_cars_main
        models[i].other_acc = acc
    end
    return models
end

function POMDPs.isterminal(mdp::GenerativeMergingMDP, o::Vector{Float64})
    s = vec_to_scene(mdp, o)
    return isterminal(mdp, s)
end

"""
    scene_to_vec(mdp::GenerativeMergingMDP, s::AugScene)

convert a scene object to a feature vector

# components:
- ego distance to merge point (start negative, 0.0 is the merge point)
- ego speed
- ego acceleration
for each other car
    - car position on main lane
    - car speed
- collision indicator
"""
function scene_to_vec(mdp::GenerativeMergingMDP, s::AugScene)
    n_features = 2*mdp.n_cars_main + 4
    features = zeros(n_features)
    ego = get_by_id(s.scene, EGO_ID)
    features[1] = dist_to_merge(mdp.env, ego)
    features[2] = ego.state.v
    features[3] = s.ego_info.acc
    @assert s.scene[1].id == EGO_ID
    for i=2:length(s.scene)
        veh = s.scene[i]
        features[2*i] = veh.state.posF.s
        features[2*i+1] = veh.state.v
    end
    features[end] = collision_checker(s.scene, EGO_ID)
    return features
end

"""
    vec_to_scene(mdp::GenerativeMergingMDP, features::Vector{Float64})

converts a feature vector to a AugScene object
"""
function vec_to_scene(mdp::GenerativeMergingMDP, features::Vector{Float64})
    scene = Scene()
    if features[1] < 0.0
        lane_ego = merge_lane(mdp.env)
        s_ego = get_end(lane_ego) + features[1]
    else
        lane_ego = main_lane(mdp.env)
        s_ego = mdp.env.roadway[mdp.env.merge_index].s + features[1]
    end

    v_ego = features[2]
    acc_ego = features[3]
    ego = Vehicle(vehicle_state(s_ego, lane_ego, v_ego, mdp.env.roadway), VehicleDef(), EGO_ID)
    push!(scene, ego)
    for i=2:mdp.n_cars_main+1
        veh = Vehicle(vehicle_state(features[2*i], main_lane(mdp.env), features[2*i+1], mdp.env.roadway),
                      VehicleDef(), i)
        push!(scene, veh)
    end
    return AugScene(scene, (acc=acc_ego,))
end

## Rendering

"""
    AutoViz.render(mdp::GenerativeMergingMDP, features::Vector{Float64}, filename=nothing)
render a scene as a CairoSurface. If filename is passed in, write the surface to a png file.
"""
function AutoViz.render(mdp::GenerativeMergingMDP, features::Vector{Float64}, filename=nothing)
    s = vec_to_scene(mdp, features)
    c = super_render(mdp, s)
    if filename != nothing
        write_to_png(c, filename)
    end
    return c
end

"""
    make_gif(mdp::GenerativeMergingMDP, feature_list::Vector{Vector{Float64}}, filename="out.gif", fps=6)
create "filename.gif" that visualizes the list of scenes.
"""
function make_gif(mdp::GenerativeMergingMDP, feature_list::Vector{Vector{Float64}}, filename="out.gif", fps=6)
    frames = Frames(MIME("image/png"), fps=fps);
    for step in 1:length(feature_list)
        s = vec_to_scene(mdp, feature_list[step])
        c = super_render(mdp, s)
        push!(frames, c)
    end
    write(filename, frames)
end

ext: jl


text:
mutable struct IterationChannel{D} <: AbstractChannel{D}
    data::Dict{Int,D}
    cond_take::Condition
    IterationChannel(data::Dict{Int,D}) where D = new{D}(data, Condition())
end

function put!(channel::IterationChannel, t, x)
    channel.data[t] = copy(x)
    notify(channel.cond_take)
    return channel
end

function take!(channel::IterationChannel, t)
    x = fetch(channel, t)
    delete!(channel.data, t)
    return x
end

isready(channel::IterationChannel) = length(channel.data) > 1
isready(channel::IterationChannel, t) = haskey(channel.data, t)

function fetch(channel::IterationChannel, t)
    wait(channel, t)
    return channel.data[t]
end

function wait(channel::IterationChannel, t)
    while !isready(channel, t)
        wait(channel.cond_take)
    end
end

mutable struct RunningAverage{D}
    average::D
    data::Vector{D}
    buffer::Dict{Int,D}

    RunningAverage(average::D, data::Vector{D}) where D = new{D}(average, data, Dict{Int,D}())
end

function average(running_average::RunningAverage)
    return running_average.average
end

function subtract!(running_average::RunningAverage, i::Integer)
    running_average.buffer[i] = copy(running_average.data[i])
end

function add!(running_average::RunningAverage, i::Integer, π::AbstractFloat)
    running_average.average -= π * running_average.buffer[i]
    running_average.average += π * running_average.data[i]
    return running_average
end

function add!(running_average::RunningAverage{D}, i::Integer, x::D, π::AbstractFloat) where D
    running_average.average -= π * running_average.buffer[i]
    running_average.average += π * x
    running_average.data[i] = copy(x)
    return running_average
end

IteratedValue{T <: AbstractFloat} = RemoteChannel{IterationChannel{T}}
RemoteRunningAverage{D} = RemoteChannel{Channel{RunningAverage{D}}}
RemoteIterates{A <: AbstractArray} = RemoteChannel{IterationChannel{A}}


text:
using Flux,CUDA,Images,BSON

loadModelsGPU = (path;top="model-up.bson",down="model-down.bson") -> (
    filename = path*top;
    up = BSON.load(filename) ;
    model_up = up ;
    filename = path*down;
    down = BSON.load(filename) ;
    model_down = down;
    return model_up,model_down)

function calc_one(image,model_up,model_down)
    image = image #|> gpu
    y_up= model_up(image)
    y_down= model_down(image)
    y_up = y_up #|> cpu
    y_down = y_down #|> cpu
    return y_up,y_down
end

function NN_eval(path_NN, path_im)
    x = [0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,0.96,0.97,0.98,0.99,1]
    a = loadModelsGPU(path_NN)
    CNN_C5_up = a[1][:md] #|> gpu
    CNN_C5_down = a[2][:md] #|> gpu

    #---- Image Loading ----
    FF=Array{Float32}(undef,(32,32,3,1))
    #--- image loading-----
    beep = path_im
    image = load(beep)
    image = channelview(image) #CWH not WHC ie. 3X32X32
    Channels,Width,Height=size(image[:,:,:])
    for C = 1:Channels-1,W = 1:Width,H = 1:Height
        FF[W,H,C,1] = image[C,W,H] #Switching from CWH to WHCN
    end
    up,down=calc_one(FF,CNN_C5_up,CNN_C5_down)
    up = hcat(x,up)
    down = hcat(x,down)
    return up, down
end


text:
module ASCIIPlots
    using Compat

    export ASCIIPlot
    export imagesc, lineplot, scatterplot

    immutable ASCIIPlot
        s::ASCIIString
    end

    Base.show(io::IO, p::ASCIIPlot) = print(io, p.s)

    include("imagesc.jl")
    include("lineplot.jl")
    include("scatterplot.jl")
end


text:
module AsmMacro
using MLStyle
export @asm

genlabel(sym::Symbol) = string(".", string(sym)[2:end], "\${:uid}")

genarg(args, arg) = @match arg begin
    :($value[$ii]) => begin
            ii′ = ii isa Number ? ii : eval(ii)
            string(ii′, "(", genarg(args, value), ")")
        end
    ::Symbol => begin
            idx = findfirst(isequal(arg), args)
            idx === nothing && return string("%", arg)
            "\$$(idx-1)"
        end
    Expr(:macrocall, sym, ::LineNumberNode) => genlabel(sym)
    _ => error("?!! $arg")
end

# TODO add more of those
typ2llvm(::Type{Float64}) = "double"
typ2llvm(::Type{Float32}) = "float"
typ2llvm(::Type{Int32}) = "i32"
typ2llvm(::Type{Int64}) = "i64"
typ2llvm(::Type{Ptr{T}}) where {T} = typ2llvm(Int)

const DEBUG_ASM = Ref(false)

function genasm(args, xs)
    io = IOBuffer()
    argnames = Symbol[]
    typs = []
    for arg in args; @match arg begin
        :($(argname :: Symbol) :: $typ) => begin
                push!(argnames, argname)
                push!(typs, eval(typ))
            end
        a => error("invalid arg sig $a")
    end end
    println(io, "call void asm \"")
    for ex in xs; @match ex begin
        ::LineNumberNode || Expr(:line, _...) => nothing
        :($op($(opargs...))) => begin
                str_opargs = join(map(a -> genarg(argnames, a), opargs), ", ")
                println(io, op, " ", str_opargs)
            end
        Expr(:macrocall, sym, ::LineNumberNode) =>
            println(io, genlabel(sym), ":")
        ::Expr => begin
            dump(ex)
            error("unknown expr $ex")
        end
        _ => error("??? $(typeof(ex))")
    end end
    llvmtypes = map(typ2llvm, typs)
    for i = 1:length(llvmtypes)
        llvmtypes[i] = string(llvmtypes[i], " %", i-1)
    end
    constr = map(_ -> "r", llvmtypes)
    println(io, "\",\"", join(constr, ","), "\"(", join(llvmtypes, ", "), ")")
    println(io, "ret void")
    asm = String(take!(io))
    DEBUG_ASM[] && println(asm)
    Expr(:call, GlobalRef(Base, :llvmcall), asm, Cvoid, Tuple{typs...}, args...)
end

macro asm(f)
    @match f begin
        :(function $f($(args...)) $(body...) end) =>
            :(function $f($(args...))
                $(genasm(args, body))
              end) |> esc
        _ => error("invalid function form $f")
    end
end

end # module


text:
export applyMatrix, popMatrix, pushMatrix, resetMatrix
export rotate, Scale, shearX, shearY, translate
export rotateX, rotateY, rotateZ

function applyMatrix(n00, n01, n02, n03, n10, n11, n12, n13, n20, n21, n22, n23, n30, n31, n32, n33)
	m = GLfloat[n00 n01 n02 n03;
		n10 n11 n12 n13;
		n20 n21 n22 n23;
		n30 n31 n32 n33]

	GLmatState.currMatrix = m * GLmatState.currMatrix

	glUniformMatrix4fv(glGetUniformLocation(shaderBank["basicShapes"], "MVP"), 1, false, GLmatState.currMatrix)
	glUniformMatrix4fv(glGetUniformLocation(shaderBank["texturedShapes"], "MVP"), 1, false, GLmatState.currMatrix)
end

function popMatrix()
	GLmatState.currMatrix = pop!(GLmatState.matrixStack)

	glUniformMatrix4fv(glGetUniformLocation(shaderBank["basicShapes"], "MVP"), 1, false, GLmatState.currMatrix)
	glUniformMatrix4fv(glGetUniformLocation(shaderBank["texturedShapes"], "MVP"), 1, false, GLmatState.currMatrix)
end

function pushMatrix()
	push!(GLmatState.matrixStack, GLmatState.currMatrix)
end

function resetMatrix()
	GLmatState.currMatrix = GLfloat[1 0 0 0;
		0 1 0 0;
		0 0 1 0;
		0 0 0 1]

	glUniformMatrix4fv(glGetUniformLocation(shaderBank["basicShapes"], "MVP"), 1, false, GLmatState.currMatrix)
	glUniformMatrix4fv(glGetUniformLocation(shaderBank["texturedShapes"], "MVP"), 1, false, GLmatState.currMatrix)
end

function rotate(angle, x, y, z)
	m = GLfloat[cos(angle) -sin(angle) 0 0;
		sin(angle) cos(angle) 0 0;
		0 0 1 0;
		0 0 0 1]

	GLmatState.currMatrix = m * GLmatState.currMatrix

	glUniformMatrix4fv(glGetUniformLocation(shaderBank["basicShapes"], "MVP"), 1, false, GLmatState.currMatrix)
	glUniformMatrix4fv(glGetUniformLocation(shaderBank["texturedShapes"], "MVP"), 1, false, GLmatState.currMatrix)
end

function rotateX(angle)
	m = rotationmatrix_x(angle)
	m = map(Float32, convert(Array, m)) # GLfloat = float32

	GLmatState.currMatrix = m * GLmatState.currMatrix

	glUniformMatrix4fv(glGetUniformLocation(shaderBank["basicShapes"], "MVP"), 1, false, GLmatState.currMatrix)
	glUniformMatrix4fv(glGetUniformLocation(shaderBank["texturedShapes"], "MVP"), 1, false, GLmatState.currMatrix)
end

function rotateY(angle)
	m = rotationmatrix_y(angle)
	m = map(Float32, convert(Array, m)) # GLfloat = float32

	GLmatState.currMatrix = m * GLmatState.currMatrix

	glUniformMatrix4fv(glGetUniformLocation(shaderBank["basicShapes"], "MVP"), 1, false, GLmatState.currMatrix)
	glUniformMatrix4fv(glGetUniformLocation(shaderBank["texturedShapes"], "MVP"), 1, false, GLmatState.currMatrix)
end

function rotateZ(angle)
	m = rotationmatrix_z(angle)
	m = map(Float32, convert(Array, m)) # GLfloat = float32

	GLmatState.currMatrix = m * GLmatState.currMatrix

	glUniformMatrix4fv(glGetUniformLocation(shaderBank["basicShapes"], "MVP"), 1, false, GLmatState.currMatrix)
	glUniformMatrix4fv(glGetUniformLocation(shaderBank["texturedShapes"], "MVP"), 1, false, GLmatState.currMatrix)
end

function Scale(x, y, z)
	m = GLfloat[x 0 0 0;
				0 y 0 0;
				0 0 z 0;
				0 0 0 1]

	GLmatState.currMatrix = m * GLmatState.currMatrix

	glUniformMatrix4fv(glGetUniformLocation(shaderBank["basicShapes"], "MVP"), 1, false, GLmatState.currMatrix)
	glUniformMatrix4fv(glGetUniformLocation(shaderBank["texturedShapes"], "MVP"), 1, false, GLmatState.currMatrix)
end

function shearX(angle)
	m = [1 tan(angle) 0 0;
		0 1 0 0;
		0 0 1 0;
		0 0 0 1]

	GLmatState.currMatrix = m * GLmatState.currMatrix

	glUniformMatrix4fv(glGetUniformLocation(shaderBank["basicShapes"], "MVP"), 1, false, GLmatState.currMatrix)
	glUniformMatrix4fv(glGetUniformLocation(shaderBank["texturedShapes"], "MVP"), 1, false, GLmatState.currMatrix)
end

function shearY(angle)
	m = [1 0 0 0;
		tan(angle) 1 0 0;
		0 0 1 0;
		0 0 0 1]

	GLmatState.currMatrix = m * GLmatState.currMatrix

	glUniformMatrix4fv(glGetUniformLocation(shaderBank["basicShapes"], "MVP"), 1, false, GLmatState.currMatrix)
	glUniformMatrix4fv(glGetUniformLocation(shaderBank["texturedShapes"], "MVP"), 1, false, GLmatState.currMatrix)
end

function translate(x, y, z)
	m = GLfloat[1 0 0 x;
				0 1 0 y;
				0 0 1 z;
				0 0 0 1]

	GLmatState.currMatrix = m * GLmatState.currMatrix

	glUniformMatrix4fv(glGetUniformLocation(shaderBank["basicShapes"], "MVP"), 1, false, GLmatState.currMatrix)
	glUniformMatrix4fv(glGetUniformLocation(shaderBank["texturedShapes"], "MVP"), 1, false, GLmatState.currMatrix)
end

function printMatrix()
	print(GLmatState.currMatrix)
end


text:
function kepler_init(time::T,mass::T,elements::Array{T,1}) where {T <: Real}
    # Takes orbital elements of a single Keplerian; returns positions & velocities.
    # This is 3D), so 6 orbital elements specified, the code returns 3D.  For
    # Inclination = pi/2, motion is in X-Z plane; sky plane is X-Y.
    # Elements are given by: period, t0, e*cos(omega), e*sin(omega), Inclination, Omega
    period = elements[1]
    # Compute the semi-major axis in AU (or other units specified by GNEWT):
    semi = cbrt(GNEWT*mass*period^2/4/pi^2)
    # Convert to eccentricity & longitude of periastron:
    ecc2=elements[3]^2+elements[4]^2
    ecc=sqrt(ecc2)
    #omega = atan2(elements[4],elements[3])
    omega = atan(elements[4],elements[3])
    # The true anomaly at the time of transit:
    f1 = 1.5*pi-omega
    # Compute the time of periastron passage:
    sqrt1mecc2 = sqrt(1.0-ecc^2)
    tp=(elements[2]+period*sqrt1mecc2/2.0/pi*(ecc*sin(f1)/(1.0+ecc*cos(f1))
                                              #    -2.0/sqrt1mecc2*atan2(sqrt1mecc2*tan(0.5*f1),1.0+ecc)))
    -2.0/sqrt1mecc2*atan(sqrt1mecc2*tan(0.5*f1),1.0+ecc)))
    # Compute the mean anomaly
    n = 2pi/period
    m=n*(time-tp)
    # Kepler solver:
    if ecc > 0.0
        f = kepler(m,ecc)
    else
        f = m
    end
    ecosfp1 = 1.0+ecc*cos(f)
    fdot = n*ecosfp1^2/sqrt1mecc2^3
    # Compute the radial distance:
    r=semi*(1.0-ecc^2)/ecosfp1
    rdot = semi*n/sqrt1mecc2*ecc*sin(f)
    # For now assume plane-parallel:
    #inc = pi/2
    #capomega = pi
    inc = elements[5]
    capomega = elements[6]
    # Now, compute the positions
    x = zeros(eltype(elements),3)
    v = zeros(eltype(elements),3)
    if abs(capomega-pi) > 1e-15
        coscapomega = cos(capomega) ; sincapomega = sin(capomega)
    else
        coscapomega = -1.0 ; sincapomega = 0.0
    end
    cosomegapf = cos(omega+f) ; sinomegapf = sin(omega+f)
    if abs(inc-pi/2) > 1e-15
        cosinc = cos(inc) ; sininc = sin(inc)
    else
        cosinc = 0.0 ; sininc = 1.0
    end
    x[1]=r*(coscapomega*cosomegapf-sincapomega*sinomegapf*cosinc)
    x[2]=r*(sincapomega*cosomegapf+coscapomega*sinomegapf*cosinc)
    x[3]=r*sinomegapf*sininc
    rdotonr = rdot/r
    # Compute the velocities:
    rfdot = r*fdot
    v[1]=x[1]*rdotonr+rfdot*(-coscapomega*sinomegapf-sincapomega*cosomegapf*cosinc)
    v[2]=x[2]*rdotonr+rfdot*(-sincapomega*sinomegapf+coscapomega*cosomegapf*cosinc)
    v[3]=x[3]*rdotonr+rfdot*cosomegapf*sininc
    return x,v
end

function kepler_init(time::T,mass::T,elements::Array{T,1},jac_init::Array{T,2}) where {T <: Real}
    # Takes orbital elements of a single Keplerian; returns positions & velocities.
    # This is 3D), so 6 orbital elements specified, the code returns 3D.  For
    # Inclination = pi/2, motion is in X-Z plane; sky plane is X-Y.
    # Elements are given by: period, t0, e*cos(omega), e*sin(omega), Inclination, Omega
    # Returns the position & velocity of Keplerian.
    # jac_init is the derivative of (x,v,m) with respect to (elements,m).
    period = elements[1]
    n = 2pi/period
    t0 = elements[2]
    # Compute the semi-major axis in AU (or other units specified by GNEWT):
    semi = cbrt(GNEWT*mass*period^2/4/pi^2)
    dsemidp = 2third*semi/period
    dsemidm = third*semi/mass
    # Convert to eccentricity & longitude of periastron:
    ecosomega = elements[3]
    esinomega = elements[4]
    ecc=sqrt(esinomega^2+ecosomega^2)
    deccdecos = ecosomega/ecc
    deccdesin = esinomega/ecc
    #omega = atan2(esinomega,ecosomega)
    # The true anomaly at the time of transit:
    #f1 = 1.5*pi-omega
    # Compute the time of periastron passage:
    sqrt1mecc2 = sqrt(1.0-ecc^2)
    #tp=(t0+period*sqrt1mecc2/2.0/pi*(ecc*sin(f1)/(1.0+ecc*cos(f1))
    #    -2.0/sqrt1mecc2*atan2(sqrt1mecc2*tan(0.5*f1),1.0+ecc)))
    den1 = esinomega-ecosomega-ecc
    tp = (t0 - sqrt1mecc2/n*ecosomega/(1.0-esinomega)-2/n*
          #     atan2(sqrt(1.0-ecc)*(esinomega+ecosomega+ecc),sqrt(1.0+ecc)*den1))
    atan(sqrt(1.0-ecc)*(esinomega+ecosomega+ecc),sqrt(1.0+ecc)*den1))
    dtpdp = (tp-t0)/period
    fac = sqrt((1.0-ecc)/(1.0+ecc))
    den2 = 1.0/den1^2
    theta = fac*(esinomega+ecosomega+ecc)/den1
    #println("theta: ",theta," tan(atan2()): ",tan(atan2(sqrt(1.0-ecc)*(esinomega+ecosomega+ecc),sqrt(1.0+ecc)*den1)))
    dthetadecc = ((ecc+ecosomega)^2+2*(1.0-ecc^2)*esinomega-esinomega^2)/(sqrt1mecc2*(1.0+ecc))*den2
    dthetadecos = 2fac*esinomega*den2
    dthetadesin = -2fac*(ecosomega+ecc)*den2
    dtpdecc = ecc/sqrt1mecc2/n*ecosomega/(1.0-esinomega)-2/n/(1.0+theta^2)*dthetadecc
    dtpdecos = dtpdecc*deccdecos -sqrt1mecc2/n/(1.0-esinomega)-2/n/(1.0+theta^2)*dthetadecos
    dtpdesin = dtpdecc*deccdesin -sqrt1mecc2/n*ecosomega/(1.0-esinomega)^2-2/n/(1.0+theta^2)*dthetadesin
    dtpdt0 = 1.0
    # Compute the mean anomaly
    m=n*(time-tp)
    dmdp = -m/period
    dmdtp = -n
    # Kepler solver: instead of true anomaly, return eccentric anomaly:
    ekep=ekepler(m,ecc)
    cosekep = cos(ekep); sinekep = sin(ekep)
    # Compute the radial distance:
    r=semi*(1.0-ecc*cosekep)
    drdekep = semi*ecc*sinekep
    drdecc = -semi*cosekep
    denom = semi/r
    dekepdecos = sinekep*denom*deccdecos
    dekepdesin = sinekep*denom*deccdesin
    dekepdm   = denom
    inc = elements[5]
    capomega = elements[6]
    # Now, compute the positions
    coscapomega = cos(capomega) ; sincapomega = sin(capomega)
    cosomega = ecosomega/ecc; sinomega = esinomega/ecc
    cosinc = cos(inc) ; sininc = sin(inc)
    # Define rotation matrices (M&D 2.119-2.120):
    P1 = [cosomega -sinomega 0.0; sinomega cosomega 0.0; 0.0 0.0 1.0]
    P2 = [1.0 0.0 0.0; 0.0 cosinc -sininc; 0.0 sininc cosinc]
    P3 = [coscapomega -sincapomega 0.0; sincapomega coscapomega 0.0; 0.0 0.0 1.0]
    P321 = P3*P2*P1
    # Express x & v in terms of eccentric anomaly (2.121, 2.41, 2.68):
    xplane = semi*[cosekep-ecc;  sqrt1mecc2*sinekep;  0.0] # position vector in orbital plane
    vplane = [-sinekep; sqrt1mecc2*cosekep; 0.0]
    x = P321*xplane
    # Now derivatives:
    dxda = x/semi
    dxdekep = P321*semi*vplane
    P32 = P3*P2
    #dxdecc  = -x/ecc + P321*semi*[-1.0; -ecc/sqrt1mecc2*sinekep; 0.0]
    # Get rid of cancellations in the prior expression:
    dxdecc  = -P321*semi/ecc*[cosekep; sinekep/sqrt1mecc2; zero(T)]
    # These may need to be rewritten to flag ecc = 0.0 case:
    dxdecos = dxdecc*deccdecos + P32/ecc*xplane
    dxdesin = dxdecc*deccdesin + P32/ecc*[-xplane[2]; xplane[1]; 0.0]
    dxdinc  = P3*[0.0 0.0 0.0; 0.0 -sininc -cosinc; 0.0 cosinc -sininc]*P1*xplane
    dxdcom  = [-sincapomega -coscapomega 0.0; coscapomega -sincapomega 0.0; 0.0 0.0 0.0]*P2*P1*xplane
    # Compute the velocities:
    v = P321*n*semi*denom*vplane
    #decc = 1e-5
    #dvp = P321*ecc/(ecc+decc)*n*semi/(1.0-(ecc+decc)*cosekep)*[-sinekep; sqrt(1.0-(ecc+decc)^2)*cosekep; 0.0]
    #dvm = P321*ecc/(ecc-decc)*n*semi/(1.0-(ecc-decc)*cosekep)*[-sinekep; sqrt(1.0-(ecc-decc)^2)*cosekep; 0.0]
    dvda    = v/semi
    dvdp    = -v/period
    dvdekep = -v*ecc*sinekep*denom+P321*n*semi*denom*[-cosekep; -sqrt1mecc2*sinekep; 0.0]
    dvdecc  = -v/ecc + v*cosekep*denom + P321*n*semi*denom*[0.0; -ecc/sqrt1mecc2*cosekep; 0.0]
    #dvdecc_num = .5*(dvp-dvm)/decc
    #println("dvdecc: ",dvdecc," dvdecc_num: ",dvdecc_num)
    #dvdecc = dvdecc_num
    dvdecos = dvdecc*deccdecos + P32*n*semi*denom/ecc*vplane
    dvdesin = dvdecc*deccdesin + P32*n*semi*denom/ecc*[-vplane[2]; vplane[1]; 0.0]
    dvdinc  = P3*[0.0 0.0 0.0; 0.0 -sininc -cosinc; 0.0 cosinc -sininc]*P1*n*semi*denom*vplane
    dvdcom  = [-sincapomega -coscapomega 0.0; coscapomega -sincapomega 0.0; 0.0 0.0 0.0]*P2*P1*n*semi*denom*vplane
    # Now, take derivatives (11/15/2017 notes):
    # Elements are given by: period, t0, e*cos(omega), e*sin(omega), Inclination, Omega
    fill!(jac_init,0.0)
    jac_init[1:3,1] = dxda*dsemidp + dxdekep*dekepdm*(dmdp+dmdtp*dtpdp)
    #if  T != BigFloat
    #  println("position vs. period: term 1 ",dxda*dsemidp," term 2: ",dxdekep*dekepdm*dmdp," term 3: ",dxdekep*dekepdm*dmdtp*dtpdp," sum: ",jac_init[1:3,1])
    #end
    jac_init[1:3,2] = dxdekep*dekepdm*dmdtp*dtpdt0
    jac_init[1:3,3] = dxdecos + dxdekep*(dekepdm*dmdtp*dtpdecos + dekepdecos)
    #if  T != BigFloat
    #  println("position vs. ecosom : term 1 ",dxdecos," term 2: ",dxdekep*dekepdm*dmdtp*dtpdecos," term 3: ",dxdekep*dekepdm*dekepdecos," sum: ",jac_init[1:3,3])
    #end
    jac_init[1:3,4] = dxdesin + dxdekep*(dekepdm*dmdtp*dtpdesin + dekepdesin)
    #if  T != BigFloat
    #  println("position vs. esinom : term 1 ",dxdesin," term 2: ",dxdekep*dekepdm*dmdtp*dtpdesin," term 3: ",dxdekep*dekepdm*dekepdesin," sum: ",jac_init[1:3,4])
    #end
    jac_init[1:3,5] = dxdinc
    jac_init[1:3,6] = dxdcom
    jac_init[1:3,7] = dxda*dsemidm
    jac_init[4:6,1] = dvdp + dvda*dsemidp +  dvdekep*dekepdm*(dmdp+dmdtp*dtpdp)
    #if  T != BigFloat
    #  println("velocity vs. period: term 1 ",dvdp," term 2: ",dvda*dsemidp," term 3: ",dvdekep*dekepdm*dmdp," term 4: ",dvdekep*dekepdm*dmdtp*dtpdp," sum: ",jac_init[4:6,1])
    #end
    jac_init[4:6,2] = dvdekep*dekepdm*dmdtp*dtpdt0
    jac_init[4:6,3] = dvdecos + dvdekep*(dekepdm*dmdtp*dtpdecos + dekepdecos)
    #if  T != BigFloat
    #  println("velocity vs. ecosom : term 1 ",dvdecos," term 2: ",dvdekep*dekepdm*dmdtp*dtpdecos," term 3: ",dvdekep*dekepdm*dekepdecos," sum: ",jac_init[4:6,3])
    #end
    jac_init[4:6,4] = dvdesin + dvdekep*(dekepdm*dmdtp*dtpdesin + dekepdesin)
    #if  T != BigFloat
    #  println("velocity vs. esinom : term 1 ",dvdesin," term 2: ",dvdekep*dekepdm*dmdtp*dtpdesin," term 3: ",dvdekep*dekepdm*dekepdesin," sum: ",jac_init[4:6,4])
    #end
    jac_init[4:6,5] = dvdinc
    jac_init[4:6,6] = dvdcom
    jac_init[4:6,7] = dvda*dsemidm
    jac_init[7,7] = 1.0
    #return x,v,tp,dtpdecos,dtpdesin
    #return x,v,ekep,dekepdm*dmdtp*dtpdecos + dekepdecos,dekepdm*dmdtp*dtpdesin + dekepdesin
    return x,v
end


text:
export tone

"""
𝑓 = tone(name; A=440.0)
  convert tone name to frequency, equal temperament
    tone("A4") == A
    tone("C♭0") or tone("Cb0")
    tone("B♯9") or tone("B#9")
"""
function tone(name; A=440.0)
  𝑔=2^(1/12)
  wholetone=Int(name[1])-65   # 'A'=0
  octave=Int(name[end])-48 # '0'=0
  semitone = 2*wholetone
  if wholetone >= 2 #C
    semitone -= 1 # no B♯
    octave -= 1 # octave at C
  end
  if wholetone >= 5 #F
    semitone -= 1 # no E♯
  end
  if length(name) == 3
    if name[2] in ['♭', 'b']
      semitone -= 1
    elseif name[2] in ['♯', '#']
      semitone += 1
    end
  end
  return round(A/16.0 * 𝑔^semitone * 2.0^octave; digits=6)
end

max_stars_repo_name: isaacwhanson/scad-flute

max_issues_repo_name: isaacwhanson/scad-flute

max_forks_repo_name: isaacwhanson/scad-flute


text:
import HTTP
import URIs
import Sockets
import Dates
import JSON2
import MsgPack

module QuickAPI

include("util.jl")
include("errors.jl")

export
    APP,
    text, json, msgpack,
    @route, @get, @post, @put, @delete,
    serve, serve_async

const APP = HTTP.Router()

const LOCALHOST = Sockets.localhost

const CONTENT_TYPES = Dict{Symbol,String}(
  :text       => "text/plain",
  :json       => "application/json",
  :msgpack    => "application/x-msgpack",
)

#========================================
# Functions for building HTTP responses #
========================================#

function text(data::Any ; status::Int = 200, header::Array{Pair{String,String},1} = Pair{String,String}[])
    HTTP.Response(
        status,
        ["Content-Type" => CONTENT_TYPES[:text], header...];
        body = string(data) * "\n"
    )
end

function json(data::Any ; status::Int = 200, header::Array{Pair{String,String},1} = Pair{String,String}[])
    HTTP.Response(
        status,
        ["Content-Type" => CONTENT_TYPES[:json], header...];
        body = JSON2.write(data) * "\n"
    )
end

function msgpack(data::Any ; status::Int = 200, header::Array{Pair{String,String},1} = Pair{String,String}[])
    HTTP.Response(
        status,
        ["Content-Type" => CONTENT_TYPES[:msgpack], header...];
        body = MsgPack.pack(data)
    )
end

#=================================
# Macros for handling GET, POST, #
# PUT, DELETE requests           #
=================================#

macro route(method::Symbol,path::String,response)
    HTTP.@register(
        APP,
        method,
        path,
        response
    )
end

macro get(path::String, response)
    quote
        @route(
            :GET,
            $path,
            $response
        )
    end
end

macro post(path::String, response)
    quote
        @route(
            :POST,
            $path,
            $response
        )
    end
end

macro put(path::String, response)
    quote
        @route(
            :PUT,
            $path,
            $response
        )
    end
end

macro delete(path::String, response)
    quote
        @route(
            :DELETE,
            $path,
            $response
        )
    end
end

#=============================
# Functions to Serve the App #
=============================#

function serve(host=LOCALHOST, port=8081; kw...)
    @trycatch HTTP.serve(
        APP,
        host,
        port;
        kw...
    )
end

function serve_async(host=LOCALHOST, port=8081; kw...)
    @async HTTP.serve(APP,host,port,kw...)
end

end # module


text:
using PkgEval
using Test

# determine which Julia to use
const julia_spec = get(ENV, "JULIA_SPEC", string(VERSION))
@info "Testing with Julia $julia_spec"
const julia_version = PkgEval.obtain_julia(julia_spec)::VersionNumber
@info "Resolved to Julia v$julia_version"
const julia_install = PkgEval.prepare_julia(julia_version)

@testset "sandbox" begin
    mktemp() do path, io
        try
            PkgEval.run_sandboxed_julia(julia_install, `-e 'print(1337)'`; stdout=io)
            close(io)
            @test read(path, String) == "1337"
        catch
            # if we failed to spawn a container, make sure to print the reason
            flush(io)
            @error read(path, String)
            rethrow()
        end
    end

    # print versioninfo so we can verify in CI logs that the correct version is used
    PkgEval.run_sandboxed_julia(julia_install, `-e 'using InteractiveUtils; versioninfo()'`)
end

const pkgnames = ["TimerOutputs", "Crayons", "Example", "Gtk"]

@testset "low-level interface" begin
    PkgEval.prepare_registry()

    pkgs = PkgEval.read_pkgs(pkgnames)

    # timeouts
    results = PkgEval.run([Configuration(julia=julia_version)], pkgs; time_limit = 0.1)
    @test all(results.status .== :kill) && all(results.reason .== :time_limit)

    # log limit
    results = PkgEval.run([Configuration(julia=julia_version)], pkgs; log_limit = 1)
    @test all(results.status .== :kill) && all(results.reason .== :log_limit)
end

@testset "main entrypoint" begin
    results = PkgEval.run([Configuration(julia=julia_version)], pkgnames)
    if !(julia_spec == "master" || julia_spec == "nightly")
        @test all(results.status .== :ok)
        for result in eachrow(results)
            @test occursin("Testing $(result.name) tests passed", result.log)
        end
    end
end

@testset "PackageCompiler" begin
    results = PkgEval.run([Configuration(julia=julia_version, compiled=true)], ["Example"])
    if !(julia_spec == "master" || julia_spec == "nightly")
        @test all(results.status .== :ok)
        for result in eachrow(results)
            @test occursin("Testing $(result.name) tests passed", result.log)
        end
    end
end

@testset "reporting" begin
    lts = Configuration(julia=v"1.0.5")
    stable = Configuration(julia=v"1.2.0")
    results = PkgEval.run([lts, stable], ["Example"])
    PkgEval.compare(results)
end

PkgEval.purge()
rm(julia_install; recursive=true)


text:
export FisherBRCLearner

mutable struct EntropyBC{A<:NeuralNetworkApproximator}
    policy::A
    α::Float32
    lr_alpha::Float32
    target_entropy::Float32
    # Logging
    policy_loss::Float32
end

mutable struct FisherBRCLearner{
    BA<:NeuralNetworkApproximator,
    BC1<:NeuralNetworkApproximator,
    BC2<:NeuralNetworkApproximator,
    R<:AbstractRNG,
} <: AbstractLearner
    policy::BA
    behavior_policy::EntropyBC
    qnetwork1::BC1
    qnetwork2::BC2
    target_qnetwork1::BC1
    target_qnetwork2::BC2
    γ::Float32
    τ::Float32
    α::Float32
    f_reg::Float32
    reward_bonus::Float32
    batch_size::Int
    pretrain_step::Int
    update_freq::Int
    update_step::Int
    lr_alpha::Float32
    target_entropy::Float32
    rng::R
    # Logging
    qnetwork_loss::Float32
    policy_loss::Float32
end

"""
    FisherBRCLearner(;kwargs...)

See paper: [Offline reinforcement learning with fisher divergence critic regularization](https://arxiv.org/abs/2103.08050).

# Keyword arguments

- `policy`, used to get action.
- `behavior_policy::EntropyBC`, used to estimate log μ(a|s).
- `qnetwork1`, used to get Q-values.
- `qnetwork2`, used to get Q-values.
- `target_qnetwork1`, used to estimate the target Q-values.
- `target_qnetwork2`, used to estimate the target Q-values.
- `γ::Float32 = 0.99f0`, reward discount rate.
- `τ::Float32 = 0.005f0`, the speed at which the target network is updated.
- `α::Float32 = 0.0f0`, entropy term.
- `f_reg::Float32 = 1.0f0`, the weight of gradient penalty regularizer.
- `reward_bonus::Float32 = 5.0f0`, add extra value to the reward.
- `batch_size::Int = 32`
- `pretrain_step::Int = 1000`, the number of pre-training rounds.
- `update_freq::Int = 50`, the frequency of updating the `approximator`.
- `lr_alpha::Float32 = 0.003f0`, learning rate of tuning entropy.
- `action_dims::Int = 0`, the dimensionality of the action.
- `update_step::Int = 0`
- `rng = Random.GLOBAL_RNG`

`policy` is expected to output a tuple `(μ, logσ)` of mean and
log standard deviations for the desired action distributions, this
can be implemented using a `GaussianNetwork` in a `NeuralNetworkApproximator`.
"""
function FisherBRCLearner(;
    policy,
    behavior_policy,
    qnetwork1,
    qnetwork2,
    target_qnetwork1,
    target_qnetwork2,
    γ = 0.99f0,
    τ = 0.005f0,
    α = 0.0f0,
    f_reg = 1.0f0,
    reward_bonus = 5.0f0,
    batch_size = 32,
    pretrain_step = 1000,
    update_freq = 50,
    lr_alpha = 0.003f0,
    behavior_lr_alpha = 0.001f0,
    action_dims = 0,
    update_step = 0,
    rng = Random.GLOBAL_RNG,
)
    copyto!(qnetwork1, target_qnetwork1)  # force sync
    copyto!(qnetwork2, target_qnetwork2)  # force sync
    entropy_behavior_policy = EntropyBC(behavior_policy, 0.0f0, behavior_lr_alpha, Float32(-action_dims), 0.0f0)
    FisherBRCLearner(
        policy,
        entropy_behavior_policy,
        qnetwork1,
        qnetwork2,
        target_qnetwork1,
        target_qnetwork2,
        γ,
        τ,
        α,
        f_reg,
        reward_bonus,
        batch_size,
        pretrain_step,
        update_freq,
        update_step,
        lr_alpha,
        Float32(-action_dims),
        rng,
        0f0,
        0f0,
    )
end

function (l::FisherBRCLearner)(env)
    D = device(l.policy)
    s = send_to_device(D, state(env))
    s = Flux.unsqueeze(s, ndims(s) + 1)
    action = dropdims(l.policy(l.rng, s; is_sampling=true), dims=2)
end

function RLBase.update!(l::FisherBRCLearner, batch::NamedTuple{SARTS})
    if l.update_step == 0
        update_behavior_policy!(l.behavior_policy, batch)
    else
        update_learner!(l, batch)
    end
end

function update_behavior_policy!(l::EntropyBC, batch::NamedTuple{SARTS})
    D = device(l.policy)
    s, a, r, t, s′ = (send_to_device(D, batch[x]) for x in SARTS)
    # Update behavior policy with entropy
    gs = gradient(Flux.params(l.policy)) do
        log_π = l.policy.model(s, a)
        _, entropy = l.policy.model(s; is_sampling=true, is_return_log_prob=true)
        loss = mean(l.α .* entropy .- log_π)
        # Update entropy
        ignore() do
            l.policy_loss = loss
            l.α -= l.lr_alpha * mean(-entropy .- l.target_entropy)
        end
        loss
    end
    update!(l.policy, gs)
end

function update_learner!(l::FisherBRCLearner, batch::NamedTuple{SARTS})
    D = device(l.policy)
    s, a, r, t, s′ = (send_to_device(D, batch[x]) for x in SARTS)
    r .+= l.reward_bonus
    γ, τ, α = l.γ, l.τ, l.α

    a′ = l.policy(l.rng, s′; is_sampling=true)
    q′_input = vcat(s′, a′)
    target_q′ = min.(l.target_qnetwork1(q′_input), l.target_qnetwork2(q′_input))

    y = r .+ γ .* (1 .- t) .* vec(target_q′)

    # Train Q Networks
    a = reshape(a, :, l.batch_size)
    q_input = vcat(s, a)
    log_μ = l.behavior_policy.policy.model(s, a) |> vec
    a_policy = l.policy(l.rng, s; is_sampling=true)

    q_grad_1 = gradient(Flux.params(l.qnetwork1)) do
        q1 = l.qnetwork1(q_input) |> vec
        q1_grad_norm = gradient(Flux.params([a_policy])) do
            q1_reg = mean(l.qnetwork1(vcat(s, a_policy)))
        end
        reg = mean(q1_grad_norm[a_policy] .^ 2)
        loss = mse(q1 .+ log_μ, y) + l.f_reg * reg
        ignore() do
            l.qnetwork_loss = loss
        end
        loss
    end
    update!(l.qnetwork1, q_grad_1)

    q_grad_2 = gradient(Flux.params(l.qnetwork2)) do
        q2 = l.qnetwork2(q_input) |> vec
        q2_grad_norm = gradient(Flux.params([a_policy])) do
            q2_reg = mean(l.qnetwork2(vcat(s, a_policy)))
        end
        reg = mean(q2_grad_norm[a_policy] .^ 2)
        loss = mse(q2 .+ log_μ, y) + l.f_reg * reg
        ignore() do
            l.qnetwork_loss += loss
        end
        loss
    end
    update!(l.qnetwork2, q_grad_2)

    # Train Policy
    p_grad = gradient(Flux.params(l.policy)) do
        a, log_π = l.policy(l.rng, s; is_sampling=true, is_return_log_prob=true)
        q_input = vcat(s, a)
        q = min.(l.qnetwork1(q_input), l.qnetwork2(q_input)) .+ log_μ
        policy_loss = mean(α .* log_π .- q)
        ignore() do
            l.policy_loss = policy_loss
            # Tune entropy automatically
            l.α -= l.lr_alpha * mean(-log_π .- l.target_entropy)
        end
        policy_loss
    end
    update!(l.policy, p_grad)

    # polyak averaging
    for (dest, src) in zip(
        Flux.params([l.target_qnetwork1, l.target_qnetwork2]),
        Flux.params([l.qnetwork1, l.qnetwork2]),
    )
        dest .= (1 - τ) .* dest .+ τ .* src
    end
end


text:
# VERA Core Physics Benchmark Progression Problem Specifications
# Revision 4, August 29, 2014
# CASL-U-2012-0131-004
using MOCNeutronTransport

add_timestamps_to_logger()
gmsh.initialize()
gmsh.option.set_number("General.NumThreads", 0) # 0 uses system default, i.e. OMP_NUM_THREADS)
gmsh.option.set_number("Geometry.OCCParallel", 1) # use parallel OCC boolean operations
gmsh.option.set_number("General.Verbosity", 2) # 1: +errors, 2: +warnings, 3: +direct, 4: +information, 5: +status, 99: +debug

#
# Model
# ----------------------------------------------------------------------------------------------
uo2_entities = Int64[]; gap_entities = Int64[]; clad_entities = Int64[]; h2o_entities = Int64[]
r_fuel = 0.4096 # Pellet radius = 0.4096 cm (pg. 4)
r_gap = 0.418   # Inner clad radius = 0.418 cm (pg.4)
r_clad = 0.475  # Outer clad radius = 0.475 cm (pg.4)
pitch = 1.26    # Pitch = 1.26 cm (pg. 4)
half_gap = 0.04 # Inter-Assembly Half Gap  = 0.04 cm (pg. 7 or deduction from assembly/pin pitch)

# Guide tube locations (i, j) (pg. 5)
coords_gt = [ (6, 15), (9, 15), (12, 15), (4, 14), (14, 14), (3, 12), (6, 12), (9, 12), (12, 12),
              (15, 12), (3, 9), (6, 9), (12, 9),  (15, 9), (3, 6),  (6, 6),  (9, 6),  (12, 6),
              (15, 6), (4, 4), (14, 4), (6, 3), (9, 3), (12, 3)]

# Instrument tube locations (pg. 5)
coords_it = [(9,9)]

# Place UO₂ pins
for i = 1:17, j = 1:17
    if (i, j) ∈  coords_gt || (i,j) ∈  coords_it
        continue
    end
    x = half_gap + i*pitch - pitch/2
    y = half_gap + j*pitch - pitch/2
    push!(uo2_entities, gmsh.model.occ.add_disk(x, y, 0, r_fuel, r_fuel))
    push!(gap_entities, gmsh.model.occ.add_disk(x, y, 0, r_gap, r_gap))
    push!(clad_entities, gmsh.model.occ.add_disk(x, y, 0, r_clad, r_clad))
end

# Place guide tubes
r_gt_inner = 0.561 # Inner Guide Tube Radius = 0.561 cm (pg. 5)
r_gt_outer = 0.602 # Outer Guide Tube Radius = 0.602 cm (pg. 5)
for (i, j) in coords_gt
    x = half_gap + i*pitch - pitch/2
    y = half_gap + j*pitch - pitch/2
    push!(h2o_entities, gmsh.model.occ.add_disk(x, y, 0, r_gt_inner, r_gt_inner))
    push!(clad_entities, gmsh.model.occ.add_disk(x, y, 0, r_gt_outer, r_gt_outer))
end

# Place instrument tube
r_it_inner = 0.559 # Inner Instrument Tube Radius = 0.559 cm (pg. 5)
r_it_outer = 0.605 # Outer Instrument Tube Radius = 0.605 cm (pg. 5)
for (i, j) in coords_it
    x = half_gap + i*pitch - pitch/2
    y = half_gap + j*pitch - pitch/2
    push!(h2o_entities, gmsh.model.occ.add_disk(x, y, 0, r_it_inner, r_it_inner))
    push!(clad_entities, gmsh.model.occ.add_disk(x, y, 0, r_it_outer, r_it_outer))
end
gmsh.model.occ.synchronize()

# Materials
# ----------------------------------------------------------------------------------------------
# http://juliagraphics.github.io/Colors.jl/stable/namedcolors/
materials = Material[
                Material(name="UO2", color="forestgreen"),
                Material(name="Gap", color="darkorange"),
                Material(name="H2O", color="royalblue"),
                Material(name="Clad", color="slategrey")
            ]
safe_add_physical_group("Material: UO2",  [(2,i) for i in uo2_entities])
safe_add_physical_group("Material: Gap",  [(2,i) for i in gap_entities])
safe_add_physical_group("Material: H2O",  [(2,i) for i in h2o_entities])
safe_add_physical_group("Material: Clad", [(2,i) for i in clad_entities])
ents = gmsh.model.get_entities(2)
# Fragment the many disks, prioritizing UO2 > Gap > H2O > Clad to fill space
safe_fragment(ents, ents, material_hierarchy = materials)

# Overlay Grid
# ---------------------------------------------------------------------------------------------
coarse_divs = [0.0]
append!(coarse_divs, [pitch*i + 0.04 for i = 1:16])
push!(coarse_divs, 21.5) # Assembly pitch = 21.50 cm (pg. 7)
scoarse_divs = Vec{18,Float64}(coarse_divs)
coarse_grid = RectilinearGrid(scoarse_divs, scoarse_divs)
mpact_grid = MPACTGridHierarchy(coarse_grid)
# We now want H2O to fill empty space, preserving all other materials,
# so we need to swap H2O to the bottom of the materials hierarchy
materials[3], materials[4] = materials[4], materials[3]
overlay_mpact_grid_hierarchy(mpact_grid, materials)

# Mesh
# ------------------------------------------------------------------------------------------------
lc = 0.25 # cm
for mat in materials
    mat.mesh_size = lc
end
set_mesh_field_using_materials(materials)
generate_mesh(order = 1, faces = "Triangle", opt_iters = 2)
gmsh.write("2a.inp")
mesh_error = get_cad_to_mesh_error()
for i in eachindex(mesh_error)
    println(mesh_error[i])
end
# gmsh.fltk.run()
gmsh.finalize()

mesh = import_mesh("2a.inp")
statistics(mesh)
# Partition mesh according to mpact grid hierarchy and write as an xdmf file
mpt = MeshPartitionTree(mesh)
export_mesh(mpt, "2a.xdmf")

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using CMark
using Test
import Markdown
using cmark_jll: cmark_jll

@testset "CMark.jl" begin
    # Write your own tests here.
    @test CMark.cmark_version() == 1900547
    @test CMark.cmark_version_string() == "0.29.0.gfm.3"

    # simple empty Markdown document
    let md = "" |> CMark.parse_document |> CMark.markdown_withsiblings
        @test length(md) == 1
        @test md[1] isa Markdown.MD
        @test isempty(md[1])
    end

    # single string
    let md = "foo" |> CMark.parse_document |> CMark.markdown_withsiblings
        @test length(md) == 1
        @test md[1] isa Markdown.MD
        @test length(md[1]) == 1
        @test md[1][1] isa Markdown.Paragraph
        @test length(md[1][1].content) == 1
        @test md[1][1].content[1] == "foo"
    end

    # some basic inline elements
    let md = "foo **strong** _emphasis_" |> markdown
        @test md isa Markdown.MD
        @test length(md) == 1
        @test md[1] isa Markdown.Paragraph
        @test length(md[1].content) == 4
        @test md[1].content[1] == "foo "
        @test md[1].content[2] isa Markdown.Bold
        @test md[1].content[3] == " "
        @test md[1].content[4] isa Markdown.Italic
    end

    # Headings and inline code
    let md = """
        # Heading `1`
        ## Heading `2`
        ### Heading `3`
        #### Heading `4`
        ##### Heading `5`
        ###### Heading `6`

        Setext `1`
        ==========

        Setext `2`
        ----------
        """ |> markdown
        @test md isa Markdown.MD
        @test length(md) == 8
        for h in md.content
            @test h isa Markdown.Header
            @test length(h.text) == 2
            @test h.text[2] isa Markdown.Code
            level = parse(Int, h.text[2].code)
            @test h isa Markdown.Header{level}
        end
    end

    # code, blockquote and thematic break blocks
    let md = """
        ```foo bar baz
        code
        example()
        ```
        ---
        > Paragraph 1
        ---
        > Paragraph 1
        >
        > Paragraph 2
        ---
            Indented code
        """ |> CMark.parse_document |> CMark.markdown
        @test md isa Markdown.MD
        @test length(md) == 7
        @test md[1] isa Markdown.Code
        @test md[1].language == "foo bar baz"
        @test md[1].code == "code\nexample()\n"
        @test md[2] isa Markdown.HorizontalRule
        @test md[3] isa Markdown.BlockQuote
        @test length(md[3].content) == 1
        @test md[4] isa Markdown.HorizontalRule
        @test md[5] isa Markdown.BlockQuote
        @test length(md[5].content) == 2
        @test md[5].content[1] isa Markdown.Paragraph
        @test md[6] isa Markdown.HorizontalRule
        @test md[7] isa Markdown.Code
        @test md[7].language == ""
        @test md[7].code == "Indented code\n"
    end

    # lists
    let md = """
        * A
        * B
        ---
        * A

          B
        * C
        ---
        1. Foo
        2. Bar
        """ |> markdown
        @test md isa Markdown.MD
        @test length(md) == 5

        @test md[1] isa Markdown.List
        @test md[1].ordered == -1
        @test length(md[1].items) == 2
        @test md[1].items[1] isa Vector
        @test length(md[1].items[1]) == 1
        @test md[1].items[1][1] isa Markdown.Paragraph

        let list = md[3]
            @test list isa Markdown.List
            @test list.ordered == -1
            @test length(list.items) == 2
            @test list.items[1] isa Vector
            @test length(list.items[1]) == 2
            @test list.items[1][1] isa Markdown.Paragraph
            @test list.items[1][2] isa Markdown.Paragraph
            @test list.items[2] isa Vector
            @test length(list.items[2]) == 1
            @test list.items[2][1] isa Markdown.Paragraph
        end

        let list = md[5]
            @test list isa Markdown.List
            @test list.ordered == 1
            @test length(list.items) == 2
            @test list.items[1] isa Vector
            @test length(list.items[1]) == 1
            @test list.items[1][1] isa Markdown.Paragraph
            @test list.items[2] isa Vector
            @test length(list.items[2]) == 1
            @test list.items[2][1] isa Markdown.Paragraph
        end
    end

    # linebreaks
    let md = "Soft\nBreak\n\n---\n\nHard  \nBreak\n" |> markdown
        @test md isa Markdown.MD
        @test length(md) == 3
        let p = md[1]
            @test p isa Markdown.Paragraph
            @test length(p.content) == 3
            @test p.content[1] == "Soft"
            @test p.content[2] == " "
            @test p.content[3] == "Break"
        end
        let p = md[3]
            @test p isa Markdown.Paragraph
            @test length(p.content) == 3
            @test p.content[1] == "Hard"
            @test p.content[2] == "\n"
            @test p.content[3] == "Break"
        end
    end

    # links and link definitions
    let md = """
        [link text](url)

        [link text](url "title")

        <http://foo.bar.baz>

        [text][link1]

        [text][link2]

        [link1]: foo
        [link2]: bar "title"
        """ |> markdown
        @test md isa Markdown.MD
        @test length(md) == 5
        @test all(isa.(md.content, Markdown.Paragraph))
        @test all(length(p.content) == 1 for p in md.content)

        let link = md[1].content[1]
            @test link isa Markdown.Link
            @test link.text == ["link text"]
            @test link.url == "url"
        end
        let link = md[2].content[1]
            @test link isa Markdown.Link
            @test link.text == ["link text"]
            @test link.url == "url"
        end
        let link = md[3].content[1]
            @test link isa Markdown.Link
            @test link.text == ["http://foo.bar.baz"]
            @test link.url == "http://foo.bar.baz"
        end
        let link = md[4].content[1]
            @test link isa Markdown.Link
            @test link.text == ["text"]
            @test link.url == "foo"
        end
        let link = md[5].content[1]
            @test link isa Markdown.Link
            @test link.text == ["text"]
            @test link.url == "bar"
        end
    end

    # image
    let md = "![alt text](url \"title\")" |> markdown
        @test md isa Markdown.MD
        @test length(md) == 1
        @test md[1] isa Markdown.Paragraph
        let image = md[1].content[1]
            @test image isa Markdown.Image
            @test image.url == "url"
            @test image.alt == ""
        end
    end

    # HTML
    let md = """
        Inline HTML!

HTML blocks

        """ |> markdown
        @test md isa Markdown.MD
        @test length(md) == 2
        let p = md[1]
            @test p isa Markdown.Paragraph
            @test length(p.content) == 5
            @test p.content[2] isa Markdown.Code
            @test p.content[2].language == "@raw html"
            @test p.content[2].code == ""
            @test p.content[4] isa Markdown.Code
            @test p.content[4].language == "@raw html"
            @test p.content[4].code == ""
        end
        let code = md[2]
            @test code isa Markdown.Code
            @test code.language == "@raw html"
            @test code.code == "
HTML blocks
\n"
        end
    end

    # GFM extensions
    let md = """
        ~foo bar~

        | foo | bar       |
        | --- | --------- |
        | baz | *bim* bim |
        """ |> markdown
        @test md isa Markdown.MD
        @test length(md) == 2
        let p = md[1]
            @test p isa Markdown.Paragraph
            @test length(p.content) == 1
            @test p.content[1] isa Vector
            @test length(p.content[1]) == 3
            @test p.content[1][1] isa Markdown.Code
            @test p.content[1][2] isa String
            @test p.content[1][3] isa Markdown.Code
        end
        let t = md[2]
            @test t isa Markdown.Table
            @test length(t.rows) == 2
            @test length(t.rows[1]) == 2
            @test length(t.rows[1][1]) == 1
            @test length(t.rows[1][2]) == 1
            @test length(t.rows[2]) == 2
            @test length(t.rows[2][1]) == 1
            @test length(t.rows[2][2]) == 2
        end
    end

    @testset "Custom CMarkLibrary" begin
        libcmark = CMark.CMarkLibrary(cmark_jll.libcmark)
        @test CMark.cmark_version(libcmark) == 7682
        @test CMark.cmark_version_string(libcmark) == "0.30.2"
    end
end


text:
@testset "InteractionMatrix tests" begin
    @test_throws ErrorException SpinMC.InteractionMatrix(ones(4,3))

    m1 = SpinMC.InteractionMatrix()
    @test m1.m11 == 0.0
    @test m1.m12 == 0.0
    @test m1.m13 == 0.0
    @test m1.m21 == 0.0
    @test m1.m22 == 0.0
    @test m1.m23 == 0.0
    @test m1.m31 == 0.0
    @test m1.m32 == 0.0
    @test m1.m33 == 0.0

    M = rand(3,3)
    m2 = SpinMC.InteractionMatrix(M)
    @test [m2.m11 m2.m12 m2.m13; m2.m21 m2.m22 m2.m23; m2.m31 m2.m32 m2.m33] == M
end


text:
function filterOutIterations(Wall, Hall, genomeErrors, numberProcessesToExtract, genomesReconstructed, removeLastPercentage)

    totalIterations = div(size(Wall, 2), numberProcessesToExtract)
    totalRemoveIter = convert(Int, round(removeLastPercentage * totalIterations ))

    closenessGenomes = zeros(totalIterations)
    for i = 1 : totalIterations
        closenessGenomes[i] = norm(genomeErrors[:, :, i])
    end

    indexClosenessGenomes = sortperm(closenessGenomes, rev=true)
    removeIterations = indexClosenessGenomes[1:totalRemoveIter]

    removeIterationSets = zeros(Int, numberProcessesToExtract * totalRemoveIter)

    for i = 1 : totalRemoveIter
        iStart = convert(Int, numberProcessesToExtract * ( removeIterations[i] - 1) + 1)
        iEnd = convert(Int, numberProcessesToExtract * removeIterations[i])
        removeIterationSets[(numberProcessesToExtract*(i-1)+1):(numberProcessesToExtract*i)] = iStart:iEnd
    end

    Wall = Wall[:, setdiff(1:end, removeIterationSets)]
    Hall = Hall[setdiff(1:end, removeIterationSets), :]
    genomeErrors = genomeErrors[:, :, setdiff(1:end, removeIterations)]
    genomesReconstructed = genomesReconstructed[:, :, setdiff(1:end, removeIterations)]

    return Wall, Hall, genomeErrors, genomesReconstructed
end

max_stars_repo_name: ypradat/sigprofilerjulia

max_issues_repo_name: ypradat/sigprofilerjulia

max_forks_repo_name: ypradat/sigprofilerjulia


text:
"""
Adds simplistic checkpointing to a DFTK self-consistent field calculation.
Requires JLD2 to be loaded.
"""
function ScfSaveCheckpoints end  # implementation in src/jld2io.jl

"""
Plot the trace of an SCF, i.e. the absolute error of the total energy at
each iteration versus the converged energy in a semilog plot. By default
a new plot canvas is generated, but an existing one can be passed and reused
along with `kwargs` for the call to `plot!`. Requires Plots to be loaded.
"""
function ScfPlotTrace end  # implementation in src/plotting.jl

"""
Default callback function for `self_consistent_field` and `newton`, which prints a convergence table.
"""
function ScfDefaultCallback()
    prev_energy = NaN
    function callback(info)
        show_magn = info.basis.model.spin_polarization == :collinear
        show_diag = hasproperty(info, :diagonalization)
        show_damp = hasproperty(info, :α)

        if show_diag
            # Gather MPI-distributed information
            # Average number of diagonalizations per k-Point needed for this SCF step
            # Note: If two Hamiltonian diagonalizations have been used (e.g. adaptive damping),
            # the per k-Point values are summed.
            diagiter = mpi_mean(sum(mean(diag.iterations) for diag in info.diagonalization),
                                info.basis.comm_kpts)
        end

        !mpi_master() && return info  # Rest is printing => only do on master
        if info.stage == :finalize
            info.converged || @warn "$(info.algorithm) not converged."
            return info
        end

        # TODO We should really do this properly ... this is really messy
        if info.n_iter == 1
            E_label = haskey(info.energies, "Entropy") ? "Free energy" : "Energy"
            label_magn = show_magn ? ("   Magnet", "   ------") : ("", "")
            label_damp = show_damp ? ("   α   ", "   ----") : ("", "")
            label_diag = show_diag ? ("   Diag", "   ----") : ("", "")
            @printf "n     %-12s      Eₙ-Eₙ₋₁     ρout-ρin" E_label
            println(label_magn[1], label_damp[1], label_diag[1])
            @printf "---   ---------------   ---------   --------"
            println(label_magn[2], label_damp[2], label_diag[2])
        end
        E    = isnothing(info.energies) ? Inf : info.energies.total
        Δρ   = norm(info.ρout - info.ρin) * sqrt(info.basis.dvol)
        magn = sum(spin_density(info.ρout)) * info.basis.dvol

        Estr    = (@sprintf "%+15.12f" round(E, sigdigits=13))[1:15]
        ΔE      = isnan(prev_energy) ? "      NaN" : @sprintf "% 3.2e" E - prev_energy
        Mstr    = show_magn ? "   $((@sprintf "%6.3f" round(magn, sigdigits=4))[1:6])" : ""
        diagstr = show_diag ? "  $(@sprintf "% 5.1f" diagiter)" : ""

        αstr = ""
        show_damp && (αstr = isnan(info.α) ? "    NaN" : @sprintf "  % 4.2f" info.α)

        @printf "% 3d   %s   %s   %2.2e" info.n_iter Estr ΔE Δρ
        println(Mstr, αstr, diagstr)
        prev_energy = info.energies.total

        flush(stdout)
        info
    end
end

"""
Flag convergence as soon as total energy change drops below tolerance
"""
function ScfConvergenceEnergy(tolerance)
    energy_total = NaN

    function is_converged(info)
        info.energies === nothing && return false # first iteration

        # The ρ change should also be small, otherwise we converge if the SCF is just stuck
        if norm(info.ρout - info.ρin) * sqrt(info.basis.dvol) > 10sqrt(tolerance)
            return false
        end

        etot_old = energy_total
        energy_total = info.energies.total
        abs(energy_total - etot_old) < tolerance
    end
    return is_converged
end

"""
Flag convergence by using the L2Norm of the change between
input density and unpreconditioned output density (ρout)
"""
function ScfConvergenceDensity(tolerance)
    info -> (norm(info.ρout - info.ρin) * sqrt(info.basis.dvol) < tolerance)
end

"""
Determine the tolerance used for the next diagonalization. This function takes
``|ρnext - ρin|`` and multiplies it with `ratio_ρdiff` to get the next `diagtol`,
ensuring additionally that the returned value is between `diagtol_min` and `diagtol_max`
and never increases.
"""
function ScfDiagtol(;ratio_ρdiff=0.2, diagtol_min=nothing, diagtol_max=0.03)
    function determine_diagtol(info)
        isnothing(diagtol_min) && (diagtol_min = 100eps(real(eltype(info.ρin))))
        info.n_iter ≤ 1 && return diagtol_max
        info.n_iter == 2 && (diagtol_max /= 5)  # Enforce more accurate Bloch wave

        ρnext = hasproperty(info, :ρnext) ? info.ρnext : info.ρout
        diagtol = (norm(ρnext - info.ρin)
                   * sqrt(info.basis.dvol)
                   * ratio_ρdiff)

        # TODO Quantum espresso divides diagtol by the number of electrons
        diagtol = clamp(diagtol, diagtol_min, diagtol_max)
        @assert isfinite(diagtol)

        # Adjust maximum to ensure diagtol may only shrink during an SCF
        diagtol_max = min(diagtol, diagtol_max)
        diagtol
    end
end


text:
using StatsFuns: normcdf, normccdf, normlogcdf, normlogccdf, normlogpdf, normpdf, logsumexp

@marginalrule Probit(:out_in) (m_out::PointMass, m_in::UnivariateNormalDistributionsFamily, meta::ProbitMeta) = begin

    # extract parameters
    mz, vz = mean_cov(m_in)
    p = mean(m_out)
    @assert p >= zero(p) && p <= one(p) "The Probit node only accepts messages on its output with values between 0 and 1."

    # calculate auxiliary variables
    γ = mz/sqrt(1+vz)

    # calculate moments of g
    if γ > 0 && p > 0.5
        log_mom0_pz = logsumexp((log(1-p), log(2*p-1) + normlogccdf(-γ)))
    elseif γ <= 0 && p > 0.5
        log_mom0_pz = logsumexp((log(1-p), log(2*p-1) + normlogcdf(γ)))
    elseif γ > 0 && p <= 0.5
        log_mom0_pz = logsumexp((log(1-p) + normlogcdf(-γ), log(p) + normlogcdf(γ)))
    else
        log_mom0_pz = logsumexp((log(1-p) + normlogccdf(γ), log(p) + normlogcdf(γ)))
    end
    tmp = log(vz) + normlogpdf(γ) - log(1+vz)/2 - log_mom0_pz
    mom1_pz = mz + (2*p-1)*exp(tmp)
    mom2_pz = vz + mz^2 + (2*p-1)*2*mz*exp(tmp) - (2p-1)*γ*exp(log(vz) - log(1 + vz)/2 + tmp)

    # calculate parameters of posterior
    mpz = mom1_pz
    vpz = clamp(mom2_pz - mom1_pz^2, tiny, vz) # ensure variance of marginal is not larger than the variance of the cavity distribution.

    return (out = m_out, in = NormalMeanVariance(mpz, vpz))

end


text:
import Random # TODO: Move `rng` field up to SpinSystem

"""
    mutable struct WangLandau{F<:Function}

Wang-Landau sampler. All parameters have default values that can be overwritten,
but a SpinSystem must be passed during construction.
"""
Base.@kwdef mutable struct WangLandau{F<:Function}
    # adaptive binned histogram
    hist::BinnedArray{Float64, Int64} = BinnedArray{Float64, Int64}()

    # binning resolution for the energy values
    bin_size::Float64 = 0.1

    # optional bounding of state space
    bounds::Vector{Float64} = Vector{Float64}()

    # interval for checking histogram criterion
    hcheck_interval::Int64 = 10000

    # flatness criterion: all hist values must ≥ this fraction of the average value
    flatness::Float64 = 0.8

    # natural log of binned density of states stored adaptively
    ln_g::BinnedArray{Float64, Float64} = BinnedArray{Float64, Float64}()

    # modification factor for accumulating ln_g
    ln_f::Float64 = 1.0

    # function for reduction schedule on ln_f
    ln_f_sched::F = (x, t)->(x/2.0)

    # termination criterion: stop when ln_f ≤ ln_f_final
    ln_f_final::Float64 = 1e-6

    # termination criterion: stop if MC budget exceeded
    max_mcsweeps::Int64 = typemax(Int64)

    # random number generator
    rng::Random.AbstractRNG = Random.MersenneTwister(
        round(Int64, time()*1000)
    )

    # spin system
    system::SpinSystem

    # minimum energy (not binned) found in simulation
    E_min::Float64 = Inf

    # divide quantities by system size
    per_spin::Bool = false

    # max cone radius for spherical cap MC spin move
    mc_step_size::Float64 = 0.1
end

"""
Update a spin `S` by applying a random rotation matrix that has an angle between
0 and ``θ_{max}``, where ``cos(θ_{max})`` is given by the parameter
`cos_max_angle`. Within the constrained space, the probability distribution of
new spin is uniform (with respect to the standard measure on the 2-sphere).

Algorithm adapted from Eqs. 3.60--3.62 of the PhD dissertation "On Classical and
Quantum Mechanical Energy Spectra of Finite Heisenberg Spin Systems", Matthias
Exler https://www.msuq.physik.uni-osnabrueck.de/ps/Dissertation-Exler.pdf
"""

function spherical_cap_update(S::Vec3, cos_max_angle::Float64, rng::Random.AbstractRNG)::Vec3
    # Step 1: Generate a normalized unit vector [x, y, z] from uniform
    # distribution, subject to the constraint that the polar angle θ is less
    # than `max_angle`. Remarkably, this can be achieved by drawing z from a
    # uniform distribution subject to the polar angle constraint.

    # Draw random numbers uniformly from [0,1]
    ξ1 = rand(rng)
    ξ2 = rand(rng)

    # Randomized z component subject to constraint on polar angle
    min_z = cos_max_angle
    z′ = 1 - ξ1 * (1 - min_z)

    # Random azimuthal angle
    ϕ = 2π*ξ2
    sinϕ, cosϕ = sincos(ϕ)

    # Resulting random x and y components
    ρ = sqrt(1 - z′^2)
    x′ = ρ * cosϕ
    y′ = ρ * sinϕ

    # Step 2: Select a reference frame in which S points in the e direction (we
    # will use either e=z or e=y). Specifically, find some rotation matrix R
    # that satisfies `S = R e`. Randomly select a new spin S′ (perturbed from e)
    # in this new reference frame. Then the desired spin update is given by `R
    # S′`.
    x, y, z = S
    if z^2 < 1/2
        # Spin S is not precisely aligned in the z direction, so we can use e=z
        # as our reference frame, and there will not be numerical difficulties
        # in defining the rotation matrix R.
        r = sqrt(x^2 + y^2)
        R = SA[ x*z/r   -y/r    x
                y*z/r    x/r    y
                   -r      0    z]
        S′ = Vec3(x′, y′, z′)
    else
        # Spin S may be precisely aligned with z, and it is safer to pick a
        # different reference frame. We can arbitrarily select e=y, effectively
        # swapping y ↔ z in the code above (also permuting matrix indices).
        r = sqrt(x^2 + z^2)
        R = SA[ x*y/r   x  -z/r
                   -r   y     0
                z*y/r   z   x/r ]
        S′ = Vec3(x′, z′, y′)
    end
    return R*S′
end

"""
Generate a random unit spin that is normally distributed about the direction
of the existing spin S. The parameter σ determines the size of the update
from the spin S. Pass a normal random vector nv so that WLS.rng is used.
"""
function gaussian_spin_update(S::Vec3, σ::Float64, rng::Random.AbstractRNG)::Vec3
    S += σ * randn(rng, Vec3)
    return S/norm(S)
end

"""
Check histogram using the average flatness criterion
"""
function check_hist(hist::BinnedArray{Float64, Int64}; p::Float64=0.8)
    # calculate average of visited bins
    avg = 0.0
    vacancies = 0
    for i in 1:hist.size
        if hist.visited[i]
            avg += hist.vals[i]
        else
            vacancies += 1
        end
    end
    avg /= (hist.size - vacancies)

    # check flatness
    for i in 1:hist.size
        if hist.visited[i] && hist.vals[i] < p*avg
            return false
        end
    end

    return true
end

"""
For new bins, shift ln_g to minimum existing and reset histogram
"""
function add_new!(WLS::WangLandau, key::Float64)
    ln_g_min = Inf
    for i in WLS.hist.size
        if WLS.ln_g.visited[i]
            # find min of ln_g
            if WLS.ln_g.vals[i] < ln_g_min
                ln_g_min = WLS.ln_g.vals[i]
            end
            # reset histogram
            WLS.hist.vals[i] = 0
        end
    end
    # shift new ln_g to min value and add to histogram
    # these will be updated again after acceptance
    WLS.ln_g[key] = ln_g_min - WLS.ln_f
    WLS.hist[key] = 0

    return nothing
end

"""
Initialize system to bounded range of states using throw-away WL sampling run
see run!(...) for comments explaining code in init loop.
"""
function init_bounded!(WLS::WangLandau)
    println("begin init.")

    ln_g_tmp = BinnedArray{Float64, Float64}(bin_size=WLS.bin_size)

    system_size = length(WLS.system)
    ps = (WLS.per_spin) ? system_size : 1

    E_curr = energy(WLS.system) / ps

    ln_g_tmp[E_curr] = 1.0

    lim_curr = Inf
    pfac = (WLS.bounds[1] < E_curr) ? 1 : -1

    # start init with finite length
    for mcsweeps_tmp in 1 : WLS.max_mcsweeps

        for pos in CartesianIndices(WLS.system)
            new_spin = gaussian_spin_update(WLS.system[pos], WLS.mc_step_size, WLS.rng)
            #new_spin = spherical_cap_update(WLS.system[pos], 1.0-WLS.mc_step_size, WLS.rng)

            E_next = E_curr + local_energy_change(WLS.system, pos, new_spin) / ps

            Δln_g = ln_g_tmp[E_curr] - ln_g_tmp[E_next]

            if (Δln_g >= 0) || ( rand(WLS.rng) <= exp(Δln_g) )

                WLS.system[pos] = new_spin
                E_curr = E_next

                if pfac*E_curr < lim_curr
                    lim_curr = pfac*E_curr
                    print("new E = ", E_curr, "\r")
                    flush(stdout)
                end

                if (E_curr >= WLS.bounds[1]) && (E_curr <= WLS.bounds[2])
                    println("\nfinish init with E = ", E_curr)
                    return :SUCCESS
                end
            end

            ln_g_tmp[E_curr] += 1.0
        end
    end
    println("init failed.")
    return :MCSLIMIT
end

"""
Run a Wang-Landau simulation.
"""
function run!(WLS::WangLandau)
    # initialize system if bounded - must supply [min, max]
    bounded = false
    if length(WLS.bounds) == 2
        bounded = true
        if init_bounded!(WLS) == :MCSLIMIT
            return :INITFAIL
        end
    end

    println("begin WL sampling.")

    system_size = length(WLS.system)
    ps = (WLS.per_spin) ? system_size : 1

    iteration = 1
    mcs = 0
    mcsweeps = 0

    # set bin sizes
    WLS.hist.bin_size = WLS.bin_size
    WLS.ln_g.bin_size = WLS.bin_size

    # initial state
    E_curr = energy(WLS.system) / ps

    # record initial state
    WLS.ln_g[E_curr] = WLS.ln_f
    WLS.hist[E_curr] = 1

    # start sampling with finite length
    while (WLS.ln_f > WLS.ln_f_final) && (mcsweeps < WLS.max_mcsweeps)

        # use MC *sweep* as unit time for histogram check interval
        for i in 1 : WLS.hcheck_interval
            for pos in CartesianIndices(WLS.system)
                # propose single spin move - random rotation on spherical cap about spin
                new_spin = gaussian_spin_update(WLS.system[pos], WLS.mc_step_size, WLS.rng)
                #new_spin = spherical_cap_update(WLS.system[pos], 1.0-WLS.mc_step_size, WLS.rng)

                E_next = E_curr + local_energy_change(WLS.system, pos, new_spin) / ps
                mcs += 1

                # enforce bounds if applicable - still update state if rejecting move
                if bounded && ( (E_next < WLS.bounds[1]) || (E_next > WLS.bounds[2]) )
                else
                    # add new bin to ln_g, histogram
                    if WLS.ln_g[E_next] <= eps()
                        add_new!(WLS, E_next)
                    end

                    # calculate ratio of ln_g for transition probability
                    Δln_g = WLS.ln_g[E_curr] - WLS.ln_g[E_next]

                    # accept move
                    if (Δln_g >= 0) || ( rand(WLS.rng) <= exp(Δln_g) )

                        WLS.system[pos] = new_spin
                        E_curr = E_next

                        # record minimum energies
                        if E_curr < WLS.E_min
                            WLS.E_min = E_curr
                            println("mcs = ", mcs, ",  E_min = ", WLS.E_min)
                        end
                    end
                end
                # update ln_g, hist
                WLS.ln_g[E_curr] += WLS.ln_f
                WLS.hist[E_curr] += 1
            end
        end
        mcsweeps += WLS.hcheck_interval

        # check histogram criterion - start new iteration if satisfied
        if check_hist(WLS.hist; p=WLS.flatness)

            # print histogram and ln_g to files for each iteration
            output = open(@sprintf("hist_iteration_%02d.dat", iteration), "w")
            print(output, WLS.hist)
            close(output)

            output = open(@sprintf("ln_g_iteration_%02d.dat", iteration), "w")
            print(output, WLS.ln_g)
            close(output)

            # reset histogram
            reset!(WLS.hist)

            # reduce modification factor by some schedule
            WLS.ln_f = WLS.ln_f_sched(WLS.ln_f, iteration)

            @printf("iteration %d complete: mcs = %d, ln_f = %.8f\n", iteration, mcs, WLS.ln_f)
            iteration += 1
        end
    end

    println("WL sampling complete. mcsweeps = ", mcsweeps)

    return ((mcsweeps < WLS.max_mcsweeps) ? :SUCCESS : :MCSLIMIT)
end


text:
const categoricalplottypes = [BarPlot, Heatmap, Volume]

to_weights(v) = weights(v)
to_weights(v::AbstractWeights) = v

function compute_edges(data, extrema, bins::Tuple{Vararg{Integer}}, closed)
    return map(extrema, bins) do (min, max), n
        histrange(min, max, n, closed)
    end
end
compute_edges(data, extrema, bins::Tuple{Vararg{AbstractArray}}, closed) = bins

function midpoints(edges::AbstractRange)
    min, s, l = minimum(edges), step(edges), length(edges)
    return range(min + s / 2, step=s, length=l - 1)
end

function _histogram(data...; bins=sturges(length(data[1])), weights=automatic,
    normalization=:none, extrema, closed=:left)

    bins_tuple = bins isa Tuple ? bins : map(_ -> bins, data)
    edges = compute_edges(data, extrema, bins_tuple, closed)
    w = weights === automatic ? () : (to_weights(weights),)
    h = fit(Histogram, data, w..., edges)
    return normalize(h, mode=normalization)
end

struct HistogramAnalysis
    options::Dict{Symbol, Any}
end

function (h::HistogramAnalysis)(le::Entry)
    options = copy(h.options)
    get!(options, :extrema, map(extrema, le.positional))

    return splitapply(le) do entry
        hist = _histogram(entry.positional...; entry.named..., options...)
        normalization = get(options, :normalization, :none)
        N = length(entry.positional)
        labels, attributes = copy(entry.labels), copy(entry.attributes)
        labels[N + 1] = normalization == :none ? "count" : string(normalization)
        default_plottype = categoricalplottypes[N]
        plottype = AbstractPlotting.plottype(entry.plottype, default_plottype)
        N == 1 && (attributes[:width] = step(hist.edges[1]))
        positional, named = (map(midpoints, hist.edges)..., hist.weights), (;)
        return Entry(entry; plottype, positional, named, labels, attributes)
    end
end

"""
    histogram(; bins=automatic, weights=automatic, normalization=:none)

Compute a histogram. `bins` can be an `Int` to create that
number of equal-width bins over the range of `values`.
Alternatively, it can be a sorted iterable of bin edges. The histogram
can be normalized by setting `normalization`. Possible values are:
*  `:pdf`: Normalize by sum of weights and bin sizes. Resulting histogram
   has norm 1 and represents a PDF.
* `:density`: Normalize by bin sizes only. Resulting histogram represents
   count density of input and does not have norm 1.
* `:probability`: Normalize by sum of weights only. Resulting histogram
   represents the fraction of probability mass for each bin and does not have
   norm 1.
*  `:none`: Do not normalize.
Weighted data is supported via the keyword `weights`.

!!! note

    Normalizations are computed withing groups. For example, in the case of
    `normalization=:pdf`, sum of weights *within each group* will be equal to `1`.
"""
histogram(; options...) = Layer((HistogramAnalysis(Dict{Symbol, Any}(options)),))


text:
@everywhere using DistributedArrays
include("../utils.jl")

# Get the cell id for a given x. Note cells are labeled 1,2,3,...
@everywhere function getCellID1d(x::Float64, nCells::Int64, nCols::Int64)
  id::Int64 = floor(Int64, x * nCells / nCols)+1
  id<0 ? 0 : (id>nCells ? nCells : id) # make it between 0 and nCells
end

# Returns (xTess, velTess)
@everywhere function getTess(As::Array{Float64,3}, nCells::Int64, nCols::Int64)
  colPerCell::Float64 = nCols/nCells
  velTess::Array{Float64,1} = zeros(nCells+1)
  xTess::Array{Float64,1} = zeros(nCells+1)
  for i::Int64=1:nCells
    xTess[i+1]::Float64 = i*colPerCell
    velTess[i+1]::Float64 = getVelocity1d(xTess[i], As, i)
  end
  xTess, velTess
end

# Assuming 0 boundary
@everywhere function getVelocity1d(x::Float64, As::Array{Float64,3}, cellId::Int64)
  As[cellId,1,1]*x + As[cellId,1,2]
end

# loop though all points to get velocity field
@everywhere function getVelocityField1d(pts_src::Array{Float64,1}, As::Array{Float64,3}, nCells::Int64, nCols::Int64)
  nPts::Int64 = length(pts_src)
  v::Array{Float64,1} = zeros(nPts)
  for i::Int64=1:nPts
    x::Float64 = pts_src[i]
    cellId::Int64 = getCellID1d(x, nCells, nCols)
    v[i]::Float64 = getVelocity1d(x, As, cellId)
  end
  v
end

@everywhere function calculateAffinePhi(x::Float64, t::Float64, cellId::Int64, As::Array{Float64,3})
  a::Float64, b::Float64 = As[cellId,1,1], As[cellId,1,2]

  xNext::Float64
  if a == 0
    xNext = x + t * b
  else
    xNext = exp(t*a)*x + b*(exp(t*a)-1)/a
  end
  xNext
end

@everywhere function findCrossingTime(x::Float64, cellId::Int64, As::Array{Float64,3}, xTess::Array{Float64,1}, nCols::Int64)
  v::Float64 = getVelocity1d(x, As, cellId)
  a::Float64, b::Float64 = As[cellId,1,1], As[cellId,1,2]

  # Calculate next cell, velocity
  cellIdNext::Int64 = cellId + round(Int64, sign(v))
  if !(1 <= cellIdNext <= nCells)
    return Inf, Union{}
  end

  xNext::Float64
  if v > 0
    xNext = xTess[cellIdNext] + 1e-16 # pick the right endpoint
  else
    xNext = xTess[cellId] - 1e-16 # pick the left endpoint
  end
  vNext::Float64 = getVelocity1d(xNext, As, cellIdNext)

  # Cases where we will never make it to the next cell
  if allclose(vNext, 0) || sign(v) != sign(vNext)
    #println("b")
    return Inf, Union{}
  end

  t::Float64
  if allclose(xNext, 0) || allclose(xNext, nCols)
    t = Inf
  elseif allclose(a, 0)
    t = (xNext - x)/b
  else
    t = 1/a * log( (xNext + b/a)/(x + b/a) )
  end
  t, cellIdNext
end

@everywhere function calculatePhi(x::Float64, t::Float64, tess::Tuple{Array{Float64,1},Array{Float64,1}}, As::Array{Float64,3}, nCells::Int64, nCols::Int64)
  xTess::Array{Float64,1} = tess[1]
  velTess::Array{Float64,1} = tess[2]

  xPrev::Float64 = x
  xNext::Float64 = 0.0
  #counter::Int64 = 0
  tRemainder::Float64 = t
  v::Float64

  #
  if 0.0<xPrev<nCols && in(xPrev, xTess)
    v = velTess[findfirst(xTess, xPrev)]
    if v>0.0
        pass
    else
      #println("line 105")
      xPrev = xPrev-1e-16
    end
  end

  # Do the function compositions
  while true
    cellId::Int64 = getCellID1d(xPrev, nCells, nCols)
    tCell::Float64, cellIdNext::Any = findCrossingTime(xPrev, cellId, As, xTess, nCols)

    # Whether or not we stay in the same cell
    crossed::Bool
    if tCell >= tRemainder
      tCell = tRemainder
      crossed = false
    else
      crossed = true
    end

    if crossed
      if cellIdNext > cellId
        # Moving to next cell
        xNext = xTess[cellIdNext] + 1e-12
      else
        # going to previous cell
        xNext = xTess[cellId] - 1e-12
      end
    else
      xNext = calculateAffinePhi(xPrev, tCell, cellId, As)
    end

    # Making it stay within 0 and nCols
    xNext = min(nCols, max(0, xNext))

    tRemainder::Float64 -= tCell

    if !crossed
      break
    end

    xPrev = xNext
  end
  xNext
end

@everywhere function calculateNumericalPhi(
    x::Float64,
    t::Float64,
    tess::Tuple{Array{Float64,1},Array{Float64,1}},
    As::Array{Float64,3},
    nCells::Int64,
    nCols::Int64,
    nSteps::Int64
  )

  yn::Float64 = x
  Δt::Float64 = t/nSteps
  for i::Int64 = 1:nSteps
    midpoint::Float64 = yn+Δt/2*getVelocity1d(yn, As, getCellID1d(yn, nCells, nCols))
    cellId::Int64 = getCellID1d(midpoint, nCells, nCols)
    yn = yn + Δt*getVelocity1d(midpoint, As, cellId)
  end
  yn
end

# loop though all points to get integration
function calculateNumericalIntegration(
    pts_src::Array{Float64,1},
    t::Float64,
    tess::Tuple{Array{Float64,1},Array{Float64,1}},
    As::Array{Float64,3},
    nCells::Int64,
    nCols::Int64,
    nSteps::Int64
  )

  nPts::Int64 = length(pts_src)
  pts_fwd::Array{Float64,1} = zeros(nPts)
  x::Float64 = 0.
  for i::Int64=1:nPts
    x = pts_src[i]
    pts_fwd[i] = calculateNumericalPhi(x, t, tess, As, nCells, nCols, nSteps)
  end
  pts_fwd
end

# loop though all points to get integration
function calculateNumericalIntegrationParallel(
    pts_src::Array{Float64,1},
    t::Float64,
    tess::Tuple{Array{Float64,1},Array{Float64,1}},
    As::Array{Float64,3},
    nCells::Int64,
    nCols::Int64,
    nSteps::Int64
  )

  nPts::Int64 = length(pts_src)
  pts_fwd::DArray{Float64,1} = dzeros(nPts)
  threads::Array{Int64,1} = workers()
  nThreads::Int64 = length(threads)
  refs::Array{RemoteRef,1} = Array(RemoteRef, nThreads)

  #spawns the tasks at each worker
  for index::Int64=1:nThreads
    p = threads[index]
    refs[index] = @spawnat(p, (
      indicesOnThread::Array{Int64,1} = localindexes(pts_fwd)[1];
      localArray = localpart(pts_fwd);
      for j::Int64 = 1:length(indicesOnThread);
        localArray[j] = calculateNumericalPhi(pts_src[indicesOnThread[j]], t, tess, As, nCells, nCols, nSteps);
      end;
      ))
  end

  # wait for the workers to finish
  map(wait, refs)

  pts_fwd
end

# loop though all points to get integration
function calculateIntegration(
    pts_src::Array{Float64,1},
    t::Float64,
    tess::Tuple{Array{Float64,1},Array{Float64,1}},
    As::Array{Float64,3},
    nCells::Int64,
    nCols::Int64
  )

  nPts::Int64 = length(pts_src)
  pts_fwd::Array{Float64,1} = zeros(nPts)
  for i::Int64=1:nPts
    x::Float64 = pts_src[i]
    pts_fwd[i] = calculatePhi(x, t, tess, As, nCells, nCols)
  end
  pts_fwd
end

function calculateIntegrationParallel(
    pts_src::Array{Float64,1},
    t::Float64,
    tess::Tuple{Array{Float64,1},Array{Float64,1}},
    As::Array{Float64,3},
    nCells::Int64,
    nCols::Int64
  )

  nPts::Int64 = length(pts_src)
  pts_fwd::DArray{Float64,1} = dzeros(nPts)
  threads::Array{Int64,1} = workers()
  nThreads::Int64 = length(threads)
  refs::Array{RemoteRef,1} = Array(RemoteRef, nThreads)

  # spawns the tasks at each worker
  for index::Int64=1:nThreads
    p = threads[index]
    refs[p] = @spawnat(p, (
      indicesOnThread::Array{Int64,1} = localindexes(pts_fwd)[1];
      localArray = localpart(pts_fwd);
      for j::Int64 = 1:length(indicesOnThread);
        localArray[j] = calculatePhi(pts_src[indicesOnThread[j]], t, tess, As, nCells, nCols);
      end;
      ))
  end

  # wait for the workers to finish
  map(wait, refs)

  pts_fwd
end

@everywhere function solver1d(
    x::Float64,
    t::Float64,
    tess::Tuple{Array{Float64,1},Array{Float64,1}},
    As::Array{Float64,3},
    nCells::Int64,
    nCols::Int64,
    NSteps::Int64,
    Δt::Float64,
    nSteps::Int64
  )

  xPrev::Float64 = x
  cellId::Int64 = getCellID1d(xPrev, nCells, nCols)
  for j::Int64 = 1:NSteps
    xTemp::Float64 = calculateAffinePhi(xPrev, Δt, cellId, As)
    if cellId == getCellID1d(xTemp, nCells, nCols)
      xPrev = xTemp
    else
      xPrev = calculateNumericalPhi(xPrev, Δt, tess, As, nCells, nCols, nSteps)
      cellId == getCellID1d(xPrev, nCells, nCols)
    end
  end
  xPrev
end

function solver1dIntegration(
    pts_src::Array{Float64,1},
    t::Float64,
    tess::Tuple{Array{Float64,1},Array{Float64,1}},
    As::Array{Float64,3},
    nCells::Int64,
    nCols::Int64,
    NSteps::Int64,
    nSteps::Int64
  )

  nPts::Int64 = length(pts_src)
  pts_fwd::Array{Float64,1} = zeros(nPts)
  Δt::Float64 = t/NSteps
  for i::Int64=1:nPts
    pts_fwd[i] = solver1d(pts_src[i], t, tess, As, nCells, nCols, NSteps, Δt, nSteps)
  end
  pts_fwd
end

function solver1dParallelIntegration(
    pts_src::Array{Float64,1},
    t::Float64,
    tess::Tuple{Array{Float64,1},Array{Float64,1}},
    As::Array{Float64,3},
    nCells::Int64,
    nCols::Int64,
    NSteps::Int64,
    nSteps::Int64
  )

  nPts::Int64 = length(pts_src)
  pts_fwd::DArray{Float64,1} = dzeros(nPts)
  Δt::Float64 = t/NSteps
  δt::Float64 = Δt/nSteps
  threads::Array{Int64,1} = workers()
  nThreads::Int64 = length(threads)
  refs::Array{RemoteRef,1} = Array(RemoteRef, nThreads)

  # spawns the tasks at each worker
  for index::Int64=1:nThreads
    p = threads[index]
    refs[index] = @spawnat(p, (
      indicesOnThread::Array{Int64,1} = localindexes(pts_fwd)[1];
      localArray = localpart(pts_fwd);
      for i::Int64 = 1:length(indicesOnThread);
        localArray[i] = solver1d(pts_src[indicesOnThread[i]], t, tess, As, nCells, nCols, NSteps, Δt, nSteps);
      end;
      ))
  end

  # wait for the workers to finish
  map(wait, refs)

  pts_fwd
end


text:
module MochaTheano
import Mocha: Config
export TheanoLossLayer
include("theano-layer.jl")

if Config.use_cuda
  include("cuda/theano-layer.jl")
end
end

max_stars_repo_name: benmoran/MochaTheano

max_issues_repo_name: benmoran/MochaTheano

max_forks_repo_name: benmoran/MochaTheano


text:
using StackOverflow

"""
    getquestions(sort::string = "creation", order::String = "desc")

Returns 30 most recent answers given a certain sort parameter.
"""
function getquestions(sort::String = "creation", order::String = "desc")

    r = HTTP.request("GET", "https://api.stackexchange.com/2.2/answers?order=$(order)&sort=$(sort)&site=stackoverflow")

    json_obj = StackOverflow.convert_HTTP_Response_To_JSON(r)

    for (k, v) in json_obj
        if occursin("items", k)
            for item in v
                for (key, value) in item
                    println("$(key): $(value)")
                end
                println("\n")
            end
        end
    end

end


text:
import UUIDs

if VERSION < v"1.4"
    stdlibs = Pkg.Types.stdlib
else
    # renamed to stdlibs
    # https://github.com/JuliaLang/Pkg.jl/pull/1559
    stdlibs = Pkg.Types.stdlibs
end

# This function is based off of a similar function here:
# https://github.com/JuliaRegistries/RegistryCI.jl/blob/master/src/RegistryCI.jl
function gather_stdlib_uuids()
    return Set{UUIDs.UUID}(x for x in keys(stdlibs()))
end


text:
using CUDArt

include("init.jl")
include("defs.jl")
include("wrappers.jl")
include("hl.jl")

# println("At the end of module loading: $_rng")