text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

SUITE = BenchmarkGroup()
model = "ISS"
cases = ["ISSF01-ISS01", "ISSF01-ISU01",
         "ISSF01-ISS01-discrete", "ISSF01-ISU01-discrete",
         "ISSC01-ISS02", "ISSC01-ISU02",
         "ISSC01-ISS02-discrete", "ISSC01-ISU02-discrete"]

SUITE[model] = BenchmarkGroup()

include("ISS.jl")
validation = []

# ----------------------------------------
#  ISSF01 (dense time)
# ----------------------------------------
prob_ISSF01 = ISSF01()

dirs = CustomDirections([C3, -C3])

# ISSF01-ISS01
sol_ISSF01 = solve(prob_ISSF01, T=20.0, alg=LGG09(δ=6e-4, template=dirs, sparse=true, cache=false))
property = (-ρ(-C3, sol_ISSF01) >= -0.0007) && (ρ(C3, sol_ISSF01) <= 0.0007)
push!(validation, Int(property))
SUITE[model][cases[1]] = @benchmarkable solve($prob_ISSF01, T=20.0, alg=LGG09(δ=6e-4, template=$dirs, sparse=true, cache=false))

# ISSF01-ISU01
sol_ISSF01 = solve(prob_ISSF01, T=20.0, alg=LGG09(δ=6e-4, template=dirs, sparse=true, cache=false))
property = (-ρ(-C3, sol_ISSF01) >= -0.0005) && (ρ(C3, sol_ISSF01) <= 0.0005)
push!(validation, Int(property))
SUITE[model][cases[2]] = @benchmarkable solve($prob_ISSF01, T=20.0, alg=LGG09(δ=6e-4, template=$dirs, sparse=true, cache=false))

# ----------------------------------------
#  ISSF01 (discrete time)
# ----------------------------------------

dirs = CustomDirections([C3, -C3])

# ISSF01-ISS01
sol_ISSF01 = solve(prob_ISSF01, T=20.0, alg=LGG09(δ=0.01, template=dirs, sparse=true, cache=false, approx_model=NoBloating()))
property = (-ρ(-C3, sol_ISSF01) >= -0.0007) && (ρ(C3, sol_ISSF01) <= 0.0007)
push!(validation, Int(property))
SUITE[model][cases[3]] = @benchmarkable solve($prob_ISSF01, T=20.0, alg=LGG09(δ=6e-4, template=$dirs, sparse=true, cache=false))

# ISSF01-ISU01
sol_ISSF01 = solve(prob_ISSF01, T=20.0, alg=LGG09(δ=0.01, template=dirs, sparse=true, cache=false, approx_model=NoBloating()))
property = (-ρ(-C3, sol_ISSF01) >= -0.0005) && (ρ(C3, sol_ISSF01) <= 0.0005)
push!(validation, Int(property))
SUITE[model][cases[4]] = @benchmarkable solve($prob_ISSF01, T=20.0, alg=LGG09(δ=6e-4, template=$dirs, sparse=true, cache=false))

# ----------------------------------------
#  ISSC01 (dense time)
# ----------------------------------------
prob_ISSC01 = ISSC01()

dirs = CustomDirections([C3_ext, -C3_ext])

# ISSC01-ISS02
sol_ISSC01 = solve(prob_ISSC01, T=20.0, alg=LGG09(δ=0.01, template=dirs, sparse=true, cache=true))
property = (-ρ(-C3_ext, sol_ISSC01) >= -0.0005) && (ρ(C3_ext, sol_ISSC01) <= 0.0005)
push!(validation, Int(property))
SUITE[model][cases[5]] = @benchmarkable solve($prob_ISSC01, T=20.0, alg=LGG09(δ=0.01, template=$dirs, sparse=true, cache=true))

# ISSC01-ISU02
sol_ISSC01 = solve(prob_ISSC01, T=20.0, alg=LGG09(δ=0.01, template=dirs, sparse=true, cache=true))
property = (-ρ(-C3_ext, sol_ISSC01) >= -0.00017) && (ρ(C3_ext, sol_ISSC01) <= 0.00017)
push!(validation, Int(property))
SUITE[model][cases[6]] = @benchmarkable solve($prob_ISSC01, T=20.0, alg=LGG09(δ=0.01, template=$dirs, sparse=true, cache=true))

# ----------------------------------------
#  ISSC01 (discrete time)
# ----------------------------------------

prob_ISSC01 = ISSC01()

dirs = CustomDirections([C3_ext, -C3_ext])

# ISSC01-ISS02
sol_ISSC01 = solve(prob_ISSC01, T=20.0, alg=LGG09(δ=0.01, template=dirs, sparse=true, cache=true, approx_model=NoBloating()))
property = (-ρ(-C3_ext, sol_ISSC01) >= -0.0005) && (ρ(C3_ext, sol_ISSC01) <= 0.0005)
push!(validation, Int(property))
SUITE[model][cases[7]] = @benchmarkable solve($prob_ISSC01, T=20.0, alg=LGG09(δ=0.01, template=$dirs, sparse=true, cache=true))

# ISSC01-ISU02
sol_ISSC01 = solve(prob_ISSC01, T=20.0, alg=LGG09(δ=0.01, template=dirs, sparse=true, cache=true, approx_model=NoBloating()))
property = (-ρ(-C3_ext, sol_ISSC01) >= -0.00017) && (ρ(C3_ext, sol_ISSC01) <= 0.00017)
push!(validation, Int(property))
SUITE[model][cases[8]] = @benchmarkable solve($prob_ISSC01, T=20.0, alg=LGG09(δ=0.01, template=$dirs, sparse=true, cache=true))

sol_ISSF01 = nothing
sol_ISSC01 = nothing
GC.gc()

# ==============================================================================
# Execute benchmarks and save benchmark results
# ==============================================================================

# tune parameters
tune!(SUITE)

# run the benchmarks
results = run(SUITE, verbose=true)

# return the sample with the smallest time value in each test
println("minimum time for each benchmark:\n", minimum(results))

# return the median for each test
println("median time for each benchmark:\n", median(results))

# export runtimes
runtimes = Dict()
for (i, c) in enumerate(cases)
    t = median(results[model][c]).time * 1e-9
    runtimes[c] = round(t, digits=4)
end

for (i, c) in enumerate(cases)
    print(io, "JuliaReach, $model, $c, $(validation[i]), $(runtimes[c])\n")
end

# ==============================================================================
# Plot
# ==============================================================================

dirs = CustomDirections([C3, -C3])
prob_ISSF01 = ISSF01()
sol_ISSF01 = solve(prob_ISSF01, T=20.0, alg=LGG09(δ=6e-4, template=dirs, sparse=true, cache=false))
πsol_ISSF01 = flatten(sol_ISSF01)

fig = Plots.plot()
Plots.plot!(fig, πsol_ISSF01, vars=(0, 1),
            linecolor=:blue, color=:blue, alpha=0.8,
            tickfont=font(30, "Times"), guidefontsize=45,
            xlab=L"t",
            ylab=L"y_{3}",
            xtick=[0, 5, 10, 15, 20.], ytick=[-0.00075, -0.0005, -0.00025, 0, 0.00025, 0.0005, 0.00075],
            xlims=(0., 20.), ylims=(-0.00075, 0.00075),
            yformatter = y -> string(floor(Int, y/1e-4)), annotate=(0.5, 7.5e-4, text("⋅10⁻⁴", :left, 30)),
            bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
            size=(1000, 1000))

savefig(fig, joinpath(@__DIR__, "ARCH-COMP21-JuliaReach-ISS-ISSF01"))

# -------------

dirs = CustomDirections([C3_ext, -C3_ext])
prob_ISSC01 = ISSC01()
sol_ISSC01 = solve(prob_ISSC01, T=20.0, alg=LGG09(δ=0.01, template=dirs, sparse=true, cache=false))
πsol_ISSC01 = flatten(sol_ISSC01)

fig = Plots.plot()
Plots.plot!(fig, πsol_ISSC01, vars=(0, 1),
            linecolor=:blue, color=:blue, alpha=0.8,
            tickfont=font(30, "Times"), guidefontsize=45,
            xlab=L"t",
            ylab=L"y_{3}",
            xtick=[0, 5, 10, 15, 20.], ytick=[-0.0002, -0.0001, 0.0, 0.0001, 0.0002],
            xlims=(0., 20.), ylims=(-0.0002, 0.0002),
            yformatter = y -> string(floor(Int, y/1e-4)), annotate=(0.5, 7.5e-4, text("⋅10⁻⁴", :left, 30)),
            bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
            size=(1000, 1000))

savefig(fig, joinpath(@__DIR__, "ARCH-COMP21-JuliaReach-ISS-ISSC01"))


text:
using StanSample

ProjDir = @__DIR__
cd(ProjDir)

shared_file = "shared_funcs.stan"

bernoulli_model = "
  functions{

    #include $(shared_file)
    //#include shared_funcs.stan // a comment
    //#include shared_funcs.stan // a comment
    //#include /Users/rob/shared_funcs.stan // a comment

    void model_specific_function(){
        real x = 1.0;
        return;
    }

  }
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    model_specific_function();
    theta ~ beta(my_function(),1);
    y ~ bernoulli(theta);
  }
";

tmpdir = joinpath(@__DIR__, "tmp")

stanmodel = SampleModel("bernoulli", bernoulli_model, tmpdir=tmpdir)

observeddata = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

(sample_file, log_file) = stan_sample(stanmodel, data=observeddata)

if sample_file !== Nothing
  # Convert to an MCMCChains.Chains object
  chns = read_samples(stanmodel)

  # Describe the MCMCChains using MCMCChains statistics
  cdf = describe(chns)
  display(cdf)

  # Show cmdstan summary in DataFrame format
  sdf = read_summary(stanmodel)
  display(sdf)
  println()

  # Retrieve mean value of theta from the summary
  sdf[:theta, :mean]
end

ext: jl


text:
module StippleUI

import Genie
import Stipple

using Stipple.Reexport

const DEFAULT_WRAPPER = Genie.Renderer.Html.template
const NO_WRAPPER = f->f()

#===#

function theme() :: String
  Genie.Router.route("/css/stipple/quasar.min.css") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "css", "quasar.min.css"), String),
      :css) |> Genie.Renderer.respond
  end

  Stipple.Elements.stylesheet("$(Genie.config.base_path)css/stipple/quasar.min.css")
end

#===#

function deps() :: String
  Genie.Router.route("/js/stipple/quasar.umd.min.js") do
    Genie.Renderer.WebRenderable(
      read(joinpath(@__DIR__, "..", "files", "js", "quasar.umd.min.js"), String),
      :javascript) |> Genie.Renderer.respond
  end

  Genie.Renderer.Html.script(src="$(Genie.config.base_path)js/stipple/quasar.umd.min.js")
end

#===#

include("API.jl")
include("Badge.jl")
include("Banner.jl")
include("BigNumber.jl")
include("Button.jl")
include("Card.jl")
include("Checkbox.jl")
include("Chip.jl")
include("Dashboard.jl")
include("Dialog.jl")
include("Drawer.jl")
include("Form.jl")
include("FormInput.jl")
include("Heading.jl")
include("Icon.jl")
include("Layout.jl")
include("List.jl")
include("Radio.jl")
include("Range.jl")
include("Select.jl")
include("Separator.jl")
include("Space.jl")
include("Table.jl")
include("Toggle.jl")
include("Uploader.jl")

#===#

import .API: quasar, quasar_pure, vue, vue_pure

export quasar, quasar_pure, vue, vue_pure, @click

@reexport using .Badge
@reexport using .Banner
@reexport using .BigNumber
@reexport using .Button
@reexport using .Card
@reexport using .Checkbox
@reexport using .Chip
@reexport using .Dashboard
@reexport using .Dialog
@reexport using .Form
@reexport using .FormInput
@reexport using .Heading
@reexport using .Icon
@reexport using .List
@reexport using .Radio
@reexport using .Range
@reexport using .Select
@reexport using .Separator
@reexport using .Space
@reexport using .Table
@reexport using .Toggle
@reexport using .Uploader

#===#

"""
    `@click(expr)`

Defines a js routine that is called by a click of the quasar component.
If a symbol argument is supplied, `@click` sets this value to true.

`@click("savefile = true")` or `@click("myjs_func();")` or `@click(:button)`

Modifers can be appended:
```
@click(:me, :native)
# "@click.native='me = true'"
```
"""
macro click(expr, mode="")
  quote
    x = $(esc(expr))
    m = $(esc(mode))
    if x isa Symbol
      """@click$(m == "" ? "" : ".$m")='$x = true'"""
    else
      "@click='$(replace(x, "'" => raw"\'"))'"
    end
  end
end

#===#

function __init__()
  push!(Stipple.Layout.THEMES, theme)
  push!(Stipple.DEPS, deps)
end

end


text:
# Unit tests for api/md_result.jl

@testset "api/md_result.jl" begin
    @testset "Unimplemented Functions" begin
        struct MockResult1 <: MolecularDynamicsResult end

        target = MockResult1()

        @test_throws UnimplementedError{MockResult1} get_time_range(target)
        @test_throws UnimplementedError{MockResult1} get_num_bodies(target)
        @test_throws UnimplementedError{MockResult1} get_bounding_box(target)
        @test_throws UnimplementedError{MockResult1} get_boundary_conditions(target)

        @test_throws UnimplementedError{MockResult1} get_positions(target)
        @test_throws UnimplementedError{MockResult1} get_positions(target, 1)
        @test_throws UnimplementedError{MockResult1} get_velocities(target)
        @test_throws UnimplementedError{MockResult1} get_velocities(target, 1)
        @test_throws UnimplementedError{MockResult1} get_particles(target)
        @test_throws UnimplementedError{MockResult1} get_particles(target, 1)

        @test_throws UnimplementedError{MockResult1} Atomistic.temperature(target)
        @test_throws UnimplementedError{MockResult1} Atomistic.temperature(target, 1)
        @test_throws UnimplementedError{MockResult1} Atomistic.kinetic_energy(target)
        @test_throws UnimplementedError{MockResult1} Atomistic.kinetic_energy(target, 1)
        @test_throws UnimplementedError{MockResult1} Atomistic.potential_energy(target)
        @test_throws UnimplementedError{MockResult1} Atomistic.potential_energy(target, 1)
    end
    @testset "Default Implementations" begin
        struct MockResult2 <: MolecularDynamicsResult end
        Atomistic.get_time_range(::MockResult2) = collect(0:2:100)u"ns"
        Atomistic.get_num_bodies(::MockResult2) = 3
        Atomistic.get_bounding_box(::MockResult2) = (@SVector [(@SVector [100, 0, 0]), (@SVector [0, 100, 0]), (@SVector [0, 0, 100])])u"bohr"
        Atomistic.get_boundary_conditions(::MockResult2) = @SVector [Periodic(), Periodic(), Periodic()]
        Atomistic.get_positions(::MockResult2, t::Integer) = [
            (@SVector [0, 0, 0])u"bohr",
            (@SVector [-t, 0, t])u"bohr",
            (@SVector [t, 0, -t])u"bohr"
        ]
        Atomistic.get_particles(::MockResult2, t::Integer) = [
            AtomsBase.Atom(:Ar, (@SVector [0, 0, 0])u"bohr", (@SVector [0, 0, 0])u"bohr/ns"),
            AtomsBase.Atom(:Ar, (@SVector [-t, 0, t])u"bohr", (@SVector [-2t, 0, 2t])u"bohr/ns"),
            AtomsBase.Atom(:Ar, (@SVector [t, 0, -t])u"bohr", (@SVector [2t, 0, -2t])u"bohr/ns")
        ]
        Atomistic.kinetic_energy(::MockResult2, t::Integer) = t * u"bohr"
        Atomistic.potential_energy(::MockResult2, t::Integer) = -2t * u"bohr"

        target = MockResult2()

        @test ismissing(reference_temperature(target))

        @test get_time(target) == 100u"ns"
        @test get_time(target, 1) == 0u"ns"

        system = get_system(target)
        @test bounding_box(system) == (@SVector [(@SVector [100, 0, 0]), (@SVector [0, 100, 0]), (@SVector [0, 0, 100])])u"bohr"
        @test boundary_conditions(system) == @SVector [Periodic(), Periodic(), Periodic()]
        @test collect(system) == [
            AtomsBase.Atom(:Ar, (@SVector [0, 0, 0])u"bohr", (@SVector [0, 0, 0])u"bohr/ns"),
            AtomsBase.Atom(:Ar, (@SVector [-51, 0, 51])u"bohr", (@SVector [-102, 0, 102])u"bohr/ns"),
            AtomsBase.Atom(:Ar, (@SVector [51, 0, -51])u"bohr", (@SVector [102, 0, -102])u"bohr/ns"),
        ]

        system1 = get_system(target, 1)
        @test bounding_box(system1) == (@SVector [(@SVector [100, 0, 0]), (@SVector [0, 100, 0]), (@SVector [0, 0, 100])])u"bohr"
        @test boundary_conditions(system1) == @SVector [Periodic(), Periodic(), Periodic()]
        @test collect(system1) == [
            AtomsBase.Atom(:Ar, (@SVector [0, 0, 0])u"bohr", (@SVector [0, 0, 0])u"bohr/ns"),
            AtomsBase.Atom(:Ar, (@SVector [-1, 0, 1])u"bohr", (@SVector [-2, 0, 2])u"bohr/ns"),
            AtomsBase.Atom(:Ar, (@SVector [1, 0, -1])u"bohr", (@SVector [2, 0, -2])u"bohr/ns"),
        ]

        @test Atomistic.total_energy(target) == -51u"bohr"
        @test Atomistic.total_energy(target, 1) == -1u"bohr"
        rdf1 = Atomistic.rdf(target)
        rdf2 = Atomistic.rdf(target, 1)
        rdf3 = Atomistic.rdf(target, 1, 51)
        @test rdf1 isa NamedTuple{(:r, :g),<:Tuple{<:AbstractVector{<:Unitful.Length},<:AbstractVector{Float64}}}
        @test rdf1 == rdf2
        @test rdf1 == rdf3
    end
end


text:
abstract type AbstractTechnology{G} end

struct ThermalGenerators{G1} <: AbstractTechnology{G1}

    name::Vector{String}             # g
    owner::Vector{String}            # g

    mingen::Vector{Float64}          # (MW/unit, g)
    maxgen::Vector{Float64}          # (MW/unit, g)

    minuptime::Vector{Int}           # (hours, g)
    mindowntime::Vector{Int}         # (hours, g)
    maxrampup::Vector{Float64}       # (MW/unit/hour, g)
    maxrampdown::Vector{Float64}     # (MW/unit/hour, g)

    capacitycredit::Vector{Float64}  # (fraction, g)

    function ThermalGenerators{}(args...)
        G = length(first(args))
        @assert all(a -> length(a) == G, args)
        new{G}(args...)
    end

end

struct VariableGenerators{G2} <: AbstractTechnology{G2}

    name::Vector{String}    # g
    owner::Vector{String}   # g

    maxgen::Vector{Float64} # (MW/unit, g)

    capacitycredit::Vector{Float64} # (fraction, g)

    function VariableGenerators{}(args...)
        G = length(first(args))
        @assert all(a -> length(a) == G, args)
        new{G}(args...)
    end

end

struct StorageDevices{G3} <: AbstractTechnology{G3}

    name::Vector{String}  # g
    owner::Vector{String} # g

    maxcharge::Vector{Float64}    # (MW, g)
    maxdischarge::Vector{Float64} # (MW, g)
    maxenergy::Vector{Float64}    # (MWh, g)

    chargeefficiency::Vector{Float64}    # (fraction, g)
    dischargeefficiency::Vector{Float64} # (fraction, g)
    carryoverefficiency::Vector{Float64}  # (fraction, g)

    capacitycredit::Vector{Float64}  # (fraction, g)

    function StorageDevices{}(args...)
        G = length(first(args))
        @assert all(a -> length(a) == G, args)
        new{G}(args...)
    end

end

struct Interfaces{I,R} <: AbstractTechnology{I}

    name::Vector{String}
    regions::Vector{Tuple{Int,Int}}

    function Interfaces{R}(
        name::Vector{String},
        regions::Vector{Tuple{Int,Int}}
    ) where R

        I = length(name)
        @assert I == length(regions)

        @assert all(x -> 1 <= first(x) <= R, regions)
        @assert all(x -> 1 <= last(x) <= R, regions)
        @assert allunique(tuple.(minimum.(regions), maximum.(regions)))

        new{I,R}(name, regions)

    end

end

struct Technologies{G1,G2,G3,I,R}

    thermal::ThermalGenerators{G1}
    variable::VariableGenerators{G2}
    storage::StorageDevices{G3}
    interface::Interfaces{I,R}

    function Technologies{}(
        thermal::ThermalGenerators{G1},
        variable::VariableGenerators{G2},
        storage::StorageDevices{G3},
        interface::Interfaces{I,R}
    ) where {G1,G2,G3,I,R}

        allunique(vcat(thermal.name, variable.name,
                       storage.name, interface.name)) ||
            error("Resource class names must be unique")

        new{G1,G2,G3,I,R}(thermal, variable, storage, interface)

    end

end


text:
using Distributions
using DynamicHMC
using LogDensityProblems
using TransformVariables
using Random

struct Schaefer{F}
    r::F
    K::F
end
(s::Schaefer)(p) = p + s.r * p * (1 - p / s.K)

struct SchaeferProblem{T}
    p0::Vector{T}
    p1::Vector{T}
    σ::T
end
function (sp::SchaeferProblem)(θ)
    ll = zero(θ.r)
    ll += logpdf(Normal(100.0, 50.0), θ.K)
    ll += logpdf(LogNormal(log(0.1) - 0.1^2 / 2, 0.1), θ.r)
    schaefer = Schaefer(θ.r, θ.K)
    errdist = LogNormal(-sp.σ^2 / 2, sp.σ)
    for idx in eachindex(sp.p0)
        res = schaefer(sp.p0[idx]) / sp.p1[idx]
        ll += logpdf(errdist, res)
    end
    ll
end

schaefer = Schaefer(0.1, 100.0)
σ = 0.1
errdist = LogNormal(σ^2 / 2, σ)
p0 = 100rand(30)
p1 = zeros(100)
for idx in eachindex(p0)
    p1[idx] = schaefer(p0[idx])
end

sprob = SchaeferProblem(p0, p1, σ)

sprob((r = 0.1, K = 100.0))
t = as((r = asℝ₊, K = asℝ₊))
ℓ = TransformedLogDensity(t, sprob)

∇ℓ = ADgradient(:ForwardDiff, ℓ)

## Fails for non-finiteness somewhere?
post = mcmc_with_warmup(Random.GLOBAL_RNG, ∇ℓ, 100_000;
                        initialization = (q = [0.1, 100.0],),
                        reporter = NoProgressReport())

histogram2d(first.(post.chain), last.(post.chain), bins = 100)

using Optim

optimize(θ -> -sprob((r = θ[1], K = θ[2])), [0.1, 100.0])

struct Rosenbrock{T}
    a::T
    b::T
end
(rosen::Rosenbrock)(θ) = -(rosen.a - θ.x)^2 - rosen.b * (θ.y - θ.x^2)^2

rosen = Rosenbrock(1.0, 10.0)
rosen_t = as((x = asℝ, y = asℝ))
rosen_ll = TransformedLogDensity(rosen_t, rosen)
∇rosen_ll = ADgradient(:ForwardDiff, rosen_ll)
rosen_post = mcmc_with_warmup(Random.GLOBAL_RNG, ∇rosen_ll, 10_000;
                              reporter = NoProgressReport())

histogram2d(first.(rosen_post.chain), last.(rosen_post.chain), bins = 100)

RNG = Random.GLOBAL_RNG

results = DynamicHMC.mcmc_keep_warmup(RNG, ∇rosen_ll, 0; reporter = NoProgressReport())
steps = DynamicHMC.mcmc_steps(results.sampling_logdensity, results.final_warmup_state)
Q = results.final_warmup_state.Q

# a single update step
Q, tree_stats = DynamicHMC.mcmc_next_step(steps, Q)

grid = range(-1.5, 3.5; length = 101)

# extract the position
Q.q

DynamicHMC.Diagnostics.leapfrog_trajectory(∇rosen_ll, Q.q, -5:5)


text:
using Documenter
using SQLDataFrameTools
using Dates

using Documenter, SQLDataFrameTools

makedocs(
    modules = [SQLDataFrameTools],
    sitename="SQLDataFrameTools.jl",
    authors = "Matt Lawless",
    format = Documenter.HTML(),
)

deploydocs(
    repo = "github.com/lawless-m/SQLDataFrameTools.jl.git",
    devbranch = "main",
    push_preview = true,
)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "GRASS"
version = v"7.8.5"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://grass.osgeo.org/grass78/source/grass-$version.tar.gz", "a359bb665524ecccb643335d70f5436b1c84ffb6a0e428b78dffebacd983ff37"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""

cd $WORKSPACE/srcdir/grass-*/

export LDFLAGS="-L${libdir}"

#patch remove g.proj from general/Makefile for right now, throws a segmentation fault
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done

if [[ "${target}" == "${MACHTYPE}" ]]; then
    # Remove system libexpat to avoid confusion
    rm /usr/lib/libexpat.so*
fi

# Force linking to libiconv with `-liconv`
sed -i 's/ICONVLIB=.*/ICONVLIB=-liconv/' configure

# Need a posix regex for Windows
if [[ "${target}" == *-mingw* ]]; then
    cp "${includedir}/pcreposix.h" "${includedir}/regex.h"
    # Force linking to libregex with `-lpcreposix-0`
    sed -i 's/-lregex/-lpcreposix-0/' configure
fi

./configure \
    --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --enable-shared \
    --with-python=no \
    --with-cxx \
    --with-lapack --with-lapack-includes=${includedir} \
    --with-tiff --with-tiff-includes=${includedir} \
    --with-png --with-png-includes=${includedir} \
    --with-sqlite --with-sqlite-includes=${includedir} \
    --with-opengl=no \
    --with-fftw --with-fftw-includes=${includedir} \
    --with-cairo --with-cairo-includes=${includedir}/cairo --with-cairo-ldflags=-lfontconfig \
    --with-freetype --with-freetype-includes=${includedir}/freetype2 \
    --with-regex \
    --with-zstd --with-zstd-includes=${includedir} \
    --with-geos=${bindir}/geos-config \
    --with-gdal \
    --with-proj --with-proj-share=${prefix}/share/proj --with-proj-includes=${includedir}

# Build only the libraries
DIRS=(
    include
    tools
    lib/external/shapelib
    lib/datetime
    lib/gis
    lib/linkm
    lib/db
    lib/btree2
    lib/vector
    db/drivers
    lib
    imagery/i.ortho.photo/lib
    vector/v.lrs/lib
    raster/r.li/r.li.daemon
    raster/r.sim/simlib
)

for dir in "${DIRS[@]}"; do
    make -j${nproc} -C "${dir}"
done

# Manually install libraries and header files
cp dist.*/lib/*.${dlext}* ${libdir}/.
cp -r dist.*/include/grass ${includedir}/.

if [[ "${target}" == *-mingw* ]]; then
    # Cover up the traces of the hack
    rm "${includedir}/regex.h"
fi

install_license COPYING GPL.TXT
"""

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

# The products that we will ensure are always built
products = [

    LibraryProduct(["libgrass_lidar", "libgrass_lidar.7.8"], :libgrass_lidar),
    LibraryProduct(["libgrass_neta", "libgrass_neta.7.8"], :libgrass_neta),
    LibraryProduct(["libgrass_htmldriver", "libgrass_htmldriver.7.8"], :libgrass_htmldriver),
    LibraryProduct(["libgrass_cluster", "libgrass_cluster.7.8"], :libgrass_cluster),
    LibraryProduct(["libgrass_gmath", "libgrass_gmath.7.8"], :libgrass_gmath),
    LibraryProduct(["libgrass_interpdata", "libgrass_interpdata.7.8"], :libgrass_interpdata),
    LibraryProduct(["libgrass_psdriver", "libgrass_psdriver.7.8"], :libgrass_psdriver),
    LibraryProduct(["libgrass_pngdriver", "libgrass_pngdriver.7.8"], :libgrass_pngdriver),
    LibraryProduct(["libgrass_dspf", "libgrass_dspf.7.8"], :libgrass_dspf),
    LibraryProduct(["libgrass_bitmap", "libgrass_bitmap.7.8"], :libgrass_bitmap),
    LibraryProduct(["libgrass_shape", "libgrass_shape.7.8"], :libgrass_shape),
    LibraryProduct(["libgrass_rtree", "libgrass_rtree.7.8"], :libgrass_rtree),
    LibraryProduct(["libgrass_cairodriver", "libgrass_cairodriver.7.8"], :libgrass_cairodriver),
    LibraryProduct(["libgrass_rowio", "libgrass_rowio.7.8"], :libgrass_rowio),
    LibraryProduct(["libgrass_lrs", "libgrass_lrs.7.8"], :libgrass_lrs),
    LibraryProduct(["libgrass_qtree", "libgrass_qtree.7.8"], :libgrass_qtree),
    LibraryProduct(["libgrass_gis", "libgrass_gis.7.8"], :libgrass_gis),
    LibraryProduct(["libgrass_gpde", "libgrass_gpde.7.8"], :libgrass_gpde),
    LibraryProduct(["libgrass_vector", "libgrass_vector.7.8"], :libgrass_vector),
    LibraryProduct(["libgrass_segment", "libgrass_segment.7.8"], :libgrass_segment),
    LibraryProduct(["libgrass_g3d", "libgrass_g3d.7.8"], :libgrass_g3d),
    LibraryProduct(["libgrass_display", "libgrass_display.7.8"], :libgrass_display),
    LibraryProduct(["libgrass_iortho", "libgrass_iortho.7.8"], :libgrass_iortho),
    LibraryProduct(["libgrass_btree2", "libgrass_btree2.7.8"], :libgrass_btree2),
    LibraryProduct(["libgrass_imagery", "libgrass_imagery.7.8"], :libgrass_imagery),
    LibraryProduct(["libgrass_datetime", "libgrass_datetime.7.8"], :libgrass_datetime),
    LibraryProduct(["libgrass_linkm", "libgrass_linkm.7.8"], :libgrass_linkm),
    LibraryProduct(["libgrass_ccmath", "libgrass_ccmath.7.8"], :libgrass_ccmath),
    LibraryProduct(["libgrass_driver", "libgrass_driver.7.8"], :libgrass_driver),
    LibraryProduct(["libgrass_gproj", "libgrass_gproj.7.8"], :libgrass_gproj),
    LibraryProduct(["libgrass_dgl", "libgrass_dgl.7.8"], :libgrass_dgl),
    LibraryProduct(["libgrass_dbmibase", "libgrass_dbmibase.7.8"], :libgrass_dbmibase),
    LibraryProduct(["libgrass_stats", "libgrass_stats.7.8"], :libgrass_stats),
    LibraryProduct(["libgrass_temporal", "libgrass_temporal.7.8"], :libgrass_temporal),
    LibraryProduct(["libgrass_sim", "libgrass_sim.7.8"], :libgrass_sim),
    LibraryProduct(["libgrass_btree", "libgrass_btree.7.8"], :libgrass_btree),
    LibraryProduct(["libgrass_dbstubs", "libgrass_dbstubs.7.8"], :libgrass_dbstubs),
    LibraryProduct(["libgrass_vedit", "libgrass_vedit.7.8"], :libgrass_vedit),
    LibraryProduct(["libgrass_calc", "libgrass_calc.7.8"], :libgrass_calc),
    LibraryProduct(["libgrass_raster", "libgrass_raster.7.8"], :libgrass_raster),
    LibraryProduct(["libgrass_cdhc", "libgrass_cdhc.7.8"], :libgrass_cdhc),
    LibraryProduct(["libgrass_interpfl", "libgrass_interpfl.7.8"], :libgrass_interpfl),
    LibraryProduct(["libgrass_dbmidriver", "libgrass_dbmidriver.7.8"], :libgrass_dbmidriver),
    LibraryProduct(["libgrass_rli", "libgrass_rli.7.8"], :libgrass_rli),
    LibraryProduct(["libgrass_arraystats", "libgrass_arraystats.7.8"], :libgrass_arraystats),
    LibraryProduct(["libgrass_sqlp", "libgrass_sqlp.7.8"], :libgrass_sqlp),
    LibraryProduct(["libgrass_dig2", "libgrass_dig2.7.8"], :libgrass_dig2),
    LibraryProduct(["libgrass_symb", "libgrass_symb.7.8"], :libgrass_symb),
    LibraryProduct(["libgrass_dbmiclient", "libgrass_dbmiclient.7.8"], :libgrass_dbmiclient),
    LibraryProduct(["libgrass_manage", "libgrass_manage.7.8"], :libgrass_manage),

    # ExecutableProduct("g.mkfontcap", :g_mkfontcap, "grass78/bin"),
    # ExecutableProduct("g.region", :g_region, "grass78/bin"),
    # ExecutableProduct("g.ppmtopng", :g_ppmtopng, "grass78/bin"),
    # ExecutableProduct("g.dirseps", :g_dirseps, "grass78/bin"),
    # ExecutableProduct("g.parser", :g_parser, "grass78/bin"),
    # ExecutableProduct("g.message", :g_message, "grass78/bin"),
    # ExecutableProduct("g.copy", :g_copy, "grass78/bin"),
    # ExecutableProduct("g.access", :g_access, "grass78/bin"),
    # ExecutableProduct("g.mapsets", :g_mapsets, "grass78/bin"),
    # ExecutableProduct("g.tempfile", :g_tempfile, "grass78/bin"),
    # ExecutableProduct("g.rename", :g_rename, "grass78/bin"),
    # ExecutableProduct("g.findfile", :g_findfile, "grass78/bin"),
    # ExecutableProduct("g.remove", :g_remove, "grass78/bin"),
    # ExecutableProduct("g.gui", :g_gui, "grass78/bin"),
    # ExecutableProduct("g.findetc", :g_findetc, "grass78/bin"),
    # ExecutableProduct("g.version", :g_version, "grass78/bin"),
    # ExecutableProduct("g.list", :g_list, "grass78/bin"),
    # ExecutableProduct("g.filename", :g_filename, "grass78/bin"),
    # ExecutableProduct("g.mapset", :g_mapset, "grass78/bin"),
    # ExecutableProduct("g.echo", :g_echo, "grass78/tools"),
    # ExecutableProduct("g.pnmcomp", :g_pnmcomp, "grass78/bin"),
    # ExecutableProduct("g.gisenv", :g_gisenv, "grass78/bin"),

    # ExecutableProduct("db.createdb", :db_createdb, "grass78/bin"),
    # ExecutableProduct("db.login", :db_login, "grass78/bin"),
    # ExecutableProduct("db.connect", :db_connect, "grass78/bin"),
    # ExecutableProduct("db.execute", :db_execute, "grass78/bin"),
    # ExecutableProduct("db.databases", :db_databases, "grass78/bin"),
    # ExecutableProduct("db.tables", :db_tables, "grass78/bin"),
    # ExecutableProduct("db.describe", :db_describe, "grass78/bin"),
    # ExecutableProduct("db.columns", :db_columns, "grass78/bin"),
    # ExecutableProduct("db.copy", :db_copy, "grass78/bin"),
    # ExecutableProduct("db.dropdb", :db_dropdb, "grass78/bin"),
    # ExecutableProduct("db.select", :db_select, "grass78/bin"),
    # ExecutableProduct("db.drivers", :db_drivers, "grass78/bin"),

    # ExecutableProduct("r.statistics", :r_statistics, "grass78/bin"),
    # ExecutableProduct("r.sim.water", :r_sim_water, "grass78/bin"),
    # ExecutableProduct("r.volume", :r_volume, "grass78/bin"),
    # ExecutableProduct("r.carve", :r_carve, "grass78/bin"),
    # ExecutableProduct("r.basins.fill", :r_basins_fill, "grass78/bin"),
    # ExecutableProduct("r.param.scale", :r_params_scale, "grass78/bin"),
    # ExecutableProduct("r.to.rast3elev", :r_to_rast3elev, "grass78/bin"),
    # ExecutableProduct("r.out.vrml", :r_out_vrml, "grass78/bin"),
    # ExecutableProduct("r.mapcalc", :r_mapcalc, "grass78/bin"),
    # ExecutableProduct("r.resample", :r_resample, "grass78/bin"),
    # ExecutableProduct("r.stream.extract", :r_stream_extract, "grass78/bin"),
    # ExecutableProduct("r.out.mat", :r_out_mat, "grass78/bin"),
    # ExecutableProduct("r.external", :r_external, "grass78/bin"),
    # ExecutableProduct("r.cost", :r_cost, "grass78/bin"),
    # ExecutableProduct("r.series", :r_series, "grass78/bin"),
    # ExecutableProduct("r.cross", :r_cross, "grass78/bin"),
    # ExecutableProduct("r.category", :r_category, "grass78/bin"),
    # ExecutableProduct("r.buffer", :r_buffer, "grass78/bin"),
    # ExecutableProduct("r.surf.contour", :r_surf_contour, "grass78/bin"),
    # ExecutableProduct("r.topmodel", :r_topmodel, "grass78/bin"),
    # ExecutableProduct("r.watershed", :r_watershed, "grass78/bin"),
    # ExecutableProduct("r.support", :r_support, "grass78/bin"),
    # ExecutableProduct("r.li.padsd", :r_li_padsd, "grass78/bin"),
    # ExecutableProduct("r.random.cells", :r_random_cells, "grass78/bin"),
    # ExecutableProduct("r.sunhours", :r_sunhours, "grass78/bin"),
    # ExecutableProduct("r.out.ppm3", :r_out_ppm3, "grass78/bin"),
    # ExecutableProduct("r.li.richness", :r_li_richness, "grass78/bin"),
    # ExecutableProduct("r.li.dominance", :r_li_dominance, "grass78/bin"),
    # ExecutableProduct("r.patch", :r_patch, "grass78/bin"),
    # ExecutableProduct("r.water.outlet", :r_water_outlet, "grass78/bin"),
    # ExecutableProduct("r.regression.multi", :r_regression_multi, "grass78/bin"),
    # ExecutableProduct("r.out.ascii", :r_out_ascii, "grass78/bin"),
    # ExecutableProduct("r.what.color", :r_what_color, "grass78/bin"),
    # ExecutableProduct("r.profile", :r_profile, "grass78/bin"),
    # ExecutableProduct("r.rescale.eq", :r_rescale_eq, "grass78/bin"),
    # ExecutableProduct("r.terraflow", :r_terraflow, "grass78/bin"),
    # ExecutableProduct("r.colors", :r_colors, "grass78/bin"),
    # ExecutableProduct("r.li.cwed", :r_li_cwed, "grass78/bin"),
    # ExecutableProduct("r.univar", :r_univar, "grass78/bin"),
    # ExecutableProduct("r.mode", :r_mode, "grass78/bin"),
    # ExecutableProduct("r.distance", :r_distance, "grass78/bin"),
    # ExecutableProduct("r.thin", :r_thin, "grass78/bin"),
    # ExecutableProduct("r.viewshed", :r_viewshed, "grass78/bin"),
    # ExecutableProduct("r.null", :r_null, "grass78/bin"),
    # ExecutableProduct("r.rescale", :r_rescale, "grass78/bin"),
    # ExecutableProduct("r.surf.area", :r_surf_area, "grass78/bin"),
    # ExecutableProduct("r.reclass", :r_reclass, "grass78/bin"),
    # ExecutableProduct("r.in.bin", :r_in_bin, "grass78/bin"),
    # ExecutableProduct("r.geomorphon", :r_geomorphon, "grass78/bin"),
    # ExecutableProduct("r.his", :r_his, "grass78/bin"),
    # ExecutableProduct("r.stats.quantile", :r_stats_quantile, "grass78/bin"),
    # ExecutableProduct("r.out.vtk", :r_out_vtk, "grass78/bin"),
    # ExecutableProduct("r.walk", :r_walk, "grass78/bin"),
    # ExecutableProduct("r.gwflow", :r_gwflow, "grass78/bin"),
    # ExecutableProduct("r.li.edgedensity", :r_li_edgedensity, "grass78/bin"),
    # ExecutableProduct("r.texture", :r_texture, "grass78/bin"),
    # ExecutableProduct("r.resamp.stats", :r_resamp_stats, "grass78/bin"),
    # ExecutableProduct("r.fill.stats", :r_fill_stats, "grass78/bin"),
    # ExecutableProduct("r.li.shannon", :r_li_shannon, "grass78/bin"),
    # ExecutableProduct("r.in.gdal", :r_in_gdal, "grass78/bin"),
    # ExecutableProduct("r.buildvrt", :r_buildvrt, "grass78/bin"),
    # ExecutableProduct("r.grow.distance", :r_grow_distance, "grass78/bin"),
    # ExecutableProduct("r.circle", :r_circle, "grass78/bin"),
    # ExecutableProduct("r.in.png", :r_in_png, "grass78/bin"),
    # ExecutableProduct("r.li.padrange", :r_li_padrange, "grass78/bin"),
    # ExecutableProduct("r.spread", :r_spread, "grass78/bin"),
    # ExecutableProduct("r.out.gridatb", :r_out_gridatb, "grass78/bin"),
    # ExecutableProduct("r.info", :r_info, "grass78/bin"),
    # ExecutableProduct("r.clump", :r_clump, "grass78/bin"),
    # ExecutableProduct("r.li.pielou", :r_li_pielou, "grass78/bin"),
    # ExecutableProduct("r.resamp.filter", :r_resamp_filter, "grass78/bin"),
    # ExecutableProduct("r.surf.fractal", :r_surf_fractal, "grass78/bin"),
    # ExecutableProduct("r.timestamp", :r_timestamp, "grass78/bin"),
    # ExecutableProduct("r.in.xyz", :r_in_xyz, "grass78/bin"),
    # ExecutableProduct("r.li.padcv", :r_li_padcv, "grass78/bin"),
    # ExecutableProduct("r.solute.transport", :r_solute_transport, "grass78/bin"),
    # ExecutableProduct("r.in.poly", :r_in_poly, "grass78/bin"),
    # ExecutableProduct("r.out.bin", :r_out_bin, "grass78/bin"),
    # ExecutableProduct("r.compress", :r_compress, "grass78/bin"),
    # ExecutableProduct("r.region", :r_region, "grass78/bin"),
    # ExecutableProduct("r.surf.gauss", :r_surf_gauss, "grass78/bin"),
    # ExecutableProduct("r.recode", :r_recode, "grass78/bin"),
    # ExecutableProduct("r.ros", :r_ros, "grass78/bin"),
    # ExecutableProduct("r.latlong", :r_latlong, "grass78/bin"),
    # ExecutableProduct("r.tile", :r_tile, "grass78/bin"),
    # ExecutableProduct("r.quantile", :r_quantile, "grass78/bin"),
    # ExecutableProduct("r.resamp.rst", :r_resamp_rst, "grass78/bin"),
    # ExecutableProduct("r.topidx", :r_topidx, "grass78/bin"),
    # ExecutableProduct("r.lake", :r_lake, "grass78/bin"),
    # ExecutableProduct("r.composite", :r_composite, "grass78/bin"),
    # ExecutableProduct("r.kappa", :r_kappa, "grass78/bin"),
    # ExecutableProduct("r.li.mpa", :r_li_mpa, "grass78/bin"),
    # ExecutableProduct("r.out.ppm", :r_out_ppm, "grass78/bin"),
    # ExecutableProduct("r.sunmask", :r_sunmask, "grass78/bin"),
    # ExecutableProduct("r.to.vect", :r_to_vect, "grass78/bin"),
    # ExecutableProduct("r.random.surface", :r_random_surface, "grass78/bin"),
    # ExecutableProduct("r.li.simpson", :r_li_simpson, "grass78/bin"),
    # ExecutableProduct("r.out.pov", :r_out_pov, "grass78/bin"),
    # ExecutableProduct("r.what", :r_what, "grass78/bin"),
    # ExecutableProduct("r.li.patchdensity", :r_li_patchdensity, "grass78/bin"),
    # ExecutableProduct("r.flow", :r_flow, "grass78/bin"),
    # ExecutableProduct("r.surf.idw", :r_surf_idw, "grass78/bin"),
    # ExecutableProduct("r.spreadpath", :r_spreadpath, "grass78/bin"),
    # ExecutableProduct("r.fill.dir", :r_fill_dir, "grass78/bin"),
    # ExecutableProduct("r.in.ascii", :r_in_ascii, "grass78/bin"),
    # ExecutableProduct("r.resamp.bspline", :r_resamp_bspline, "grass78/bin"),
    # ExecutableProduct("r.uslek", :r_uslek, "grass78/bin"),
    # ExecutableProduct("r.mfilter", :r_mfilter, "grass78/bin"),
    # ExecutableProduct("r.series.accumulate", :r_series_accumulate, "grass78/bin"),
    # ExecutableProduct("r.stats.zonal", :r_stats_zonal, "grass78/bin"),
    # ExecutableProduct("r.li.shape", :r_li_shape, "grass78/bin"),
    # ExecutableProduct("r.regression.line", :r_regression_line, "grass78/bin"),
    # ExecutableProduct("r.sim.sediment", :r_sim_sediment, "grass78/bin"),
    # ExecutableProduct("r.transect", :r_transect, "grass78/bin"),
    # ExecutableProduct("r.in.gridatb", :r_in_gridatb, "grass78/bin"),
    # ExecutableProduct("r.coin", :r_coin, "grass78/bin"),
    # ExecutableProduct("r.slope.aspect", :r_slope_aspect, "grass78/bin"),
    # ExecutableProduct("r.covar", :r_covar, "grass78/bin"),
    # ExecutableProduct("r.horizon", :r_horizon, "grass78/bin"),
    # ExecutableProduct("r.out.png", :r_out_png, "grass78/bin"),
    # ExecutableProduct("r.neighbors", :r_neighbors, "grass78/bin"),
    # ExecutableProduct("r.li.mps", :r_li_mps, "grass78/bin"),
    # ExecutableProduct("r.contour", :r_contour, "grass78/bin"),
    # ExecutableProduct("r.in.mat", :r_in_mat, "grass78/bin"),
    # ExecutableProduct("r.relief", :r_relief, "grass78/bin"),
    # ExecutableProduct("r.path", :r_path, "grass78/bin"),
    # ExecutableProduct("r.series.interp", :r_series_interp, "grass78/bin"),
    # ExecutableProduct("r.proj", :r_proj, "grass78/bin"),
    # ExecutableProduct("r.li.patchnum", :r_li_patchnum, "grass78/bin"),
    # ExecutableProduct("r.usler", :r_usler, "grass78/bin"),
    # ExecutableProduct("r.to.rast3", :r_to_rast3, "grass78/bin"),
    # ExecutableProduct("r.quant", :r_quant, "grass78/bin"),
    # ExecutableProduct("r.report", :r_report, "grass78/bin"),
    # ExecutableProduct("r.colors.out", :r_colors_out, "grass78/bin"),
    # ExecutableProduct("r.support.stats", :r_support_stats, "grass78/bin"),
    # ExecutableProduct("r.li.renyi", :r_li_renyi, "grass78/bin"),
    # ExecutableProduct("r.external.out", :r_external_out, "grass78/bin"),
    # ExecutableProduct("r.random", :r_random, "grass78/bin"),
    # ExecutableProduct("r.stats", :r_stats, "grass78/bin"),
    # ExecutableProduct("r.surf.random", :r_surf_random, "grass78/bin"),
    # ExecutableProduct("r.resamp.interp", :r_resamp_interp, "grass78/bin"),
    # ExecutableProduct("r.sun", :r_sun, "grass78/bin"),
    # ExecutableProduct("r.describe", :r_describe, "grass78/bin"),
    # ExecutableProduct("r.out.gdal", :r_out_gdal, "grass78/bin"),
    # ExecutableProduct("r.out.mpeg", :r_out_mpeg, "grass78/bin"),

    # ExecutableProduct("v.in.db", :v_in_db, "grass78/bin"),
    # ExecutableProduct("v.lidar.growing", :v_lidar_growing, "grass78/bin"),
    # ExecutableProduct("v.out.vtk", :v_out_vtk, "grass78/bin"),
    # ExecutableProduct("v.surf.bspline", :v_surf_bspline, "grass78/bin"),
    # ExecutableProduct("v.net.distance", :v_net_distance, "grass78/bin"),
    # ExecutableProduct("v.support", :v_support, "grass78/bin"),
    # ExecutableProduct("v.univar", :v_univar, "grass78/bin"),
    # ExecutableProduct("v.out.ascii", :v_out_ascii, "grass78/bin"),
    # ExecutableProduct("v.reclass", :v_reclass, "grass78/bin"),
    # ExecutableProduct("v.surf.idw", :v_surf_idw, "grass78/bin"),
    # ExecutableProduct("v.net.visibility", :v_net_visibility, "grass78/bin"),
    # ExecutableProduct("v.build.polylines", :v_build_polylines, "grass78/bin"),
    # ExecutableProduct("v.outlier", :v_outlier, "grass78/bin"),
    # ExecutableProduct("v.split", :v_split, "grass78/bin"),
    # ExecutableProduct("v.vol.rst", :v_vol_rst, "grass78/bin"),
    # ExecutableProduct("v.parallel", :v_parallel, "grass78/bin"),
    # ExecutableProduct("v.net.steiner", :v_net_steiner, "grass78/bin"),
    # ExecutableProduct("v.net.path", :v_net_path, "grass78/bin"),
    # ExecutableProduct("v.transform", :v_transform, "grass78/bin"),
    # ExecutableProduct("v.kcv", :v_kcv, "grass78/bin"),
    # ExecutableProduct("v.category", :v_category, "grass78/bin"),
    # ExecutableProduct("v.profile", :v_profile, "grass78/bin"),
    # ExecutableProduct("v.drape", :v_drape, "grass78/bin"),
    # ExecutableProduct("v.out.ogr", :v_out_ogr, "grass78/bin"),
    # ExecutableProduct("v.db.select", :v_db_select, "grass78/bin"),
    # ExecutableProduct("v.what.rast", :v_what_rast, "grass78/bin"),
    # ExecutableProduct("v.neighbors", :v_neighbors, "grass78/bin"),
    # ExecutableProduct("v.net.salesman", :v_net_salesman, "grass78/bin"),
    # ExecutableProduct("v.normal", :v_normal, "grass78/bin"),
    # ExecutableProduct("v.in.dxf", :v_in_dxf, "grass78/bin"),
    # ExecutableProduct("v.colors", :v_colors, "grass78/bin"),
    # ExecutableProduct("v.perturb", :v_perturb, "grass78/bin"),
    # ExecutableProduct("v.decimate", :v_decimate, "grass78/bin"),
    # ExecutableProduct("v.extrude", :v_extrude, "grass78/bin"),
    # ExecutableProduct("v.to.db", :v_to_db, "grass78/bin"),
    # ExecutableProduct("v.kernel", :v_kernel, "grass78/bin"),
    # ExecutableProduct("v.random", :v_random, "grass78/bin"),
    # ExecutableProduct("v.net.flow", :v_net_flow, "grass78/bin"),
    # ExecutableProduct("v.lrs.segment", :v_lrs_segment, "grass78/bin"),
    # ExecutableProduct("v.label", :v_label, "grass78/bin"),
    # ExecutableProduct("v.lrs.create", :v_lrs_create, "grass78/bin"),
    # ExecutableProduct("v.mkgrid", :v_mkgrid, "grass78/bin"),
    # ExecutableProduct("v.generalize", :v_generalize, "grass78/bin"),
    # ExecutableProduct("v.net.alloc", :v_net_alloc, "grass78/bin"),
    # ExecutableProduct("v.vect.stats", :v_vect_stats, "grass78/bin"),
    # ExecutableProduct("v.in.ogr", :v_in_ogr, "grass78/bin"),
    # ExecutableProduct("v.lrs.label", :v_lrs_label, "grass78/bin"),
    # ExecutableProduct("v.to.points", :v_to_points, "grass78/bin"),
    # ExecutableProduct("v.lidar.edgedetection", :v_lidar_edgedetection, "grass78/bin"),
    # ExecutableProduct("v.out.pov", :v_out_pov, "grass78/bin"),
    # ExecutableProduct("v.db.connect", :v_db_connect, "grass78/bin"),
    # ExecutableProduct("v.colors.out", :v_colors_out, "grass78/bin"),
    # ExecutableProduct("v.distance", :v_distance, "grass78/bin"),
    # ExecutableProduct("v.sample", :v_sample, "grass78/bin"),
    # ExecutableProduct("v.timestamp", :v_timestamp, "grass78/bin"),
    # ExecutableProduct("v.what", :v_what, "grass78/bin"),
    # ExecutableProduct("v.extract", :v_extract, "grass78/bin"),
    # ExecutableProduct("v.buffer", :v_buffer, "grass78/bin"),
    # ExecutableProduct("v.net.iso", :v_net_iso, "grass78/bin"),
    # ExecutableProduct("v.net.components", :v_net_components, "grass78/bin"),
    # ExecutableProduct("v.net.centrality", :v_net_centrality, "grass78/bin"),
    # ExecutableProduct("v.surf.rst", :v_surf_rst, "grass78/bin"),
    # ExecutableProduct("v.to.3d", :v_to_3d, "grass78/bin"),
    # ExecutableProduct("v.qcount", :v_qcount, "grass78/bin"),
    # ExecutableProduct("v.edit", :v_edit, "grass78/bin"),
    # ExecutableProduct("v.rectify", :v_rectify, "grass78/bin"),
    # ExecutableProduct("v.label.sa", :v_label_sa, "grass78/bin"),
    # ExecutableProduct("v.net.allpairs", :v_net_allpairs, "grass78/bin"),
    # ExecutableProduct("v.external.out", :v_external_out, "grass78/bin"),
    # ExecutableProduct("v.what.rast3", :v_what_rast3, "grass78/bin"),
    # ExecutableProduct("v.in.ascii", :v_in_ascii, "grass78/bin"),
    # ExecutableProduct("v.patch", :v_patch, "grass78/bin"),
    # ExecutableProduct("v.clean", :v_clean, "grass78/bin"),
    # ExecutableProduct("v.info", :v_info, "grass78/bin"),
    # ExecutableProduct("v.net.timetable", :v_net_timetable, "grass78/bin"),
    # ExecutableProduct("v.net.bridge", :v_net_bridge, "grass78/bin"),
    # ExecutableProduct("v.out.dxf", :v_out_dxf, "grass78/bin"),
    # ExecutableProduct("v.to.rast", :v_to_rast, "grass78/bin"),
    # ExecutableProduct("v.lrs.where", :v_lrs_where, "grass78/bin"),
    # ExecutableProduct("v.delaunay", :v_delaunay, "grass78/bin"),
    # ExecutableProduct("v.in.region", :v_in_region, "grass78/bin"),
    # ExecutableProduct("v.external", :v_external, "grass78/bin"),
    # ExecutableProduct("v.lidar.correction", :v_lidar_correction, "grass78/bin"),
    # ExecutableProduct("v.class", :v_class, "grass78/bin"),
    # ExecutableProduct("v.segment", :v_segment, "grass78/bin"),
    # ExecutableProduct("v.proj", :v_proj, "grass78/bin"),
    # ExecutableProduct("v.cluster", :v_cluster, "grass78/bin"),
    # ExecutableProduct("v.overlay", :v_overlay, "grass78/bin"),
    # ExecutableProduct("v.type", :v_type, "grass78/bin"),
    # ExecutableProduct("v.voronoi", :v_voronoi, "grass78/bin"),
    # ExecutableProduct("v.hull", :v_hull, "grass78/bin"),
    # ExecutableProduct("v.out.svg", :v_out_svg, "grass78/bin"),
    # ExecutableProduct("v.net.connectivity", :v_net_connectivity, "grass78/bin"),
    # ExecutableProduct("v.build", :v_build, "grass78/bin"),
    # ExecutableProduct("v.net.spanningtree", :v_net_spanningtree, "grass78/bin"),
    # ExecutableProduct("v.net", :v_net, "grass78/bin"),
    # ExecutableProduct("v.to.rast3", :v_to_rast3, "grass78/bin"),
    # ExecutableProduct("v.select", :v_select, "grass78/bin"),

    # ExecutableProduct("i.atcorr", :i_atcorr, "grass78/bin"),
    # ExecutableProduct("i.eb.hsebal01", :i_eb_hsebal01, "grass78/bin"),
    # ExecutableProduct("i.ifft", :i_ifft, "grass78/bin"),
    # ExecutableProduct("i.biomass", :i_biomass, "grass78/bin"),
    # ExecutableProduct("i.gensig", :i_gensig, "grass78/bin"),
    # ExecutableProduct("i.eb.soilheatflux", :i_eb_soilheatflux, "grass78/bin"),
    # ExecutableProduct("i.albedo", :i_albedo, "grass78/bin"),
    # ExecutableProduct("i.evapo.mh", :i_evapo_mh, "grass78/bin"),
    # ExecutableProduct("i.eb.eta", :i_eb_eta, "grass78/bin"),
    # ExecutableProduct("i.aster.toar", :i_aster_toar, "grass78/bin"),
    # ExecutableProduct("i.target", :i_target, "grass78/bin"),
    # ExecutableProduct("i.rgb.his", :i_rgb_his, "grass78/bin"),
    # ExecutableProduct("i.evapo.pt", :i_evapo_pt, "grass78/bin"),
    # ExecutableProduct("i.find", :i_find, "grass78/etc"),
    # ExecutableProduct("i.ortho.camera", :i_ortho_camera, "grass78/bin"),
    # ExecutableProduct("i.evapo.pm", :i_evapo_pm, "grass78/bin"),
    # ExecutableProduct("i.eb.evapfr", :i_eb_evapfr, "grass78/bin"),
    # ExecutableProduct("i.ortho.elev", :i_ortho_elev, "grass78/bin"),
    # ExecutableProduct("i.group", :i_group, "grass78/bin"),
    # ExecutableProduct("i.ortho.rectify", :i_ortho_rectify, "grass78/bin"),
    # ExecutableProduct("i.ortho.init", :i_ortho_init, "grass78/bin"),
    # ExecutableProduct("i.eb.netrad", :i_eb_netrad, "grass78/bin"),
    # ExecutableProduct("i.segment", :i_segment, "grass78/bin"),
    # ExecutableProduct("i.vi", :i_vi, "grass78/bin"),
    # ExecutableProduct("i.maxlik", :i_maxlik, "grass78/bin"),
    # ExecutableProduct("i.evapo.time", :i_evapo_time, "grass78/bin"),
    # ExecutableProduct("i.ortho.transform", :i_ortho_transform, "grass78/bin"),
    # ExecutableProduct("i.topo.corr", :i_topo_corr, "grass78/bin"),
    # ExecutableProduct("i.modis.qc", :i_modis_qc, "grass78/bin"),
    # ExecutableProduct("i.cca", :i_cca, "grass78/bin"),
    # ExecutableProduct("i.his.rgb", :i_his_rgb, "grass78/bin"),
    # ExecutableProduct("i.fft", :i_fft, "grass78/bin"),
    # ExecutableProduct("i.landsat.toar", :i_landsat_toar, "grass78/bin"),
    # ExecutableProduct("i.gensigset", :i_gensigset, "grass78/bin"),
    # ExecutableProduct("i.rectify", :i_rectify, "grass78/bin"),
    # ExecutableProduct("i.pca", :i_pca, "grass78/bin"),
    # ExecutableProduct("i.emissivity", :i_emissivity, "grass78/bin"),
    # ExecutableProduct("i.landsat.acca", :i_landsat_acca, "grass78/bin"),
    # ExecutableProduct("i.zc", :i_zc, "grass78/bin"),
    # ExecutableProduct("i.ortho.photo", :i_ortho_photo, "grass78/bin"),
    # ExecutableProduct("i.smap", :i_smap, "grass78/bin"),
    # ExecutableProduct("i.cluster", :i_cluster, "grass78/bin"),
    # ExecutableProduct("i.ortho.target", :i_ortho_target, "grass78/bin"),

    # ExecutableProduct("r3.in.ascii", :r3_in_ascii, "grass78/bin"),
    # ExecutableProduct("r3.in.bin", :r3_in_bin, "grass78/bin"),
    # ExecutableProduct("r3.to.rast", :r3_to_rast, "grass78/bin"),
    # ExecutableProduct("r3.timestamp", :r3_timestamp, "grass78/bin"),
    # ExecutableProduct("r3.support", :r3_support, "grass78/bin"),
    # ExecutableProduct("r3.retile", :r3_retile, "grass78/bin"),
    # ExecutableProduct("r3.colors.out", :r3_colors_out, "grass78/bin"),
    # ExecutableProduct("r3.out.vtk", :r3_out_vtk, "grass78/bin"),
    # ExecutableProduct("r3.in.v5d", :r3_in_v5d, "grass78/bin"),
    # ExecutableProduct("r3.mapcalc", :r3_mapcalc, "grass78/bin"),
    # ExecutableProduct("r3.info", :r3_info, "grass78/bin"),
    # ExecutableProduct("r3.flow", :r3_flow, "grass78/bin"),
    # ExecutableProduct("r3.neighbors", :r3_neighbors, "grass78/bin"),
    # ExecutableProduct("r3.out.ascii", :r3_out_ascii, "grass78/bin"),
    # ExecutableProduct("r3.gwflow", :r3_gwflow, "grass78/bin"),
    # ExecutableProduct("r3.mask", :r3_mask, "grass78/bin"),
    # ExecutableProduct("r3.null", :r3_null, "grass78/bin"),
    # ExecutableProduct("r3.out.v5d", :r3_out_v5d, "grass78/bin"),
    # ExecutableProduct("r3.out.bin", :r3_out_bin, "grass78/bin"),
    # ExecutableProduct("r3.colors", :r3_colors, "grass78/bin"),
    # ExecutableProduct("r3.mkdspf", :r3_mkdspf, "grass78/bin"),
    # ExecutableProduct("r3.gradient", :r3_gradient, "grass78/bin"),
    # ExecutableProduct("r3.univar", :r3_univar, "grass78/bin"),
    # ExecutableProduct("r3.stats", :r3_stats, "grass78/bin"),
    # ExecutableProduct("r3.cross.rast", :r3_cross_rast, "grass78/bin"),

    # ExecutableProduct("d.info", :d_info, "grass78/bin"),
    # ExecutableProduct("d.rast.num", :d_rast_num, "grass78/bin"),
    # ExecutableProduct("d.text", :d_text, "grass78/bin"),
    # ExecutableProduct("d.barscale", :d_barscale, "grass78/bin"),
    # ExecutableProduct("d.rast.arrow", :d_rast_arrow, "grass78/bin"),
    # ExecutableProduct("d.font", :d_font, "grass78/bin"),
    # ExecutableProduct("d.geodesic", :d_geodesic, "grass78/bin"),
    # ExecutableProduct("d.rhumbline", :d_rhumbline, "grass78/bin"),
    # ExecutableProduct("d.vect", :d_vect, "grass78/bin"),
    # ExecutableProduct("d.northarrow", :d_northarrow, "grass78/bin"),
    # ExecutableProduct("d.linegraph", :d_linegraph, "grass78/bin"),
    # ExecutableProduct("d.path", :d_path, "grass78/bin"),
    # ExecutableProduct("d.where", :d_where, "grass78/bin"),
    # ExecutableProduct("d.vect.thematic", :d_vect_thematic, "grass78/bin"),
    # ExecutableProduct("d.fontlist", :d_fontlist, "grass78/bin"),
    # ExecutableProduct("d.vect.chart", :d_vect_chart, "grass78/bin"),
    # ExecutableProduct("d.mon", :d_mon, "grass78/bin"),
    # ExecutableProduct("d.legend.vect", :d_legend_vect, "grass78/bin"),
    # ExecutableProduct("d.erase", :d_erase, "grass78/bin"),
    # ExecutableProduct("d.labels", :d_labels, "grass78/bin"),
    # ExecutableProduct("d.histogram", :d_histogram, "grass78/bin"),
    # ExecutableProduct("d.colorlist", :d_colorlist, "grass78/bin"),
    # ExecutableProduct("d.his", :d_his, "grass78/bin"),
    # ExecutableProduct("d.title", :d_title, "grass78/bin"),
    # ExecutableProduct("d.profile", :d_profile, "grass78/bin"),
    # ExecutableProduct("d.grid", :d_grid, "grass78/bin"),
    # ExecutableProduct("d.rgb", :d_rgb, "grass78/bin"),
    # ExecutableProduct("d.graph", :d_graph, "grass78/bin"),
    # ExecutableProduct("d.legend", :d_legend, "grass78/bin"),
    # ExecutableProduct("d.rast", :d_rast, "grass78/bin"),
    # ExecutableProduct("d.colortable", :d_colortable, "grass78/bin"),

    # ExecutableProduct("ps.map", :ps_map, "grass78/bin"),

    # ExecutableProduct("test.r3flow", :test_r3flow, "grass78/bin"),
    # ExecutableProduct("test.raster3d.lib", :test_raster3d_lib, "grass78/bin"),

    # ExecutableProduct("echo", :echo, "grass78/etc"),
    # ExecutableProduct("run", :grass_run, "grass78/etc"),
    # ExecutableProduct("lock", :grass_lock, "grass78/etc"),
    # ExecutableProduct("current_time_s_ms", :current_time_s_ms, "grass78/etc"),
    # ExecutableProduct("clean_temp", :clean_temp, "grass78/etc"),

    # ExecutableProduct("m.cogo", :m_cogo, "grass78/bin"),
    # ExecutableProduct("m.transform", :m_transform, "grass78/bin"),
    # ExecutableProduct("m.measure", :m_measure, "grass78/bin"),
    # ExecutableProduct("m.nviz.script", :m_nviz_script, "grass78/bin"),

    # ExecutableProduct("t.connect", :t_connect, "grass78/bin"),

    # ExecutableProduct("dbf", :dbf, "grass78/driver/db"),
    # ExecutableProduct("sqlite", :sqlite, "grass78/driver/db"),
    # ExecutableProduct("ogr", :ogr, "grass78/driver/db"),
    # ExecutableProduct("seg", :seg, "grass78/etc/r.watershed"),
    # ExecutableProduct("ram", :ram, "grass78/etc/r.watershed"),
    # ExecutableProduct("vector", :vector, "grass78/etc/lister"),
    # ExecutableProduct("cell", :cell, "grass78/etc/lister")

]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GDAL_jll", uuid="a7073274-a066-55f0-b90d-d619367d196c"))
    Dependency(PackageSpec(name="PROJ_jll", uuid="58948b4f-47e0-5654-a9ad-f609743f8632"))
    Dependency(PackageSpec(name="Cairo_jll", uuid="83423d85-b0ee-5818-9007-b63ccbeb887a"))
    Dependency(PackageSpec(name="FreeType2_jll", uuid="d7e528f0-a631-5988-bf34-fe36492bcfd7"))
    Dependency(PackageSpec(name="LAPACK_jll", uuid="51474c39-65e3-53ba-86ba-03b1b862ec14"))
    Dependency(PackageSpec(name="FFTW_jll", uuid="f5851436-0d7a-5f13-b9de-f02708fd171a"))
    Dependency(PackageSpec(name="Libiconv_jll", uuid="94ce4f54-9a6c-5748-9c1c-f9c7231a4531"))
    BuildDependency(PackageSpec(name="Xorg_xorgproto_jll", uuid="c4d99508-4286-5418-9131-c86396af500b"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"7.1.0", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).
@recipe function f(
    parshapes::NamedTupleShape,
    samples::DensitySampleVector,
    param::Symbol;
    intervals = standard_confidence_vals,
    bins=200,
    normalize = true,
    colors = standard_colors,
    intervallabels = [],
    mean = false,
    std_dev = false,
    globalmode = false,
    localmode = true,
    filter = true
)
    i = findfirst(x -> x == param, keys(parshapes))

    @series begin
        intervals --> intervals
        bins --> bins
        normalize --> normalize
        colors --> colors
        intervallabels --> intervallabels
        mean --> mean
        std_dev -->  std_dev
        globalmode --> globalmode
        localmode --> localmode
        filter --> filter

        samples, i
    end

end

@recipe function f(
    posterior::PosteriorDensity,
    samples::DensitySampleVector,
    param::Symbol;
    intervals = standard_confidence_vals,
    bins=200,
    normalize = true,
    colors = standard_colors,
    intervallabels = [],
    mean = false,
    std_dev = false,
    globalmode = false,
    localmode = true,
    filter = true
)
    i = findfirst(x -> x == param, keys(params_shape(posterior)))

    @series begin
        intervals --> intervals
        bins --> bins
        normalize --> normalize
        colors --> colors
        intervallabels --> intervallabels
        mean --> mean
        std_dev -->  std_dev
        globalmode --> globalmode
        localmode --> localmode
        filter --> filter

        samples, i
    end

end

@recipe function f(
    samples::DensitySampleVector,
    param::Integer;
    intervals = standard_confidence_vals,
    bins=200,
    normalize = true,
    colors = standard_colors,
    intervallabels = [],
    mean = false,
    std_dev = false,
    globalmode = false,
    localmode = true,
    filter = true
)
    if filter
        samples = BAT.drop_low_weight_samples(samples)
    end

    #seriestype = get(plotattributes, :seriestype, :smallest_intervals)

    hist = fit(Histogram, flatview(samples.params)[param, :], FrequencyWeights(samples.weight), nbins = bins, closed = :left)
    normalize ? hist=StatsBase.normalize(hist) : nothing
    line_height = maximum(hist.weights)*1.03

    orientation = get(plotattributes, :orientation, :vertical)
    (orientation != :vertical) ? swap=true : swap = false

    stats = MCMCBasicStats(samples)

    if swap
        yguide --> "\$\\theta_$(param)\$"
        xguide --> "\$p(\\theta_$(param))\$"
    else
        yguide --> "\$p(\\theta_$(param))\$"
        xguide --> "\$\\theta_$(param)\$"
    end

    @series begin
        seriestype --> :smallest_intervals

        intervals --> intervals
        bins --> bins
        normalize --> normalize
        colors --> colors
        intervallabels --> intervallabels

        hist, param
    end

    mean_options = convert_to_options(mean)
    globalmode_options = convert_to_options(globalmode)
    localmode_options = convert_to_options(localmode)
    stddev_options = convert_to_options(std_dev)

    if stddev_options != ()
        Σ_all = stats.param_stats.cov
        Σ = Σ_all[param, param]
        dev = sqrt(Σ)
        meanvalue = stats.param_stats.mean[param]
        @series begin
            seriestype := :shape
            label := get(stddev_options, "label", "std. dev.")
            linewidth := 0
            fillcolor := get(stddev_options, "fillcolor", :grey)
            fillalpha := get(stddev_options, "fillalpha", 0.5)
            uncertaintyband(meanvalue, dev, line_height, swap=swap)
        end
    end

    if mean_options != ()
        meanvalue = stats.param_stats.mean[param]
        @series begin
            seriestype := :line
            label := get(mean_options, "label", "mean") # : $(@sprintf("%.2f", meanvalue))
            linestyle := get(mean_options, "linestyle", :solid)
            linecolor := get(mean_options, "linecolor", :dimgrey)
            linewidth := get(mean_options, "linewidth", 1)
            alpha := get(mean_options, "alpha", 1)

            line(meanvalue, line_height, swap=swap)
        end
    end

   if globalmode_options != ()
        globalmode_value = stats.mode[param]
         @series begin
            seriestype := :line
            label := get(globalmode_options, "label", "global mode") #: $(@sprintf("%.2f", globalmode_value))
            linestyle := get(globalmode_options, "linestyle", :dash)
            linecolor := get(globalmode_options, "linecolor", :black)
            linewidth := get(globalmode_options, "linewidth", 1)
            alpha := get(globalmode_options, "alpha", 1)

            line(globalmode_value, line_height, swap=swap)
        end
    end

    if localmode_options != ()
        localmode_values = calculate_localmode(hist)
        for (i, l) in enumerate(localmode_values)
         @series begin
            seriestype := :line
            if i==1 && length(localmode_values)==1
                label := get(localmode_options, "label", "local mode") #: $(@sprintf("%.2f", l))
            elseif i ==1
                label := get(localmode_options, "label", "local modes")
            else
                label :=""
            end
            linestyle := get(localmode_options, "linestyle", :dot)
            linecolor := get(localmode_options, "linecolor", :black)
            linewidth := get(localmode_options, "linewidth", 1)
            alpha := get(localmode_options, "alpha", 1)
            line(l, line_height, swap=swap)
            end
        end
    end

end

function line(pos, height; swap=false)
    if swap
        return  [(0, pos), (height, pos)]
    else
        return  [(pos, 0), (pos, height)]
    end
end

function convert_to_options(arg::Dict)
    return arg
end

function convert_to_options(arg::Bool)
    arg ? Dict() : ()
end

function uncertaintyband(m, u, h; swap=false)
    if _plots_module() != nothing
        if swap
            return _plots_module().Shape([0,0,h,h], [m-u,m+u,m+u,m-u])
        else
            return _plots_module().Shape([m-u,m+u,m+u,m-u], [0,0,h,h])
        end
    else
        ()
    end
end


text:
# # Minimization of simple loss functions
#
# First we load the required packages.

using
    Distributions,
    LinearAlgebra,
    Random,
    Plots

using
    EnsembleKalmanProcesses.EnsembleKalmanProcessModule,
    EnsembleKalmanProcesses.ParameterDistributionStorage

# ## Loss function with single minimum
#
# Here, we minimize the loss function
# ```math
# G₁(u) = \|u - u_*\| ,
# ```
# where ``u`` is a 2-vector of parameters and ``u_*`` is given; here ``u_* = (-1, 1)``.
u_min = [0.2, 0.05] #minimizer is 0,0
function G₁(u)

    if (u[1]>=0) && (u[2] >=0)
        return sqrt((u[1]-u_min[1])^2 + (u[2]-u_min[2])^2)
    else
        return NaN
    end
end

# We set the seed for pseudo-random number generator for reproducibility.
rng_seed = 41
Random.seed!(rng_seed)

# We set a stabilization level, which can aid the algorithm convergence
dim_output = 1
stabilization_level = 1e-4
Γ_stabilization = stabilization_level * Matrix(I, dim_output, dim_output)

# The functional is positive so to minimize it we may set the target to be 0,
G_target  = [0]

# ### Prior distributions
#
# As we work with a Bayesian method, we define a prior. This will behave like an "initial guess"
# for the likely region of parameter space we expect the solution to live in.
prior_distributions = [Parameterized(Normal(0, 1)), Parameterized(Normal(0, 1))]

constraints = [[no_constraint()], [no_constraint()]]

parameter_names = ["u1", "u2"]

prior = ParameterDistribution(prior_distributions, constraints, parameter_names)

# ### Calibration
#
# We choose the number of ensemble members and the number of iterations of the algorithm
N_ensemble   = 30
N_iterations = 10
nothing # hide

# The initial ensemble is constructed by sampling the prior
initial_ensemble =
    EnsembleKalmanProcessModule.construct_initial_ensemble(prior, N_ensemble;
                                                           rng_seed=rng_seed)
# We then initialize the Ensemble Kalman Process algorithm, with the initial ensemble, the
# target, the stabilization and the process type (for EKI this is `Inversion`, initialized
# with `Inversion()`).
ensemble_kalman_process =
    EnsembleKalmanProcessModule.EnsembleKalmanProcess(initial_ensemble, G_target,
                                                      Γ_stabilization, Inversion())

# Then we calibrate by *(i)* obtaining the parameters, *(ii)* calculate the loss function on
# the parameters (and concatenate), and last *(iii)* generate a new set of parameters using
# the model outputs:
for i in 1:N_iterations
    params_i = get_u_final(ensemble_kalman_process)

    g_ens = hcat([G₁(params_i[:, i]) for i in 1:N_ensemble]...)
    println("Iteration: ",i," number of particle failures: ", sum([1 for i = 1:size(g_ens,2) if isnan(g_ens[1,i])]), "/", size(g_ens,2))

    EnsembleKalmanProcessModule.update_ensemble!(ensemble_kalman_process, g_ens)

end

# and visualize the results:
u_init = get_u_prior(ensemble_kalman_process)

anim_unique_minimum = @animate for i in 1:N_iterations
    u_i = get_u(ensemble_kalman_process, i)
    g_i = get_g(ensemble_kalman_process, i)
    failed_particles = [i for i = 1:size(g_i,2) if isnan(g_i[1,i])]
    successful_particles = filter(x-> !(x in failed_particles), collect(1:size(g_i,2)))
    u_i_succ = u_i[:,successful_particles]
    u_i_fail = u_i[:,failed_particles]

    plot([u_min[1]], [u_min[2]],
          seriestype = :scatter,
         markershape = :star5,
          markersize = 11,
         markercolor = :green,
               label = "optimum u⋆"
         )

    plot!(u_i_succ[1, :], u_i_succ[2, :],
             seriestype = :scatter,
                  xlims = extrema(u_init[1, :]),
                  ylims = extrema(u_init[2, :]),
                 xlabel = "u₁",
                 ylabel = "u₂",
             markersize = 5,
            markeralpha = 0.6,
            markercolor = :blue,
                  label = "successes",
                  title = "EKI iteration = " * string(i)
          )
        plot!(u_i_fail[1, :], u_i_fail[2, :],
             seriestype = :scatter,
                  xlims = extrema(u_init[1, :]),
                  ylims = extrema(u_init[2, :]),
             markersize = 5,
            markeralpha = 0.6,
              markercolor = :red,
              markershape = :x,
                  label = "failures",
         )

end
# The results show that the minimizer of ``G_1`` is ``u=u_*``.

gif(anim_unique_minimum, "unique_minimum.gif", fps = 1) # hide


text:
"""
function euler_fluxes(UL,UR)
    UL, UR = states = (rho,u,v,beta)
    2D version.  Figure out multiple dispatch version for 3D?

"""

# assumes field ordering: (rhoL,uL,...,betaL) = UL, (rhoR,uR,...,betaR) = UR
function euler_fluxes(UL,UR)
    rhoL = first(UL); betaL = last(UL)
    rhoR = first(UR); betaR = last(UR)
    logL = (log.(rhoL),log.(betaL))
    logR = (log.(rhoR),log.(betaR))
    return euler_fluxes(UL,UR,logL,logR)
end

# dispatches to 2D/3D
function euler_fluxes(UL,UR,logL,logR)
    return euler_fluxes(UL...,UR...,logL...,logR...)
end

# 2d version
function euler_fluxes(rhoL,uL,vL,betaL,rhoR,uR,vR,betaR,
                      rhologL,betalogL,rhologR,betalogR)

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))
    vavg   = (@. .5*(vL+vR))

    unorm = (@. uL*uR + vL*vR)
    pa    = (@. rhoavg/(betaL+betaR))
    f4aux = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. FxS1*vavg)
    FxS4 = (@. f4aux*uavg)

    FyS1 = (@. rholog*vavg)
    FyS2 = (@. FxS3)
    FyS3 = (@. FyS1*vavg + pa)
    FyS4 = (@. f4aux*vavg)
    return (FxS1,FxS2,FxS3,FxS4),(FyS1,FyS2,FyS3,FyS4)
end

# 3d version
function euler_fluxes(rhoL,uL,vL,wL,betaL,rhoR,uR,vR,wR,betaR,
                      rhologL,betalogL,rhologR,betalogR)

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))
    vavg   = (@. .5*(vL+vR))
    wavg   = (@. .5*(wL+wR))

    unorm = (@. uL*uR + vL*vR + wL*wR)
    pa    = (@. rhoavg/(betaL+betaR))
    E_plus_p  = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. FxS1*vavg) # rho * u * v
    FxS4 = (@. FxS1*wavg) # rho * u * w
    FxS5 = (@. E_plus_p*uavg)

    FyS1 = (@. rholog*vavg)
    FyS2 = (@. FxS3) # rho * u * v
    FyS3 = (@. FyS1*vavg + pa)
    FyS4 = (@. FyS1*wavg) # rho * v * w
    FyS5 = (@. E_plus_p*vavg)

    FzS1 = (@. rholog*wavg)
    FzS2 = (@. FxS4) # rho * u * w
    FzS3 = (@. FyS4) # rho * v * w
    FzS4 = (@. FzS1*wavg + pa) # rho * w^2 + p
    FzS5 = (@. E_plus_p*wavg)

    Fx = (FxS1,FxS2,FxS3,FxS4,FxS5)
    Fy = (FyS1,FyS2,FyS3,FyS4,FyS5)
    Fz = (FzS1,FzS2,FzS3,FzS4,FzS5)
    return Fx,Fy,Fz
end

max_stars_repo_name: jlchan/ESDG-CNS

max_issues_repo_name: jlchan/ESDG-CNS

max_forks_repo_name: jlchan/ESDG-CNS


text:
module GLMModels

using SparseArrays, LinearAlgebra, Statistics

using ..RelayUtils

export Poisson, Binomial, Gaussian, Identity, Exponential, Logistic,
    nobs, nvar, activate!, get_objective, hessian, init_glm, predictors, response,
    bin_size

export Distribution, ActivationFunction, AbstractGLM, GLM, RegularizedGLM

const SPMat = SparseMatrixCSC{Float64,Int64}
const StorageVec = Union{Vector{Float64},Nothing}
const StorageMat = Union{Matrix{Float64},Nothing}
const ObjFlag = Union{Float64,Nothing}

abstract type Distribution end
struct Poisson <: Distribution end
struct Binomial <: Distribution end
struct Gaussian <: Distribution end

abstract type ActivationFunction end
struct Exponential <: ActivationFunction end
struct Logistic <: ActivationFunction end
struct Identity <: ActivationFunction end
# ============================================================================ #
abstract type AbstractGLM{D<:Distribution, A<:ActivationFunction} end
# ---------------------------------------------------------------------------- #
struct GLM{D,A} <: AbstractGLM{D,A}
    x::Matrix{Float64}
    y::Vector{Bool}
    r::Vector{Float64}
    dr::Vector{Float64}
    bin_size::Float64
end
# ---------------------------------------------------------------------------- #
struct RegularizedGLM{D,A,P<:AbstractMatrix} <: AbstractGLM{D,A}
    glm::GLM{D,A}
    prior::P
end
# ---------------------------------------------------------------------------- #
# main GLM constructor
function GLM(::Type{D}, ::Type{A}, x::Matrix{<:Real}, y::Vector{Bool}, bin_size::Real) where {D,A}
    return GLM{D,A}(x, y, zeros(Float64, length(y)), zeros(Float64, length(y)), bin_size)
end
# ---------------------------------------------------------------------------- #
# GLM copy constructor for easy partitioning
function GLM(glm::GLM{D,A}, x::Matrix{<:Real}, y::Vector{Bool}, bin_size::Real) where {D,A}
    return GLM(D, A, x, y, bin_size)
end
# ---------------------------------------------------------------------------- #
# construct a GLM from a RegularizedGLM, for use when training utilizes
# regularization but testing does not (which is always the case)
function GLM(rg::RegularizedGLM, x::Matrix{<:Real}, y::Vector{Bool}, bin_size::Real)
    return GLM(rg.glm, x, y, bin_size)
end
# ---------------------------------------------------------------------------- #
function RegularizedGLM(::Type{D}, ::Type{A}, x::Matrix{<:Real}, y::Vector{Bool}, bin_size::Real, pr::P) where {D,A,P}
    return RegularizedGLM{D,A,P}(GLM(D, A, x, y, bin_size), pr)
end
# ---------------------------------------------------------------------------- #
# alternate copy constructors for convienence use in GLMFit.cross_validate
function init_glm(glm::GLM{D,A}, x::Matrix{<:Real}, y::Vector{Bool}, bin_size::Real) where {D,A}
    return GLM(D, A, x, y, bin_size)
end
# ---------------------------------------------------------------------------- #
function init_glm(glm::RegularizedGLM{D,A}, x::Matrix{<:Real}, y::Vector{Bool}, bin_size::Real) where {D,A}
    return RegularizedGLM(D, A, x, y, bin_size, glm.prior)
end
# ============================================================================ #
@inline nobs(d::GLM) = length(d.y)
@inline nvar(d::GLM) = size(d.x, 2)
@inline nobs(d::RegularizedGLM) = length(d.glm.y)
@inline nvar(d::RegularizedGLM) = size(d.glm.x, 2)
@inline predictors(d::GLM) = d.x
@inline response(d::GLM) = d.y
@inline bin_size(d::GLM) = d.bin_size
@inline predictors(d::RegularizedGLM) = d.glm.x
@inline response(d::RegularizedGLM) = d.glm.y
@inline bin_size(d::RegularizedGLM) = d.glm.bin_size
# ============================================================================ #
get_objective(d::AbstractGLM) = (f::ObjFlag, g::StorageVec, p::Vector{Float64}) -> nlli!(d, f, p, g)

# function get_objective(d::GLM, pr::AbstractMatrix)
#     return (f::ObjFlag, g::StorageVec, p::Vector{Float64}) -> nlli_prior!(d, f, p, g, pr)
# end

get_objective(d::AbstractGLM{Binomial,Logistic}) =
    (f::ObjFlag, g::StorageVec, h::StorageMat, p::Vector{Float64}) -> nlli!(d, f, p, g, h)

# get_objective(d::GLM{Binomial,Logistic}, pr::AbstractMatrix) =
#     (f::ObjFlag, g::StorageVec, h::StorageMat, p::Vector{Float64}) -> nlli_prior!(d, f, p, g, h, pr)
# ============================================================================ #
function nlli!(d::RegularizedGLM{D,A,SPMat}, f::ObjFlag, p::Vector{Float64}, g::Vector{Float64}) where {D,A}
    nli = nlli!(d.glm, f, p, g)

    # is Optim requesting a gradient?
    if g != nothing
        g .+= d.prior * p
    end

    # Optim requested gradient only, don't bother calculating full objective
    f == nothing && return +Inf

    nli += .5 * (p' * d.prior * p)
    return nli
end
# ============================================================================ #
function nlli!(d::RegularizedGLM{D,A,SPMat}, f::ObjFlag, p::Vector{Float64}, g::StorageVec, h::StorageMat) where {D,A}
    nli = nlli!(d.glm, f, p, g, h)

    # is Optim requesting a gradient?
    if g != nothing
        g .+= d.prior * p
    end

    if h != nothing
        h .+= d.prior
    end

    # Optim requested gradient only, don't bother calculating full objective
    f == nothing && return +Inf

    nli += .5 * (p' * d.prior * p)
    return nli
end
# ============================================================================ #
function nlli!(d::RegularizedGLM{D,A,Matrix{Float64}}, f::ObjFlag, p::Vector{Float64}, g::Vector{Float64}) where {D,A}
    nli = nlli!(d.glm, f, p, g)

    # is Optim requesting a gradient?
    if g != nothing
        # g .+= cin * p
        LinearAlgebra.BLAS.gemv!('N', 1.0, d.prior, p, 1.0, g)
    end

    # Optim requested gradient only, don't bother calculating full objective
    f == nothing && return +Inf

    nli += .5 * (p' * d.prior * p)
    return nli
end
# ============================================================================ #
@inline function predict(d::RegularizedGLM, p::Vector{Float64})
    return predict(d.glm, p)
end
# ---------------------------------------------------------------------------- #
@inline function predict(d::GLM{D,L}, p::Vector{Float64}) where {D<:Distribution,L<:ActivationFunction}
    mul!(d.r, d.x, p)
    return activate!(L, d.r, d.dr, d.bin_size)
end
# ---------------------------------------------------------------------------- #
@inline function predict(d::GLM{D,Identity}, p::Vector) where D<:Distribution
    mul!(d.r, d.x, p)
    return d.r, d.dr
end
# ============================================================================ #
function nlli!(d::GLM{Gaussian,L}, p::Vector, g::StorageVec) where L <: ActivationFunction
    error("Not implemented!")
end
# ============================================================================ #
# source: https://www.comp.nus.edu.sg/~cs5240/lecture/matrix-diff.pdf slide 32 - 33
function nlli!(d::GLM{Gaussian,Identity}, f::ObjFlag, p::Vector{Float64}, g::StorageVec)
    r, _ = predict(d, p)

    # is Optim requesting a gradient?
    if g != nothing
        # Gaussian-Identity maximum likelyhood gradient: 2X'Xp - 2X'y
        # our approach is to rearrange the gradient as:
        # 2X'Xp - 2X'y = 2X'(Xp - y)

        # compute the term in the parenthese
        # d.r = Xp - y
        d.r .= d.y
        LinearAlgebra.BLAS.gemv!('N', 1.0, d.x, p, -1.0, d.r)

        # update the output variable <g> as
        # g = 2X'd.r
        LinearAlgebra.BLAS.gemv!('T', 2.0, d.x, d.r, 0.0, g)
    end

    # Optim requested gradient only, don't bother calculating full objective
    f == nothing && return +Inf

    # -log(exp(-(y .- r).^2))
    return sum(abs2, d.y .- r)
end
# ============================================================================ #
null_nlli(d::GLM{Gaussian,Identity}) = sum(abs2, d.y .- mean(d.y))
# ============================================================================ #
@inline function dnlli!(d::GLM{Poisson,Exponential}, r::Vector{Float64},
    dr::Vector{Float64}, g::Vector{Float64})

    # b/c exponential is it's own derivitive, r is likely to be aliased with
    # dr, thus we can simplify the gradient calculation that is used for
    # Poisson with any other activation function (see next impl. of dnlli!() below)
    # (i.e. dr / r == 1)

    # less efficent procedure that we implement below
    # g .= (d.x' * r) .- (d.x' * d.y)

    # second term g .= (d.x' * d.y)
    mul!(g, d.x', d.y)

    # first term and update (g .= d.x' * r -. g)
    LinearAlgebra.BLAS.gemv!('T', 1.0, d.x, r, -1.0, g)

    return g
end
# ============================================================================ #
@inline function dnlli!(d::GLM{Poisson,L}, r::Vector{Float64},
    dr::Vector{Float64}, g::Vector{Float64}) where L <: ActivationFunction

    # less efficent procedure that we implement below
    # g .= (d.x' * dr) .- (d.x' * ((dr ./ r) .* d.y))

    # second term g .= (d.x' * ((dr ./ r) .* d.y))
    d.dr .= ((dr ./ r) .* d.y) # avoid allocating a temporary for element-wise part
    mul!(g, d.x', d.dr)

    # first term and update (g .= d.x' * r -. g)
    LinearAlgebra.BLAS.gemv!('T', 1.0, d.x, r, -1.0, g)

    return g
end
# ============================================================================ #
function nlli!(d::GLM{Poisson,L}, f::ObjFlag, p::Vector{Float64}, g::StorageVec) where L <: ActivationFunction
    r, dr = predict(d, p)

    # is Optim requesting a gradient?
    if g != nothing
        dnlli!(d, r, dr, g)
    end

    # Optim requested gradient only, don't bother calculating full objective
    f == nothing && return +Inf

    return -dot(log.(r), d.y) + sum(r)
end
# ============================================================================ #
function null_nlli(d::GLM{Poisson,L}) where L <: ActivationFunction
    n = sum(d.y)
    return -n * log(n / length(d.y)) + n
end
# ============================================================================ #
@inline function dnlli!(d::GLM{Binomial,Logistic}, r::Vector{Float64},
    dr::Vector{Float64}, g::Vector{Float64}) where L <: ActivationFunction

    # as for Poisson GLM, for Binomial-Logistic we can simplify the gradient
    # compared to the case where we have an arbitrary, unknown activation
    # function (see next implmentation of dnlli!() below)

    # avoid temporary allocation
    d.dr .= d.y .- r

    # -X' (y - r)
    LinearAlgebra.BLAS.gemv!('T', -1.0, d.x, d.dr, 0.0, g)
    return g
end
# ---------------------------------------------------------------------------- #
@inline function dnlli!(d::GLM{Binomial,L}, r::Vector{Float64},
    dr::Vector{Float64}, g::Vector{Float64}) where L <: ActivationFunction

    # less efficent procedure that we implement below
    # g .= d.x' * (.-((r .- d.y) ./ (r.^2 .- r)) .* dr)

    # avoid allocating a temporary for element-wise part
    d.dr .= .-((r .- d.y) ./ (r.^2 .- r)) .* dr
    mul!(g, d.x', d.dr)

    return g
end
# ---------------------------------------------------------------------------- #
# sources:
# * https://web.stanford.edu/class/archive/cs/cs109/cs109.1178/lectureHandouts/220-logistic-regression.pdf
# * http://www.stat.cmu.edu/~cshalizi/uADA/12/lectures/ch12.pdf
function nlli!(d::GLM{Binomial,L}, f::ObjFlag, p::Vector{Float64}, g::StorageVec) where L <: ActivationFunction
    r, dr = predict(d, p)

    # is Optim requesting a gradient?
    if g != nothing
        dnlli!(d, r, dr, g)
    end

    # Optim requested gradient only, don't bother calculating full objective
    f == nothing && return +Inf

    return -RelayUtils.binomial_lli_turbo(r, d.y)
end
# ---------------------------------------------------------------------------- #
function nlli!(d::GLM{Binomial,Logistic}, f::ObjFlag, p::Vector{Float64}, g::StorageVec, h::StorageMat)
    r, dr = predict(d, p)

    # is Optim requesting a gradient?
    if g != nothing
        dnlli!(d, r, dr, g)
    end

    # is Optim requesting a hessian?
    if h != nothing
        # h .= d.x' * (d.x .* (r .* (1.0 .- r)))
        d.dr .= (r .* (1.0 .- r))
        # h .= d.x' * (d.x .* d.dr)
        LinearAlgebra.BLAS.gemm!('T', 'N', 1.0, d.x, d.x .* d.dr, 0.0, h)
    end

    # Optim requested gradient only, don't bother calculating full objective
    f == nothing && return +Inf

    return -RelayUtils.binomial_lli_turbo(r, d.y)
end
# ============================================================================ #
function null_nlli(d::GLM{Binomial,L}) where L <: ActivationFunction
    return -RelayUtils.binomial_lli(d.y)

    # if we normalize by the number of spikes (length(d.y)) we get:
    # ef = sum(d.y) / length(y)
    # return -(ef * log(ef) + (1.0 - ef) * log(1.0 - ef))
end
# ============================================================================ #
@inline function activate!(::Type{Identity}, xp::Vector{Float64}, dr::Vector{Float64}, bin_size::Real)
    dr .= 1.0
    return xp, dr
end
# ============================================================================ #
function activate!(::Type{Exponential}, xp::Vector{Float64}, dr::Vector{Float64}, bin_size::Real)
    # exponential is its own derivitive
    xp .= exp.(xp) .* bin_size
    return xp, xp
end
# ============================================================================ #
# see also: https://web.stanford.edu/class/archive/cs/cs109/cs109.1178/lectureHandouts/220-logistic-regression.pdf slide 5
function activate!(::Type{Logistic}, xp::Vector{Float64}, dr::Vector{Float64}, bin_size::Real)
    # # logistic function, f(xp)
    # xp .= 1.0 ./ (1.0 .+ exp.(.-xp))
    #
    # # logistic derivitive, f'(xp))
    # dr .= xp .* (1.0 .- xp)

    return RelayUtils.logistic_derivative_turbo!(xp, dr)
end
# ============================================================================ #
function hessian(d::GLM{Poisson,Exponential}, p::Vector)
    # for computing the stderr of parameter estimates

    r, _ = predict(d, p)

    # NOTE: see also GLMspiketools/glmtools_fitting/Loss_GLM_logli.m
    # exponential is it's own derivitive (and 2nd), so the "weights" for all
    # "spiking terms" are 0 (line 101), thus the whole complicated mess just
    # becomes the following (seems crazy that it simplifies this much, but I've
    # checked the result agains Pillow's original code - for 1 example pair -
    # and it comes out correct)

    return d.x' * (d.x .* r)
end
# ============================================================================ #
function hessian(d::GLM{Binomial,Logistic}, p::Vector)
    # for computing the stderr of parameter estimates
    r, _ = predict(d, p)

    #= NOTE: see also
        https://stats.stackexchange.com/questions/68391/hessian-of-logistic-function
        http://personal.psu.edu/jol2/course/stat597e/notes2/logit.pdf slide 19
        https://gist.github.com/faisal-w/9acdeefbddf2bbc6eb72
    =#

    return d.x' * (d.x .* (r .* (1.0 .- r)))
end
# ============================================================================ #
end

# nlli -(y log(r) + (1 - y)log(1 - r))
# d/dr = -((r - y) / (r^2 - r))


text:
function mul!(C::T, A::T, B::T, idx_A::String, idx_B::String, idx_C::String; nt=C_NULL) where T <: TTensor{T2} where T2 <: AbstractFloat
    tblis_tmul = dlsym(tblis,:tblis_tensor_mult)
    ccall(tblis_tmul,Cvoid,(Ptr{Nothing},Ptr{Nothing},
                            Ref{TTensor{T2}},Cstring,
                            Ref{TTensor{T2}},Cstring,
                            Ref{TTensor{T2}},Cstring),

          C_NULL,C_NULL,
          A,idx_A,
          B,idx_B,
          C,idx_C)
end

function add!(A::T, B::T, idx_A::String, idx_B::String) where T <: TTensor{T2} where T2 <: AbstractFloat
    tblis_tadd = dlsym(tblis,:tblis_tensor_add)
    ccall(tblis_tadd,Cvoid,(Ptr{Nothing},Ptr{Nothing},
                            Ref{TTensor{T2}},Cstring,
                            Ref{TTensor{T2}},Cstring),
          C_NULL,C_NULL,
          A,idx_A,
          B,idx_B)
end


text:
# This file includes the plugin for calculation of Lyapunov exponents.

using ChaosTools

"""
    Lyapunov(;m::Int=15, J::Int=5, ni::Int=300, ts::Float64=0.01)

Constructs a `Lyapunov` plugin. The [`process(plg::Lyapunov, x)`](@ref) function computes the maximum numerical Lyapunov exponents. `m` is the reconstruction dimension, `J` is the amount of delay in reconstruction, `ni` is the number of steps during transient steps and `ts` is the sampling time between samples of the input data `x`. See also: (https://juliadynamics.github.io/DynamicalSystems.jl/latest/chaos/nlts/)
"""
struct Lyapunov <: AbstractPlugin
    m::Int
    J::Int
    ni::Int
    ts::Float64
end
Lyapunov(;m::Int=15, J::Int=5, ni::Int=300, ts::Float64=0.01) = Lyapunov(m, J, ni, ts)

show(io::IO, plg::Lyapunov) = print(io,
    "Lyapunov(embeddingdim:$(plg.m), numlags:$(plg.J), numiteration:$(plg.ni), samplingtime:$(plg.ts)")

"""
    process(plg::Lyapunov, x)

Computes the maximum Lyapunov exponent of the input data `x`.

# Example
```julia
julia> using Random

julia> rng = MersenneTwister(1234);

julia> x = rand(rng, 1000);

julia> plg = Plugins.Lyapunov()
Lyapunov(embeddingdim:15, numlags:5, numiteration:300, samplingtime:0.01

julia> process(plg, x)
-0.42032928176193973
```
"""
function process(plg::Lyapunov, x)
    ndims(x) == 1 || (x = x[:])
    ntype = FixedMassNeighborhood(5)
    ks = 1 : 4 : plg.ni
    R = reconstruct(x, plg.m, plg.J)
    E = numericallyapunov(R, ks, ntype=ntype)
    val = linear_region(plg.ts .* ks, E)[2]
    return val
end

# using NearestNeighbors
# using LinearAlgebra
# using Statistics
# using LsqFit
# import Base.log
# process(plg::Lyapunov, x) = lyapunov(x, plg.m, plg.J, plg.ni, plg.ts)[2]

# function reconstruct(x, m, J)
#     N = length(x)
#     M = N - (m - 1) * J
#     X = zeros(m, M)
#     for i = 1 : M
#         data =  x[i : J : i + (m - 1)* J]
#         X[:, i] = data
#     end
#     X
# end

# function knneighbours(X)
#     tree = KDTree(X)
#     [knn(tree, X[:, j], 2)[1][1] for j = 1 : size(X, 2)]
# end

# distance(Xj, Xjhat) = norm(Xj - Xjhat)

# Base.log(a::Vector) = log.(a)

# function lyapunov(x, m, J, ni, ts)
#     X = reconstruct(x, m, J)
#     M = size(X, 2)
#     js = collect(1:M)
#     jbars = knneighbours(X)
#     m1 = js .+ ni .<= M
#     m2 = jbars .+ ni .<= M
#     m = m1 .& m2
#     mjs = js[m]
#     mjbars = jbars[m]
#     y = mean(log.([distance.(eachcol(X[:, j : j + ni]), eachcol(X[:, jbar : jbar + ni]))
#         for (j, jbar) in zip(mjs, mjbars)])) / ts
#     @. model(x, p) = p[1]*x
#     ydata = y[round(Int, length(y) * 0.25) : end]   # Discard first transient region and go to the linear region.
#     lambda = coef(curve_fit(model, collect(1:length(ydata)), ydata, rand(1)))[1]
#     return y, lambda
# end


text:
export sgbonds

function solve(::Type{TT}, sg::AbstractSpinglass{LT}; usecuda=false) where {LT<:SquareLattice,TT}
    # Yao gates
    lt = sg.lattice
    Lx, Ly = lt.Nx, lt.Ny
    reg = _init_reg(TT, lt.Ny, Val(usecuda))
    k = 0

    for i=1:Lx
        println("Layer $i/$Lx")
        i!=1 && for j=1:Ly
            k += 1
            reg |> put(Ly, j=>Ghb(bondtensor(TT, sg, k)))
        end
        for j=1:Ly
            reg |> put(Ly, j=>Gh(vertextensor(TT, sg, (i-1)*Ly + j)))
        end
        for j=1:Ly-1
            k += 1
            reg |> put(Ly, (j,j+1)=>Gvb(bondtensor(TT, sg, k)))
        end
    end
    sum(state(reg))
end

function solve(lattice::Viznet.AbstractLattice, Js::Vector{T}, hs::Vector{T}; usecuda=false) where {T}
    sg = Spinglass(lattice, Js, hs)
    solve(sg; usecuda=usecuda)
end

# to index `h`.
function sgvertices(lt::SquareLattice)
    v = Int[]
    LI = LinearIndices(lt)
    for i=1:lt.Nx
        for j=1:lt.Ny
            push!(v, LI[i,j])
        end
    end
    return v
end

# to index `J`.
function sgbonds(lt::SquareLattice)
    edges  =Tuple{Int,Int}[]
    append!(edges, v_bonds(lt, 1))
    for i=2:lt.Nx
        append!(edges, h_bonds(lt, i-1))
        append!(edges, v_bonds(lt, i))
    end
    edges
end

function v_bonds(lt::SquareLattice, i::Int)
    LI = LinearIndices(size(lt))
    map(j -> (LI[i,j], LI[i,j+1]), 1:lt.Ny-1)
end

function h_bonds(lt::SquareLattice, i::Int)
    LI = LinearIndices(size(lt))
    map(j -> (LI[i,j], LI[i+1,j]), 1:lt.Ny)
end

regsize(lt::SquareLattice) = lt.Ny

function _init_reg(::Type{T}, lt::SquareLattice, usecuda) where T<:TropicalTypes
    _init_reg(T, lt.Ny, usecuda)
end


text:
# Types

export
    Interaction,
    GeneralInteraction,
    SpecificInteraction,
    AbstractCutoff,
    Simulator,
    Thermostat,
    NeighborFinder,
    Logger,
    Atom,
    AtomData,
    Simulation

const DefaultFloat = Float64

"An interaction between atoms that contributes to forces on the atoms."
abstract type Interaction end

"""
A general interaction that will apply to all or most atom pairs.
Custom general interactions should sub-type this type.
"""
abstract type GeneralInteraction <: Interaction end

"""
A specific interaction between sets of specific atoms, e.g. a bond angle.
Custom specific interactions should sub-type this type.
"""
abstract type SpecificInteraction <: Interaction end

"""
A general type of cutoff encoding the approximation used for a potential.
Interactions can be parameterized by the cutoff behavior.
"""
abstract type AbstractCutoff end

"""
A type of simulation to run, e.g. leap-frog integration or energy minimisation.
Custom simulators should sub-type this type.
"""
abstract type Simulator end

"""
A way to keep the temperature of a simulation constant.
Custom thermostats should sub-type this type.
"""
abstract type Thermostat end

"""
A way to find near atoms to save on simulation time.
Custom neighbor finders should sub-type this type.
"""
abstract type NeighborFinder end

"""
A way to record a property, e.g. the temperature, throughout a simulation.
Custom loggers should sub-type this type.
"""
abstract type Logger end

"""
    Atom(; <keyword arguments>)

An atom and its associated information.
Properties unused in the simulation or in analysis can be left with their
default values.
The types used should be bits types if the GPU is going to be used.

# Arguments
- `index::Int`: the index of the atom in the system.
- `charge::C=0.0u"q"`: the charge of the atom, used for electrostatic interactions.
- `mass::M=0.0u"u"`: the mass of the atom.
- `σ::S=0.0u"nm"`: the Lennard-Jones finite distance at which the inter-particle
    potential is zero.
- `ϵ::E=0.0u"kJ * mol^-1"`: the Lennard-Jones depth of the potential well.
"""
struct Atom{C, M, S, E}
    index::Int
    charge::C
    mass::M
    σ::S
    ϵ::E
end

function Atom(;
                index=1,
                charge=0.0u"q",
                mass=0.0u"u",
                σ=0.0u"nm",
                ϵ=0.0u"kJ * mol^-1")
    return Atom(index, charge, mass, σ, ϵ)
end

function Base.show(io::IO, a::Atom)
    print(io, "Atom with index ", a.index, ", charge=", a.charge,
            ", mass=", a.mass, ", σ=", a.σ, ", ϵ=", a.ϵ)
end

"""
    AtomData(atom_type, atom_name, res_number, res_name)

Data associated with an atom.
Storing this separately allows the atom types to be bits types and hence
work on the GPU.
"""
struct AtomData
    atom_type::String
    atom_name::String
    res_number::Int
    res_name::String
    element::String
end

function AtomData(;
                    atom_type="?",
                    atom_name="?",
                    res_number=1,
                    res_name="???",
                    element="?")
    return AtomData(atom_type, atom_name, res_number, res_name, element)
end

"""
    Simulation(; <keyword arguments>)

The data needed to define and run a molecular simulation.
Properties unused in the simulation or in analysis can be left with their
default values.

# Arguments
- `simulator::Simulator`: the type of simulation to run.
- `atoms::A`: the atoms, or atom equivalents, in the simulation. Can be
    of any type.
- `atoms_data::AD`: data associated with the atoms, allowing the atoms to be
    bits types and hence work on the GPU.
- `specific_inter_lists::SI=()`: the specific interactions in the simulation,
    i.e. interactions between specific atoms such as bonds or angles. Typically
    a `Tuple`.
- `general_inters::GI=()`: the general interactions in the simulation, i.e.
    interactions between all or most atoms such as electrostatics. Typically a
    `Tuple`.
- `coords::C`: the coordinates of the atoms in the simulation. Typically a
    `Vector` of `SVector`s of any dimension and type `T`, where `T` is an
    `AbstractFloat` type.
- `velocities::V=zero(coords)`: the velocities of the atoms in the simulation,
    which should be the same type as the coordinates. The meaning of the
    velocities depends on the simulator used, e.g. for the `VelocityFreeVerlet`
    simulator they represent the previous step coordinates for the first step.
- `temperature::T=0.0`: the temperature of the simulation.
- `box_size::B`: the size of the box in which the simulation takes place.
- `neighbor_finder::NeighborFinder=NoNeighborFinder()`: the neighbor finder
    used to find close atoms and save on computation.
- `thermostat::Thermostat=NoThermostat()`: the thermostat which applies during
    the simulation.
- `loggers::Dict{String, <:Logger}=Dict()`: the loggers that record properties
    of interest during the simulation.
- `timestep::S=0.0`: the timestep of the simulation.
- `n_steps::Integer=0`: the number of steps in the simulation.
- `force_unit::F=u"kJ * mol^-1 * nm^-1"`: the unit of force used in the
    simulation.
- `energy_unit::E=u"kJ * mol^-1"`: the unit of energy used in the simulation.
- `gpu_diff_safe::Bool`: whether to use the GPU implementation. Defaults to
    `isa(coords, CuArray)`.
"""
struct Simulation{D, T, A, AD, C, V, GI, SI, B, S, F, E}
    simulator::Simulator
    atoms::A
    atoms_data::AD
    specific_inter_lists::SI
    general_inters::GI
    coords::C
    velocities::V
    temperature::T
    box_size::B
    neighbors::Vector{Tuple{Int, Int, Bool}}
    neighbor_finder::NeighborFinder
    thermostat::Thermostat
    loggers::Dict{String, <:Logger}
    timestep::S
    n_steps::Int
    n_steps_made::Vector{Int}
    force_unit::F
    energy_unit::E
end

Simulation{D}(args...) where {D, T, A, AD, C, V, GI, SI, B, S, F, E} = Simulation{D, T, A, AD, C, V, GI, SI, B, S, F, E}(args...)

function Simulation(;
                    simulator=VelocityVerlet(),
                    atoms,
                    atoms_data=[],
                    specific_inter_lists=(),
                    general_inters=(),
                    coords,
                    velocities=zero(coords),
                    temperature=0.0u"K",
                    box_size,
                    neighbors=Tuple{Int, Int, Bool}[],
                    neighbor_finder=NoNeighborFinder(),
                    thermostat=NoThermostat(),
                    loggers=Dict{String, Logger}(),
                    timestep=0.0u"ps",
                    n_steps=0,
                    n_steps_made=[0],
                    force_unit=u"kJ * mol^-1 * nm^-1",
                    energy_unit=u"kJ * mol^-1",
                    gpu_diff_safe=isa(coords, CuArray))
    T = typeof(temperature)
    A = typeof(atoms)
    AD = typeof(atoms_data)
    C = typeof(coords)
    V = typeof(velocities)
    GI = typeof(general_inters)
    SI = typeof(specific_inter_lists)
    B = typeof(box_size)
    S = typeof(timestep)
    F = typeof(force_unit)
    E = typeof(energy_unit)
    return Simulation{gpu_diff_safe, T, A, AD, C, V, GI, SI, B, S, F, E}(
                simulator, atoms, atoms_data, specific_inter_lists, general_inters,
                coords, velocities, temperature, box_size, neighbors, neighbor_finder,
                thermostat, loggers, timestep, n_steps, n_steps_made, force_unit,
                energy_unit)
end

function Base.show(io::IO, s::Simulation)
    print(io, "Simulation with ", length(s.coords), " atoms, ",
                typeof(s.simulator), " simulator, ", s.timestep, " timestep, ",
                s.n_steps, " steps, ", first(s.n_steps_made), " steps made")
end