text:
using Yao, YaoFlux, Zygote, Flux.Optimise

t = ArrayReg(bit"0") + ArrayReg(bit"1")
normalize!(t)

function fid(xs)
    r = zero_state(1)
    U = mat(chain(Rx(xs[1]), Rz(xs[2])))
    return abs(statevec(t)' * U * statevec(r))
end

function train()
    opt = ADAM()
    xs = rand(2)
    for _ in 1:1000
        println(fid(xs))
        Optimise.update!(opt, xs, -fid'(xs))
    end
    return xs
end

# train()

using Makie, GLMakie, Observables, Colors

AbstractPlotting.inline!(true)

function bloch_arrow(θ, ϕ)
    x = sin(2θ) * cos(ϕ)
    y = sin(2θ) * sin(ϕ)
    z = cos(2θ)
    return (x, y, z)
end

function bloch_arrow(r::ArrayReg)
    @assert nqubits(r) == 1 "invalid quantum state, expect only one qubit"
    st = statevec(r)
    global_phase = exp(im * angle(st[1]))
    st = st ./ global_phase
    θ = acos(real(st[1])) # st[1] is real
    ϕ = iszero(θ) ? zero(θ) : angle(st[2] / sin(θ))
    return bloch_arrow(θ, ϕ)
end

function plot_arrow!(scene, r::ArrayReg; kwargs...)
    x, y, z = bloch_arrow(r)
    lines!(scene, [0.0, x], [0.0, y], [0.0, z]; kwargs...)
    return scene
end

function plot_block!(scene, x::AbstractBlock; linewidth=5.0, kwargs...)
    r = ArrayReg(bit"0")
    plot_arrow!(scene, r; color=:red, linewidth=linewidth, kwargs...)
    rand_color = RGB(rand(0:255)/255, rand(0:255)/255, rand(0:255)/255)
    plot_arrow!(scene, apply!(r, x); color=rand_color ,linewidth=linewidth, kwargs...)
    return scene
end

function bloch_sphere!(scene)
    wireframe!(scene, Sphere(Point3f0(0), 1f0), color=:cyan3, linewidth=0.1, thickness = 0.6f0, transparency = true)
    return scene
end

function plotapply(x::AbstractBlock)
    AbstractPlotting.inline!(false)
    scene = Scene(resolution=(1000, 1000), center=false)
    bloch_sphere!(scene)
    plot_block!(scene, x)
    rotate_cam!(scene, 0.5, 0.0, 0.0)
    return scene
end

function animate_rot()
    scene = Scene(resolution=(1000, 1000))
    bloch_sphere!(scene)
    scene.center = false

    time = Node(0.0)
    scene = lift()
    bloch_arrow!(scene, x)
    rotate_cam!(scene, 0.5, 0.0, 0.0)
end

# Pauli Gates
function animate_rot(blk::AbstractBlock)
    AbstractPlotting.inline!(false)
    screen = GLMakie.Screen()
    GLMakie.GLFW.set_visibility!(GLMakie.to_native(screen), AbstractPlotting.use_display[])

    AbstractPlotting.inline!(true)
    scene = Scene(resolution=(1000, 1000), center=false)

    bloch_sphere!(scene)

    r = ArrayReg(bit"0")
    plot_arrow!(scene, r; color=:red, linewidth=5.0)
    time = Node(0.0)

    lifted = lift(time) do t
        x, y, z = bloch_arrow(apply!(r, rot(blk, t)))
        [0.0, x], [0.0, y], [0.0, z]
    end

    lines!(scene, lifted; color=:blue, linewidth=5.0)
    rotate_cam!(scene, 0.1, 0.0, 0.0)

    for _ in 1:1000
        push!(time, to_value(time) + 0.001)
        display(screen, scene)
    end
    return scene
end

function animate_train()
    AbstractPlotting.inline!(false)
    screen = GLMakie.Screen()
    GLMakie.GLFW.set_visibility!(GLMakie.to_native(screen), AbstractPlotting.use_display[])

    AbstractPlotting.inline!(true)
    scene = Scene(resolution=(1000, 1000), center=false)

    bloch_sphere!(scene)

    plot_arrow!(scene, zero_state(1); color=:red, linewidth=5.0)
    GHZ = normalize!(ArrayReg(bit"0") + ArrayReg(bit"1"))
    plot_arrow!(scene, GHZ; color=:green, linewidth=5.0)

    θ, ϕ = Node(0.0), Node(0.0)

    lifted = lift(θ, ϕ) do θ, ϕ
        x, y, z = bloch_arrow(apply!(zero_state(1), chain(Rx(θ), Rz(ϕ))))
        [0.0, x], [0.0, y], [0.0, z]
    end

    lines!(scene, lifted; color=:blue, linewidth=5.0)
    # scatter!(scene, lifted; color=:black, linewidth=5.0)
    rotate_cam!(scene, -1.1, 0.0, 0.0)

    opt = ADAM()
    xs = [0.01, 0.01]
    for k in 1:1000
        # sleep(0.001)
        println(fid(xs))
        Optimise.update!(opt, xs, -fid'(xs))

        r = apply!(zero_state(1), chain(Rx(xs[1]), Rz(xs[2])))
        x, y, z = bloch_arrow(r)

        if k % 10 == 0
            scatter!(scene, [x], [y], [z]; color=:blue, markersize=0.02)
        end
        push!(θ, xs[1])
        push!(ϕ, xs[2])
        display(screen, scene)
    end
    return scene
end

function make_gif()
    AbstractPlotting.inline!(false)
    screen = GLMakie.Screen()
    GLMakie.GLFW.set_visibility!(GLMakie.to_native(screen), AbstractPlotting.use_display[])

    AbstractPlotting.inline!(true)
    scene = Scene(resolution=(1000, 1000), center=false)

    bloch_sphere!(scene)

    plot_arrow!(scene, zero_state(1); color=:red, linewidth=5.0)
    GHZ = normalize!(ArrayReg(bit"0") + ArrayReg(bit"1"))
    plot_arrow!(scene, GHZ; color=:green, linewidth=5.0)

    θ, ϕ = Node(0.0), Node(0.0)

    lifted = lift(θ, ϕ) do θ, ϕ
        x, y, z = bloch_arrow(apply!(zero_state(1), chain(Rx(θ), Rz(ϕ))))
        [0.0, x], [0.0, y], [0.0, z]
    end

    lines!(scene, lifted; color=:blue, linewidth=5.0)
    # scatter!(scene, lifted; color=:black, linewidth=5.0)
    rotate_cam!(scene, -1.1, 0.0, 0.0)

    opt = ADAM()
    xs = [0.01, 0.01]

    record(scene, "bloch.gif", 1:1000) do i
        Optimise.update!(opt, xs, -fid'(xs))

        r = apply!(zero_state(1), chain(Rx(xs[1]), Rz(xs[2])))
        x, y, z = bloch_arrow(r)

        if k % 10 == 0
            scatter!(scene, [x], [y], [z]; color=:blue, markersize=0.02)
        end
        push!(θ, xs[1])
        push!(ϕ, xs[2])
    end
end

s = animate_train()
make_gif()

ext: jl


text:
using JUGRNModelGenerator

# setup paths -
path_to_input_file = "/Users/jeffreyvarner/Desktop/julia_work/JUGRNModelGenerator.jl/test/data/Venus.json"
path_to_output_dir = "/Users/jeffreyvarner/Desktop/julia_work/JUGRNModelGenerator.jl/test/generated_code"

# build the model object -
julia_model_object = build_julia_model_object(path_to_input_file, path_to_output_dir) |> check

# try to generate code -
result = make_julia_model(julia_model_object)


text:
# Utilities to support the testing of distributions and samplers

import Test: @test

# auxiliary functions

# to workaround issues of Base.linspace
function _linspace(a::Float64, b::Float64, n::Int)
    intv = (b - a) / (n - 1)
    r = Vector{Float64}(undef, n)
    @inbounds for i = 1:n
        r[i] = a + (i-1) * intv
    end
    r[n] = b
    return r
end

#################################################
#
#   Driver functions
#
#################################################

# testing the implementation of a discrete univariate distribution
#
function test_distr(distr::DiscreteUnivariateDistribution, n::Int; testquan::Bool=true)

    test_range(distr)
    vs = get_evalsamples(distr, 0.00001)

    test_support(distr, vs)
    test_evaluation(distr, vs, testquan)
    test_range_evaluation(distr)

    test_stats(distr, vs)
    test_samples(distr, n)
    test_params(distr)
end

# testing the implementation of a continuous univariate distribution
#
function test_distr(distr::ContinuousUnivariateDistribution, n::Int; testquan::Bool=true)

    test_range(distr)
    vs = get_evalsamples(distr, 0.01, 2000)

    test_support(distr, vs)
    test_evaluation(distr, vs, testquan)

    xs = test_samples(distr, n)
    allow_test_stats(distr) && test_stats(distr, xs)
    test_params(distr)
end

#################################################
#
#   Core testing functions
#
#################################################

#### Testing sampleable objects (samplers)

# for discrete samplers
#
function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable instance
                      distr::DiscreteUnivariateDistribution,    # corresponding distribution
                      n::Int;                                   # number of samples to generate
                      q::Float64=1.0e-7,                        # confidence interval, 1 - q as confidence
                      verbose::Bool=false)                      # show intermediate info (for debugging)

    # The basic idea
    # ------------------
    #   Generate n samples, and count the occurences of each value within a reasonable range.
    #   For each distinct value, it computes an confidence interval of the counts
    #   and checks whether the count is within this interval.
    #
    #   If the distribution has a bounded range, it also checks whether
    #   the samples are all within this range.
    #
    #   By setting a small q, we ensure that failure of the tests rarely
    #   happen in practice.
    #

    verbose && println("test_samples on $(typeof(s))")

    n > 1 || error("The number of samples must be greater than 1.")
    0.0 < q < 0.1 || error("The value of q must be within the open interval (0.0, 0.1).")

    # determine the range of values to examine
    vmin = minimum(distr)
    vmax = maximum(distr)

    rmin = floor(Int,quantile(distr, 0.00001))::Int
    rmax = floor(Int,quantile(distr, 0.99999))::Int
    m = rmax - rmin + 1  # length of the range
    p0 = pdf.(Ref(distr), rmin:rmax)  # reference probability masses
    @assert length(p0) == m

    # determine confidence intervals for counts:
    # with probability q, the count will be out of this interval.
    #
    clb = Vector{Int}(undef, m)
    cub = Vector{Int}(undef, m)
    for i = 1:m
        bp = Binomial(n, p0[i])
        clb[i] = floor(Int,quantile(bp, q/2))
        cub[i] = ceil(Int,cquantile(bp, q/2))
        @assert cub[i] >= clb[i]
    end

    # generate samples
    samples = rand(s, n)
    @assert length(samples) == n

    # scan samples and get counts
    cnts = zeros(Int, m)
    for i = 1:n
        @inbounds si = samples[i]
        if rmin <= si <= rmax
            cnts[si - rmin + 1] += 1
        else
            vmin <= si <= vmax ||
                error("Sample value out of valid range.")
        end
    end

    # check the counts
    for i = 1:m
        verbose && println("v = $(rmin+i-1) ==> ($(clb[i]), $(cub[i])): $(cnts[i])")
        clb[i] <= cnts[i] <= cub[i] ||
            error("The counts are out of the confidence interval.")
    end
    return samples
end

test_samples(distr::DiscreteUnivariateDistribution, n::Int; q::Float64=1.0e-6, verbose::Bool=false) =
    test_samples(distr, distr, n; q=q, verbose=verbose)

# for continuous samplers
#
function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable instance
                      distr::ContinuousUnivariateDistribution,  # corresponding distribution
                      n::Int;                                   # number of samples to generate
                      nbins::Int=50,                            # divide the main interval into nbins
                      q::Float64=1.0e-6,                        # confidence interval, 1 - q as confidence
                      verbose::Bool=false)                      # show intermediate info (for debugging)

    # The basic idea
    # ------------------
    #   Generate n samples, divide the interval [q00, q99] into nbins bins, where
    #   q01 and q99 are 0.01 and 0.99 quantiles, and count the numbers of samples
    #   falling into each bin. For each bin, we will compute a confidence interval
    #   of the number, and see whether the actual number is in this interval.
    #
    #   If the distribution has a bounded range, it also checks whether
    #   the samples are all within this range.
    #
    #   By setting a small q, we ensure that failure of the tests rarely
    #   happen in practice.
    #

    verbose && println("test_samples on $(typeof(s))")

    n > 1 || error("The number of samples must be greater than 1.")
    nbins > 1 || error("The number of bins must be greater than 1.")
    0.0 < q < 0.1 || error("The value of q must be within the open interval (0.0, 0.1).")

    # determine the range of values to examine
    vmin = minimum(distr)
    vmax = maximum(distr)

    local rmin::Float64
    local rmax::Float64
    if applicable(quantile, distr, 0.5)
        rmin = quantile(distr, 0.01)
        rmax = quantile(distr, 0.99)
    elseif isfinite(vmin) && isfinite(vmax)
        rmin = vmin
        rmax = vmax
    end
    edges = _linspace(rmin, rmax, nbins + 1)

    # determine confidence intervals for counts:
    # with probability q, the count will be out of this interval.
    #
    clb = Vector{Int}(undef, nbins)
    cub = Vector{Int}(undef, nbins)
    cdfs = cdf.(Ref(distr), edges)

    for i = 1:nbins
        pi = cdfs[i+1] - cdfs[i]
        bp = Binomial(n, pi)
        clb[i] = floor(Int,quantile(bp, q/2))
        cub[i] = ceil(Int,cquantile(bp, q/2))
        @assert cub[i] >= clb[i]
    end

    # generate samples
    samples = rand(s, n)
    @assert length(samples) == n

    # check whether all samples are in the valid range
    for i = 1:n
        @inbounds si = samples[i]
        vmin <= si <= vmax ||
            error("Sample value out of valid range.")
    end

    # get counts
    cnts = fit(Histogram, samples, edges; closed=:right).weights
    @assert length(cnts) == nbins

    # check the counts
    for i = 1:nbins
        if verbose
            @printf("[%.4f, %.4f) ==> (%d, %d): %d\n", edges[i], edges[i+1], clb[i], cub[i], cnts[i])
        end
        clb[i] <= cnts[i] <= cub[i] ||
            error("The counts are out of the confidence interval.")
    end
    return samples
end

test_samples(distr::ContinuousUnivariateDistribution, n::Int; nbins::Int=50, q::Float64=1.0e-6, verbose::Bool=false) =
    test_samples(distr, distr, n; nbins=nbins, q=q, verbose=verbose)

#### Testing range & support methods

function test_range(d::UnivariateDistribution)
    vmin = minimum(d)
    vmax = maximum(d)
    @test vmin <= vmax

    is_lb = islowerbounded(d)
    is_ub = isupperbounded(d)

    @test isfinite(vmin) == is_lb
    @test isfinite(vmax) == is_ub
    @test isbounded(d) == (is_lb && is_ub)
end

function get_evalsamples(d::DiscreteUnivariateDistribution, q::Float64)
    # samples for testing evaluation functions (even spacing)

    T = eltype(d)
    lv = (islowerbounded(d) ? minimum(d) : floor(T,quantile(d, q/2)))::T
    hv = (isupperbounded(d) ? maximum(d) : ceil(T,cquantile(d, q/2)))::T
    @assert lv <= hv
    return lv:hv
end

function get_evalsamples(d::ContinuousUnivariateDistribution, q::Float64, n::Int)
    # samples for testing evaluation functions (even spacing)

    lv = quantile(d, q/2)
    hv = cquantile(d, q/2)
    @assert isfinite(lv) && isfinite(hv) && lv <= hv
    return _linspace(lv, hv, n)
end

function test_support(d::UnivariateDistribution, vs::AbstractVector)
    for v in vs
        @test insupport(d, v)
    end
    @test all(insupport(d, vs))

    if islowerbounded(d)
        @test isfinite(minimum(d))
        @test insupport(d, minimum(d))
        @test !insupport(d, minimum(d)-1)
    end
    if isupperbounded(d)
        @test isfinite(maximum(d))
        @test insupport(d, maximum(d))
        @test !insupport(d, maximum(d)+1)
    end

    @test isbounded(d) == (isupperbounded(d) && islowerbounded(d))

    if isbounded(d)
        if isa(d, DiscreteUnivariateDistribution)
            s = support(d)
            @test isa(s, AbstractUnitRange)
            @test first(s) == minimum(d)
            @test last(s) == maximum(d)
        end
    end
end

#### Testing evaluation methods

function test_range_evaluation(d::DiscreteUnivariateDistribution)
    # check the consistency between range-based and ordinary pdf
    vmin = minimum(d)
    vmax = maximum(d)
    @test vmin <= vmax
    if islowerbounded(d)
        @test isa(vmin, Integer)
    end
    if isupperbounded(d)
        @test isa(vmax, Integer)
    end

    rmin = round(Int, islowerbounded(d) ? vmin : quantile(d, 0.001))::Int
    rmax = round(Int, isupperbounded(d) ? vmax : quantile(d, 0.999))::Int

    p0 = pdf.(Ref(d), collect(rmin:rmax))
    @test pdf.(Ref(d), rmin:rmax) ≈ p0
    if rmin + 2 <= rmax
        @test pdf.(Ref(d), rmin+1:rmax-1) ≈ p0[2:end-1]
    end

    if isbounded(d)
        @test pdf.(Ref(d), support(d)) ≈ p0
        @test pdf.(Ref(d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0)
        @test pdf.(Ref(d), rmin:rmax+3) ≈ vcat(p0, 0.0, 0.0, 0.0)
        @test pdf.(Ref(d), rmin-2:rmax+3) ≈ vcat(0.0, 0.0, p0, 0.0, 0.0, 0.0)
    elseif islowerbounded(d)
        @test pdf.(Ref(d), rmin-2:rmax) ≈ vcat(0.0, 0.0, p0)
    end
end

function test_evaluation(d::DiscreteUnivariateDistribution, vs::AbstractVector, testquan::Bool=true)
    nv  = length(vs)
    p   = Vector{Float64}(undef, nv)
    c   = Vector{Float64}(undef, nv)
    cc  = Vector{Float64}(undef, nv)
    lp  = Vector{Float64}(undef, nv)
    lc  = Vector{Float64}(undef, nv)
    lcc = Vector{Float64}(undef, nv)
    ci  = 0.

    for (i, v) in enumerate(vs)
        p[i] = pdf(d, v)
        c[i] = cdf(d, v)
        cc[i] = ccdf(d, v)
        lp[i] = logpdf(d, v)
        lc[i] = logcdf(d, v)
        lcc[i] = logccdf(d, v)

        @assert p[i] >= 0.0
        @assert (i == 1 || c[i] >= c[i-1])

        ci += p[i]
        @test ci ≈ c[i]
        @test isapprox(c[i] + cc[i], 1.0       , atol=1.0e-12)
        @test isapprox(lp[i]       , log(p[i]) , atol=1.0e-12)
        @test isapprox(lc[i]       , log(c[i]) , atol=1.0e-12)
        @test isapprox(lcc[i]      , log(cc[i]), atol=1.0e-12)

        if testquan
            ep = 1.0e-8
            if p[i] > 2 * ep   # ensure p[i] is large enough to guarantee a reliable result
                @test quantile(d, c[i] - ep) == v
                @test cquantile(d, cc[i] + ep) == v
                @test invlogcdf(d, lc[i] - ep) == v
                if 0.0 < c[i] < 1.0
                    @test invlogccdf(d, lcc[i] + ep) == v
                end
            end
        end
    end

    # consistency of scalar-based and vectorized evaluation
    @test pdf.(Ref(d), vs)  ≈ p
    @test cdf.(Ref(d), vs)  ≈ c
    @test ccdf.(Ref(d), vs) ≈ cc

    @test logpdf.(Ref(d), vs)  ≈ lp
    @test logcdf.(Ref(d), vs)  ≈ lc
    @test logccdf.(Ref(d), vs) ≈ lcc
end

function test_evaluation(d::ContinuousUnivariateDistribution, vs::AbstractVector, testquan::Bool=true)
    nv  = length(vs)
    p   = Vector{Float64}(undef, nv)
    c   = Vector{Float64}(undef, nv)
    cc  = Vector{Float64}(undef, nv)
    lp  = Vector{Float64}(undef, nv)
    lc  = Vector{Float64}(undef, nv)
    lcc = Vector{Float64}(undef, nv)

    for (i, v) in enumerate(vs)
        p[i] = pdf(d, v)
        c[i] = cdf(d, v)
        cc[i] = ccdf(d, v)
        lp[i] = logpdf(d, v)
        lc[i] = logcdf(d, v)
        lcc[i] = logccdf(d, v)

        @assert p[i] >= 0.0
        @assert (i == 1 || c[i] >= c[i-1])

        @test isapprox(c[i] + cc[i], 1.0       , atol=1.0e-12)
        @test isapprox(lp[i]       , log(p[i]) , atol=1.0e-12)
        @test isapprox(lc[i]       , log(c[i]) , atol=1.0e-12)
        @test isapprox(lcc[i]      , log(cc[i]), atol=1.0e-12)

        if testquan
            # TODO: remove this line when we have more accurate implementation
            # of quantile for InverseGaussian
            qtol = isa(d, InverseGaussian) ? 1.0e-4 : 1.0e-10
            if p[i] > 1.0e-6
                @test isapprox(quantile(d, c[i])    , v, atol=qtol * (abs(v) + 1.0))
                @test isapprox(cquantile(d, cc[i])  , v, atol=qtol * (abs(v) + 1.0))
                @test isapprox(invlogcdf(d, lc[i])  , v, atol=qtol * (abs(v) + 1.0))
                @test isapprox(invlogccdf(d, lcc[i]), v, atol=qtol * (abs(v) + 1.0))
            end
        end
    end

    # check: pdf should be the derivative of cdf
    for i = 2:(nv-1)
        if p[i] > 1.0e-6
            v = vs[i]
            ap = (cdf(d, v + 1.0e-6) - cdf(d, v - 1.0e-6)) / (2.0e-6)
            @test isapprox(p[i], ap, atol=p[i] * 1.0e-3)
        end
    end

    # consistency of scalar-based and vectorized evaluation
    @test pdf.(Ref(d), vs)  ≈ p
    @test cdf.(Ref(d), vs)  ≈ c
    @test ccdf.(Ref(d), vs) ≈ cc

    @test logpdf.(Ref(d), vs)  ≈ lp
    @test logcdf.(Ref(d), vs)  ≈ lc
    @test logccdf.(Ref(d), vs) ≈ lcc
end

#### Testing statistics methods

function test_stats(d::DiscreteUnivariateDistribution, vs::AbstractVector)
    # using definition (or an approximation)

    vf = Float64[v for v in vs]
    p = pdf.(Ref(d), vf)
    xmean = dot(p, vf)
    xvar = dot(p, abs2.(vf .- xmean))
    xstd = sqrt(xvar)
    xentropy = entropy(p)
    xskew = dot(p, (vf .- xmean).^3) / (xstd.^3)
    xkurt = dot(p, (vf .- xmean).^4) / (xvar.^2) - 3.0

    if isbounded(d)
        @test isapprox(mean(d), xmean, atol=1.0e-8)
        @test isapprox(var(d) , xvar , atol=1.0e-8)
        @test isapprox(std(d) , xstd , atol=1.0e-8)

        if applicable(skewness, d) && isfinite(skewness(d))
            @test isapprox(skewness(d), xskew   , atol=1.0e-8)
        end
        if applicable(kurtosis, d) && isfinite(kurtosis(d))
            @test isapprox(kurtosis(d), xkurt   , atol=1.0e-8)
        end
        if applicable(entropy, d)
            @test isapprox(entropy(d) , xentropy, atol=1.0e-8)
        end
    else
        @test isapprox(mean(d), xmean, atol=1.0e-3 * (abs(xmean) + 1.0))
        @test isapprox(var(d) , xvar , atol=0.01 * xvar)
        @test isapprox(std(d) , xstd , atol=0.01 * xstd)
    end
end

allow_test_stats(d::UnivariateDistribution) = true
allow_test_stats(d::NoncentralBeta) = false

function test_stats(d::ContinuousUnivariateDistribution, xs::AbstractVector{Float64})
    # using Monte Carlo methods

    if !(isfinite(mean(d)) && isfinite(var(d)))
        return
    end
    vd = var(d)

    n = length(xs)
    xmean = mean(xs)
    xvar = var(xs)
    xstd = sqrt(xvar)

    # we utilize central limit theorem, and consider xmean as (approximately) normally
    # distributed with std = std(d) / sqrt(n)
    #
    # For a normal distribution, it is extremely rare for a sample to deviate more
    # than 5 * std.dev, (chance < 6.0e-7)
    mean_tol = 5.0 * (sqrt(vd / n))
    @test isapprox(mean(d), xmean, atol=mean_tol)

    # test variance
    if applicable(kurtosis, d)
        kd = kurtosis(d)
        # when the excessive kurtosis is sufficiently large (i.e. > 2)
        # the sample variance has a finite variance, given by
        #
        #   (sigma^4 / n) * (k + 3 - (n-3)/(n-1))
        #
        # where k is the excessive kurtosis
        #
        if isfinite(kd) && kd > -2.0
            @test isapprox(var(d), xvar, atol=5.0 * vd * (kd + 2) / sqrt(n))
        end
    end
end

function test_params(d::Distribution)
    # simply test that params returns something sufficient to
    # reconstruct d
    D = typeof(d)
    pars = params(d)
    d_new = D(pars...)
    @test d_new == d
end

function test_params(d::Truncated)
    # simply test that params returns something sufficient to
    # reconstruct d
    d_unt = d.untruncated
    D = typeof(d_unt)
    pars = params(d_unt)
    d_new = Truncated(D(pars...), d.lower, d.upper)
    @test d_new == d
end


text:
struct CardinalDict{V,K} <: AbstractCardinalDictionary{V,K}

    present::BitArray{1}    # does this key hold content
    content::Vector{V}      # the keyed content (onlyif present)

    idxmax::K

    function CardinalDict(::Type{V}, maxidx::K) where {K<:Signed, V}
        present = falses(maxidx)
        content = Vector{V}(maxidx)
        T = type_for_indexing(maxidx)
        idxmax = T(maxidx)
        return new{V,T}(present, content, idxmax)
    end
end

function CardinalDict{V}(maxidx::K) where {K<:Signed, V}
   return CardinalDict(V, maxidx)
end

function CardinalDict(values::A) where A<:AbstractArray{T,1} where T
    n = length(values)
    dict = CardinalDict{T}(n)
    for (idx,val) in enumerate(values)
        dict[idx] = val
    end
    return dict
end


text:
import Base: size
import Base: stride
import Base: eachindex
import Base: getindex
import Base: setindex!
import Base: similar
import Base: copy
import Base: eltype
import Base: convert
import Base: collect
import Base.Broadcast: broadcast
import Base: hcat
import Base: vcat

"""
    Base.convert(::Type{SimpleSDMResponse}, layer::T) where {T <: SimpleSDMPredictor}

Returns a response with the same grid and bounding box as the predictor.
"""
function Base.convert(::Type{SimpleSDMResponse}, layer::T) where {T <: SimpleSDMPredictor}
   return copy(SimpleSDMResponse(layer.grid, layer.left, layer.right, layer.bottom, layer.top))
end

"""
    Base.convert(::Type{SimpleSDMPredictor}, layer::T) where {T <: SimpleSDMResponse}

Returns a predictor with the same grid and bounding box as the response.
"""
function Base.convert(::Type{SimpleSDMPredictor}, layer::T) where {T <: SimpleSDMResponse}
   return copy(SimpleSDMPredictor(layer.grid, layer.left, layer.right, layer.bottom, layer.top))
end

"""
    Base.convert(::Type{Matrix}, layer::T) where {T <: SimpleSDMLayer}

Returns the grid as an array.
"""
function Base.convert(::Type{Matrix}, layer::T) where {T <: SimpleSDMLayer}
   return copy(layer.grid)
end

"""
    Base.eltype(layer::SimpleSDMLayer{T}) where {T}

Returns the type of the values stored in the grid, where the `Nothing` type is
omitted.
"""
Base.eltype(::SimpleSDMResponse{T}) where {T} = T
Base.eltype(::SimpleSDMPredictor{T}) where {T} = T

"""
    Base.size(layer::T) where {T <: SimpleSDMLayer}

Returns the size of the grid.
"""
Base.size(layer::T) where {T <: SimpleSDMLayer} = size(layer.grid)

"""
    Base.size(layer::T, i...) where {T <: SimpleSDMLayer}

Returns the size of the grid alongside a dimension.
"""
Base.size(layer::T, i...) where {T <: SimpleSDMLayer} = size(layer.grid, i...)

"""
    Base.stride(layer::T; dims::Union{Nothing,Integer}=nothing) where {T <: SimpleSDMLayer}

Returns the stride, *i.e.* half the length, of cell dimensions, possibly
alongside a side of the grid. The first position is the length of the
*longitude* cells, the second the *latitude*.
"""
function Base.stride(layer::T; dims::Union{Nothing,Integer}=nothing) where {T <: SimpleSDMLayer}
   lon_stride = (layer.right-layer.left)/2.0size(layer, 2)
   lat_stride = (layer.top-layer.bottom)/2.0size(layer, 1)
   isnothing(dims) && return (lon_stride, lat_stride)
   dims == 1 && return lon_stride
   dims == 2 && return lat_stride
end
Base.stride(layer::T, i::Int) where {T<:SimpleSDMLayer} = stride(layer; dims=i)

"""
    Base.eachindex(layer::T) where {T <: SimpleSDMLayer}

Returns the index of the grid.
"""
Base.eachindex(layer::T) where {T <: SimpleSDMLayer} = eachindex(layer.grid)

"""
Extracts a  value from a layer by its grid position.
"""
Base.getindex(layer::T, i::Int64) where {T <: SimpleSDMLayer} = layer.grid[i]

"""
    Base.getindex(layer::T, i::R, j::R) where {T <: SimpleSDMLayer, R <: UnitRange}

Extracts a series of positions in a layer, and returns a layer corresponding to
the result. This is essentially a way to rapidly crop a layer to a given subset
of its extent. The `i` and `j` arguments are `UnitRange`s (of `Integer`).

The layer returned by this function will have the same type as the layer passed
as its argument, but this can be changed using `convert`. Note that this function
performs additional checks to ensure that the range is not empty, and to also
ensure that it does not overflows from the size of the layer.
"""
function Base.getindex(layer::T, i::R, j::R) where {T <: SimpleSDMLayer, R <: UnitRange}
   i_min = isempty(i) ? max(i.start-1, 1) : i.start
   i_max = isempty(i) ? max(i.stop+2, size(layer, 1)) : i.stop
   j_min = isempty(j) ? max(j.start-1, 1) : j.start
   j_max = isempty(j) ? max(j.stop+2, size(layer, 2)) : j.stop
   i_fix = i_min:i_max
   j_fix = j_min:j_max
   RT = T <: SimpleSDMResponse ? SimpleSDMResponse : SimpleSDMPredictor
   return RT(
            layer.grid[i_fix,j_fix],
            minimum(longitudes(layer)[j_fix])-stride(layer,1),
            maximum(longitudes(layer)[j_fix])+stride(layer,1),
            minimum(latitudes(layer)[i_fix])-stride(layer,2),
            maximum(latitudes(layer)[i_fix])+stride(layer,2)
           )
end

"""
Given a layer and a latitude, returns `nothing` if the latitude is outside the
range, or the grid index containing this latitude if it is within range
"""
function _match_latitude(layer::T, lat::K) where {T <: SimpleSDMLayer, K <: AbstractFloat}
   lat > layer.top && return nothing
   lat < layer.bottom && return nothing
   return last(findmin(abs.(lat .- latitudes(layer))))
end

"""
Given a layer and a longitude, returns `nothing` if the longitude is outside the
range, or the grid index containing this longitude if it is within range
"""
function _match_longitude(layer::T, lon::K) where {T <: SimpleSDMLayer, K <: AbstractFloat}
   lon > layer.right && return nothing
   lon < layer.left && return nothing
   return last(findmin(abs.(lon .- longitudes(layer))))
end

"""
    Base.getindex(layer::T, longitude::K, latitude::K) where {T <: SimpleSDMLayer, K <: AbstractFloat}

Extracts the value of a layer at a given latitude and longitude. If values
outside the range are requested, will return `nothing`.
"""
function Base.getindex(layer::T, longitude::K, latitude::K) where {T <: SimpleSDMLayer, K <: AbstractFloat}
   i = _match_longitude(layer, longitude)
   j = _match_latitude(layer, latitude)
   isnothing(i) && return nothing
   isnothing(j) && return nothing
   return layer.grid[j, i]
end

"""
    Base.getindex(layer::T; left=nothing, right=nothing, top=nothing, bottom=nothing) where {T <: SimpleSDMLayer, K <: Union{Nothing,AbstractFloat}}

Returns a subset of the argument layer, where the new limits are given by
`left`, `right`, `top`, and `bottom`. Up to three of these can be omitted, and
if so these limits will not be affected.
"""
function Base.getindex(layer::T; left=nothing, right=nothing, top=nothing, bottom=nothing) where {T <: SimpleSDMLayer}
   for limit in [left, right, top, bottom]
      if !isnothing(limit)
         @assert typeof(limit) <: AbstractFloat
      end
   end
   imax = _match_longitude(layer, isnothing(right) ? layer.right : right)
   imin = _match_longitude(layer, isnothing(left) ? layer.left : left)
   jmax = _match_latitude(layer, isnothing(top) ? layer.top : top)
   jmin = _match_latitude(layer, isnothing(bottom) ? layer.bottom : bottom)
   any(isnothing.([imin, imax, jmin, jmax])) && throw(ArgumentError("Unable to extract, coordinates outside of range"))
   # Note that this is LATITUDE first
   return layer[jmin:jmax, imin:imax]
end

"""
    Base.getindex(layer::T, n::NT) where {T <: SimpleSDMLayer, NT <: NamedTuple}

Returns a subset of the argument layer, where the new limits are given in
a NamedTuple by `left`, `right`, `top`, and `bottom`, in any order. Up to
three of these can be omitted, and if so these limits will not be affected.
"""
function Base.getindex(layer::T, n::NT) where {T <: SimpleSDMLayer, NT <: NamedTuple}
    l = isdefined(n, :left) ? n.left : nothing
    r = isdefined(n, :right) ? n.right : nothing
    t = isdefined(n, :top) ? n.top : nothing
    b = isdefined(n, :bottom) ? n.bottom : nothing
    Base.getindex(layer; left=l, right=r, top=t, bottom=b)
end

"""
    Base.getindex(layer1::T1, layer2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}

Extract a layer based on a second layer. Note that the two layers must be
*compatible*, which is to say they must have the same bounding box and grid
size.
"""
function Base.getindex(layer1::T1, layer2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}
   SimpleSDMLayers._layers_are_compatible(layer1, layer2)
   return layer1[left=layer2.left, right=layer2.right, bottom=layer2.bottom, top=layer2.top]
end

"""
     Base.setindex!(layer::SimpleSDMResponse{T}, v::T, i...) where {T}

Changes the value of a cell, or a range of cells, as indicated by their grid
positions.
"""
function Base.setindex!(layer::SimpleSDMResponse{T}, v::T, i...) where {T}
   typeof(v) <: eltype(layer.grid) || throw(ArgumentError("Impossible to set a value to a non-matching type"))
   layer.grid[i...] = v
end

"""
    Base.setindex!(layer::T, v, lon::Float64, lat::Float64) where {T <: SimpleSDMResponse}

Changes the values of the cell including the point at the requested latitude and
longitude.
"""
function Base.setindex!(layer::SimpleSDMResponse{T}, v::T, lon::Float64, lat::Float64) where {T}
   i = _match_longitude(layer, lon)
   j = _match_latitude(layer, lat)
   layer[j,i] = v
end

"""
    Base.similar(l::T) where {T <: SimpleSDMLayer}

Returns a `SimpleSDMResponse` of the same dimensions as the original layer, with
`nothing` in the same positions. The rest of the values are replaced by the
output of `zero(eltype(layer.grid))`, which implies that there must be a way to
get a zero for the type. If not, the same result can always be achieved through
the use of `copy`, manual update, and `convert`.
"""
function Base.similar(layer::T) where {T <: SimpleSDMLayer}
   emptygrid = similar(layer.grid)
   for i in eachindex(emptygrid)
      emptygrid[i] = isnothing(layer.grid[i]) ? nothing : zero(eltype(layer.grid[i]))
   end
   return SimpleSDMResponse(emptygrid, layer.left, layer.right, layer.bottom, layer.top)
end

"""
    Base.similar(::Type{TC}, l::T) where {TC <: Any, T <: SimpleSDMLayer}

Returns a `SimpleSDMResponse` of the same dimensions as the original layer, with
`nothing` in the same positions. The rest of the values are replaced by the
output of `zero(TC)`, which implies that there must be a way to get a zero for
the type. If not, the same result can always be achieved through the use of
`copy`, manual update, and `convert`.
"""
function Base.similar(::Type{TC}, layer::T) where {TC <: Any, T <: SimpleSDMLayer}
   emptygrid = convert(Matrix{Union{Nothing,TC}}, zeros(TC, size(layer)))
   emptygrid[findall(isnothing, layer.grid)] .= nothing
   return SimpleSDMResponse(emptygrid, layer.left, layer.right, layer.bottom, layer.top)
end

"""
    Base.copy(l::T) where {T <: SimpleSDMLayer}

Returns a new copy of the layer, which has the same type.
"""
function Base.copy(layer::T) where {T <: SimpleSDMLayer}
   copygrid = copy(layer.grid)
   RT = T <: SimpleSDMResponse ? SimpleSDMResponse : SimpleSDMPredictor
   return RT(copygrid, copy(layer.left), copy(layer.right), copy(layer.bottom), copy(layer.top))
end

"""
   Broadcast.broadcast(f, L::LT) where {LT <: SimpleSDMLayer}

TODO
"""
function Base.Broadcast.broadcast(f, L::LT) where {LT <: SimpleSDMLayer}
    newgrid = Array{Any}(nothing, size(L))
    N = SimpleSDMResponse(newgrid, L)
    v = filter(!isnothing, L.grid)
    fv = f.(v)
    N.grid[findall(!isnothing, L.grid)] .= fv

    internal_types = unique(typeof.(N.grid))

    RT = LT <: SimpleSDMResponse ? SimpleSDMResponse : SimpleSDMPredictor
    return RT(convert(Matrix{Union{internal_types...}}, N.grid), N)
end

"""
    Base.collect(l::T) where {T <: SimpleSDMLayer}

Returns the non-`nothing` values of a layer.
"""
function Base.collect(l::T) where {T <: SimpleSDMLayer}
    v = filter(!isnothing, l.grid)
    return convert(Vector{typeof(v[1])}, v)
end

"""
    Base.vcat(l1::T, l2::T) where {T <: SimpleSDMLayers}

Adds the second layer *under* the first one, assuming the strides and left/right
coordinates match. This will automatically re-order the layers if the second is
above the first.
"""
function Base.vcat(l1::T, l2::T) where {T <: SimpleSDMLayer}
    (l1.left == l2.left) || throw(ArgumentError("The two layers passed to vcat must have the same left coordinate"))
    (l1.right == l2.right) || throw(ArgumentError("The two layers passed to vcat must have the same right coordinate"))
    all(stride(l1) .≈ stride(l2)) || throw(ArgumentError("The two layers passed to vcat must have the same stride"))
    (l2.top == l1.bottom) && return vcat(l2, l1)
    new_grid = vcat(l1.grid, l2.grid)
    RT = T <: SimpleSDMPredictor ? SimpleSDMPredictor : SimpleSDMResponse
    return RT(new_grid, l1.left, l1.right, l1.top, l2.bottom)
end

"""
    Base.hcat(l1::T, l2::T) where {T <: SimpleSDMLayers}

Adds the second layer *to the right of* the first one, assuming the strides and
left/right coordinates match. This will automatically re-order the layers if the
second is to the left the first.
"""
function Base.hcat(l1::T, l2::T) where {T <: SimpleSDMLayer}
    (l1.top == l2.top) || throw(ArgumentError("The two layers passed to hcat must have the same top coordinate"))
    (l1.bottom == l2.bottom) || throw(ArgumentError("The two layers passed to hcat must have the same bottom coordinate"))
    all(stride(l1) .≈ stride(l2)) || throw(ArgumentError("The two layers passed to hcat must have the same stride"))
    (l2.right == l1.left) && return hcat(l2, l1)
    new_grid = hcat(l1.grid, l2.grid)
    RT = T <: SimpleSDMPredictor ? SimpleSDMPredictor : SimpleSDMResponse
    return RT(new_grid, l1.left, l2.right, l1.top, l1.bottom)
end

"""
    quantile(layer::T, p) where {T <: SimpleSDMLayer}

Returns the quantiles of `layer` at `p`, using `Statistics.quantile`.
"""
function Statistics.quantile(layer::T, p) where {T <: SimpleSDMLayer}
    return quantile(collect(layer), p)
end


text:
using JuMP, JuMPeR, Gurobi, Random, Distributions, LinearAlgebra, Plots

n_moves = 2
n_players = 3
# Some payoff matrix
payoffM = Dict(1 => zeros(2,2,2), 2 => zeros(2,2,2), 3 => zeros(2,2,2))
# Jean's payoff
payoffM[1][:,:,1] = [3 2; 1 5]
payoffM[1][:,:,2] = [2 4; 0 7]
# Enrico's payoff
payoffM[2][:,:,1] = [5 4; 3 2]
payoffM[2][:,:,2] = [4 6; 2 7]
# Guy's payoff
payoffM[3][:,:,1] = [8 1; 2 4]
payoffM[3][:,:,2] = [5 3; 1 7]

function collaborative_game()
    m = Model(solver = GurobiSolver())
    @variable(m, 0 <= strategy[1:n_moves, 1:n_moves, 1:n_moves] <= 1)
    @variable(m, payoff[1:n_players])
    @constraint(m, sum(strategy) == 1)
    for i = 1:n_players
        @constraint(m, payoff[i] == sum(payoffM[i] .* strategy))
    end
    @objective(m, Max, sum(payoff))
    return m
end

function invididual_game(ind=1)
    m = Model(solver = GurobiSolver());
    @variable(m, 0 <= move <= 1);
    @variable(m, payoff >= 0);
    all_opposing_strats = reshape(collect(Base.product(zip([1,2], [1,2])...)), (1, n_moves^(n_players-1)))
    @objective(m, Max, sum(dot(payoffM[ind][:, all_opposing_strats[i][1], all_opposing_strats[i][2]], [move, 1-move])
                for i=1:length(all_opposing_strats))); # Average return
    solve(m)
    println("Average return fixed move: $(getvalue(move)) with payoff $(getobjectivevalue(m)).")

    for i = 1:length(all_opposing_strats)
        @constraint(m, payoff <= sum(dot(payoffM[ind][:, all_opposing_strats[i][1], all_opposing_strats[i][2]],
                    [move, 1-move]))) # Minimum return
    end
    @objective(m, Max, payoff)
    solve(m)
    println("Max-min return, fixed move: $(getvalue(move)) with payoff $(getvalue(payoff)).")

    # Adaptive strategy
    m = Model(solver = GurobiSolver());
    @variable(m, 0 <= move <= 1);
    @variable(m, strat[1:n_players-1])
    @variable(m, strat0)
    @variable(m, 0 <= opp_moves[1:n_players-1] <= 1)
    @constraint(m, move == dot(strat, opp_moves) + strat0)
    @variable(m, payoff >= 0);

    # TODO: clean-up!
    @variable(m, zero_payoffs[1:n_players-1] >= 0)
    @variable(m, one_payoffs[1:n_players-1] >= 0)
    @constraint(m, zero_payoffs .<= opp_moves' * payoffM[ind][:,:,1] * opp_moves)
    @constraint(m, one_payoffs .<= payoffM[ind][:,:,2] * opp_moves)
    @objective(m, Max, sum(zero_payoffs) + sum(one_payoffs))
    @constraint(m, payoff <= )
    @objective(m, Max, payoff)


text:
"""based on Base.binomial, but without negative values for n and without overflow checks
(index calculations here should not overflow if the array does not have more elements than an Int64 can represent)"""
@inline function binomial_simple(n::T, k::T) where T<:Integer
    (k < 0 || k > n) && return zero(T)
    (k == 0 || k == n) && return one(T)
    if k > (n>>1)
        k = n - k
    end
    k == 1 && return n
    x::T = nn = n - k + 1
    nn += 1
    rr = 2
    while rr <= k
        x = div(x*nn, rr)
        rr += 1
        nn += 1
    end
    x
end

@inline @generated function binomial_unrolled(n::T, ::Val{k}) where {k,T<:Integer}
    terms = [:(n - $(k-j)) for j=1:k]
    if k < 10
        # for small k, just return the unrolled calculation directly
        # (n k) = prod(n+i-k, i=1:k)/k!
        # typemax(Int)/factorial(9) ~ 23*10^12
        # so for k<10, this does not overflow for index calculation of up to ~100TB arrays.
        # so we do not worry about it
        # use the precomputed factorial
        binom = :( *($(terms...)) ÷ $(factorial(k)) )
    else
        binom = terms[1] # j=1
        for j=2:k
            # careful about operation order:
            # first multiply, the product is then always divisible by j
            binom = :( ($binom * $(terms[j])) ÷ $j )
        end
        # (n k) == (n n-k) -> should be faster for n-k < k
        # but when we do this replacement, we cannot unroll the loop explicitly,
        # so heuristically use n-k < k/2 to ensure that it wins
        :( n < $(3k÷2) ? binomial_simple(n,n-$k) : $binom )
    end
end

"""calculate binomial(ii+n+offset,n)
This shows up in size and index calculations for arrays with symmetric indices."""
#@inline symind_binomial(ii,::Val{n},::Val{offset}) where {n,offset} = binomial_unrolled(ii+(n+offset),Val(n))

# binary search for finding an integer m such that func(m) <= ind < func(m+1), with low <= m <= high
function searchlast_func(x,func,low::T,high::T) where T<:Integer
    high += one(T)
    while low < high-one(T)
        mid = (low+high) >> 1
        if x < func(mid)
            high = mid
        else
            low = mid
        end
    end
    return low
end


text:
using BinaryBuilder, Pkg

include("../../fancy_toys.jl")

name = "CUTENSOR"
version = v"1.3.0"#.3

platforms_and_sources = Dict(
    Platform("x86_64", "linux") => [
        ArchiveSource("https://developer.download.nvidia.com/compute/cutensor/1.3.0/local_installers/libcutensor-linux-x86_64-1.3.0.3.tar.gz",
                      "80c3124497bd8692fd4e98376589372eee616a55f875a2a28961087322d45d00")],
    Platform("powerpc64le", "linux") => [
        ArchiveSource("https://developer.download.nvidia.com/compute/cutensor/1.3.0/local_installers/libcutensor-linux-ppc64le-1.3.0.3.tar.gz",
                      "5e60569be31be0377fec94392be3c7b844181dbf9b0731dd87e87d0b9ff38258")],
    Platform("x86_64", "windows") => [
        ArchiveSource("https://developer.download.nvidia.com/compute/cutensor/1.3.0/local_installers/libcutensor-windows-x86_64-1.3.0.3.zip",
                      "599ad469846d9a22d608a768d45c29138c6ba03247593aee56a999306489e840")],
)

products = [
    LibraryProduct(["libcutensor", "cutensor"], :libcutensor, dont_dlopen = true),
]

# XXX: CUDA_loader_jll's CUDA tag should match the library's CUDA version compatibility.
#      lacking that, we can't currently dlopen the library

dependencies = [
    Dependency(PackageSpec(name="CUDA_loader_jll")),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll",
                           uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

cuda_versions = [v"10.2", v"11.0", v"11.1", v"11.2", v"11.3"]
for cuda_version in cuda_versions
    cuda_tag = "$(cuda_version.major).$(cuda_version.minor)"
    include("build_$(cuda_tag).jl")

    for (platform, sources) in platforms_and_sources
        should_build_platform(triplet(platform)) || continue
        platform.tags["cuda"] = cuda_tag

        build_tarballs(ARGS, name, version, sources, script, [platform], products, dependencies;
                       lazy_artifacts=true)
    end
end


text:
abstract type AbstractLine <: AbstractShape end

struct VerticalLine{I <: Integer} <: AbstractLine
    i_min::I
    i_max::I
    j::I
end

struct HorizontalLine{I <: Integer} <: AbstractLine
    i::I
    j_min::I
    j_max::I
end

struct Line{I <: Integer} <: AbstractLine
    point1::Point{I}
    point2::Point{I}
end

struct ThickLine{I <: Integer} <: AbstractLine
    point1::Point{I}
    point2::Point{I}
    thickness::I
end

#####
##### VerticalLine
#####

is_valid(shape::VerticalLine) = shape.i_min <= shape.i_max

get_i_min(shape::VerticalLine) = shape.i_min
get_i_max(shape::VerticalLine) = shape.i_max

get_j_min(shape::VerticalLine) = shape.j
get_j_max(shape::VerticalLine) = shape.j

clip(image::AbstractMatrix, shape::VerticalLine) = VerticalLine(max(get_i_min(shape), get_i_min(image)), min(get_i_max(shape), get_i_max(image)), shape.j)

get_drawing_optimization_style(::VerticalLine) = CLIP

function draw!(f::F, image::AbstractMatrix, shape::VerticalLine, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    i_min = shape.i_min
    i_max = shape.i_max
    j = shape.j

    for i in i_min:i_max
        f(image, i, j, color)
    end

    return nothing
end

#####
##### HorizontalLine
#####

is_valid(shape::HorizontalLine) = shape.j_min <= shape.j_max

get_i_min(shape::HorizontalLine) = shape.i
get_i_max(shape::HorizontalLine) = shape.i

get_j_min(shape::HorizontalLine) = shape.j_min
get_j_max(shape::HorizontalLine) = shape.j_max

clip(image::AbstractMatrix, shape::HorizontalLine) = HorizontalLine(shape.i, max(get_j_min(shape), get_j_min(image)), min(get_j_max(shape), get_j_max(image)))

get_drawing_optimization_style(::HorizontalLine) = CLIP

function draw!(f::F, image::AbstractMatrix, shape::HorizontalLine, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    i = shape.i
    j_min = shape.j_min
    j_max = shape.j_max

    for j in j_min:j_max
        f(image, i, j, color)
    end

    return nothing
end

#####
##### Line
#####

get_i_min(shape::Line) = min(shape.point1.i, shape.point2.i)
get_i_max(shape::Line) = max(shape.point1.i, shape.point2.i)
get_i_extrema(shape::Line) = minmax(shape.point1.i, shape.point2.i)

get_j_min(shape::Line) = min(shape.point1.j, shape.point2.j)
get_j_max(shape::Line) = max(shape.point1.j, shape.point2.j)
get_j_extrema(shape::Line) = minmax(shape.point1.j, shape.point2.j)

get_drawing_optimization_style(::Line) = CHECK_BOUNDS

function draw!(f::F, image::AbstractMatrix, shape::Line, color) where {F <: Function}
    point1 = shape.point1
    point2 = shape.point2

    i1 = point1.i
    j1 = point1.j
    i2 = point2.i
    j2 = point2.j

    I = typeof(i1)

    di = abs(i2 - i1)
    dj = -abs(j2 - j1)
    si = i1 < i2 ? one(I) : -one(I)
    sj = j1 < j2 ? one(I) : -one(I)
    err = di + dj

    while true
        f(image, i1, j1, color)

        if (i1 == i2 && j1 == j2)
            break
        end

        e2 = convert(I, 2) * err

        if (e2 >= dj)
            err += dj
            i1 += si
        end

        if (e2 <= di)
            err += di
            j1 += sj
        end
    end

    return nothing
end

#####
##### ThickLine
#####

is_valid(shape::ThickLine) = shape.thickness > zero(shape.thickness)

get_i_min(shape::ThickLine) = min(shape.point1.i, shape.point2.i) - shape.thickness ÷ 2

function get_i_max(shape::ThickLine)
    point1 = shape.point1
    point2 = shape.point2
    thickness = shape.thickness

    i_max = max(point1.i, point2.i)

    I = typeof(thickness)

    i_max = max(point1.j, point2.j)

    half_thickness = thickness ÷ 2

    return i_max - half_thickness + thickness - one(I)
end

get_j_min(shape::ThickLine) = min(shape.point1.j, shape.point2.j) - shape.thickness ÷ 2

function get_j_max(shape::ThickLine)
    point1 = shape.point1
    point2 = shape.point2
    thickness = shape.thickness

    I = typeof(thickness)

    j_max = max(point1.j, point2.j)

    half_thickness = thickness ÷ 2

    return j_max - half_thickness + thickness - one(I)
end

get_drawing_optimization_style(::ThickLine) = CHECK_BOUNDS

function draw!(f::F, image::AbstractMatrix, shape::ThickLine, color) where {F <: Function}
    @assert is_valid(shape) "Cannot draw invalid shape $(shape)"

    point1 = shape.point1
    point2 = shape.point2
    thickness = shape.thickness

    I = typeof(thickness)

    i1 = point1.i
    j1 = point1.j
    i2 = point2.i
    j2 = point2.j

    half_thickness = thickness ÷ 2

    draw!(image, Line(point1, point2), color) do image, i, j, color
        draw!(f, image, FilledRectangle(Point(i - half_thickness, j - half_thickness), thickness, thickness), color)
    end

    return nothing
end


text:
using Compat

# Implement user-defined MPI reduction operations, by passing Julia
# functions as callbacks to MPI.

# for defining thread-local variables; can use Compat
# once JuliaLang/Compat.jl#223 is resolved.
if isdefined(Base, :Threads)
    import Base.Threads: nthreads, threadid
else
    nthreads() = 1
    threadid() = 1
end

# Unfortunately, MPI_Op_create takes a function that does not accept
# a void* "thunk" parameter, making it impossible to fully simulate
# a closure.  So, we have to use a global variable instead.   (Since
# the reduction functions are barriers, being re-entrant is probably
# not important in practice, fortunately.)   For MPI_THREAD_MULTIPLE
# using Julia native threading, however, we do make this global thread-local
const _user_functions = Array{Function}(undef, 1) # resized to nthreads() at runtime
const _user_op = Op(MPI_OP_NULL) # _mpi_user_function operation, initialized below

# C callback function corresponding to MPI_User_function
function _mpi_user_function(_a::Ptr{Nothing}, _b::Ptr{Nothing}, _len::Ptr{Cint}, t::Ptr{Cint})
    len = unsafe_load(_len)
    T = mpitype_dict_inverse[unsafe_load(t)]
    a = Ptr{T}(_a)
    b = Ptr{T}(_b)
    f = _user_functions[threadid()]
    for i = 1:len
        unsafe_store!(b, f(unsafe_load(a,i), unsafe_load(b,i)), i)
    end
    return nothing
end

function user_op(opfunc::Function)
    # we must initialize these at runtime, but it can't be done in __init__
    # since MPI.Init is not called yet.  So we do it lazily here:
    if _user_op.val == MPI_OP_NULL
        # FIXME: to be thread-safe, there should really be a mutex lock
        # of some sort so that this initialization only occurs once.
        # To do when native threading in Julia stabilizes (and is documented).
        resize!(_user_functions, nthreads())
        user_function = @cfunction(_mpi_user_function, Nothing, (Ptr{Nothing}, Ptr{Nothing}, Ptr{Cint}, Ptr{Cint}))
        opnum = Ref{Cint}()
        ccall(MPI_OP_CREATE, Nothing, (Ptr{Nothing}, Ref{Cint}, Ref{Cint}, Ref{Cint}),
             user_function, false, opnum, 0)
        _user_op.val = opnum[]
    end

    _user_functions[threadid()] = opfunc
    return _user_op
end

# use function types in Julia 0.7 to automatically use built-in
# MPI operations for the corresponding Julia functions.
for (f,op) in ((+,SUM), (*,PROD),
               (min,MIN), (max,MAX),
               (&, BAND), (|, BOR), (⊻, BXOR))
    @eval user_op(::$(typeof(f))) = $op
end

Allreduce!(sendbuf::MPIBuffertype{T}, recvbuf::MPIBuffertype{T},
           count::Integer, opfunc::Function, comm::Comm) where {T} =
    Allreduce!(sendbuf, recvbuf, count, user_op(opfunc), comm)

Reduce(sendbuf::MPIBuffertype{T}, count::Integer,
       opfunc::Function, root::Integer, comm::Comm) where {T} =
    Reduce(sendbuf, count, user_op(opfunc), root, comm)


text:
function compute_derivatives(obj::FieldOrdinaryDifferentialEquation, arg0::RealFieldElement, arg1::Vector{RealFieldElement})
    return jcall(obj, "computeDerivatives", Vector{RealFieldElement}, (RealFieldElement, Vector{RealFieldElement}), arg0, arg1)
end

function get_dimension(obj::FieldOrdinaryDifferentialEquation)
    return jcall(obj, "getDimension", jint, ())
end

function init(obj::FieldOrdinaryDifferentialEquation, arg0::RealFieldElement, arg1::Vector{RealFieldElement}, arg2::RealFieldElement)
    return jcall(obj, "init", void, (RealFieldElement, Vector{RealFieldElement}, RealFieldElement), arg0, arg1, arg2)
end


text:
rotate_plus_90(vec::Vector2D) = typeof(vec)(-vec[2], vec[1])
rotate_minus_90(vec::Vector2D) = typeof(vec)(vec[2], -vec[1])
rotate_180(vec::Vector2D) = -vec

rotate(x, y, c, s) = Vector2D(c * x - s * y, s * x + c * y)
rotate(vec::Vector2D, theta) = rotate(vec[1], vec[2], cos(theta), sin(theta))
rotate(vec::Vector2D, dir::Vector2D) = rotate(vec[1], vec[2], dir[1], dir[2])

get_relative_direction(dir1::Vector2D, dir2::Vector2D) = typeof(dir1)(dir2[1] * dir1[1] + dir2[2] * dir1[2], dir2[2] * dir1[1] - dir2[1] * dir1[2])
invert_relative_direction(dir::Vector2D) = typeof(dir)(dir[1], -dir[2])

function invert(pos::Vector2D, dir::Vector2D)
    inv_dir = invert_relative_direction(dir)
    inv_pos = -rotate(pos, inv_dir)
    return inv_pos, inv_dir
end


text:
using Mill, Flux

root_node = BagNode(BagNode(ArrayNode(rand(3,4)),[1:2, 3:4]), [1:2])
root_model = Mill.reflectinmodel(root_node, d -> Dense(d, 4), d -> SegmentedMeanMax(d))
a, m = root_node.data, root_model.im
ϕ = (x...) -> Mill.onoff(x..., 1)

explain(root_node, root_model, (x...) -> Mill.onoff(x..., 3), explaining_fun)


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

using RLESUtils, FileUtils

const REFDIR = "../quicktest_ref"
const RESULTDIR = "./results"

include("runtests.jl")

"""
checks if quicktest results are same as reference set in repo
returns true if all files are identical
"""
function reftest(doquicktest::Bool=true)
    if doquicktest
        runscript("config_2ac_quicktest.ini")
    end
    fs = readdir(REFDIR)

    for f in fs
        if !isidentical(joinpath(REFDIR,f), joinpath(RESULTDIR,f))
            return false
        end
    end
    return true
end

#reftest()

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
using ..SearchProblem: AbstractSearchProblem, successors, cost

mutable struct TreeSearchNode{S, A}
    state::S
    path_cost::Float64
    depth::Int
    action::Union{A,Nothing}
    parent::Union{TreeSearchNode{S,A},Nothing}
    popped_using_w::Float64
    function TreeSearchNode{S, A}(state::S, path_cost::Real, depth::Integer, action::Union{A, Nothing}, parent::Union{TreeSearchNode,Nothing}) where {S, A}
        return new{S, A}(state, path_cost, depth, action, parent, 0.0)
    end
end

TreeSearchNode{S,A}(state) where {S,A} = TreeSearchNode{S,A}(state, 0, 0, nothing, nothing)

function get_solution(node::TreeSearchNode{S,A})::Tuple{Vector{A}, Vector{Float64}} where {S,A}
    if node.depth == 1
        return A[node.action], Float64[node.popped_using_w]
    else
        sol, w_trace = get_solution(node.parent)
        push!(sol, node.action)
        push!(w_trace, node.popped_using_w)
        return sol, w_trace
    end
end

function get_children_nodes(sp::AbstractSearchProblem{S, A}, node::TreeSearchNode{S, A}) where {S, A}
    return Iterators.map(successors(sp, node.state)) do successor
        successor_state::S, successor_action::A = successor
        path_cost = node.path_cost + cost(sp, node.state, successor_action, successor_state)
        child_node = TreeSearchNode{S, A}(successor_state, path_cost, node.depth + 1, successor_action, node)
        return child_node
    end
end


text:
export Specification, Goal
export is_goal, is_violated, get_cost, get_reward

"Abstract specification for a planning problem"
abstract type Specification end

"Check if `state` is a goal state according to the specification."
is_goal(spec::Specification, domain::Domain, state::State) =
    error("Not implemented.")

"Check if `state` violates specified constraints."
is_violated(spec::Specification, domain::Domain, state::State) =
    error("Not implemented.")

"Returns the cost of going from state `s1` to state `s2` via action `a`."
get_cost(spec::Specification, domain::Domain, s1::State, a::Term, s2::State) =
    error("Not implemented.")

"Returns the reward of going from state `s1` to state `s2` via action `a`."
get_reward(spec::Specification, domain::Domain, s1::State, a::Term, s2::State) =
    error("Not implemented.")

"Returns the reward discount factor."
get_discount(spec::Specification) =
    error("Not implemented.")

"Null specification."
struct NullSpecification <: Specification end

is_goal(::NullSpecification, ::Domain, ::State) = false
is_violated(::NullSpecification, ::Domain, ::State) = false
get_cost(::NullSpecification, ::Domain, ::State, ::Term, ::State) = 0.0
get_reward(::NullSpecification, ::Domain, ::State, ::Term, ::State) = 0.0
get_discount(::NullSpecification) = 1.0

"Abstract type for goal-based specifications."
abstract type Goal <: Specification end

"Return goal terms."
get_goal_terms(spec::Goal) = error("Not implemented.")

# No discounting for goal specifications by default
get_discount(spec::Goal) = 1.0

"Null goal specification."
struct NullGoal <: Goal end

is_goal(::NullGoal, ::Domain, ::State) = true
is_violated(::NullGoal, ::Domain, ::State) = false
get_cost(::NullGoal, ::Domain, ::State, ::Term, ::State) = 0.0
get_reward(::NullGoal, ::Domain, ::State, ::Term, ::State) = 0.0
get_discount(::NullGoal) = 1.0
get_goal_terms(::NullGoal) = Term[]

include("specifications/min_steps.jl")
include("specifications/min_metric.jl")
include("specifications/max_metric.jl")
include("specifications/state_constrained.jl")
include("specifications/discounted.jl")
include("specifications/goal_reward.jl")
include("specifications/backward.jl")

# Convenience constructors
Specification(problem::Problem) = problem.metric === nothing ?
    MinStepsGoal(problem) : Specification(problem)
Specification(goals::AbstractVector{<:Term}) = MinStepsGoal(goals)
Specification(goal::Term) = MinStepsGoal(goal)


text:
import MLBase
import StatsBase

"""
"""
function getallrocnums(
        ytrue::AbstractVector{<:Integer},
        yscore::AbstractVector{<:Real};
        additionalthreshold::AbstractFloat = 0.5,
        )
    allthresholds = get_binary_thresholds(
        yscore;
        additionalthreshold = additionalthreshold,
        )
    allrocnums = MLBase.roc(
        ytrue,
        yscore,
        allthresholds,
        )
    return allrocnums, allthresholds
end

"""
"""
accuracy(x::MLBase.ROCNums) = (x.tp + x.tn)/(x.p + x.n)

"""
"""
true_positive_rate(x::MLBase.ROCNums) = (x.tp)/(x.p)

"""
"""
true_negative_rate(x::MLBase.ROCNums) = (x.tn)/(x.n)

"""
"""
false_positive_rate(x::MLBase.ROCNums) = (x.fp)/(x.n)

"""
"""
false_negative_rate(x::MLBase.ROCNums) = (x.fn)/(x.p)

"""
"""
function positive_predictive_value(x::MLBase.ROCNums)
    if (x.tp == 0) && (x.tp + x.fp == 0)
        result = 1
    elseif (x.tp != 0) && (x.tp + x.fp == 0)
        error("x.tp != 0) && (x.tp + x.fp == 0)")
    else
        result = (x.tp) / (x.tp + x.fp)
    end
    return result
end

"""
"""
function negative_predictive_value(x::MLBase.ROCNums)
    if (x.tn == 0) && (x.tn + x.fn ==0)
        result = 1
    elseif (x.tn != 0) && (x.tn + x.fn == 0)
        error("(x.tn != 0) && (x.tn + x.fn == 0)")
    else
        result = (x.tn) / (x.tn + x.fn)
    end
    return result
end

"""
"""
sensitivity(x::MLBase.ROCNums) = true_positive_rate(x)

"""
"""
specificity(x::MLBase.ROCNums) = true_negative_rate(x)

"""
"""
precision(x::MLBase.ROCNums) = positive_predictive_value(x)

"""
"""
recall(x::MLBase.ROCNums) = true_positive_rate(x)

"""
"""
function fbetascore(
        x::MLBase.ROCNums,
        beta::Real,
        )
    p = precision(x)
    r = recall(x)
    result = ( 1 + beta^2 ) * ( p*r ) / ( ((beta^2) * p) + r )
    return result
end

"""
"""
f1score(x::MLBase.ROCNums) = fbetascore(x, 1)


text:
using Symbolics
using Test
using Latexify
using ReferenceTests

@variables x y z u(x)
Dx = Differential(x)

# issue 260
@test Symbolics._toexpr(3*x/y) == :((3x) / y)
@test Symbolics._toexpr(3*x^y) == :(3x^y)

@test_reference "latexify_refs/inverse.txt" latexify(x^-1)

@test_reference "latexify_refs/frac1.txt" latexify((z + x*y^-1) / sin(z))
@test_reference "latexify_refs/frac2.txt"  latexify((3x - 7y*z^23) * (z - z^2) / x)

@test_reference "latexify_refs/minus1.txt" latexify(x - y)
@test_reference "latexify_refs/minus2.txt" latexify(x - y * z)
@test_reference "latexify_refs/minus3.txt" latexify(sin(x+y-z))

@test_reference "latexify_refs/unary_minus1.txt" latexify(-y)
@test_reference "latexify_refs/unary_minus2.txt" latexify(-y * z)

@test_reference "latexify_refs/derivative1.txt" latexify(Dx(y))
@test_reference "latexify_refs/derivative2.txt" latexify(Dx(u))
@test_reference "latexify_refs/derivative3.txt" latexify(Dx(x^2 + y^2 + z^2))

@test_reference "latexify_refs/stable_mul_ordering1.txt" latexify(x * y)
@test_reference "latexify_refs/stable_mul_ordering2.txt" latexify(y * x)

@test_reference "latexify_refs/equation1.txt" latexify(x ~ y + z)
@test_reference "latexify_refs/equation2.txt" latexify(x ~ Dx(y + z))

@test_reference "latexify_refs/equation_vec1.txt" latexify([
    x ~   y +  z
    y ~   x - 3z
    z ~ -5x + 4y
])
@test_reference "latexify_refs/equation_vec2.txt" latexify([
    Dx(u) ~   z
    Dx(y) ~   y*x
])


text:
@inline anytagged() = false
@inline anytagged(x::Tagged, args...) = true
@inline anytagged(x, args...) = anytagged(args...)

macro proptagcontext(name)
    quote
        Cassette.@context($name)

        ## promote(x,y) should not tag the output tuple with the union of tags
        ## of x and y. So here we recurse into promote, and then tag each
        ## element of the result with the original tag
        function Cassette.overdub(ctx::$name, f::typeof(promote), args...)
            promoted = Cassette.fallback(ctx, f, args...)

            # put the tags back on:
            tagged_promoted = map(args, promoted) do orig, prom
                if Cassette.hasmetadata(orig, ctx)
                    tag(prom, ctx, metadata(orig, ctx))
                else
                    prom
                end
            end
            Cassette.recurse(ctx, tuple, tagged_promoted...)
        end

        function Cassette.overdub(ctx::$name, f, args...)
            # this check can be inferred (in theory)
            if anytagged(args...)
                # This is a slower check
                if all(x->metatype(x, ctx) <: Cassette.NoMetaData, args)
                    if isleaf(f)
                        return Cassette.fallback(ctx, f, args...)
                    else
                        # there maybe closures closing over tagged values
                        return Cassette.recurse(ctx, f, args...)
                    end
                end
                if isleaf(f)
                    val = Cassette.fallback(ctx, f, args...)
                    return propagate_tags(ctx, f, val, args...)
                else
                    #@show f, args
                    val = Cassette.recurse(ctx, f, args...)
                    # Inputs were tagged but the output wasn't
                    # So just leave the input tags on.
                    if !(val isa Tagged) || metatype(val, ctx) <: Cassette.NoMetaData
                        return propagate_tags(ctx, f, val, args...)
                    else
                        return val
                    end
                end
            else
                if isleaf(f)
                    return Cassette.fallback(ctx, f, args...)
                else
                    return Cassette.recurse(ctx, f, args...)
                end
            end
        end
    end |> esc
end

"""
`propagate_tags(ctx, f, result, args...)`

Called only if any of the `args` are Tagged.
must return `result` or a tagged version of `result`.
"""
@inline function propagate_tags(ctx, f, result, args...)
    result
end


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 = "MKL_Headers"
version = v"2022.0.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource(
        "https://anaconda.org/intel/mkl-include/2022.0.1/download/linux-64/mkl-include-2022.0.1-intel_117.tar.bz2",
        "e664eb639faf03d5a4e76b435deb28a67ef7f59f3c55adaeeef771a7c94b56e3";
        unpack_target = "mkl-include-x86_64-linux-gnu"
    ),
    ArchiveSource(
        "https://anaconda.org/intel/mkl-include/2022.0.1/download/linux-32/mkl-include-2022.0.1-intel_117.tar.bz2",
        "411f636d3e95f75453de80a5574273454b67181c7c30d3276c6cb18a91ef3eeb";
        unpack_target = "mkl-include-i686-linux-gnu"
    ),
    ArchiveSource(
        "https://anaconda.org/intel/mkl-include/2022.0.0/download/osx-64/mkl-include-2022.0.0-intel_105.tar.bz2",
        "d7aa8f0c7904045c6f569462c0d04cc34250ee991855f5a3ab2fba8bbc5a0d81";
        unpack_target = "mkl-include-x86_64-apple-darwin14"
    ),
    ArchiveSource(
        "https://anaconda.org/intel/mkl-include/2022.0.0/download/win-32/mkl-include-2022.0.0-intel_115.tar.bz2",
        "2d9999a92959be25b5dabbd5b69a4b7f738b87215254c776d59c6b2f8b2addf3";
        unpack_target = "mkl-include-i686-w64-mingw32"
    ),
    ArchiveSource(
        "https://anaconda.org/intel/mkl-include/2022.0.0/download/win-64/mkl-include-2022.0.0-intel_115.tar.bz2",
        "eb8f6187d8e4ded8c48b701d00658ef49adc2180e42f1ff52f43c1fecc79a870";
        unpack_target = "mkl-include-x86_64-w64-mingw32"
    ),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/mkl-include-$target
if [[ $target == *-mingw* ]]; then
    rsync -av Library/include/ ${includedir}
else
    rsync -av include/ ${includedir}
fi
install_license info/licenses/*.txt
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
    Platform("i686", "windows"),
    Platform("x86_64", "windows"),
]

# The products that we will ensure are always built
products = [
    FileProduct("include/mkl.h", :include_mkl_h),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
@get_units @with_kw struct SurfaceFlow{T,R,L}
    β::T | "-"                              # constant in Manning's equation
    sl::Vector{T} | "m m-1"                 # Slope [m m⁻¹]
    n::Vector{T} | "s m-1/3"                # Manning's roughness [s m⁻⅓]
    dl::Vector{T} | "m"                     # Drain length [m]
    q::Vector{T} | "m3 s-1"                 # Discharge [m³ s⁻¹]
    qin::Vector{T} | "m3 s-1"               # Inflow from upstream cells [m³ s⁻¹]
    q_av::Vector{T} | "m3 s-1"              # Average discharge [m³ s⁻¹]
    qlat::Vector{T} | "m2 s-1"              # Lateral inflow per unit length [m² s⁻¹]
    inwater::Vector{T} | "m3 s-1"           # Lateral inflow [m³ s⁻¹]
    volume::Vector{T} | "m3"                # Kinematic wave volume [m³] (based on water level h)
    h::Vector{T} | "m"                      # Water level [m]
    h_av::Vector{T} | "m"                   # Average water level [m]
    Δt::T | "s"                             # Model time step [s]
    its::Int | "-"                          # Number of fixed iterations
    width::Vector{T} | "m"                  # Flow width [m]
    alpha_pow::T | "-"                      # Used in the power part of α
    alpha_term::Vector{T} | "-"             # Term used in computation of α
    α::Vector{T} | "s3/5 m1/5"              # Constant in momentum equation A = αQᵝ, based on Manning's equation
    eps::T | "s3/5 m1/5"                    # Maximum allowed change in α, if exceeded cross sectional area and h is recalculated
    cel::Vector{T} | "m s-1"                # Celerity of the kinematic wave
    to_river::Vector{T} | "m3 s-1"          # Part of overland flow [m³ s⁻¹] that flows to the river
    rivercells::Vector{Bool} | "-"          # Location of river cells (0 or 1)
    wb_pit::Vector{Bool} | "-"              # Boolean location (0 or 1) of a waterbody (wb, reservoir or lake).
    reservoir_index::Vector{Int} | "-"      # map cell to 0 (no reservoir) or i (pick reservoir i in reservoir field)
    lake_index::Vector{Int} | "-"           # map cell to 0 (no lake) or i (pick lake i in lake field)
    reservoir::R                            # Reservoir model struct of arrays
    lake::L                                 # Lake model struct of arrays

    # TODO unclear why this causes a MethodError
    # function SurfaceFlow{T,R,L}(args...) where {T,R,L}
    #     equal_size_vectors(args)
    #     return new(args...)
    # end
end

statevars(::SurfaceFlow) = (:q, :h, :h_av)

function update(
    sf::SurfaceFlow,
    network;
    frac_toriver = nothing,
    inflow_wb = nothing,
    do_iter = false,
    doy = 0,
)
    @unpack graph,
    order,
    subdomain_order,
    topo_subdomain,
    indices_subdomain,
    upstream_nodes = network

    n = length(order)
    ns = length(subdomain_order)
    # two options for iteration, fixed or based on courant number.
    if do_iter
        if sf.its > 0
            its = sf.its
        else
            # calculate celerity
            courant = zeros(n)
            for k = 1:ns
                @threads for m in subdomain_order[k]
                    for v in topo_subdomain[m]
                        if sf.q[v] > 0.0
                            sf.cel[v] = 1.0 / (sf.α[v] * sf.β * pow(sf.q[v], (sf.β - 1.0)))
                            courant[v] = (sf.Δt / sf.dl[v]) * sf.cel[v]
                        end
                    end
                end
            end
            filter!(x -> x ≠ 0.0, courant)
            its = isempty(courant) ? 1 : ceil(Int, (1.25 * quantile!(courant, 0.95)))
        end
    else
        its = 1
    end

    # sub time step
    adt = sf.Δt / its

    @. sf.alpha_term = pow(sf.n / sqrt(sf.sl), sf.β)
    @. sf.α = sf.alpha_term * pow(sf.width + 2.0 * sf.h, sf.alpha_pow)

    sf.q_av .= 0.0
    sf.h_av .= 0.0
    sf.to_river .= 0.0
    # because of possible iterations set reservoir and lake inflow and total outflow at
    # start to zero, the total sum of inflow and outflow at each sub time step is calculated
    if !isnothing(sf.reservoir)
        sf.reservoir.inflow .= 0.0
        sf.reservoir.totaloutflow .= 0.0
    end
    if !isnothing(sf.lake)
        sf.lake.inflow .= 0.0
        sf.lake.totaloutflow .= 0.0
    end

    for _ = 1:its
        sf.qin .= 0.0
        for k = 1:ns
            @threads for m in subdomain_order[k]
                for (n, v) in zip(indices_subdomain[m], topo_subdomain[m])

                    # for overland flow frac_toriver needs to be defined
                    if frac_toriver !== nothing # run kinematic wave for land domain
                        # for a river cell without a reservoir or lake (wb_pit is false) part of the
                        # upstream surface flow goes to the river (frac_toriver) and part goes to
                        # the surface flow reservoir (1.0 - frac_toriver), upstream nodes with a
                        # reservoir or lake are excluded
                        sf.to_river[v] += sum_at(
                            i -> sf.q[i] * frac_toriver[i],
                            upstream_nodes[n],
                            eltype(sf.to_river),
                        )
                        if sf.width[v] > 0.0
                            sf.qin[v] = sum_at(
                                i -> sf.q[i] * (1.0 - frac_toriver[i]),
                                upstream_nodes[n],
                                eltype(sf.q),
                            )
                        end
                    else # run kinematic wave for river domain (including reservoirs and lakes)
                        # sf.qin by outflow from upstream reservoir or lake location is added
                        sf.qin[v] += sum_at(sf.q, upstream_nodes[n])

                        if !isnothing(sf.reservoir) && sf.reservoir_index[v] != 0
                            # run reservoir model and copy reservoir outflow to inflow (qin) of
                            # downstream river cell
                            i = sf.reservoir_index[v]
                            update(sf.reservoir, i, sf.q[v] + inflow_wb[v], adt)
                            j = only(outneighbors(graph, v))
                            sf.qin[j] = sf.reservoir.outflow[i]

                        elseif !isnothing(sf.lake) && sf.lake_index[v] != 0
                            # run lake model and copy lake outflow to inflow (qin) of downstream river
                            # cell
                            i = sf.lake_index[v]
                            update(sf.lake, i, sf.q[v] + inflow_wb[v], doy, adt)
                            j = only(outneighbors(graph, v))
                            sf.qin[j] = sf.lake.outflow[i]
                        end
                    end

                    sf.q[v] = kinematic_wave(
                        sf.qin[v],
                        sf.q[v],
                        sf.qlat[v],
                        sf.α[v],
                        sf.β,
                        adt,
                        sf.dl[v],
                    )

                    # update alpha
                    crossarea = sf.α[v] * pow(sf.q[v], sf.β)
                    sf.h[v] = crossarea / sf.width[v]
                    wetper = sf.width[v] + (2.0 * sf.h[v]) # wetted perimeter
                    α = sf.α[v]
                    sf.α[v] = sf.alpha_term[v] * pow(wetper, sf.alpha_pow)

                    if abs(α - sf.α[v]) > sf.eps
                        crossarea = sf.α[v] * pow(sf.q[v], sf.β)
                        sf.h[v] = crossarea / sf.width[v]
                    end

                    sf.q_av[v] += sf.q[v]
                    sf.h_av[v] += sf.h[v]
                end

            end
        end
    end
    sf.q_av ./= its
    sf.h_av ./= its
    sf.to_river ./= its
    sf.volume .= sf.dl .* sf.width .* sf.h
end

@get_units @with_kw struct LateralSSF{T}
    kh₀::Vector{T} | "mm Δt-1"              # Horizontal hydraulic conductivity at soil surface [mm Δt⁻¹]
    f::Vector{T} | "mm-1"                   # A scaling parameter [mm⁻¹] (controls exponential decline of kh₀)
    soilthickness::Vector{T} | "mm"         # Soil thickness [mm]
    θₛ::Vector{T} | "mm mm-1"               # Saturated water content (porosity) [-]
    θᵣ::Vector{T} | "mm mm-1"               # Residual water content [-]
    t::T | "Δt s"                           # time step [Δt s]
    Δt::T | "s"                             # model time step [s]
    βₗ::Vector{T} | "m m-1"                  # Slope [m m⁻¹]
    dl::Vector{T}                           # Drain length [mm]
    dw::Vector{T}                           # Flow width [mm]
    zi::Vector{T}                           # Pseudo-water table depth [mm] (top of the saturated zone)
    exfiltwater::Vector{T}                  # Exfiltration [mm]  (groundwater above surface level, saturated excess conditions)
    recharge::Vector{T}                     # Net recharge to saturated store [mm]
    ssf::Vector{T} | "mm3 Δt-1"             # Subsurface flow [mm³ Δt⁻¹]
    ssfin::Vector{T} | "mm3 Δt-1"           # Inflow from upstream cells [mm³ Δt⁻¹]
    ssfmax::Vector{T} | "mm2 Δt-1"          # Maximum subsurface flow [mm² Δt⁻¹]
    to_river::Vector{T} | "mm3 Δt-1"        # Part of subsurface flow [mm³ Δt⁻¹] that flows to the river
    wb_pit::Vector{Bool} | "-"              # Boolean location (0 or 1) of a waterbody (wb, reservoir or lake).

    function LateralSSF{T}(args...) where {T}
        equal_size_vectors(args)
        return new(args...)
    end
end

statevars(::LateralSSF) = (:ssf,)

function update(ssf::LateralSSF, network, frac_toriver)
    @unpack subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes = network

    ns = length(subdomain_order)
    for k = 1:ns
        @threads for m in subdomain_order[k]
            for (n, v) in zip(indices_subdomain[m], topo_subdomain[m])
                # for a river cell without a reservoir or lake (wb_pit is false) part of the
                # upstream subsurface flow goes to the river (frac_toriver) and part goes to the
                # subsurface flow reservoir (1.0 - frac_toriver) upstream nodes with a reservoir or
                # lake are excluded
                ssf.ssfin[v] = sum_at(
                    i -> ssf.ssf[i] * (1.0 - frac_toriver[i]),
                    upstream_nodes[n],
                    eltype(ssf.ssfin),
                )
                ssf.to_river[v] = sum_at(
                    i -> ssf.ssf[i] * frac_toriver[i],
                    upstream_nodes[n],
                    eltype(ssf.to_river),
                )

                ssf.ssf[v], ssf.zi[v], ssf.exfiltwater[v] = kinematic_wave_ssf(
                    ssf.ssfin[v],
                    ssf.ssf[v],
                    ssf.zi[v],
                    ssf.recharge[v],
                    ssf.kh₀[v],
                    ssf.βₗ[v],
                    ssf.θₛ[v] - ssf.θᵣ[v],
                    ssf.f[v],
                    ssf.soilthickness[v],
                    ssf.t,
                    ssf.dl[v],
                    ssf.dw[v],
                    ssf.ssfmax[v],
                )
            end
        end
    end
end

@get_units @with_kw struct GroundwaterExchange{T}
    Δt::T | "s"                         # model time step [s]
    exfiltwater::Vector{T}              # Exfiltration [mm]  (groundwater above surface level, saturated excess conditions)
    zi::Vector{T}                       # Pseudo-water table depth [mm] (top of the saturated zone)
    to_river::Vector{T} | "mm3 Δt-1"    # Part of subsurface flow [mm³ Δt⁻¹] that flows to the river
    ssf::Vector{T} | "mm3 Δt-1"         # Subsurface flow [mm³ Δt⁻¹]
end


text:
include("../examples/ellipticintegral.jl")
include("../examples/erfc.jl")
include("../examples/Wiener-Hopf.jl")
include("../examples/Fracture_a.jl")
include("../examples/Fracture_b.jl")
include("../examples/Fracture_c.jl")
include("../examples/Laplace.jl")
include("../examples/HelmholtzDirichlet.jl")
include("../examples/HelmholtzNeumann.jl")


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

@recipe function f(
    samples::DensitySampleVector;
    vsel=collect(1:5),
    mean=false,
    std=false,
    bins = 200,
    globalmode=false,
    marginalmode=false,
    diagonal = Dict(),
    upper = Dict(),
    lower = Dict(),
    vsel_label = []
)
    vsel = collect(vsel)
    if isa(vsel, Vector{<:Integer}) == false
        vsel = asindex.(Ref(BAT.varshape(samples)), vsel)
        vsel = reduce(vcat, vsel)
    end
    vsel = vsel[vsel .<=  totalndof(BAT.varshape(samples))]

    xlabel = ["v$i" for i in vsel]
    ylabel = ["p(v$i)" for i in vsel]

    if isshaped(samples)
        xlabel = getstring.(Ref(samples), vsel)
        ylabel = ["p($l)" for l in xlabel]
    end

    if length(vsel_label) > 0
        xlabel = [vsel_label[i] for i in 1:length(vsel_label)]
        ylabel = ["p("*vsel_label[i]*")" for i in 1:length(vsel_label)]
    end

    nparams = length(vsel)
    layout --> nparams^2
    size --> (1000, 600)

    partition_tree = get(upper, "partition_tree", false)
    min_samples = minimum(flatview(unshaped.(samples.v)), dims=2)[:,1]
    max_samples = maximum(flatview(unshaped.(samples.v)), dims=2)[:,1]

    # get bin edges
    bins = if isa(bins, Integer)
        fill(bins, length(vsel))
    elseif isa(bins, Tuple{Vararg{Union{Integer, StepRangeLen}}})
        bins
    elseif isa(bins, NamedTuple)
        tmp_bins::Vector{Any} = fill(200, length(vsel))
        for k in keys(bins)
            idx = findfirst(isequal(string(k)), getstring.(Ref(samples), vsel))
            tmp_bins[idx] = bins[idx]
        end
        tmp_bins
    end

    for i in 1:nparams
        # diagonal
        @series begin
            subplot := i + (i-1)*nparams

            seriestype --> get(diagonal, "seriestype", :smallest_intervals)
            bins := bins[i]
            colors --> get(diagonal, "colors", default_colors)
            intervals --> get(diagonal, "intervals", default_credibilities)
            interval_labels --> get(diagonal, "interval_labels", [])
            legend --> get(diagonal, "legend", false)
            mean --> get(diagonal, "mean", mean)
            std --> get(diagonal, "std", std)
            globalmode --> get(diagonal, "globalmode", globalmode)
            marginalmode --> get(diagonal, "marginalmode", marginalmode)
            xlims --> get(diagonal, "xlims", :auto)
            ylims --> get(diagonal, "ylims", :auto)
            xguide --> xlabel[i]
            yguide --> ylabel[i]

            unshaped.(samples), (vsel[i])
        end

        # upper right plots
        for j in i+1:nparams
            @series begin
                subplot := j + (i-1)*nparams

                seriestype --> get(upper, "seriestype", :histogram)
                bins := (bins[i], bins[j])
                colors --> get(upper, "colors", default_colors)
                intervals --> get(upper, "intervals", default_credibilities)
                interval_labels --> get(upper, "interval_labels", [])
                legend --> get(upper, "legend", false)
                colorbar --> get(upper, "colorbar", false)
                mean --> get(upper, "mean", mean)
                std --> get(upper, "std", std)
                globalmode --> get(upper, "globalmode", globalmode)
                marginalmode --> get(upper, "marginalmode", marginalmode)
                xlims --> get(upper, "xlims", :auto)
                ylims --> get(upper, "ylims", :auto)

                if partition_tree != false
                    xguide --> xlabel[j]
                    yguide --> xlabel[i]

                    unshaped.(samples), (vsel[j], vsel[i])
                else
                    xguide --> xlabel[i]
                    yguide --> xlabel[j]

                    unshaped.(samples), (vsel[i], vsel[j])
                end
            end

            if partition_tree != false
                subspaces_rect_bounds = get_tree_par_bounds(partition_tree)
                axes = [vsel[j], vsel[i]]
                for rect in subspaces_rect_bounds
                    for ind in Base.OneTo(size(rect)[1])
                        rect[ind,:] = replace(rect[ind,:], Pair(Inf, max_samples[ind]), Pair(-Inf, min_samples[ind]))
                    end
                    box = _rect_box(rect[axes[1],1], rect[axes[1],2], rect[axes[2],1], rect[axes[2],2])

                    @series begin
                        subplot := j + (i-1)*nparams
                        seriestype --> :path
                        linewidth --> get(upper, "linewidth", 1)
                        linealpha --> get(upper, "linealpha", 1)
                        linecolor -->  get(upper, "linecolor", :darkgray)
                        legend --> false

                        xlims --> (min_samples[vsel[j]], max_samples[vsel[j]])
                        ylims --> (min_samples[vsel[i]], max_samples[vsel[i]])

                        box[1], box[2]
                    end
                end
            end

            # lower left plots
             @series begin
                 subplot := i + (j-1)*nparams

                 seriestype --> get(lower, "seriestype", :smallest_intervals)
                 bins := (bins[i], bins[j])
                 colors --> get(lower, "colors", default_colors)
                 colorbar --> get(lower, "colorbar", false)
                 intervals --> get(lower, "intervals", default_credibilities)
                 interval_labels --> get(lower, "interval_labels", [])
                 legend --> get(lower, "legend", false)
                 mean --> get(lower, "mean", mean)
                 std --> get(lower, "std", std)
                 globalmode --> get(lower, "globalmode", globalmode)
                 marginalmode --> get(lower, "marginalmode", marginalmode)
                 xlims --> get(lower, "xlims", :auto)
                 ylims --> get(lower, "ylims", :auto)
                 xguide --> xlabel[i]
                 yguide --> xlabel[j]

                 unshaped.(samples), (vsel[i], vsel[j])
             end

        end
    end

end

@recipe function f(x::Union{StepRangeLen, Vector},
        model::Function,
        sample_from::Union{DensitySampleVector};
        n_samples = 10^4,
        conf_intervals = default_credibilities,
        colors = default_colors,
        global_mode = true,
        marginal_mode = false)

    if typeof(sample_from) <: DensitySampleVector
        samples = bat_sample(sample_from, OrderedResampling(nsamples = n_samples)).result
    else
        samples = bat_sample(sample_from.prior.dist, IIDSampling(nsamples = n_samples)).result
    end

    y_ribbons = zeros(Float64, length(x), 2*length(conf_intervals))
    y_median = zeros(Float64, length(x))
    quantile_values = zeros(Float64, 2*length(conf_intervals))

    quantile_values[1:2:end] .= 0.5*(1 .- conf_intervals)
    quantile_values[2:2:end] .= 1 .- 0.5*(1 .- conf_intervals)

    for x_ind in Base.OneTo(length(x))
        y_samples = model.(samples.v, x[x_ind])
        y_median[x_ind] = quantile(y_samples, weights(samples.weight), 0.5)
        y_ribbons[x_ind,:] .= [quantile(y_samples, weights(samples.weight), quantile_tmp) for quantile_tmp in quantile_values]
        y_ribbons[x_ind,:] .= abs.(y_ribbons[x_ind,:] .- y_median[x_ind])
    end

    xguide --> "x"
    yguide --> "f(x)"
    title --> ""
    legend --> :topleft
    size --> (600, 400)

    for interval_ind in length(conf_intervals):-1:1
        @series begin
            ribbon --> (y_ribbons[:,interval_ind*2 - 1],y_ribbons[:,interval_ind*2])
            fillcolor --> colors[interval_ind]
            linecolor --> colors[interval_ind]
            seriesalpha --> 1
            linealpha --> 1
            fillalpha --> 1
            label --> "$(conf_intervals[interval_ind])"
            x, y_median
        end
    end

    @series begin
        linecolor --> :black
        linestyle --> :solid
        linewidth --> 1.5
        label --> "Median"
        x, y_median
    end

    if global_mode
        global_mode_params = mode(samples)
        @series begin
            linecolor --> :black
            linestyle --> :dot
            linewidth --> 1.5
            label --> "Global Mode"
            x, broadcast(x -> model(global_mode_params, x), x)
        end
    end

    if marginal_mode
        marginal_mode_params = bat_marginalmode(samples).result
        @series begin
            linecolor --> :black
            linestyle --> :dash
            linewidth --> 1.5
            label --> "Marginal Mode"
            x, broadcast(x -> model(marginal_mode_params, x), x)
        end
    end
end

function _rect_box(x_s::F, x_f::F, y_s::F, y_f::F) where {F<:AbstractFloat}
    x = [x_s, x_f, x_f, x_s, x_s]
    y = [y_s, y_s, y_f, y_f, y_s]
    return (x, y)
end


text:
"""
    save_localapproxvi_policy_to_jld2(policy_fn::String, policy::LocalApproximationValueIterationPolicy,
                                           mdp::Union{MDP,POMDP}, mersenne_seed::Int)

Custom saving method for storing a LocalApproximationValueIterationPolicy as a JLD2 object. It saves the individual
components based on the fields of the object. NOTE - Can only work with load_localapproxvi_policy_from_jld2
"""
function save_localapproxvi_policy_to_jld2(policy_fn::String, policy::LocalApproximationValueIterationPolicy,
                                           mdp::MDP)

    save(policy_fn, "interp", policy.interp,
                    "mdp", mdp,
                    "is_mdp_generative",policy.is_mdp_generative,
                    "n_generative_samples",policy.n_generative_samples,
                    "seed",policy.rng.seed)
end

"""
    load_localapproxvi_policy_from_jld2(policy_fn::String)

Custom loading method for retrieving a LocalApproximationValueIterationPolicy from a JLD2 file. It loads the individual
fields for the object. NOTE - Can only work with save_localapproxvi_policy_to_jld2
"""
function load_localapproxvi_policy_from_jld2(policy_fn::String, rng_type::Type{RNG} = MersenneTwister) where {RNG <: AbstractRNG}

    policy_interp = load(policy_fn, "interp")
    policy_mdp = load(policy_fn, "mdp")
    policy_isgen = load(policy_fn,"is_mdp_generative")
    policy_n_gen_samples = load(policy_fn,"n_generative_samples")
    policy_seed = load(policy_fn, "seed")

    return LocalApproximationValueIterationPolicy(policy_interp, ordered_actions(policy_mdp), policy_mdp, policy_isgen,
                                                  policy_n_gen_samples, rng_type(policy_seed))
end

function load_modal_horizon_policy_localapprox(inhor_policy_fn::String, outhor_policy_fn::String,
                                       rng_type::Type{RNG} = MersenneTwister) where {RNG <: AbstractRNG}

    inhor_policy = load_localapproxvi_policy_from_jld2(inhor_policy_fn, rng_type)
    outhor_policy = load_localapproxvi_policy_from_jld2(outhor_policy_fn, rng_type)

    @assert inhor_policy.action_map == outhor_policy.action_map "Action Maps for in and out horizon must match!"

    return ModalHorizonPolicy(inhor_policy, outhor_policy)
end

function save_modal_horizon_policy_localapprox(modal_policy::ModalHorizonPolicy, inhor_policy_fn::String, outhor_policy_fn::String,
                                       mdp::MDP)
    save_localapproxvi_policy_to_jld2(inhor_policy_fn, modal_policy.in_horizon_policy, mdp)
    save_localapproxvi_policy_to_jld2(outhor_policy_fn, modal_policy.out_horizon_policy, mdp)
end

"""
    log2space_symmetric(symm_val::Float64, n::Int64, base_val::Int64=2)

Returns a vector of logarithmically spaced numbers from -x to +x, symmetrically
spaced around 0. The number of values must be odd to reflect 0 and a symmetric
number of arguments around it.
"""
function log2space_symmetric(symm_val::Float64, n::Int64, base_val::Int64=2)

    # Ensure that the number of terms is odd (so around 0.0)
    @assert n%2==1

    # Ensure that the value itself is positive
    @assert symm_val > 0

    vals = Vector{Float64}(undef,0)
    idx = 1
    midpt = convert(Int64,round((n-1)/2))

    for i=1:midpt
        push!(vals, -(symm_val/idx))
        idx = idx*base_val
    end

    # Add 0.0 to axis
    symm_vect = reverse(-vals)
    push!(vals, 0.0)
    append!(vals, symm_vect)

    return vals
end

"""
    polyspace_symmetric(symm_val::Float64, n::Int64, exp_val::Int64=3)

Returns a vector of polynomially spaced numbers from -x to +x, symmetrically
spaced around 0. The number of values must be odd to reflect 0 and a symmetric
number of arguments around it.
"""
function polyspace_symmetric(symm_val::Float64, n::Int64, exp_val::Int64=3)

    # Ensure that the number of terms is odd (so around 0.0)
    @assert n%2==1
    # Ensure that the value itself is positive
    @assert symm_val > 0

    vals = Vector{Float64}(undef,0)
    idx = 1
    midpt = convert(Int64,round((n-1)/2))

    x = (symm_val/(midpt^exp_val))^(1/exp_val)

    for i=midpt:-1:1
        val = -1*(i*x)^exp_val
        push!(vals,val)
    end

    symm_vect = reverse(-vals)
    push!(vals, 0.0)
    append!(vals, symm_vect)

    return vals
end


text:
export iauCal2jd
"""
Gregorian Calendar to Julian Date.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   iy,im,id  int     year, month, day in Gregorian calendar (Note 1)

Returned:
   djm0      double  MJD zero-point: always 2400000.5
   djm       double  Modified Julian Date for 0 hrs

Returned (function value):
            int     status:
                        0 = OK
                        -1 = bad year   (Note 3: JD not computed)
                        -2 = bad month  (JD not computed)
                        -3 = bad day    (JD computed)

Notes:

   1. The algorithm used is valid from -4800 March 1, but this
   implementation rejects dates before -4799 January 1.

   2. The Julian Date is returned in two pieces, in the usual SOFA
   manner, which is designed to preserve time resolution.  The
   Julian Date is available as a single number by adding djm0 and
   djm.

   3. In early eras the conversion is from the "Proleptic Gregorian
   Calendar";  no account is taken of the date(s) of adoption of
   the Gregorian Calendar, nor is the AD/BC numbering convention
   observed.

Reference:

   Explanatory Supplement to the Astronomical Almanac,
   P. Kenneth Seidelmann (ed), University Science Books (1992),
   Section 12.92 (p604).

This revision:  2013 August 7

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauCal2jd(int iy, int im, int id, double *djm0, double *djm)

function iauCal2jd(iy::Real, im::Real, id::Real)

   # Allocate return value
   ref_djm0 = Ref{Float64}(0.0)
   ref_djm  = Ref{Float64}(0.0)

   status = ccall((:iauCal2jd, libsofa_c), Cint,
         (Cint, Cint, Cint,
         Ref{Cdouble}, Ref{Cdouble}),
         convert(Int32, iy), convert(Int32, im), convert(Int32, id),
         ref_djm0, ref_djm)

   # Warn on non-zero status code
   if status == -1
      @warn "bad year (Note 3: JD not computed"
   elseif status == -2
      @warn "bad month"
   elseif status == -3
      @warn "bad day"
   end

   return status, ref_djm0[], ref_djm[]
end