text:
"""
    stfr! = STFRBlochSim(Tfree, Tg, TE, rftipdown, rftipup, [spoiling], [nTR], [save_transients])
    stfr! = STFRBlochSim(Tfree, Tg, TE, rftipdown, β, ϕ, [spoiling], [nTR], [save_transients])
    stfr!(spin, [workspace])

Simulate a small-tip fast recovery (STFR) scan on `spin`, overwriting the spin's
magnetization vector. The resultant magnetization is stored in `spin.M`. If
`nTR > 0` and `save_transients === true`, then `stfr!(...)` returns a `Vector`
with the magnetization vectors at the echo time for each of the `nTR` simulated
TRs.

# Arguments
- `Tfree::Real`: Time from the center of the tip-down RF pulse to the center of
  the tip-up RF pulse (ms)
- `Tg::Real`: Time from the center of the tip-up RF pulse to the center of the
  tip-down RF pulse (ms); gradient spoiling occurs during this time
- `TE::Real`: Echo time (ms)
- `rftipdown`:
  - `::Real`: Excitation (tip-down) flip angle (rad) (assumes an instantaneous
    RF pulse)
  - `::AbstractRF`: Excitation (tip-down) RF pulse
- One of the following:
  - Either
    - `β::Real`: Tip-up flip angle (rad) (assumes an instantaneous RF pulse)
    - `ϕ::Real`: Tip-up phase (rad)
  - or
    - `rftipup::AbstractRF`: Tip-up RF pulse
- `spoiling::AbstractSpoiling = IdealSpoiling()`: Type of spoiling to apply
- `nTR::Val = Val(0)`: Number of TRs to simulate; `Val(0)` indicates to simulate
  a steady-state scan
- `save_transients::Val = Val(false)`: Whether or not to return the
  magnetization vectors at the TE for each of the `nTR` simulated TRs; does
  nothing if `nTR == 0`

`workspace isa STFRBlochSimWorkspace`.
"""
struct STFRBlochSim{T<:Real,T1<:AbstractRF,T2<:AbstractRF,T3<:AbstractSpoiling,nTR,save_transients}
    Tfree::T
    Tg::T
    TE::T
    rftipdown::T1
    rftipup::T2
    spoiling::T3

    function STFRBlochSim(
        Tfree::Real,
        Tg::Real,
        TE::Real,
        rftipdown::T1,
        rftipup::T2,
        spoiling::T3,
        nTR::Val{T4},
        save_transients::Val{T5}
    ) where {T1<:AbstractRF,T2<:AbstractRF,T3<:AbstractSpoiling,T4,T5}

        Tfree >= TE + duration(rftipup) / 2 ||
            error("Tfree must be greater than or equal to TE + duration(rftipup) / 2")
        TE >= duration(rftipdown) / 2 || error("TE must not be during the tip-down pulse")
        Tg >= duration(rftipup) / 2 + spoiler_gradient_duration(spoiling) + duration(rftipdown) / 2 ||
            error("Tg must be greater than or equal to duration(rftipup) / 2 + " *
                  "duration(rftipdown) / 2 + spoiler_gradient_duration")
        (T4 isa Int && T4 >= 0) || error("nTR must be a nonnegative Int")
        T3 <: Union{<:RFSpoiling,<:RFandGradientSpoiling} && (T4 > 0 ||
            error("nTR must be positive when simulating RF spoiling"))
        T5 isa Bool || error("save_transients must be a Bool")
        T4 == 0 && T5 &&
            @warn("save_transients is true, but nTR = 0; no transients will be saved")
        args = promote(Tfree, Tg, TE)
        T = typeof(args[1])
        new{T,T1,T2,T3,T4,T5}(args..., rftipdown, rftipup, spoiling)

    end
end

STFRBlochSim(Tfree, Tg, TE, rftipdown, rftipup, spoiling::AbstractSpoiling, nTR::Val = Val(0)) = STFRBlochSim(Tfree, Tg, TE, rftipdown, rftipup, spoiling, nTR, Val(false))
STFRBlochSim(Tfree, Tg, TE, rftipdown, rftipup, nTR::Val, save_transients::Val = Val(false)) = STFRBlochSim(Tfree, Tg, TE, rftipdown, rftipup, IdealSpoiling(), nTR, save_transients)
STFRBlochSim(Tfree, Tg, TE, rftipdown, rftipup) = STFRBlochSim(Tfree, Tg, TE, rftipdown, rftipup, IdealSpoiling(), Val(0), Val(false))
STFRBlochSim(Tfree, Tg, TE, α::Real, rftipup, spoiling::AbstractSpoiling, nTR, save_transients) = STFRBlochSim(Tfree, Tg, TE, InstantaneousRF(α), rftipup, spoiling, nTR, save_transients)
STFRBlochSim(Tfree, Tg, TE, rftipdown, β::Real, ϕ::Real, args...) = STFRBlochSim(Tfree, Tg, TE, rftipdown, InstantaneousRF(-β, ϕ), args...)

Base.show(io::IO, stfr::STFRBlochSim{<:Real,<:AbstractRF,<:AbstractRF,<:AbstractSpoiling,nTR,save}) where {nTR,save} =
    print(io, "STFRBlochSim(", stfr.Tfree, ", ", stfr.Tg, ", ", stfr.TE, ", ", stfr.rftipdown, ", ", stfr.rftipup, ", ", stfr.spoiling, ", Val(", nTR, "), Val(", save, "))")

function Base.show(io::IO, ::MIME"text/plain", stfr::STFRBlochSim{<:Real,<:AbstractRF,<:AbstractRF,<:AbstractSpoiling,nTR,save}) where {nTR,save}

    print(io, "Small-Tip Fast Recovery (STFR) Bloch Simulation:")
    print(io, "\n Tfree = ", stfr.Tfree, " ms")
    print(io, "\n Tg = ", stfr.Tg, " ms")
    print(io, "\n TE = ", stfr.TE, " ms")
    print(io, "\n rftipdown (tip-down excitation pulse) = ")
    show(io, "text/plain", stfr.rftipdown)
    print(io, "\n rftipup (tip-up RF pulse) = ")
    show(io, "text/plain", stfr.rftipup)
    print(io, "\n spoiling = ")
    show(io, "text/plain", stfr.spoiling)
    if nTR == 0
        print(io, "\n steady-state")
    else
        print(io, "\n nTR = ", nTR)
        print(io, "\n save transients = ", save)
    end

end

struct STFRBlochSimWorkspace{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14}
    Atd::T1
    Btd::T2
    Ate::T3
    Bte::T4
    Atf::T3
    Btf::T4
    Atu::T5
    Btu::T6
    Atg::T3
    Btg::T4
    As::T7
    Bs::T8
    tmpA1::T9
    tmpB1::T4
    tmpA2::T9
    tmpB2::T4
    mat::T10
    vec::T11
    bm_workspace::T12
    td_workspace::T13
    tu_workspace::T14
end

function STFRBlochSimWorkspace(
    spin::AbstractSpin,
    scan::STFRBlochSim,
    bm_workspace = spin isa Spin ? nothing : BlochMcConnellWorkspace(spin)
)

    STFRBlochSimWorkspace(typeof(spin), typeof(scan), bm_workspace)

end

function STFRBlochSimWorkspace(
    spin::Union{Type{Spin{T}},Type{SpinMC{T,N}}},
    scan::Type{<:STFRBlochSim{<:Real,T1,T2,T3}},
    bm_workspace = spin <: Spin ? nothing : BlochMcConnellWorkspace(spin)
) where {T,N,T1,T2,T3}

    if T1 <: InstantaneousRF
        Atd = ExcitationMatrix{T}()
        Btd = nothing
        td_workspace = nothing
    else
        if spin <: Spin
            Atd = BlochMatrix{T}()
            Btd = Magnetization{T}()
        else
            Atd = BlochMcConnellMatrix{T}(N)
            Btd = MagnetizationMC{T}(N)
        end
        td_workspace = ExcitationWorkspace(spin, bm_workspace)
    end
    if T2 <: InstantaneousRF
        Atu = ExcitationMatrix{T}()
        Btu = nothing
        tu_workspace = nothing
    else
        if spin <: Spin
            Atu = BlochMatrix{T}()
            Btu = Magnetization{T}()
        else
            Atu = BlochMcConnellMatrix{T}(N)
            Btu = MagnetizationMC{T}(N)
        end
        tu_workspace = ExcitationWorkspace(spin, bm_workspace)
    end
    if T3 <: IdealSpoiling
        As = idealspoiling
        Bs = nothing
    elseif T3 <: RFSpoiling
        As = nothing
        Bs = nothing
    elseif spin <: Spin
        As = FreePrecessionMatrix{T}()
        Bs = Magnetization{T}()
    else
        As = BlochMcConnellMatrix{T}(N)
        Bs = MagnetizationMC{T}(N)
    end
    if spin <: Spin
        Ate = FreePrecessionMatrix{T}()
        Bte = Magnetization{T}()
        Atf = FreePrecessionMatrix{T}()
        Btf = Magnetization{T}()
        Atg = FreePrecessionMatrix{T}()
        Btg = Magnetization{T}()
        tmpA1 = BlochMatrix{T}()
        tmpB1 = Magnetization{T}()
        tmpA2 = BlochMatrix{T}()
        tmpB2 = Magnetization{T}()
        mat = Matrix{T}(undef, 3, 3)
        vec = Vector{T}(undef, 3)
    else
        Ate = BlochMcConnellMatrix{T}(N)
        Bte = MagnetizationMC{T}(N)
        Atf = BlochMcConnellMatrix{T}(N)
        Btf = MagnetizationMC{T}(N)
        Atg = BlochMcConnellMatrix{T}(N)
        Btg = MagnetizationMC{T}(N)
        tmpA1 = BlochMcConnellMatrix{T}(N)
        tmpB1 = MagnetizationMC{T}(N)
        tmpA2 = BlochMcConnellMatrix{T}(N)
        tmpB2 = MagnetizationMC{T}(N)
        mat = Matrix{T}(undef, 3N, 3N)
        vec = Vector{T}(undef, 3N)
    end
    STFRBlochSimWorkspace(Atd, Btd, Ate, Bte, Atf, Btf, Atu, Btu, Atg, Btg, As,
                          Bs, tmpA1, tmpB1, tmpA2, tmpB2, mat, vec,
                          bm_workspace, td_workspace, tu_workspace)

end

# Case when nTR = 0
function (scan::STFRBlochSim{<:Real,<:AbstractRF,<:AbstractRF,<:AbstractSpoiling,0})(
    spin::AbstractSpin,
    workspace::STFRBlochSimWorkspace = STFRBlochSimWorkspace(spin, scan)
)

    rfduration = (duration(scan.rftipdown) + duration(scan.rftipup)) / 2

    excite!(workspace.Atd, workspace.Btd, spin, scan.rftipdown, workspace.td_workspace)
    freeprecess!(workspace.Atf, workspace.Btf, spin, scan.Tfree - rfduration, workspace.bm_workspace)
    excite!(workspace.Atu, workspace.Btu, spin, scan.rftipup, workspace.tu_workspace)
    freeprecess!(workspace.Atg, workspace.Btg, spin, scan.Tg - rfduration - spoiler_gradient_duration(scan.spoiling), workspace.bm_workspace)
    spoil!(workspace.As, workspace.Bs, spin, scan.spoiling, workspace.bm_workspace)

    combine!(workspace.tmpA1, workspace.tmpB1, workspace.Atf, workspace.Btf, workspace.Atu, workspace.Btu)
    combine!(workspace.tmpA2, workspace.tmpB2, workspace.tmpA1, workspace.tmpB1, workspace.Atg, workspace.Btg)
    combine!(workspace.tmpA1, workspace.tmpB1, workspace.tmpA2, workspace.tmpB2, workspace.As, workspace.Bs)
    combine!(workspace.tmpA2, workspace.tmpB2, workspace.tmpA1, workspace.tmpB1, workspace.Atd, workspace.Btd)
    subtract!(workspace.mat, I, workspace.tmpA2)
    copyto!(workspace.vec, workspace.tmpB2)
    F = lu!(workspace.mat)
    ldiv!(F, workspace.vec)
    copyto!(spin.M, workspace.vec)

    freeprecess!(workspace.Ate, workspace.Bte, spin, scan.TE - duration(scan.rftipdown) / 2, workspace.bm_workspace)
    applydynamics!(spin, workspace.tmpB1, workspace.Ate, workspace.Bte)

end

# Case when nTR > 0
function (scan::STFRBlochSim{<:Real,<:AbstractRF,<:AbstractRF,T,nTR,save})(
    spin::AbstractSpin,
    workspace::STFRBlochSimWorkspace = STFRBlochSimWorkspace(spin, scan)
) where {T,nTR,save}

    rftd = scan.rftipdown
    rftu = scan.rftipup
    rfspoiling = T <: Union{<:RFSpoiling,<:RFandGradientSpoiling}
    rfduration = (duration(rftd) + duration(rftu)) / 2

    if rfspoiling
        rftd isa RF && (rftd.Δθ[] = rftd.Δθ_initial)
        rftu isa RF && (rftu.Δθ[] = rftu.Δθ_initial)
        Δθinc = rfspoiling_increment(scan.spoiling)
        θ = zero(Δθinc) # For knowing how much phase to remove when recording signal
        Δθ = Δθinc
    else
        excite!(workspace.Atd, workspace.Btd, spin, rftd, workspace.td_workspace)
        excite!(workspace.Atu, workspace.Btu, spin, rftu, workspace.tu_workspace)
    end

    if save
        M = Vector{typeof(spin.M)}(undef, nTR)
        freeprecess!(workspace.Atf, workspace.Btf, spin, scan.Tfree - scan.TE - duration(rftu) / 2, workspace.bm_workspace)
    else
        freeprecess!(workspace.Atf, workspace.Btf, spin, scan.Tfree - rfduration, workspace.bm_workspace)
    end
    freeprecess!(workspace.Atg, workspace.Btg, spin, scan.Tg - rfduration - spoiler_gradient_duration(scan.spoiling), workspace.bm_workspace)
    spoil!(workspace.As, workspace.Bs, spin, scan.spoiling, workspace.bm_workspace)
    combine!(workspace.tmpA1, workspace.tmpB1, workspace.Atg, workspace.Btg, workspace.As, workspace.Bs)
    freeprecess!(workspace.Ate, workspace.Bte, spin, scan.TE - duration(rftd) / 2, workspace.bm_workspace)

    for rep = 1:nTR-1

        rfspoiling && excite!(workspace.Atd, workspace.Btd, spin, rftd, workspace.td_workspace)
        applydynamics!(spin, workspace.tmpB2, workspace.Atd, workspace.Btd)
        if save
            applydynamics!(spin, workspace.tmpB2, workspace.Ate, workspace.Bte)
            M[rep] = copy(spin.M)
            if rfspoiling
                modulation = exp(im * θ)
                if spin isa Spin
                    tmp = signal(spin.M) * modulation
                    M[rep].x = real(tmp)
                    M[rep].y = imag(tmp)
                else
                    for i = 1:spin.N
                        tmp = signal(spin.M[i]) * modulation
                        M[rep][i].x = real(tmp)
                        M[rep][i].y = imag(tmp)
                    end
                end
            end
        end
        applydynamics!(spin, workspace.tmpB2, workspace.Atf, workspace.Btf)
        rfspoiling && excite!(workspace.Atu, workspace.Btu, spin, rftu, workspace.tu_workspace)
        applydynamics!(spin, workspace.tmpB2, workspace.Atu, workspace.Btu)
        applydynamics!(spin, workspace.tmpB2, workspace.tmpA1, workspace.tmpB1)

        if rfspoiling
            if rftd isa InstantaneousRF
                rftd = InstantaneousRF(rftd.α, rftd.θ + Δθ)
            else
                rftd.Δθ[] += Δθ
            end
            if rftu isa InstantaneousRF
                rftu = InstantaneousRF(rftu.α, rftu.θ + Δθ)
            else
                rftu.Δθ[] += Δθ
            end
            θ += Δθ
            Δθ += Δθinc
        end

    end

    rfspoiling && excite!(workspace.Atd, workspace.Btd, spin, rftd, workspace.td_workspace)
    applydynamics!(spin, workspace.tmpB2, workspace.Atd, workspace.Btd)
    applydynamics!(spin, workspace.tmpB2, workspace.Ate, workspace.Bte)
    if rfspoiling
        modulation = exp(im * θ)
        if spin isa Spin
            tmp = signal(spin.M) * modulation
            spin.M.x = real(tmp)
            spin.M.y = imag(tmp)
        else
            for i = 1:spin.N
                tmp = signal(spin.M[i]) * modulation
                spin.M[i].x = real(tmp)
                spin.M[i].y = imag(tmp)
            end
        end
    end
    if save
        M[nTR] = copy(spin.M)
        return M
    end

end


text:
struct_types = Type[]
# Export everything
for name in names(@__MODULE__, all=true)
    name in [Symbol("#eval"), Symbol("#include"), :include, :eval] && continue
    @eval begin
        export $name
        $name isa Type && isstructtype($name) && push!(struct_types, $name)
    end
end

function Base.getproperty(x::Ptr{<:Union{struct_types...}}, f::Symbol)
    T = eltype(x)
    fieldpos = Base.fieldindex(T, f)
    field_pointer = convert(Ptr{fieldtype(T, fieldpos)}, x + fieldoffset(T, fieldpos))
    return field_pointer
end

ext: jl


text:
module TestFits
    verbose = false

    import IntertemporalChoiceHeuristics:
        load,
        fit,
        models,
        training_losses

    inputs, weights = load()

    for model in models
        for loss in training_losses
            res = fit(model, inputs, weights, loss)
            if verbose
                println(model)
                println(loss)
                println(res)
                println()
            end
        end
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct RadiusServer <: SwaggerModel
    radiusServerAddress::Any # spec type: Union{ Nothing, String } # spec name: radiusServerAddress
    radiusServerScore::Any # spec type: Union{ Nothing, Int64 } # spec name: radiusServerScore
    radiusServerSecret::Any # spec type: Union{ Nothing, String } # spec name: radiusServerSecret

    function RadiusServer(;radiusServerAddress=nothing, radiusServerScore=nothing, radiusServerSecret=nothing)
        o = new()
        validate_property(RadiusServer, Symbol("radiusServerAddress"), radiusServerAddress)
        setfield!(o, Symbol("radiusServerAddress"), radiusServerAddress)
        validate_property(RadiusServer, Symbol("radiusServerScore"), radiusServerScore)
        setfield!(o, Symbol("radiusServerScore"), radiusServerScore)
        validate_property(RadiusServer, Symbol("radiusServerSecret"), radiusServerSecret)
        setfield!(o, Symbol("radiusServerSecret"), radiusServerSecret)
        o
    end
end # type RadiusServer

const _property_map_RadiusServer = Dict{Symbol,Symbol}(Symbol("radiusServerAddress")=>Symbol("radiusServerAddress"), Symbol("radiusServerScore")=>Symbol("radiusServerScore"), Symbol("radiusServerSecret")=>Symbol("radiusServerSecret"))
const _property_types_RadiusServer = Dict{Symbol,String}(Symbol("radiusServerAddress")=>"String", Symbol("radiusServerScore")=>"Int64", Symbol("radiusServerSecret")=>"String")
Base.propertynames(::Type{ RadiusServer }) = collect(keys(_property_map_RadiusServer))
Swagger.property_type(::Type{ RadiusServer }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_RadiusServer[name]))}
Swagger.field_name(::Type{ RadiusServer }, property_name::Symbol) =  _property_map_RadiusServer[property_name]

function check_required(o::RadiusServer)
    (getproperty(o, Symbol("radiusServerAddress")) === nothing) && (return false)
    true
end

function validate_property(::Type{ RadiusServer }, name::Symbol, val)
end


text:
module MonadTest

using Funs, Monad
using Base.Test

function main()
    # Array{T,0}
    let
        u = munit(Array{Int,0}, 42)
        @test length(u) == 1
        @test isa(u, Array{Int,0})
        @test u == fill(42)

        uu = munit(Array{Array{Int,0},0}, u)
        @test length(uu) == 1
        @test isa(uu, Array{Array{Int,0},0})
        @test uu[1] == fill(42)

        j = mjoin(uu)
        @test length(j) == 1
        @test isa(j, Array{Int,0})
        @test j == fill(42)

        b = mbind(x::Int->fill(x+1), fill(1), R=Array{Int,0})
        @test b == fill(2)

        b = mbind(Fun{Array{Int,0}}(x::Int->fill(x+1)), fill(1))
        @test b == fill(2)
    end

    # Array{T,1}
    let
        u = munit(Array{Int,1}, 42)
        @test length(u) == 1
        @test isa(u, Array{Int,1})
        @test u == [42]

        uu = munit(Array{Array{Int,1},1}, u)
        @test length(uu) == 1
        @test isa(uu, Array{Array{Int,1},1})
        @test uu[1] == [42]

        j = mjoin(uu)
        @test length(j) == 1
        @test isa(j, Array{Int,1})
        @test j == [42]

        a = Array(Array{Int,1}, 2)
        a[1] = [1,2]
        a[2] = [3,4]
        @test isa(a, Array{Array{Int,1},1})
        j2 = mjoin(a)
        @test j2 == [1,2,3,4]

        b = mbind(x::Int->Int[x,x+1], [1,2], R=Array{Int,1})
        @test b == [1,2,2,3]

        b = mbind(Fun{Array{Int,1}}(x::Int->Int[x,x+1]), [1,2])
        @test b == [1,2,2,3]

        z = mzero(Array{Int,1})
        @test isempty(z)
        @test isa(z, Array{Int,1})
        @test size(z) == (0,)

        p1 = mplus([1,2])
        p2 = mplus([1,2], [3,4])
        p3 = mplus([1,2], [3,4], [5,6])
        @test p1 == [1,2]
        @test p2 == [1,2,3,4]
        @test p3 == [1,2,3,4,5,6]
    end

    # Array{T,2}
    let
        u = munit(Array{Int,2}, 42)
        @test size(u) == (1,1)
        @test u[1] == 42

        uu = munit(Array{Array{Int,2},2}, u)
        @test size(uu) == (1,1)
        @test uu[1] == u

        j = mjoin(uu)
        @test size(j) == (1,1)
        @test j[1] == 42

        a = Array(Array{Int,2}, 2,2)
        a[1] = [11 12; 13 14]
        a[2] = [21 22; 23 24]
        a[3] = [31 32; 33 34]
        a[4] = [41 42; 43 44]
        @test isa(a, Array{Array{Int,2},2})
        j2 = mjoin(a)
        @test j2 == [[11 12; 13 14] [31 32; 33 34]; [21 22; 23 24] [41 42; 43 44]]

        b = mbind(x::Int->Int[x x+1; 2*x 2*x+1], [1 2; 3 4], R=Array{Int,2})
        @test b == [[1 2; 2 3] [2 3; 4 5]; [3 4; 6 7] [4 5; 8 9]]

        z = mzero(Array{Int,2})
        @test isempty(z)
        @test isa(z, Array{Int,2})
        @test size(z) == (0,0)

        p1 = mplus([1 2; 3 4])
        p2 = mplus([1 2; 3 4], [5 6; 7 8])
        p3 = mplus([1 2; 3 4], [5 6; 7 8], [9 10; 11 12])
        @test p1 == [1 2; 3 4]
        @test p2 == [1 2; 3 4; 5 6; 7 8]
        @test p3 == [1 2; 3 4; 5 6; 7 8; 9 10; 11 12]
    end

    # Fun
    let
        u = munit(Fun{Int}, 42)
        @test isa(u, Fun{Int})
        @test fcall(u,nothing) == 42

        uu = munit(Fun{Fun{Int}}, u)
        @test isa(uu, Fun{Fun{Int}})
        @test fcall(fcall(uu,nothing),nothing) == 42

        j = mjoin(uu)
        @test isa(j, Fun{Int})
        @test fcall(j,nothing) == 42

        inc = Fun{Int}(x->x+1)
        mul = Fun{Fun{Float64}}(x->Fun{Float64}(y->float(x*y)))
        b = mbind(mul, inc)
        @test isa(b, Fun{Float64})
        @test fcall(b,1) == 2.0
        @test fcall(b,2) == 6.0
        @test fcall(b,3) == 12.0

        b = mbind(x->Fun{Float64}(y->float(x*y)), inc, R=Fun{Float64})
        @test isa(b, Fun{Float64})
        @test fcall(b,1) == 2.0
        @test fcall(b,2) == 6.0
        @test fcall(b,3) == 12.0
    end

    # Set
    let
        u = munit(Set{Int}, 42)
        @test length(u) == 1
        @test isa(u, Set{Int})
        @test u == Set([42])

        uu = munit(Set{Set{Int}}, u)
        @test length(uu) == 1
        @test isa(uu, Set{Set{Int}})
        @test uu == Set([Set([42])])

        j = mjoin(uu)
        @test length(j) == 1
        @test isa(j, Set{Int})
        @test j == Set([42])

        s = Set([Set([1,2]), Set([3,4])])
        @test isa(s, Set{Set{Int}})
        j2 = mjoin(s)
        @test j2 == Set([1,2,3,4])

        b = mbind(x::Int->Set(Int[x,x+1]), Set([1,2]), R=Set{Int})
        @test b == Set([1,2,2,3])

        b = mbind(Fun{Set{Int}}(x::Int->Set(Int[x,x+1])), Set([1,2]))
        @test b == Set([1,2,2,3])

        z = mzero(Set{Int})
        @test isempty(z)
        @test isa(z, Set{Int})

        p1 = mplus(Set([1,2]))
        p2 = mplus(Set([1,2]), Set([3,4]))
        p3 = mplus(Set([1,2]), Set([3,4]), Set([5,6]))
        @test p1 == Set([1,2])
        @test p2 == Set([1,2,3,4])
        @test p3 == Set([1,2,3,4,5,6])
    end
end

end


text:
#======================================================================================================================
Reach-Avoid Specification
======================================================================================================================#
dfa_states = collect(1:2)
dfa_props = ["a"]
dfa_transitions = [(1, "!a", 1),
                    (1, "a", 2),
                    (2, "true", 2)]
dfa_accepting_state = 2
dfa_sink_state = nothing
dfa_initial_state = 1


text:
@testset "Statistics" begin

      # Test Statistics
      stats = Statistics()
      T = Float64
      N = 10
      dt = 0.1
      p = 3
      model = UnicycleGame(p=p)
      probsize = ProblemSize(N,model)
      x0 = rand(SVector{probsize.n})

      k = 1
      t_elap = 0.2
      res_norm = 0.3
      Δ_traj = 0.1
      dyn_vio = DynamicsViolation(N)
      con_vio = ControlViolation(N)
      sta_vio = StateViolation(N)
      opt_vio = OptimalityViolation(N)
      record!(stats, t_elap, res_norm, Δ_traj, dyn_vio, con_vio, sta_vio, opt_vio, k)
      @test stats.iter == 1

      game_con = GameConstraintValues(probsize)
      u_max =  0.1*ones(model.m)
      u_min = -0.1*ones(model.m)
      add_control_bound!(game_con, u_max, u_min)
      walls = [Wall([0.,1], [1,0], [1,1]/sqrt(2))]
      add_wall_constraint!(game_con, walls)
      pdtraj = PrimalDualTraj(probsize, dt)
      core = NewtonCore(probsize)

      opts = Options()
      Q = [Diagonal(10*ones(SVector{model.ni[i],T})) for i=1:p] # Quadratic state cost
      R = [Diagonal(0.1*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
      # Desrired state
      xf = [SVector{model.ni[1],T}([2,+0.4,0,0]),
            SVector{model.ni[2],T}([2, 0.0,0,0]),
            SVector{model.ni[3],T}([3,-0.4,0,0]),
            ]
      # Desired control
      uf = [zeros(SVector{model.mi[i],T}) for i=1:p]
      # Objectives of the game
      game_obj = GameObjective(Q,R,xf,uf,N,model)
      prob = GameProblem(N,dt,x0,model,opts,game_obj,game_con)

      k = 2
      record!(stats, prob, model, game_con, pdtraj, t_elap, Δ_traj, k)
      @test stats.iter == 2
      @test stats.outer_iter == [1,2]
      @test stats.Δ_traj == [0.1,0.1]
      @test stats.t_elap == [0.2,0.2]
      stats.t_elap

      k = 3
      i = 2
      record!(stats, prob, model, game_con, pdtraj, t_elap, Δ_traj, k, i)
      @test stats.iter == 3
      @test stats.outer_iter == [1,2,3]
      @test stats.Δ_traj == [0.1,0.1,0.1]
      @test stats.t_elap == [0.2,0.2,0.2]
      stats.t_elap

      reset!(stats)
      @test stats.iter == 0
      @test stats.outer_iter == Vector{Int}()

end


text:
module ProductionFunctionsLH

using DocStringExtensions

export output, mproducts, validate_prod_fct, n_inputs, productivities, constant_returns
export CES, ces_output, ces_mproducts
export output_cobb_douglas

abstract type AbstractProductionFunction{F1<: AbstractFloat} end

include("cobb_douglas.jl");
include("ces.jl");

end


text:
using Distributions
using LinearAlgebra
using Distances
using Statistics

##

function rbfkernel(X::AbstractArray{Float64}, a::Float64)
    exp.(- a .* pairwise(SqEuclidean(), X, dims=2))
end

function rbfkernel(
    X::AbstractArray{Float64},
    Z::AbstractArray{Float64},
    a::Float64
)
    exp.(- a .* pairwise(SqEuclidean(), X, Z, dims=2))
end

mutable struct GaussianProcessSampler
    X::AbstractMatrix{Float64}
    y::AbstractVector{Float64}
    tested::AbstractVector{Bool}
    a::Float64
    σ::Float64
    K::Matrix{Float64}
    b::Float64
    offset::Float64

    function GaussianProcessSampler(
            X::Array{Float64};
            dim::Int=ndims(X),
            σ::Float64=0.001,
            a::Float64=0.5,
            b::Float64=1.0,
            offset::Float64=0.0)
        N = size(X, dim)
        y = zeros(N)
        tested = zeros(Bool, N)
        K = b * rbfkernel(X, a)
        new(X, y, tested, a, σ, K, b, offset)
    end
end

function addobs!(
        gp::GaussianProcessSampler,
        pos::Vector{Int},
        ynew::Vector{Float64})
    for (i, yi) in zip(pos, ynew)
        gp.tested[i] = true
        gp.y[i] = yi
    end
end

function gpsample(gp::GaussianProcessSampler, n::Int=1)
    # gp regression mean and variance
    @assert n <= length(gp.y) - sum(gp.tested)
    if sum(gp.tested) == 0  # choose at random
        M = length(gp.y)
        idx = sample(1:M, n, replace=false)
        vals = [-Inf for _ in 1:n]
        Xidx = gp.X[:, idx]
        return idx, vals, Xidx
    end
    y = gp.y
    σ2 = gp.σ^2
    s = gp.tested
    K = gp.K
    offset = gp.offset
    #
    A = K[s, s] + (σ2 + 1e-12)*I
    β = y[s] .- offset
    f = K[:, s] * (A \ β) .+ offset
    Σ = K - K[:, s] * (A \ K[s, :])
    # sample from multivariate normal
    # and take elements with higher sampled values
    Σ = Symmetric(Σ + 1e-12I)
    vals = Float64[]
    idx = Int[]
    idx_ = Set{Int}()
    for _ in 1:n
        W = rand(MvNormal(f, Σ))
        Z = sortperm(W, rev=true)
        i = 1
        while Z[i] ∈ idx_
            i += 1
        end
        push!(idx, Z[i])
        push!(idx_, Z[i])
        push!(vals, W[Z[i]])
    end
    Xidx = gp.X[:, idx]
    predidx = f[idx]
    bandidx = sqrt.(diag(Σ)[idx])
    return idx, vals, Xidx, predidx, bandidx
end

function gpeval(gp::GaussianProcessSampler)
    # gp regression mean and variance
    @assert sum(gp.tested) > 0
    σ2 = gp.σ^2
    s = gp.tested
    K = gp.K
    y = gp.y
    N = length(gp.y)
    offset = gp.offset
    #
    A = K[s, s] + (σ2 + 1e-12)*I
    β = y[s] .- offset
    f = K[:, s] * (A \ β) .+ offset
    Σ = K - K[:, s] * (A \ K[s, :])
    band = sqrt.(diag(Σ))
    return f, band
end

##

mutable struct RandomGaussianProcessSampler
    dists::Vector{T} where T <: Distribution
    X::Vector{Vector{Float64}}
    y::Vector{Float64}
    a::Float64
    σ::Float64
    b::Float64
    offset::Float64
    adaptive_normalization::Bool

    function RandomGaussianProcessSampler(
        dists::Vector{T};
        σ::Float64=0.001,
        a::Float64=0.5,
        b::Float64=1.0,
        offset::Float64=0.0,
        adaptive_normalization::Bool=false
    ) where T <: Distribution
        X = Vector{Float64}[]
        y = Float64[]
        new(dists, X, y, a, σ, b, offset, adaptive_normalization)
    end
end
Base.length(gp::RandomGaussianProcessSampler) = length(gp.X)

function addobs!(
    gp::RandomGaussianProcessSampler,
    xnew::Vector{Float64},
    ynew::Float64
)
    # update X, y
    push!(gp.X, xnew)
    push!(gp.y, ynew)
    # todo: eval kernel smartly on the fly
end

function gpsample(
    gp::RandomGaussianProcessSampler;
    batch_size::Int = 100
)
    # gp regression mean and variance
    n = batch_size
    candidates = [rand(d, n) for d in gp.dists]
    if length(gp) == 0  # choose at random
        pars = vcat([x' for x in candidates]...)
        pred = fill(-Inf, n)
        fhat = fill(-Inf, n)
        band = zeros(n)
        return pars, pred, fhat, band
    else
        y = gp.y
        if gp.adaptive_normalization
            μ_seen = mean(y)
            σ_seen = std(y)
            y = (y .- μ_seen) ./ (σ_seen + eps())
        end
        Xoldmat = hcat([x for x in gp.X]...)
        Xnewmat = vcat([x' for x in candidates]...)
        K00 = gp.b * rbfkernel(Xoldmat, gp.a)
        K11 = gp.b * rbfkernel(Xnewmat, gp.a)
        K01 = gp.b * rbfkernel(Xoldmat, Xnewmat, gp.a)
        K10 = K01'
        s = length(gp)
        #
        σ2 = gp.σ^2
        A = K00 + (σ2 + 1e-12)I
        fhat = K10 * (A \ y)
        Σ = K11 - K10 * (A \ K01)
        Σ = Symmetric(Σ + 1e-12I)
        pred = rand(MvNormal(fhat, Σ))
        band = sqrt.(diag(Σ))
        order = sortperm(pred, rev=true)
        pred = pred[order]
        band = band[order]
        fhat = fhat[order]
        Xnewmat = Xnewmat[:, order]
        if gp.adaptive_normalization
            pred .= (σ_seen + eps()) .* pred .+ μ_seen
            fhat .= (σ_seen + eps()) .* fhat .+ μ_seen
            band .= (σ_seen + eps()) .* band
        end
        return Xnewmat, pred, fhat, band
    end
end

function gpsample(
    gp::RandomGaussianProcessSampler,
    n::Int;
    kwargs...
)
    xs = AbstractVector{Float64}[]
    fs = Float64[]
    ps = Float64[]
    bs = Float64[]
    for i in 1:n
        x, p, f, b = gpsample(gp; kwargs...)
        push!(xs, x[:, 1])
        push!(fs, f[1])
        push!(ps, p[1])
        push!(bs, b[1])
    end
    xs = hcat(xs...)
    xs, ps, fs, bs
end

function gpeval(gp::RandomGaussianProcessSampler)
    # gp regression mean and variance
    @assert length(gp) > 0  "add data first"
    y = gp.y
    if gp.adaptive_normalization
        μ_seen = mean(y)
        std_seen = std(y)
        y = (y .- μ_seen) ./ (std_seen + eps())
    end
    σ2 = gp.σ^2
    X = hcat([x for x in gp.X]...)
    (size(X, 2) ≥ 1_000) && warn("kernel matrix too large, try sampling less candidates")
    K = gp.b * rbfkernel(X, gp.a)
    s = length(gp)
    #
    A = K + (σ2 + 1e-12)I
    f = K * (A \ y)
    Σ = K - K * (A \ K)
    Σ = Symmetric(Σ + 1e-12I)
    band = sqrt.(diag(Σ))
    if gp.adaptive_normalization
        f = (σ_seen + eps()) .* f .+ μ_seen
        band = (σ_seen + eps()) .* band
    end
    return f, band
end

function gpeval(
    gp::RandomGaussianProcessSampler,
    Xnew::Vector{Vector{Float64}}
)
    # gp regression mean and variance
    if length(gp) == 0  # choose at random
        pars = hcat([x for x in Xnew]...)
        n = length(Xnew)
        pred = fill(-Inf, n)
        band = zeros(n)
        return pred, band
    else
        y = gp.y
        if gp.adaptive_normalization
            μ_seen = mean(y)
            σ_seen = std(y)
            y = (y .- μ_seen) ./ (σ_seen + eps())
        end
        σ2 = gp.σ^2
        Xoldmat = hcat([x for x in gp.X]...)
        Xnewmat = hcat([x for x in Xnew]...)
        K00 = gp.b * rbfkernel(Xoldmat, gp.a)
        K11 = gp.b * rbfkernel(Xnewmat, gp.a)
        K01 = gp.b * rbfkernel(Xoldmat, Xnewmat, gp.a)
        K10 = K01'
        offset = gp.offset
        s = length(gp)
        #
        A = K00 + (σ2 + 1e-12)I
        β = y .- offset
        pred = K10 * (A \ β) .+ offset
        Σ = K11 - K10 * (A \ K01)
        Σ = Symmetric(Σ + 1e-12I)
        band = sqrt.(diag(Σ))
        if gp.adaptive_normalization
            pred = (σ_seen + eps()) .* pred .+ μ_seen
            band = (σ_seen + eps()) .* band
        end
        return pred, band
    end
end

## test random gp
# using Plots
# N = 11
# a = 0.2
# σ = 0.05
# b = 1.0
# scale = 100.0
# offset = -10.0 / scale
# x = collect(range(0., stop=2π, length=N))
# y = zeros(N)
# draw(x::Float64, σ::Float64) = (sin(2x) + σ * randn()) / scale + offset
# for i in 1:N
#     y[i] = draw(x[i], σ)
# end
# dists = [Uniform(0.0, 2π)]
# gp = RandomGaussianProcessSampler(dists, a=a, σ=σ, b=b, adaptive_normalization=true)

# for i in 1:N
#     addobs!(gp, [x[i]], y[i])
# end

# @time pars, pred, fhat, band = gpsample(gp, 16; batch_size=500)
# x_sample = vec(pars)
# y_sample = pred

# xseq = [[xi] for xi in collect(range(0., stop=2π, length=50))]
# pred_eval, band_eval = gpeval(gp, xseq)
# xseq = vcat(xseq...)

# ytruth = sin.(2xseq) ./ scale .+ offset
# p1 = plot(xseq, ytruth, label="truth", linestyle=:dash, color="black")
# plot!(p1, vec(xseq), pred_eval, ribbon=1.96 * band_eval, fillalpha=.2, label="fit_0", color="blue")
# plot!(p1, x, gp.y, st = :scatter, label="data", color="blue", alpha=0.4)
# plot!(p1, x_sample, y_sample, st = :scatter, label="candidate_1", color=:red)

# ylims!(p1, -2.5, 2.5)

##

# #  make sure it works fine!
# using Plots
# function test()
#     N = 200
#     a = 0.5
#     σ = 0.1
#     b = 1.0
#     x = collect(range(0., stop=2π, length=N))
#     y = zeros(N)
#     seen = zeros(Bool, N)
#     seen[[1, div(N, 2), N]] .= true
#     draw(x::Float64, σ::Float64) = sin(x) + σ * randn()
#     for i in findall(seen)
#         y[i] = draw(x[i], σ)
#     end
#     ytruth = sin.(x)
#     gp = GaussianProcessSampler(x, a=a, σ=σ, b=b)
#     addobs!(gp, findall(seen), y[seen])
#     f, band = gpeval(gp)
#     p1 = plot(x, ytruth, label="truth", linestyle=:dash, color="black")
#     plot!(p1, x, f, ribbon=1.96 * band, fillalpha=.2, label="fit_0", color="blue")
#     plot!(p1, x[seen], gp.y[seen], st = :scatter, label="candidate_0", color="blue", alpha=0.4)
#     ylims!(p1, -2.5, 2.5)

#     iu, yu, _, _, _ = gpsample(gp, 5)
#     p2 = plot(x[iu], yu, st = :scatter, label="candidate_1", color="red", alpha=0.4)
#     obs = draw.(x[iu], σ)
#     plot!(p2, x[iu], obs, st = :scatter, label="obs_1", color="blue", alpha=0.4)
#     addobs!(gp, iu, obs)
#     f, band = gpeval(gp)
#     plot!(p2, x, ytruth, label="truth", linestyle=:dash, color="black")
#     plot!(p2, x, f, ribbon=1.96 * band, fillalpha=.2, label="fit_1", color="red")
#     ylims!(p2, -2.5, 2.5)

#     iu, yu, _, _, _ = gpsample(gp, 5)
#     p3 = plot(x[iu], yu, st = :scatter, label="candidate_2", color="red", alpha=0.4)
#     obs = draw.(x[iu], σ)
#     plot!(p3, x[iu], obs, st = :scatter, label="obs_2", color="blue", alpha=0.4)
#     addobs!(gp, iu, obs)
#     f, band = gpeval(gp)
#     plot!(p3, x, ytruth, label="truth", linestyle=:dash, color="black")
#     plot!(p3, x, f, ribbon=1.96 * band, fillalpha=.2, label="fit_2", color="red")
#     ylims!(p3, -2.5, 2.5)

#     iu, yu, _, _, _ = gpsample(gp, 5)
#     p4 = plot(x[iu], yu, st = :scatter, label="candidate_3", color="red", alpha=0.4)
#     obs = draw.(x[iu], σ)
#     plot!(p4, x[iu], obs, st = :scatter, label="obs_3", color="blue", alpha=0.4)
#     addobs!(gp, iu, obs)
#     f, band = gpeval(gp)
#     plot!(p4, x, ytruth, label="truth", linestyle=:dash, color="black")
#     plot!(p4, x, f, ribbon=1.96 * band, fillalpha=.2, label="fit_3", color="red")
#     ylims!(p4, -2.5, 2.5)

#     return plot(p1, p2, p3, p4, layout=(2, 2), size=(1000, 600))
# end
# test()
# savefig(test(), "bayesian_optimization_example.png")


text:
module segment
import Term
import Term: remove_markup, remove_ansi, unescape_brackets
import ..style: apply_style, MarkupStyle
import ..measure: Measure

export Segment

# ---------------------------------------------------------------------------- #
#                                    SEGMENT                                   #
# ---------------------------------------------------------------------------- #

"""
    Segment

stores one piece of text with all the styles applied to it.
"""
struct Segment
    text::String   # text with ANSI codes injected
    measure::Measure       # measure of plain text
end

# ------------------------------- constructors ------------------------------- #
"""
    Segment(text::AbstractString)

Construct a Segment out of a string with markup.
"""
function Segment(text)
    return Segment(apply_style(text), Measure(text))
end

"""
    Segment(text::Union{Segment, AbstractString}, markup::AbstractString)

Construct a Segment out of a plain string and a markup string with style info
"""
function Segment(text, markup::String)
    Segment("[$markup]" * text * "[/$markup]")
end

Segment(seg::Segment) = seg

# --------------------------------- printing --------------------------------- #
"""print styled in stdout, info otherwise"""
function Base.show(io::IO, seg::Segment)
    if io isa IOContext{Base.TTY}
        print(io, "Segment \e[2m(size: $(seg.measure))\e[0m")
    else
        print(io, unescape_brackets(seg.text))

    end
end

# ---------------------------------------------------------------------------- #
#                                    LAYOUT                                    #
# ---------------------------------------------------------------------------- #

"""
    Term.fillin(segments::Vector{Segment})::Vector{Segment}

Ensure that for each segment the text has the same width
"""
function Term.fillin(segments::Vector{Segment})::Vector{Segment}
    widths = [seg.measure.w for seg in segments]
    w = max(widths...)

    filled::Vector{Segment} = []
    for seg in segments
        push!(filled, Segment(seg.text * " "^(w-seg.measure.w)))
    end
    return filled
end

# -------------------------------- concatenate ------------------------------- #
"""
concatenate strings and segments
"""
Base.:*(seg::Segment, str::AbstractString) = Segment(seg.text * str)
Base.:*(str::AbstractString, seg::Segment) = Segment(str * seg.text)
Base.:*(seg1::Segment, seg2::Segment) = Segment(seg1.text * seg2.text)

end


text:
function StopWatches.time_table()
    table = Tables.columntable(append_times!(StructVector{Time}(undef, 0)))
    total_time = uconvert.(u"s", float.(table.total_time_ns) .* u"ns")
    average_time = total_time ./ table.ncalls
    return StructVector(
        _time_table(;
            average_time = average_time,
            total_time = total_time,
            # Append the rest:
            table...,
        ),
    )
end

_time_table(;
    # Removing from the table:
    total_time_ns,
    # Keeping the rest:
    table...,
) = values(table)

StopWatches.time_dataframe() = DataFrame(StopWatches.time_table(), copycols = false)

function StopWatches.summary()
    display(StopWatches.time_dataframe())
    return
end


text:
function step(model, scheme::ModifiedBridge, t0, x0, Δw)
  x1 = scheme.s
  t1 = scheme.t
  μ = (x1 - x0) / (t1 - t0)
  σ = sqrt((t1 - (t0 + scheme.Δt)) / (t1 - t0)) * diffusion(model, t0, x0)
  x = x0 + μ * scheme.Δt + σ * Δw
  x
end

function _normal_transition_params(model, scheme::ModifiedBridge, t0, x0, s1)
  x1 = scheme.s
  t1 = scheme.t
  μ = x0 + (x1 - x0) / (t1 - t0) * scheme.Δt
  σ = diffusion(model, t0, x0)
  Σ = (t1 - (t0 + scheme.Δt)) / (t1 - t0) * scheme.Δt * σ * σ'
  μ, Σ
end


text:
using Documenter, Colors

abstract type SVG end
function Base.show(io::IO, ::MIME"text/html", svg::SVG)
    write(io, "")
    write(io, take!(svg.buf))
    write(io, "")
    flush(io)
end

include("crosssectionalcharts.jl")
include("colordiffcharts.jl")
include("colormaps.jl")
include("namedcolorcharts.jl")
include("sampleimages.jl")

makedocs(
    clean = false,
    modules = [Colors],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
                             assets = ["assets/resize_svg.js", "assets/favicon.ico"]),
    checkdocs = :exports,
    sitename = "Colors",
    pages    = Any[
        "Introduction"                => "index.md",
        "Construction and Conversion" => "constructionandconversion.md",
        "Color Differences"           => "colordifferences.md",
        "Colormaps and Colorscales"   => "colormapsandcolorscales.md",
        "Named Colors"                => "namedcolors.md",
        "Advanced Functions"          => "advancedfunctions.md",
        "References"                  => "references.md",
        "Index"                       => "functionindex.md",
        ]
    )

deploydocs(
    repo = "github.com/JuliaGraphics/Colors.jl.git",
    target = "build")


text:
export  ApproximateValueIteration,
        solve,
        NearestNeighborValueFunction,
        fit!,
        LocallyWeightedValueFunction,
        MultilinearValueFunction,
        SimplexValueFunction,
        LinearRegressionValueFunction

"""
    struct ApproximateValueIteration
        Uθ # initial parameterized value function that supports fit!
        S # set of discrete states for performing backups
        k_max # maximum number of iterations
    end
"""
struct ApproximateValueIteration
    Uθ # initial parameterized value function that supports fit!
    S # set of discrete states for performing backups
    k_max # maximum number of iterations
end

"""
    function solve(M::ApproximateValueIteration, 𝒫::MDP)
"""
function solve(M::ApproximateValueIteration, 𝒫::MDP)
    Uθ, S, k_max = M.Uθ, M.S, M.k_max
    for k in 1:k_max
        U = [backup(𝒫, Uθ, s) for s in S]
        fit!(Uθ, S, U)
    end

    return ValueFunctionPolicy(𝒫, Uθ)
end

"""
    mutable struct NearestNeighborValueFunction
        k # number of neighbors
        d # distance function d(s, s′)
        S # set of discrete states
        θ # vector of values at states in S
    end
"""
mutable struct NearestNeighborValueFunction
    k # number of neighbors
    d # distance function d(s, s′)
    S # set of discrete states
    θ # vector of values at states in S
end

"""
    function (Uθ::NearestNeighborValueFunction)(s)

Approximate value function at s as average of value function over k nearest states (whose value function value is already known).
"""
function (Uθ::NearestNeighborValueFunction)(s)
    # distances from s to each state s′ in state subset S
    dists = [Uθ.d(s,s′) for s′ in Uθ.S]
    # indicies of shortest 1:k distances (ind[1] == index of closest neighbor)
    ind = sortperm(dists)[1:Uθ.k]
    # average value of k nearest neighbors
    return mean(Uθ.θ[i] for i in ind)
end

"""
    function fit!(Uθ::NearestNeighborValueFunction, S, U)
"""
function fit!(Uθ::NearestNeighborValueFunction, S, U)
    Uθ.θ = U
    return Uθ
end

"""
    mutable struct LocallyWeightedValueFunction
        k # kernel function k(s, s′)
        S # set of discrete states
        θ # vector of values at states in S
    end
"""
mutable struct LocallyWeightedValueFunction
    k # kernel function k(s, s′)
    S # set of discrete states
    θ # vector of values at states in S
end

"""
    function (Uθ::LocallyWeightedValueFunction)(s)
"""
function (Uθ::LocallyWeightedValueFunction)(s)
    w = normalize([Uθ.k(s,s′) for s′ in Uθ.S], 1)
    return Uθ.θ ⋅ w
end

"""
    function fit!(Uθ::LocallyWeightedValueFunction, S, U)
"""
function fit!(Uθ::LocallyWeightedValueFunction, S, U)
    Uθ.θ = U
    return Uθ
end

# ------------------ MultilinearValueFunction ---------------------

# Also implemented in Interpolations.jl

"""
    mutable struct MultilinearValueFunction
        o # position of lower-left corner
        δ # vector of widths
        θ # vector of values at states in S
    end
"""
mutable struct MultilinearValueFunction
    o # position of lower-left corner
    δ # vector of widths
    θ # vector of values at states in S
end

"""
    function (Uθ::MultilinearValueFunction)(s)

Use multilinear interpolation to estimate the value of state vector `s` for known state values `θ` over a grid defined by lower-left vertex `o` and vector of widths `δ`.  Verticies of the grid can all be written as `o + δ.*i` for some non-negative integer vector `i`.
"""
function (Uθ::MultilinearValueFunction)(s)
    o, δ, θ = Uθ.o, Uθ.δ, Uθ.θ
    Δ = (s - o)./δ
    # Multidimensional index of lower-left cell
    i = min.(floor.(Int, Δ) .+ 1, size(θ) .- 1)
    vertex_index = similar(i)
    d = length(s)
    u = 0.0
    for vertex in 0:2^d-1
        weight = 1.0
        for j in 1:d
            # Check whether jth bit is set
            if vertex & (1 << (j-1)) > 0
                vertex_index[j] = i[j] + 1
                weight *= Δ[j] - i[j] + 1
            else
                weight *= i[j] - Δ[j]
            end
        end
        u += θ[vertex_index...]*weight
    end

    return u
end

"""
    function fit!(Uθ::MultilinearValueFunction, S, U)
"""
function fit!(Uθ::MultilinearValueFunction, S, U)
    Uθ.θ = U
    return Uθ
end

# --------------- SimplexValueFunction ------------------
#
# Also implemented in GridInterpolations.jl

"""
    mutable struct SimplexValueFunction
        o # position of lower-left corner
        δ # vector of widths
        θ # vector of values at states in S
    end
"""
mutable struct SimplexValueFunction
    o # position of lower-left corner
    δ # vector of widths
    θ # vector of values at states in S
end

"""
    function (Uθ::SimplexValueFunction)(s)

Estimate value of state vector `s` for known values `θ` over a grid defined by a lower-left vertex `o` and vector of widths `δ`.  Vertices of the grid can all be written `o + δ.*i` for some non-negative integral vector `i`.
"""
function (Uθ::SimplexValueFunction)(s)
    Δ = (s - Uθ.o)./Uθ.δ
    # Multidimensional index of upper-right cell
    i = min.(floor.(Int, Δ) .+ 1, size(Uθ.θ) .- 1) .+ 1
    u = 0.0
    s′ = (s - (Uθ.o + Uθ.δ.*(i.-2))) ./ Uθ.δ
    p = sortperm(s′) # increasing order
    w_tot = 0.0
    for j in p
        w = s′[j] - w_tot
        u += w*Uθ.θ[i...]
        i[j] -= 1
        w_tot += w
    end

    u += (1 - w_tot)*Uθ.θ[i...]
    return u
end

"""
    function fit!(Uθ::SimplexValueFunction, S, U)
"""
function fit!(Uθ::SimplexValueFunction, S, U)
    Uθ.θ = U
    return Uθ
end

# -------------- Linear Regression Value Function ---------------

"""
    mutable struct LinearRegressionValueFunction
        β # basis vector function
        θ # vector of paramters
    end
"""
mutable struct LinearRegressionValueFunction
    β # basis vector function
    θ # vector of paramters
end

"""
    function (Uθ::LinearRegressionValueFunction)(s)
"""
function (Uθ::LinearRegressionValueFunction)(s)
    return Uθ.β(s) ⋅ Uθ.θ
end

"""
    function fit!(Uθ::LinearRegressionValueFunction, S, U)
"""
function fit!(Uθ::LinearRegressionValueFunction, S, U)
    X = hcat([Uθ.β(s) for s in S]...)'
    Uθ.θ = pinv(X)*U
    return Uθ
end


text:
function diffeq(du,u,p,t)
    CycE = p[C.CycElevel] - u[V.CycEp27]

    Vdp27 = p[C.kd27] + (p[C.kd27e]*CycE)*u[V.Skp2]
    Vdcyce = p[C.kdcyce] + p[C.kdcycee]*CycE
    Vdskp2 = p[C.kdskp2] + p[C.kdskp2c1]*u[V.Cdh1]

    Vicdh1 = p[C.kicdh1e]*CycE

    du[V.p27T] = p[C.ks27] - Vdp27*u[V.p27T]

    du[V.Skp2] = p[C.ksskp2] - Vdskp2*u[V.Skp2]

    du[V.CycEp27] = p[C.kasse]*(p[C.CycElevel]-u[V.CycEp27])*(u[V.p27T]-u[V.CycEp27])-(p[C.kdise]+Vdp27+Vdcyce)*u[V.CycEp27]

    du[V.EmiC] = p[C.kasec]*(p[C.Cdh1T]-u[V.EmiC])*(p[C.Emi1T]-u[V.EmiC]) - (p[C.kdiec]+p[C.kdemi1])*u[V.EmiC]

    du[V.Cdh1dp] = p[C.kacdh1]*(p[C.Cdh1T]-u[V.Cdh1dp]) - Vicdh1*u[V.Cdh1dp]

    du[V.Cdh1] = (p[C.kdiec]+p[C.kdemi1])*(u[V.Cdh1dp]-u[V.Cdh1]) - p[C.kasec]*u[V.Cdh1]*(p[C.Emi1T]-u[V.EmiC])+p[C.kacdh1]*(p[C.Cdh1T]-u[V.EmiC]-u[V.Cdh1])-Vicdh1*u[V.Cdh1]

end

function get_steady_state(p::Vector{Float64})
    tspan::Tuple{Float64,Float64} = (0.0,Inf)
    u0::Vector{Float64} = zeros(V.NUM)

    prob = ODEProblem(diffeq,u0,tspan,p)
    prob = SteadyStateProblem(prob)
    sol = solve(prob,DynamicSS(CVODE_BDF()),dt=1.0)

    return sol.u
end

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


text:
#HDF5 Plots: Save/replay plots to/from HDF5
#-------------------------------------------------------------------------------

#==Usage
===============================================================================
Write to .hdf5 file using:
    p = plot(...)
   Plots.hdf5plot_write(p, "plotsave.hdf5")

Read from .hdf5 file using:
    pyplot() #Must first select backend
    pread = Plots.hdf5plot_read("plotsave.hdf5")
    display(pread)
==#

#==TODO
===============================================================================
 1. Support more features.
    - GridLayout known not to be working.
 2. Improve error handling.
    - Will likely crash if file format is off.
 3. Save data in a folder parallel to "plot".
    - Will make it easier for users to locate data.
    - Use HDF5 reference to link data?
 4. Develop an actual versioned file format.
    - Should have some form of backward compatibility.
    - Should be reliable for archival purposes.
 5. Fix construction of plot object with hdf5plot_read.
    - Not building object correctly when backends do not natively support
      a certain feature (ex: :steppre)
==#

import FixedPointNumbers: N0f8 #In core Julia

#Dispatch types:
struct HDF5PlotNative; end #Indentifies a data element that can natively be handled by HDF5
struct HDF5CTuple; end #Identifies a "complex" tuple structure

mutable struct HDF5Plot_PlotRef
	ref::Union{Plot, Nothing}
end

#==Useful constants
===============================================================================#
const _hdf5_nullable{T} = Union{T, Nothing}
const _hdf5_plotroot = "plot"
const _hdf5_dataroot = "data" #TODO: Eventually move data to different root (easier to locate)?
const _hdf5plot_datatypeid = "TYPE" #Attribute identifying type
const _hdf5plot_countid = "COUNT" #Attribute for storing count

#Dict has problems using "Types" as keys.  Initialize in "_initialize_backend":
const HDF5PLOT_MAP_STR2TELEM = Dict{String, Type}()
const HDF5PLOT_MAP_TELEM2STR = Dict{Type, String}()

#Don't really like this global variable... Very hacky
const HDF5PLOT_PLOTREF = HDF5Plot_PlotRef(nothing)

#Simple sub-structures that can just be written out using _hdf5plot_gwritefields:
const HDF5PLOT_SIMPLESUBSTRUCT = Union{Font, BoundingBox,
	GridLayout, RootLayout, ColorGradient, SeriesAnnotations, PlotText,
	Shape,
}

#==
===============================================================================#
is_marker_supported(::HDF5Backend, shape::Shape) = true

if length(HDF5PLOT_MAP_TELEM2STR) < 1
    #Possible element types of high-level data types:
    telem2str = Dict{String, Type}(
        "NATIVE" => HDF5PlotNative,
        "VOID" => Nothing,
        "BOOL" => Bool,
        "SYMBOL" => Symbol,
        "TUPLE" => Tuple,
        "CTUPLE" => HDF5CTuple, #Tuple of complex structures
        "RGBA" => ARGB{N0f8},
        "EXTREMA" => Extrema,
        "LENGTH" => Length,
        "ARRAY" => Array, #Dict won't allow Array to be key in HDF5PLOT_MAP_TELEM2STR

        #Sub-structure types:
        "FONT" => Font,
        "BOUNDINGBOX" => BoundingBox,
        "GRIDLAYOUT" => GridLayout,
        "ROOTLAYOUT" => RootLayout,
        "SERIESANNOTATIONS" => SeriesAnnotations,
        "PLOTTEXT" => PlotText,
        "SHAPE" => Shape,
        "COLORGRADIENT" => ColorGradient,
        "AXIS" => Axis,
        "SURFACE" => Surface,
        "SUBPLOT" => Subplot,
        "NULLABLE" => _hdf5_nullable,
    )
    merge!(HDF5PLOT_MAP_STR2TELEM, telem2str)
    merge!(HDF5PLOT_MAP_TELEM2STR, Dict{Type, String}(v=>k for (k,v) in HDF5PLOT_MAP_STR2TELEM))
end

#==Helper functions
===============================================================================#

_hdf5_plotelempath(subpath::String) = "$_hdf5_plotroot/$subpath"
_hdf5_datapath(subpath::String) = "$_hdf5_dataroot/$subpath"
_hdf5_map_str2telem(k::String) = HDF5PLOT_MAP_STR2TELEM[k]
_hdf5_map_str2telem(v::Vector) = HDF5PLOT_MAP_STR2TELEM[v[1]]

function _hdf5_merge!(dest::Dict, src::Dict)
    for (k, v) in src
        if isa(v, Axis)
            _hdf5_merge!(dest[k].plotattributes, v.plotattributes)
        else
            dest[k] = v
        end
    end
    return
end

#==
===============================================================================#

# Create the window/figure for this backend.
function _create_backend_figure(plt::Plot{HDF5Backend})
    #Do nothing
end

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

# # this is called early in the pipeline, use it to make the plot current or something
# function _prepare_plot_object(plt::Plot{HDF5Backend})
# end

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

# Set up the subplot within the backend object.
function _initialize_subplot(plt::Plot{HDF5Backend}, sp::Subplot{HDF5Backend})
    #Do nothing
end

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

# Add one series to the underlying backend object.
# Called once per series
# NOTE: Seems to be called when user calls plot()... even if backend
#       plot, sp.o has not yet been constructed...
function _series_added(plt::Plot{HDF5Backend}, series::Series)
    #Do nothing
end

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

# When series data is added/changed, this callback can do dynamic updates to the backend object.
# note: if the backend rebuilds the plot from scratch on display, then you might not do anything here.
function _series_updated(plt::Plot{HDF5Backend}, series::Series)
    #Do nothing
end

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

# called just before updating layout bounding boxes... in case you need to prep
# for the calcs
function _before_layout_calcs(plt::Plot{HDF5Backend})
    #Do nothing
end

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

# Set the (left, top, right, bottom) minimum padding around the plot area
# to fit ticks, tick labels, guides, colorbars, etc.
function _update_min_padding!(sp::Subplot{HDF5Backend})
    #Do nothing
end

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

# Override this to update plot items (title, xlabel, etc), and add annotations (plotattributes[:annotations])
function _update_plot_object(plt::Plot{HDF5Backend})
    #Do nothing
end

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

# Display/show the plot (open a GUI window, or browser page, for example).
function _display(plt::Plot{HDF5Backend})
    msg = "HDF5 interface does not support `display()` function."
    msg *= "\nUse `Plots.hdf5plot_write(::String)` method to write to .HDF5 \"plot\" file instead."
    @warn(msg)
    return
end

#==HDF5 write functions
===============================================================================#

function _hdf5plot_writetype(grp, k::String, tstr::Array{String})
    d = HDF5.d_open(grp, k)
    HDF5.a_write(d, _hdf5plot_datatypeid, tstr)
end
function _hdf5plot_writetype(grp, k::String, T::Type)
    tstr = HDF5PLOT_MAP_TELEM2STR[T]
    d = HDF5.d_open(grp, k)
    HDF5.a_write(d, _hdf5plot_datatypeid, tstr)
end
function _hdf5plot_overwritetype(grp, k::String, T::Type)
    tstr = HDF5PLOT_MAP_TELEM2STR[T]
    d = HDF5.d_open(grp, k)
    HDF5.a_delete(d, _hdf5plot_datatypeid)
    HDF5.a_write(d, _hdf5plot_datatypeid, tstr)
end
function _hdf5plot_writetype(grp, T::Type) #Write directly to group
    tstr = HDF5PLOT_MAP_TELEM2STR[T]
    HDF5.a_write(grp, _hdf5plot_datatypeid, tstr)
end
function _hdf5plot_overwritetype(grp, T::Type) #Write directly to group
    tstr = HDF5PLOT_MAP_TELEM2STR[T]
    HDF5.a_delete(grp, _hdf5plot_datatypeid)
    HDF5.a_write(grp, _hdf5plot_datatypeid, tstr)
end
function _hdf5plot_writetype(grp, ::Type{Array{T}}) where T<:Any
    tstr = HDF5PLOT_MAP_TELEM2STR[Array] #ANY
    HDF5.a_write(grp, _hdf5plot_datatypeid, tstr)
end
function _hdf5plot_writetype(grp, ::Type{T}) where T<:BoundingBox
    tstr = HDF5PLOT_MAP_TELEM2STR[BoundingBox]
    HDF5.a_write(grp, _hdf5plot_datatypeid, tstr)
end
function _hdf5plot_writecount(grp, n::Int) #Write directly to group
    HDF5.a_write(grp, _hdf5plot_countid, n)
end
function _hdf5plot_gwritefields(grp, k::String, v)
    grp = HDF5.g_create(grp, k)
    for _k in fieldnames(typeof(v))
        _v = getfield(v, _k)
        kstr = string(_k)
        _hdf5plot_gwrite(grp, kstr, _v)
    end
    _hdf5plot_writetype(grp, typeof(v))
    return
end

# Write data
# ----------------------------------------------------------------

function _hdf5plot_gwrite(grp, k::String, v) #Default
    T = typeof(v)
    if !(T <: Number || T <: String)
        tstr = string(T)
        path = HDF5.name(grp) * "/" * k
        @info("Type not supported: $tstr\npath: $path")
#        @show v
        return
    end
    grp[k] = v
    _hdf5plot_writetype(grp, k, HDF5PlotNative)
end
function _hdf5plot_gwrite(grp, k::String, v::Array{T}) where T<:Number #Default for arrays
    grp[k] = v
    _hdf5plot_writetype(grp, k, HDF5PlotNative)
end
#=
function _hdf5plot_gwrite(grp, k::String, v::Array{Any})
#    @show grp, k
    @warn("Cannot write Array: $k=$v")
end
=#
function _hdf5plot_gwrite(grp, k::String, v::Nothing)
    grp[k] = 0
    _hdf5plot_writetype(grp, k, Nothing)
end
function _hdf5plot_gwrite(grp, k::String, v::Bool)
    grp[k] = Int(v)
    _hdf5plot_writetype(grp, k, Bool)
end
function _hdf5plot_gwrite(grp, k::String, v::Symbol)
    grp[k] = string(v)
    _hdf5plot_writetype(grp, k, Symbol)
end
function _hdf5plot_gwrite(grp, k::String, v::Tuple)
    varr = [v...]
    elt = eltype(varr)
#    if isleaftype(elt)

    _hdf5plot_gwrite(grp, k, varr)
    if elt <: Number
        #We just wrote a simple dataset
        _hdf5plot_overwritetype(grp, k, Tuple)
    else #Used a more complex scheme (using subgroups):
        _hdf5plot_overwritetype(grp[k], HDF5CTuple)
    end
    #NOTE: _hdf5plot_overwritetype overwrites "Array" type with "Tuple".
end
function _hdf5plot_gwrite(grp, k::String, plotattributes::Dict)
#    @warn("Cannot write dict: $k=$plotattributes")
end
function _hdf5plot_gwrite(grp, k::String, v::AbstractRange)
    _hdf5plot_gwrite(grp, k, collect(v)) #For now
end
function _hdf5plot_gwrite(grp, k::String, v::ARGB{N0f8})
    grp[k] = [v.r.i, v.g.i, v.b.i, v.alpha.i]
    _hdf5plot_writetype(grp, k, ARGB{N0f8})
end
function _hdf5plot_gwrite(grp, k::String, v::Colorant)
    _hdf5plot_gwrite(grp, k, ARGB{N0f8}(v))
end
#Custom vector (when not using simple numeric type):
function _hdf5plot_gwritearray(grp, k::String, v::Array{T}) where T
    vgrp = HDF5.g_create(grp, k)
    _hdf5plot_writetype(vgrp, Array) #ANY
    sz = size(v)
    lidx = LinearIndices(sz)

    for iter in eachindex(v)
        coord = lidx[iter]
        elem = v[iter]
        idxstr = join(coord, "_")
        _hdf5plot_gwrite(vgrp, "v$idxstr", elem)
    end

    _hdf5plot_gwrite(vgrp, "dim", [sz...])
    return
end
_hdf5plot_gwrite(grp, k::String, v::Array) =
	_hdf5plot_gwritearray(grp, k, v)
function _hdf5plot_gwrite(grp, k::String, v::Extrema)
    grp[k] = [v.emin, v.emax]
    _hdf5plot_writetype(grp, k, Extrema)
end
function _hdf5plot_gwrite(grp, k::String, v::Length{T}) where T
    grp[k] = v.value
    _hdf5plot_writetype(grp, k, [HDF5PLOT_MAP_TELEM2STR[Length], string(T)])
end

# Write more complex structures:
# ----------------------------------------------------------------

function _hdf5plot_gwrite(grp, k::String, v::Plot)
    #Don't write plot references
end
function _hdf5plot_gwrite(grp, k::String, v::HDF5PLOT_SIMPLESUBSTRUCT)
    _hdf5plot_gwritefields(grp, k, v)
    return
end
function _hdf5plot_gwrite(grp, k::String, v::Axis)
    grp = HDF5.g_create(grp, k)
    for (_k, _v) in v.plotattributes
        kstr = string(_k)
        _hdf5plot_gwrite(grp, kstr, _v)
    end
    _hdf5plot_writetype(grp, Axis)
    return
end
function _hdf5plot_gwrite(grp, k::String, v::Surface)
	grp = HDF5.g_create(grp, k)
	_hdf5plot_gwrite(grp, "data2d", v.surf)
	_hdf5plot_writetype(grp, Surface)
end
# #TODO: "Properly" support Nullable using _hdf5plot_writetype?
# function _hdf5plot_gwrite(grp, k::String, v::_hdf5_nullable)
#     if isnull(v)
#         _hdf5plot_gwrite(grp, k, nothing)
#     else
#         _hdf5plot_gwrite(grp, k, v.value)
#     end
#     return
# end

function _hdf5plot_gwrite(grp, k::String, v::Subplot)
    grp = HDF5.g_create(grp, k)
    _hdf5plot_gwrite(grp, "index", v[:subplot_index])
    _hdf5plot_writetype(grp, Subplot)
    return
end
function _hdf5plot_write(grp, plotattributes::Dict)
    for (k, v) in plotattributes
        kstr = string(k)
        _hdf5plot_gwrite(grp, kstr, v)
    end
    return
end

# Write main plot structures:
# ----------------------------------------------------------------

function _hdf5plot_write(sp::Subplot{HDF5Backend}, subpath::String, f)
    f = f::HDF5.HDF5File #Assert
    grp = HDF5.g_create(f, _hdf5_plotelempath("$subpath/attr"))
    _hdf5plot_write(grp, sp.attr)
    grp = HDF5.g_create(f, _hdf5_plotelempath("$subpath/series_list"))
    _hdf5plot_writecount(grp, length(sp.series_list))
    for (i, series) in enumerate(sp.series_list)
        grp = HDF5.g_create(f, _hdf5_plotelempath("$subpath/series_list/series$i"))
        _hdf5plot_write(grp, series.plotattributes)
    end

    return
end

function _hdf5plot_write(plt::Plot{HDF5Backend}, f)
    f = f::HDF5.HDF5File #Assert

    grp = HDF5.g_create(f, _hdf5_plotelempath("attr"))
    _hdf5plot_write(grp, plt.attr)

    grp = HDF5.g_create(f, _hdf5_plotelempath("subplots"))
    _hdf5plot_writecount(grp, length(plt.subplots))

    for (i, sp) in enumerate(plt.subplots)
        _hdf5plot_write(sp, "subplots/subplot$i", f)
    end

    return
end
function hdf5plot_write(plt::Plot{HDF5Backend}, path::AbstractString)
    HDF5.h5open(path, "w") do file
        _hdf5plot_write(plt, file)
    end
end
hdf5plot_write(path::AbstractString) = hdf5plot_write(current(), path)

#==HDF5 playback (read) functions
===============================================================================#

function _hdf5plot_readcount(grp) #Read directly from group
    return HDF5.a_read(grp, _hdf5plot_countid)
end

_hdf5plot_convert(T::Type{HDF5PlotNative}, v) = v
_hdf5plot_convert(T::Type{Nothing}, v) = nothing
_hdf5plot_convert(T::Type{Bool}, v) = (v!=0)
_hdf5plot_convert(T::Type{Symbol}, v) = Symbol(v)
_hdf5plot_convert(T::Type{Tuple}, v) = tuple(v...) #With Vector{T<:Number}
function _hdf5plot_convert(T::Type{ARGB{N0f8}}, v)
    r, g, b, a = reinterpret(N0f8, v)
    return Colors.ARGB{N0f8}(r, g, b, a)
end
_hdf5plot_convert(T::Type{Extrema}, v) = Extrema(v[1], v[2])

# Read data structures:
# ----------------------------------------------------------------

function _hdf5plot_read(grp, k::String, T::Type, dtid)
    v = HDF5.d_read(grp, k)
    return _hdf5plot_convert(T, v)
end
function _hdf5plot_read(grp, k::String, T::Type{Length}, dtid::Vector)
    v = HDF5.d_read(grp, k)
    TU = Symbol(dtid[2])
    T = typeof(v)
    return Length{TU,T}(v)
end

# Read more complex data structures:
# ----------------------------------------------------------------
function _hdf5plot_read(grp, k::String, T::Type{Font}, dtid)
    grp = HDF5.g_open(grp, k)

    family = _hdf5plot_read(grp, "family")
    pointsize = _hdf5plot_read(grp, "pointsize")
    halign = _hdf5plot_read(grp, "halign")
    valign = _hdf5plot_read(grp, "valign")
    rotation = _hdf5plot_read(grp, "rotation")
    color = _hdf5plot_read(grp, "color")
    return Font(family, pointsize, halign, valign, rotation, color)
end
function _hdf5plot_read(grp, k::String, T::Type{Array}, dtid) #ANY
    grp = HDF5.g_open(grp, k)
    sz = _hdf5plot_read(grp, "dim")
    if [0] == sz; return []; end
    sz = tuple(sz...)
    result = Array{Any}(undef, sz)
    lidx = LinearIndices(sz)

    for iter in eachindex(result)
        coord = lidx[iter]
        idxstr = join(coord, "_")
        result[iter] = _hdf5plot_read(grp, "v$idxstr")
    end

    #Hack: Implicitly make Julia detect element type.
    #      (Should probably write it explicitly to file)
    result = [result[iter] for iter in eachindex(result)] #Potentially make more specific
    return reshape(result, sz)
end
function _hdf5plot_read(grp, k::String, T::Type{HDF5CTuple}, dtid)
    v = _hdf5plot_read(grp, k, Array, dtid)
    return tuple(v...)
end
function _hdf5plot_read(grp, k::String, T::Type{PlotText}, dtid)
    grp = HDF5.g_open(grp, k)

    str = _hdf5plot_read(grp, "str")
    font = _hdf5plot_read(grp, "font")
    return PlotText(str, font)
end
function _hdf5plot_read(grp, k::String, T::Type{SeriesAnnotations}, dtid)
    grp = HDF5.g_open(grp, k)

    strs = _hdf5plot_read(grp, "strs")
    font = _hdf5plot_read(grp, "font")
    baseshape = _hdf5plot_read(grp, "baseshape")
    scalefactor = _hdf5plot_read(grp, "scalefactor")
    return SeriesAnnotations(strs, font, baseshape, scalefactor)
end
function _hdf5plot_read(grp, k::String, T::Type{Shape}, dtid)
    grp = HDF5.g_open(grp, k)

    x = _hdf5plot_read(grp, "x")
    y = _hdf5plot_read(grp, "y")
    return Shape(x, y)
end
function _hdf5plot_read(grp, k::String, T::Type{ColorGradient}, dtid)
    grp = HDF5.g_open(grp, k)

    colors = _hdf5plot_read(grp, "colors")
    values = _hdf5plot_read(grp, "values")
    return ColorGradient(colors, values)
end
function _hdf5plot_read(grp, k::String, T::Type{BoundingBox}, dtid)
    grp = HDF5.g_open(grp, k)

    x0 = _hdf5plot_read(grp, "x0")
    a = _hdf5plot_read(grp, "a")
    return BoundingBox(x0, a)
end
_hdf5plot_read(grp, k::String, T::Type{RootLayout}, dtid) = RootLayout()
function _hdf5plot_read(grp, k::String, T::Type{GridLayout}, dtid)
    grp = HDF5.g_open(grp, k)

#    parent = _hdf5plot_read(grp, "parent")
parent = RootLayout()
    minpad = _hdf5plot_read(grp, "minpad")
    bbox = _hdf5plot_read(grp, "bbox")
    grid = _hdf5plot_read(grp, "grid")
    widths = _hdf5plot_read(grp, "widths")
    heights = _hdf5plot_read(grp, "heights")
    attr = KW() #TODO support attr: _hdf5plot_read(grp, "attr")

    return GridLayout(parent, minpad, bbox, grid, widths, heights, attr)
end
function _hdf5plot_read(grp, k::String, T::Type{Axis}, dtid)
    grp = HDF5.g_open(grp, k)
    kwlist = KW()
    _hdf5plot_read(grp, kwlist)
    return Axis([], kwlist)
end
function _hdf5plot_read(grp, k::String, T::Type{Surface}, dtid)
    grp = HDF5.g_open(grp, k)
    data2d = _hdf5plot_read(grp, "data2d")
    return Surface(data2d)
end
function _hdf5plot_read(grp, k::String, T::Type{Subplot}, dtid)
    grp = HDF5.g_open(grp, k)
    idx = _hdf5plot_read(grp, "index")
    return HDF5PLOT_PLOTREF.ref.subplots[idx]
end
function _hdf5plot_read(grp, k::String)
    dtid = HDF5.a_read(grp[k], _hdf5plot_datatypeid)
    T = _hdf5_map_str2telem(dtid) #expect exception
    return _hdf5plot_read(grp, k, T, dtid)
end

#Read in values in group to populate plotattributes:
function _hdf5plot_read(grp, plotattributes::Dict)
    gnames = names(grp)
    for k in gnames
        try
            v = _hdf5plot_read(grp, k)
            plotattributes[Symbol(k)] = v
        catch e
            @show e
            @show grp
            @warn("Could not read field $k")
        end
    end
    return
end

# Read main plot structures:
# ----------------------------------------------------------------

function _hdf5plot_read(sp::Subplot, subpath::String, f)
    f = f::HDF5.HDF5File #Assert

    grp = HDF5.g_open(f, _hdf5_plotelempath("$subpath/series_list"))
    nseries = _hdf5plot_readcount(grp)

    for i in 1:nseries
        grp = HDF5.g_open(f, _hdf5_plotelempath("$subpath/series_list/series$i"))
        kwlist = KW()
        _hdf5plot_read(grp, kwlist)
        plot!(sp, kwlist[:x], kwlist[:y]) #Add data & create data structures
        _hdf5_merge!(sp.series_list[end].plotattributes, kwlist)
    end

    #Perform after adding series... otherwise values get overwritten:
    grp = HDF5.g_open(f, _hdf5_plotelempath("$subpath/attr"))
    kwlist = KW()
    _hdf5plot_read(grp, kwlist)
    _hdf5_merge!(sp.attr, kwlist)

    return
end

function _hdf5plot_read(plt::Plot, f)
    f = f::HDF5.HDF5File #Assert
    #Assumpltion: subplots are already allocated (plt.subplots)

    HDF5PLOT_PLOTREF.ref = plt #Used when reading "layout"
    grp = HDF5.g_open(f, _hdf5_plotelempath("attr"))
    _hdf5plot_read(grp, plt.attr)

    for (i, sp) in enumerate(plt.subplots)
        _hdf5plot_read(sp, "subplots/subplot$i", f)
    end

    return
end

function hdf5plot_read(path::AbstractString)
    plt = nothing
    HDF5.h5open(path, "r") do file
        grp = HDF5.g_open(file, _hdf5_plotelempath("subplots"))
        n = _hdf5plot_readcount(grp)
        plt = plot(layout=n) #Get reference to a new plot
        _hdf5plot_read(plt, file)
    end
    return plt
end

#Last line


text:
using SimpleSolvers
using Test

n = 2
x = rand(n)
g = 2x
T = eltype(x)

function F(x::Vector)
    1 + sum(x.^2)
end

function ∇F!(g::Vector, x::Vector)
    g .= 0
    for i in eachindex(x,g)
        g[i] = 2x[i]
    end
end

∇PAD = GradientParameters{T}(F, n; mode = :autodiff, diff_type = :forward)
∇PFD = GradientParameters{T}(F, n; mode = :autodiff, diff_type = :finite)
∇PUS = GradientParameters{T}(∇F!, n; mode = :user)

@test typeof(∇PAD) <: GradientParametersAD
@test typeof(∇PFD) <: GradientParametersFD
@test typeof(∇PUS) <: GradientParametersUser

function test_grad(g1, g2, atol)
    for i in eachindex(g1,g2)
        @test g1[i] ≈ g2[i] atol=atol
    end
end

gad = zero(g)
gfd = zero(g)
gus = zero(g)

compute_gradient!(gad, x, ∇PAD)
compute_gradient!(gfd, x, ∇PFD)
compute_gradient!(gus, x, ∇PUS)

test_grad(gad, g, eps())
test_grad(gfd, g, 1E-7)
test_grad(gus, g, 0)

gad1 = zero(g)
gfd1 = zero(g)
gus1 = zero(g)

compute_gradient!(gad1, x, F; mode = :autodiff, diff_type = :forward)
compute_gradient!(gfd1, x, F; mode = :autodiff, diff_type = :finite)
compute_gradient!(gus1, x, ∇F!; mode = :user)

test_grad(gad, gad1, 0)
test_grad(gfd, gfd1, 0)
test_grad(gus, gus1, 0)

gad2 = zero(g)
gfd2 = zero(g)

compute_gradient_ad!(gad2, x, F)
compute_gradient_fd!(gfd2, x, F)

test_grad(gad, gad2, 0)
test_grad(gfd, gfd2, 0)


text:
export DDPGPolicy

using Random
using Flux

mutable struct DDPGPolicy{
    BA<:NeuralNetworkApproximator,
    BC<:NeuralNetworkApproximator,
    TA<:NeuralNetworkApproximator,
    TC<:NeuralNetworkApproximator,
    P,
    R<:AbstractRNG,
} <: AbstractPolicy

    behavior_actor::BA
    behavior_critic::BC
    target_actor::TA
    target_critic::TC
    γ::Float32
    ρ::Float32
    batch_size::Int
    start_steps::Int
    start_policy::P
    update_after::Int
    update_every::Int
    act_limit::Float64
    act_noise::Float64
    step::Int
    rng::R
end

"""
    DDPGPolicy(;kwargs...)

# Keyword arguments

- `behavior_actor`,
- `behavior_critic`,
- `target_actor`,
- `target_critic`,
- `start_policy`,
- `γ = 0.99f0`,
- `ρ = 0.995f0`,
- `batch_size = 32`,
- `start_steps = 10000`,
- `update_after = 1000`,
- `update_every = 50`,
- `act_limit = 1.0`,
- `act_noise = 0.1`,
- `step = 0`,
- `seed = nothing`,
"""
function DDPGPolicy(;
    behavior_actor,
    behavior_critic,
    target_actor,
    target_critic,
    start_policy,
    γ = 0.99f0,
    ρ = 0.995f0,
    batch_size = 32,
    start_steps = 10000,
    update_after = 1000,
    update_every = 50,
    act_limit = 1.0,
    act_noise = 0.1,
    step = 0,
    seed = nothing,
)
    rng = MersenneTwister(seed)
    copyto!(behavior_actor, target_actor)  # force sync
    copyto!(behavior_critic, target_critic)  # force sync
    DDPGPolicy(
        behavior_actor,
        behavior_critic,
        target_actor,
        target_critic,
        γ,
        ρ,
        batch_size,
        start_steps,
        start_policy,
        update_after,
        update_every,
        act_limit,
        act_noise,
        step,
        rng,
    )
end

# TODO: handle Training/Testing mode
function (p::DDPGPolicy)(obs)
    p.step += 1

    if p.step <= p.start_steps
        p.start_policy(obs)
    else
        D = device(p.behavior_actor)
        s = get_state(obs)
        s = Flux.unsqueeze(s, ndims(s) + 1)
        action = p.behavior_actor(send_to_device(D, s)) |> vec |> send_to_host
        clamp(action[] + randn(p.rng) * p.act_noise, -p.act_limit, p.act_limit)
    end
end

function RLBase.update!(p::DDPGPolicy, t::CircularCompactSARTSATrajectory)
    length(t) > p.update_after || return
    p.step % p.update_every == 0 || return

    inds = rand(p.rng, 1:(length(t)-1), p.batch_size)
    SARTS = (:state, :action, :reward, :terminal, :next_state)
    s, a, r, t, s′ = map(x -> select_last_dim(get_trace(t, x), inds), SARTS)

    A = p.behavior_actor
    C = p.behavior_critic
    Aₜ = p.target_actor
    Cₜ = p.target_critic

    γ = p.γ
    ρ = p.ρ

    # !!! we have several assumptions here, need revisit when we have more complex environments
    # state is vector
    # action is scalar
    a′ = Aₜ(s′)
    qₜ = Cₜ(vcat(s′, a′)) |> vec
    y = r .+ γ .* (1 .- t) .* qₜ
    a = Flux.unsqueeze(a, 1)

    gs1 = gradient(Flux.params(C)) do
        q = C(vcat(s, a)) |> vec
        mean((y .- q) .^ 2)
    end

    update!(C, gs1)

    gs2 = gradient(Flux.params(A)) do
        -mean(C(vcat(s, A(s))))
    end

    update!(A, gs2)

    # polyak averaging
    for (dest, src) in zip(Flux.params([Aₜ, Cₜ]), Flux.params([A, C]))
        dest .= ρ .* dest .+ (1 - ρ) .* src
    end
end


text:
"""

# `check_open`

$(SIGNATURES)

Internal.
"""
function check_open(d::DAG, path::Vector{Symbol}, conditioning_set::SymbolList; debug=false)
  res = :open
  path_closed = false
  if typeof(conditioning_set) == Symbol
    conditioning_set = [conditioning_set]
  end
  r = range(2, stop=length(path)-1)
  for i in r

    # If this vertex a collider, i.e. path[i-1] -> path[i] <- path[i+1], and
    # 1. not in the conditioning_set or
    # 2. has does not have a causal descendent that is in the conditioning set
    # then the path is blocked

    vertex_type = :unknown

    if d.a[path[i], path[i-1]] == 1 && d.a[path[i], path[i+1]] == 1
      vertex_type = :collider
      if length(conditioning_set) > 0
        path_closed = path[i] in conditioning_set ? false : true
      else
        path_closed = true
      end
    elseif d.a[path[i-1], path[i]] == 1 && d.a[path[i+1], path[i]] == 1
      vertex_type = :fork
      if length(conditioning_set) > 0
        path_closed = path[i] in conditioning_set ? true : false
      else
        path_closed = false
      end
    elseif d.a[path[i-1], path[i]] == 1 && d.a[path[i], path[i+1]] == 1 ||
        d.a[path[i], path[i-1]] == 1 && d.a[path[i+1], path[i]] == 1
      vertex_type = :chain
      if length(conditioning_set) > 0
        path_closed = path[i] in conditioning_set ? true : false
      else
        path_closed = false
      end
    else
      @info "Vertex_type of $(path[i]) = $(vertex_type)\n"*
        ""
    end
    debug && println("Vertex_type of $(path[i]) in $path = $(vertex_type)\n")

    # If path_closed == true, check that no causal descendent
    # of path[i] is in conditioning set. If that is the case,
    # flip the state. Max_length of decendent path? What if forks
    # in the path?

    res = path_closed ? :closed : :open
    res == :closed && break
  end
  res
end

"""

# `open_paths`

$(SIGNATURES)

Internal.
"""
function open_paths(d::DAG, paths::Vector{Vector{Symbol}}, cond::Vector{Symbol};debug=false)
  open = fill(:closed, length(paths))
  for (ind, path) in enumerate(paths)
    if length(path) == 2
      if d.a[path[2], path[1]] == 1
        open[ind] = :open
      end
    else
      open[ind] = check_open(d, path, cond)
    end
  end
  ind = findall(x -> x == :open, open)
  paths[ind]
end

export
  open_paths,
  check_open


text:
export allelefreq, allelefreq_vec, avg_allelefreq

"""
    allelefreq(data::PopData)
Return a NamedTuple of `Dicts` of allele frequencies of all
loci in a `PopData` object. These are global allele frequencies.
"""
function allelefreq(data::PopData)
    NamedTuple{Tuple(Symbol.(loci(data)))}(
            DataFrames.combine(
            groupby(data.genodata, :locus),
            :genotype => allelefreq => :frq
        )[:, :frq] |> Tuple
    )
end

"""
    allelefreq(locus::GenoArray)
Return a `Dict` of allele frequencies of a GenoArray (typically a single locus) in
a `PopData` object.
"""
@inline function allelefreq(locus::GenoArray)
    all(ismissing.(locus)) == true && return Dict{eltype(nonmissingtype(eltype(locus))), Float64}()
    proportionmap(alleles(locus))
end

"""
    allelefreq(geno::Genotype)
Return a `Dict` of allele frequencies of the alleles within a single Genotype in a `PopData`
object.
"""
@inline function allelefreq(geno::Genotype)
    d = Dict{eltype(geno),Float32}()
    len = length(geno)
    @inbounds @simd for allele in geno
        d[allele] = @inbounds get!(d, allele, 0.0) + 1.0/len
    end
    return d
end

"""
    allelefreq_vec(locus::GenoArray)
Return a Vector of allele frequencies of a single locus in a `PopData`
object. Similar to `allelefreq()`, except it returns only the frequencies,
without the allele names, meaning they can be in any order. This is useful
for getting the expected genotype frequencies.
"""
@inline function allelefreq_vec(locus::GenoArray)
    flat_alleles = alleles(locus)
    len = length(flat_alleles)
    d = [count(i -> i == j, flat_alleles) for j in unique(flat_alleles)]
    return d ./ len
end

@inline function allelefreq_vec(::Missing)
    return missing
end

"""
    avg_allelefreq(allele_dicts::AbstractVector{Dict{T, Float64}}, power::Int = 1)
Takes a Vector of Dicts generated by `allelefreq` and returns a Dict of the average
allele frequencies raised to the `power` (exponent) specified (default: `1`).
This is typically done to calculate average allele frequencies across populations.

**Example**
```
cats = @nancycats;
alleles_df = DataFrames.combine(
    groupby(cats.genodata, [:locus, :population]),
    :genotype => allelefreq => :alleles
);
DataFrames.combine(
    groupby(alleles_df, :locus),
    :alleles => (i -> sum(avg_allelefreq(i, 2))) => :avg_freq
)
```
"""
function avg_allelefreq(allele_dicts::AbstractVector{Dict{T, Float64}}, power::Int = 1) where T<:Integer
   sum_dict = Dict{Int16, Tuple{Float32, Int}}()
   # remove any dicts with no entries (i.e. from a group without that locus)
   allele_dicts = allele_dicts[findall(i -> length(i) > 0, allele_dicts)]
   # create a list of all the alleles
   all_alleles = keys.(allele_dicts) |> Base.Iterators.flatten |> collect |> unique
   # populate the sum dict with allele frequency and n for each allele
   @inbounds for allele in all_alleles
       for allele_dict in allele_dicts
           @inbounds sum_dict[allele] = get!(sum_dict, allele, (0., 0)) .+ (get!(allele_dict, allele, 0.), 1)
       end
   end
   avg_dict = Dict{Int16, Float32}()
   # collapse the sum dict into a dict of averages
   @inbounds for (key, value) in sum_dict
       freq_sum, n = value
       avg = (freq_sum / n) ^ power
       # drop zeroes
       if !iszero(avg)
           @inbounds avg_dict[key] = avg
       end
   end
   return avg_dict
end

# method for nei_fst (pairwise)
function avg_allelefreq(allele_dicts::T, power::Int = 1) where T<:Tuple
    sum_dict = Dict{Int16, Tuple{Float32, Int}}()
    # remove any dicts with no entries (i.e. from a group without that locus)
    allele_dicts = allele_dicts[findall(i -> length(i) > 0, allele_dicts)]
    # create a list of all the alleles
    all_alleles = keys.(allele_dicts) |> Base.Iterators.flatten |> collect |> unique
    # populate the sum dict with allele frequency and n for each allele
    @inbounds for allele in all_alleles
        for allele_dict in allele_dicts
            @inbounds sum_dict[allele] = get!(sum_dict, allele, (0., 0)) .+ (get!(allele_dict, allele, 0.), 1)
        end
    end
    avg_dict = Dict{Int16, Float32}()
    # collapse the sum dict into a dict of averages
    @inbounds for (key, value) in sum_dict
        freq_sum, n = value
        avg = (freq_sum / n) ^ power
        # drop zeroes
        if !iszero(avg)
            @inbounds avg_dict[key] = avg
        end
    end
    return avg_dict
 end

"""
    allelefreq(data::PopData, locus::String; population::Bool = false)
Return a `Dict` of allele frequencies of a single locus in a `PopData`
object. Use `population = true` to return a table of allele frequencies
of that locus per population.
### Example
```
cats = @nancycats;
allelefreq(cats, "fca8")
allelefreq(cats, "fca8", population = true)
```
"""
function allelefreq(data::PopData, locus::String; population::Bool=false)
    if !population
        data.genodata[data.genodata.locus .== locus, :genotype] |> allelefreq
    else
        tmp = groupby(data.genodata[data.genodata.locus .== locus, :], :population)
        DataFrames.combine(tmp, :genotype => allelefreq => :frequency)
    end
end

#TODO swtich order of args do it's data, allele?
# Does doing that break anything?
"""
    allelefreq(allele::Int, genos::GenoArray)
Return the frequency of a specific `allele` from a vector of `genotypes`.

### Example
```
using DataFramesMeta
ncats = @nancycats;
ncats_sub @where(ncats.genodata, :locus .== "fca8", :genotype .!== missing)
pop_grp = groupby(ncats_sub, :population)
DataFrames.combine(pop_grp, :genotype => (geno -> allelefreq(137, geno)) => :freq_137)
```
"""
function allelefreq(allele::Int, genos::T) where T<:GenoArray
    tmp = allelefreq(genos)
    haskey(tmp, allele) ? getindex(tmp, allele) : 0.
end


text:
type CBFData
    name::String
    sense::Symbol
    var::Vector{Tuple{String,Int}}
    psdvar::Vector{Int}
    con::Vector{Tuple{String,Int}}
    psdcon::Vector{Int}
    objacoord::Vector{Tuple{Int,Float64}}
    objfcoord::Vector{Tuple{Int,Int,Int,Float64}}
    objoffset::Float64
    fcoord::Vector{Tuple{Int,Int,Int,Int,Float64}}
    acoord::Vector{Tuple{Int,Int,Float64}} # linear coefficients
    bcoord::Vector{Tuple{Int,Float64}} # linear offsets
    hcoord::Vector{Tuple{Int,Int,Int,Int,Float64}}
    dcoord::Vector{Tuple{Int,Int,Int,Float64}}
    intlist::Vector{Int}
    nvar::Int
    nconstr::Int
end

CBFData() = CBFData("",:xxx,[],[],[],[],[],[],0.0,[],[],[],[],[],[],0,0)

function parse_matblock(fd,outputmat,num_indices)
    nextline = readline(fd)
    nnz = parse(Int,strip(nextline))
    for k in 1:nnz
        nextline = readline(fd)
        tup = split(strip(nextline))
        push!(outputmat, (map(s->parse(Int,s)+1,tup[1:num_indices])...,parse(Float64,tup[end])))
    end
end

function readcbfdata(filename)
    if endswith(filename,"cbf.gz")
        fd = gzopen(filename,"r")
    else
        @assert endswith(filename, "cbf")
        fd = open(filename,"r")
    end

    dat = CBFData()
    dat.name = split(basename(filename),".")[1]

    while !eof(fd)
        line = readline(fd)
        startswith(line,"#") && continue # comments
        length(line) == 1 && continue # blank lines

        # new block

        if startswith(line,"VER")
            nextline = readline(fd)
            @assert startswith(nextline,"1") || startswith(nextline,"2")
            continue
        end

        if startswith(line,"OBJSENSE")
            nextline = readline(fd)
            if strip(nextline) == "MIN"
                dat.sense = :Min
            else
                dat.sense = :Max
            end
            continue
        end

        if startswith(line,"VAR")
            nextline = readline(fd)
            totalvars, lines = split(nextline)
            totalvars = parse(Int,strip(totalvars))
            lines = parse(Int,strip(lines))
            varcnt = 0

            for k in 1:lines
                nextline = readline(fd)
                cone, sz = split(nextline)
                sz = parse(Int,strip(sz))
                push!(dat.var, (cone, sz))
                varcnt += sz
            end
            @assert totalvars == varcnt
            dat.nvar = varcnt
            continue
        end

        if startswith(line, "INT")
            nextline = readline(fd)
            intvar = parse(Int,strip(nextline))
            for k in 1:intvar
                nextline = readline(fd)
                idx = parse(Int,strip(nextline))
                push!(dat.intlist,idx+1)
            end
            continue
        end

        if startswith(line,"CON")
            nextline = readline(fd)
            totalconstr, lines = split(nextline)
            totalconstr = parse(Int,strip(totalconstr))
            lines = parse(Int,strip(lines))
            constrcnt = 0

            for k in 1:lines
                nextline = readline(fd)
                cone, sz = split(nextline)
                sz = parse(Int,strip(sz))
                push!(dat.con, (cone, sz))
                constrcnt += sz
            end
            @assert totalconstr == constrcnt
            dat.nconstr = constrcnt
            continue
        end

        if startswith(line,"PSDVAR")
            nextline = readline(fd)
            lines = parse(Int,strip(nextline))

            for k in 1:lines
                nextline = readline(fd)
                sz = parse(Int,strip(nextline))
                push!(dat.psdvar, sz)
            end
            continue
        end

        if startswith(line,"PSDCON")
            nextline = readline(fd)
            lines = parse(Int,strip(nextline))

            for k in 1:lines
                nextline = readline(fd)
                sz = parse(Int,strip(nextline))
                push!(dat.psdcon, sz)
            end
            continue
        end

        if startswith(line,"OBJACOORD")
            parse_matblock(fd,dat.objacoord,1)
        end

        if startswith(line,"OBJBCOORD")
            nextline = readline(fd)
            dat.objoffset = float(strip(nextline))
            warn("Instance has objective offset")
        end

        if startswith(line,"BCOORD")
            parse_matblock(fd,dat.bcoord,1)
        end

        if startswith(line,"ACOORD")
            parse_matblock(fd,dat.acoord,2)
        end

        if startswith(line,"OBJFCOORD")
            parse_matblock(fd,dat.objfcoord,3)
        end

        if startswith(line,"FCOORD")
            parse_matblock(fd,dat.fcoord,4)
        end

        if startswith(line,"HCOORD")
            parse_matblock(fd,dat.hcoord,4)
        end

        if startswith(line,"DCOORD")
            parse_matblock(fd,dat.dcoord,3)
        end
    end

    return dat
end


text:
function read_input(filename::String)
    d =  Dict{Symbol,Any}()
    open(filename) do f
        # read meshfilename and outfilename options
        for i=1:2
            r = split(readline(f))
            d[Symbol(r[1])] = r[2]
        end
        # read maxsteps and printskip options
        for i=1:2
            r = split(readline(f))
            d[Symbol(r[1])] = parse(Int,r[2])
        end
        # read pause(1=true) option
        r = split(readline(f))
        d[Symbol(r[1])] = Bool(parse(Int,r[2]))

        #skip blank line
        readline(f)

        #read dt and source
        for i=1:2
            r = split(readline(f))
            d[Symbol(r[1])] = parse(Float64,r[2])
        end

        #read gradtype
        r = split(readline(f))
        d[Symbol(r[1])] = r[2]

        #skip blank line
        readline(f)

        #read condudctivity and rho*C
        for i=1:2
            r = split(readline(f))
            d[Symbol(r[1])] = parse(Float64,r[2])
        end

        #read velocity
        r = split(readline(f))
        d[Symbol(r[1])] = Vec2D(parse(Float64,r[2]), parse(Float64,r[3]))

        #read viscosity and density
        for i=1:2
            r = split(readline(f))
            d[Symbol(r[1])] = parse(Float64,r[2])
        end

        #skip blank line
        readline(f)

        #read velocity
        r = split(readline(f))
        d[Symbol(r[1])] = parse(Bool,r[2])
        #read boundary conditions
        r = split(readline(f))
        r[1] == "BCs" || error("Error reading $filename file line 15. Expected parameter \"BCs\", got \"$(r[1])\" instead.")
        nbcs = parse(UInt,r[2])
        bcs = ()
        for i=1:nbcs
            r = split(readline(f))
            bc = Symbol(r[1])=>Vec2D(parse(Float64,r[2]), parse(Float64,r[3]))
            bcs = (bcs...,bc)
        end
        d[:uBCs] = bcs
        #skip blank line
        readline(f)

        #read temperature
        r = split(readline(f))
        d[Symbol(r[1])] = parse(Bool,r[2])
        #read boundary conditions
        r = split(readline(f))
        r[1] == "BCs" || error("Error reading $filename file line 15. Expected parameter \"BCs\", got \"$(r[1])\" instead.")
        nbcs = parse(UInt,r[2])
        bcs2 = ()
        for i=1:nbcs
            r = split(readline(f))
            bc2 = Symbol(r[1])=>parse(Float64,r[2])
            bcs2 = (bcs2...,bc2)
        end
        d[:BCs] = bcs2
    end
    return d
end

read_input() = read_input("input.txt")

function show_parameters(d::Dict)
    msg = "Parameters from input:\n"
    for el in d
        msg *= string(el[1],": ",el[2],"\n")
    end
    msg*="\n"
    return msg
end

const available_laplacians = ("Method-B","Method-Gauss","Corr-Gauss","Mimetic-Cell","Implicit-B")

function get_laplacian_method(Tc,bcond,k,mesh,d)
    method = d[:gradtype]
    if method == "Method-B"
        return MethodB()
    elseif method == "Method-Gauss"
        return Gauss(Tc,mesh)
    elseif method == "Corr-Gauss"
        return CorrectedGauss(Tc,mesh)
    elseif method == "Mimetic-Cell" || method == "Mimetic_Cell"
        return CellMimetic(k,bcond,mesh)
    elseif method == "Implicit-B" || method == "Implicit_B"
        return ImplicitDiff(k,bcond,mesh,d)
    else
        warn("Gradtype \"",method,"\" not implemented. Using \"Corr-Gauss\". Available options are: ",available_laplacians)
        return CorrectedGauss(Tc,mesh)
    end
end

const available_advections = ("FullUpWind","UpWind2ndOrder")

function get_advection_method(d)
    if haskey(d,:advectype)
        method = d[:advectype]
        if method == "FullUpWind"
            return FullUpWind()
        elseif method == "UpWind2ndOrder"
            return UpWind2ndOrder()
        else
            warn("Advectype \"",method,"\" not implemented. Using \"UpWind2ndOrder\". Available options are: ",available_advections)
            return FullUpWind()
        end
    else
        return UpWind2ndOrder()
    end
end


text:
using MacroTools

macro internal_patch!(assert_existence, _expr)
    expr = prettify(_expr)
    mname = Symbol[]
    check = Expr[]
    fixed = Expr[]

    function extract_param_expr(ex::Expr)
        if  (ex.head == :.)         &&
            (length(ex.args) == 2)
            cname = ex.args[1]
            # pname = ex.args[2]  unused...

            if  (cname.head == :ref)       &&
                (length(cname.args) == 2)
                if isa(cname.args[1], Symbol)
                    return (cname.args[1], :(haskey($(cname.args[1]), $(cname.args[2]))), :($(ex).fixed = true))
                elseif isa(cname.args[1], Expr)  &&
                       (cname.args[1].head == :ref)       &&
                       (length(cname.args[1].args) == 2)  &&
                       isa(cname.args[1].args[1], Symbol)
                    return (cname.args[1].args[1], :(haskey($(cname.args[1]), $(cname.args[2]))), :($(ex).fixed = true))
                end
            end
        end
        return nothing
    end

    function recursive_search(ex)
        if isa(ex, Expr)
            if  (ex.head == :(=))       &&
                (length(ex.args) == 2)
                dd = extract_param_expr(ex.args[1])
                if !isnothing(dd)
                    push!(mname, dd[1])
                    push!(check, dd[2])
                    push!(fixed, dd[3])
                end
                for i in 2:length(ex.args)
                    recursive_search(ex.args[i])
                end
            else
                dd = extract_param_expr(ex)
                if !isnothing(dd)
                    push!(mname, dd[1])
                    push!(check, dd[2])
                else
                    for i in 1:length(ex.args)
                        recursive_search(ex.args[i])
                    end
                end
            end
        end
    end

    recursive_search(expr)
    @assert length(unique(mname)) == 1 "More than one Model or MultiModel name was found, but only one was expected: $mname"
    mname = mname[1]
    check = unique(check);  # for i in 1:length(check);  println(check[i]);  end
    fixed = unique(fixed);  # for i in 1:length(fixed);  println(fixed[i]);  end
    check = Expr(:&&, check...)
    fixed = Expr(:block, fixed...)
    if assert_existence
        out = prettify(:(
            @assert $check;
            $fixed;
            GFit.patch!($mname, GFit.@exprfunc $mname -> $expr)
        ))
    else
        out = prettify(:(
            if $check
            $fixed;
            GFit.patch!($mname, GFit.@exprfunc $mname -> $expr)
            end
        ))
    end
    # println(out)
    return esc(out)
end

macro try_patch!(_expr)
    return esc(:(GFit.@internal_patch!(false, $_expr)))
end

macro patch!(_expr)
    return esc(:(GFit.@internal_patch!(true, $_expr)))
end


text:
# This is code for the figures in the introduction of the paper

using Revise

using BridgeLandmarks
const BL=BridgeLandmarks
using RCall
using Plots
using Random
using Distributions
using DataFrames
using DelimitedFiles
using CSV
using StaticArrays
using LinearAlgebra
using StaticArrays
using LinearAlgebra

Random.seed!(3)

workdir = @__DIR__
cd(workdir)
include(joinpath(BL.dir(),"scripts", "postprocessing.jl"))
outdir = workdir
#mkpath(joinpath(outdir, "forward"))

model = [:ms, :ahs][1]
T = 1.0; dt = 0.005; t = 0.0:dt:T; tt_ =  tc(t,T)

n = 10
nshapes = 1

q0 = [PointF(2.0cos(t), sin(t))  for t in collect(0:(2pi/n):2pi)[2:end]]
fac = model==:ms ? 1.0 : 45.0
p0 = fac*[PointF(1.0, -3.0) for i in 1:n]
x0 = State(q0, p0)

ainit = mean(norm.([x0.q[i]-x0.q[i-1] for i in 2:n]))/2.0   # Let op: door 2 gedeeld
if model == :ms
    cinit = 0.2
    γinit = .2
    P = MarslandShardlow(ainit, cinit, γinit, 0.0, n)
elseif model == :ahs
    cinit = 0.02
    γinit = 0.2
    stdev = 0.75
    nfsinit = construct_nfs(2.5, stdev, γinit)
    P = Landmarks(ainit, cinit, n, 2.5, stdev, nfsinit)
end

# simulate forward
Wf, Xf = landmarksforward(tt_, x0, P)
xobsT =[ BL.q(Xf.yy[end])] # [Xf.yy[end].q[i] for i in 1:n ]

Xfsave = []# typeof(zeros(length(tt_) * P.n * 2 * d * nshapes))[]
push!(Xfsave, BL.convert_samplepath(Xf))
write_mcmc_iterates(Xfsave, tt_, n, nshapes, [1], joinpath(outdir,"forward"))

################################# start settings #################################
ITER = 200
subsamples = 0:1:ITER

fixinitmomentato0 = false
obs_atzero = true
updatescheme =  [:innov, :mala_mom]

σobs = 0.01   # noise on observations
Σobs = fill([σobs^2 * one(UncF) for i in 1:n],2)

################################# MCMC tuning pars #################################
ρinit = 0.9              # pcN-step
covθprop =   [0.04 0. 0.; 0. 0.04 0.; 0. 0. 0.04]
if model==:ms
    δinit = [0.001, 0.02] # first comp is not used
else
    δinit = [0.1, 0.1] # first comp is not used
end
η(n) = min(0.2, 10/n)  # adaptation rate for adjusting tuning pars
adaptskip = 20  # adapt mcmc tuning pars every adaptskip iters
maxnrpaths = 10 # update at most maxnrpaths Wiener increments at once
tp = MCMC_TuningPars(ρinit, maxnrpaths, δinit,covθprop,η,adaptskip)

################################# initialise P #################################
xobs0 = BL.q(x0)
################## prior specification with θ = (a, c, γ) ########################
priorθ = product_distribution([Exponential(ainit), Exponential(cinit), Exponential(γinit)])
κ = 100.0
priormom = prior_momenta = MvNormalCanon(gramkernel(xobs0,P)/κ)

#########################
xinit = State(x0.q, zeros(PointF,n))
mT = zeros(PointF,n)

start = time() # to compute elapsed time
    Xsave, parsave, objvals, accpcn, accinfo, δ, ρ, covθprop =
    lm_mcmc(tt_, (xobs0,xobsT), Σobs, mT, P,
              obs_atzero, fixinitmomentato0, ITER, subsamples,
              xinit, tp, priorθ, priormom, updatescheme,
            outdir)
elapsed = time() - start

#----------- post processing -------------------------------------------------
println("Elapsed time: ",round(elapsed/60;digits=2), " minutes")
perc_acc_pcn = mean(accpcn)*100
println("Acceptance percentage pCN step: ", round(perc_acc_pcn;digits=2))
write_mcmc_iterates(Xsave, tt_, n, nshapes, subsamples, outdir)
write_info(model,ITER, n, tt_, updatescheme, Σobs, tp, ρ, δ, perc_acc_pcn, elapsed, outdir)
write_acc(accinfo, accpcn, nshapes, outdir)

v0 = Vector(xobs0)
vT = [Vector(xobsT[1])]
write_observations(v0, vT, n, nshapes, outdir)
write_noisefields(P, outdir)


text:
# Copyright (c) 2021 Thibaut Cuvelier and contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module Chuffed

import ConstraintProgrammingExtensions
using Chuffed_jll

const CP = ConstraintProgrammingExtensions

function run_chuffed(args)
    return String(read(`$(Chuffed_jll.fznchuffed()) $(args)`))
end

function Optimizer()
    return CP.FlatZinc.Optimizer(Chuffed_jll.fznchuffed())
end

end # module


text:
# VERSION >= v"0.4" && __precompile__()

module VegaLite

    using JSON, Compat

    import Base: +, *, scale, show

    export VegaLiteVis, scale, axis, legend

    export svg, buttons

    export data_values

    export config_grid, config_facet_axis, config_facet_cell,
           config_facet_scale

    export config, config_cell, config_mark, config_scale,
           config_axis, config_legend

    export mark_bar, mark_circle, mark_square, mark_tick,
           mark_line, mark_area, mark_point, mark_text

    export encoding_x_quant, encoding_x_temp, encoding_x_ord, encoding_x_nominal,
           encoding_y_quant, encoding_y_temp, encoding_y_ord, encoding_y_nominal,
           encoding_color_quant, encoding_color_temp, encoding_color_ord, encoding_color_nominal,
           encoding_shape_quant, encoding_shape_temp, encoding_shape_ord, encoding_shape_nominal,
           encoding_path_quant, encoding_path_temp, encoding_path_ord, encoding_path_nominal,
           encoding_size_quant, encoding_size_temp, encoding_size_ord, encoding_size_nominal,
           encoding_text_quant, encoding_text_temp, encoding_text_ord, encoding_text_nominal,
           encoding_detail_quant, encoding_detail_temp, encoding_detail_ord, encoding_detail_nominal,
           encoding_order_quant, encoding_order_temp, encoding_order_ord, encoding_order_nominal,
           encoding_row_quant, encoding_row_temp, encoding_row_ord, encoding_row_nominal,
           encoding_column_quant, encoding_column_temp, encoding_column_ord, encoding_column_nominal

    include("utils.jl")
    include("render.jl")
    include("axis.jl")
    include("scale.jl")
    include("legend.jl")
    include("config.jl")
    include("data_values.jl")
    include("mark.jl")
    include("encoding.jl")

    # include("render.jl")

    ### Integration with Escher
    Pkg.installed("Escher") != nothing && include("escher_integration.jl")

    ### Integration with DataFrames
    Pkg.installed("DataFrames") != nothing && include("dataframes_integration.jl")

    ### Integration with IJulia - Jupyter
    Pkg.installed("IJulia") != nothing && include("ijulia_integration.jl")

end