text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: February 2021

function getindex(J::HDF5.Dataset, ::Colon, ::Colon, ::Colon, a::Array{Int64,1})
    X = zeros(Float32, size(J, 1), size(J, 2), size(J, 3), length(a))

    for (i, idx) in enumerate(a)
        X[:, :, :, i] = J[:, :, :, idx]
    end
    return X
end


text:
function func2(x)
return x^2
end

function func1(x)
    function func3(x)
        return func4(x)
    end
return func3(x)
end

function anotherfun(x)
    return awayfun(x)
end

function nonanonfun(x)
    return anotherfun(x)
end

anonfun = (x::Float64) -> nonanonfun(x)
anonfun = (x::Int64) -> nonanonfun(x)

ext: jl

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
import PowerSystems: LazyDictFromIterator

@testset "PowerSystemTableData parsing" begin
    resolutions = (
        (resolution = Dates.Minute(5), len = 288),
        (resolution = Dates.Minute(60), len = 24),
    )

    for (resolution, len) in resolutions
        create_func = [create_rts_system, create_rts_multistart_system]
        for f in create_func
            sys = f(resolution)
            for time_series in get_time_series_multiple(sys)
                @test length(time_series) == len
            end
        end
    end
end

@testset "PowerSystemTableData parsing invalid directory" begin
    baddir = abspath(joinpath(RTS_GMLC_DIR, "../../test"))
    @test_throws ErrorException PowerSystemTableData(baddir, 100.0, DESCRIPTORS)
end

@testset "consistency between PowerSystemTableData and standardfiles" begin
    mpsys = PSB.build_system(PSB.MatPowerTestSystems, "matpower_RTS_GMLC_sys")
    cdmsys = PSB.build_system(PSB.PSITestSystems, "test_RTS_GMLC_sys")
    mp_iter = get_components(HydroGen, mpsys)
    mp_generators = LazyDictFromIterator(String, HydroGen, mp_iter, get_name)
    for cdmgen in get_components(HydroGen, cdmsys)
        mpgen = get(mp_generators, uppercase(get_name(cdmgen)))
        if isnothing(mpgen)
            error("did not find $cdmgen")
        end
        @test cdmgen.available == mpgen.available
        @test lowercase(cdmgen.bus.name) == lowercase(mpgen.bus.name)
        gen_dat = (
            structname = nothing,
            fields = (
                :active_power,
                :reactive_power,
                :rating,
                :active_power_limits,
                :reactive_power_limits,
                :ramp_limits,
            ),
        )
        function check_fields(chk_dat)
            for field in chk_dat.fields
                n = get(chk_dat, :structname, nothing)
                (cdmd, mpd) =
                    isnothing(n) ? (cdmgen, mpgen) :
                    (getfield(cdmgen, n), getfield(mpgen, n))
                cdmgen_val = getfield(cdmd, field)
                mpgen_val = getfield(mpd, field)
                if isnothing(cdmgen_val) || isnothing(mpgen_val)
                    @warn "Skip value with nothing" repr(cdmgen_val) repr(mpgen_val)
                    continue
                end
                @test cdmgen_val == mpgen_val
            end
        end
        check_fields(gen_dat)
    end

    mp_iter = get_components(ThermalGen, mpsys)
    mp_generators = LazyDictFromIterator(String, ThermalGen, mp_iter, get_name)
    for cdmgen in get_components(ThermalGen, cdmsys)
        mpgen = get(mp_generators, uppercase(get_name(cdmgen)))
        @test cdmgen.available == mpgen.available
        @test lowercase(cdmgen.bus.name) == lowercase(mpgen.bus.name)
        for field in (:active_power_limits, :reactive_power_limits, :ramp_limits)
            cdmgen_val = getfield(cdmgen, field)
            mpgen_val = getfield(mpgen, field)
            if isnothing(cdmgen_val) || isnothing(mpgen_val)
                @warn "Skip value with nothing" repr(cdmgen_val) repr(mpgen_val)
                continue
            end
            @test cdmgen_val == mpgen_val
        end

        if length(mpgen.operation_cost.variable) == 4
            @test [
                isapprox(
                    cdmgen.operation_cost.variable[i][1],
                    mpgen.operation_cost.variable[i][1],
                    atol = 0.1,
                ) for i in 1:4
            ] == [true, true, true, true]
            #@test compare_values_without_uuids(cdmgen.operation_cost, mpgen.operation_cost)
        end
    end

    mp_iter = get_components(RenewableGen, mpsys)
    mp_generators = LazyDictFromIterator(String, RenewableGen, mp_iter, get_name)
    for cdmgen in get_components(RenewableGen, cdmsys)
        mpgen = get(mp_generators, uppercase(get_name(cdmgen)))
        # Disabled since data is inconsisten between sources
        #@test cdmgen.available == mpgen.available
        @test lowercase(cdmgen.bus.name) == lowercase(mpgen.bus.name)
        for field in (:rating, :power_factor)
            cdmgen_val = getfield(cdmgen, field)
            mpgen_val = getfield(mpgen, field)
            if isnothing(cdmgen_val) || isnothing(mpgen_val)
                @warn "Skip value with nothing" repr(cdmgen_val) repr(mpgen_val)
                continue
            end
            @test cdmgen_val == mpgen_val
        end
        #@test compare_values_without_uuids(cdmgen.operation_cost, mpgen.operation_cost)
    end

    cdm_ac_branches = collect(get_components(ACBranch, cdmsys))
    @test get_rate(cdm_ac_branches[2]) == get_rate(get_branch(mpsys, cdm_ac_branches[2]))
    @test get_rate(cdm_ac_branches[6]) == get_rate(get_branch(mpsys, cdm_ac_branches[6]))
    @test get_rate(cdm_ac_branches[120]) ==
          get_rate(get_branch(mpsys, cdm_ac_branches[120]))

    cdm_dc_branches = collect(get_components(DCBranch, cdmsys))
    @test get_active_power_limits_from(cdm_dc_branches[1]) ==
          get_active_power_limits_from(get_branch(mpsys, cdm_dc_branches[1]))
end

@testset "Test reserve direction" begin
    @test PSY.get_reserve_direction("Up") == ReserveUp
    @test PSY.get_reserve_direction("Down") == ReserveDown

    for invalid in ("up", "down", "right", "left")
        @test_throws PSY.DataFormatError PSY.get_reserve_direction(invalid)
    end
end

@testset "Test consistency between variable cost and heat rate parsing" begin
    fivebus_dir = joinpath(DATA_DIR, "5-bus-hydro")
    rawsys_hr = PowerSystemTableData(
        fivebus_dir,
        100.0,
        joinpath(fivebus_dir, "user_descriptors_var_cost.yaml"),
        generator_mapping_file = joinpath(fivebus_dir, "generator_mapping.yaml"),
    )
    rawsys = PowerSystemTableData(
        fivebus_dir,
        100.0,
        joinpath(fivebus_dir, "user_descriptors_var_cost.yaml"),
        generator_mapping_file = joinpath(fivebus_dir, "generator_mapping.yaml"),
    )
    sys_hr = System(rawsys_hr)
    sys = System(rawsys)

    g_hr = get_components(ThermalStandard, sys_hr)
    g = get_components(ThermalStandard, sys)
    @test get_variable.(get_operation_cost.(g)) == get_variable.(get_operation_cost.(g))
end


text:
using Turing
using DifferentialEquations

# Load StatsPlots for visualizations and diagnostics.
using StatsPlots

using LinearAlgebra

# Set a seed for reproducibility.
using Random
Random.seed!(14);

# Define Lotka-Volterra model.
function lotka_volterra(du, u, p, t)
    # Model parameters.
    α, β, γ, δ = p
    # Current state.
    x, y = u

    # Evaluate differential equations.
    du[1] = (α - β * y) * x # prey
    du[2] = (δ * x - γ) * y # predator

    return nothing
end

# Define initial-value problem.
u0 = [1.0, 1.0]
p = [1.5, 1.0, 3.0, 1.0]
tspan = (0.0, 10.0)
prob = ODEProblem(lotka_volterra, u0, tspan, p)

# Plot simulation.
plot(solve(prob, Tsit5()))

sol = solve(prob, Tsit5(); saveat=0.1)
odedata = Array(sol) + 0.8 * randn(size(Array(sol)))

# Plot simulation and noisy observations.
plot(sol; alpha=0.3)
scatter!(sol.t, odedata'; color=[1 2], label="")

@model function fitlv(data, prob)
    # Prior distributions.
    σ ~ InverseGamma(2, 3)
    α ~ truncated(Normal(1.5, 0.5), 0.5, 2.5)
    β ~ truncated(Normal(1.2, 0.5), 0, 2)
    γ ~ truncated(Normal(3.0, 0.5), 1, 4)
    δ ~ truncated(Normal(1.0, 0.5), 0, 2)

    # Simulate Lotka-Volterra model.
    p = [α, β, γ, δ]
    predicted = solve(prob, Tsit5(); p=p, saveat=0.1)

    # Observations.
    for i in 1:length(predicted)
        data[:, i] ~ MvNormal(predicted[i], σ^2 * I)
    end

    return nothing
end

model = fitlv(odedata, prob)

# Sample 3 independent chains with forward-mode automatic differentiation (the default).
chain = sample(model, NUTS(0.65), MCMCSerial(), 1000, 3; progress=false)

plot(chain)

plot(; legend=false)
posterior_samples = sample(chain[[:α, :β, :γ, :δ]], 300; replace=false)
for p in eachrow(Array(posterior_samples))
    sol_p = solve(prob, Tsit5(); p=p, saveat=0.1)
    plot!(sol_p; alpha=0.1, color="#BBBBBB")
end

# Plot simulation and noisy observations.
plot!(sol; color=[1 2], linewidth=1)
scatter!(sol.t, odedata'; color=[1 2])

@model function fitlv2(data::AbstractVector, prob)
    # Prior distributions.
    σ ~ InverseGamma(2, 3)
    α ~ truncated(Normal(1.5, 0.5), 0.5, 2.5)
    β ~ truncated(Normal(1.2, 0.5), 0, 2)
    γ ~ truncated(Normal(3.0, 0.5), 1, 4)
    δ ~ truncated(Normal(1.0, 0.5), 0, 2)

    # Simulate Lotka-Volterra model but save only the second state of the system (predators).
    p = [α, β, γ, δ]
    predicted = solve(prob, Tsit5(); p=p, saveat=0.1, save_idxs=2)

    # Observations of the predators.
    data ~ MvNormal(predicted.u, σ^2 * I)

    return nothing
end

model2 = fitlv2(odedata[2, :], prob)

# Sample 3 independent chains.
chain2 = sample(model2, NUTS(0.45), MCMCSerial(), 5000, 3; progress=false)

plot(; legend=false)
posterior_samples = sample(chain2[[:α, :β, :γ, :δ]], 300; replace=false)
for p in eachrow(Array(posterior_samples))
    sol_p = solve(prob, Tsit5(); p=p, saveat=0.1)
    plot!(sol_p; alpha=0.1, color="#BBBBBB")
end

# Plot simulation and noisy observations.
plot!(sol; color=[1 2], linewidth=1)
scatter!(sol.t, odedata'; color=[1 2])

function delay_lotka_volterra(du, u, h, p, t)
    # Model parameters.
    α, β, γ, δ = p

    # Current state.
    x, y = u
    # Evaluate differential equations
    du[1] = α * h(p, t - 1; idxs=1) - β * x * y
    du[2] = -γ * y + δ * x * y

    return nothing
end

# Define initial-value problem.
p = (1.5, 1.0, 3.0, 1.0)
u0 = [1.0; 1.0]
tspan = (0.0, 10.0)
h(p, t; idxs::Int) = 1.0
prob_dde = DDEProblem(delay_lotka_volterra, u0, h, tspan, p);

sol_dde = solve(prob_dde; saveat=0.1)
ddedata = Array(sol_dde) + 0.5 * randn(size(sol_dde))

# Plot simulation and noisy observations.
plot(sol_dde)
scatter!(sol_dde.t, ddedata'; color=[1 2], label="")

@model function fitlv_dde(data, prob)
    # Prior distributions.
    σ ~ InverseGamma(2, 3)
    α ~ Truncated(Normal(1.5, 0.5), 0.5, 2.5)
    β ~ Truncated(Normal(1.2, 0.5), 0, 2)
    γ ~ Truncated(Normal(3.0, 0.5), 1, 4)
    δ ~ Truncated(Normal(1.0, 0.5), 0, 2)

    # Simulate Lotka-Volterra model.
    p = [α, β, γ, δ]
    predicted = solve(prob, MethodOfSteps(Tsit5()); p=p, saveat=0.1)

    # Observations.
    for i in 1:length(predicted)
        data[:, i] ~ MvNormal(predicted[i], σ^2 * I)
    end
end

model_dde = fitlv_dde(ddedata, prob_dde)

# Sample 3 independent chains.
chain_dde = sample(model_dde, NUTS(0.65), MCMCSerial(), 300, 3; progress=false)

plot(chain_dde)

plot(; legend=false)
posterior_samples = sample(chain_dde[[:α, :β, :γ, :δ]], 300; replace=false)
for p in eachrow(Array(posterior_samples))
    sol_p = solve(prob_dde, MethodOfSteps(Tsit5()); p=p, saveat=0.1)
    plot!(sol_p; alpha=0.1, color="#BBBBBB")
end

# Plot simulation and noisy observations.
plot!(sol_dde; color=[1 2], linewidth=1)
scatter!(sol_dde.t, ddedata'; color=[1 2])

using Zygote, DiffEqSensitivity

# Sample a single chain with 1000 samples using Zygote.
setadbackend(:zygote)
sample(model, NUTS(0.65), 1000; progress=false)

@model function fitlv_sensealg(data, prob)
    # Prior distributions.
    σ ~ InverseGamma(2, 3)
    α ~ truncated(Normal(1.5, 0.5), 0.5, 2.5)
    β ~ truncated(Normal(1.2, 0.5), 0, 2)
    γ ~ truncated(Normal(3.0, 0.5), 1, 4)
    δ ~ truncated(Normal(1.0, 0.5), 0, 2)

    # Simulate Lotka-Volterra model and use a specific algorithm for computing sensitivities.
    p = [α, β, γ, δ]
    predicted = solve(
        prob;
        p=p,
        saveat=0.1,
        sensealg=InterpolatingAdjoint(; autojacvec=ReverseDiffVJP(true)),
    )

    # Observations.
    for i in 1:length(predicted)
        data[:, i] ~ MvNormal(predicted[i], σ^2 * I)
    end

    return nothing
end;

model_sensealg = fitlv_sensealg(odedata, prob)

# Sample a single chain with 1000 samples using Zygote.
setadbackend(:zygote)
sample(model_sensealg, NUTS(0.65), 1000; progress=false)

u0 = [1.0, 1.0]
tspan = (0.0, 10.0)
function multiplicative_noise!(du, u, p, t)
    x, y = u
    du[1] = p[5] * x
    return du[2] = p[6] * y
end
p = [1.5, 1.0, 3.0, 1.0, 0.1, 0.1]

function lotka_volterra!(du, u, p, t)
    x, y = u
    α, β, γ, δ = p
    du[1] = dx = α * x - β * x * y
    return du[2] = dy = δ * x * y - γ * y
end

prob_sde = SDEProblem(lotka_volterra!, multiplicative_noise!, u0, tspan, p)

ensembleprob = EnsembleProblem(prob_sde)
data = solve(ensembleprob, SOSRI(); saveat=0.1, trajectories=1000)
plot(EnsembleSummary(data))

@model function fitlv_sde(data, prob)
    # Prior distributions.
    σ ~ InverseGamma(2, 3)
    α ~ truncated(Normal(1.3, 0.5), 0.5, 2.5)
    β ~ truncated(Normal(1.2, 0.25), 0.5, 2)
    γ ~ truncated(Normal(3.2, 0.25), 2.2, 4.0)
    δ ~ truncated(Normal(1.2, 0.25), 0.5, 2.0)
    ϕ1 ~ truncated(Normal(0.12, 0.3), 0.05, 0.25)
    ϕ2 ~ truncated(Normal(0.12, 0.3), 0.05, 0.25)

    # Simulate stochastic Lotka-Volterra model.
    p = [α, β, γ, δ, ϕ1, ϕ2]
    predicted = solve(prob, SOSRI(); p=p, saveat=0.1)

    # Early exit if simulation could not be computed successfully.
    if predicted.retcode !== :Success
        Turing.@addlogprob! -Inf
        return nothing
    end

    # Observations.
    for i in 1:length(predicted)
        data[:, i] ~ MvNormal(predicted[i], σ^2 * I)
    end

    return nothing
end;

model_sde = fitlv_sde(odedata, prob_sde)

setadbackend(:forwarddiff)
chain_sde = sample(
    model_sde,
    NUTS(0.25),
    5000;
    init_params=[1.5, 1.3, 1.2, 2.7, 1.2, 0.12, 0.12],
    progress=false,
)
plot(chain_sde)

max_stars_repo_name: TuringLang/TuringTutorialsOutput

max_issues_repo_name: TuringLang/TuringTutorialsOutput

max_forks_repo_name: TuringLang/TuringTutorialsOutput


text:
@testset "lexicon plain text IO" begin
    # test basic lexicon reading
    @testset "basic lexicon reading" begin
        lex = Lexicon("res/lexicon-0.txt")
        for elem in ["word", "tagless", "target", "lexicon", "pencil", "pen", "manytags"]
            @test elem in lex
        end
        for elem in ["word", "target", "pencil", "pen"]
            @test elem in tagged_with(lex, "tag2")
        end
        for  elem in ["tagless", "lexicon", "manytags"]
            @test !(elem in tagged_with(lex, "tag2"))
        end
        for tag in ["tag1", "tag2", "meta", "unique", "tag", "another", "more", "evenmore", "lots",
        "toomany"]
            @test tag in lex.tags
        end
        for tag in ["tag", "another", "more", "evenmore", "lots", "toomany"]
            @test tag in lex["manytags"]
        end
    end

    # words and tags with special characters
    @testset "reading lexicons with special chars" begin
        lex = Lexicon("res/lexicon-1.txt")
        for elem in ["with space", "with,comma", "multipart,sentence with both", "word"]
            @test elem in lex
        end
        for tag in ["symbol", "whitespace"]
            @test tag in lex["multipart,sentence with both"]
        end
        @test "spaced tag" in lex["word"]
        @test !("spaced" in lex["word"])
        @test !("tag" in lex["word"])
    end

    # writing and reading a lexicon from file
    @testset "lexicon writing and reading" begin
        lex = Lexicon()
        push!(lex, "command", "exception", "ide", tags=("computer",))
        push!(lex, "don't know where this word comes from")
        push!(lex, "pencil", tags=("stationary",))
        push!(lex, "two category item", tags=("category1", "category2"))

        file = IOBuffer()
        println(file, lex)
        read_lex = Lexicon(file)

        for elem in ["command", "exception", "ide", "pencil"]
            @test elem in lex
        end
        for tag in ["category1", "category2"]
            @test tag in lex["two category item"]
        end
        @test "stationary" in lex["pencil"]
        for elem in ["command", "exception", "ide"]
            @test "computer" in lex[elem]
        end
    end

end

@testset "lexicon IO (other formats)" begin

end

@testset "coverage under lexicon" begin

end


text:
"""
sample_theta()

Description of sample_theta
------------------------------
Returns θ in radians Uniformly distributed. θ is the polar angle. See Distributions.Uniform

"""
function sample_theta()
    return acos(1 - 2 * rand(Uniform())) ## polar angle in radians
end

"""
    sample_phi()

Description of sample_phi
------------------------------
Returns ϕ in radians Uniformly distributed. phi is the azimuthal angle. See Distributions.Uniform

Sampling equation is:
    + First θ1 is sampled uniformly from 2π uniform pdf.
    + Second angle θdif (the angle between the two electrons) is sampled from pdf = 0.5 - 0.5x. This pdf represents the 1 - cos(θ) angular distribution.
        It is sampled through inverse CDF method. Integrating pdf and finding the inverse function we get ``` x = 1 - 2\\sqrt{1 - \\gamma}```
    + θ2 is determined as θ2 = θdif + θ1.

"""
function sample_phi()
    return 2*π*rand(Uniform()) ## azimuthal angle in radians
end

"""
    sample_theta_dif()

Description of sample_theta_dif
------------------------------
Returns θdif in radians distrtibuted according to pdf = 0.5 - 0.5x (where x = cos(θdif)). θdif is the angle between the two electrons.

Sampling equation is:
    + θdif (the angle between the two electrons) is sampled from pdf = 0.5 - 0.5x. This pdf represents the 1 - cos(θ) angular distribution.
        It is sampled through inverse CDF method. Integrating pdf and finding the inverse function we get ``` x = 1 - 2\\sqrt{1 - \\gamma}```
    + returned angle is the acos(x)

"""
function sample_theta_dif()
    a = 1.
    b = -2.
    c = -3 + 4*rand(Uniform())
    cosθdif = solve_quadratic(a,b,c)[1]
    return acos(cosθdif)
end

"""
#### function ```sample_energies(df::DataFrame, thickness = 0.001)```

    Description of ```sample_energies```
    ------------------------------
Returns a tuple of energies (E1, E2) given by the distribution specified in dataframe ```df```.

```df``` must have fields: ```:E1, :minE, :maxE, :a, :b, :cdf (given by get_cdf!(df))```

Sampling is done in the following way.
E2 is sampled using inverse CDF method - where cdf is obtained from ```get_cdf!(df)```.
E1 is dependent on E2 as it can only be sampled in the given row where E2 was chosen.
E1 is sampled as a random uniform number between (row of E1, next row).

It takes roughly 100s to sample 1e5 energies.

"""
function sample_energies(df::DataFrame, thickness = 0.001)
    gamma = rand(Uniform(0.0, maximum(df.cdf))) # random uniform number from 0 to maximum of cdf column (it is not exactly 1.0 due to discrete numbers used in PC)

    r_id = findfirst(x -> x .>= gamma,df.cdf)

    r_id > 1 ? Pᵢ₋₁ = df[r_id-1, 8] : Pᵢ₋₁ = 0

    df_Pᵢ= df[r_id, :]

    a = 0.5* df_Pᵢ.a * thickness
    b = df_Pᵢ.b * thickness
    c = Pᵢ₋₁ - gamma - a*df_Pᵢ.minE^2 - b*df_Pᵢ.minE

    x1, x2 = solve_quadratic(a, b, c)

    if (df_Pᵢ.minE < x1) && (df_Pᵢ.maxE > x1)
        return rand(Uniform(df_Pᵢ.E1, df_Pᵢ.E1+thickness)), x1  #E1 is sampled as uniform number between steps,
                                                                #E2 is sampled by CDF
    elseif (df_Pᵢ.minE < x2) && (df_Pᵢ.maxE > x2)
        return rand(Uniform(df_Pᵢ.E1, df_Pᵢ.E1+thickness)), x2
    else
        sample_energies(df)
    end

end

function solve_quadratic(a, b, c)
    d  = sqrt(b^2 - 4*a*c)
    return (-b - d) / (2*a), (-b + d) / (2*a)
end


text:
"""
```
    SubpixelNonmaximaSuppression <: AbstractEdgeThinningAlgorithm
    SubpixelNonmaximaSuppression(; threshold::Union{Number, Percentile} = Percentile(20))

    f = SubpixelNonmaximaSuppression()
    f(out₁::AbstractArray, out₂::Matrix{<:AbstractArray}, mag::AbstractArray, g₁::AbstractArray, g₂::AbstractArray, f::SubpixelNonmaximaSuppression)
```

Isolates local maxima of gradient magnitude `mag` along the local gradient
direction to subpixel precision.  The arguments `g₁` and `g₂` represent the
gradient in the first spatial dimension (y), and the second spatial dimension
(x), respectively.

The integer components of the local maxima correspond to non-zero row and column
entries `out₁`. The accompanying subpixel offsets  are stored in a 2-D array
`out₂` as length-2 vectors. One can recover the  subpixel coordinates by adding
the subpixel offsets to the integer components.

# Details

TODO

# Example

```julia

using TestImages, FileIO, ImageView, ImageEdgeDetection, ImageFiltering

img =  testimage("mandril_gray")

# Gradient in the first and second spatial dimension
g₁, g₂ = imgradients(img, KernelFactors.scharr)

# Gradient magnitude
mag = hypot.(g₁, g₂)

nms = zeros(eltype(mag), axes(mag))
subpixel_offsets = zeros(SVector{2,Float64}, axes(mag))

# Instantiate the NonmaximaSuppression functor.
f = SubpixelNonmaximaSuppression()

# Suppress the non-maximal gradient magnitudes and store the result in `nms`.
f(nms, subpixel_offsets, mag, g₁, g₂)

imshow(img)
imshow(mag)
imshow(nms)
```

# References
1. J. Canny, "A Computational Approach to Edge Detection," in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. PAMI-8, no. 6, pp. 679-698, Nov. 1986, doi: 10.1109/TPAMI.1986.4767851.
2. F. Devernay, "A non-maxima suppression method for edge detection with sub-pixel accuracy", Tech. Report RR-2724, INRIA, 1995.
"""
@with_kw struct SubpixelNonmaximaSuppression{T <: Union{Real,AbstractGray, Percentile}} <: AbstractEdgeThinningAlgorithm
    threshold::T = Percentile(20)
end

function (f::SubpixelNonmaximaSuppression)(out₁::AbstractArray, out₂::AbstractArray{<:StaticVector}, mag::AbstractArray, g₁::AbstractArray, g₂::AbstractArray)
    @unpack threshold = f

    low_threshold =  typeof(threshold) <: Percentile ? StatsBase.percentile(vec(mag), threshold.p) : threshold

    # Isolate local maxima of gradient magnitude by “non-maximum suppression”
    # along the local gradient direction.
    suppress_subpixel_non_maxima!(out₁, out₂, mag, g₁, g₂, low_threshold)

    return out₁, out₂
end

"""
```
    suppress_subpixel_non_maxima!(out₁::AbstractArray, out₂::AbstractArray{<:StaticVector}, mag::AbstractArray, gx::AbstractArray, gy::AbstractArray, threshold::Number)
```

Isolates local maxima of gradient magnitude by “non-maximum suppression” along the local gradient direction.

"""
function suppress_subpixel_non_maxima!(out₁::AbstractArray, out₂::AbstractArray{<:StaticVector},  mag::AbstractArray, g₁::AbstractArray, g₂::AbstractArray, threshold::Number)
    itp = interpolate(mag, BSpline(Linear()))
    etp = extrapolate(itp, Flat())
    m₁ = zero(eltype(mag))
    m₂ = zero(eltype(mag))
    @inbounds for i in CartesianIndices(mag)
        r, c = i.I
        d₁ = gray(g₁[i])
        d₂ = gray(g₂[i])
        mc = mag[i]
        if mc < threshold || mc == 0 || isnan(mc)
            out₁[r,c] = zero(eltype(mag))
        else
            # Ensure the vector 𝐝 = [d₁, d₂] has unit norm.
            d₁ = d₁ / mc
            d₂ = d₂ / mc
            m₁ = etp(r + d₁, c + d₂)
            m₂ = etp(r - d₁, c - d₂)
            if ((m₁ <= mc) && (m₂ < mc)) || ((m₁ < mc) && (m₂ <= mc))
                out₁[r,c] = mag[r,c]
                # To obtain subpixel precision we fit a parabola through the
                # three points:
                # (x₁ = -1, y₁ = m₂) , (x₂ = 0, y₂ = mc), (x₃ = 1, y₃ = m₂)
                # and then solve for the vertex of the parabola.
                λ = (m₂ - m₁) / (2*(m₂ - 2*mc + m₁))
                out₂[r,c] = SVector(λ*d₁, λ*d₂)
            else
                out₁[r,c] = zero(eltype(mag))
                out₂[r,c] = zero(eltype(out₂))
            end
        end
    end
    return nothing
end


text:
# DirichletSpaces

export ChebyshevDirichlet

#TODO: Type of Interval
immutable ChebyshevDirichlet{left,right,D} <: PolynomialSpace{D}
    domain::D
    ChebyshevDirichlet(d)=new(d)
    ChebyshevDirichlet()=new(D())
end

for TYP in (:Number,:AbstractArray,:Fun)
    @eval evaluate(f::AbstractVector,S::ChebyshevDirichlet,x::$TYP) =
        evaluate(Fun(Fun(f,S),canonicalspace(S)),x)
end
@compat (::Type{ChebyshevDirichlet{l,r}}){l,r}() =
    ChebyshevDirichlet{l,r,Interval{Float64}}()
@compat (::Type{ChebyshevDirichlet{l,r}}){l,r}(d::Domain) =
    ChebyshevDirichlet{l,r,typeof(d)}(d)

spacescompatible{l,r,D}(a::ChebyshevDirichlet{l,r,D},b::ChebyshevDirichlet{l,r,D}) =
    domainscompatible(a,b)

ChebyshevDirichlet() = ChebyshevDirichlet{1,1,Interval{Float64}}()
ZeroChebyshevDirichlet(d) =
    ChebyshevDirichlet{1,1,Interval{Float64}}(d)|(3:∞)
ZeroChebyshevDirichlet() =
    ChebyshevDirichlet{1,1,Interval{Float64}}()|(3:∞)

canonicalspace(S::ChebyshevDirichlet) = Chebyshev(domain(S))

setdomain{l,r}(S::ChebyshevDirichlet{l,r},d::Domain) = ChebyshevDirichlet{l,r}(d)

# These are used to make sure Chebyshev comes first

Base.isless(a::Chebyshev,b::ChebyshevDirichlet)=true
<(a::Chebyshev,b::ChebyshevDirichlet)=true
<=(a::Chebyshev,b::ChebyshevDirichlet)=true
>(a::Chebyshev,b::ChebyshevDirichlet)=false
>=(a::Chebyshev,b::ChebyshevDirichlet)=false

Base.isless(a::ChebyshevDirichlet,b::Chebyshev)=false
<(a::ChebyshevDirichlet,b::Chebyshev)=false
<=(a::ChebyshevDirichlet,b::Chebyshev)=false
>(a::ChebyshevDirichlet,b::Chebyshev)=true
>=(a::ChebyshevDirichlet,b::Chebyshev)=true

## coefficients

# converts to f_0 T_0 + f_1 T_1 +  \sum f_k (T_k - T_{k-2})

function dirichlettransform!(w::Vector)
    for k=length(w)-2:-1:1
        @inbounds w[k] += w[k+2]
    end

    w
end

function idirichlettransform!(w::Vector)
    for k=3:length(w)
        @inbounds w[k-2]-= w[k]
    end

    w
end

# converts to f_0 T_0 + \sum f_k (T_k ± T_{k-1})

function idirichlettransform!(s::Bool,w::Vector)
    for k=2:length(w)
        @inbounds w[k-1]+= (s?-1:1)*w[k]
    end

    w
end

function dirichlettransform!(s::Bool,w::Vector)
    for k=length(w)-1:-1:1
        @inbounds w[k] += (s?1:-1)*w[k+1]
    end

    w
end

coefficients(v::Vector,::Chebyshev,::ChebyshevDirichlet{1,1})=dirichlettransform!(copy(v))
coefficients(v::Vector,::Chebyshev,::ChebyshevDirichlet{0,1})=dirichlettransform!(true,copy(v))
coefficients(v::Vector,::Chebyshev,::ChebyshevDirichlet{1,0})=dirichlettransform!(false,copy(v))

coefficients(v::Vector,::ChebyshevDirichlet{1,1},::Chebyshev)=idirichlettransform!(copy(v))
coefficients(v::Vector,::ChebyshevDirichlet{0,1},::Chebyshev)=idirichlettransform!(true,copy(v))
coefficients(v::Vector,::ChebyshevDirichlet{1,0},::Chebyshev)=idirichlettransform!(false,copy(v))

## Dirichlet Conversion operators

Conversion(D::ChebyshevDirichlet,C::Chebyshev)=ConcreteConversion(D,C)

getindex{D,CC<:Chebyshev,T}(C::ConcreteConversion{ChebyshevDirichlet{1,0,D},CC,T},k::Integer,j::Integer) =
    j==k || j==k+1 ? one(T) : zero(T)

getindex{D,CC<:Chebyshev,T}(C::ConcreteConversion{ChebyshevDirichlet{0,1,D},CC,T},k::Integer,j::Integer) =
    j==k ? one(T) : ( j==k+1? -one(T) : zero(eltype(C)))

getindex{D,CC<:Chebyshev,T}(C::ConcreteConversion{ChebyshevDirichlet{1,1,D},CC,T},k::Integer,j::Integer) =
    j==k ? one(T) : ( j==k+2? -one(T) : zero(eltype(C)))

function getindex{D,CC<:Chebyshev,T}(C::ConcreteConversion{ChebyshevDirichlet{2,2,D},CC,T},k::Integer,j::Integer)
    if j==k
        one(T)
    elseif j==k+4
        2one(T)*(k+1)/k-1
    elseif k≥3 && j==k+2
        -2one(T)*(k-1)/(k-2)
    else
        zero(T)
    end
end

bandinds{D,C<:Chebyshev}(::ConcreteConversion{ChebyshevDirichlet{1,0,D},C})=0,1
bandinds{D,C<:Chebyshev}(::ConcreteConversion{ChebyshevDirichlet{0,1,D},C})=0,1
bandinds{D,C<:Chebyshev}(::ConcreteConversion{ChebyshevDirichlet{1,1,D},C})=0,2
bandinds{D,C<:Chebyshev}(::ConcreteConversion{ChebyshevDirichlet{2,2,D},C})=0,4

conversion_rule(b::ChebyshevDirichlet,a::Chebyshev)=b

# return the space that has banded Conversion to the other
# function conversion_rule(a::ChebyshevDirichlet,b::Ultraspherical)
#     @assert domainscompatible(a,b)
#
#     a
# end

## Evaluation Functional

bandinds{D}(B::ConcreteEvaluation{ChebyshevDirichlet{1,0,D},Bool}) = 0,B.x?∞:0
bandinds{D}(B::ConcreteEvaluation{ChebyshevDirichlet{0,1,D},Bool}) = 0,B.x?0:∞
bandinds{D}(B::ConcreteEvaluation{ChebyshevDirichlet{1,1,D},Bool}) = 0,B.x?0:1

function getindex{D}(B::ConcreteEvaluation{ChebyshevDirichlet{1,0,D},Bool},kr::Range)
    d = domain(B)

    if B.x == false && B.order == 0
        Float64[k==1?1.0:0.0 for k=kr]
    elseif B.x == true && B.order == 0
        Float64[k==1?1.0:2.0 for k=kr]
    else
        (Evaluation(d,B.x,B.order)*Conversion(domainspace(B)))[kr]
    end
end

function getindex{D}(B::ConcreteEvaluation{ChebyshevDirichlet{0,1,D},Bool},kr::Range)
    S = Space(domain(B))

    if B.x == true && B.order == 0
        Float64[k==1?1.0:0.0 for k=kr]
    elseif B.x == false && B.order == 0
        Float64[k==1?1.0:-(-1)^k*2.0 for k=kr]
    else
        (Evaluation(S,B.x,B.order)*Conversion(domainspace(B),S))[kr]
    end
end

function getindex{D}(B::ConcreteEvaluation{ChebyshevDirichlet{1,1,D},Bool},kr::Range)
    S = Space(domain(B))

    if B.x == false && B.order == 0
        Float64[k==1?1.0:(k==2?-1.0:0.0) for k=kr]
    elseif B.x == true && B.order == 0
        Float64[k==1||k==2?1.0:0.0 for k=kr]
    else
        getindex(Evaluation(S,B.x,B.order)*Conversion(domainspace(B),S),kr)
    end
end

function Evaluation(sp::ChebyshevDirichlet,x::Float64,ord::Integer)
    S=Space(domain(sp))
    EvaluationWrapper(sp,x,ord,Evaluation(S,x,ord)*Conversion(sp,S))
end


text:
############################################################################################
"""
$(SIGNATURES)
Construct a Tempering tuning struct with the appropriate type for model parameter.

# Examples
```julia
```

"""
function construct(
    tempering::BaytesCore.IterationTempering,
    model::ModelWrapper,
    chains::Integer
)
    MType = model.info.reconstruct.default.output
    @unpack L, k, x₀ = tempering.parameter
    ## Adjust types for parameter and initial value
    parameter = BaytesCore.TemperingParameter(MType(L), MType(k), MType(x₀))
    val₀ = BaytesCore.ValueHolder(MType(tempering.initial.current))
    ## Return correct Tempering struct
    tempertune = BaytesCore.IterationTempering(tempering.adaption, val₀, parameter)
    return tempertune, val₀.current
end
function construct(
    tempering::BaytesCore.JointTempering,
    model::ModelWrapper,
    chains::Integer
)
    tempertune = BaytesCore.JointTempering(
        model.info.reconstruct.default.output,
        tempering.adaption,
        tempering.val.current,
        tempering.ESSTarget.current,
        chains
    )
    return tempertune, tempertune.val.current
end

############################################################################################
"""
$(SIGNATURES)
Construct separate models and data tuning container for each chain.

# Examples
```julia
```

"""
function construct(model::ModelWrapper, datatune::DataTune, chains::Integer, args...)
    ## Initialize separate models that work with each chain
    modelᵛ = map(
        iter -> ModelWrapper(deepcopy(model.val), model.info, model.id), Base.OneTo(chains)
    )
    datatuneᵛ = map(iter -> deepcopy(datatune), Base.OneTo(chains))
    return modelᵛ, datatuneᵛ
end
function construct(model::ModelWrapper, datatune::DataTune, chains::Integer, smc::S) where {S<:Union{SMC, SMCConstructor}}
    return model, datatune
end

############################################################################################
"""
$(SIGNATURES)
Construct `chains` MCMC chains and initiate all algorithms stated in `args` separately. A separate model for each chain is provided to avoid pointer issues.
If `args` is a single `SMC` algorithm, chains will not be separately allocated but instead used within the algorithm.
Note that smc still works as intended if used alongside other mcmc sampler in `args`.

# Examples
```julia
```

"""
function construct(
    _rng::Random.AbstractRNG,
    model::ModelWrapper,
    data::D,
    default::BaytesCore.SampleDefault,
    tempering::BaytesCore.TemperingMethod,
    datatune::DataTune,
    chains::Integer,
    args...;
) where {D}
    ## Initiate Tempering Tune with correct model parameter type.
    temperingtune, temperature₀ = construct(tempering, model, default.chains)
    ## Initialize separate models and data tunes that work with each chain
    modelᵛ, datatuneᵛ = construct(model, datatune, chains, args...)
    ## Initialize sampler combination of first chain to deduce algorithm types
    algorithms = map(
        algorithm -> algorithm(
            _rng,
            modelᵛ[begin],
            BaytesCore.adjust(datatuneᵛ[begin], data),
            temperature₀,
            default
        ),
        args,
    )
    algorithmsᵛ = Vector{typeof(algorithms)}(undef, chains)
    algorithmsᵛ[1] = algorithms
    ## Initialize algorithms from args... chains-1 times
    Threads.@threads for Nchain in 2:chains
        algorithmsᵛ[Nchain] = map(
            algorithm -> algorithm(
                _rng, modelᵛ[Nchain],
                BaytesCore.adjust(datatuneᵛ[Nchain], data),
                temperature₀,
                default
            ),
            args,
        )
    end
    return modelᵛ, algorithmsᵛ, temperingtune, datatuneᵛ
end
function construct(
    _rng::Random.AbstractRNG,
    model::ModelWrapper,
    data::D,
    default::BaytesCore.SampleDefault,
    tempering::BaytesCore.TemperingMethod,
    datatune::DataTune,
    chains::Integer,
    smc::SMCConstructor;
) where {D}
    ## Assign tempering struct with correct temperature type
    temperingtune, temperature₀ = construct(tempering, model, default.chains)
    ## Assign model for chains
    modelᵛ, datatuneᵛ = construct(model, datatune, chains, smc)
    ## Assign algorithm
    algorithmsᵛ = smc(_rng, model, BaytesCore.adjust(datatune, data), temperature₀, default)
    return modelᵛ, algorithmsᵛ, temperingtune, datatuneᵛ
end

############################################################################################
#export
export construct


text:
include("main_Approx_V_f_Pred.jl")

io = open("V_f.txt", "w");
plot_on = 1

##

for Δt_replan_update in [1.0, 2.0, 5.0, 10.0, 20.0]
    for Approx_mode in 0:2
        s_guidance.Approx_mode = Approx_mode
        s_guidance.n           = 1.5

## Series 1: Initial Planning Only Simulation
        Δt_plan_update = 300.0
        df_plansim, t_plansim, x_plansim, h_plansim, γ_plansim, V_plansim, A_L_plansim, R_plansim, c_plansim =
main(x₀, h₀, γ₀, V₀, Δt_plan_update, s_guidance)

## Series 2: Online Replanning Update Simulation
        Δt_plan_update = Δt_replan_update
        df_replansim, t_replansim, x_replansim, h_replansim, γ_replansim, V_replansim, A_L_replansim, R_replansim, c_replansim =
main(x₀, h₀, γ₀, V₀, Δt_plan_update, s_guidance)

## Series 3: Initial Planner Output
x_plan, h_plan, γ_plan, V_plan, A_L_plan, R_plan, c_plan = plan_path(x₀, h₀, γ₀, V₀, s_guidance)

## Plotting
        if plot_on == 1
            label_string_xplot = ["PlanSim" "RePlanSim" "Plan" ]
            label_string_tplot = ["PlanSim" "RePlanSim"]

            f_x_h   = fig_print([x_plansim, x_replansim, x_plan], [h_plansim, h_replansim, h_plan],         "x_h_$(round(Int, Δt_plan_update))_$(Approx_mode)", label_string_xplot, "\$x ~[\\textrm{m}]\$", "\$h ~[\\textrm{m}]\$"; save_file = 0, lgnd_val = :bottomleft) # ar_val = :equal
            f_x_γ   = fig_print([x_plansim, x_replansim, x_plan], [γ_plansim, γ_replansim, γ_plan],     "x_gamma_$(round(Int, Δt_plan_update))_$(Approx_mode)", label_string_xplot, "\$x ~[\\textrm{m}]\$", "\$\\gamma ~[\\textrm{deg}]\$")
            f_x_V   = fig_print([x_plansim, x_replansim, x_plan], [V_plansim, V_replansim, V_plan],         "x_V_$(round(Int, Δt_plan_update))_$(Approx_mode)", label_string_xplot, "\$x ~[\\textrm{m}]\$", "\$V ~[\\textrm{m}/\\textrm{s}]\$")
            f_x_A_L = fig_print([x_plansim, x_replansim, x_plan], [A_L_plansim, A_L_replansim, A_L_plan], "x_A_L_$(round(Int, Δt_plan_update))_$(Approx_mode)", label_string_xplot, "\$x ~[\\textrm{m}]\$", "\$A_{L} ~[\\textrm{m}/\\textrm{s}^{2}]\$")

            f_t_x   = fig_print([t_plansim, t_replansim], [x_plansim, x_replansim],         "t_x_$(round(Int, Δt_plan_update))_$(Approx_mode)", label_string_tplot, "\$t ~[\\textrm{s}]\$", "\$x ~[\\textrm{m}]\$")
            f_t_h   = fig_print([t_plansim, t_replansim], [h_plansim, h_replansim],         "t_h_$(round(Int, Δt_plan_update))_$(Approx_mode)", label_string_tplot, "\$t ~[\\textrm{s}]\$", "\$h ~[\\textrm{m}]\$")
            f_t_γ   = fig_print([t_plansim, t_replansim], [γ_plansim, γ_replansim],     "t_gamma_$(round(Int, Δt_plan_update))_$(Approx_mode)", label_string_tplot, "\$t ~[\\textrm{s}]\$", "\$\\gamma ~[\\textrm{deg}]\$")
            f_t_V   = fig_print([t_plansim, t_replansim], [V_plansim, V_replansim],         "t_V_$(round(Int, Δt_plan_update))_$(Approx_mode)", label_string_tplot, "\$t ~[\\textrm{s}]\$", "\$V ~[\\textrm{m}/\\textrm{s}]\$")
            f_t_A_L = fig_print([t_plansim, t_replansim], [A_L_plansim, A_L_replansim],   "t_A_L_$(round(Int, Δt_plan_update))_$(Approx_mode)", label_string_tplot, "\$t ~[\\textrm{s}]\$", "\$A_{L} ~[\\textrm{m}/\\textrm{s}^{2}]\$")
            f_t_R   = fig_print([t_plansim, t_replansim], [R_plansim, R_replansim],         "t_R_$(round(Int, Δt_plan_update))_$(Approx_mode)", label_string_tplot, "\$t ~[\\textrm{s}]\$", "\$R\$")

            f_t_c = plot()
            fig_print(t_plansim,   c_plansim,     [], "PlanSim",   "\$t ~[\\textrm{s}]\$", "\$\\mathbf{c}\$", f_t_c; save_file = 0)
            fig_print(t_replansim, c_replansim,   [], "RePlanSim", "\$t ~[\\textrm{s}]\$", "\$\\mathbf{c}\$", f_t_c; save_file = 0)
            f_t_c = fig_print([], [], "t_c_$(round(Int, Δt_plan_update))_$(Approx_mode)", [], [], [], f_t_c)
        end

        write(io, "Δt_replan_update: $(round(Int, Δt_plan_update)) / Approx_mode: $(Approx_mode) \n")
        write(io, "x_f: plansim = $(x_plansim[end]), replansim = $(x_replansim[end]), plan = $(x_plan[end])\n")
        write(io, "h_f: plansim = $(h_plansim[end]), replansim = $(h_replansim[end]), plan = $(h_plan[end])\n")
        write(io, "γ_f: plansim = $(γ_plansim[end]), replansim = $(γ_replansim[end]), plan = $(γ_plan[end])\n")
        write(io, "V_f: plansim = $(V_plansim[end]), replansim = $(V_replansim[end]), plan = $(V_plan[end])\n\n\n")
    end
end
close(io)

## -----------------------------------------------------------------------------------------
# Test
# -----------------------------------------------------------------------------------------
# ##
# status, J_opt, c_opt = Path_Generator(s_guidance, param.R, x₀, h₀, γ₀)
# V̂_f = Speed_Predictor(s_guidance, x₀, V₀, c_opt)

# ##
# R_list = -5:1:5
# V̂_f_list = R_list |> Map(R -> begin
#     status, J_opt, c_opt = Path_Generator(s_guidance, R, x₀, h₀, γ₀)
#     V̂_f = Speed_Predictor(s_guidance, x₀, V₀, c_opt)
# end) |> tcollect

##
# s₀     = State(env_model)(x₀, h₀, γ₀, V₀)
# R₀     = 1.0
# c₀     = zeros(p)
# param₀ = Params(env_model)(R₀, c₀)

# s_guidance.Approx_mode = 0
# param_exact  = IASCG_Planner(s_guidance, s₀, param₀)
# s_guidance.Approx_mode = 1
# s_guidance.n = 1.5
# param_approx = IASCG_Planner(s_guidance, s₀, param₀)

# ##
# x = x₀:1:x_f

# order = p-1:-1:0
# h_exact = x |> Map(x -> h_ref * ((x / x_ref).^order)' * param_exact.c) |> collect
# h_exact=hcat(h_exact...)'
# h_approx = x |> Map(x -> h_ref * ((x / x_ref).^order)' * param_approx.c) |> collect
# h_approx=hcat(h_approx...)'

# f_x_h = fig_print(x, [h_exact h_approx], "x_h_plan_comparison", ["exact" "approx"], "\$x ~[\\textrm{m}]\$", "\$h ~[\\textrm{m}]\$"; save_file=1)
# # fig_print([], [], "x_h", [], [], [], f_x_h; ar_val=:equal, lfs_val=6)


text:
"""
    all_bound(atom)

Returns `true` if all orbitals in `atom` are bound orbitals.
"""
all_bound(atom::Atom) = all(isbound.(atom.orbitals))


text:
#! /u1/vian294/julia/julia-1.1.0/bin/julia -p 8
#PBS -N jHexaG1e3p08
#PBS -l select=1:ncpus=8
#PBS -j oe
#PBS -k oe

# qsub /u1/vian294/METMAT/AutoDiff/rev_paper1/example_3DSpring_r1.jl

using Printf, WriteVTK, AbaqusReader, JLD, Dates, LinearAlgebra
@everywhere cd("/u1/vian294/METMAT/AutoDiff/rev_paper1")
@show pwd()

@everywhere using AD4SM

@show Elements.setp(8)

mean(x)   = sum(x)/length(x)

sMeshFile = "3DSpringHexaj.inp"
mat       = Materials.NeoHooke(10, 1e3)
bisinc    = true
sVTKpath  = "./vtk_files/"
bVTKall   = true
Δz        = 250
sPosFix   = "_1e3w250"
N         = 300
LF        = vcat(range(0.0, 0.7, length=N÷3),
                 range(0.7, 1.0, length=2N÷3+1))
nSteps    = length(LF)

sFileName = splitext(sMeshFile)[1]*sPosFix
mymodel   = AbaqusReader.abaqus_read_mesh(sMeshFile)
nodes     = [mymodel["nodes"][ii] for ii in 1:mymodel["nodes"].count]
el_nodes  = [item[2]              for item in mymodel["elements"]]
id_b      = mymodel["node_sets"]["BTM"]
id_t      = mymodel["node_sets"]["TOP"]
n_t, n_b  = length(id_t), length(id_b)

points    = hcat(nodes...)
cells     = [if length(item)==4
               MeshCell(VTKCellTypes.VTK_TETRA, item)
             else
               MeshCell(VTKCellTypes.VTK_HEXAHEDRON, item)
             end for item in el_nodes]
elems     = [if length(item)==4
               Elements.Tet04(item, nodes[item], mat=mat)
             else
               Elements.Hex08(item, nodes[item], mat=mat)
             end for item  in el_nodes]

@show nNodes, nElems   = length(nodes), length(elems)
@printf("starting %s \n\n", sFileName)

ifree = trues(3, nNodes)
unew  = zeros(3, nNodes)
λnew  = zeros(6)
idxes = LinearIndices(unew)
fnew  = zeros(3, nNodes)
allus = []
t0    = Base.time_ns()

for (ii,LF) in enumerate(LF)
  global unew, λnew, fnew
  w0    = LF*Δz
  @printf("doing step %3i/%i, LF = %.4f, w0 = %.3f",
          ii, nSteps, LF, w0); flush(stdout)

  eqns  = [Elements.ConstEq(x->sum(x),          idxes[1,id_t][:], adiff.D1),
           Elements.ConstEq(x->sum(x),          idxes[2,id_t][:], adiff.D1),
           Elements.ConstEq(x->sum(x)/n_t-w0/2, idxes[3,id_t][:], adiff.D1),
           Elements.ConstEq(x->sum(x),          idxes[1,id_b][:], adiff.D1),
           Elements.ConstEq(x->sum(x),          idxes[2,id_b][:], adiff.D1),
           Elements.ConstEq(x->sum(x)/n_b+w0/2, idxes[3,id_b][:], adiff.D1) ]

  lastu = copy(unew)
  lastλ = copy(λnew)
  fnew  = zeros(3, nNodes)
  T     = @elapsed (bfailed, normr, iter) =
  Elements.solvestep!(elems, lastu, unew, ifree,
                      eqns      = eqns,
                      λ         = λnew,
                      fe        = fnew,
                      dTolu     = 1e-4,
                      dTole     = 1e-3,
                      dNoise    = 1e-9,
                      bprogress = false,
                      bpredict  = false,
                      becho     = true)
  if bfailed
    @printf("!! failed at LF: %.3f, with normr: %.3e\n\n", LF, normr)
    unew = lastu
    λnew = lastλ
    break
  else
    push!(allus, (copy(unew), copy(fnew), copy(λnew)))
    @printf("step %2i done in %2i iter, after %.2f sec.\n\n", ii, iter, T)
  end
  flush(stdout)
end
@printf("completed in %s\n",(Base.time_ns()-t0)÷1e9|>
        Dates.Second|>Dates.CompoundPeriod|>Dates.canonicalize)
flush(stdout)

Δu_tot = [mean(item[1][3,id_b]) for item in allus]
rf_tot = [item[3][3]            for item in allus]

JLD.save(sFileName*".jld", "nodes", nodes, "elems", elems, "allus", allus,
         "points", points, "cells", cells, "LF", LF, "mat", mat,
         "Nsteps", length(LF), "id_b", id_b, "id_t", id_t,
         "rf_tot", rf_tot, "Δu_tot", Δu_tot, "sFileName", sFileName)
@printf("results written to %s\n", sFileName)

paraview_collection(sVTKpath*sFileName) do pvd
  for (ii, item) in enumerate(allus)
    (u0,rf) = item[1:2]
    F   = Elements.getinfo(elems,u0,info=:F)
    LE  = Elements.getinfo(elems,u0,info=:LE)
    J   = [det(F) for F in F]
    E   = [0.5*(transpose(F)*F-I) for F in F]
    σ   = [Materials.getinfo(F,mat,info=:σ)        for F in F]  # Cauchy stress
    S   = [Materials.getinfo(F,mat,info=:S)        for F in F]  # 2nd PK
    Ii  = [Materials.getInvariants(transpose(F)F)  for F in F]  # 2nd PK
    σP  = [sort(eigvals(σ))  for σ  in σ]
    σVM = [sqrt((x[1]-x[2])^2+(x[2]-x[3])^2+(x[3]-x[1])^2)/sqrt(2)  for x in σP]

    vtkobj  = vtk_grid(@sprintf("%s_%03i", sVTKpath*sFileName, ii), points+u0, cells)
    vtk_point_data(vtkobj, (u0[1,:],u0[2,:],u0[3,:]),                   "u")
    vtk_point_data(vtkobj, (rf[1,:],rf[2,:],rf[3,:]),                   "rf")
    vtk_cell_data(vtkobj,  tuple([getindex.(F,ii)  for ii in 1:9]...),  "F")
    vtk_cell_data(vtkobj,  tuple([getindex.(LE,ii)  for ii in 1:9]...), "LE")
    vtk_cell_data(vtkobj,  J,   "J")
    vtk_cell_data(vtkobj,  tuple([getindex.(σ,ii)   for ii in 1:9]...), "\$\\sigma\$")
    vtk_cell_data(vtkobj,  tuple([getindex.(σP,ii)  for ii in 1:3]...), "\$\\sigma_p\$")
    vtk_cell_data(vtkobj,  tuple([getindex.(Ii,ii)  for ii in 1:3]...), "I")
    vtk_cell_data(vtkobj,  σVM, "\$\\sigma_{VM}\$")

    collection_add_timestep(pvd, vtkobj, ii)
    vtk_save(vtkobj)
  end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    Ref{T}

An object that safely references data of type `T`. This type is guaranteed to point to
valid, Julia-allocated memory of the correct type. The underlying data is protected from
freeing by the garbage collector as long as the `Ref` itself is referenced.

In Julia, `Ref` objects are dereferenced (loaded or stored) with `[]`.

Creation of a `Ref` to a value `x` of type `T` is usually written `Ref(x)`.
Additionally, for creating interior pointers to containers (such as Array or Ptr),
it can be written `Ref(a, i)` for creating a reference to the `i`-th element of `a`.

When passed as a `ccall` argument (either as a `Ptr` or `Ref` type), a `Ref` object will be
converted to a native pointer to the data it references.

There is no invalid (NULL) `Ref` in Julia, but a `C_NULL` instance of `Ptr` can be passed to
a `ccall` Ref argument.

# Use in broadcasting

Broadcasting with `Ref(x)` treats `x` as a scalar:
```jldoctest
julia> isa.(Ref([1,2,3]), [Array, Dict, Int])
3-element BitArray{1}:
 1
 0
 0
```
"""
Ref

# C NUL-terminated string pointers; these can be used in ccall
# instead of Ptr{Cchar} and Ptr{Cwchar_t}, respectively, to enforce
# a check for embedded NUL chars in the string (to avoid silent truncation).
if Int === Int64
    primitive type Cstring  64 end
    primitive type Cwstring 64 end
else
    primitive type Cstring  32 end
    primitive type Cwstring 32 end
end

### General Methods for Ref{T} type

eltype(x::Type{<:Ref{T}}) where {T} = @isdefined(T) ? T : Any
convert(::Type{Ref{T}}, x::Ref{T}) where {T} = x
size(x::Ref) = ()
axes(x::Ref) = ()
length(x::Ref) = 1
ndims(x::Ref) = 0
ndims(::Type{<:Ref}) = 0
iterate(r::Ref) = (r[], nothing)
iterate(r::Ref, s) = nothing
IteratorSize(::Type{<:Ref}) = HasShape{0}()

# create Ref objects for general object conversion
unsafe_convert(::Type{Ref{T}}, x::Ref{T}) where {T} = unsafe_convert(Ptr{T}, x)
unsafe_convert(::Type{Ref{T}}, x) where {T} = unsafe_convert(Ptr{T}, x)

convert(::Type{Ref{T}}, x) where {T} = RefValue{T}(x)

### Methods for a Ref object that is backed by an array at index i
struct RefArray{T,A<:AbstractArray{T},R} <: Ref{T}
    x::A
    i::Int
    roots::R # should be either ::Nothing or ::Any
    RefArray{T,A,R}(x,i,roots=nothing) where {T,A<:AbstractArray{T},R} = new(x,i,roots)
end
RefArray(x::AbstractArray{T}, i::Int, roots::Any) where {T} = RefArray{T,typeof(x),Any}(x, i, roots)
RefArray(x::AbstractArray{T}, i::Int=1, roots::Nothing=nothing) where {T} = RefArray{T,typeof(x),Nothing}(x, i, nothing)
convert(::Type{Ref{T}}, x::AbstractArray{T}) where {T} = RefArray(x, 1)

function unsafe_convert(P::Type{Ptr{T}}, b::RefArray{T}) where T
    if datatype_pointerfree(RefValue{T})
        p = pointer(b.x, b.i)
    elseif isconcretetype(T) && T.mutable
        p = pointer_from_objref(b.x[b.i])
    else
        # see comment on equivalent branch for RefValue
        p = pointerref(Ptr{Ptr{Cvoid}}(pointer(b.x, b.i)), 1, Core.sizeof(Ptr{Cvoid}))
    end
    return convert(P, p)
end
function unsafe_convert(P::Type{Ptr{Any}}, b::RefArray{Any})
    return convert(P, pointer(b.x, b.i))
end
unsafe_convert(::Type{Ptr{Cvoid}}, b::RefArray{T}) where {T} = convert(Ptr{Cvoid}, unsafe_convert(Ptr{T}, b))

###
if is_primary_base_module
    Ref(x::Any) = RefValue(x)
    Ref{T}() where {T} = RefValue{T}() # Ref{T}()
    Ref{T}(x) where {T} = RefValue{T}(x) # Ref{T}(x)

    Ref(x::Ref, i::Integer) = (i != 1 && error("Ref only has one element"); x)
    Ref(x::Ptr{T}, i::Integer) where {T} = x + (i - 1) * Core.sizeof(T)

    # convert Arrays to pointer arrays for ccall
    function Ref{P}(a::Array{<:Union{Ptr,Cwstring,Cstring}}) where P<:Union{Ptr,Cwstring,Cstring}
        return RefArray(a) # effectively a no-op
    end
    function Ref{P}(a::Array{T}) where P<:Union{Ptr,Cwstring,Cstring} where T
        if (!isbitstype(T) && T <: eltype(P))
            # this Array already has the right memory layout for the requested Ref
            return RefArray(a,1,false) # root something, so that this function is type-stable
        else
            ptrs = Vector{P}(undef, length(a)+1)
            roots = Vector{Any}(undef, length(a))
            for i = 1:length(a)
                root = cconvert(P, a[i])
                ptrs[i] = unsafe_convert(P, root)::P
                roots[i] = root
            end
            ptrs[length(a)+1] = C_NULL
            return RefArray(ptrs,1,roots)
        end
    end
    Ref(x::AbstractArray, i::Integer) = RefArray(x, i)
end

cconvert(::Type{Ptr{P}}, a::Array{<:Ptr}) where {P<:Ptr} = a
cconvert(::Type{Ref{P}}, a::Array{<:Ptr}) where {P<:Ptr} = a
cconvert(::Type{Ptr{P}}, a::Array) where {P<:Union{Ptr,Cwstring,Cstring}} = Ref{P}(a)
cconvert(::Type{Ref{P}}, a::Array) where {P<:Union{Ptr,Cwstring,Cstring}} = Ref{P}(a)

###

getindex(b::RefArray) = b.x[b.i]
setindex!(b::RefArray, x) = (b.x[b.i] = x; b)

###

max_stars_repo_name: milesfrain/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
using Roots, EngEconomics

# Given
A1 = 9000
N1 = 4
N2 = 5
i = 0.05

# Find
A2 = A1 * N1 * capitalRecoveryFactor(i, N2)
NPWi(x) = A1 + A1 * seriesPresentAmountFactor(x, N1 - 1) - A2 * seriesPresentAmountFactor(x, N2) * presentWorthFactor(x, 4)
find_zero(NPWi, 0.05)


text:
using ProjAssistant
@quickactivate

## ------------------------------------------------------
@time begin
    import Chemostat_InSilico
    const Dyn = Chemostat_InSilico.Dynamic
    import Chemostat_InSilico.Dynamic:
        SimD2, run_simD2!, check_stst, hist,
        n_ave_conv!, tail_ave, maxent,
        Container, vec!, Vi,
        normalizeP!,
        save_simdat, load_simdat, simdat_file,
        set_status, get_status,
        load_batch, save_batch,
        UNDONE_SIM_STATUS, DEAD_SIM_STATUS,
        EXPLODED_SIM_STATUS, NITERS_SIM_STATUS,
        STST_SIM_STATUS, collect_ts,
        STST, is_steady

    using Plots
    import GR
    !isinteractive() && GR.inline("png")

    using Base.Threads
    using ProgressMeter
    using BenchmarkTools
    using ExtractMacro

end

## ------------------------------------------------------
# globlas
# batch = (; push_frec, Xts, sgts, z_avts, ug_avts, cgD_Xts, Pzts, Pugts)
params = lglob(Dyn, :dyn, :params, :finite_cg)
@extract params: simid Ds ϵs cg

simparams = let
    iter = collect(Iterators.product(Ds, ϵs, cg))
    filter(iter) do tp
        D, ϵ, cg = tp
        status = Dyn.get_status(simid, (;D, ϵ, cg))
        (status != STST_SIM_STATUS) && return false
        return true
    end
end

## ----------------------------------------------------------------------------
# ME_AVZ_EXPECTED_AVUG_BOUNDED
let

    WLOCK = ReentrantLock()
    MEmode = "ME_AVZ_EXPECTED_AVUG_BOUNDED"

    sim_tot = length(simparams)
    @threads for (simi, (D, ϵ, cg)) in collect(enumerate(simparams))
        thid = threadid()

        simparams = (;D, ϵ, cg)

        # ---------------------------------------------------------------
        # Welcome
        @info("At", simid, simparams, thid)

        # ---------------------------------------------------------------
        datfile = Dyn.maxent_file(simid, MEmode, simparams)
        isfile(datfile) && continue

        # ---------------------------------------------------------------
        # Sim globals
        S = Dyn.load_simdat(simid, simparams)
        X = S.X
        cgD_X = S.cg * S.D / S.X
        V = S.V
        S = nothing; GC.gc()

        # ---------------------------------------------------------------
        # maxent
        me_out = maxent(V, D, cgD_X;
            # Top loop
            top_maxiter = 100,

            # GD globals
            gd_maxiter = 1000,
            gdth = 1e-2,
            verbose = true,

            # stst
            stst_th = 1e-2,
            stst_w = 10
        )

        sdat(Dyn, me_out, datfile)

    end # for (simi, (D, ϵ, cg))
end


text:
# Engine collecting cash flows along a market-model simulation for doing pathwise computation of Deltas and vegas using Giles--Glasserman smoking adjoints method
# note only works with displaced LMM.

# The method is intimately connected with log-normal Euler evolution

# We must work with discretely compounding MM account
# To compute a vega means changing the pseudo-square root at each time step
# So for each vega, we have a vector of matrices. So we need a vector of vectors of matrices to compute all the vegas.
# We do the outermost vector by time step and inner one by which vega.
# This implementation is different in that all the linear combinations by the bumps are done as late as possible,
# whereas PathwiseVegasAccountingEngine does them as early as possible.

mutable struct PathwiseVegasOuterAccountingEngine{P <: MarketModelPathwiseMultiProduct, M <: AbstractMarketModel}
  evolver::LogNormalFwdRateEuler
  product::P
  pseudoRootStructure::M
  vegaBumps::Vector{Vector{Matrix{Float64}}}
  numeraires::Vector{Int}
  initialNumeraireValue::Float64
  numberProducts::Int
  numberRates::Int
  numberCashFlowTimes::Int
  numberSteps::Int
  factors::Int
  numberBumps::Int
  numberElementaryVegas::Int
  jacobianComputers::Vector{RatePseudoRootJacobianAllElements}
  doDeflation::Bool
  currentForwards::Vector{Float64}
  lastForwards::Vector{Float64}
  numerairesHeld::Vector{Float64}
  numberCashFlowsThisStep::Vector{Int}
  cashFlowsGenerated::Vector{Vector{MarketModelPathWiseCashFlow}}
  discounters::Vector{MarketModelPathwiseDiscounter}
  V::Vector{Matrix{Float64}}
  LIBORRatios::Matrix{Float64}
  Discounts::Matrix{Float64}
  StepsDiscountsSquared::Matrix{Float64}
  stepsDiscounts::Vector{Float64}
  LIBORRates::Matrix{Float64}
  partials::Matrix{Float64}
  elementaryVegasThisPath::Vector{Vector{Matrix{Float64}}}
  jacobiansThisPaths::Vector{Vector{Matrix{Float64}}}
  deflatorAndDerivatives::Vector{Float64}
  fullDerivatives::Vector{Float64}
  numberCashFlowsThisIndex::Vector{Vector{Int}}
  totalCashFlowsThisIndex::Vector{Matrix{Float64}}
  cashFlowIndicesThisStep::Vector{Vector{Int}}
end

function PathwiseVegasOuterAccountingEngine(evolver::LogNormalFwdRateEuler,
                                            product::MarketModelPathwiseMultiProduct,
                                            pseudoRootStructure::AbstractMarketModel,
                                            vegaBumps::Vector{Vector{Matrix{Float64}}},
                                            initialNumeraireValue::Float64)
  product = clone(product)
  numberProducts = number_of_products(product)
  doDeflation = ~already_deflated(product)
  numerairesHeld = zeros(numberProducts)
  numberCashFlowsThisStep = zeros(Int, numberProducts)
  cashFlowsGenerated = Vector{Vector{MarketModelPathWiseCashFlow}}(undef, numberProducts)
  stepsDiscounts = zeros(pseudoRootStructure.numberOfRates + 1)
  elementaryVegasThisPath = Vector{Vector{Matrix{Float64}}}(undef, numberProducts)
  deflatorAndDerivatives = zeros(pseudoRootStructure.numberOfRates + 1)

  stepsDiscounts[1] = 1.0

  numberRates = pseudoRootStructure.numberOfRates
  numberSteps = pseudoRootStructure.numberOfSteps
  factors = pseudoRootStructure.numberOfFactors
  fullDerivatives = zeros(numberRates)

  evolution = pseudoRootStructure.evolution
  numeraires = money_market_measure(evolution)

  length(vegaBumps) == numberSteps || error("we need precisely one vector of vega bumps for each step")

  numberBumps = length(vegaBumps[1])

  jacobiansThisPathsModel = Vector{Matrix{Float64}}(undef, numberRates)

  @simd for i in eachindex(jacobiansThisPathsModel)
    @inbounds jacobiansThisPathsModel[i] = zeros(numberRates, factors)
  end

  jacobianComputers = Vector{RatePseudoRootJacobianAllElements}(undef, numberSteps)
  jacobiansThisPaths = Vector{Vector{Matrix{Float64}}}(undef, numberSteps)
  @simd for i in eachindex(jacobianComputers)
    @inbounds jacobianComputers[i] = RatePseudoRootJacobianAllElements(pseudoRootStructure.pseudoRoots[i], evolution.firstAliveRate[i],
                            numeraires[i], evolution.rateTaus, pseudoRootStructure.displacements)
    @inbounds jacobiansThisPaths[i] = deepcopy(jacobiansThisPathsModel)
  end

  VModel = zeros(numberSteps + 1, numberRates)

  Discounts = zeros(numberSteps + 1, numberRates + 1)

  @simd for i = 1:numberSteps+1
    @inbounds Discounts[i, 1] = 1.0
  end

  V = Vector{Matrix{Float64}}(undef, numberProducts)

  modelCashFlowIndex = zeros(length(possible_cash_flow_times(product)), numberRates + 1)

  numberCashFlowsThisIndex = Vector{Vector{Int}}(undef, numberProducts)
  V = Vector{Matrix{Float64}}(undef, numberProducts)
  totalCashFlowsThisIndex = Vector{Matrix{Float64}}(undef, numberProducts)

  @inbounds @simd for i in eachindex(numberCashFlowsThisIndex)
    p = max_number_of_cashflows_per_product_per_step(product)
    cashFlowsGenerated[i] = MarketModelPathWiseCashFlow[MarketModelPathWiseCashFlow(numberRates + 1) for i = 1:p]

    numberCashFlowsThisIndex[i] = Vector{Int}(undef, length(possible_cash_flow_times(product)))
    V[i] = copy(VModel)
    totalCashFlowsThisIndex[i] = copy(modelCashFlowIndex)
  end

  LIBORRatios = copy(VModel)
  StepsDiscountsSquared = copy(VModel)
  LIBORRates = copy(VModel)

  cashFlowTimes = possible_cash_flow_times(product)
  numberCashFlowTimes = length(cashFlowTimes)

  rateTimes = product.evolution.rateTimes
  evolutionTimes = product.evolution.evolutionTimes

  discounters = MarketModelPathwiseDiscounter[MarketModelPathwiseDiscounter(cashFlowTimes[j], rateTimes) for j = 1:numberCashFlowTimes]

  ## need to check that we are in money market measure

  # we need to allocate cash-flow times to steps, i.e. what is the last step completed before a flow occurs
  # what we really need for each step, what cash flow time indices to look at

  cashFlowIndicesThisStep = Vector{Int}[Vector{Int}() for i = 1:numberSteps]
  @inbounds @simd for i in eachindex(cashFlowTimes)
    # idx = upper_bound(evolutionTimes, cashFlowTimes[i])
    # if idx != 1
    #   idx -= 1
    # end
    idx = findlast(isequal(cashFlowTimes[i]), evolutionTimes)
    if idx == nothing
      if cashFlowTimes[i] > evolutionTimes[end]
        idx = length(evolutionTimes)
      else
        idx = 1
      end
    end
    push!(cashFlowIndicesThisStep[idx], i)
  end

  partials = zeros(factors, numberRates)

  begin
    modelVegaMatrix = zeros(numberRates, factors)
    @simd for i in eachindex(elementaryVegasThisPath)
      @inbounds elementaryVegasThisPath[i] = Matrix[copy(modelVegaMatrix) for j = 1:numberSteps]
    end
  end

  numberElementaryVegas = numberSteps * numberRates * factors

  return PathwiseVegasOuterAccountingEngine(evolver, product, pseudoRootStructure, vegaBumps, numeraires, initialNumeraireValue, numberProducts, numberRates,
                                            numberCashFlowTimes, numberSteps, factors, numberBumps, numberElementaryVegas, jacobianComputers, doDeflation,
                                            Vector{Float64}(), Vector{Float64}(), numerairesHeld, numberCashFlowsThisStep, cashFlowsGenerated, discounters,
                                            V, LIBORRatios, Discounts, StepsDiscountsSquared, stepsDiscounts, LIBORRates, partials, elementaryVegasThisPath,
                                            jacobiansThisPaths, deflatorAndDerivatives, fullDerivatives, numberCashFlowsThisIndex, totalCashFlowsThisIndex,
                                            cashFlowIndicesThisStep)
end

function single_path_values!(pwEng::PathwiseVegasOuterAccountingEngine, values::Vector{Float64})
  initialForwards = copy(pwEng.pseudoRootStructure.initialRates)
  pwEng.currentForwards = copy(initialForwards)

  @inbounds @simd for i = 1:pwEng.numberProducts
    pwEng.numerairesHeld[i] = 0.0
    for j = 1:pwEng.numberCashFlowTimes
      pwEng.numberCashFlowsThisIndex[i][j] = 0
      for k = 1:pwEng.numberRates + 1
        pwEng.totalCashFlowsThisIndex[i][j, k] = 0.0
      end
    end
    for l = 1:pwEng.numberRates
      for m = 1:pwEng.numberSteps + 1
        pwEng.V[i][m, l] = 0.0
      end
    end
  end

  weight = start_new_path!(pwEng.evolver)
  reset!(pwEng.product)

  isDone = false
  thisStep = 1
  while ~isDone
    thisStep = pwEng.evolver.currentStep
    storeStep = thisStep + 1
    weight *= advance_step!(pwEng.evolver)

    isDone = next_time_step!(pwEng.product, current_state(pwEng.evolver), pwEng.numberCashFlowsThisStep, pwEng.cashFlowsGenerated)

    pwEng.lastForwards = copy(pwEng.currentForwards)

    pwEng.currentForwards = copy(current_state(pwEng.evolver).forwardRates)

    @inbounds @simd for i = 1:pwEng.numberRates
      x = discount_ratio(current_state(pwEng.evolver), i+1, i)
      pwEng.stepsDiscounts[i+1] = x
      pwEng.StepsDiscountsSquared[storeStep, i] = x*x

      pwEng.LIBORRatios[storeStep, i] = pwEng.currentForwards[i] / pwEng.lastForwards[i]
      pwEng.LIBORRates[storeStep, i] = pwEng.currentForwards[i]
      pwEng.Discounts[storeStep, i+1] = discount_ratio(current_state(pwEng.evolver), i+1, 1)
    end

    get_bumps!(pwEng.jacobianComputers[thisStep], pwEng.lastForwards, pwEng.stepsDiscounts, pwEng.currentForwards, brownians_this_step(pwEng.evolver),
              pwEng.jacobiansThisPaths[thisStep])

    # for each product and each cash flow
    @inbounds for i = 1:pwEng.numberProducts, j = 1:pwEng.numberCashFlowsThisStep[i]
      k = pwEng.cashFlowsGenerated[i][j].timeIndex
      pwEng.numberCashFlowsThisIndex[i][k] += 1

      @simd for l = 1:pwEng.numberRates + 1
        pwEng.totalCashFlowsThisIndex[i][k, l] += pwEng.cashFlowsGenerated[i][j].amount[l] * weight
      end
    end

    # if done
    #   break
    # end
  end

  # ok we've gathered cash-flows, still have to do backwards computation
  factors = pwEng.pseudoRootStructure.numberOfFactors
  taus = pwEng.pseudoRootStructure.evolution.rateTaus

  flowsFound = false
  finalStepDone = thisStep

  @inbounds for currentStep = pwEng.numberSteps:-1:1
    stepToUse = min(currentStep, finalStepDone) + 1
    for k in eachindex(pwEng.cashFlowIndicesThisStep[currentStep])
      cashFlowIndex = pwEng.cashFlowIndicesThisStep[currentStep][k]

      # first check to see if anything actually happened before spending time on computing stuff
      noFlows = true
      for l = 1:pwEng.numberProducts
        noFlows = noFlows && (pwEng.numberCashFlowsThisIndex[l][cashFlowIndex] == 0)
        if ~noFlows
          break
        end
      end

      flowsFound = flowsFound || ~noFlows

      if ~noFlows
        if pwEng.doDeflation
          pwEng.deflatorAndDerivatives = get_factors(pwEng.discounters[cashFlowIndex], pwEng.LIBORRates, pwEng.Discounts, stepToUse)
        end

        for j = 1:pwEng.numberProducts
          if pwEng.numberCashFlowsThisIndex[j][cashFlowIndex] > 0
            deflatedCashFlow = pwEng.totalCashFlowsThisIndex[j][cashFlowIndex, 1]
            if pwEng.doDeflation
              deflatedCashFlow *= pwEng.deflatorAndDerivatives[1]
            end

            pwEng.numerairesHeld[j] += deflatedCashFlow

            for i = 2:pwEng.numberRates + 1
              thisDerivative = pwEng.totalCashFlowsThisIndex[j][cashFlowIndex, i]
              if pwEng.doDeflation
                thisDerivative *= pwEng.deflatorAndDerivatives[1]
                thisDerivative += pwEng.totalCashFlowsThisIndex[j][cashFlowIndex, 1] * pwEng.deflatorAndDerivatives[i]
              else
                pwEng.fullDerivatives[i-1] = thisDerivative
              end

              pwEng.V[j][stepToUse, i-1] += thisDerivative
            end # end of for i = 2:numberRates + 1
          end # end of (numberCashFlowsThisIndex[j][cashFlowIndex] > 0)
        end # end of j = 1:numberProducts
      end # end of ~noFlows
    end # end of k in eachindex(cashFlowIndicesThisStep)

    # need to do backwards updating
    if flowsFound
      nextStepToUse = min(currentStep - 1, finalStepDone)
      nextStepIndex = nextStepToUse + 1

      if nextStepIndex != stepToUse # then we need to update V
        thisPseudoRoot = pwEng.pseudoRootStructure.pseudoRoots[currentStep]

        for i = 1:pwEng.numberProducts
          # compute partials
          for f = 1:factors
            libor = pwEng.LIBORRates[stepToUse, pwEng.numberRates]
            V = pwEng.V[i][stepToUse, pwEng.numberRates]
            pseudo = thisPseudoRoot[pwEng.numberRates, f]
            thisPartialTerm = libor * V * pseudo
            pwEng.partials[f, pwEng.numberRates] = thisPartialTerm

            for r = pwEng.numberRates-1:-1:1
              thisPartialTermr = pwEng.LIBORRates[stepToUse, r] * pwEng.V[i][stepToUse, r] * thisPseudoRoot[r, f]
              pwEng.partials[f, r] = pwEng.partials[f, r+1] + thisPartialTermr
            end
          end # end of f = 1:factors

          for j = 1:pwEng.numberRates
            nextV = pwEng.V[i][stepToUse, j] * pwEng.LIBORRatios[stepToUse, j]
            pwEng.V[i][nextStepIndex, j] = nextV

            summandTerm = 0.0
            for f = 1:factors
              summandTerm += thisPseudoRoot[j, f] * pwEng.partials[f, j]
            end

            summandTerm *= taus[j] * pwEng.StepsDiscountsSquared[stepToUse, j]
            pwEng.V[i][nextStepIndex, j] += summandTerm
          end # end of for j = 1:numberRates
        end # end of i = 1:numberProducts
      end # end of if nextStepIndex != stepToUse
    end # end of if flowsFound
  end # end of for currentStep = numberSteps:-1:1

  # all V matricies computed, we now compute the elementary vegas for this path
  @inbounds for i = 1:pwEng.numberProducts, j = 1:pwEng.numberSteps
    nextIndex = j+1

    # we know V, we need to pair against the sensitivity of the rate to the elementary vega
    # note the simplification here arising from the fact that the elementary vega affects the evolution on precisely one step
    for k = 1:pwEng.numberRates, f = 1:pwEng.factors
      sensitivity = 0.0
      @simd for r = 1:pwEng.numberRates
        sensitivity += pwEng.V[i][nextIndex, r] * pwEng.jacobiansThisPaths[j][r][k, f]
      end

      pwEng.elementaryVegasThisPath[i][j][k, f] = sensitivity
    end
  end

  # write answer into values
  entriesPerProduct = 1 + pwEng.numberRates + pwEng.numberElementaryVegas

  @inbounds @simd for i = 1:pwEng.numberProducts
    values[((i - 1) * entriesPerProduct) + 1] = pwEng.numerairesHeld[i] * pwEng.initialNumeraireValue

    for j = 1:pwEng.numberRates
      values[(i - 1) * entriesPerProduct + 1 + j] = pwEng.V[i][1, j] * pwEng.initialNumeraireValue
    end

    for k = 1:pwEng.numberSteps, l = 1:pwEng.numberRates, m = 1:pwEng.factors
      values[(i-1) * entriesPerProduct + pwEng.numberRates + 1 + m + (l - 1) * factors + (k - 1) * pwEng.numberRates * pwEng.factors] =
            pwEng.elementaryVegasThisPath[i][k][l, m] * pwEng.initialNumeraireValue
    end
  end

  return 1.0 # we have put the weight in already, this results in lower variance since weight changes along the path
end

function multiple_path_values_elementary!(pwEng::PathwiseVegasOuterAccountingEngine, means::Vector{Float64}, errors::Vector{Float64}, numberOfPaths::Int)
  numberOfElementaryVegas = pwEng.numberRates * pwEng.numberSteps * pwEng.factors
  values = zeros(number_of_products(pwEng.product) * (1 + pwEng.numberRates + numberOfElementaryVegas))
  resize!(means, length(values))
  resize!(errors, length(values))

  sums = zeros(length(values))
  sumsqs = zeros(length(values))

  @inbounds @simd for i = 1:numberOfPaths
    single_path_values!(pwEng, values)
    for j in eachindex(values)
      sums[j] += values[j]
      sumsqs[j] += values[j] * values[j]
    end
  end

  @inbounds @simd for j in eachindex(values)
    means[j] = sums[j] / numberOfPaths
    meanSq = sumsqs[j] / numberOfPaths
    variance = meanSq - means[j] * means[j]
    errors[j] = sqrt(variance / numberOfPaths)
  end

  return means, errors
end

function multiple_path_values!(pwEng::PathwiseVegasOuterAccountingEngine, means::Vector{Float64}, errors::Vector{Float64}, numberOfPaths::Int)
  allMeans = Vector{Float64}()
  allErrors = Vector{Float64}()

  multiple_path_values_elementary!(pwEng, allMeans, allErrors, numberOfPaths)

  outDataPerProduct = 1 + pwEng.numberRates + pwEng.numberBumps
  inDataPerProduct = 1 + pwEng.numberRates + pwEng.numberElementaryVegas

  resize!(means, (1 + pwEng.numberRates + pwEng.numberBumps) * pwEng.numberProducts)
  resize!(errors, (1 + pwEng.numberRates + pwEng.numberBumps) * pwEng.numberProducts)

  @inbounds for p = 1:pwEng.numberProducts, i = 1:pwEng.numberRates + 1
    means[i + (p - 1) * outDataPerProduct] = allMeans[i + (p - 1) * outDataPerProduct]
    errors[i + (p - 1) * outDataPerProduct] = allErrors[i + (p - 1) * outDataPerProduct]

    @simd for bump = 1:pwEng.numberBumps
      thisVega = 0.0
      for t = 1:pwEng.numberSteps, r = 1:pwEng.numberRates, f = 1:pwEng.factors
        thisVega += pwEng.vegaBumps[t][bump][r, f] * allMeans[(p - 1) * inDataPerProduct + 1 + pwEng.numberRates + (t - 1) * pwEng.numberRates * pwEng.factors + (r - 1) * pwEng.factors + f]
      end

      means[(p - 1) * outDataPerProduct + 1 + pwEng.numberRates + bump] = thisVega
    end
  end

  return means, errors
end


text:
mux(f) = f
mux(m, f) = (xs...) -> m(f, xs...)
mux(ms...) = foldr(mux, ms)

type Context{T}
  interp::T
  cache::ObjectIdDict
  stack::Vector{Any}
  data::Dict{Symbol,Any}
end

Context(interp; kws...) = Context(interp, ObjectIdDict(), [], Dict{Symbol,Any}(kws))

Base.getindex(ctx::Context, k::Symbol) = ctx.data[k]
Base.setindex!(ctx::Context, v, k::Symbol) = ctx.data[k] = v

function stack(c::Context)
  stk = []
  isempty(c.stack) && return stk
  frame = nothing
  for i = 1:length(c.stack)
    isa(c.stack[i], Frame) || continue
    i > 1 && isa(c.stack[i-1], Line) && unshift!(stk, (frame, c.stack[i-1]))
    frame = c.stack[i].f
  end
  isa(c.stack[end], Line) && unshift!(stk, (frame, c.stack[end]))
  return stk
end

function interpv(ctx::Context, graph::IVertex)
  haskey(ctx.cache, graph) && return ctx.cache[graph]
  ctx.cache[graph] = ctx.interp(ctx, value(graph), inputs(graph)...)
end

interpv(ctx::Context, xs::Tuple) = map(x -> interpv(ctx, x), xs)

function interpret(ctx::Context, graph::IVertex, args::IVertex...)
  graph = spliceinputs(graph, args...)
  interpv(ctx, graph)
end

interpret(ctx::Context, graph::IVertex, args...) =
  interpret(ctx, graph, map(constant, args)...)

# The `ifoo` convention denotes a piece of interpreter middleware

iconst(f, ctx::Context, x::Constant) = x.value

function iline(f, ctx::Context, l::Union{Line,Frame}, v)
  push!(ctx.stack, l)
  val = interpv(ctx, v)
  pop!(ctx.stack)
  return val
end

ilinev(f, ctx::Context, l::Union{Line,Frame}, v) = vertex(l, iline(f, ctx, l, v))

ilambda(f, ctx::Context, ::Flosure, body, vars...) =
  (xs...) -> interpret(ctx, flopen(body), vars..., xs...)

iargs(cb, ctx::Context, f, xs...) = cb(f, interpv(ctx, xs)...)

function ituple(f, s::Split, xs)
  isa(xs, Vertex) && value(xs) == tuple ? inputs(xs)[s.n] :
  isa(xs, Tuple) ? xs[s.n] :
    f(s, constant(xs))
end

for m in :[iconst, iline, ilinev, ilambda, ituple].args
  @eval $m(f, args...) = f(args...)
end

interpeval = mux(iline, ilambda, iconst, iargs, ituple, (f, xs...) -> f(xs...))

interpret(graph::IVertex, args...) =
  interpret(Context(interpeval), graph, args...)

# Error Handling

import Juno: errmsg, errtrace

framename(f::Function) = typeof(f).name.mt.name
framename(f::Void) = Symbol("<none>")
framename(x) = symbol(string(typeof(x)))

totrace(stack) = [StackFrame(framename(f), Symbol(line.file), line.line)
                  for (f, line) in stack]

type Exception{T}
  err::T
  trace::StackTrace
end

errmsg(e::Exception) = errmsg(e.err)
errtrace(e::Exception, bt) = errtrace(e.err, [e.trace..., bt...])


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/05/2014

module Encounter

export
    AbstractEncounterModel,
    AddObserver,

    update,

    generateEncounter,
    getInitialState,
    getNextCommand,
    getTrajectory,

    CorrAEM,
    CorrAEMInitialState,
    CorrAEMCommand,
    generateEncountersToFile,
    setInitialDistributions,
    validate,
    getInitialSample,

    LLAEM,
    LLAEMInitialState,
    LLAEMCommand

import CommonInterfaces: initialize, update
using AbstractEncounterModelImpl

using CorrAEMImpl
using LLAEMImpl

end


text:
using Petri
import Petri: evaluate
using LabelledArrays
using ModelingToolkit
using MacroTools
import MacroTools: postwalk, striplines
using Test
import GeneralizedGenerated: mk_function

@testset "Stochastic Simulation" begin
    @variables S,E,I,R, β,γ,μ
    @testset "SIR" begin
        u0 = @LArray [100, 1, 0] (:S, :I, :R)
        params = @LArray [0.15, 0.55/101] (:β, :μ)
        sir = Petri.Model([S,I,R], [(S+I, 2*I), (I,R)])
        sir′ = funckit(sir)
        m′ = evaluate(sir′)
        p = Petri.ParamProblem(m′, u0, params, 250)
        soln = Petri.solve(p)
        @test soln.S <= 10
        @test sum(soln) == 101
    end
    @testset "SEIRS" begin
        seirs = Petri.Model([S,E,I,R],[(S+I, E+I), (E,I), (I,R), (R,S)])

        # Δ1 = @show Δ(seirs)
        # Λ1 = @show Λ(seirs)
        m′ = funckit(seirs)

        @show m′
        # @test all(Δ1 .== m′.Δ)
        # @test all(Λ1 .== m′.Λ)

        m′′ = Petri.evaluate(m′)
        u0 = @LArray [100, 1, 0, 0] (:S, :E, :I, :R)
        params = @LArray [0.55/101, 0.15, 0.15, 0.15] (:β, :γ, :μ, :η)
        p = Petri.ParamProblem(m′′, u0, params, 250)
        seirs_soln = Petri.solve(p)
        @test seirs_soln.S >= 10
        @test sum(seirs_soln) == 101
    end
end

function stepfunc(m::Model)
    transitionnames = map(m.Δ) do f
        f.args[1].args[1]
    end
    n = length(m.Λ)
    ratecalc = map(enumerate(m.Λ)) do (i,f)
        body = f.args[2] |> striplines
        if length(body.args) == 1
            body = body.args[1]
        end

        quote λ[$i] = $body end
    end
    actionchain = map(enumerate(m.Δ)) do (i, f)
        body = f.args[2] |> striplines
        if length(body.args) == 1
            body = body.args[1]
        end
        # push!(body.args, :(return state))
        quote if nexti == $i
            $body
        end
        end
    end
    fdef = quote
        (p, state) -> begin
            params = p.param
            λ = zeros(Float64, $n)
            $(ratecalc...)
            nexti = Petri.sample(λ)
            if nexti == nothing
                return state
            end
            $(actionchain...)
            return state
        end
    end
    fdef
end

@testset "Unrolled Stochastic Solve" begin
    @variables S,E,I,R, β,γ,μ
    seirs = Petri.Model([S,E,I,R],[(S+I, E+I), (E,I), (I,R), (R,S)])
    m′ = funckit(seirs)
    u0 = @LArray [100, 1, 0, 0] (:S, :E, :I, :R)
    params = @LArray [0.55/101, 0.15, 0.15, 0.15] (:β, :γ, :μ, :η)
    p = Petri.ParamProblem(m′, u0, params, 250)
    fexp = stepfunc(p.m)
    arglist = [:p,:state]
    body = fexp.args[2]
    # f1 = mk_function(@__MODULE__, arglist, [], body)
    # f = (p, state) -> f1(p, state)
    f = eval(body)
    @show state′ = f(p, p.initial)
    sol = Petri.solve(p, f)
    @show sol
end


text:
module Arrays
#####
##### PersistentArray
#####

export PersistentArray

using ..Lib
using ..Transaction
using ..Persist

import Base: size, sizeof, getindex, setindex!, IndexStyle, pointer, unsafe_convert, similar, elsize

struct ArrayHandle{T,N}
    size::NTuple{N, Int}
    base::PersistentOID{T}
end

mutable struct PersistentArray{T,N} <: DenseArray{T,N}
    # A handle to the underlying base for the array
    handle::Persistent{ArrayHandle{T,N}}

    # Cache the size and pointer to the base, these elements are volatile.
    size::NTuple{N, Int}
    base::Ptr{T}
end

function PersistentArray{T}(pool, size::NTuple{N, <:Integer}) where {T,N}
    handle = transaction(pool) do
        # Allocate space for the array in persistent memory
        allocsize = prod(size) * sizeof(T)
        baseoid = Lib.tx_alloc(allocsize, T)

        persist(pool, ArrayHandle{T,N}(size, baseoid))
    end
    return PersistentArray(handle)
end

function PersistentArray(handle::Persistent{ArrayHandle{T,N}}) where {T,N}
    _handle = retrieve(handle)

    # We can instantiate the array directly from the handle
    size = _handle.size
    base = Lib.direct(_handle.base)
    P = PersistentArray{T,N}(handle, size, base)
    finalizer(_free, P)
    return P
end

function _free(P::PersistentArray)
    # Free the array than free the handle
    Lib.free(Lib.oid(P.base))
    Persist.destroy(P.handle)
end

#####
##### customizations
#####

function Transaction.transaction(f, P::PersistentArray{T}) where {T}
    pool = getpool(P.handle)
    return transaction(pool) do
        Lib.add_range_direct(P.base, prod(size(P)) * sizeof(T))
        f()
    end
end

#####
##### Array Interface Methods
#####

# Conversion
pointer(P::PersistentArray) = P.base
unsafe_convert(::Type{Ptr{T}}, P::PersistentArray{T}) where {T} = pointer(P)

# Should maintain invariant that the unsafe_load is always valid.
size(P::PersistentArray) = P.size
sizeof(P::PersistentArray{T}) where {T} = prod(P.size) * sizeof(T)
elsize(P::PersistentArray{T}) where {T} = sizeof(T)

getindex(P::PersistentArray, i::Integer) = unsafe_load(P.base, i)
setindex!(P::PersistentArray, v, i::Integer) = _unsafe_setindex!(P, v, i)
IndexStyle(::Type{<:PersistentArray}) = IndexLinear()

function _safe_setindex!(P::PersistentArray{T}, v, i::Integer) where {T}
    # Get the object pool backing this array
    pool = Lib.pool(P.base)

    # Store the result in a transaction.
    transaction(pool) do
        Lib.add_range_direct(P.base + i - 1, sizeof(T))
        unsafe_store!(P.base, v, i)
    end
    return nothing
end

_unsafe_setindex!(P::PersistentArray{T}, v, i::Integer) where {T} = unsafe_store!(P.base, v, i)

# Similar
function similar(array::PersistentArray{T,N}, ::Type{U}, dims::NTuple{D,Int64}) where {T,N,U,D}
    # Store this in the same pool as the original array
    pool = Persist.getpool(array.handle)
    return PersistentArray{U}(pool, dims)
end

end


text:
@doc raw"""
    ChiSquared()

The Chi-Squared base function is given by:

```math
f(\mathbf{x}, \mathbf{y}) = \sum_i \frac{(x_i - y_i)^2}{x_i + y_i}
```
"""
struct ChiSquared <: PreMetric end

@inline function base_aggregate(::ChiSquared, s::T, x::T, y::T) where {T}
    x == y == zero(T) ? s : s + (x-y)^2/(x+y)
end