text:
struct SingleIntegrator <: AbstractEnv
end

"""
Even for scalar integrator, you should put an array of length 1;
due to the limitation of in-place method `Dynamics!`.
"""
function State(env::SingleIntegrator)
    return function (x::AbstractArray)
        x
    end
end

function Dynamics!(env::SingleIntegrator)
    @Loggable function dynamics!(dx, x, p, t; u)
        @assert length(dx) == length(u)
        @log integral = x
        @log integrand = u
        dx .= u
    end
end


text:
# bcf.jl
# ======
#
# A submodule for reading, writing, and working with BCF formatted files.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/GeneticVariation.jl/blob/master/LICENSE

module BCF

import BioCore: BioCore, isfilled
import GeneticVariation.VCF
import BGZFStreams
import BufferedStreams

include("record.jl")
include("reader.jl")
include("writer.jl")

end

ext: jl


text:
function ψ₁(cache, solution, k, args...)
    @↓ U, T = cache
    # if k > 1
    #     V = getU(solution[k-1])
    #     # @↓ V ← U = getU(solution[k-1])
    # else
    @↓ V ← F = cache
    # end
    r = 0.0
    N = length(U)
    for n = 1:N
        r += norm(U[n] - V[n]) / norm(U[n])
    end
    return r / N
end

function ψ₂(cache, solution, k, weights)
    @↓ U, T, F = cache
    @↓ w = weights
    w = max(1.0, w)
    r = 0.0
    N = length(U)
    for n = 1:N
        Wₙ = w^(T[1] - T[n])
        r += norm(Wₙ * (U[n] - F[n]))
    end
    return r / N
end

"""
    ErrorControl <: AbstractTimeParallelParameters

A composite type for the error control mechanism of an [`AbstractTimeParallelSolver`](@ref).

# Constructors
```julia
ErrorControl(; ϵ=1e-12, ψ=ψ₁, weights=ErrorWeights())
```

# Arguments
- `ϵ :: Real` : tolerance
- `ψ :: Function` : error function.
- `weights :: ErrorWeights` : weights for ψ.

# Functions
- [`show`](@ref) : shows name and contents.
- [`summary`](@ref) : shows name.
"""
struct ErrorControl{ϵ_T, ψ_T, weights_T} <: AbstractTimeParallelParameters
    ϵ::ϵ_T
    ψ::ψ_T
    weights::weights_T
end

ErrorControl(; ϵ=1e-12, ψ=ψ₁, weights=ErrorWeights()) = ErrorControl(ϵ, ψ, weights)


text:
@testset "data modification tests" begin
    @testset "30-bus case file incremental" begin
        data = PowerModels.parse_file("../test/data/case30.m")

        result = run_opf(data, ACPPowerModel, ipopt_solver)
        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 204.96; atol = 1e-1)

        data["branch"]["6"]["br_status"] = 0

        result = run_opf(data, ACPPowerModel, ipopt_solver)
        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 195.25; atol = 1e-1)

        data["gen"]["6"]["gen_status"] = 0

        result = run_opf(data, ACPPowerModel, ipopt_solver)
        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 195.896; atol = 1e-1)

        data["bus"]["5"]["pd"] = 0
        data["bus"]["5"]["qd"] = 0

        result = run_opf(data, ACPPowerModel, ipopt_solver)
        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 104.428; atol = 1e-1)
    end

    @testset "30-bus case file batch" begin
        data = PowerModels.parse_file("../test/data/case30.m")

        data_delta = JSON.parse("
        {
            \"branch\":{
                \"6\":{
                    \"br_status\":0
                }
            },
            \"gen\":{
                \"6\":{
                    \"gen_status\":0
                }
            },
            \"bus\":{
                \"5\":{
                    \"pd\":0,
                    \"qd\":0
                }
            }
        }")

        PowerModels.update_data(data, data_delta)

        result = run_opf(data, ACPPowerModel, ipopt_solver)
        @test result["status"] == :LocalOptimal
        @test isapprox(result["objective"], 104.428; atol = 1e-1)
    end
end


text:
function peakfinder(source::AbstractArray; σ::Real=2., threshold::Real=10., backgroundRemove::Bool=true, deconIterations::Int=3, markov::Bool=true, averWindow::Int=3)
    sigma::Float64 = σ
    PEAK_WINDOW::Int = 1024 # const of TSpectrum
    if (sigma < 1) throw(ArgumentError("Invalid sigma, must be greater than or equal to 1")) end
    if (threshold<=0 || threshold>=100) throw(ArgumentError("Invalid threshold, must be between 0 and 100")) end
    if markov if (averWindow <= 0) throw(ArgumentError("Averaging Window must be positive")) end end
    ssize = length(source)
    numberIterations = round(Int, 7*sigma + 0.5, RoundDown)
    if backgroundRemove if (ssize < 2*numberIterations + 1) throw(ArgumentError("Too large clipping window")) end end
    j = round(Int, 5*sigma + 0.5, RoundDown)
    if (j >= PEAK_WINDOW/2) throw(ArgumentError("Too large sigma")) end
    peak_index::Int=0; size_ext::Int=ssize+2*numberIterations; shift::Int=numberIterations; bw::Int=2
    plocha::Float64 = 0
    m0low, m1low, m2low, l0low, l1low = zeros(Float64, 5)

    k::Int = round(Int, 2*sigma+0.5, RoundDown)
    if k >= 2
        for i in 0:k-1
            a = i; b=source[i+1]
            m0low += 1; m1low += a; m2low += a^2; l0low += b; l1low += a*b
        end
        detlow = m0low*m2low - m1low*m1low
        l1low = (detlow != 0) ? (-l0low*m1low + l1low*m0low)/detlow : l1low = 0.
        if (l1low > 0) l1low=0. end
    else
        l1low = 0.
    end

    i::Int = 2*round(Int, 7*sigma+0.5, RoundDown)
    working_space::Array{Float64, 1} = zeros(Float64, 7*(ssize+i))
    for i in 0:size_ext-1
        if i < shift
            a = i - shift
            working_space[i+size_ext+1] = source[1]+l1low*a
            if (working_space[i+size_ext+1]<0) working_space[i+size_ext+1] = 0. end
        elseif i >= ssize+shift
            a = i - (ssize - 1 + shift)
            working_space[i+size_ext+1] = source[ssize]
            if (working_space[i+size_ext+1]<0)  working_space[i+size_ext+1] = 0. end
        else
            working_space[i+size_ext+1] = source[i-shift+1]
        end
    end

    if backgroundRemove
        for i in 1:numberIterations
            for j in i:size_ext-i-1
                if !markov
                    a = Float64(working_space[size_ext+j+1])
                    b = (working_space[size_ext+j-i+1] + working_space[size_ext+j+i+1]) / 2.0
                    if (b < a) a=b end
                    working_space[j+1]=a
                else
                    a = working_space[size_ext+j+1]
                    av = 0.
                    men = 0.
                    for w = j-bw:j+bw
                        if ((w>=0) && (w<size_ext))
                            av += working_space[size_ext+w+1]
                            men += 1
                        end
                    end
                    av = av/men
                    b = 0.
                    men = 0.
                    for w in j-i-bw:j-i+bw
                        if ((w>=0) && (w<size_ext))
                             b += working_space[size_ext+w+1]
                             men +=1
                        end
                    end
                    b = b/men
                    c = 0.
                    men = 0.
                    for w in j+i-bw:j+i+bw
                        if ((w>=0) && (w<size_ext))
                             c += working_space[size_ext+w+1]
                             men += 1
                        end
                    end
                    c = c/men
                    b = (b+c)/2.
                    if (b<a) av=b end
                    working_space[j+1]=av;
                end
            end
            for j in i:size_ext-i-1
                working_space[size_ext + j + 1] = working_space[j+1];
            end
        end

        for j in 0:size_ext-1
            if j<=shift
                a = j-shift
                b = source[1] + l1low*a
                if (b<0) b=0. end
                working_space[size_ext+j+1] = b - working_space[size_ext+j+1]
            elseif (j >= ssize + shift)
                a = j - (ssize-1+shift)
                b = source[ssize]
                if (b<0) b=0. end
                working_space[size_ext+j+1] = b - working_space[size_ext+j+1];
            else
                working_space[size_ext+j+1] = source[j-shift+1] - working_space[size_ext+j+1];
            end
        end
        for j in 0:size_ext-1
            if (working_space[size_ext+j+1] < 0) working_space[size_ext+j+1] = 0. end
        end
    end

    for i in 0:size_ext-1
        working_space[i + 6*size_ext+1] = working_space[i + size_ext + 1]
    end

    if markov
        for j in 0:size_ext-1
            working_space[2*size_ext+j+1] = working_space[size_ext + j + 1]
        end
        xmin::Int = 0
        xmax::Int = size_ext-1
        maxch::Float64 = 0.
        for i in 0:size_ext-1
            working_space[i+1] = 0.
            if (maxch < working_space[2*size_ext+i+1]) maxch = working_space[2*size_ext+i+1] end
            plocha += working_space[2*size_ext+i+1]
        end
        if (maxch == 0) error("'maxch' = 0") end

        nom::Float64 = 1.
        working_space[xmin+1]=1.
        for i in xmin:xmax-1
            nip::Float64 = working_space[2 * size_ext + i + 1] / maxch
            nim::Float64 = working_space[2 * size_ext + i + 1 + 1] / maxch
            sp::Float64 = 0.
            sm::Float64 = 0.
            for l in 1:averWindow
                a::Float64 = (i+l>xmax) ? (working_space[2 * size_ext + xmax + 1] / maxch) : (working_space[2 * size_ext + i + l + 1] / maxch)
                b::Float64 = a-nip
                a = (a+nip <= 0) ? 1 : sqrt(a+nip)

                b = exp(b/a)
                sp += b
                a = ((i-l+a)<xmin) ? (working_space[2 * size_ext + xmin + 1] / maxch) : (working_space[2 * size_ext + i - l + 1 + 1] / maxch)
                b = a-nim
                a = (a+nim<=0) ? 1 : sqrt(a+nim)
                b = exp(b/a)
                sm +=b
            end
            a = sp/sm
            a = working_space[i + 1 + 1] = working_space[i + 1] * a
            nom += a
        end
        for i in xmin:xmax
            working_space[i+1] = working_space[i+1] / nom
        end
        for j in 0:size_ext-1
            working_space[size_ext + j + 1] = working_space[j + 1] * plocha
        end
        for j in 0:size_ext-1
            working_space[2 * size_ext + j + 1] = working_space[size_ext + j + 1]
        end
        if backgroundRemove
            for i in 1:numberIterations
                for j in i:size_ext-i-1
                    a = working_space[size_ext + j + 1]
                    b = (working_space[size_ext + j - i + 1] + working_space[size_ext + j + i + 1]) / 2.0
                    if (b<a) a=b end
                    working_space[j+1] = a;
                end
                for j in i:size_ext-i-1
                    working_space[size_ext + j + 1] = working_space[j + 1]
                end
            end
            for j in 0:size_ext-1
                working_space[size_ext + j + 1] = working_space[2 * size_ext + j + 1] - working_space[size_ext + j + 1]
            end
         end
    end

    # deconvolution starts
    area::Float64 = 0
    lh_gold::Int = -1
    posit::Int = 0
    maximum::Float64 = 0
    maximum_decon::Float64 = 0
    lda::Float64 = 0
    ldb::Float64 = 0
    ldc::Float64 = 0
    # generate response vector
    for i in 0:size_ext-1
        lda = i - 3 * sigma
        lda = lda^2 / (2 * sigma^2)
        j = round(Int, 1000.0*exp(-lda), RoundDown)
        lda = j
        if (lda != 0) lh_gold = i+1. end
        working_space[i+1] = lda
        area += lda
        if lda > maximum
            maximum = lda
            posit = i
        end
    end
    # read source vector
    for i in 0:size_ext-1
        working_space[2*size_ext+i+1] = abs(working_space[size_ext + i + 1])
    end

    # create matrix at a*vector(transponse(b))
    i = lh_gold-1
    if (i > size_ext) i=size_ext end

    imin::Int = -i
    imax::Int = i
    for i in imin:imax
        lda = 0.
        jmin::Int = 0
        if (i<0) jmin = -i end
        jmax::Int = lh_gold-i
        if (jmax > lh_gold-1) jmax = lh_gold-1 end
        for j in jmin:jmax
            ldb = working_space[j+1]
            ldc = working_space[i+j+1]
            lda += ldb*ldc
        end
        working_space[size_ext + i - imin + 1] = lda;
    end

    # # create vector p
    i = lh_gold-1
    imin = -i
    imax = size_ext+i-1
    for i in imin:imax
        lda = 0.
        for j in 1:lh_gold-1
            ldb = working_space[j+1]
            k = i + j
            if ((k >= 0) && (k<size_ext))
                ldc = working_space[2*size_ext+k+1]
                lda += ldb*ldc
            end
        end
        working_space[4*size_ext+i-imin+1] = lda
    end

    # # move vector p
    for i in imin:imax
        working_space[2 * size_ext + i - imin + 1] = working_space[4 * size_ext + i - imin + 1]
    end
    # # initialization of resulting vector
    for i in 0:size_ext-1
        working_space[i+1] = 1.
    end
    # # START OF ITERATIONS
    for lindex = 0:deconIterations-1
        for i in 0:size_ext-1
            if ((abs(working_space[2*size_ext+i+1]) > 0.00001) && (abs(working_space[i+1]) > 0.00001))
                lda = 0.
                jmin = lh_gold-1
                if (jmin>i) jmin = i end
                jmin = -jmin
                jmax = lh_gold-1
                if (jmax>(size_ext-i-1)) jmax = size_ext-i-1 end
                for j in jmin:jmax
                    ldb = working_space[j+lh_gold+size_ext]
                    ldc = working_space[i+j+1]
                    lda += ldb*ldc
                end
                ldb = working_space[2*size_ext+i+1]
                if (lda!=0)
                    lda = ldb/lda
                else
                    lda = 0.
                end
                ldb = working_space[i+1]
                lda = lda*ldb
                working_space[3*size_ext+i+1] = lda
            end
        end
        for i in 0:size_ext-1
            working_space[i+1] = working_space[3*size_ext+i+1]
        end
    end

    # # shift resulting spectrum
    for i in 0:size_ext-1
        lda = working_space[i+1]
        j = i+posit
        j = j % size_ext
        working_space[size_ext+j+1] = lda
    end
    # # write back resulting spectrum
    maximum = 0.
    maximum_decon = 0.
    j = lh_gold-1
    for i in 0:size_ext-1
        if ((i >=shift) && (i < ssize+shift))
            working_space[i+1] = area*working_space[size_ext+i+j + 1]
            if (maximum_decon < working_space[i+1]) maximum_decon = working_space[i+1] end
            if (maximum < working_space[6 * size_ext + i + 1]) maximum = working_space[6 * size_ext + i + 1] end
        else
            working_space[i+1] = 0
        end
    end
    lda = 1.
    if (lda > threshold) lda = threshold end
    lda = lda/100.

    fMaxPeaks::Int = 10000
    fPosition = zeros(Float64, fMaxPeaks)
    fPositionX = zeros(Float64, fMaxPeaks)
    fPositionY = zeros(Float64, fMaxPeaks)
    fResolution = 1
    fHistogram = 0
    fNPeaks = 0
    # peak_index = 0# already defined at the start
    # # searching for peaks in deconvolved spectrum
    for i in 1:size_ext-1-1
        if ((working_space[i+1] > working_space[i]) && working_space[i+1] > working_space[i+1+1])
            if ((i >= shift) && i<ssize+shift)
                if ((working_space[i+1]>lda*maximum_decon) && (working_space[6*size_ext+i+1] > threshold*maximum / 100.))
                    a = 0.
                    b = 0.
                    for j in i-1:i+1
                        a += (j-shift)*working_space[j+1]
                        b += working_space[j+1]
                    end
                    a = a/b
                    if (a<0) a = 0. end
                    if (a>=ssize) a=ssize-1 end
                    if (peak_index == 0)
                        fPositionX[0+1] = a
                        peak_index = 1
                    else
                        priz = 0
                        j = 0
                        while ((priz == 0) && (j<peak_index))
                        # for j in 0:peak_index-1
                            if (working_space[6*size_ext+shift+round(Int, a, RoundDown)+1] > working_space[6*size_ext+shift+round(Int, fPositionX[j+1], RoundDown)+1])
                                priz = 1
                            end
                            j+=1
                        end
                        if (priz == 0)
                            if (j<fMaxPeaks) fPositionX[j+1] = a end
                        else
                            for k in peak_index:-1:j
                                if (k<fMaxPeaks) fPositionX[k+1] = fPositionX[k-1+1] end
                            end
                            fPositionX[j] = a
                        end
                        if (peak_index < fMaxPeaks) peak_index += 1 end
                    end
                end
            end
        end
    end

    # return working_space
    destVector = zeros(Float64, ssize)
    for i in 0:ssize-1
        destVector[i+1] = working_space[i + shift + 1]
    end
    fNPeaks = peak_index
    if (peak_index == fMaxPeaks) warn("Peak buffer full", prefix="Warning: SearchHighRes -> ") end

    return destVector, fPositionX[1:fNPeaks]
end

"""
    peakfinder(h::Histogram; <keyword arguments>)::Tuple{Histogram, Array{Float64, 1}}

Returns a deconvoluted spectrum and an array of peak positions.

# Keywords
- `σ::Real=2.0`: The expected sigma of a peak in the spectrum. In units of bins.
- `threshold::Real=10.0`: Threshold for being identified as a peak in the deconvoluted spectrum. A single bin is identified as a peak when its weight exceeds the `threshold` and the previous bin was not identified as an peak.
- `backgroundRemove::Bool=true`
- `deconIterations::Int=3`
- `markov::Bool=true`
- `averWindow::Int=3`

# Source
This function is basically a copy of `TSpectrum::SearchHighRes` from [ROOT](https://root.cern.ch/doc/master/classTSpectrum.html).

- M.A. Mariscotti: A method for identification of peaks in the presence of background and its application to spectrum analysis. NIM 50 (1967), 309-320.
- M. Morhac;, J. Kliman, V. Matouoek, M. Veselsky, I. Turzo.:Identification of peaks in multidimensional coincidence gamma-ray spectra. NIM, A443 (2000) 108-125.
- Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376 (1996), 451.
"""
function peakfinder(h::Histogram; kw...)::Tuple{Histogram, Array{Float64, 1}}
    destVector, fPositionX = peakfinder(h.weights; kw...)
    rh = Histogram(h.edges[1], destVector, :left)
    rh.weights = destVector

    mp = StatsBase.midpoints(h.edges[1]);
    fPositionX_idcs = Int.(round.(1 .+ fPositionX, RoundUp));
    weight_in_bin = (fPositionX .+ 2) .% fPositionX_idcs .- 0.5
    binvolumes = [StatsBase.binvolume(h, i) for i in fPositionX_idcs];

    fPositionX = mp[fPositionX_idcs] + binvolumes .* weight_in_bin - 0.5 * binvolumes

    return rh, fPositionX
end


text:
result = [2.32551,
-0.155293,
0.469214,
-0.0870457,
0.271364,
-0.222398,
1.01423,
-0.179009,
1.01724,
-0.152972,
0.216039,
-0.251637,
1.05034,
-0.231289,
0.245388,
-0.231289,
1.03371,
-0.258872,
1.03256,
-0.262519,
1.04748,
-0.259143,
1.0535,
-0.266484,
1.04711,
-0.267177,
1.02131,
-0.280381]

include(joinpath(BASE_DIR, "test", "data_5bus_pu.jl"))
include(joinpath(BASE_DIR, "test", "data_14bus_pu.jl"))

nodes_5 = nodes5()
nodes_14 = nodes14()
c_sys14 = System(nodes_14, thermal_generators14(nodes_14), loads14(nodes_14), branches14(nodes_14), nothing, 100.0, nothing, nothing);

c_sys5_re = System(nodes_5, vcat(thermal_generators5(nodes_5), renewable_generators5(nodes_5)), loads5(nodes_5),
                nothing, nothing,  100.0, nothing, nothing)

@testset begin
    using NLsolve
    pf!, x0 = make_pf(c_sys14);
    res = NLsolve.nlsolve(pf!, x0)
    @test res.zero ≈ result rtol=1e-3

    solve_powerflow!(c_sys14, nlsolve, method = :newton)
    pf!, x0 = make_pf(c_sys5_re);
    res = NLsolve.nlsolve(pf!, x0)
    @test !res.f_converged
end


text:
"""
    renyi(p, q)

    Computes the Renyi heterogeneity of order q for the vector p
"""
function renyi(p::Vector{Float64}, q::Union{Float64, Int64}=1)
    p = filter(x -> x > 0, p)
    p = p ./ sum(p)
    if q == 0
        Π = length(p)
    elseif q == 1
        Π = exp.(-p'*log.(p))
    elseif q == Inf
        Π = 1/maximum(p)
    else
        Π = sum(p.^q)^(1/(1-q))
    end
    return Π
end

"""
    renyi_pooled(P, W, q)

    Computes the pooled Renyi heterogeneity of order q for the NxM matrix P, given the N-dimensional
    weight vector W.
"""
function renyi_pooled(P::Matrix{Float64}, W::Vector{Float64}, q::Union{Float64, Int64}=1)
    return renyi(P'*W, q)
end

"""
    renyi_within(P, W, q)

    Computes the within-sample Renyi heterogeneity of order q for the NxM matrix P, given the
    N-dimensional weight vector W.
"""
function renyi_within(P::Matrix{Float64}, W::Vector{Float64}, q::Union{Float64, Int64}=1)
    P = P ./ repeat(sum(P, dims=2), 1, size(P)[2])  # Normalize rows
    if q == 1
        out = exp(-W'*map(p -> filter(x-> x > 0, p[:]) |> a -> a'*log.(a), eachrow(P)))
    elseif (q > 0) & (q < Inf) & (q != 1)
        Wq = W.^q
        w = Wq ./ sum(Wq)
        p  = sum(P.^q, dims=2)[:]
        out = (w'*p)^(1/(1-q))
    end
    return out
end

"""
    renyi_between(P, W, q)

    Computes the within-sample Renyi heterogeneity of order q for the NxM matrix P, given the
    N-dimensional weight vector W.
"""
function renyi_between(P::Matrix{Float64}, W::Vector{Float64}, q::Union{Float64, Int64}=1)
    return renyi_pooled(P,W,q)/renyi_within(P,W,q)
end


text:
"""
    struct RegularPrism{T,CO,N,TR} <: AbstractVolumePrimitive{T, CO}

Volume primitive describing a [Prism](@ref) with base plates are regular polygons
which are parallel to the `xy` plane. If the regular polygon base plate is projected to
the `xy` plane, one of the vertices lays on the `x` axis.

## Parametric types
* `T`: Precision type.
* `CO`: Describes whether the surface belongs to the primitive.
    It can be `ClosedPrimitive`, i.e. the surface points belong to the primitive,
    or `OpenPrimitive`, i.e. the surface points do not belong to the primitive.
* `N`: Number of vertices of the regular polygon that defines the base of the prism.
* `TR`: Type of `r`.
    * `TR == T`: Regular polygon base (all vertices have the same distance to the center).

## Fields
* `r::TR`: Distance of the vertices to the center of the regular polygon base (in m).
* `hZ::T`: Half of the width in `z` dimension (in m).
* `origin::CartesianPoint{T}`: The position of the center of the `RegularPrism`.
* `rotation::SMatrix{3,3,T,9}`: Matrix that describes a rotation of the `RegularPrism` around its `origin`.

## Definition in Configuration File

So far, only `HexagonalPrism` can be defined in the configuration files.
A `HexagonalPrism` is defined in the configuration file as part of the `geometry` field
of an object through the field `HexagonalPrism`.

An example definition of a `HexagonalPrism` looks like this:
```yaml
HexagonalPrism:
  r: 1.0 # => r = 1.0
  h: 2.0 # => hZ = 1.0
```

See also [Constructive Solid Geometry (CSG)](@ref).
"""
@with_kw struct RegularPrism{T,CO,N,TR} <: AbstractVolumePrimitive{T, CO}
    r::TR = 1
    hZ::T = 1

    origin::CartesianPoint{T} = zero(CartesianPoint{T})
    rotation::SMatrix{3,3,T,9} = one(SMatrix{3, 3, T, 9})
end

RegularPrism{T, CO, N, TR}( rp::RegularPrism{T, CO, N, TR}; COT = CO,
            origin::CartesianPoint{T} = rp.origin,
            rotation::SMatrix{3,3,T,9} = rp.rotation) where {T, CO<:Union{ClosedPrimitive, OpenPrimitive}, N, TR} =
    RegularPrism{T, COT, N, TR}(rp.r, rp.hZ, origin, rotation)

const TriangularPrism{T,CO,TR} = RegularPrism{T,CO,3,TR}
const QuadranglePrism{T,CO,TR} = RegularPrism{T,CO,4,TR}
const PentagonalPrism{T,CO,TR} = RegularPrism{T,CO,5,TR}
const HexagonalPrism{T,CO,TR}  = RegularPrism{T,CO,6,TR}

function Geometry(::Type{T}, ::Type{P}, dict::AbstractDict, input_units::NamedTuple, transformations::Transformations{T}
            ) where {T, P <: Union{TriangularPrism, QuadranglePrism, PentagonalPrism, HexagonalPrism}}
    length_unit = input_units.length
    angle_unit = input_units.angle
    origin = get_origin(T, dict, length_unit)
    rotation = get_rotation(T, dict, angle_unit)

    r = parse_r_of_primitive(T, dict, length_unit)
    hZ = if haskey(dict, "h")
        _parse_value(T, dict["h"], length_unit) / 2
    elseif haskey(dict, "z")
        z = parse_height_of_primitive(T, dict, length_unit)
        if z isa Real
            @warn "Deprecation warning: Field `z` for `RegularPrism` is deprecated.
                Use `h` instead to specify the height of the primitive."
            z
        else # z isa Tuple
            @warn "Deprecation warning: Field `z` for `RegularPrism` is deprecated.
                Use `h` instead to specify the height of the primitive.
                There might be a conflict with the possible field `origin`:
                The `z` component of the origin of the primitive is overwritten by the `z`."
            origin = CartesianPoint{T}(origin[1], origin[2], mean(z))
            (z[2] - z[1])/2
       end
    end

    g = if r isa Tuple # lazy workaround for now
        P{T,ClosedPrimitive,T}(r = r[2], hZ = hZ, origin = origin, rotation = rotation) -
        P{T,ClosedPrimitive,T}(r = r[1], hZ = T(1.1)*hZ, origin = origin, rotation = rotation)
    else
        P{T,ClosedPrimitive,T}(r = r, hZ = hZ, origin = origin, rotation = rotation)
    end
    transform(g, transformations)
end

const PrismAliases = Dict{Int, String}(
    3 => "TriangularPrism",
    4 => "QuadranglePrism",
    5 => "PentagonalPrism",
    6 => "HexagonalPrism"
)

function Dictionary(rp::RegularPrism{T, <:Any, N})::OrderedDict{String, Any} where {T, N}
    dict = OrderedDict{String, Any}()
    dict["r"] = rp.r # always a Real
    dict["h"] = rp.hZ*2
    if rp.origin != zero(CartesianVector{T}) dict["origin"] = rp.origin end
    if rp.rotation != one(SMatrix{3,3,T,9}) dict["rotation"] = Dictionary(rp.rotation) end
    OrderedDict{String, Any}(PrismAliases[N] => dict)
end

function vertices(rp::RegularPrism{T,ClosedPrimitive,N,T}) where {T,N}
    xys = [rp.r .* sincos(T(2π)*(n-1)/N) for n in 1:N]
    pts = [CartesianPoint{T}(xy[2], xy[1], z) for z in (-rp.hZ, rp.hZ) for xy in xys]
    _transform_into_global_coordinate_system(pts, rp)
end
function vertices(rp::RegularPrism{T,OpenPrimitive,N,T}) where {T,N}
    xys = [rp.r .* sincos(T(2π)*(n-1)/N) for n in N:-1:1]
    pts = [CartesianPoint{T}(xy[2], xy[1], z) for z in (-rp.hZ, rp.hZ) for xy in xys]
    _transform_into_global_coordinate_system(pts, rp)
end

function surfaces(rp::RegularPrism{T,<:Any,N,T}) where {T,N}
    vs = (vertices(rp))
    p_bot = Polygon{N,T}(vs[1:N])
    p_top = Polygon{N,T}(reverse(vs[N+1:end]))
    quads = Vector{Quadrangle{T}}(undef, N)
    for i in 1:N-1
        quads[i] = Quadrangle{T}((vs[i], vs[N+i], vs[N+i+1], vs[i+1]))
    end
    quads[N] = Quadrangle{T}((vs[N], vs[2N], vs[N+1], vs[1]))
    p_bot, p_top, quads...
end

function sample(rp::RegularPrism{T,<:Any,N,T})::Vector{CartesianPoint{T}} where {T,N}
    [vertices(rp)...]
end

function _in(pt::CartesianPoint{T}, rp::RegularPrism{T,ClosedPrimitive,N,T}; csgtol::T = csg_default_tol(T)) where {T,N}
    abs(pt.z) <= rp.hZ + csgtol && begin
        r, φ = hypot(pt.x, pt.y), atan(pt.y, pt.x)
        α = T(π/N)
        _r = r * cos(α - mod(φ, 2α)) / cos(α)
        _r <= rp.r + csgtol
    end
end
function _in(pt::CartesianPoint{T}, rp::RegularPrism{T,OpenPrimitive,N,T}; csgtol::T = csg_default_tol(T)) where {T,N}
    abs(pt.z) < rp.hZ - csgtol && begin
        r, φ = hypot(pt.x, pt.y), atan(pt.y, pt.x)
        α = T(π/N)
        _r = r * cos(α - mod(φ, 2α)) / cos(α)
        _r < rp.r - csgtol
    end
end

# # Convenience functions
# const TriangularPrism{T,TR,TZ} = RegularPrism{3,T,TR,TZ}
# const SquarePrism{T,TR,TZ}  = RegularPrism{4,T,TR,TZ}
# const PentagonalPrism{T,TR,TZ} = RegularPrism{5,T,TR,TZ}
# const HexagonalPrism{T,TR,TZ}  = RegularPrism{6,T,TR,TZ}

# TriangularPrism(args...) = RegularPrism(3, args...)
# SquarePrism(args...)  = RegularPrism(4, args...)
# PentagonalPrism(args...) = RegularPrism(5, args...)
# HexagonalPrism(args...)  = RegularPrism(6, args...)

# print(io::IO, rp::TriangularPrism{T, TR, TZ}) where {T,TR,TZ} = print(io, "TriangularPrism{$(T), $(TR), $(TZ)}($(rp.r), $(rp.z))")
# print(io::IO, rp::SquarePrism{T, TR, TZ})  where {T,TR,TZ} = print(io, "SquarePrism{$(T), $(TR), $(TZ)}($(rp.r), $(rp.z))")
# print(io::IO, rp::PentagonalPrism{T, TR, TZ}) where {T,TR,TZ} = print(io, "PentagonalPrism{$(T), $(TR), $(TZ)}($(rp.r), $(rp.z))")
# print(io::IO, rp::HexagonalPrism{T, TR, TZ})  where {T,TR,TZ} = print(io, "HexagonalPrism{$(T), $(TR), $(TZ)}($(rp.r), $(rp.z))")


text:
const VectorLPSet = Union{MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives}

scalar_set_type(::Type{<:MOI.Zeros}, T::Type) = MOI.EqualTo{T}
scalar_set_type(::Type{<:MOI.Nonpositives}, T::Type) = MOI.LessThan{T}
scalar_set_type(::Type{<:MOI.Nonnegatives}, T::Type) = MOI.GreaterThan{T}

__constant(f::Union{MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction}, T::DataType) = MOI._constant(f)
__constant(f::Union{MOI.VectorAffineFunction, MOI.VectorQuadraticFunction}, T::DataType) = MOI._constant(f)
__constant(f::MOI.SingleVariable, T::DataType) = zero(T)
__constant(f::MOI.VectorOfVariables, T::DataType) = zeros(T, length(f.variables))

"""
    ScalarizeBridge{T}

Transforms a constraint `AbstractVectorFunction`-in-`S` where `S <: LPCone` to
`AbstactVectorFunction`-in-`scalar_set_type(S)`.
"""
mutable struct ScalarizeBridge{T, F, S} <: AbstractBridge
    scalar_constraints::Vector{CI{F, S}}
    constants::Vector{T}
end
function ScalarizeBridge{T, F, S}(model::MOI.ModelLike,
                                  f::MOI.AbstractVectorFunction,
                                  set::VectorLPSet) where {T, F, S}
    dimension = MOI.output_dimension(f)
    constants = __constant(f, T)
    new_f = MOIU.scalarize(f, true)
    constraints = Vector{CI{F, S}}(undef, dimension)
    for i in 1:dimension
        constraints[i] = MOIU.add_scalar_constraint(model, new_f[i], S(-constants[i]))
    end
    return ScalarizeBridge{T, F, S}(constraints, constants)
end

function MOI.supports_constraint(::Type{ScalarizeBridge{T}},
                                 ::Type{<:MOI.AbstractVectorFunction},
                                 ::Type{<:VectorLPSet}) where T
    return true
end
function added_constraint_types(::Type{ScalarizeBridge{T, F, S}}) where {T, F, S}
    return [(F, S)]
end
function concrete_bridge_type(::Type{<:ScalarizeBridge{T}},
                              F::Type{<:MOI.AbstractVectorFunction},
                              S::Type{<:VectorLPSet}) where T
    return ScalarizeBridge{T, MOIU.scalar_type(F), scalar_set_type(S, T)}
end

# Attributes, Bridge acting as an model
function MOI.get(bridge::ScalarizeBridge{T, F, S},
                 ::MOI.NumberOfConstraints{F, S}) where {T, F, S}
    return length(bridge.scalar_constraints)
end
function MOI.get(bridge::ScalarizeBridge{T, F, S},
                 ::MOI.ListOfConstraintIndices{F, S}) where {T, F, S}
    return bridge.scalar_constraints
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::ScalarizeBridge)
    MOI.delete.(model, bridge.scalar_constraints)
    nothing
end

# Attributes, Bridge acting as a constraint

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal,
                 bridge::ScalarizeBridge)
    return MOI.get.(model, attr, bridge.scalar_constraints) .+ bridge.constants
end
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::ScalarizeBridge)
    return MOI.get.(model, attr, bridge.scalar_constraints)
end
function MOI.modify(model::MOI.ModelLike, bridge::ScalarizeBridge{T,F,S},
    change::MOI.VectorConstantChange{T}) where {T,F,S}
    bridge.constants = change.new_constant
    MOI.set.(model, MOI.ConstraintSet(), bridge.scalar_constraints,
             S.(-change.new_constant))
end
function MOI.modify(model::MOI.ModelLike, bridge::ScalarizeBridge,
                    change::MOI.MultirowChange{T}) where T
    for (index, value) in change.new_coefficients
        MOI.modify(model, bridge.scalar_constraints[index],
               MOI.ScalarCoefficientChange{T}(change.variable, value))
    end
    nothing
end
function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintFunction,
    bridge::ScalarizeBridge{T}, func) where T
    old_constants = bridge.constants
    bridge.constants = __constant(func, T)
    new_func = MOIU.scalarize(func, true)
    MOI.set.(model, MOI.ConstraintFunction(), bridge.scalar_constraints,
             new_func)
    for i in eachindex(bridge.constants)
        if bridge.constants[i] != old_constants[i]
            MOI.set(model, MOI.ConstraintSet(), bridge.scalar_constraints,
                     S(-bridge.constants[i]))
        end
    end
end

# TODO implement transform


text:
"""
    like(state, agent_idx, config)

Generate a like for a post.

# Arguments
- `state`: A tuple of the current graph and agent_list
- `agent_idx`: Agent index
- `config`: Config object as provided by Config

See also: [`Config`](@ref), [`Agent`](@ref)
"""
function like!(
    state::Tuple{AbstractGraph, AbstractArray},
    agent_idx::Integer,
    config::Config
)

    agent_list = state[2]
    this_agent = agent_list[agent_idx]

    for post in this_agent.feed
        if (
            (
                abs(post.opinion - this_agent.opinion)
                < config.opinion_threshs.like
            )
            && !(post in this_agent.liked_posts)
            && !(post in this_agent.disliked_posts)
        )
            post.like_count += 1
            post.weight *= 1.01
            push!(this_agent.liked_posts, post)
        end
    end

    return state

end

# suppress output of include()
;


text:
@testset "NeuralNetworkQ" begin
    q = NeuralNetworkQ(;model=Dense(2, 2, initW = ones, initb = zeros), optimizer = Flux.Optimise.Descent(0.1))
    x = [1, 2]
    y = [1, 0]
    gs = Flux.gradient(q.params) do
        Flux.crossentropy(RL.batch_estimate(q, x), y)
    end
    update!(q, gs)

    @test q(x) == [3.2, 3.0]
end


text:
StandaloneLightningEigenfunction(
    domain::Triangle{T},
    i::Integer;
    outside = false,
    l = 1,
    σ = 4,
    even::Bool = false,
    odd::Bool = false,
    reversed::Bool = false,
) where {T<:Union{arb,fmpq}} =
    StandaloneLightningEigenfunction{arb,T}(domain, i; outside, l, σ, even, odd, reversed)

function StandaloneLightningEigenfunction{T,S}(
    domain::Triangle{S},
    i::Integer;
    outside::Bool = false,
    l = 1,
    σ = 4,
    even::Bool = false,
    odd::Bool = false,
    reversed::Bool = false,
) where {T,S}
    θ = domain.angles[i]

    if i == 1
        orientation = zero(θ)
    elseif i == 2
        orientation = ifelse(S == arb, domain.parent(π), 1) - θ
    elseif i == 3
        orientation = 2ifelse(S == arb, domain.parent(π), 1) - domain.angles[2] - θ
    end

    if outside
        orientation += θ
        θ = ifelse(T == arb, 2domain.parent(π), 2) - θ
    end

    return StandaloneLightningEigenfunction{T,S}(
        vertex(domain, i),
        orientation,
        θ;
        l,
        σ,
        even,
        odd,
        reversed,
    )
end

# TODO: Currently the orientation is computed from the placement of
# the vertices and will therefore always be an arb and never and fmpq.
# The value for θ is also converted to an arb for that reason.
function StandaloneLightningEigenfunction(
    domain::Polygon{T},
    i::Integer;
    outside = false,
    l::arb = domain.parent(1),
    σ::arb = domain.parent(4),
    even::Bool = false,
    odd::Bool = false,
    reversed::Bool = false,
) where {T<:Union{arb,fmpq}}
    θ = angle(domain, i)

    v = vertex(domain, mod1(i + 1, length(vertices(domain)))) - vertex(domain, i)
    orientation = atan(v[2], v[1])

    if contains_zero(v[2])
        @warn "orientation could not be computed accurately, orientation = $orientation"
    end

    if outside
        orientation += θ
        θ = 2domain.parent(π) - θ
    end

    return StandaloneLightningEigenfunction(
        vertex(domain, i),
        orientation,
        θ,
        domain.parent;
        l,
        σ,
        even,
        odd,
        reversed,
    )
end

function StandaloneLightningEigenfunction{T,arb}(
    domain::Polygon,
    i::Integer;
    outside = false,
    l = 1,
    σ = 4,
    even::Bool = false,
    odd::Bool = false,
    reversed::Bool = false,
) where {T}
    θ = angle(domain, i)

    v = vertex(domain, mod1(i + 1, length(vertices(domain)))) - vertex(domain, i)
    orientation = atan(v[2], v[1])

    if contains_zero(v[2])
        @warn "orientation could not be computed accurately, orientation = $orientation"
    end

    if outside
        orientation += θ
        θ = 2domain.parent(π) - θ
    end

    return StandaloneLightningEigenfunction{T,arb}(
        vertex(domain, i),
        orientation,
        θ;
        l,
        σ,
        even,
        odd,
        reversed,
    )
end

function StandaloneLightningEigenfunction(
    domain::TransformedDomain,
    i::Integer;
    outside = false,
    l::arb = domain.parent(1),
    σ::arb = domain.parent(4),
    even::Bool = false,
    odd::Bool = false,
    reversed::Bool = false,
)
    u = StandaloneLightningEigenfunction(
        domain.original,
        i;
        outside,
        l,
        σ,
        even,
        odd,
        reversed,
    )
    if typeof(u.orientation) == typeof(domain.rotation)
        orientation = u.orientation + domain.rotation
    else
        if u.orientation isa fmpq
            orientation = domain.parent(π) * u.orientation + domain.rotation
        else
            orientation = u.orientation + domain.parent(π) * domain.rotation
        end
    end

    u = StandaloneLightningEigenfunction(
        domain.map(u.vertex),
        orientation,
        u.θ,
        domain.parent,
        l = u.l,
        σ = u.σ,
        even = u.even,
        odd = u.odd,
        reversed = u.reversed,
    )

    return u
end

function StandaloneLightningEigenfunction{T,S}(
    domain::TransformedDomain,
    i::Integer;
    outside = false,
    l = 1,
    σ = 4,
    even::Bool = false,
    odd::Bool = false,
    reversed::Bool = false,
) where {T,S}
    u = StandaloneLightningEigenfunction{T,S}(
        domain.original,
        i;
        outside,
        l,
        σ,
        even,
        odd,
        reversed,
    )
    if typeof(u.orientation) == typeof(domain.rotation)
        orientation = u.orientation + domain.rotation
    else
        if u.orientation isa fmpq
            orientation = domain.parent(π) * u.orientation + domain.rotation
        else
            orientation = u.orientation + domain.parent(π) * domain.rotation
        end
    end

    u = StandaloneLightningEigenfunction{T,S}(
        domain.map(u.vertex),
        orientation,
        u.θ;
        u.l,
        u.σ,
        u.even,
        u.odd,
        u.reversed,
    )

    return u
end

function Base.show(io::IO, u::StandaloneLightningEigenfunction{T}) where {T}
    println(
        io,
        "StandaloneLightningEigenfunction{$T}" *
        ifelse(u.even, " - even", "") *
        ifelse(u.odd, " - odd", "") *
        ifelse(u.reversed, " - reversed", ""),
    )
    if !haskey(io, :compact) || !io[:compact]
        println(io, "vertex: $(u.vertex)")
        println(io, "orientation: $(u.orientation), θ: $(u.θ)")
        println(io, "l = $(u.l), σ = $(u.σ)")
        print(io, "number of set coefficients: $(length(u.coefficients))")
    end
end

function set_eigenfunction!(u::StandaloneLightningEigenfunction{arb}, coefficients::Vector)
    resize!(u.coefficients, length(coefficients))
    copy!(u.coefficients, u.parent.(coefficients))
end

function set_eigenfunction!(
    u::StandaloneLightningEigenfunction{T},
    coefficients::Vector,
) where {T}
    resize!(u.coefficients, length(coefficients))
    copy!(u.coefficients, coefficients)
end

"""
    coordinate_transformation(u::StandaloneLightningEigenfunction, xy::AbstractVector)

Takes a 2-element vector `xy` representing a point in the plane in
Cartesian coordinates and makes a (affine) change of coordinates so
that `u.vertex` is put at the origin and rotated `u.orientation + u.θ/2`
clockwise.
"""
function coordinate_transformation(
    u::StandaloneLightningEigenfunction{T,S},
    xy::AbstractVector,
) where {T,S}
    if S == fmpq
        s, c = sincospi(-u.orientation - u.θ // 2, u.parent)
    elseif S == arb
        s, c = sincos(-u.orientation - u.θ / 2)
    end

    s, c = convert.(T, (s, c))

    if eltype(xy) == arb_series
        xy = convert(SVector{2,arb_series}, xy)
    else
        xy = convert(SVector{2,T}, xy)
    end

    M = SMatrix{2,2}(c, s, -s, c)
    res = M * (xy .- u.vertex)
    if u.reversed
        return SVector(res[1], -res[2])
    else
        return res
    end
end

"""
    chargedistance(u::StandaloneLightningEigenfunction, i::Integer, n::Integer)

Compute the distance from the vertex for the `i`th charge point
assuming a total of `n` charge points are used.

We use a tapered distribution as proposed by Trefethen but we reverse
the order.

An old alternative version is `u.l*exp(-u.σ*(i - 1)/sqrt(u.parent(n)))`.
"""
chargedistance(u::StandaloneLightningEigenfunction{arb}, i::Integer, n::Integer) =
    u.l * exp(-u.σ * (sqrt(u.parent(n)) - sqrt(u.parent(n + 1 - i))))

chargedistance(u::StandaloneLightningEigenfunction{T}, i::Integer, n::Integer) where {T} =
    u.l * exp(-u.σ * (sqrt(convert(T, n)) - sqrt(convert(T, n + 1 - i))))

"""
    charge(
        u::StandaloneLightningEigenfunction,
        i::Integer,
        n::Integer,
        standard_coordinates = false,
)

Compute the location of the `i`th charge point assuming a total of `n`
points is used. The charge points are numbered from `1` to `n`.

They are returned in the coordinate system used by `u`, so they all
lie on the negative part of the x-axis. If `standard_coordinates` is
true then return them in the standard coordinate system instead.
"""
function charge(
    u::StandaloneLightningEigenfunction{T,S},
    i::Integer,
    n::Integer,
    standard_coordinates = false,
) where {T,S<:Union{arb,fmpq}}
    d = chargedistance(u, i, n)
    if standard_coordinates
        if S == fmpq
            s, c = sincospi(u.orientation + u.θ // 2, u.parent)
        else
            s, c = sincos(u.orientation + u.θ / 2)
        end

        return u.vertex - d .* SVector{2,T}(c, s)
    else
        SVector{2}(-d, zero(d))
    end
end

function (u::StandaloneLightningEigenfunction{T,S})(
    xy::AbstractVector,
    λ::Union{arb,Real},
    ks::UnitRange{Int};
    boundary = nothing,
    notransform::Bool = false,
) where {T,S}
    if T == arb && eltype(xy) == arb_series
        xy = SVector{2,arb_series}(xy[1], xy[2])
        λ = u.parent(λ)
    elseif T == arb
        xy = SVector{2,arb}(u.parent(xy[1]), u.parent(xy[2]))
        λ = u.parent(λ)
    else
        xy = convert(SVector{2,T}, xy)
        λ = convert(T, λ)
    end

    if !notransform
        xy = coordinate_transformation(u, xy)
    end

    if T == arb && eltype(xy) == arb_series
        res = similar(ks, arb_series)
    else
        res = similar(ks, T)
    end
    i = 1
    if u.even
        charge_indices = (div(ks.start - 1, 2)+1):(div(ks.stop - 1, 2)+1)
    elseif u.odd
        charge_indices = ks
    else
        charge_indices = (div(ks.start - 1, 3)+1):(div(ks.stop - 1, 3)+1)
    end
    for charge_index in charge_indices
        # Perform the change of coordinates corresponding to the charge
        # point
        xy_local = xy - charge(u, charge_index, charge_indices.stop)
        r, θ = polar_from_cartesian(xy_local)

        rsqrtλ = r * sqrt(λ)
        b = bessel_y(one(λ), rsqrtλ)
        s, c = sincos(θ)

        if u.even
            if charge_index == charge_indices.start == charge_indices.stop
                # Special case when first and last charge are the same
                if mod1(ks.start, 2) == 1
                    res[i] = bessel_y(zero(λ), rsqrtλ)
                    i += 1
                end
                if mod1(ks.stop, 2) == 2
                    res[i] = b * c
                    i += 1
                end
            elseif charge_index == charge_indices.start
                # Might not want all terms from the first charge
                if mod1(ks.start, 2) == 1
                    res[i] = bessel_y(zero(λ), rsqrtλ)
                    i += 1
                end
                res[i] = b * c
                i += 1
            elseif charge_index == charge_indices.stop
                # Might not want all terms from the last charge
                res[i] = bessel_y(zero(λ), rsqrtλ)
                i += 1
                if mod1(ks.stop, 2) == 2
                    res[i] = b * c
                    i += 1
                end
            else
                res[i] = bessel_y(zero(λ), rsqrtλ)
                res[i+1] = b * c
                i += 2
            end
        elseif u.odd
            # We only take one expansion from each charge in this
            # case. So it's a lot simpler.
            res[i] = b * s
            i += 1
        else
            if charge_index == charge_indices.start == charge_indices.stop
                # Special case when first and last charge are the same
                if mod1(ks.start, 3) == 1
                    res[i] = bessel_y(zero(λ), rsqrtλ)
                    i += 1
                end
                if mod1(ks.start, 3) <= 2 && mod1(ks.stop, 3) >= 2
                    res[i] = b * s
                    i += 1
                end
                if mod1(ks.stop, 3) == 3
                    res[i] = b * c
                    i += 1
                end
            elseif charge_index == charge_indices.start
                # Might not want all terms from the first charge
                if mod1(ks.start, 3) == 1
                    res[i] = bessel_y(zero(λ), rsqrtλ)
                    i += 1
                end
                if mod1(ks.start, 3) <= 2
                    res[i] = b * s
                    i += 1
                end
                res[i] = b * c
                i += 1
            elseif charge_index == charge_indices.stop
                # Might not want all terms from the last charge
                res[i] = bessel_y(zero(λ), rsqrtλ)
                i += 1
                if mod1(ks.stop, 3) >= 2
                    res[i] = b * s
                    i += 1
                end
                if mod1(ks.stop, 3) >= 3
                    res[i] = b * c
                    i += 1
                end
            else
                res[i] = bessel_y(zero(λ), rsqrtλ)
                res[i+1] = b * s
                res[i+2] = b * c
                i += 3
            end
        end
    end

    return res
end

# TODO: Figure out how to handle this
set_domain!(u::StandaloneLightningEigenfunction, ::AbstractDomain) = u

# TODO: Figure out how to handle this
recompute!(u::StandaloneLightningEigenfunction) = u


text:
module LegolasFlux

export write_model_row, read_model_row
export fetch_weights, load_weights!

using Legolas
using Arrow
using Arrow.ArrowTypes
using Tables
using Functors
using Base: IdSet

const LEGOLAS_SCHEMA = Legolas.Schema("legolas-flux.model@1")

#####
##### `FlatArray`
#####

struct FlatArray{T <: Number}
    vec::Vector{T}
    size::Vector{Int}
end

FlatArray{T}(x::AbstractArray{T}) where {T} = FlatArray{T}(vec(x), collect(size(x)))
FlatArray{T}(x::AbstractArray{S}) where {T, S} = FlatArray{T}(convert(Vector{T}, vec(x)), collect(size(x)))
FlatArray(x::AbstractArray{T}) where {T} = FlatArray{T}(vec(x), collect(size(x)))
Base.Array(x::FlatArray) = reshape(x.vec, x.size...)

for op in (:(==), :isequal)
    @eval Base.$(op)(f1::FlatArray, f2::FlatArray) = $op(f1.size, f2.size) && $op(f1.vec, f2.vec)
end
Base.hash(f::FlatArray, h::UInt) = hash(:FlatArray, hash(f.vec, hash(f.size, h)))

const FLATARRAY_ARROW_NAME = Symbol("JuliaLang.LegolasFlux.FlatArray")
ArrowTypes.arrowname(::Type{<:FlatArray}) = FLATARRAY_ARROW_NAME
ArrowTypes.JuliaType(::Val{FLATARRAY_ARROW_NAME}) = FlatArray

#####
##### `Weights`
#####

struct Weights{T}
    weights::Vector{FlatArray{T}}
end
Base.hash(w::Weights, h::UInt) = hash(:Weights, hash(w.weights, h))

for op in (:(==), :isequal)
    @eval Base.$(op)(w1::Weights, w2::Weights) = $op(w1.weights, w2.weights)
end

Weights(w::Weights) = w
Base.length(w::Weights) = length(w.weights)
Base.eltype(w::Weights{T}) where T = Array{T}

function Base.iterate(w::Weights, state...)
    result = iterate(w.weights, state...)
    result === nothing && return nothing
    arr, state = result
    return Array(arr), state
end

# for saving weights from the user
function Weights(v::AbstractVector)
    T = foldl(promote_type, (eltype(x) for x in v))
    return Weights(FlatArray{T}.(v))
end

# for deserializing from Arrow
# here we assume the elements are all of the same type
# but the `eltype`s may not be correct (Arrow seems to
# deserialize them as `FlatArray{Any}`).
function Weights(v::AbstractVector{<:FlatArray})
    T = typeof(first(first(v).vec))
    return Weights(convert.(FlatArray{T}, v))
end

const WEIGHTS_ARROW_NAME = Symbol("JuliaLang.LegolasFlux.Weights")
ArrowTypes.arrowname(::Type{<:Weights}) = WEIGHTS_ARROW_NAME
ArrowTypes.JuliaType(::Val{WEIGHTS_ARROW_NAME}) = Weights

#####
##### `ModelRow`
#####

const ModelRow = Legolas.@row("legolas-flux.model@1",
                              weights::Weights = Weights(weights),
                              architecture_version::Union{Missing,Int})
#####
##### Utilities
#####

"""
    write_model_row(io_or_path; kwargs...)

A light wrapper around `Legolas.write` to write a table
with a single row. `kwargs` are forwarded to an internal
invocation of `Arrow.write`.
"""
function write_model_row(io_or_path, row; kwargs...)
    return Legolas.write(io_or_path, [row], LEGOLAS_SCHEMA; validate=true, kwargs...)
end

"""
    read_model_row(io_or_path) -> ModelRow

A light wrapper around `Legolas.read` to retrieve
a `ModelRow` from a table with a single row, such
as the output of [`write_model_row`](@ref)`.
"""
function read_model_row(io_or_path)
    table = Legolas.read(io_or_path; validate=true)
    rows = ModelRow.(Tables.rows(table))
    return only(rows)
end

include("functors.jl")

end # module


text:
using BinaryProvider
using BinaryBuilder
using BinaryBuilder: preferred_runner
using ObjectFile
using Base.Test
using SHA
using Compat, Compat.Random, Compat.LibGit2, Compat.Libdl

# The platform we're running on
const platform = platform_key()

# On windows, the `.exe` extension is very important
const exe_ext = Compat.Sys.iswindows() ? ".exe" : ""

# We are going to build/install libfoo a lot, so here's our function to make sure the
# library is working properly
function check_foo(fooifier_path = "fooifier$(exe_ext)",
                   libfoo_path = "libfoo.$(Libdl.dlext)")
    # We know that foo(a, b) returns 2*a^2 - b
    result = 2*2.2^2 - 1.1

    # Test that we can invoke fooifier
    @test !success(`$fooifier_path`)
    @test success(`$fooifier_path 1.5 2.0`)
    @test parse(Float64,readchomp(`$fooifier_path 2.2 1.1`)) ≈ result

    # Test that we can dlopen() libfoo and invoke it directly
    libfoo = Libdl.dlopen_e(libfoo_path)
    @test libfoo != C_NULL
    foo = Libdl.dlsym_e(libfoo, :foo)
    @test foo != C_NULL
    @test ccall(foo, Cdouble, (Cdouble, Cdouble), 2.2, 1.1) ≈ result
    Libdl.dlclose(libfoo)
end

@testset "File Collection" begin
    temp_prefix() do prefix
        # Create a file and a link, ensure that only the one file is returned by collect_files()
        f = joinpath(prefix, "foo")
        f_link = joinpath(prefix, "foo_link")
        touch(f)
        symlink(f, f_link)

        files = collect_files(prefix)
        @test length(files) == 2
        @test f in files
        @test f_link in files

        collapsed_files = collapse_symlinks(files)
        @test length(collapsed_files) == 1
        @test f in collapsed_files
    end
end

@testset "Target properties" begin
    for t in ["i686-linux-gnu", "i686-w64-mingw32", "arm-linux-gnueabihf"]
        @test BinaryBuilder.target_nbits(t) == "32"
    end

    for t in ["x86_64-linux-gnu", "x86_64-w64-mingw32", "aarch64-linux-gnu",
              "powerpc64le-linux-gnu", "x86_64-apple-darwin14"]
        @test BinaryBuilder.target_nbits(t) == "64"
    end

    for t in ["x86_64-linux-gnu", "x86_64-apple-darwin14", "i686-w64-mingw32"]
        @test BinaryBuilder.target_proc_family(t) == "intel"
    end
    for t in ["aarch64-linux-gnu", "arm-linux-gnueabihf"]
        @test BinaryBuilder.target_proc_family(t) == "arm"
    end
    @test BinaryBuilder.target_proc_family("powerpc64le-linux-gnu") == "power"

    for t in ["aarch64-linux-gnu", "x86_64-unknown-freebsd11.1"]
        @test BinaryBuilder.target_dlext(t) == "so"
    end
    @test BinaryBuilder.target_dlext("x86_64-apple-darwin14") == "dylib"
    @test BinaryBuilder.target_dlext("i686-w64-mingw32") == "dll"
end

@testset "UserNS utilities" begin
    # Test that is_ecryptfs works for something we're certain isn't encrypted
    isecfs = (false, "/proc/")
    @test BinaryBuilder.is_ecryptfs("/proc"; verbose=true) == isecfs
    @test BinaryBuilder.is_ecryptfs("/proc/"; verbose=true) == isecfs
    @test BinaryBuilder.is_ecryptfs("/proc/not_a_file"; verbose=true) == isecfs
end

@testset "Builder Dependency" begin
    temp_prefix() do prefix
        # First, let's create a Dependency that just installs a file
        begin
            ur = preferred_runner()(prefix.path; cwd="/workspace/", platform=platform, verbose=true)

            # Our simple executable file, generated by bash
            test_exe_sandbox_path = joinpath("/workspace/bin", "test_exe")
            test_exe_path = joinpath(bindir(prefix),"test_exe")
            test_exe = ExecutableProduct(test_exe_path, :test_exe)
            results = [test_exe]

            # These commands will be run within the cross-compilation environment
            script = """
            /bin/mkdir -p $(dirname(test_exe_sandbox_path))
            printf '#!/bin/bash\necho test' > $(test_exe_sandbox_path)
            /bin/chmod 775 $(test_exe_sandbox_path)
            """
            dep = Dependency("bash_test", results, script, platform, prefix)

            @test build(ur, dep; verbose=true)
            @test satisfied(dep)
            @test String(read(`$(test_exe_path)`)) == "test\n"
        end
    end

    begin
        build_path = tempname()
        mkpath(build_path)
        prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], platform)
        cd(joinpath(dirname(@__FILE__),"build_tests","libfoo")) do
            run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)

            libfoo = LibraryProduct(prefix, "libfoo", :libfoo)
            fooifier = ExecutableProduct(prefix, "fooifier", :fooifier)
            script="""
            /usr/bin/make clean
            /usr/bin/make install
            """
            dep = Dependency("fooifier", [libfoo, fooifier], script, platform, prefix)

            # Build it
            @test build(ur, dep; verbose=true)
            @test satisfied(dep; verbose=true)

            # Test the binaries
            check_foo(locate(fooifier), locate(libfoo))

            # Also test the binaries through `activate()`
            activate(prefix)
            check_foo()
            deactivate(prefix)

            # Test that `collect_files()` works:
            all_files = collect_files(prefix)
            @test locate(libfoo) in all_files
            @test locate(fooifier) in all_files
        end
        rm(build_path, recursive = true)
    end
end

libfoo_products(prefix) = [
    LibraryProduct(prefix, "libfoo", :libfoo)
    ExecutableProduct(prefix, "fooifier", :fooifier)
]
libfoo_script = """
/usr/bin/make clean
/usr/bin/make install
"""

@testset "Builder Packaging" begin
    # Clear out previous build products
    for f in readdir(".")
        if !endswith(f, ".tar.gz") || !endswith(f, ".tar.gz.256")
            continue
        end
        rm(f; force=true)
    end

    # Gotta set this guy up beforehand
    tarball_path = nothing
    tarball_hash = nothing

    begin
        build_path = tempname()
        mkpath(build_path)
        prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], platform)
        cd(joinpath(dirname(@__FILE__),"build_tests","libfoo")) do
            run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)

            # First, build libfoo
            dep = Dependency("foo", libfoo_products(prefix), libfoo_script, platform, prefix)

            @test build(ur, dep)
        end

        # Next, package it up as a .tar.gz file
        tarball_path, tarball_hash = package(prefix, "./libfoo", v"1.0.0"; verbose=true)
        @test isfile(tarball_path)

        # Delete the build path
        rm(build_path, recursive = true)
    end

    # Test that we can inspect the contents of the tarball
    contents = list_tarball_files(tarball_path)
    @test "bin/fooifier" in contents
    @test "lib/libfoo.$(Libdl.dlext)" in contents

    # Install it within a new Prefix
    temp_prefix() do prefix
        # Install the thing
        @test install(tarball_path, tarball_hash; prefix=prefix, verbose=true)

        # Ensure we can use it
        fooifier_path = joinpath(bindir(prefix), "fooifier")
        libfoo_path = joinpath(libdir(prefix), "libfoo.$(Libdl.dlext)")
        check_foo(fooifier_path, libfoo_path)
    end

    rm(tarball_path; force=true)
    rm("$(tarball_path).sha256"; force=true)
end

if lowercase(get(ENV, "BINARYBUILDER_FULL_SHARD_TEST", "false") ) == "true"
    # Perform a sanity test on each and every shard.
    @testset "Shard sanity tests" begin
        for shard_platform in supported_platforms()
            build_path = tempname()
            mkpath(build_path)
            prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], shard_platform)
            cd(joinpath(dirname(@__FILE__),"build_tests","libfoo")) do
                run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)

                # Build libfoo, warn if we fail
                dep = Dependency("foo", libfoo_products(prefix), libfoo_script, shard_platform, prefix)
                @test build(ur, dep)
            end

            # Delete the build path
            rm(build_path, recursive = true)
        end
    end
end

@testset "environment and history saving" begin
    build_path = tempname()
    mkpath(build_path)
    prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], platform)
    dep = Dependency("foo", libfoo_products(prefix), "MARKER=1\nexit 1", platform, prefix)
    @test_throws ErrorException build(ur, dep)

    # Ensure that we get a metadir, and that our history and .env files are in there!
    metadir = joinpath(prefix.path, "..", "metadir")
    @test isdir(metadir)

    hist_file = joinpath(metadir, ".bash_history")
    env_file = joinpath(metadir, ".env")
    @test isfile(hist_file)
    @test isfile(env_file)

    # Test that exit 1 is in .bash_history
    @test contains(read(open(hist_file), String), "\nexit 1\n")

    # Test that MARKER=1 is in .env:
    @test contains(read(open(env_file), String), "\nMARKER=1\n")

    # Delete the build path
    rm(build_path, recursive = true)
end

# Testset to make sure we can build_tarballs() from a local directory
@testset "build_tarballs() local directory based" begin
    build_path = tempname()
    local_dir_path = joinpath(build_path, "libfoo")
    mkpath(local_dir_path)

    cd(build_path) do
        # Just like we package up libfoo into a tarball above, we'll just copy it
        # into a new directory and use build_tarball's ability to auto-package
        # local directories to do all the heavy lifting.
        libfoo_dir = joinpath(@__DIR__, "build_tests", "libfoo")
        run(`cp -r $(libfoo_dir)/$(readdir(libfoo_dir)) $local_dir_path`)

        build_tarballs(
            [], # fake ARGS
            "libfoo",
            v"1.0.0",
            [local_dir_path],
            libfoo_script,
            [Linux(:x86_64, :glibc)],
            libfoo_products,
            [], # no dependencies
        )

        # Make sure that worked
        @test isfile("products/libfoo.v1.0.0.x86_64-linux-gnu.tar.gz")
        @test isfile("products/build_libfoo.v1.0.0.jl")
    end
end

# Testset to make sure we can build_tarballs() from a git repository
@testset "build_tarballs() Git-Based" begin
    build_path = tempname()
    git_path = joinpath(build_path, "libfoo.git")
    mkpath(git_path)

    cd(build_path) do
        # Just like we package up libfoo into a tarball above, we'll create a fake
        # git repo for it here, then build from that.
        repo = LibGit2.init(git_path)
        LibGit2.commit(repo, "Initial empty commit")
        libfoo_dir = joinpath(@__DIR__, "build_tests", "libfoo")
        run(`cp -r $(libfoo_dir)/$(readdir(libfoo_dir)) $git_path/`)
        for file in ["fooifier.c", "libfoo.c", "Makefile"]
            LibGit2.add!(repo, file)
        end
        commit = LibGit2.commit(repo, "Add libfoo files")

        # Now build that git repository for Linux x86_64
        sources = [
            git_path =>
            LibGit2.hex(LibGit2.GitHash(commit)),
        ]

        build_tarballs(
            [], # fake ARGS
            "libfoo",
            v"1.0.0",
            sources,
            "cd libfoo\n$libfoo_script",
            [Linux(:x86_64, :glibc)],
            libfoo_products,
            [], # no dependencies
        )

        # Make sure that worked
        @test isfile("products/libfoo.v1.0.0.x86_64-linux-gnu.tar.gz")
        @test isfile("products/build_libfoo.v1.0.0.jl")
    end

    rm(build_path; force=true, recursive=true)
end

@testset "build_tarballs() --only-buildjl" begin
    build_path = tempname()
    mkpath(build_path)
    cd(build_path) do
        # Clone down OggBuilder.jl
        repo = LibGit2.clone("https://github.com/staticfloat/OggBuilder", ".")

        # Check out a known-good tag
        LibGit2.checkout!(repo, hex(LibGit2.GitHash(LibGit2.GitCommit(repo, "v1.3.3-6"))))

        # Reconstruct binaries!  We don't want it to pick up BinaryBuilder.jl information from CI,
        # so wipe out those environment variables through withenv:
        blacklist = ["CI_REPO_OWNER", "CI_REPO_NAME", "TRAVIS_REPO_SLUG", "TRAVIS_TAG", "CI_COMMIT_TAG"]
        withenv((envvar => nothing for envvar in blacklist)...) do
            m = Module(:__anon__)
            eval(m, quote
                ARGS = ["--only-buildjl"]
                include(joinpath($(build_path), "build_tarballs.jl"))
            end)
        end

        # Read in `products/build.jl` to get download_info
        m = Module(:__anon__)
        download_info = eval(m, quote
            using BinaryProvider
            # Override BinaryProvider functionality so that it doesn't actually install anything
            function install(args...; kwargs...); end
            function write_deps_file(args...; kwargs...); end

            # Include build.jl file to extract download_info
            include(joinpath($build_path, "products", "build_Ogg.v1.3.3.jl"))
            download_info
        end)

        # Test that we get the info right about some of these platforms
        bin_prefix = "https://github.com/staticfloat/OggBuilder/releases/download/v1.3.3-6"
        @test download_info[Linux(:x86_64)] == (
            "$bin_prefix/Ogg.v1.3.3.x86_64-linux-gnu.tar.gz",
            "6ef771242553b96262d57b978358887a056034a3c630835c76062dca8b139ea6",
        )
        @test download_info[Windows(:i686)] == (
            "$bin_prefix/Ogg.v1.3.3.i686-w64-mingw32.tar.gz",
            "3f6f6f524137a178e9df7cb5ea5427de6694c2a44ef78f1491d22bd9c6c8a0e8",
        )
    end
end

@testset "Auditor - ISA tests" begin
    begin
        build_path = tempname()
        mkpath(build_path)
        isa_platform = Linux(:x86_64)
        prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], isa_platform)

        main_sse = ExecutableProduct(prefix, "main_sse", :main_sse)
        main_avx = ExecutableProduct(prefix, "main_avx", :main_avx)
        main_avx2 = ExecutableProduct(prefix, "main_avx2", :main_avx2)
        products = [main_sse, main_avx, main_avx2]

        cd(joinpath(dirname(@__FILE__),"build_tests","isa_tests")) do
            run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)

            # Build isa tests
            script="""
            /usr/bin/make clean
            /usr/bin/make install
            """
            dep = Dependency("isa_tests", products, script, isa_platform, prefix)

            # Build it
            @test build(ur, dep; verbose=true)

            # Ensure it's satisfied
            @test satisfied(dep; verbose=true)
        end

        # Next, test isa of these files
        readmeta(locate(main_sse)) do oh
            isa_sse = BinaryBuilder.analyze_instruction_set(oh; verbose=true)
            @test isa_sse == :core2
        end

        readmeta(locate(main_avx)) do oh
            isa_avx = BinaryBuilder.analyze_instruction_set(oh; verbose=true)
            @test isa_avx == :sandybridge
        end

        readmeta(locate(main_avx2)) do oh
            isa_avx2 = BinaryBuilder.analyze_instruction_set(oh; verbose=true)
            @test isa_avx2 == :haswell
        end

        # Delete the build path
        rm(build_path, recursive = true)
    end
end

@testset "Auditor - .dll moving" begin
    begin
        build_path = tempname()
        mkpath(build_path)
        dll_platform = Windows(:x86_64)
        prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], dll_platform)
        cd(joinpath(dirname(@__FILE__),"build_tests","libfoo")) do
            run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)

            # First, build libfoo, but with a dumb script that doesn't know to put .dll files in bin
            dumb_script = """
            /usr/bin/make clean
            /usr/bin/make install libdir=\$prefix/lib
            """
            dep = Dependency("foo", libfoo_products(prefix), libfoo_script, dll_platform, prefix)

            @test build(ur, dep; autofix=false)
        end

        # Test that libfoo puts its .dll's into lib, even on windows:
        @test !isfile(joinpath(prefix, "bin", "libfoo.dll"))
        @test isfile(joinpath(prefix, "lib", "libfoo.dll"))

        # Test that `audit()` moves it to `bin`.
        BinaryBuilder.audit(prefix; platform=dll_platform, verbose=true, autofix=true)
        @test isfile(joinpath(prefix, "bin", "libfoo.dll"))
        @test !isfile(joinpath(prefix, "lib", "libfoo.dll"))
    end
end

@testset "Auditor - absolute paths" begin
    prefix = Prefix(tempname())
    try
        sharedir = joinpath(prefix.path, "share")
        mkpath(sharedir)
        open(joinpath(sharedir, "foo.conf"), "w") do f
            write(f, "share_dir = \"$sharedir\"")
        end

        # Test that `audit()` warns about an absolute path within the prefix
        Compat.@info("Expecting a warning about share/foo.conf:")
        BinaryBuilder.audit(prefix)
    finally
        rm(prefix.path; recursive=true)
    end
end
@testset "GitHub releases build.jl reconstruction" begin
    # Download some random release that is relatively small
    product_hashes = product_hashes_from_github_release("staticfloat/OggBuilder", "v1.3.3-0")

    # Ground truth hashes for each product
    true_product_hashes = Dict(
        "aarch64-linux-gnu" => (
            "Ogg.aarch64-linux-gnu.tar.gz",
            "4150d19fe0dc773ef3917498379a9148cd780d44322fc30e44cf4a241fb8e688"
        ),
        "i686-w64-mingw32" => (
            "Ogg.i686-w64-mingw32.tar.gz",
            "3f9940c1c8614fbb40f35ab28dac9237226e3e7fcfb45d1fe7488e0289284ff8"
        ),
        "powerpc64le-linux-gnu" => (
            "Ogg.powerpc64le-linux-gnu.tar.gz",
            "cf519a13c3b343334aed8771d50b991d7ea00d0bbbd490ee0c8f5ffbd3ba65f4"
        ),
        "x86_64-linux-gnu" => (
            "Ogg.x86_64-linux-gnu.tar.gz",
            "3952d4def1505ad5090622a50662b6e0a38d1977abb7bb61d2483a47b626c807"
        ),
        "x86_64-apple-darwin14" => (
            "Ogg.x86_64-apple-darwin14.tar.gz",
            "e6c0fec453c4f833a0364fbf92a4a6e4bc738aa8de84059e12f1ee40f02dbba1"
        ),
        "x86_64-w64-mingw32" => (
            "Ogg.x86_64-w64-mingw32.tar.gz",
            "a47c33147f7e572f40178d47de54ad3a0a04e6b17a9fe3cf15791ba11203544e"
        ),
        "arm-linux-gnueabihf" => (
            "Ogg.arm-linux-gnueabihf.tar.gz",
            "dcafb1c46b4363f84fc194c28732e3080be82227e36938d883a2d0e381cae20c"
        ),
        "i686-linux-gnu" => (
            "Ogg.i686-linux-gnu.tar.gz",
            "f5227b205ee64e03f7b030af0c6b6eb9ddb6a3175efe1c650935e71dcd3ae658"
        ),
    )

    @test length(product_hashes) == length(true_product_hashes)

    for target in keys(true_product_hashes)
        @test haskey(product_hashes, target)
        product_platform = extract_platform_key(product_hashes[target][1])
        true_product_platform = extract_platform_key(true_product_hashes[target][1])
        @test product_platform == true_product_platform
        @test product_hashes[target][2] == true_product_hashes[target][2]
    end
end

include("wizard.jl")

# Run the package tests if we ask for it
if lowercase(get(ENV, "BINARYBUILDER_PACKAGE_TESTS", "false") ) == "true"
    cd("package_tests") do
        include(joinpath(pwd(), "runtests.jl"))
    end
end


text:
t = Template(; user=me)
pkg_dir = joinpath(t.dir, test_pkg)

@testset "GitHubPages" begin
    @testset "Plugin creation" begin
        p = GitHubPages()
        @test p.gitignore == ["/docs/build/", "/docs/site/"]
        @test isempty(p.assets)
        p = GitHubPages(; assets=[test_file])
        @test p.assets == [test_file]
        @test_throws ArgumentError GitHubPages(; assets=[fake_path])
    end

    @testset "Badge generation" begin
        p = GitHubPages()
        @test badges(p, me, test_pkg) ==  [
            "[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://$me.github.io/$test_pkg.jl/stable)"
            "[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://$me.github.io/$test_pkg.jl/dev)"
        ]
    end

    @testset "File generation" begin
        p = GitHubPages()
        @test gen_plugin(p, t, test_pkg) == ["docs/"]
        @test isdir(joinpath(pkg_dir, "docs"))
        @test isfile(joinpath(pkg_dir, "docs", "make.jl"))
        make = readchomp(joinpath(pkg_dir, "docs", "make.jl"))
        @test occursin("assets=String[]", make)
        @test !occursin("deploydocs", make)
        @test isdir(joinpath(pkg_dir, "docs", "src"))
        @test isfile(joinpath(pkg_dir, "docs", "src", "index.md"))
        index = readchomp(joinpath(pkg_dir, "docs", "src", "index.md"))
        @test occursin("autodocs", index)
        rm(joinpath(pkg_dir, "docs"); recursive=true)
        p = GitHubPages(; assets=[test_file])
        @test gen_plugin(p, t, test_pkg) == ["docs/"]
        make = readchomp(joinpath(pkg_dir, "docs", "make.jl"))
        # Check the formatting of the assets list.
        @test occursin(
            strip("""
            assets=[
                    "assets/$(basename(test_file))",
                ]
            """),
            make,
        )
        @test isfile(joinpath(pkg_dir, "docs", "src", "assets", basename(test_file)))
        rm(joinpath(pkg_dir, "docs"); recursive=true)
        t.plugins[TravisCI] = TravisCI()
        @test gen_plugin(p, t, test_pkg) == ["docs/"]
        make = readchomp(joinpath(pkg_dir, "docs", "make.jl"))
        @test occursin("deploydocs", make)
        rm(joinpath(pkg_dir, "docs"); recursive=true)
    end

    @testset "Package generation with GitHubPages plugin" begin
        temp_dir = mktempdir()
        t = Template(; user=me, dir=temp_dir, plugins=[GitHubPages()])
        generate(test_pkg, t; gitconfig=gitconfig)

        # Check that the gh-pages branch exists.
        repo = LibGit2.GitRepo(joinpath(t.dir, test_pkg))
        branches = map(b -> LibGit2.shortname(first(b)), LibGit2.GitBranchIter(repo))
        @test in("gh-pages", branches)
    end
end

rm(pkg_dir; recursive=true)


text:
## MvNormal ##

"""
    TuringMvNormal{Tm<:AbstractVector, TC<:Cholesky} <: ContinuousMultivariateDistribution

A multivariate Normal distribution whose covariance is dense. Compatible with Tracker.
"""
struct TuringMvNormal{Tm<:AbstractVector, TC<:Cholesky} <: ContinuousMultivariateDistribution
    m::Tm
    C::TC
end

TuringMvNormal(m::AbstractVector, A::AbstractMatrix) = TuringMvNormal(m, cholesky(A))

Distributions.dim(d::TuringMvNormal) = length(d.m)
function Distributions.rand(rng::Random.AbstractRNG, d::TuringMvNormal)
    return d.m .+ d.C.U' * randn(rng, dim(d))
end

"""
    TuringDiagNormal{Tm<:AbstractVector, Tσ<:AbstractVector} <: ContinuousMultivariateDistribution

A multivariate normal distribution whose covariance is diagonal. Compatible with Tracker.
"""
struct TuringDiagNormal{Tm<:AbstractVector, Tσ<:AbstractVector} <: ContinuousMultivariateDistribution
    m::Tm
    σ::Tσ
end

Distributions.dim(d::TuringDiagNormal) = length(d.m)
function Distributions.rand(rng::Random.AbstractRNG, d::TuringDiagNormal)
    return d.m .+ d.σ .* randn(rng, dim(d))
end
for T in (:AbstractVector, :AbstractMatrix)
    @eval Distributions.logpdf(d::TuringDiagNormal, x::$T) = _logpdf(d, x)
    @eval Distributions.logpdf(d::TuringMvNormal, x::$T) = _logpdf(d, x)
end

function _logpdf(d::TuringDiagNormal, x::AbstractVector)
    return -(dim(d) * log(2π) + 2 * sum(log.(d.σ)) + sum(abs2, (x .- d.m) ./ d.σ)) / 2
end
function _logpdf(d::TuringDiagNormal, x::AbstractMatrix)
    return -(dim(d) * log(2π) .+ 2 * sum(log.(d.σ)) .+ sum(abs2, (x .- d.m) ./ d.σ, dims=1)') ./ 2
end
function _logpdf(d::TuringMvNormal, x::AbstractVector)
    return -(dim(d) * log(2π) + logdet(d.C) + sum(abs2, zygote_ldiv(d.C.U', x .- d.m))) / 2
end
function _logpdf(d::TuringMvNormal, x::AbstractMatrix)
    return -(dim(d) * log(2π) .+ logdet(d.C) .+ sum(abs2, zygote_ldiv(d.C.U', x .- d.m), dims=1)') ./ 2
end

# zero mean, dense covariance
MvNormal(A::TrackedMatrix) = MvNormal(zeros(size(A, 1)), A)

# zero mean, diagonal covariance
MvNormal(σ::TrackedVector) = MvNormal(zeros(length(σ)), σ)

# dense mean, dense covariance
MvNormal(m::TrackedVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvNormal(m, A)
MvNormal(m::TrackedVector{<:Real}, A::Matrix{<:Real}) = TuringMvNormal(m, A)
MvNormal(m::AbstractVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvNormal(m, A)

# dense mean, diagonal covariance
function MvNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{T, <:TrackedVector{T}} where {T<:Real},
)
    return MvNormal(m, sqrt.(D.diag))
end
function MvNormal(
    m::AbstractVector{<:Real},
    D::Diagonal{T, <:TrackedVector{T}} where {T<:Real},
)
    return MvNormal(m, sqrt.(D.diag))
end

# dense mean, diagonal covariance
MvNormal(m::TrackedVector{<:Real}, σ::TrackedVector{<:Real}) = TuringDiagNormal(m, σ)
MvNormal(m::TrackedVector{<:Real}, σ::AbstractVector{<:Real}) = TuringDiagNormal(m, σ)
MvNormal(m::TrackedVector{<:Real}, σ::Vector{<:Real}) = TuringDiagNormal(m, σ)
MvNormal(m::AbstractVector{<:Real}, σ::TrackedVector{<:Real}) = TuringDiagNormal(m, σ)

# dense mean, constant variance
MvNormal(m::TrackedVector{<:Real}, σ::TrackedReal) = MvNormal(m, fill(σ, length(m)))
MvNormal(m::TrackedVector{<:Real}, σ::Real) = MvNormal(m, fill(σ, length(m)))
MvNormal(m::AbstractVector{<:Real}, σ::TrackedReal) = MvNormal(m, fill(σ, length(m)))

# dense mean, constant variance
MvNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:TrackedReal}) = MvNormal(m, A.λ)
MvNormal(m::AbstractVector{<:Real}, A::UniformScaling{<:TrackedReal}) = MvNormal(m, A.λ)
MvNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:Real}) = MvNormal(m, A.λ)

# zero mean,, constant variance
MvNormal(d::Int, σ::TrackedReal{<:Real}) = MvNormal(zeros(d), σ)

## MvLogNormal ##

struct TuringMvLogNormal{TD} <: AbstractMvLogNormal
    normal::TD
end
MvLogNormal(d::TuringMvNormal) = TuringMvLogNormal(d)
MvLogNormal(d::TuringDiagNormal) = TuringMvLogNormal(d)
Distributions.dim(d::TuringMvLogNormal) = length(d.normal)
function Distributions.rand(rng::Random.AbstractRNG, d::TuringMvLogNormal)
    return exp!(rand(rng, d.normal))
end
for T in (:AbstractVector, :AbstractMatrix)
    @eval Distributions.logpdf(d::TuringMvLogNormal, x::$T) = _logpdf(d, x)
end
for T in (:(Tracker.TrackedVector), :(Tracker.TrackedMatrix))
    @eval Distributions.logpdf(d::TuringMvLogNormal, x::$T) = _logpdf(d, x)
end
function _logpdf(d::TuringMvLogNormal, x::AbstractVecOrMat{T}) where {T<:Real}
    return insupport(d, x) ? (_logpdf(d.normal, log.(x)) - sum(log.(x))) : -Inf
end

# zero mean, dense covariance
MvLogNormal(A::TrackedMatrix) = MvLogNormal(zeros(size(A, 1)), A)

# zero mean, diagonal covariance
MvLogNormal(σ::TrackedVector) = MvLogNormal(zeros(length(σ)), σ)

# dense mean, dense covariance
MvLogNormal(m::TrackedVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))
MvLogNormal(m::TrackedVector{<:Real}, A::Matrix{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))
MvLogNormal(m::AbstractVector{<:Real}, A::TrackedMatrix{<:Real}) = TuringMvLogNormal(TuringMvNormal(m, A))

# dense mean, diagonal covariance
function MvLogNormal(
    m::TrackedVector{<:Real},
    D::Diagonal{T, <:TrackedVector{T}} where {T<:Real},
)
    return MvLogNormal(m, sqrt.(D.diag))
end
function MvLogNormal(
    m::AbstractVector{<:Real},
    D::Diagonal{T, <:TrackedVector{T}} where {T<:Real},
)
    return MvLogNormal(m, sqrt.(D.diag))
end

# dense mean, diagonal covariance
MvLogNormal(m::TrackedVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvLogNormal(TuringDiagNormal(m, σ))
MvLogNormal(m::TrackedVector{<:Real}, σ::AbstractVector{<:Real}) = TuringMvLogNormal(TuringDiagNormal(m, σ))
MvLogNormal(m::TrackedVector{<:Real}, σ::Vector{<:Real}) = TuringMvLogNormal(TuringDiagNormal(m, σ))
MvLogNormal(m::AbstractVector{<:Real}, σ::TrackedVector{<:Real}) = TuringMvLogNormal(TuringDiagNormal(m, σ))

# dense mean, constant variance
MvLogNormal(m::TrackedVector{<:Real}, σ::TrackedReal) = MvLogNormal(m, fill(σ, length(m)))
MvLogNormal(m::TrackedVector{<:Real}, σ::Real) = MvLogNormal(m, fill(σ, length(m)))
MvLogNormal(m::AbstractVector{<:Real}, σ::TrackedReal) = MvLogNormal(m, fill(σ, length(m)))

# dense mean, constant variance
MvLogNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:TrackedReal}) = MvLogNormal(m, A.λ)
MvLogNormal(m::AbstractVector{<:Real}, A::UniformScaling{<:TrackedReal}) = MvLogNormal(m, A.λ)
MvLogNormal(m::TrackedVector{<:Real}, A::UniformScaling{<:Real}) = MvLogNormal(m, A.λ)

# zero mean,, constant variance
MvLogNormal(d::Int, σ::TrackedReal{<:Real}) = MvLogNormal(zeros(d), σ)


text:
abstract type LinearSolver end

# residual
function r!(r, z, θ, κ)
    @warn "residual not defined"
    nothing
end

# residual Jacobian wrt z
function rz!(rz, z, θ)
    @warn "residual Jacobian wrt z not defined"
    nothing
end

# residual Jacobian wrt θ
function rθ!(rθ, z, θ)
    @warn "residual Jacobian wrt θ not defined"
    nothing
end

function r_update!(r, r̄)
    r .= r̄
end

# optimization spaces
abstract type Space end

# Euclidean
struct Euclidean <: Space
    n::Int
end

function candidate_point!(z̄::Vector{T}, ::Euclidean, z::Vector{T}, Δ::Vector{T}, α::T) where T
    z̄ .= z - α .* Δ
end

function update_point!(z::Vector{T}, ::Space, z̄::Vector{T}) where T
    z .= z̄
end

function mapping!(δz, s::Euclidean, δzs, z) # TODO: make allocation free
    δz .= δzs
end

using Parameters
# interior-point solver options
@with_kw mutable struct InteriorPointOptions{T}
    r_tol::T = 1.0e-5
    κ_tol::T = 1.0e-5
    κ_init::T = 1.0
    κ_scale::T = 0.1
    ls_scale::T = 0.5
    max_iter_inner::Int = 100
    max_iter_outer::Int = 10
    max_ls::Int = 50
    max_time::T = 60.0
    diff_sol::Bool = false
    res_norm::Real = Inf
    reg::Bool = false
    reg_pr_init = 0.0
    reg_du_init = 0.0
    solver::Symbol = :lu_solver
end

mutable struct ResidualMethods
    r!
    rz!
    rθ!
end

# regularize Jacobian / Hessian
function regularize!(v_pr, v_du, reg_pr, reg_du)
    v_pr .+= reg_pr
    v_du .-= reg_du
end

mutable struct InteriorPoint{T}
    s::Space
    methods::ResidualMethods
    z::Vector{T}               # current point
    z̄::Vector{T}               # candidate point
    r#::Vector{T}               # residual
    r_norm::T                  # residual norm
    r̄#::Vector{T}               # candidate residual
    r̄_norm::T                  # candidate residual norm
    rz#::SparseMatrixCSC{T,Int} # residual Jacobian wrt z
    rθ#::SparseMatrixCSC{T,Int} # residual Jacobian wrt θ
    Δ::Vector{T}               # search direction
    idx_ineq::Vector{Int}      # indices for inequality constraints
    idx_soc::Vector{Vector{Int}} # indices for second-order cone constraints
    idx_pr::Vector{Int}        # indices for primal variables
    idx_du::Vector{Int}        # indices for dual variables
    δz::Matrix{T}              # solution gradients (this is always dense)
    δzs::Matrix{T}             # solution gradients (in optimization space; δz = δzs for Euclidean)
    θ::Vector{T}               # problem data
    κ::Vector{T}               # barrier parameter
    num_var::Int
    num_data::Int
    solver::LinearSolver
    v_pr
    v_du
    reg_pr
    reg_du
    opts::InteriorPointOptions
end

function interior_point(z, θ;
        s = Euclidean(length(z)),
        num_var = length(z),
        num_data = length(θ),
        idx_ineq = collect(1:0),
        idx_soc = Vector{Int}[],
        idx_pr = collect(1:s.n),
        idx_du = collect(1:0),
        r! = r!, rz! = rz!, rθ! = rθ!,
        r  = zeros(s.n),
        rz = spzeros(s.n, s.n),
        rθ = spzeros(s.n, num_data),
        reg_pr = [0.0], reg_du = [0.0],
        v_pr = view(rz, CartesianIndex.(idx_pr, idx_pr)),
        v_du = view(rz, CartesianIndex.(idx_du, idx_du)),
        opts = InteriorPointOptions()) where T

    rz!(rz, z, θ) # compute Jacobian for pre-factorization

    InteriorPoint(
        s,
        ResidualMethods(r!, rz!, rθ!),
        z,
        zeros(length(z)),
        r,
        0.0,
        deepcopy(r),
        0.0,
        rz,
        rθ,
        zeros(s.n),
        idx_ineq,
        idx_soc,
        idx_pr,
        idx_du,
        zeros(length(z), num_data),
        zeros(s.n, num_data),
        θ,
        zeros(1),
        num_var,
        num_data,
        eval(opts.solver)(rz),
        v_pr,
        v_du,
        reg_pr, reg_du,
        opts)
end

# interior point solver
function interior_point!(ip::InteriorPoint{T}) where T

    # space
    s = ip.s

    # methods
    r! = ip.methods.r!
    rz! = ip.methods.rz!
    rθ! = ip.methods.rθ!

    # options
    opts = ip.opts
    r_tol = opts.r_tol
    κ_tol = opts.κ_tol
    κ_init = opts.κ_init
    κ_scale = opts.κ_scale
    ls_scale = opts.ls_scale
    max_iter_inner = opts.max_iter_inner
    max_iter_outer = opts.max_iter_outer
    max_time = opts.max_time
    max_ls = opts.max_ls
    diff_sol = opts.diff_sol
    res_norm = opts.res_norm
    reg = opts.reg

    # unpack pre-allocated data
    z = ip.z
    z̄ = ip.z̄
    r = ip.r
    r_norm = ip.r_norm
    r̄ = ip.r̄
    r̄_norm = ip.r̄_norm
    rz = ip.rz
    Δ = ip.Δ
    idx_ineq = ip.idx_ineq
    idx_soc = ip.idx_soc
    θ = ip.θ
    κ = ip.κ
    v_pr = ip.v_pr
    v_du = ip.v_du
    reg_pr = ip.reg_pr
    reg_du = ip.reg_du
    solver = ip.solver

    # initialize barrier parameter
    κ[1] = κ_init

    # initialize regularization
    reg_pr[1] = opts.reg_pr_init
    reg_du[1] = opts.reg_du_init

    # compute residual, residual Jacobian
    r!(r, z, θ, κ[1])
    r_norm = norm(r, res_norm)

    for i = 1:max_iter_outer
        for j = 1:max_iter_inner
            # check for converged residual
            if r_norm < r_tol
                break
            end
            # @show r_norm

            # compute residual Jacobian
            rz!(rz, z, θ)

            # regularize (fixed, TODO: adaptive)
            reg && regularize!(v_pr, v_du, reg_pr[1], reg_du[1])

            # compute step
            linear_solve!(solver, Δ, rz, r)

            # @show Δ
            # initialize step length
            α = 1.0

            # candidate point
            candidate_point!(z̄, s, z, Δ, α)

            # @show z̄
            # check cones
            iter = 0
            while cone_check(z̄, idx_ineq, idx_soc)
                α *= ls_scale
                candidate_point!(z̄, s, z, Δ, α)

                iter += 1
                if iter > max_ls
                    @error "backtracking line search fail"
                    return false
                end
            end

            # reduce norm of residual
            r!(r̄, z̄, θ, κ[1])
            r̄_norm = norm(r̄, res_norm)
            # @show r̄_norm

            while r̄_norm >= (1.0 - 0.001 * α) * r_norm
                α *= ls_scale
                candidate_point!(z̄, s, z, Δ, α)

                r!(r̄, z̄, θ, κ[1])
                r̄_norm = norm(r̄, Inf)
                # @show r̄
                # @show r̄_norm

                iter += 1
                if iter > max_ls
                    @error "line search fail"
                    return false
                end
            end

            # update
            update_point!(z, s, z̄)
            r_update!(r, r̄)
            r_norm = r̄_norm
        end
    end

    if κ[1] <= κ_tol
        # differentiate solution
        diff_sol && differentiate_solution!(ip)
        return true
    else
        # update barrier parameter
        κ[1] *= κ_scale

        # update residual
        r!(r, z, θ, κ[1])
        r_norm = norm(r, res_norm)
    end
end

function interior_point!(ip::InteriorPoint{T}, z::AbstractVector{T}, θ::AbstractVector{T}) where T
    ip.z .= z
    ip.θ .= θ
    interior_point!(ip)
end

function differentiate_solution!(ip::InteriorPoint)
    s = ip.s
    z = ip.z
    θ = ip.θ
    rz = ip.rz
    rθ = ip.rθ
    δz = ip.δz
    δzs = ip.δzs

    κ = ip.κ

    ip.methods.rz!(rz, z, θ) #TODO: maybe not needed
    ip.methods.rθ!(rθ, z, θ)

    linear_solve!(ip.solver, δzs, rz, rθ)
    @inbounds @views @. δzs .*= -1.0
    mapping!(δz, s, δzs, z)

    nothing
end

linear_solve!(solver::LinearSolver, x::Vector{T}, A::Array{T, 2}, b::Vector{T}) where T = linear_solve!(solver, x, sparse(A), b)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# Geometry for vectors/arrows/segments

struct SegmentGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement
    arrow::Bool
    filled::Bool
    tag::Symbol
end
SegmentGeometry(default_statistic=Gadfly.Stat.identity(); arrow=false, filled=false, tag=empty_tag) =
    SegmentGeometry(default_statistic, arrow, filled, tag)

"""
    Geom.segment[(; arrow=false, filled=false)]

Draw line segments from `x`, `y` to `xend`, `yend`.  Optionally specify their
`color`.  If `arrow` is `true` a `Scale` object for both axes must be
provided.  If `filled` is `true` the arrows are drawn with a filled polygon,
otherwise with a stroked line.
"""
const segment = SegmentGeometry

# Leave this as a function, pending extra arguments e.g. arrow attributes
"""
    Geom.vector[(; filled=false)]

This geometry is equivalent to [`Geom.segment(arrow=true)`](@ref).
"""
vector(; filled::Bool=false) = SegmentGeometry(arrow=true, filled=filled)

"""
    Geom.hair[(; intercept=0.0, orientation=:vertical)]

Draw lines from `x`, `y` to y=`intercept` if `orientation` is `:vertical` or
x=`intercept` if `:horizontal`.  Optionally specify their `color`.  This geometry
is equivalent to [`Geom.segment`](@ref) with [`Stat.hair`](@ref).
"""
hair(; intercept=0.0, orientation=:vertical) =
    SegmentGeometry(Gadfly.Stat.hair(intercept, orientation))

"""
    Geom.vectorfield[(; smoothness=1.0, scale=1.0, samples=20, filled=false)]

Draw a gradient vector field of the 2D function or a matrix in the `z`
aesthetic.  This geometry is equivalent to [`Geom.segment`](@ref) with
[`Stat.vectorfield`](@ref); see the latter for more information.
"""
function vectorfield(;smoothness=1.0, scale=1.0, samples=20, filled::Bool=false)
    return SegmentGeometry(
        Gadfly.Stat.vectorfield(smoothness, scale, samples),
        arrow=true, filled=filled )
end

default_statistic(geom::SegmentGeometry) = geom.default_statistic
element_aesthetics(::SegmentGeometry) = [:x, :y, :xend, :yend, :color]

function render(geom::SegmentGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)

    Gadfly.assert_aesthetics_defined("Geom.segment", aes, :x, :y, :xend, :yend)

    function arrow(x::T, y::T, xmax::T, ymax::T, xyrange::Vector{T}) where T<:Real
        dx = xmax-x
        dy = ymax-y
        vl = 0.225*hypot(dy/xyrange[2], dx/xyrange[1])
        θ =  atan(dy/xyrange[2], dx/xyrange[1])
        ϕ = pi/15
        xr = -vl*xyrange[1]*[cos(θ+ϕ), cos(θ-ϕ)]
        yr = -vl*xyrange[2]*[sin(θ+ϕ), sin(θ-ϕ)]
        [ (xmax+xr[1],ymax+yr[1]), (xmax,ymax), (xmax+xr[2],ymax+yr[2]) ]
    end

    n = length(aes.x)
    default_aes = Gadfly.Aesthetics()
    default_aes.color = fill(RGBA{Float32}(theme.default_color), n)

    aes = inherit(aes, default_aes)

    line_style = Gadfly.get_stroke_vector(theme.line_style[1])

    # Geom.vector requires information about scales

    if geom.arrow
        check = [aes.xviewmin, aes.xviewmax, aes.yviewmin, aes.yviewmax ]
        if any( map(x -> x === nothing, check) )
            error("For Geom.vector, Scale minvalue and maxvalue must be manually provided for both axes")
        end
         xyrange = [aes.xviewmax-aes.xviewmin, aes.yviewmax-aes.yviewmin]

         arrows = [ arrow(x, y, xend, yend, xyrange)
                for (x, y, xend, yend) in zip(aes.x, aes.y, aes.xend, aes.yend) ]
    end

    segments = [ [(x,y), (xend,yend)]
        for (x, y, xend, yend) in zip(aes.x, aes.y, aes.xend, aes.yend) ]

    classes = [svg_color_class_from_label( aes.color_label([c])[1] ) for c in aes.color ]

    ctx = context()

    compose!( ctx, (context(), Compose.line(segments, geom.tag),
                stroke(aes.color), linewidth(theme.line_width),
                strokedash(line_style), svgclass( classes )),
              svgclass("geometry")  )
    if geom.arrow
        if geom.filled
            compose!(ctx, (context(), Compose.polygon(arrows), fill(aes.color), strokedash([])) )
        else
            compose!(ctx, (context(), Compose.line(arrows), stroke(aes.color), linewidth(theme.line_width),
            strokedash([]))  )
        end
    end

    return ctx
end


text:
#=
chain:
- Julia version:
- Author: Dan
- Date: 2021-06-01
=#

import ..CONCEPT.getName
export
    Chain

"Represents a Chain of Residues. See also [`AbstractComposite`](@ref)"
mutable struct Chain <: AbstractChain
    id_                                         ::Union{Char, Nothing}
    number_of_children_                         ::Union{Int64,Nothing}
    parent_                                     ::Union{AbstractComposite, Nothing}
    previous_                                   ::Union{AbstractComposite, Nothing}
    next_                                       ::Union{AbstractComposite, Nothing}
    first_child_                                ::Union{AbstractComposite, Nothing}
    last_child_                                 ::Union{AbstractComposite, Nothing}
    properties_                                 ::Vector{Tuple{String,UInt8}}
    contains_selection_                         ::Union{Bool,Nothing}
    number_of_selected_children_                ::Union{Int64,Nothing}
    number_of_children_containing_selection_    ::Union{Int64,Nothing}
    selection_stamp_                            ::Union{TimeStamp,Nothing}
    modification_stamp_                         ::Union{TimeStamp,Nothing}
    trait_                                      ::Union{AbstractComposite, Nothing}
    selected_                                   ::Bool

    Chain() = new(nothing,nothing,nothing,nothing,nothing,nothing,nothing,Vector{Tuple{String,UInt8}}(),
                    nothing,nothing,nothing,nothing,nothing,nothing,false)

    Chain(   id_                                         ::Union{Char, Nothing},
             number_of_children_                         ::Union{Int64,Nothing},
             parent_                                     ::Union{AbstractComposite, Nothing},
             previous_                                   ::Union{AbstractComposite, Nothing},
             next_                                       ::Union{AbstractComposite, Nothing},
             first_child_                                ::Union{AbstractComposite, Nothing},
             last_child_                                 ::Union{AbstractComposite, Nothing},
             properties_                                 ::Vector{Tuple{String,UInt8}},
             contains_selection_                         ::Union{Bool,Nothing},
             number_of_selected_children_                ::Union{Int64,Nothing},
             number_of_children_containing_selection_    ::Union{Int64,Nothing},
             selection_stamp_                            ::Union{TimeStamp,Nothing},
             modification_stamp_                         ::Union{TimeStamp,Nothing},
             trait_                                      ::Union{AbstractComposite, Nothing},
             selected_                                   ::Bool

         ) = new(   id_                                     ,
                    number_of_children_                     ,
                    parent_                                 ,
                    previous_                               ,
                    next_                                   ,
                    first_child_                            ,
                    last_child_                             ,
                    properties_                             ,
                    contains_selection_                     ,
                    number_of_selected_children_            ,
                    number_of_children_containing_selection_,
                    selection_stamp_                        ,
                    modification_stamp_                     ,
                    trait_                                  ,
                    selected_
            )

end

getName(chain::Chain) = !isnothing(chain.id_) ? (chain.id_) : "-"


text:
module GeometryWrapper

using JavaCall

const Class = @jimport java.lang.Class
const Enum = @jimport java.lang.Enum
const FieldPosition = @jimport java.text.FieldPosition
const JString = @jimport java.lang.String
const JVector = @jimport org.hipparchus.geometry.Vector
const Locale = @jimport java.util.Locale
const LocalizedGeometryFormats = @jimport org.hipparchus.geometry.LocalizedGeometryFormats
const NumberFormat = @jimport java.text.NumberFormat
const Object = @jimport java.lang.Object
const Optional = @jimport java.util.Optional
const ParsePosition = @jimport java.text.ParsePosition
const Point = @jimport org.hipparchus.geometry.Point
const Space = @jimport org.hipparchus.geometry.Space
const StringBuffer = @jimport java.lang.StringBuffer
const VectorFormat = @jimport org.hipparchus.geometry.VectorFormat

export JVector
export LocalizedGeometryFormats
export Point
export Space
export VectorFormat
export add
export compare_to
export describe_constable
export distance
export distance1
export distance_inf
export distance_sq
export dot_product
export equals
export format
export get_available_locales
export get_class
export get_declaring_class
export get_dimension
export get_format
export get_localized_string
export get_norm
export get_norm1
export get_norm_inf
export get_norm_sq
export get_prefix
export get_separator
export get_source_string
export get_space
export get_sub_space
export get_suffix
export get_zero
export hash_code
export is_infinite
export is_na_n
export name
export negate
export normalize
export notify
export notify_all
export ordinal
export parse
export scalar_multiply
export subtract
export to_string
export value_of
export values
export wait

include("EnclosingWrapper/EnclosingWrapper.jl")
include("EuclideanWrapper/EuclideanWrapper.jl")
include("HullWrapper/HullWrapper.jl")
include("PartitioningWrapper/PartitioningWrapper.jl")
include("SphericalWrapper/SphericalWrapper.jl")
include("j_vector.jl")
include("localized_geometry_formats.jl")
include("point.jl")
include("space.jl")
include("vector_format.jl")

end


text:
"""
    function meow()

A meowing tiger?

# Example
```jldoctest
julia> cowsay("We're kind of stupid that way", cow=Cowsay.meow)
 _______________________________
< We're kind of stupid that way >
 -------------------------------
  \\
   \\ ,   _ ___.--'''`--''//-,-_--_.
      \\`"' ` || \\\\ \\ \\\\/ / // / ,-\\\\`,_
     /'`  \\ \\ || Y  | \\|/ / // / - |__ `-,
    /@"\\  ` \\ `\\ |  | ||/ // | \\/  \\  `-._`-,_.,
   /  _.-. `.-\\,___/\\ _/|_/_\\_\\/|_/ |     `-._._)
   `-'``/  /  |  // \\__/\\__  /  \\__/ \\
        `-'  /-\\/  | -|   \\__ \\   |-' |
          __/\\ / _/ \\/ __,-'   ) ,' _|'
         (((__/(((_.' ((___..-'((__,'
```
"""
function meow(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
  $thoughts
   $thoughts ,   _ ___.--'''`--''//-,-_--_.
      \\`"' ` || \\\\ \\ \\\\/ / // / ,-\\\\`,_
     /'`  \\ \\ || Y  | \\|/ / // / - |__ `-,
    /@"\\  ` \\ `\\ |  | ||/ // | \\/  \\  `-._`-,_.,
   /  _.-. `.-\\,___/\\ _/|_/_\\_\\/|_/ |     `-._._)
   `-'``/  /  |  // \\__/\\__  /  \\__/ \\
        `-'  /-\\/  | -|   \\__ \\   |-' |
          __/\\ / _/ \\/ __,-'   ) ,' _|'
         (((__/(((_.' ((___..-'((__,'
"""
return the_cow
end


text:
# ================================================================
# Damping factor
# ================================================================
#
struct LM_Damping
    _μ::Float64
    _ν::Float64

    function LM_Damping(μ,ν)
        @assert μ≥0 "μ≥0  μ=$μ"
        @assert ν>0
        new(μ,ν)
    end

    function LM_Damping(μ)
        LM_Damping(μ,2)
    end
end

get_μ(d::LM_Damping) = d._μ

increase_μ(d::LM_Damping) = LM_Damping(d._μ * d._ν, 2* d._ν)

function update_μ(d::LM_Damping,ρ::Real)
    if ρ>0
        return LM_Damping(d._μ *max(1/3,1-(2ρ-1)^3))
    end

    increase_μ(d)
end