text:
# plot_polaron.jl

# Interactive Figures
# plotly()
# Plots.PlotlyBackend()

# Static Figures
# pyplot()
# default(fmt=:svg)

"""
plot_polaron(polaron::Polaron; N::Int)

    Creates the following plots of a given Polaron mutable struct:

        • μ_Ω: Mobility μ(Ω) (cm^2/Vs) versus Frequency Ω/ω (multiple of phonon frequency). Includes plots at N different evenly-spaced Temperatures T (K).
        • μ_T: Mobility μ(T) (cm^2/Vs) versus Temperature T (K). Includes plots at N different evenly-spaced Frequencies Ω/ω (multiple of phonon frequency).
        • Γ_Ω: Optical Absorption Γ(Ω) (cm^-1) versus Frequency Ω/ω (multiple of phonon frequency). Includes plots at N different evenly-spaced Temperatures T (K).
        • Γ_T: Optical Absorption Γ(Ω) (cm^-1) versus Temperature T (K). Includes plots at N different evenly-spaced Frequencies Ω/ω (multiple of phonon frequency).
        • κ_T: Spring Constant κ(T)/m_e (kg/s^2) (multiple of electron masses) versus Temperature T (K).
        • M_T: Fictitious mass M(T)/m_e (kg) (multiple of electron masses) versus Temperature T (K).
        • vw_T: Variational Parameters v and w (s^-1) versus Temperature T (K).
        • F_T: Free Energy (meV) versus Temperature T (K).
        ⦿ all_plots: A final plot that includes all the above as different subplots.

    polaron is a Polaron Mutable Stuct type created from make_polaron(). N is an Integer that specifies how many different temperatures or frequencies to plots on the Mobility and Optical Absorption plots.

    returns plots: all_plots, μ_Ω, μ_T, Γ_Ω, Γ_T, κ_T, M_T, vw_T, F_T
"""

function plot_polaron(polaron; N = 6)

    # Extract polaron data.
    α = polaron.α # Alpha parameter
    v = polaron.v # v variational parameter
    w = polaron.w # w variational parameter
    F = polaron.F # Free energy
    κ = polaron.κ # Spring constant
    M = polaron.M # Fictitious mass
    Ω = polaron.Ω # Electric field frequencies
    β = polaron.β # Reduced thermodynamic betas
    T = polaron.T # Temperatures
    μ = polaron.μ # Mobilities
    Γ = polaron.Γ # Optical absorptions

    # Plot Mobility versus Frequency for N different Temperatures.
    μ_Ω = Plots.plot(Ω, μ[1:Int(floor(length(T)/(N - 1))):end], label = hcat(["T = $i" for i in T][1:Int(floor(length(T)/(N - 1))):end]...), title = "Mobility \$\\mu(\\Omega)\$", ylim = (0, maximum(μ)[end] * 1.1), xlabel = "\$\\Omega\$ / \$\\omega\$", ylabel = "\$\\mu(\\Omega)\$ \$(cm^2/Vs)\$", legend = true, linewidth = 1.5, xtickfontsize = 10, ytickfontsize = 10, xguidefontsize = 10, yguidefontsize = 10, legendfontsize = 8, thickness_scaling = 1.2, size = (600, 600))

    # Plot Mobility versus Temperature for N different Frequencies.
    μ_T = Plots.plot(T, vcat.(μ...)[1:Int(floor(length(Ω)/(N - 1))):end], label = hcat(["Ω = $i" for i in Ω][1:Int(floor(length(Ω)/(N - 1))):end]...), title = "Mobility \$\\mu(T)\$", ylim = (0, maximum(μ)[end] * 1.1),  xlabel = "\$T\$ \$(K)\$", ylabel = "\$\\mu(T)\$ \$(cm^2/Vs)\$", legend = true, linewidth = 1.5, xtickfontsize = 10, ytickfontsize = 10, xguidefontsize = 10, yguidefontsize = 10, legendfontsize = 8, thickness_scaling = 1.2, size = (600, 600))

    # Plot Asbsorption versus Frequency for N different Temperatures.
    Γ_Ω = Plots.plot(Ω, Γ[1:Int(floor(length(T)/(N - 1))):end], label = hcat(["T = $i" for i in T][1:Int(floor(length(T)/(N - 1))):end]...), title = "Optical Absorption \$\\Gamma(\\Omega)\$", ylim = (0, sort(maximum(vcat.(Γ...)))[end - 1] * 1.1),  xlabel = "\$\\Omega\$ / \$\\omega\$", ylabel = "\$\\Gamma(\\Omega)\$ \$(cm^{-1})\$", legend = true, linewidth = 1.5, xtickfontsize = 10, ytickfontsize = 10, xguidefontsize = 10, yguidefontsize = 10, legendfontsize = 8, thickness_scaling = 1.2, size = (600, 600))

    # Plot Absorption versus Temperature for N different Frequencies.
    Γ_T = Plots.plot(T, vcat.(Γ...)[1:Int(floor(length(Ω)/(N - 1))):end], label = hcat(["Ω = $i" for i in Ω][1:Int(floor(length(Ω)/(N - 1))):end]...), title = "Optical Absorption \$\\Gamma(T)\$", ylim = (0, sort(maximum(vcat.(Γ...)))[end - 1] * 1.1),  xlabel = "\$T\$ \$(K)\$", ylabel = "\$\\Gamma(T)\$ \$(cm^{-1})\$", legend = true, linewidth = 1.5, xtickfontsize = 10, ytickfontsize = 10, xguidefontsize = 10, yguidefontsize = 10, legendfontsize = 8, thickness_scaling = 1.2, size = (600, 600))

    # Plot Spring Constant versus Temperature.
    κ_T = Plots.plot(T, κ, title = "Spring Constant \$\\kappa(T)\$", xlabel = "\$T\$ \$(K)\$", ylabel = "\$\\kappa(T)\$ / \$m_e\$ \$(kg/s^2)\$", legend = false, linewidth = 1.5, xtickfontsize = 10, ytickfontsize = 10, xguidefontsize = 10, yguidefontsize = 10, legendfontsize = 8, thickness_scaling = 1.2, size = (600, 600))

    # Plot Fictitious Mass versus Temperature.
    M_T = Plots.plot(T, M, title = "Fictitious Mass \$M(T)\$", xlabel = "\$T\$ \$(K)\$", ylabel = "\$M(T)\$ / \$m_e\$ \$(kg)\$", legend = false, linewidth = 1.5, xtickfontsize = 10, ytickfontsize = 10, xguidefontsize = 10, yguidefontsize = 10, legendfontsize = 8, thickness_scaling = 1.2, size = (600, 600))

    # Plot Variational Parameters v & w versus Temperature.
    vw_T = Plots.plot(T, [v, w], label = ["v" "w"], title = "Variational Parameters \$v(T)\$ & \$w(T)\$", xlabel = "\$T\$ \$(K)\$", ylabel = "\$v\$ & \$w\$ \$(s^{-1})\$", legend = true, linewidth = 1.5, xtickfontsize = 10, ytickfontsize = 10, xguidefontsize = 10, yguidefontsize = 10, legendfontsize = 8, thickness_scaling = 1.2, size = (600, 600))

    # Plot Free Energy versus Temperature.
    F_T = Plots.plot(T, F, title = "Free Energy \$F(T)\$", xlabel = "\$T\$ \$(K)\$", ylabel = "\$F(T)\$ \$(meV)\$", legend = false, linewidth = 1.5, xtickfontsize = 10, ytickfontsize = 10, xguidefontsize = 10, yguidefontsize = 10, legendfontsize = 8, thickness_scaling = 1.2, size = (600, 600))

    # Combine all plots above as subplots for combined view.
    all_plots = Plots.plot(μ_Ω, μ_T, Γ_Ω, Γ_T, κ_T, M_T, vw_T, F_T, layout = (2, 4), size = (1800, 900), linewidth = 1.5, xtickfontsize = 10, ytickfontsize = 8, xguidefontsize = 10, yguidefontsize = 10, legendfontsize = 10, thickness_scaling = 1.2)

    # Return all plots.
    all_plots, μ_Ω, μ_T, Γ_Ω, Γ_T, κ_T, M_T, vw_T, F_T
end

"""
save_polaron_plots(plots::Array{Plot}, path::String, ext::String)

    Saves all the plots created from plot_polaron() at the specified path with the specified extension.
"""

function save_polaron_plots(plots, path, ext = "svg")
    for plot in plots
        savefig(plot, path * "$(plot)." * ext)
    end
end

function plot_polaron_interactive(ϵ_optic, ϵ_static, phonon_freq, m_eff; t = 100:15:400, f = 0:1:20)

    fig = GLMakie.Figure(resolution = (1920, 1080), backgroundcolor = RGBf0(0.9, 0.9, 0.9))

    energy_axis = fig[1, 3] = Axis(fig, xlabel = "T [K]", ylabel = "F(T) [meV]", title = "Free Energy")
    vary_axis = fig[2, 3] = Axis(fig, xlabel = "T [K]", ylabel = "v/w(T) [Hz]", title = "Variational Parameters")
    mass_axis = fig[1, 4] = Axis(fig, xlabel = "T [K]", ylabel = "M(T) [m_e kg]", title = "Fictitious Mass")
    spring_axis = fig[2, 4] = Axis(fig, xlabel = "T [K]", ylabel = "κ(T) [me kg/s^2]", title = "Fictitious Spring Constant")
    mobility_axis_temp = fig[1, 2] = Axis(fig, xlabel = "T [K]", ylabel = "μ(T) [cm^2/Vs]", title = "Mobility vs Temperature")
    absorption_axis_temp = fig[2, 2] = Axis(fig, xlabel = "T [K]", ylabel = "Γ(T) [cm^-1]", title = "Optical Absorption vs Temperature")
    mobility_axis_freq = fig[1, 5] = Axis(fig, xlabel = "Ω [Hz]", ylabel = "μ(Ω) [cm^2/Vs]", title = "Mobility vs Frequency")
    absorption_axis_freq = fig[2, 5] = Axis(fig, xlabel = "Ω [Hz]", ylabel = "Γ(Ω) [cm^-1]", title = "Optical Absorption vs Frequency")

    lsgrid = labelslidergrid!(
        fig,
        ["ϵ_0", "ϵ_∞", "m", "ω"],
        [1:0.01:100, 1:0.01:100, 0.1:0.01:1, 0.01:0.01:100],
        formats = [x -> "$(round(x, digits = 3)) $s" for s in ["F/m", "F/m", "m_e", "THz"]],
        width = 1000
    )

    set_close_to!(lsgrid.sliders[1], ϵ_static)
    set_close_to!(lsgrid.sliders[2], ϵ_optic)
    set_close_to!(lsgrid.sliders[3], m_eff)
    set_close_to!(lsgrid.sliders[4], phonon_freq)

    fig[3, 1:5] = lsgrid.layout

    Ω_slider = Slider(fig[1:2, 1], range = f, startvalue = 0.0, horizontal = false, tellwidth = true, height = nothing, width = Auto())
    freq_sliders = [s.value for s in vcat(lsgrid.sliders , [Ω_slider])]
    Ω_slider_label = Label(fig[1:2, 0], lift(s -> "Ω: $(s[]) Hz", Ω_slider.value), rotation = pi/2)

    polaron_temp = lift(freq_sliders...) do slvalues...
        make_polaron(slvalues[2], slvalues[1], slvalues[4] * 1e12, slvalues[3]; temp = t, efield_freq = slvalues[5], verbose = false)
    end

    T = lift(p -> p.T, polaron_temp)
    v = lift(p -> Vector{Float64}(p.v), polaron_temp)
    w = lift(p -> Vector{Float64}(p.w), polaron_temp)
    κ = lift(p -> Vector{Float64}(p.κ), polaron_temp)
    M = lift(p -> Vector{Float64}(p.M), polaron_temp)
    F = lift(p -> Vector{Float64}(p.F), polaron_temp)
    μ_T = lift(p -> vcat.(p.μ...)[1], polaron_temp)
    Γ_T = lift(p -> vcat.(p.Γ...)[1], polaron_temp)

    T_slider = Slider(fig[1:2, 7], range = t, startvalue = 300, horizontal = false, tellwidth = true, height = nothing, width = Auto())
    T_slider_label = Label(fig[1:2, 8], lift(s -> "T: $(s[]) K", T_slider.value), rotation = -pi/2)

    temp_sliders = [s.value for s in vcat(lsgrid.sliders , [T_slider])]

    polaron_freq = lift(temp_sliders...) do slvalues...
        make_polaron(slvalues[2], slvalues[1], slvalues[4] * 1e12, slvalues[3]; temp = slvalues[5], efield_freq = f, verbose = false)
    end

    Ω = lift(p -> p.Ω, polaron_freq)
    μ_Ω = lift(p -> Vector{Float64}(p.μ[1]), polaron_freq)
    Γ_Ω = lift(p -> Vector{Float64}(p.Γ[1]), polaron_freq)

    scatterlines!(mobility_axis_temp, T, μ_T)
    scatterlines!(absorption_axis_temp, T, Γ_T)
    scatterlines!(energy_axis, T, F)
    v_plot = scatterlines!(vary_axis, T, v, label = "v")
    w_plot = scatterlines!(vary_axis, T, w, label = "w", marker = :diamond)
    scatterlines!(spring_axis, T, κ)
    scatterlines!(mass_axis, T, M)

    axislegend(vary_axis, [v_plot, w_plot], ["v", "w"], position = :rb)

    on(polaron_temp) do x
        autolimits!(mobility_axis_temp)
        autolimits!(absorption_axis_temp)
        autolimits!(energy_axis)
        autolimits!(vary_axis)
        autolimits!(spring_axis)
        autolimits!(mass_axis)
    end

    scatterlines!(mobility_axis_freq, Ω, μ_Ω)
    scatterlines!(absorption_axis_freq, Ω, Γ_Ω)

    on(polaron_freq) do x
        autolimits!(mobility_axis_freq)
        autolimits!(absorption_axis_freq)
    end

    fig
end


text:
#Turn off precompilations because of GG bug https://github.com/cscherrer/Soss.jl/issues/267
const fwhmfac = 2*sqrt(2*log(2))

using Soss
using Comrade
using MeasureTheory
import Distributions
const Dists = Distributions

gamps = @model stations, spriors begin
    #Gain amps
    σ ~ For(eachindex(spriors)) do i
        Dists.LogNormal(zero(spriors[i]), spriors[i])
    end
    g = NamedTuple{stations}(σ)
    return g
end

gphases = @model stations, spriors begin
    σ ~ For(2:(length(spriors))) do i
            Dists.truncated(Dists.Normal(zero(spriors[i]), spriors[i]), -π, π)
        end
    g = NamedTuple{stations}(((zero(eltype(σ)),)..., σ...))
    return g
end

function amp_gains(img, g1,g2, u, v)
    return g1*g2*Comrade.amplitude(img, u, v)
end

function _vamps(img, g, s1, s2, uamp, vamp)
    vm = similar(uamp, eltype(g))
    for i in eachindex(uamp)
        vm[i] = amp_gains(img,
                            g[s1[i]],
                            g[s2[i]],
                            uamp[i],
                            vamp[i]
                            )
    end
    return vm
end

va = @model image, gamps, uamp, vamp, s1, s2, erramp begin
    img ~ image
    g ~ gamps
    vm = _vamps(img, g, s1, s2, uamp, vamp)
    amp ~ Dists.MvNormal(vm, erramp)
end

vacp = @model image, gamps, uamp, vamp, s1, s2, erramp,
               u1cp, v1cp, u2cp, v2cp, u3cp, v3cp, errcp begin

    img ~ image
    g ~ gamps
    vm = _vamps(img, g, s1, s2, uamp, vamp)
    amp ~ Dists.MvNormal(vm, erramp)

    cp = Comrade.closure_phase.(Ref(img),
                            u1cp,
                            v1cp,
                            u2cp,
                            v2cp,
                            u3cp,
                            v3cp)

    cphase ~ Comrade.CPVonMises{(:μ,:σ)}(mcp, errcp)
                            #For(eachindex(cp)) do i
    #    Comrade.CPVonMises(cp[i], errcp[i])
    #end

end

function vis_gains(img, ga1,ga2, gp1,gp2, u, v)
    Δθ = gp1 - gp2
    s,c = sincos(Δθ)
    g1 = ga1
    g2 = ga2
    vis = visibility(img, u, v)
    vr = g1*g2*(real(vis)*c + imag(vis)*s)
    vi = g1*g2*(-real(vis)*s + imag(vis)*c)
    return vr+vi*im
end

function _vism(img, ga, gp, s1, s2, u, v)
    vm = similar(u,Complex{eltype(u)})
    for i in eachindex(u)
        vm[i] = vis_gains(img,
                          ga[s1[i]],
                          ga[s2[i]],
                          gp[s1[i]],
                          gp[s2[i]],
                          u[i],
                          v[i]
                        )
    end
    return vm
end

vis = @model image, gamps, gphases, u, v, s1, s2, error begin
    img ~ image
    ga ~ gamps
    gp ~ gphases
    vis = _vism(img, ga, gp, s1, s2, u, v)

    visr ~ Dists.MvNormal(real.(vis), err)
    visi ~ Dists.MvNormal(imag.(vis), err)

end

viswnoise = @model image, gamps, gphases, u, v, s1, s2, error begin
    img ~ image
    ga ~ gamps
    gp ~ gphases
    f ~ Dists.Uniform(0.0, 0.25)
    vis = _vism(img, ga, gp, s1, s2, u, v)
    vamps = abs.(vis)
    errwnoise = hypot.(err, f.*vamps)
    visr ~ Dists.MvNormal(real.(vis), errwnoise)
    visi ~ Dists.MvNormal(imag.(vis), errwnoise)
end

lcacp = @model image,
               u1a,v1a,u2a,v2a,u3a,v3a,u4a,v4a,errcamp,
               u1p,v1p,u2p,v2p,u3p,v3p,errcp begin
    img ~ image

    mlca = Comrade.logclosure_amplitude.(Ref(img), u1a, v1a,
                                               u2a, v2a,
                                               u3a, v3a,
                                               u4a, v4a,
                                    )
    lcamp ~ Dists.MvNormal(mlca, errcamp)
    #For(eachindex(mlca,errcamp)) do i
    #    Dists.Normal(mlca[i], errcamp[i])
    #end

    mcp = Comrade.closure_phase.(Ref(img),
                             u1cp,
                             v1cp,
                             u2cp,
                             v2cp,
                             u3cp,
                             v3cp)
    #cphase ~ Dists.MvNormal(mcp, errcp)
    cphase ~ Comrade.CPVonMises{(:μ,:σ)}(mcp, errcp)
end

#=
mringwb = @model N, M, fov begin
    diam ~ Dists.Uniform(25.0, 85.0)
    fwhm ~ Dists.Uniform(1.0, 40.0)
    rad = diam/2
    σ = fwhm/fwhmfac

    ma ~ Dists.Uniform(0.0, 0.5) |> iid(N)
    mp ~ Dists.Uniform(-1π, 1π) |> iid(N)
    α = ma.*cos.(mp)
    β = ma.*sin.(mp)

    #Fraction of floor flux
    floor ~ Dists.Uniform(0.0, 1.0)
    f ~ Dists.Uniform(0.8, 1.2)
    mring = renormed(Comrade.MRing{N}(rad, α, β), f-f*floor)
    rimg = smoothed(mring,σ)

    coeff ~ Dists.Dirichlet(fill(M*M, 1.0))
    ϵ ~ Dists.LogNormal(0.0, 0.1)
    bimg = renormed(stretched(Comrade.RImage(coeff, Comrade.SqExpKernel(ϵ)),fov,fov), f*floor)
    img = rimg+bimg
    return img
end
=#
mring = @model N, fmin, fmax begin
    diam ~ Dists.Uniform(25.0, 85.0)
    fwhm ~ Dists.Uniform(1.0, 40.0)
    rad = diam/2
    σ = fwhm/fwhmfac

    ma ~ Dists.Uniform(0.0, 0.5) |> iid(N)
    mp ~ Dists.Uniform(-1π, 1π) |> iid(N)
    α = ma.*cos.(mp)
    β = ma.*sin.(mp)

    f ~ Dists.Uniform(fmin, fmax)
    mring = renormed(Comrade.MRing{N}(rad, α, β), f)
    img = smoothed(mring,σ)
    return img
end

mringwfloor = @model N, fmin, fmax begin
    diam ~ Dists.Uniform(25.0, 85.0)
    fwhm ~ Dists.Uniform(1.0, 40.0)
    rad = diam/2
    σ = fwhm/fwhmfac

    ma ~ Dists.Uniform(0.0, 0.5) |> iid(N)
    mp ~ Dists.Uniform(-1π, 1π) |> iid(N)
    α = ma.*cos.(mp)
    β = ma.*sin.(mp)

    #Fraction of floor flux
    floor ~ Dists.Uniform(0.0, 1.0)
    f ~ Dists.Uniform(fmin, fmax)

    mring = renormed(Comrade.MRing{N}(rad, α, β), f*(1-floor))
    disk = renormed(stretched(Comrade.Disk(), rad, rad), f*floor)
    img = smoothed(mring+disk,σ)
    return img
end

mringwffloor = @model N, fmin, fmax begin
    diam ~ Dists.Uniform(25.0, 85.0)
    fwhm ~ Dists.Uniform(1.0, 40.0)
    rad = diam/2
    σ = fwhm/fwhmfac

    ma ~ Dists.Uniform(0.0, 0.5) |> iid(N)
    mp ~ Dists.Uniform(-1π, 1π) |> iid(N)
    α = ma.*cos.(mp)
    β = ma.*sin.(mp)

    #Fraction of floor flux
    floor ~ Dists.Uniform(0.0, 1.0)
    f ~ Dists.Uniform(fmin, fmax)
    σf ~ Dists.Uniform(1.0, 40.0)

    mring = renormed(Comrade.MRing{N}(rad, α, β), f*(1-floor))
    disk = renormed(stretched(Comrade.Disk(), rad, rad), f*floor)
    img = smoothed(mring, σ) +smoothed(disk,σf)
    return img
end

mringwgfloor = @model N, fmin, fmax begin
    diam ~ Dists.Uniform(25.0, 85.0)
    fwhm ~ Dists.Uniform(1.0, 40.0)
    rad = diam/2
    σ = fwhm/fwhmfac

    ma ~ Dists.Uniform(0.0, 0.5) |> iid(N)
    mp ~ Dists.Uniform(-1π, 1π) |> iid(N)
    α = ma.*cos.(mp)
    β = ma.*sin.(mp)

    #Fraction of floor flux
    floor ~ Dists.Uniform(0.0, 1.0)
    f ~ Dists.Uniform(fmin, fmax)
    dg ~ Dists.Uniform(40.0, 250.0)
    rg = dg/fwhmfac
    mring = smoothed(renormed(Comrade.MRing{N}(rad, α, β), f*(1-floor)), σ)
    g = renormed(stretched(Comrade.Gaussian(), rg, rg), f*floor)
    img = mring + g
    return img
end

smring = @model N, fmin, fmax begin
    diam ~ Dists.Uniform(25.0, 85.0)
    fwhm ~ Dists.Uniform(1.0, 40.0)
    rad = diam/2
    σ = fwhm/fwhmfac

    ma ~ Dists.Uniform(0.0, 0.5) |> iid(N)
    mp ~ Dists.Uniform(-1π, 1π) |> iid(N)
    α = ma.*cos.(mp)
    β = ma.*sin.(mp)

    #Stretch
    τ ~ Dists.Uniform( 0.0, 0.5)
    ξτ ~ Dists.Uniform(-π/2, π/2)
    scx = 1/sqrt(1-τ)
    scy = sqrt(1-τ)

    f ~ Dists.Uniform(fmin, fmax)

    mring = renormed(Comrade.MRing{N}(rad, α, β), f)
    img = smoothed(rotated(stretched(mring,scx,scy),ξτ),σ)
    return img
end

smringwfloor = @model N, fmin, fmax begin
    diam ~ Dists.Uniform(25.0, 85.0)
    fwhm ~ Dists.Uniform(1.0, 40.0)
    rad = diam/2
    σ = fwhm/fwhmfac

    ma ~ Dists.Uniform(0.0, 0.5) |> iid(N)
    mp ~ Dists.Uniform(-1π, 1π) |> iid(N)
    α = ma.*cos.(mp)
    β = ma.*sin.(mp)

    #Stretch
    τ ~ Dists.Uniform(0.0, 0.5)
    ξτ ~ Dists.Uniform(-π/2, π/2)
    scx = 1/sqrt(1-τ)
    scy = sqrt(1-τ)

    #Fraction of floor flux
    floor ~ Dists.Uniform(0.0, 1.0)
    f ~ Dists.Uniform(fmin, fmax)

    mring = renormed(Comrade.MRing{N}(rad, α, β), f*(1-floor))
    disk = renormed(stretched(Comrade.Disk(), rad, rad), f*floor)
    img = smoothed(rotated(stretched(mring+disk,scx,scy),ξτ),σ)
    return img
end

ext: jl


text:
"""
    CenteredAxis(origin=0, indices)

A `CenteredAxis` takes `indices` and provides a user facing set of keys centered around zero.
The `CenteredAxis` is a subtype of `AbstractOffsetAxis` and its keys are treated as the predominant indexing style.
Note that the element type of a `CenteredAxis` cannot be unsigned because any instance with a length greater than 1 will begin at a negative value.

## Examples

A `CenteredAxis` sends all indexing arguments to the keys and only maps to the indices when `to_index` is called.
```jldoctest
julia> using AxisIndices

julia> axis = AxisIndices.CenteredAxis(1:10)
center(0)(SimpleAxis(1:10))

julia> axis[10]  # the indexing goes straight to keys and is centered around zero
ERROR: BoundsError: attempt to access center(0)(SimpleAxis(1:10)) at index [10]
[...]

julia> axis[-4]
-4

```
"""
struct CenteredAxis{I,Inds,F} <: AbstractOffsetAxis{I,Inds,F}
    origin::F
    parent::Inds

    function CenteredAxis{I,Inds,F}(origin::Integer, inds::AbstractAxis) where {I,Inds,F}
        if inds isa Inds && origin isa F
            return new{I,Inds,F}(origin, inds)
        else
            return CenteredAxis(origin, convert(Inds, inds))
        end
    end
    function CenteredAxis{I,Inds,F}(inds::AbstractRange) where {I,Inds,F<:StaticInt}
        return new{I,Inds,F}(F(), inds)
    end
    function CenteredAxis{I,Inds,F}(inds::AbstractRange) where {I,Inds,F}
        return new{I,Inds,F}(F(0), inds)
    end

    function CenteredAxis{I,Inds,F}(origin::Integer, inds::AbstractRange) where {I,Inds,F}
        return CenteredAxis{I,Inds}(origin, inds)
    end

    function CenteredAxis{I,Inds}(inds::AbstractRange) where {I,Inds}
        return CenteredAxis{I,Inds}(Zero(), inds)
    end
    function CenteredAxis{I,Inds}(origin::Integer, inds::AbstractArray) where {I,Inds}
        return CenteredAxis{I,Inds}(origin, compose_axis(inds))
    end
    function CenteredAxis{I,Inds}(origin::Integer, inds::AbstractAxis) where {I,Inds}
        if inds isa Inds
            return CenteredAxis{I,Inds,typeof(origin)}(origin, inds)
        else
            return CenteredAxis{I}(origin, convert(Inds, inds))
        end
    end

    function CenteredAxis{I}(origin::Integer, inds::AbstractArray) where {I}
        return CenteredAxis{I}(origin, compose_axis(inds))
    end
    function CenteredAxis{I}(origin::Integer, inds::AbstractOffsetAxis) where {I}
        return CenteredAxis{I}(origin, parent(inds))
    end
    function CenteredAxis{I}(origin::Integer, inds::AbstractAxis) where {I}
        if eltype(inds) <: I
            return CenteredAxis{I,typeof(inds)}(origin, inds)
        else
            return CenteredAxis{I}(origin, convert(AbstractUnitRange{I}, inds))
        end
    end
    CenteredAxis{I}(inds::AbstractRange) where {I} = CenteredAxis{I}(Zero(), inds)

    CenteredAxis(inds::AbstractRange) = CenteredAxis(Zero(), inds)
    function CenteredAxis(origin::Integer, inds::AbstractOffsetAxis)
        return CenteredAxis(origin, parent(inds))
    end
    function CenteredAxis(origin::Integer, inds::AbstractArray)
        return CenteredAxis(origin, compose_axis(inds))
    end
    function CenteredAxis(origin::Integer, inds::AbstractAxis)
        return CenteredAxis{eltype(inds)}(origin, inds)
    end
end

@inline Base.getproperty(axis::CenteredAxis, k::Symbol) = getproperty(parent(axis), k)

function ArrayInterface.unsafe_reconstruct(axis::CenteredAxis, inds; kwargs...)
    return CenteredAxis(origin(axis), unsafe_reconstruct(parent(axis), inds; kwargs...))
end

ArrayInterface.known_first(::Type{T}) where {T<:CenteredAxis{<:Any,<:Any,<:Any}} = nothing
function ArrayInterface.known_first(::Type{T}) where {Inds,F,T<:CenteredAxis{<:Any,Inds,StaticInt{F}}}
    if known_length(Inds) === nothing
        return nothing
    else
        return F - div(known_length(Inds), 2)
    end
end
Base.first(axis::CenteredAxis) = origin(axis) - div(length(parent(axis)), 2)

ArrayInterface.known_last(::Type{T}) where {T<:CenteredAxis{<:Any,<:Any,<:Any}} = nothing
function ArrayInterface.known_last(::Type{T}) where {Inds,F,T<:CenteredAxis{<:Any,Inds,StaticInt{F}}}
    if known_length(Inds) === nothing
        return nothing
    else
        return F - div(known_length(Inds), 2) + known_length(Inds)
    end
end
Base.last(axis::CenteredAxis) = last(parent(axis)) + _origin_to_offset(axis)

function _origin_to_offset(axis::CenteredAxis)
    p = parent(axis)
    return _origin_to_offset(first(p), length(p), origin(axis))
end
_origin_to_offset(start, len, origin) = (origin - div(len, 2one(start))) - start

origin(axis::CenteredAxis) = getfield(axis, :origin)
@inline function ArrayInterface.offsets(axis::CenteredAxis)
    inds = parent(axis)
    return (_origin_to_offset(static_first(inds), static_length(inds), origin(axis)),)
end

"""
    center(collection, origin)
    center(collection, origin)

Shortcut for creating [`CenteredAxis`](@ref).

## Examples
```jldoctest
julia> using AxisIndices

julia> AxisArray(ones(3), center(0))
3-element AxisArray(::Vector{Float64}
  • axes:
     1 = -1:1
)
      1
  -1  1.0
  0   1.0
  1   1.0

```
"""
struct Center <: AxisInitializer end
const center = Center()
axis_method(::Center, x, inds) = CenteredAxis(x, inds)
center(collection::AbstractArray) = center(collection, Zero())

"""
    CenteredArray(A::AbstractArray)

Provides centered axes for indexing `A`.

## Examples
```jldoctest
julia> using AxisIndices

julia> AxisIndices.CenteredArray(ones(3,3))
3×3 AxisArray(::Matrix{Float64}
  • axes:
     1 = -1:1
     2 = -1:1
)
      -1    0    1
  -1   1.0  1.0  1.0
  0    1.0  1.0  1.0
  1    1.0  1.0  1.0

```
"""
const CenteredArray{T,N,P,A<:Tuple{Vararg{<:CenteredAxis}}} = AxisArray{T,N,P,A}

CenteredArray(A::AbstractArray{T,N}) where {T,N} = CenteredArray{T,N}(A)

CenteredArray(A::AbstractArray{T,0}) where {T} = AxisArray(A)

CenteredArray{T}(A::AbstractArray) where {T} = CenteredArray{T,ndims(A)}(A)

CenteredArray{T,N}(A::AbstractArray) where {T,N} = CenteredArray{T,N,typeof(A)}(A)

CenteredArray{T,N,P}(A::CenteredArray{T,N,P}) where {T,N,P} = A
function CenteredArray{T,N,P}(A::CenteredArray) where {T,N,P}
    return CenteredArray{T,N,P}(parent(A))
end

function CenteredArray{T,N,P}(A::AbstractArray) where {T,N,P<:AbstractArray{T,N}}
    return CenteredArray{T,N,P}(convert(P, A))
end

function CenteredArray{T,N,P}(A::P) where {T,N,P<:AbstractArray{T,N}}
    axs = map(center, axes(A))
    return CenteredArray{T,N,P,typeof(axs)}(A, axs)
end

function CenteredArray{T}(init::ArrayInitializer, sz::Tuple=()) where {T}
    return CenteredArray{T,length(inds)}(init, sz)
end

function CenteredArray{T,N}(init::ArrayInitializer, sz::Tuple=()) where {T,N}
    return CenteredArray{T,N}(Array{T,N}(init, sz))
end

"""
    CenteredVector(v::AbstractVector)

Provides a centered axis for indexing `v`.

## Examples
```jldoctest
julia> using AxisIndices

julia> AxisIndices.CenteredVector(ones(3))
3-element AxisArray(::Vector{Float64}
  • axes:
     1 = -1:1
)
      1
  -1  1.0
  0   1.0
  1   1.0

```
"""
const CenteredVector{T,P<:AbstractVector{T},Ax1<:CenteredAxis} = CenteredArray{T,1,P,Tuple{Ax1}}

CenteredVector(A::AbstractVector) = CenteredArray{eltype(A)}(A)
CenteredVector{T}(A::AbstractVector) where {T} = CenteredArray{T,1}(A)

"""
    CenteredVector{T}(init::ArrayInitializer, sz::Integer)

Creates a vector with elements of type `T` of size `sz` and a centered axis.

## Examples
```jldoctest
julia> using AxisIndices

julia> AxisIndices.CenteredVector{Union{Missing, Int}}(missing, 3)
3-element AxisArray(::Vector{Union{Missing, Int64}}
  • axes:
     1 = -1:1
)
      1
  -1   missing
  0    missing
  1    missing

```
"""
function CenteredVector{T}(init::ArrayInitializer, arg) where {T}
    return CenteredVector{T}(Vector{T}(init, arg))
end

function print_axis(io::IO, axis::CenteredAxis)
    print(io, "center($(Int(origin(axis))))($(parent(axis)))")
end


text:
module CounterfactualExplanations

# Dependencies:
using Flux
import Flux.Losses
using LinearAlgebra

include("data/Data.jl")
using .Data

include("data_preprocessing/DataPreprocessing.jl")
using .DataPreprocessing
export CounterfactualData, select_factual, apply_domain_constraints

include("models/Models.jl")
using .Models
export AbstractFittedModel, AbstractDifferentiableModel,
    FluxModel, LogisticModel, BayesianLogisticModel,
    RTorchModel, PyTorchModel,
    probs, logits

include("generators/Generators.jl")
using .Generators
export AbstractGenerator, AbstractGradientBasedGenerator, GenericGenerator, GreedyGenerator, CounterfactualState,
    generate_perturbations, conditions_satisified, mutability_constraints

include("counterfactuals/Counterfactuals.jl")
using .Counterfactuals
export CounterfactualExplanation, initialize!, update!,
    total_steps, converged, terminated, path, target_probs

include("generate_counterfactual.jl")
export generate_counterfactual

end


text:
# Make sure we cover all edge cases in the plotting code
using Makie.Colors: Gray
@testset "plotting" begin
    @testset "NaN color" begin
        confusion = [
            NaN 0;
            1.0 0.5
        ]
        nan_confusion = plot_confusion_matrix(confusion, ["test1", "test2"], :Row)
        @testplot nan_confusion
        nan_custom_confusion = with_theme(ConfusionMatrix = (Heatmap=(nan_color=:red,), Text=(color=:red,))) do
            plot_confusion_matrix(confusion, ["test1", "test2"], :Row)
        end
        @testplot nan_custom_confusion
    end

    @testset "Kappa placement" begin
        classes = ["class $i" for i in 1:5]
        kappa_text_placement = with_theme(Kappas = (Text=(color=Gray(0.5),),)) do
            plot_kappas((1:5) ./ 5 .- 0.1, classes, (1:5) ./ 5, color = [Gray(0.4), Gray(0.2)])
        end
        @testplot kappa_text_placement

        kappa_text_placement_single = with_theme(Kappas = (Text=(color=:red,),)) do
            plot_kappas((1:5) ./ 5, classes, color = [Gray(0.4), Gray(0.2)])
        end
        @testplot kappa_text_placement_single
    end

    @testset "binary discriminiation calibration curves" begin
        rng = StableRNG(22)
        curves = [(LinRange(0, 1, 10), range(0, stop=i/2, length=10) .+ (randn(rng, 10) .* 0.1)) for i in -1:3]
        binary_discrimination_calibration_curves_plot = with_theme(
                BinaryDiscriminationCalibrationCurves = (
                    Ideal = (
                        linewidth = 3,
                        color = (:green, 0.5)
                    ),
                    CalibrationCurve = (
                        solid_color = :green,
                        markersize = 50, # should be overwritten by user kw
                        linewidth = 5,
                    ),
                    PerExpert = (
                        solid_color = :red,
                        linewidth=1
                    ),
                )
            ) do
            Lighthouse.plot_binary_discrimination_calibration_curves(
                curves[3],
                rand(rng, 5),
                curves[[1, 2, 4, 5]],
                nothing, nothing,
                "",
                markersize=10
            )
        end
        @testplot binary_discrimination_calibration_curves_plot
    end
end


text:
module AutoHOOT

using PyCall

export set_logger

# disable loggings from the python library
function set_logger(; disabled::Bool = true)
    logging = pyimport("logging")
    logger = logging.getLogger("optimizer")
    logger.disabled = disabled
end

include("backend.jl")
include("autodiff.jl")
include("graph_ops.jl")

end


text:
module MultiAssemblers

using Gridap
using Gridap.Helpers
using Gridap.MultiFESpaces: _compute_offsets

using SparseArrays
using SparseMatricesCSR

export MultiAssembler
export MultiSparseMatrixAssembler

import Gridap.Assemblers: assemble
import Gridap.Assemblers: assemble!
import Gridap.Assemblers: SparseMatrixAssembler
import Gridap.Assemblers: SparseMatrixAssembler
using  Gridap.Assemblers: sparse_from_coo

abstract type MultiAssembler{M<:AbstractMatrix,V<:AbstractVector} end

function assemble(
  ::MultiAssembler{M,V},
  ::Vararg{Tuple{<:MultiCellVector,<:CellNumber}})::V where {M,V}
  @abstractmethod
end

function assemble(
  ::MultiAssembler{M,V},
  ::Vararg{Tuple{<:MultiCellMatrix,<:CellNumber,<:CellNumber}})::M where {M,V}
  @abstractmethod
end

function assemble!(
  ::V,
  ::MultiAssembler{M,V},
  ::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where {M,V}
  @abstractmethod
end

function assemble!(
  ::M,
  ::MultiAssembler{M,V},
  ::Vararg{Tuple{<:MultiCellMatrix,<:CellNumber,<:CellNumber}}) where {M,V}
  @abstractmethod
end

function assemble(
  a::MultiAssembler,cv::MultiCellVector)
  l = length(cv)
  ide = IdentityCellNumber(Int,l)
  assemble(a,(cv,ide))
end

function assemble(
  a::MultiAssembler,cv::MultiCellMatrix)
  l = length(cv)
  ide = IdentityCellNumber(Int,l)
  assemble(a,(cv,ide,ide))
end

function assemble!(r,a::MultiAssembler,cv::MultiCellVector)
  l = length(cv)
  ide = IdentityCellNumber(Int,l)
  assemble!(r,a,(cv,ide))
end

function assemble!(r,a::MultiAssembler,cv::MultiCellMatrix)
  l = length(cv)
  ide = IdentityCellNumber(Int,l)
  assemble!(r,a,(cv,ide,ide))
end

"""
MultiAssembler that produces SparseMatrices from the SparseArrays package
"""
struct MultiSparseMatrixAssembler{E,M} <: MultiAssembler{M,Vector{E}}
  testfesps::MultiFESpace{E}
  trialfesps::MultiFESpace{E}

  function MultiSparseMatrixAssembler(
    ::Type{M},
    testfesp::MultiFESpace{E},
    trialfesp::MultiFESpace{E}) where {M<:AbstractSparseMatrix,E}
    new{E,M}(testfesp,trialfesp)
  end
end

function MultiSparseMatrixAssembler(
  testfesp::MultiFESpace{E},
  trialfesp::MultiFESpace{E}) where {E}
  MultiSparseMatrixAssembler(SparseMatrixCSC, testfesps, trialfesp)
end

function MultiSparseMatrixAssembler(
  ::Type{M},
  testfesps::Vector{<:FESpaceWithDirichletData},
  trialfesps::Vector{<:FESpaceWithDirichletData}) where {M<:AbstractSparseMatrix}
  V = MultiFESpace(testfesps)
  U = MultiFESpace(trialfesps)
  MultiSparseMatrixAssembler(M,V,U)
end

function MultiSparseMatrixAssembler(
  testfesps::Vector{<:FESpaceWithDirichletData},
  trialfesps::Vector{<:FESpaceWithDirichletData})
  MultiSparseMatrixAssembler(SparseMatrixCSC,V,U)
end

function SparseMatrixAssembler(
  ::Type{M},
  testfesps::Vector{<:FESpaceWithDirichletData},
  trialfesps::Vector{<:FESpaceWithDirichletData}) where {M<:AbstractSparseMatrix}
  MultiSparseMatrixAssembler(M,testfesps,trialfesps)
end

function SparseMatrixAssembler(
  testfesps::Vector{<:FESpaceWithDirichletData},
  trialfesps::Vector{<:FESpaceWithDirichletData})
  MultiSparseMatrixAssembler(SparseMatrixCSC,testfesps,trialfesps)
end

function MultiSparseMatrixAssembler(
  ::Type{M},
  testfesps::MultiFESpace{E,S},
  trialfesps::MultiFESpace{E,S}) where {M<:AbstractSparseMatrix,E,S}
  @notimplementedif S != ConsequtiveMultiFieldStyle
  MultiSparseMatrixAssembler{E,M}(M,testfesps,trialfesps)
end

function MultiSparseMatrixAssembler(
  testfesps::MultiFESpace{E,S},
  trialfesps::MultiFESpace{E,S}) where {E,S}
  @notimplementedif S != ConsequtiveMultiFieldStyle
  MultiSparseMatrixAssembler(SparseMatrixCSC,testfesps,trialfesps)
end

function assemble(
  this::MultiSparseMatrixAssembler{E},
  allvals::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where {E}

  n = num_free_dofs(this.testfesps)
  vec = zeros(E,n)
  assemble!(vec,this,allvals...)
  vec
end

function assemble!(
  vec::Vector{E},
  this::MultiSparseMatrixAssembler{E},
  allmcv::Vararg{Tuple{<:MultiCellVector,<:CellNumber}}) where {E}

  vec .= zero(E)
  V = this.testfesps
  offsets = _compute_offsets(V)
  _mf_rows = celldofids(V)
  for (mcv,cellids) in allmcv
    mf_vals = apply_constraints(V,mcv,cellids)
    i_to_fieldid = mf_vals.fieldids
    mf_rows = reindex(_mf_rows,cellids)
    _assemble_vec!(vec,mf_rows,mf_vals,i_to_fieldid,offsets)
  end
end

function assemble(
  this::MultiSparseMatrixAssembler{E,M},
  allmcm::Vararg{Tuple{<:MultiCellMatrix,<:CellNumber,<:CellNumber}}) where {E,M}

  V = this.testfesps
  U = this.trialfesps
  offsets_row = _compute_offsets(V)
  offsets_col = _compute_offsets(U)
  _mf_rows = celldofids(V)
  _mf_cols = celldofids(U)
  aux_row = Int[]; aux_col = Int[]; aux_val = E[]

  for (mcm,cellids_row,cellids_col) in allmcm
    mf_vals = apply_constraints_rows(V, mcm, cellids_row)
    mf_rows = reindex(_mf_rows, cellids_row)
    mf_vals = apply_constraints_cols(U,mf_vals,cellids_col)
    mf_cols = reindex(_mf_cols, cellids_col)
    i_to_fieldid = mf_vals.fieldids
    _assemble_mat!(M,
      aux_row,aux_col,aux_val,mf_rows,mf_cols,mf_vals,
      i_to_fieldid,offsets_row,offsets_col)
  end
  num_rows = num_free_dofs(this.testfesps)
  num_cols = num_free_dofs(this.trialfesps)
  finalize_coo!(M,aux_row,aux_col,aux_val,num_rows,num_cols)
  sparse_from_coo(M,aux_row, aux_col, aux_val)
end

function assemble!(
  mat::AbstractSparseMatrix{E},
  this::MultiSparseMatrixAssembler{E},
  allvals::Vararg{Tuple{<:MultiCellMatrix,<:CellNumber,<:CellNumber}}) where E
  # This routine can be optimized a lot taking into a count the sparsity graph of mat
  # For the moment we create an intermediate matrix and then transfer the nz values
  m = assemble(this,allvals...)
  nonzeros(mat) .= nonzeros(m)
end

function _assemble_vec!(vec,mf_rows,mf_vals,i_to_fieldid,offsets)
  for (mf_rows_c,mf_vals_c) in zip(mf_rows,mf_vals)
    for (i,vals_c) in enumerate(mf_vals_c)
      ifield, = i_to_fieldid[i]
      rows_c = mf_rows_c[ifield]
      offset = offsets[ifield]
      for (lid,gid) in enumerate(rows_c)
        if gid > 0
          vec[gid+offset] += vals_c[lid]
        end
      end
    end
  end
end

function _assemble_mat!(::Type{M},
  aux_row,aux_col,aux_val,mf_rows,mf_cols,mf_vals,
  i_to_fieldid,offsets_row,offsets_col) where {M<:AbstractSparseMatrix}

  for (mf_rows_c,mf_cols_c,mf_vals_c) in zip(mf_rows,mf_cols,mf_vals)
    for (i,vals_c) in enumerate(mf_vals_c)
      ifield,jfield = i_to_fieldid[i]
      rows_c = mf_rows_c[ifield]
      cols_c = mf_cols_c[jfield]
      offset_row = offsets_row[ifield]
      offset_col = offsets_col[jfield]
      for (lidcol,gidcol) in enumerate(cols_c)
        if gidcol > 0
          for (lidrow,gidrow) in enumerate(rows_c)
            if gidrow > 0
              push_coo!(M,aux_row, aux_col, aux_val,
                    gidrow+offset_row, gidcol+offset_col, vals_c[lidrow,lidcol])
            end
          end
        end
      end
    end
  end

end

end # module MultiAssemblers


text:
using Documenter, SimpleDifferentialOperators

# Compile the online documentation.
makedocs(sitename = "SimpleDifferentialOperators.jl",
	pages = [
        "index.md",
        "notebooks.md",
        "api.md"
    ],
	doctest = :fix)

# Deploy to git
deploydocs(
    repo = "github.com/QuantEcon/SimpleDifferentialOperators.jl.git",
)


text:
module ExternalBallistics
using ..Types, ..Utils, ..InternalBallistics

include("OneD.jl")
using .OneD
include("Pmm.jl")
using .Pmm
include("PmmZeroDrag.jl")
using .PmmZeroDrag
include("Mpmm.jl")
using .Mpmm
include("SixDof.jl")
using .SixDof
include("Flight.jl")
using .Flight

end


text:
#solve u_t + div(v * u) - div(D * grad(u)) - Q = 0 on an unstructured mesh using an upwind discretization of the velocity term (or downwind if use_upwind=false)
@NonlinearEquations.equations exclude=(areasoverlengths, volumes, dirichletnodes, coords) function transport(u, vxs, vys, vzs, Ds, neighbors, areasoverlengths, dirichletnodes, dirichletus, Qs, volumes, coords; use_upwind=true)
	@assert size(coords, 1) == 3#it is assumed that this is a 3d problem
	NonlinearEquations.setnumequations(length(u))
	for i = 1:length(u)
		if i in dirichletnodes
			NonlinearEquations.addterm(i, u[i] - dirichletus[i])
		else
			NonlinearEquations.addterm(i, Qs[i] / (volumes[i]))
		end
	end
	delta = zeros(eltype(coords), size(coords, 1))
	for (i, (node_a, node_b)) in enumerate(neighbors)
		for (node1, node2) in [(node_a, node_b), (node_b, node_a)]
			#add the advection term
			delta .= coords[:, node1] .- coords[:, node2]
			is_upwind = ((vxs[i] * delta[1] + vys[i] * delta[2] + vzs[i] * delta[3]) > 0)
			if (is_upwind && use_upwind) || (!is_upwind && !use_upwind)#!use_upwind basically means use_downwind and !is_upwind basically means is_downwind
				if !(node1 in dirichletnodes) && !(node2 in dirichletnodes)
					NonlinearEquations.addterm(node1, -(vxs[i] * delta[1] + vys[i] * delta[2] + vzs[i] * delta[3]) * (u[node1] - u[node2]) * areasoverlengths[i] / (volumes[node1]))
				elseif !(node1 in dirichletnodes) && node2 in dirichletnodes
					NonlinearEquations.addterm(node1, -(vxs[i] * delta[1] + vys[i] * delta[2] + vzs[i] * delta[3]) * (u[node1] - dirichletus[node2]) * areasoverlengths[i] / (volumes[node1]))
				end
			end
			#add the diffusion term
			if !(node1 in dirichletnodes) && !(node2 in dirichletnodes)
				NonlinearEquations.addterm(node1, Ds[i] * (u[node2] - u[node1]) * areasoverlengths[i] / (volumes[node1]))
			elseif !(node1 in dirichletnodes) && node2 in dirichletnodes
				NonlinearEquations.addterm(node1, Ds[i] * (dirichletus[node2] - u[node1]) * areasoverlengths[i] / (volumes[node1]))
			end
		end
	end
end


text:
using JLD2
using FileIO
using HDF5

# NOTE: we use jld2 instead of h5 becuase h5 does not easily build in a lot of
# environments (it takes a long time to compile or cannot becuase the
# appropriate compiler isn't available). jld2 builds anywhere that julia does.
cochba = h5open(joinpath("cochba.h5")) do file
  read(file,"/real") + read(file,"/imag")*im
end

jldopen(joinpath("cochba.jld2"),"w") do file
  M = size(cochba,2)
  for ch in 1:M
    p  = floor(Int,real(cochba[1, ch]))
    B  = real(cochba[(0:p).+2, ch])
    A  = imag(cochba[(0:p).+2, ch])
    file["filters/$(string(ch))/A"] = A
    file["filters/$(string(ch))/B"] = B
  end
  file["norm"] = imag(cochba[1, M])
end


text:
using Random: AbstractRNG
using Distributions
using DynamicPPL
import AbstractMCMC
import Distributions: logpdf

abstract type AbstractEvaluator end

struct SimpleModel{F} <: AbstractMCMC.AbstractModel
    f::F
end

(model::SimpleModel)(data...) = model.f(data...)

# sampling from prior: record all latent variables and their logpdf
struct SampleFromPrior <: AbstractEvaluator end

init_state(::SampleFromPrior) = (logp = Ref(0.0), latent = Dict{VarName, Tuple{<:Real, Float64}}())

function assume!(::SampleFromPrior, state, vn, dist)
    r = rand(dist)
    ℓ = logpdf(dist, r)
    push!(state.latent, vn => (r, ℓ))
    state.logp[] += ℓ
    return r, state
end

function observe!(::SampleFromPrior, state, vn, dist, obs)
    state.logp[] += logpdf(dist, obs)
    return state
end

# evaluataing the logpdf of a model given parameters θ
const VarDict{T} = Dict{VarName, T}

struct EvalLogp <: AbstractEvaluator
    θ::VarDict{<:Real}
end

init_state(::EvalLogp) = 0.0

function assume!(evaluator::EvalLogp, state, vn, dist)
    r = evaluator.θ[vn]
    state += logpdf(dist, r)
    return r, state
end

function observe!(evaluator::EvalLogp, state, vn, dist, obs)
    state += logpdf(dist, obs)
    return state
end

logpdf(model::SimpleModel, θ) = model(EvalLogp(θ))

# example sampler, like in https://turing.ml/dev/docs/using-turing/interface
# with static proposals
struct MetropolisHastings{D} <: AbstractMCMC.AbstractSampler
    init_θ::VarDict{<:Real}
    proposal::D
end

MetropolisHastings(init_θ::VarDict) =
    MetropolisHastings(init_θ, MvNormal(length(init_θ), 1))

# Create a very basic Transition type, only stores the
# parameter draws and the log probability of the draw.
struct Transition{L}
    θ::VarDict{<:Real}
    logp::L
end

# Store the new draw and its log density.
Transition(model::SimpleModel, θ::VarDict) = Transition(θ, logpdf(model, θ))
logpdf(model::SimpleModel, t::Transition) = t.logp

function AbstractMCMC.step!(
    rng::AbstractRNG,
    model::SimpleModel,
    spl::MetropolisHastings,
    N::Integer,
    ::Nothing;
    kwargs...
)
    return Transition(model, spl.init_θ)
end

# Define a function that makes a basic proposal.
function propose(rng, spl::MetropolisHastings, model::SimpleModel, θ::VarDict)
    # proposal = VarDict{Real}(keys(θ) .=> values(θ) .+ rand(rng, spl.proposal))
    proposal = VarDict{Real}(keys(θ) .=> rand(rng, spl.proposal))
    return Transition(model, proposal)
end
propose(rng, spl::MetropolisHastings, model::SimpleModel, t::Transition) =
    propose(rng, spl, model, t.θ)

# Calculates the probability `q(θ|θcond)`, using the proposal distribution `spl.proposal`.
q(spl::MetropolisHastings, θ::VarDict, θcond::VarDict) = logpdf(spl.proposal, values(θ) .- values(θcond))
q(spl::MetropolisHastings, t1::Transition, t2::Transition) = q(spl, t1.θ, t2.θ)

# Define the other step function. Returns a Transition containing
# either a new proposal (if accepted) or the previous proposal
# (if not accepted).
function AbstractMCMC.step!(
    rng::AbstractRNG,
    model::SimpleModel,
    spl::MetropolisHastings,
    ::Integer,
    t_prev::Transition;
    kwargs...
)
    # Generate a new proposal.
    t = propose(rng, spl, model, t_prev)

    # Calculate the log acceptance probability.
    log_α = logpdf(model, t) - logpdf(model, t_prev) + q(spl, t_prev, t) - q(spl, t, t_prev)

    # Decide whether to return the previous transition or the new one.
    if log(rand(rng)) < min(log_α, 0.0)
        return t
    else
        return t_prev
    end
end

# Example Model
# @model coinflip(y) = begin
#     p ~ Beta(1, 1)
#     N = length(y)
#     for i = 1:N
#         y[i] ~ Bernoulli(p)
#     end
# end

coinflip(y) = SimpleModel() do evaluator
    state = init_state(evaluator)
    p, state = assume!(evaluator, state, @varname(p), Beta(1, 1))

    N = length(y)
    for i = 1:N
        state = observe!(evaluator, state, @varname(y[i]), Bernoulli(p), y[i])
    end

    return state
end

model = coinflip([1,1,0])
println(model(SampleFromPrior()))
println(model(EvalLogp(VarDict{Real}(@varname(p) => 0.5))))
@show logpdf(Beta(1, 1), 0.5) + sum(logpdf.(Bernoulli(0.5), [1,1,0]))

# chain = AbstractMCMC.sample(model,
                            # MetropolisHastings(VarDict{Real}(@varname(p) => rand()),
                                               # Product([Uniform(0.3, 0.7)])),
                            # 10000)
# ps = [t.θ[@varname(p)] for t in chain]
# @show mean(ps)


text:
function internalsolve(mdp::MDP, pvi::ParallelValueIteration)

  solvechunk! = function emptyfunc() end
  if length(mdp.transition.argnames) == length(mdp.statemap) + length(mdp.actionmap)
    solvechunk! = solveset_chunk!
  elseif length(mdp.transition.argnames) > length(mdp.statemap) + length(mdp.actionmap)
    solvechunk! = solveprob_chunk!
  else
    error(string(
      "unkown transition function of neither T(s,a,s') nor T(s,a) type\n",
      "number of transition arguments: ",
      length(mdp.transition.argnames), "\n",
      "number of state and action variables",
      length(mdp.statemap) + length(mdp.actionmap)))
  end

  statedim = length(mdp.statemap)
  stateargs = mdp.reward.argnames[1:statedim]
  actionargs = mdp.reward.argnames[1 + statedim:end]

  nstates = length(pvi.stategrid)
  nactions = length(pvi.actiongrid)

  vold = SharedArray(
    Float64,
    (nstates),
    init = S -> S[localindexes(S)] = 0.0,
    pids = collect(1:pvi.nthreads))

  vnew = SharedArray(
    Float64,
    (nstates),
    init = S -> S[localindexes(S)] = 0.0,
    pids = collect(1:pvi.nthreads))

  qval = SharedArray(
    Float64,
    (nactions, nstates),
    init = S -> S[localindexes(S)] = 0.0,
    pids = collect(1:pvi.nthreads))

  lst = segment(pvi.nthreads, 1:nstates)

  iter = 0
  itertime = 0.0
  cputime = 0.0

  resid = -Inf

  for iter in 1:pvi.maxiter

    tic()
    resid = -Inf

    newresids = pmap(
      x -> (idxs = x; solvechunk!(
        mdp,
        pvi,
        nstates,
        nactions,
        stateargs,
        actionargs,
        vold,
        vnew,
        qval,
        idxs)),
      lst)

    newresid = maximum(newresids)
    newresid > resid ? resid = newresid : nothing
    results = pmap(x -> (idxs = x; updatechunk!(vold, vnew, idxs)), lst)

    vtmp = vold
    vold = vnew
    vnew = vtmp

    itertime = toq()
    cputime += itertime

    if pvi.verbose
      println(string(
        "iter $iter, resid: $resid, itertime: $itertime, cputime: $cputime"))
    end

    resid < pvi.tol ? break : nothing

  end

  if iter == pvi.maxiter
    @warn(string(
      "maximum number of iterations reached; check accuracy of solutions"))
  end

  @info(string(
    "value iteration solution generated\n",
    "cputime [s] = ", cputime, "\n",
    "number of iterations = ", iter, "\n",
    "residual = ", resid))

  return ValueIterationSolution(
    shmat2mat(qval),
    pvi.stategrid,
    pvi.actiongrid,
    cputime,
    iter,
    resid)

end

function solveset_chunk!(
    mdp::MDP,
    pvi::ParallelValueIteration,
    nstates::Int64,
    nactions::Int64,
    stateargs::Vector{String},
    actionargs::Vector{String},
    vold::SharedArray{Float64, 1},
    vnew::SharedArray{Float64, 1},
    qval::SharedArray{Float64, 2},
    istates::UnitRange)

  resid = -Inf

  for istate in istates

    state = getvar(pvi.stategrid, mdp.statemap, stateargs, istate)
    qhi = -Inf

    for iaction in 1:nactions

      action = getvar(pvi.actiongrid, mdp.actionmap, actionargs, iaction)

      statepIdxs, probs = transition(mdp, pvi, state, action, stateargs)
      qnow = 0.0

      for istatep in 1:length(statepIdxs)
        qnow += probs[istatep] * vold[statepIdxs[istatep]]
      end

      qnow *= pvi.discount
      qnow += reward(mdp, state, action)

      qval[iaction, istate] = qnow

      if qnow > qhi
        qhi = qnow
        vnew[istate] = qhi
      end

    end

    # use infinity-norm
    newresid = (vold[istate] - vnew[istate])^2
    newresid > resid ? resid = newresid : nothing

  end

  return resid

end

function solveprob_chunk!(
    mdp::MDP,
    pvi::ParallelValueIteration,
    nstates::Int64,
    nactions::Int64,
    stateargs::Vector{String},
    actionargs::Vector{String},
    vold::SharedArray{Float64, 1},
    vnew::SharedArray{Float64, 1},
    qval::SharedArray{Float64, 2},
    istates::UnitRange)

  resid = -Inf

  for istate in istates

    state = getvar(pvi.stategrid, mdp.statemap, stateargs, istate)
    qhi = -Inf

    for iaction in 1:nactions

      action = getvar(pvi.actiongrid, mdp.actionmap, actionargs, iaction)
      qnow = 0.0

      for istatep in 1:nstates
        statep = getvar(pvi.stategrid, mdp.statemap, stateargs, istatep)
        prob = mdp.transition.fn(state..., action..., statep...)
        qnow += prob * vold[istatep]
      end

      qnow *= pvi.discount
      qnow += reward(mdp, state, action)

      qval[iaction, istate] = qnow

      if qnow > qhi
        qhi = qnow
        vnew[istate] = qhi
      end

    end

    # use infinity-norm
    newresid = (vold[istate] - vnew[istate])^2
    newresid > resid ? (resid = newresid) : (nothing)

  end

  return resid

end

function updatechunk!(
    vold::SharedArray{Float64, 1},
    vnew::SharedArray{Float64, 1},
    idxs::UnitRange)

  for i in idxs
    vold[i] = vnew[i]
  end
  return idxs
end


text:
# CPXCALLSOLUTIONTYPE
const CPXCALLBACKSOLUTION_NOCHECK = Cint(-1)
const CPXCALLBACKSOLUTION_CHECKFEAS = Cint(0)
const CPXCALLBACKSOLUTION_PROPAGATE = Cint(1)
const CPXCALLBACKSOLUTION_SOLVE = Cint(2)

mutable struct CallbackContext
    model::Model
    context::Ptr{Cvoid}
end

function callback_wrapper(
    context::Ptr{Cvoid}, context_id::Clong, user_handle::Ptr{Cvoid}
)
    (model, callback_func) = unsafe_pointer_to_objref(user_handle)::Tuple{Model, Function}
    # Run the `callback_func`.
    callback_func(CallbackContext(model, context), context_id)
    # Exit normally.
    return Cint(0)
end

function return_status_or_throw(model, status)
    if status != 0
        throw(CplexError(model.env, status))
    else
        return status
    end
end

 """
    cpx_callbacksetfunc(model::Model, context_mask::Clong, callback_func::Function)

Set the general callback function of model to `callback_func`.

`callback_func` is a function that takes two arguments. The first is a
`CallbackContext` instance, and the second is a `Clong` indicating where the
callback is being called from.

### Examples

     function my_callback(cb_context::CPLEX.CallbackContext, context_id::Clong)
        CPLEX.cpx_callbackabort(cb_context)
    end
    CPLEX.cpx_callbacksetfunc(model, Clong(0), my_callback)
"""
function cbsetfunc(model::Model, context_mask::Clong, callback_func::Function)
    c_callback = @cfunction(callback_wrapper, Cint, (Ptr{Cvoid}, Clong, Ptr{Cvoid}))
    user_handle = (model, callback_func)::Tuple{Model, Function}
    return_status = @cpx_ccall(
        callbacksetfunc,
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}, Clong, Ptr{Cvoid}, Any),
        model.env, model.lp, context_mask, c_callback, user_handle
    )
    # We need to keep a reference to the callback function so that it isn't
    # garbage collected.
    model.callback = user_handle
    return return_status_or_throw(model, return_status)
end

function cbcandidateispoint(callback_data::CallbackContext)
    is_point = Ref{Cint}()
    stat = @cpx_ccall(
        callbackcandidateispoint,
        Cint,
        (Ptr{Cvoid}, Ptr{Cint}),
        callback_data.context, is_point
    )
    if stat != 0
        throw(CplexError(callback_data.model.env, stat))
    end
    return is_point[]
end

function cbgetcandidatepoint(
    callback_data::CallbackContext,
    x::Vector{Float64},
    begin_index::Cint,
    end_index::Cint,
    obj::Ref{Float64}
)
    return_status = @cpx_ccall(
        callbackgetcandidatepoint,
        Cint,
        (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint, Ref{Cdouble}),
        callback_data.context, x, begin_index, end_index, obj
    )
    return return_status_or_throw(callback_data.model, return_status)
end

function cbgetrelaxationpoint(
    callback_data::CallbackContext,
    x::Vector{Float64},
    begin_index::Cint,
    end_index::Cint,
    obj::Ref{Float64}
)
    return_status = @cpx_ccall(
        callbackgetrelaxationpoint,
        Cint,
        (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}),
        callback_data.context, x, begin_index, end_index, obj
    )
    return return_status_or_throw(callback_data.model, return_status)
end

function cbpostheursoln(
    callback_data::CallbackContext,
    cnt::Cint,
    ind::Vector{Cint},
    val::Vector{Cdouble},
    obj::Cdouble,
    strat::Cint
)
    return_status = @cpx_ccall(
        callbackpostheursoln,
        Cint,
        (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}, Cdouble, Cint),
        callback_data.context, cnt, ind, val, obj, strat
    )
    return return_status_or_throw(callback_data.model, return_status)
end

function cbaddusercuts(
    callback_data::CallbackContext,
    rcnt::Cint,
    nzcnt::Cint,
    rhs::Vector{Cdouble},
    sense::Vector{Cchar},
    rmatbeg::Vector{Cint},
    rmatind::Vector{Cint},
    rmatval::Vector{Cdouble},
    purgeable::Vector{Cint},
    _local::Vector{Cint}
)
    return_status = @cpx_ccall(
        callbackaddusercuts,
        Cint,
        (Ptr{Cvoid}, Cint, Cint, Ptr{Cdouble}, Ptr{Cchar}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}),
        callback_data.context, rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval, purgeable, _local
    )
    return return_status_or_throw(callback_data.model, return_status)
end

function cbrejectcandidate(
    callback_data::CallbackContext,
    rcnt::Cint,
    nzcnt::Cint,
    rhs::Vector{Cdouble},
    sense::Vector{Cchar},
    rmatbeg::Vector{Cint},
    rmatind::Vector{Cint},
    rmatval::Vector{Cdouble}
)
    return_status = @cpx_ccall(
        callbackrejectcandidate,
        Cint,
        (Ptr{Cvoid}, Cint, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}),
        callback_data.context, rcnt, nzcnt, rhs, sense, rmatbeg, rmatind, rmatval
    )
    return return_status_or_throw(callback_data.model, return_status)
end

function cbabort(callback_data::CallbackContext)
    return_status = @cpx_ccall(
        callbackabort, Cint, (Ptr{Cvoid},), callback_data.context
    )
    return return_status_or_throw(callback_data.model, return_status)
end

# # Check if the callback for the context CPX_CALLBACKCONTEXT_CANDIDATE is invoked by an unbounded relaxation
# function cbcandidateisray(callback_data::CallbackContext, is_ray::Int)
#     return_status = @cpx_ccall(callbackcandidateisray,
#         Cint, (Ptr{Cvoid}, Ptr{Cint}), callback_data.context, is_ray)
#     return return_status_or_throw(callback_data.model, return_status)
# end

# function cbgetcandidateray(
#         callback_data::CallbackContext, x::Vector{Cdouble}, begin_index::Int,
#         end_index::Int)
#     return_status = @cpx_ccall(callbackgetcandidateray,
#         Cint,
#         (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint),
#         callback_data.context, x, begin_index - 1, end_index - 1)
#     return return_status_or_throw(callback_data.model, return_status)
# end

# # This function may have bug, should be tested.
# # function cbgetfunc(model::Model, context_id::Clong, callback_func::Function)
# #     c_callback = @cfunction(callback_wrapper, Cint, (Ptr{Cvoid}, Clong, Ptr{Cvoid}))
# #     user_handle = (model, callback_func)::Tuple{Model, Function}
# #     return_status = @cpx_ccall(callbackgetfunc,
# #                     Cint,
# #                     (Ptr{Cvoid}, Ptr{Cvoid}, Ref{Clong}, Ptr{Cvoid}, Any),
# #                     model.env, model.lp, context_mask, c_callback, user_handle)
# #     if return_status != 0
# #         throw(CplexError(model.env, return_status))
# #     end
# #     model.callback = user_handle
# #     return return_status
# # end

# function cbgetincumbent(
#         callback_data::CallbackContext, x::Vector{Cdouble}, begin_index::Int,
#         end_index::Int, obj::Ref{Float64})
#     return_status = @cpx_ccall(callbackgetincumbent,
#         Cint,
#         (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint, Ref{Cdouble}),
#         callback_data.context, x, begin_index - 1, end_index - 1, obj)
#     return return_status_or_throw(callback_data.model, return_status)
# end

# function cbgetinfodbl(
#         callback_data::CallbackContext, cbinfo::CbInfo, data_p::Cdouble)
#     return_status = @cpx_ccall(callbackgetinfodbl,
#         Cint,
#         (Ptr{Cvoid}, Cint, Ref{Cdouble}),
#         callback_data.context, cbinfo, data_p)
#     return return_status_or_throw(callback_data.model, return_status)
# end

# function cbgetinfoint(
#         callback_data::CallbackContext, cbinfo::CbInfo, data_p::Int)
#     return_status = @cpx_ccall(callbackgetinfoint,
#         Cint,
#         (Ptr{Cvoid}, Cint, Ref{Cint}),
#         callback_data.context, cbinfo, data_p)
#     return return_status_or_throw(callback_data.model, return_status)
# end

# function cbgetinfolong(
#         callback_data::CallbackContext, cbinfo::CbInfo, data_p::Clong)
#     return_status = @cpx_ccall(callbackgetinfolong,
#         Cint,
#         (Ptr{Cvoid}, Cint, Ref{Clong}),
#         callback_data.context, cbinfo, data_p)
#     return return_status_or_throw(callback_data.model, return_status)
# end

# function cbgetlocallb(
#         callback_data::CallbackContext, lb::Vector{Cdouble}, begin_index::Int,
#         end_index::Int)
#     return_status = @cpx_ccall(callbackgetlocallb,
#         Cint,
#         (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint),
#         callback_data.context, lb, begin_index - 1, end_index - 1)
#     return return_status_or_throw(callback_data.model, return_status)
# end

# function cbgetlocalub(
#         callback_data::CallbackContext, ub::Vector{Cdouble}, begin_index::Int,
#         end_index::Int)
#     return_status = @cpx_ccall(callbackgetlocalub,
#         Cint,
#         (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint),
#         callback_data.context, ub, begin_index - 1, end_index - 1)
#     return return_status_or_throw(callback_data.model, return_status)
# end

# function cbgetgloballb(
#         callback_data::CallbackContext, lb::Vector{Cdouble}, begin_index::Int,
#         end_index::Int)
#     return_status = @cpx_ccall(callbackgetgloballb,
#         Cint,
#         (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint),
#         callback_data.context, lb, begin_index - 1, end_index - 1)
#     return return_status_or_throw(callback_data.model, return_status)
# end

# function cbgetglobalub(
#         callback_data::CallbackContext, ub::Vector{Cdouble}, begin_index::Int,
#         end_index::Int)
#     return_status = @cpx_ccall(callbackgetglobalub,
#         Cint,
#         (Ptr{Cvoid}, Ptr{Cdouble}, Cint, Cint),
#         callback_data.context, ub, begin_index - 1, end_index - 1)
#     return return_status_or_throw(callback_data.model, return_status)
# end


text:
@testset "Test utility functions" begin
    concrete_types = IS.get_all_concrete_subtypes(IS.InfrastructureSystemsType)
    @test length([x for x in concrete_types if isconcretetype(x)]) == length(concrete_types)
end

@testset "Test strip_module_name" begin
    @test IS.strip_module_name("PowerSystems.HydroDispatch") == "HydroDispatch"

    @test IS.strip_module_name(
        "InfrastructureSystems.Deterministic{PowerSystems.HydroDispatch}") ==
        "Deterministic{PowerSystems.HydroDispatch}"

    @test IS.strip_module_name(
        "Deterministic{PowerSystems.HydroDispatch}") ==
        "Deterministic{PowerSystems.HydroDispatch}"
end

@testset "Test strip_parametric_type" begin
    @test IS.strip_parametric_type(
        "Deterministic{PowerSystems.HydroDispatch}") == "Deterministic"

    @test IS.strip_parametric_type(
        "InfrastructureSystems.Deterministic{PowerSystems.HydroDispatch}") ==
        "InfrastructureSystems.Deterministic"
end

@testset "Test serialization utility functions" begin
    text = "SomeType{ParameterType1, ParameterType2}"
    type_str, parameters = IS.separate_type_and_parameter_types(text)
    @test type_str == "SomeType"
    @test parameters == ["ParameterType1", "ParameterType2"]

    text = "SomeType"
    type_str, parameters = IS.separate_type_and_parameter_types(text)
    @test type_str == "SomeType"
    @test parameters == []
end

@testset "Test exported names" begin
    @test IS.validate_exported_names(IS)
end


text:
##################### Define DeepPolyFreshVars Solver ##########################

@with_kw struct DeepPolyFreshVars <: Solver
    max_vars::Int64 = 10 # maximum number of fresh variables to introduce
    var_frac::Float64 = 0.5 # at most var_frac * n_neurons are introduced per layer
end

######################## Forward Propagation ###################################

function forward_linear(solver::DeepPolyFreshVars, L::LayerNegPosIdx, input::SymbolicIntervalFV)
    output_Low, output_Up = interval_map(L.W_neg, L.W_pos, input.sym.Low, input.sym.Up)
    output_Up[:, end] += L.bias
    output_Low[:, end] += L.bias
    sym = SymbolicInterval(output_Low, output_Up, domain(input))
    return SymbolicIntervalFV(sym, input.lbs, input.ubs, input.var_los, input.var_his, input.max_vars)
end

function forward_act(solver::DeepPolyFreshVars, L::LayerNegPosIdx{ReLU}, input::SymbolicIntervalFV)
    n_node = n_nodes(L)
    n_sym = size(input.sym.Low, 2) - 1 # last column is constant term
    n_in = dim(domain(input))
    current_n_vars = n_sym - n_in

    subs_sym_lo, subs_sym_hi = substitute_variables(input.sym.Low, input.sym.Up,
                                                    input.var_los, input.var_his,
                                                    n_in, current_n_vars)
    # TODO: write better function for bounds instead of this monstrosity!!!
    tmp_sym = SymbolicInterval(subs_sym_lo, subs_sym_hi, domain(input))
    los, his = bounds(tmp_sym, input.lbs[L.index], input.ubs[L.index])

    n_vars = min(input.max_vars - current_n_vars, floor(Int, solver.var_frac * n_node))
    fv_idxs = fresh_var_largest_range(los, his, n_vars)
    n_vars = length(fv_idxs)

    out_Low, out_Up = zeros(n_node, n_sym + 1 + n_vars), zeros(n_node, n_sym + 1 + n_vars)

    slopes = relaxed_relu_gradient.(los, his)
    out_Up[:, (1:n_sym) ∪ [end]] .= input.sym.Up .* slopes
    out_Up[:, end] .+= slopes .* max.(-los, 0)

    # also apply relaxation to substituted variables
    subs_sym_hi .= subs_sym_hi .* slopes
    subs_sym_hi[:, end] .+= slopes .* max.(-los, 0)

    out_Low[:, (1:n_sym) ∪ [end]] .= input.sym.Low .* relaxed_relu_gradient_lower.(los, his)

    # also apply relaxation to substituted variables
    subs_sym_lo .= subs_sym_lo .* relaxed_relu_gradient_lower.(los, his)

    for (i, v) in enumerate(fv_idxs)
        # store symbolic bounds on fresh variables
        input.var_los[current_n_vars + i, :] .= subs_sym_lo[v, :]
        input.var_his[current_n_vars + i, :] .= subs_sym_hi[v, :]

        # set corresponding entry to unit-vec
        out_Low[v,:] .= unit_vec(n_sym + i, n_sym + 1 + n_vars)
        out_Up[v,:]  .= unit_vec(n_sym + i, n_sym + 1 + n_vars)
    end

    sym = SymbolicInterval(out_Low, out_Up, domain(input))
    output = SymbolicIntervalFV(sym, input.lbs, input.ubs, input.var_los, input.var_his, input.max_vars)
    output.lbs[L.index] .= los
    output.ubs[L.index] .= his
    return output
end

function forward_act(solver::DeepPolyFreshVars, L::LayerNegPosIdx{Id}, input::SymbolicIntervalFV)
    return input
end

########################### Helper Functions ###################################

"""
Returns n indices of entries in lbs, ubs with largest range (ubs[i] - lbs[i])
"""
function fresh_var_largest_range(lbs::Vector{Float64}, ubs::Vector{Float64}, n::Int64)
    n <= 0 && return Int64[]  # early termination, if no variables are introduced

    ranges = ubs - lbs
    p = sortperm(-ranges) # sort descending
    p = p[(lbs[p] .< 0) .& (ubs[p] .> 0)] # want only crossing ReLUs
    return p[1:min(length(p), n)]
end

# TODO: remove after debugging!!!
function fresh_var_first(lbs::Vector{Float64}, ubs::Vector{Float64}, n::Int64)
    crossing_mask = ((lbs .< 0) .& (ubs .> 0))
    crossing = (1:length(lbs))[crossing_mask]
    return crossing[1:min(length(crossing), n)]
end

function unit_vec(i, n)
    e_i = zeros(n)
    e_i[i] = 1.
    return e_i
end


text:
function test_kernel_function(k)
    default_args, alt_args = get(kernel_functions_arguments, k, ((), ()))

    local n = length(default_args)
    K = (k)()

    @testset "Testing constructors" begin
        for j = 1:n
            @test getfield(K,j) == default_args[j]
        end
        @test eltype(K) == Float64

        for T in FloatingPointTypes, i = 1:n
            K = (k)([T(a) for a in alt_args[1:i]]...)
            for j = 1:n
                @test getfield(K,j) == (j <= i ? alt_args[j] : default_args[j])
            end
            @test eltype(K) == T
        end
    end

    @testset "Testing properties" begin
        @test MLK.ismercer(K) == isa(K, MLK.MercerKernel)
        @test MLK.isnegdef(K) == isa(K, MLK.NegativeDefiniteKernel)
        @test MLK.isstationary(K) == MLK.isstationary(MLK.basefunction(K))
        @test MLK.isisotropic(K) == MLK.isisotropic(MLK.basefunction(K))
    end

    @testset "Testing conversions" begin
        psi = k
        while psi != Any
            for T1 in FloatingPointTypes, T2 in FloatingPointTypes
                K1 = k{T1}()
                K2 = convert(psi{T2}, K1)
                @test eltype(K2) == T2
                @test isa(K2, k)
            end
            psi = supertype(psi)
        end
    end

    @testset "Testing kappa function" begin
        f = get(kernel_functions_kappa, k, x->error(""))
        args1, args2 = get(kernel_functions_arguments, k, ((), ()))
        for i = 0:length(args1)
            args = [j <= i ? args1[j] : args2[j] for j in eachindex(args1)]
            for T in FloatingPointTypes
                K = k{T}(args...)
                for z in [T(0), T(1), T(2)]
                    v1 = MLK.kappa(K, z)
                    v2 = (f)(z, args...)
                    @test isapprox(v1, v2)
                end
            end
        end
    end

    @testset "Testing show function" begin
        @test show(devnull, K) == nothing
    end
end

for k in kernel_functions
    @testset "Testing $k" begin
        test_kernel_function(k)
    end
end


text:
type Errors
    a::Float64 # parameter distribution smoothing
    h::Float64

    odev::Vector{Float64} # σ for observation, param distributions
    pdev::Vector{Float64}

    lb::Vector{Float64} # lower/upper bounds for truncated normal param distrs.
    ub::Vector{Float64}

    function Errors(nparams=9)
        this = new();
        δ = 0.99; # ∃ (0, 1, typically within 0.95-0.99), lower values = less parameter dispersion
        this.a = (3*δ-1)/(2*δ);
        this.h = sqrt.(1-this.a^2);
        this.odev = [2.5*mmHgToPa,2.5*cm3Tom3,2.5*cm3Tom3];
        this.pdev = [1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12]; # no params
        # this.pdev = [2e4,3e6,1e-9,1e-4,0.025,0.04,0.2,2,1.4e7]; # all params
        # this.pdev = [2e4,2e6,1e-9,1e-12,0.1,0.02,0.2,2,7e6]; # all params - V0
        # this.pdev = [2e4,3e6,1e-9,1e-4,0.025,1e-12,0.2,2,1.4e7]; # all params - t2
        # this.pdev = [2e4,2e6,1e-12,1e-12,0.1,0.1,0.1,2.7,1e-12]; # all params - C - V0 - Emax
        # this.pdev = [2e4,2e6,1e-12,1e-12,0.1,0.1,0.1,2.7,7e6]; # all params - C - V0
        # this.pdev = [2e4,2e6,1e-9,1e-12,1e-12,1e-12,1e-12,1e-12,1e-5]; # Windkessel only
        # this.pdev = [2e4,2e6,1e-9,1e-12,1e-12,1e-12,1e-12,1e-12,1.4e7]; # Windkessel + Emax
        # this.pdev = [2e4,1e-12,1e-9,1e-12,1e-12,1e-12,1e-12,1e-12,7e6]; # Windkessel - R + Emax
        # this.pdev = [2e4,2e6,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,7e6]; # Windkessel - C + Emax
        # this.pdev = [2e4,2e6,1e-9,1e-12,1e-12,0.1,1e-12,1e-12,1.4e7]; # Windkessel + Emax + τ2
        # this.pdev = [2e4,2e6,1e-8,1e-12,1e-12,1e-12,1e-12,2.7,1e-12]; # Windkessel + m2
        # this.pdev = [2e4,2e6,1e-8,1e-12,1e-12,1e-12,1e-12,2.7,7e6]; # Windkessel + Emax + m2
        # this.pdev = [2e4,2e6,1e-12,1e-12,1e-12,1e-12,1e-12,2.7,7e6]; # Windkessel - C + Emax + m2
        # this.pdev = [2e4,2e6,1e-9,1e-12,1e-12,1e-12,0.2,2,1.4e7]; # Windkessel + Emax + m1 + m2
        # this.pdev = [2e4,2e6,1e-9,1e-12,0.1,1e-12,0.2,2,7e6]; # Windkessel + Emax + t1 + m1 + m2
        # this.pdev = [2e4,1e-12,1e-9,1e-12,1e-12,1e-12,0.1,2.7,7e6]; # Windkessel - R + Emax + m1 + m2
        # this.pdev = [2e4,1e-12,1e-9,1e-12,1e-2,1e-12,0.1,2.7,7e6]; # Windkessel - R + Emax + t1 + m1 + m2
        # this.pdev = [2e4,2e6,1e-12,1e-12,1e-12,1e-12,0.1,2.7,7e6]; # Windkessel - C + Emax + m1 + m2
        # this.pdev = [2e4,2e6,1e-12,1e-12,0.2,1e-12,0.1,2.7,7e6]; # Windkessel - C + Emax + m1 + m2 + t1
        # this.pdev = [2e4,2e6,1e-12,1e-12,1e-12,0.2,0.1,2.7,7e6]; # Windkessel - C + Emax + m1 + m2 + Δt
        # this.pdev = [2e4,2e6,1e-12,1e-12,1e-12,1e-12,0.1,2.7,1e-12]; # Windkessel - C + m1 + m2
        # this.pdev = [1e-12,2e6,1e-9,1e-12,1e-12,1e-12,0.1,1e-12,7e6]; # Windkessel - Zc + Emax + m1
        # this.pdev = [1e-12,2e6,1e-12,1e-12,1e-12,1e-12,0.1,2.7,7e6]; # R + Emax + m1 + m2
        # this.pdev = [1e-12,2e6,1e-12,1e-12,1e-12,3e-1,0.1,2.7,7e6]; # R + Emax + m1 + m2 + Δτ
        # this.pdev = [1e-12,1e-12,1e-12,1e-12,0.1,0.1,0.2,2.7,1e-12]; # t1 + t2 + m1 + m2
        # this.pdev = [2e4,1e-12,1e-12,1e-12,1e-12,1e-12,0.1,2.7,7e6]; # Zc + Emax + m1 + m2
        # this.pdev = [1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,0.1,2.7,7e6]; # Emax + m1 + m2
        # this.pdev = [2e4,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,7e6]; # Zc + Emax
        # this.pdev = [1e-12,2e6,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,7e6]; # R + Emax
        # this.pdev = [1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,0.1,1e-12,7e6]; # Emax + m1
        # this.pdev = [1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,2.7,7e6]; # Emax + m2
        # this.pdev = [1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,7e6]; # Emax
        # this.pdev = [1e-12,1e-12,1e-12,1e-4,1e-12,1e-12,1e-12,1e-12,1e-12]; # V0
        this.lb = -Inf*ones(nparams);
        this.ub = Inf*ones(nparams);
        this.lb[1] = 1e3; # Zc
        this.lb[2] = 1e5; # R
        this.lb[3] = 5e-9; # C
        this.lb[5] = 0.01; # τ1
        this.lb[6] = 0.01; # τ2
        this.lb[7] = 1.; # m1
        this.lb[8] = 1; # m2
        this.lb[9] = 1e6; # Emax
        # this.ub[5] = 0.4; # τ1
        # this.ub[6] = 0.55; # Δτ
        this.ub[9] = 2e8; # Emax
        return this
    end
end


text:
using GeometricalPredicates
using Test

@testset "2D orientation" begin
    @testset "Points" begin
        p = Point(1,1)
        @test getx(p) === 1.0
        @test gety(p) === 1.0

        p = Point(1, 1, 1)
        @test getx(p) === 1.0
        @test gety(p) === 1.0
        @test getz(p) === 1.0
    end

    @testset "Lines" begin
        a = Point(1.1, 1.1)
        b = Point(1.5, 1.5)

        l = Line(a, b)
        @test orientation(l, Point(1.4, 1.6)) == 1
        @test orientation(l, Point(1.4,1.05)) == -1
        @test orientation(l, Point(1.4,1.4)) == 0

        l = Line(b, a)
        @test orientation(l, Point(1.4, 1.6)) == -1
        @test orientation(l, Point(1.4,1.05)) == 1
        @test orientation(l, Point(1.4,1.4)) == 0

        @test abs(length2(l)-0.4*0.4*2) < 1e-5

        a = Point(1, 1)
        b = Point(2, 2)

        l = Line(a, b)
        @test orientation(l, Point(1.4, 1.6)) == 1
        @test orientation(l, Point(1.4,1.05)) == -1
        @test orientation(l, Point(1.4,1.4)) == 0

        l = Line(b, a)
        @test orientation(l, Point(1.4, 1.6)) == -1
        @test orientation(l, Point(1.4,1.05)) == 1
        @test orientation(l, Point(1.4,1.4)) == 0

        @test abs(length2(l)-2) < 1e-5
    end

    @testset "Triangles" begin
        ax = 1.1; ay = 1.1
        bx = 1.2; by = 1.2
        cx = 1.3; cy = 1.3
        @test orientation(ax,ay,bx,by,cx,cy) == 0
        @test orientation(bx,by,ax,ay,cx,cy) == 0
        @test orientation(bx,by,cx,cy,ax,ay) == 0

        ax = 1.1; ay = 1.1
        bx = 1.2; by = 1.199999999999
        cx = 1.3; cy = 1.3
        @test orientation(ax,ay,bx,by,cx,cy) == 1
        @test orientation(bx,by,ax,ay,cx,cy) == -1
        @test orientation(bx,by,cx,cy,ax,ay) == 1

        ax = 1.1; ay = 1.1
        bx = 1.2; by = 1.199999999999999
        cx = 1.3; cy = 1.3
        @test orientation(ax,ay,bx,by,cx,cy) == 1
        @test orientation(bx,by,ax,ay,cx,cy) == -1
        @test orientation(bx,by,cx,cy,ax,ay) == 1

        ax = 1.1; ay = 1.1
        bx = 1.2; by = 1.21
        cx = 1.3; cy = 1.3
        @test orientation(ax,ay,bx,by,cx,cy) == -1
        @test orientation(bx,by,ax,ay,cx,cy) == 1
        @test orientation(bx,by,cx,cy,ax,ay) == -1

        ax = 1.1; ay = 1.1
        bx = 1.2; by = 1.20000000000001
        cx = 1.3; cy = 1.3
        @test orientation(ax,ay,bx,by,cx,cy) == -1
        @test orientation(bx,by,ax,ay,cx,cy) == 1
        @test orientation(bx,by,cx,cy,ax,ay) == -1

        p1 = Point(1, 1)
        p2 = Point2D(1.9, 1.5)
        p3 = Point2D(1.45, 1.25)
        tr = Triangle(p1, p2, p3)
        @test orientation(tr) == 0

        p3 = Point2D(getx(p3), 1.3)
        tr = Triangle(p1, p2, p3)
        @test orientation(tr) == 1

        p3 = Point2D(getx(p3), 1.2)
        tr = Triangle(p1, p2, p3)
        @test orientation(tr) == -1

        tr = Triangle(p1, p2, p3, positivelyoriented)
        @test orientation(tr) == 1

        tr = Triangle(p1, p3, p2, negativelyoriented)
        @test orientation(tr) == -1
    end
end

@testset "2D incircle" begin
    ax = 1.1; ay = 1.1
    bx = 1.2; by = 1.2
    cx = 1.3; cy = 1.3
    dx = 1.4; dy = 1.7
    @test incircle(ax,ay,bx,by,cx,cy,dx,dy) == 2
    ax = 1; ay = 1
    @test incircle(ax,ay,bx,by,cx,cy,dx,dy) == 2

    ax = 1.1; ay = 1.1
    bx = 1.3; by = 1.1
    cx = 1.3; cy = 1.3
    dx = 1.1; dy = 1.3
    @test incircle(ax,ay,bx,by,cx,cy,dx,dy) == 0
    @test incircle(bx,by,ax,ay,cx,cy,dx,dy) == 0
    @test incircle(bx,by,cx,cy,ax,ay,dx,dy) == 0
    ax = 1; ay = 1
    bx = 1.3; by = 1
    cx = 1.3; cy = 1.3
    dx = 1; dy = 1.3
    @test incircle(ax,ay,bx,by,cx,cy,dx,dy) == 0
    @test incircle(bx,by,ax,ay,cx,cy,dx,dy) == 0
    @test incircle(bx,by,cx,cy,ax,ay,dx,dy) == 0

    ax = 1.1; ay = 1.1
    bx = 1.3; by = 1.1
    cx = 1.3; cy = 1.3
    dx = 1.1; dy = 1.3000001
    @test incircle(ax,ay,bx,by,cx,cy,dx,dy) == -1
    @test incircle(bx,by,ax,ay,cx,cy,dx,dy) == -1
    @test incircle(bx,by,cx,cy,ax,ay,dx,dy) == -1

    ax = 1.1; ay = 1.1
    bx = 1.3; by = 1.1
    cx = 1.3; cy = 1.3
    dx = 1.1; dy = 1.30000000000001
    @test incircle(ax,ay,bx,by,cx,cy,dx,dy) == -1
    @test incircle(bx,by,ax,ay,cx,cy,dx,dy) == -1
    @test incircle(bx,by,cx,cy,ax,ay,dx,dy) == -1

    ax = 1.1; ay = 1.1
    bx = 1.3; by = 1.1
    cx = 1.3; cy = 1.3
    dx = 1.1; dy = 1.29999
    @test incircle(ax,ay,bx,by,cx,cy,dx,dy) == 1
    @test incircle(bx,by,ax,ay,cx,cy,dx,dy) == 1
    @test incircle(bx,by,cx,cy,ax,ay,dx,dy) == 1

    ax = 1.1; ay = 1.1
    bx = 1.3; by = 1.1
    cx = 1.3; cy = 1.3
    dx = 1.1; dy = 1.29999999999999
    @test incircle(bx,by,ax,ay,cx,cy,dx,dy) == 1
    @test incircle(bx,by,ax,ay,cx,cy,dx,dy) == 1
    @test incircle(bx,by,cx,cy,ax,ay,dx,dy) == 1

    p1 = Point(1, 1)
    p2 = Point2D(1.0625, 1.0)
    p3 = Point2D(1.0625, 1.0625)
    tr = Triangle(p1, p2, p3)

    p4 = Point2D(1.03, 1.003)
    p5 = Point2D(1.99, 1.99)
    p6 = Point2D(1.0, 1.0625)
    @test incircle(tr, p4) == 1
    @test incircle(tr, p5) == -1
    @test incircle(tr, p6) == 0
    @test incircle(tr, p2) == 0
    @test incircle(tr, p3) == 0
    @test incircle(tr, p1) == 0
end

@testset "3D orientation" begin
    ax = 1.1; ay = 1.1; az=1.1
    bx = 1.2; by = 1.2; bz=1.1
    cx = 1.3; cy = 1.3; cz=1.1
    dx = 1.4; dy = 1.4; dz=1.1
    @test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == 0
    @test orientation(dx,dy,dz,bx,by,bz,cx,cy,cz,ax,ay,az) == 0
    @test orientation(dx,dy,dz,cx,cy,cz,bx,by,bz,ax,ay,az) == 0

    ax = 1.1; ay = 1.1; az=1.1
    bx = 1.2; by = 1.2; bz=1.1
    cx = 1.3; cy = 1.3; cz=1.1
    dx = 1.4; dy = 1.4; dz=1.100001
    @test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == 0
    @test orientation(dx,dy,dz,bx,by,bz,cx,cy,cz,ax,ay,az) == 0
    @test orientation(dx,dy,dz,cx,cy,cz,bx,by,bz,ax,ay,az) == 0

    ax = 1.1; ay = 1.1; az=1.1
    bx = 1.2; by = 1.2; bz=1.1
    cx = 1.3; cy = 1.3; cz=1.1
    dx = 1.4; dy = 1.4; dz=1.100000000001
    @test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == 0
    @test orientation(dx,dy,dz,bx,by,bz,cx,cy,cz,ax,ay,az) == 0
    @test orientation(dx,dy,dz,cx,cy,cz,bx,by,bz,ax,ay,az) == 0

    ax = 1.1; ay = 1.1; az=1.1
    bx = 1.2; by = 1.2; bz=1.1
    cx = 1.3; cy = 1.15; cz=1.1
    dx = 1.4; dy = 1.17; dz=1.1
    @test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == 0
    @test orientation(dx,dy,dz,bx,by,bz,cx,cy,cz,ax,ay,az) == 0
    @test orientation(dx,dy,dz,cx,cy,cz,bx,by,bz,ax,ay,az) == 0

    ax = 1.1; ay = 1.1; az=1.1
    bx = 1.3; by = 1.1; bz=1.1
    cx = 1.3; cy = 1.3; cz=1.1
    dx = 1.2; dy = 1.2; dz=1.5
    @test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == 1
    @test orientation(dx,dy,dz,bx,by,bz,cx,cy,cz,ax,ay,az) == -1
    @test orientation(dx,dy,dz,cx,cy,cz,bx,by,bz,ax,ay,az) == 1

    ax = 1.1; ay = 1.1; az=1.1
    bx = 1.3; by = 1.1; bz=1.1
    cx = 1.3; cy = 1.3; cz=1.1
    dx = 1.2; dy = 1.2; dz=1.10000000000001
    @test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == 1
    @test orientation(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz) == -1

    ax = 1.1; ay = 1.1; az=1.1
    bx = 1.3; by = 1.1; bz=1.1
    cx = 1.3; cy = 1.3; cz=1.1
    dx = 1.2; dy = 1.2; dz=1.09
    @test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == -1
    @test orientation(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz) == 1
    @test orientation(ax,ay,az,bx,by,bz,dx,dy,dz,cx,cy,cz) == 1

    ax = 1.1; ay = 1.1; az=1.1
    bx = 1.3; by = 1.1; bz=1.1
    cx = 1.3; cy = 1.3; cz=1.1
    dx = 1.2; dy = 1.2; dz=1.099999999999999
    @test orientation(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz) == -1
    @test orientation(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz) == 1
    @test orientation(ax,ay,az,bx,by,bz,dx,dy,dz,cx,cy,cz) == 1

    p1 = Point3D(ax, ay, az)
    p2 = Point3D(bx, by, bz)
    p3 = Point3D(cx, cy, cz)
    p4 = Point3D(dx, dy, dz)

    t = Tetrahedron(p1, p2, p3, p4)
    @test orientation(t) == -1

    t = Tetrahedron(p1, p2, p3, p4, positivelyoriented)
    @test orientation(t) == 1

    t = Tetrahedron(p1, p2, p4, p3, negativelyoriented)
    @test orientation(t) == -1
end

@testset "3D incircle" begin
    ax = 1.1; ay = 1.1; az=1.1
    bx = 1.3; by = 1.1; bz=1.1
    cx = 1.3; cy = 1.3; cz=1.1
    dx = 1.4; dy = 1.4; dz=1.2
    ex = 1.1; ey = 1.1; ez=1.1
    @test incircle(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,ex,ey,ez) == 0
    @test incircle(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz,ex,ey,ez) == 0
    @test incircle(ax,ay,az,bx,by,bz,dx,dy,dz,cx,cy,cz,ex,ey,ez) == 0

    ax = 1.1; ay = 1.1; az=1.1
    bx = 1.3; by = 1.1; bz=1.1
    cx = 1.3; cy = 1.3; cz=1.1
    dx = 1.2; dy = 1.2; dz=1.5
    ex = 1.2; ey = 1.15; ez=1.13
    @test incircle(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,ex,ey,ez) == 1
    @test incircle(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz,ex,ey,ez) == 1
    @test incircle(ax,ay,az,bx,by,bz,dx,dy,dz,cx,cy,cz,ex,ey,ez) == 1

    ax = 1.1; ay = 1.1; az=1.1
    bx = 1.3; by = 1.1; bz=1.1
    cx = 1.3; cy = 1.3; cz=1.1
    dx = 1.2; dy = 1.2; dz=1.5
    ex = 1.2; ey = 1.15; ez=1.1000000000001
    @test incircle(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,ex,ey,ez) == 1
    @test incircle(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz,ex,ey,ez) == 1
    @test incircle(ax,ay,az,bx,by,bz,dx,dy,dz,cx,cy,cz,ex,ey,ez) == 1

    ax = 1.1; ay = 1.1; az=1.1
    bx = 1.3; by = 1.1; bz=1.1
    cx = 1.3; cy = 1.3; cz=1.1
    dx = 1.2; dy = 1.2; dz=1.5
    ex = 1.9; ey = 1.15; ez=1.01
    @test incircle(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz,ex,ey,ez) == -1
    @test incircle(bx,by,bz,ax,ay,az,cx,cy,cz,dx,dy,dz,ex,ey,ez) == -1
    @test incircle(ax,ay,az,bx,by,bz,dx,dy,dz,cx,cy,cz,ex,ey,ez) == -1
end

@testset "intriangle (2D)!" begin
    ax = 1.1; ay = 1.1
    bx = 1.4; by = 1.1
    cx = 1.1; cy = 1.4
    dx = 1.11; dy = 1.11
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == 1
    @test intriangle(bx,by,ax,ay,cx,cy,dx,dy) == 1
    @test intriangle(bx,by,cx,cy,ax,ay,dx,dy) == 1

    ax = 1.1; ay = 1.1
    bx = 1.4; by = 1.1
    cx = 1.1; cy = 1.4
    dx = 1.01; dy = 1.01
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) < 0
    @test intriangle(bx,by,ax,ay,cx,cy,dx,dy) < 0
    @test intriangle(bx,by,cx,cy,ax,ay,dx,dy) < 0

    ax = 1.1; ay = 1.1
    bx = 1.4; by = 1.1
    cx = 1.1; cy = 1.4
    dx = 1.31; dy = 1.01
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) < 0
    @test intriangle(bx,by,ax,ay,cx,cy,dx,dy) < 0
    @test intriangle(bx,by,cx,cy,ax,ay,dx,dy) < 0

    ax = 1.1; ay = 1.1
    bx = 1.4; by = 1.1
    cx = 1.1; cy = 1.4
    dx = 1.1; dy = 1.2
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == 3
    @test intriangle(bx,by,ax,ay,cx,cy,dx,dy) == 2
    @test intriangle(bx,by,cx,cy,ax,ay,dx,dy) == 2

    ax = 1.1; ay = 1.1
    bx = 1.4; by = 1.1
    cx = 1.1; cy = 1.4
    dx = 1.2; dy = 1.1
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == 4
    @test intriangle(bx,by,ax,ay,cx,cy,dx,dy) == 4
    @test intriangle(bx,by,cx,cy,ax,ay,dx,dy) == 3

    ax = 1.1; ay = 1.1
    bx = 1.4; by = 1.1
    cx = 1.1; cy = 1.4
    dx = 1.25; dy = 1.25
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == 2
    @test intriangle(bx,by,ax,ay,cx,cy,dx,dy) == 3
    @test intriangle(bx,by,cx,cy,ax,ay,dx,dy) == 4

    ax = 1.1; ay = 1.1
    bx = 1.4; by = 1.1
    cx = 1.1; cy = 1.4
    dx = 1.25; dy = 1.250000000000001
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) < 0
    @test intriangle(bx,by,ax,ay,cx,cy,dx,dy) < 0
    @test intriangle(bx,by,cx,cy,ax,ay,dx,dy) < 0

    ax = 1.1; ay = 1.1
    bx = 1.4; by = 1.1
    cx = 1.1; cy = 1.4
    dx = 1.25; dy = 1.249999999999999
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == 1
    @test intriangle(bx,by,ax,ay,cx,cy,dx,dy) == 1
    @test intriangle(bx,by,cx,cy,ax,ay,dx,dy) == 1

    ax = 1.1; ay = 1.1
    bx = 1.4; by = 1.1
    cx = 1.1; cy = 1.4
    dx = 1.7; dy = 1.1
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -1
    dx = 1.1; dy = 1.7
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -1
    dx = 1.01; dy = 1.1
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -2
    dx = 1.1; dy = 1.01
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -3

    ax = 1.1; ay = 1.1
    bx = 1.4; by = 1.1
    cx = 1.1; cy = 1.4
    dx = 1.55; dy = 1.55
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -1
    dx = 1.09; dy = 1.11
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -2
    dx = 1.11; dy = 1.09
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -3

    ax = min_coord; ay = min_coord
    bx = 1.1; by = 1.1
    cx = max_coord; cy = min_coord
    dx = 1.2; dy = 1.2
    @test intriangle(ax,ay,bx,by,cx,cy,dx,dy) == -1
end

@testset "intriangle (3D) !" begin
    ax = 1.1; ay = 1.1; az = 1.1
    bx = 1.4; by = 1.1; bz = 1.1
    cx = 1.1; cy = 1.4; cz = 1.1
    dx = 1.1; dy = 1.1; dz = 1.4
    px = 1.11; py = 1.11; pz = 1.11
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == 1
    @test intriangle(bx, by, bz, ax, ay, az, cx, cy, cz, dx, dy, dz, px, py, pz) == 1
    @test intriangle(ax, ay, az, cx, cy, cz, bx, by, bz, dx, dy, dz, px, py, pz) == 1

    ax = 1.1; ay = 1.1; az = 1.1
    bx = 1.4; by = 1.1; bz = 1.1
    cx = 1.1; cy = 1.4; cz = 1.1
    dx = 1.1; dy = 1.1; dz = 1.4
    px = 1.02; py = 1.2; pz = 1.2
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) < 0
    @test intriangle(bx, by, bz, ax, ay, az, cx, cy, cz, dx, dy, dz, px, py, pz) < 0
    @test intriangle(ax, ay, az, cx, cy, cz, bx, by, bz, dx, dy, dz, px, py, pz) < 0

    ax = 1.1; ay = 1.1; az = 1.1
    bx = 1.4; by = 1.1; bz = 1.1
    cx = 1.1; cy = 1.4; cz = 1.1
    dx = 1.1; dy = 1.1; dz = 1.4
    px = 1.09999999999999; py = 1.2; pz = 1.2
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) < 0
    @test intriangle(bx, by, bz, ax, ay, az, cx, cy, cz, dx, dy, dz, px, py, pz) < 0
    @test intriangle(ax, ay, az, cx, cy, cz, bx, by, bz, dx, dy, dz, px, py, pz) < 0

    ax = 1.1; ay = 1.1; az = 1.1
    bx = 1.4; by = 1.1; bz = 1.1
    cx = 1.1; cy = 1.4; cz = 1.1
    dx = 1.1; dy = 1.1; dz = 1.4
    px = 1.10000000000001; py = 1.2; pz = 1.2
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == 1
    @test intriangle(bx, by, bz, ax, ay, az, cx, cy, cz, dx, dy, dz, px, py, pz) == 1
    @test intriangle(ax, ay, az, cx, cy, cz, bx, by, bz, dx, dy, dz, px, py, pz) == 1

    ax = 1.1; ay = 1.1; az = 1.1
    bx = 1.4; by = 1.1; bz = 1.1
    cx = 1.1; cy = 1.4; cz = 1.1
    dx = 1.1; dy = 1.1; dz = 1.4
    px = 1.11; py = 1.11; pz = 1.1
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == 5
    @test intriangle(bx, by, bz, ax, ay, az, cx, cy, cz, dx, dy, dz, px, py, pz) == 5
    @test intriangle(ax, ay, az, cx, cy, cz, bx, by, bz, dx, dy, dz, px, py, pz) == 5
    @test intriangle(dx, dy, dz, bx, by, bz, cx, cy, cz, ax, ay, az, px, py, pz) == 2
    @test intriangle(ax, ay, az, dx, dy, dz, cx, cy, cz, bx, by, bz, px, py, pz) == 3
    @test intriangle(ax, ay, az, bx, by, bz, dx, dy, dz, cx, cy, cz, px, py, pz) == 4

    ax = 1.1; ay = 1.1; az = 1.1
    bx = 1.4; by = 1.1; bz = 1.1
    cx = 1.1; cy = 1.4; cz = 1.1
    dx = 1.1; dy = 1.1; dz = 1.4
    px = 1.55; py = 1.55; pz = 1.55
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
    px = 1.09; py = 1.11; pz = 1.11
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -2
    px = 1.11; py = 1.09; pz = 1.11
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -3
    px = 1.11; py = 1.11; pz = 1.09
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -4

    ax = 1.1; ay = 1.1; az = 1.1
    bx = 1.4; by = 1.1; bz = 1.1
    cx = 1.1; cy = 1.4; cz = 1.1
    dx = 1.1; dy = 1.1; dz = 1.4
    px = 1.7; py = 1.1; pz = 1.1
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
    px = 1.7; py = 1.1; pz = 1.1
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
    px = 1.1; py = 1.7; pz = 1.1
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
    px = 1.1; py = 1.1; pz = 1.7
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
    px = 1.9; py = 1.9; pz = 1.1
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
    px = 1.9; py = 1.1; pz = 1.9
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
    px = 1.1; py = 1.9; pz = 1.9
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -1
    px = 1.01; py = 1.1; pz = 1.1
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -2
    px = 1.1; py = 1.01; pz = 1.1
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -3
    px = 1.1; py = 1.1; pz = 1.01
    @test intriangle(ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, px, py, pz) == -4
end

@testset "Scale dependednt Peano-Hilbert stuff" begin
    @testset "2D" begin
        @test peanokey(Point2D(1.01,1.01), 2) == 0
        @test peanokey(Point2D(1.9901,1.01), 2) == 4*4-1
        @test Point2D(15, 2) == Point2D(1.75, 1.0)
        @test Point2D(0, 2) == Point(1, 1)
        for x in range(1.0,stop=1.999999,length=100), y in range(1.0,stop=1.999999,length=100)
            p = Point2D(x, y)
            d = peanokey(p)
            pp= Point2D(d)
            @test abs(getx(p) - getx(pp)) < 1e-7
            @test abs(gety(p) - gety(pp)) < 1e-7
        end
    end

    @testset "3D" begin
        @test peanokey(Point3D(1.01,1.01,1.01), 2) == 0
        @test peanokey(Point3D(1.01,1.01,1.9901), 2) == 4*4*4-1
        @test Point3D(15, 2) == Point3D(1.25,1.5,1.0)
        @test Point3D(0, 2) == Point3D(1.0, 1.0, 1.0)
        for x in range(1.0,stop=1.99999,length=30), y in range(1.0,stop=1.99999,length=30), z in range(1.0,stop=1.999999,length=30)
            p = Point3D(x, y, z)
            d = peanokey(p)
            pp= Point3D(d)
            @test abs(getx(p) - getx(pp)) < 1e-5
            @test abs(gety(p) - gety(pp)) < 1e-5
            @test abs(getz(p) - getz(pp)) < 1e-5
        end
    end
end

@testset "Scale independent Peano-Hilbert stuff" begin
    a2 = Point2D[]
    N=1000
    for i in 1:N
        push!(a2, Point2D(1.0+rand(), 1.0+rand()))
    end
    # TODO: not much of a test, I need to betteer think how to test this
    # At least this runs the code...
    @test length(mssort!(a2)) == N

    a3 = Point3D[]
    for i in 1:N
        push!(a3, Point3D(1.0+rand(), 1.0+rand(), 1.0+rand()))
    end
    # TODO: not much of a test, I need to betteer think how to test this
    # At least this runs the code...
    @test length(mssort!(a3)) == N
end

@testset "Qttys functions" begin
    @testset "2D" begin
        a = Point(1, 1)
        b = Point2D(1.0, 1.5)
        c = Point2D(1.5, 1.0)
        t = Triangle(a, b, c)
        @test abs(area(t)-1/8) < 1e-7
        c = centroid(t)
        @test abs(getx(c)-1-1/6) < 1e-7
        @test abs(gety(c)-1-1/6) < 1e-7
        c = circumcenter(t)
        @test abs(getx(c)-1.25) < 1e-7
        @test abs(gety(c)-1.25) < 1e-7
        @test abs(circumradius2(t)-1/8) < 1e-7
    end

    @testset "3D" begin
        a = Point3D(1.0, 1.0, 1.0)
        b = Point3D(1.0, 1.0, 1.5)
        c = Point3D(1.0, 1.5, 1.0)
        d = Point3D(1.5, 1.0, 1.0)
        t = Tetrahedron(a, b, c, d)
        @test abs(volume(t)-1/16) < 1e-7
        c = centroid(t)
        @test abs(getx(c)-1.125) < 1e-7
        @test abs(gety(c)-1.125) < 1e-7
        @test abs(getz(c)-1.125) < 1e-7
        c = circumcenter(t)
        @test abs(getx(c)-1.25) < 1e-7
        @test abs(gety(c)-1.25) < 1e-7
        @test abs(getz(c)-1.25) < 1e-7
        @test abs(circumradius2(t)-0.25^2*3) < 1e-7
    end
end

@testset "2D Polygons" begin
    ll = Point(1.0,1.0)
    lr = Point(1.2,1.0)
    ur = Point(1.2,1.2)
    ul = Point(1.0,1.2)
    poly = Polygon(ll, lr, ur, ul)
    @test isempty( setdiff(getpoints(poly), [ll,lr,ur,ul]) )
    @test isempty( setdiff(getlines(poly), [Line(ll,lr), Line(lr,ur), Line(ur,ul), Line(ul,ll)]) )
    @test inpolygon(poly, Point(1.1,1.1)) == true
    @test inpolygon(poly, Point(1.1,1.3)) == false
end
# that's it for today


text:
using Plots,DifferentialEquations,LaTeXStrings

##
V(x,y) = (1/2)*(x^2 + y^2 + 2*x^2*y - (2/3)*y^3)
x1 = y1 = range(-1,stop = 1, length = 1000)

##
V(x,y) = (1/2)*(x^2 + y^2 + 2*x^2*y - (2/3)*y^3)
x = y = range(-1,stop = 1, length = 1000)
contourf(x,y,V)
contour!(x,y,V,levels=0:1/6:1,title = "Henon-Heiles potential: Intensity map",linecolor = :white)
# savefig("potential_together.png")

##
function Hénon_Heiles(du,u,p,t)
    x  = u[1]
    y  = u[2]
    dx = u[3]
    dy = u[4]
    du[1] = dx
    du[2] = dy
    du[3] = -x - 2x*y
    du[4] = y^2 - y -x^2
end

plot()
e = 0.0833
for y in -0.2:0.1:0.4
    for ydot in -0.1:0.1:0.2
        # y = 0.02
        # ydot = -0.08
        xdot = sqrt(2*e - ydot^2 - y^2 + (2//3)*y^3)
        initial = [0.0,y,xdot,ydot]
        tspan = (0,500.)
        prob = ODEProblem(Hénon_Heiles, initial, tspan)
        sol = solve(prob, saveat = 0.001,Rodas5());
        Y = [i[2] for i in sol.u if (abs(i[1]) < 9.e-4)];
        Ydot = [i[4] for i in sol.u if (abs(i[1]) < 9.e-4)];
        scatter!(Y,Ydot,xlabel = "y",ylabel = "ẏ",label = false,title = "Poincare section: E = $(e)", titlefont=10,markersize = 2,markercolor = :black)
    end
end
plot!()
##
function Hénon_Heiles(du,u,p,t)
    x  = u[1]
    y  = u[2]
    dx = u[3]
    dy = u[4]
    du[1] = dx
    du[2] = dy
    du[3] = -x - 2x*y
    du[4] = y^2 - y -x^2
end
e = 0.125
y = 0.02
ydot = -0.08
xdot = sqrt(2*e - ydot^2 - y^2 + (2//3)*y^3)
initial = [0.0,y,xdot,ydot]
tspan = (0,500.)
prob = ODEProblem(Hénon_Heiles, initial, tspan)
sol = solve(prob, saveat = 0.001);
plot(sol,vars = (1,2,4))

##

max_stars_repo_name: Ved-Mahajan/Henon-Heiles-System

max_issues_repo_name: Ved-Mahajan/Henon-Heiles-System

max_forks_repo_name: Ved-Mahajan/Henon-Heiles-System


text:
"""

Basically:

if I ∈ [3,5,7,9]
    c[(I - 1) >> 1
else
    b * I
end

because

c = b .* [3, 5, 7, 9]

"""
@generated function vmulnoth(::Static{I}, b, c::NTuple{N,Int}) where {I,N}
    Is = (I - 1) >> 1
    if isodd(I) && 1 ≤ Is ≤ N
        Expr(:block, Expr(:meta,:inline), Expr(:call, :getindex, :c, Is))
    elseif I ∈ (6, 10)
        Expr(:block, Expr(:meta,:inline), Expr(:call, :vmul, Expr(:call, Expr(:curly, :Static, 2)), Expr(:call, :getindex, :c, I >> 2)))
    else
        Expr(:block, Expr(:meta,:inline), Expr(:call, :vmul, Expr(:call, Expr(:curly, :Static, I)), :b))
    end
end
@inline vmulnoth(a, b, c) = vmul(a, b)
@inline vmulnoth(a::Static, b, ::Nothing) = vmul(a, b)
@inline tdot(a::Tuple{Static{I1},Vararg}, b::Tuple{I2,Vararg}, c::Tuple{I3,Vararg}) where {I1,I2,I3} = vmulnoth(Static{I1}(),b[1],c[1])
@inline tdot(a::Tuple{Static{I1},Vararg}, b::Tuple{I2,Vararg}, c::Tuple{Nothing,Vararg}) where {I1,I2} = vmul(Static{I1}(),b[1])
@inline tdot(a::Tuple{I1,Vararg}, b::Tuple{I2,Vararg}, c::Tuple{I3,Vararg}) where {I1,I2,I3} = vmul(a[1],b[1])

@inline tdot(a::Tuple{Static{I1},I4,Vararg}, b::Tuple{I2,I5,Vararg}, c::Tuple{I3,I6,Vararg}) where {I1,I2,I3,I4,I5,I6} = vadd(vmulnoth(Static{I1}(),b[1],c[1]), tdot(Base.tail(a), Base.tail(b), Base.tail(c)))
@inline tdot(a::Tuple{Static{I1},I4,Vararg}, b::Tuple{I2,I5,Vararg}, c::Tuple{Nothing,I6,Vararg}) where {I1,I2,I4,I5,I6} = vadd(vmul(Static{I1}(),b[1]), tdot(Base.tail(a), Base.tail(b), Base.tail(c)))
@inline tdot(a::Tuple{I1,I4,Vararg}, b::Tuple{I2,I5,Vararg}, c::Tuple{I3,I6,Vararg}) where {I1,I2,I3,I4,I5,I6} = vadd(vmul(a[1],b[1]), tdot(Base.tail(a), Base.tail(b), Base.tail(c)))

struct OffsetPrecalc{T,P<:AbstractStridedPointer{T},I} <: AbstractStridedPointer{T}
    ptr::P
    precalc::I
end
@inline Base.pointer(ptr::OffsetPrecalc) = pointer(ptr.ptr)
@inline Base.similar(ptr::OffsetPrecalc, p::Ptr) = OffsetPrecalc(similar(ptr.ptr, p), ptr.precalc)
@inline pointerforcomparison(p::OffsetPrecalc) = pointerforcomparison(p.ptr)
@inline pointerforcomparison(p::OffsetPrecalc, i::Tuple) = pointerforcomparison(p.ptr, i)
@inline offsetprecalc(x, ::Any) = x
@inline offsetprecalc(x::PackedStridedBitPointer, ::Val{<:Any}) = x
@inline offsetprecalc(x::RowMajorStridedBitPointer, ::Val{<:Any}) = x
# descript is a tuple of (unrollfactor) for ech ind; if it shouldn't preallocate, unrollfactor may be set to 1
function precalc_quote_from_descript(descript)
    precalc = Expr(:tuple)
    anyprecalcs = false
    pstrides = Expr(:block, Expr(:(=), :pstride, Expr(:(.), :p, QuoteNode(:strides))))
    for (i,uf) ∈ enumerate(descript)
        if uf < 3
            push!(precalc.args, nothing)
        else
            t = Expr(:tuple)
            anyprecalcs = true
            pstride_i = Symbol(:pstride_, i)
            push!(pstrides.args, Expr(:(=), pstride_i, Expr(:ref, :pstride, i)))
            for u ∈ 3:uf
                if isodd(u)
                    push!(t.args, Expr(:call, :vmul, u, pstride_i))
                end
            end
            push!(precalc.args, t)
        end
    end
    if anyprecalcs
        Expr(:block, Expr(:meta,:inline), pstrides, Expr(:call, :OffsetPrecalc, :p, precalc))
    else
        Expr(:block, Expr(:meta,:inline), :p)
    end
end
@inline offsetprecalc(p::AbstractColumnMajorStridedPointer{T,0}, ::Val{<:Any}) where {T} = p
@generated function offsetprecalc(p::AbstractColumnMajorStridedPointer, ::Val{descript}) where {descript}
    precalc_quote_from_descript(Base.tail(descript))
end

@inline stride1offset(ptr::OffsetPrecalc, i) = stride1offset(ptr.ptr, i)

@inline stridedoffset(ptr::OffsetPrecalc{T,<:AbstractColumnMajorStridedPointer}, ::Tuple{}) where {T} = Zero()
@inline stridedoffset(ptr::OffsetPrecalc{T,<:AbstractColumnMajorStridedPointer{T}}, i::Tuple{I}) where {I,T} = Zero()
# @inline stridedoffset(ptr::OffsetPrecalc{T,<:AbstractColumnMajorStridedPointer{T}}, i::Integer) where {T} = vmulnp(static_sizeof(T), i)
@inline function stridedoffset(ptr::OffsetPrecalc{T,<:AbstractColumnMajorStridedPointer{T}}, i::Tuple{I1,I2,Vararg}) where {I1,I2,T}
    @inbounds tdot(Base.tail(i), ptr.ptr.strides, ptr.precalc)
end

@inline offsetprecalc(p::AbstractRowMajorStridedPointer{T,0}, ::Val{<:Any}) where {T} = p
@generated function offsetprecalc(p::AbstractRowMajorStridedPointer, ::Val{descript}) where {descript}
    precalc_quote_from_descript(Base.tail(reverse(descript)))
end

@inline function stridedoffset(ptr::OffsetPrecalc{T,<:AbstractRowMajorStridedPointer{T,1}}, i::Tuple{I1,I2}) where {I1,I2,T}
    @inbounds vmulnoth(i[1], ptr.ptr.strides[1], ptr.precalc[1])
end
@inline function stridedoffset(ptr::OffsetPrecalc{T,<:AbstractRowMajorStridedPointer{T}}, i::Tuple{I1,I2,I3,Vararg}) where {I1,I2,I3,T}
    ri = reverse(i)
    @inbounds tdot(Base.tail(ri), ptr.ptr.strides, ptr.precalc)
end

@generated function offsetprecalc(p::AbstractSparseStridedPointer, ::Val{descript}) where {descript}
    precalc_quote_from_descript(descript)
end
@inline stridedoffset(ptr::OffsetPrecalc{T,<:AbstractSparseStridedPointer{T}}, i::Tuple{I}) where {I,T} = @inbounds vmulnoth(i[1], ptr.ptr.strides[1], ptr.precalc[1])
# @inline stridedoffset(ptr::OffsetPrecalc{T,<:AbstractSparseStridedPointer{T}}, i::Integer) where {T} = @inbounds vmulnoth(i, ptr.ptr.strides[1], ptr.precalc[1])
@inline function stridedoffset(ptr::OffsetPrecalc{T,<:AbstractSparseStridedPointer{T}}, i::Tuple{I1,I2,Vararg}) where {I1,I2,T}
    @inbounds tdot(i, ptr.ptr.strides, ptr.precalc)
end


text:
#=
Author:= Shreeda Bhat
License:= MIT
UTF-8 =#

#=
While in Julia holds this syntax

while *condition*
    *loop body*
end

=#

println("*********************** 1st Example *************************")

n = 0
c = 0
while n<10
    global n+=1
    println("number is now $n")
    global c += 1
end
println("""count = $c""")

println("***********************  2nd Example*************************")

myfriends = ["a","b","c","d"]

i = 1
while i<=length(myfriends)
    friend = myfriends[i]
    println("Hi $friend")
    global i+=1
end