text:
#const TableauHSPARKsecondary = AbstractTableauSPARK{:hspark_secondary}
const TableauHSPARKsecondary = TableauVSPARKsecondary
const ParametersHSPARKsecondary = AbstractParametersSPARK{:hspark_secondary}

@doc raw"""
Specialised Partitioned Additive Runge-Kutta integrator for Hamiltonian systems
subject to Dirac constraints with projection on secondary constraint *EXPERIMENTAL*.

This integrator solves the following system of equations for the internal stages,
```math
\begin{aligned}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\
P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\
\tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\
\tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\
0 &= \sum \limits_{j=1}^{r} \omega_{ij} \tilde{\Phi}_{n,j} , & i &= 1, ..., r-1 ,
\end{aligned}
```
with definitions
```math
\begin{aligned}
V_{n,i} &= \hphantom{-} \frac{\partial H}{\partial p} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\
F_{n,i} &=           -  \frac{\partial H}{\partial q} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\
U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
G_{n,i} &=           -  \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r ,
\end{aligned}
```
and update rule
```math
\begin{aligned}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} , \\
0 &= \phi (q_{n+1}, p_{n+1}) .
\end{aligned}
```
"""
struct IntegratorHSPARKsecondary{DT, TT, D, S, R, PT <: ParametersHSPARKsecondary{DT,TT,D,S,R},
                                                  ST <: NonlinearSolver{DT},
                                                  IT <: InitialGuessPODE{TT}} <: AbstractIntegratorHSPARK{DT,TT,D,S,R}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorHSPARKsecondary(params::ParametersHSPARKsecondary{DT,TT,D,S,R}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,R,ST,IT}
        new{DT, TT, D, S, R, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorHSPARKsecondary{DT,D}(equations::NamedTuple, tableau::TableauHSPARKsecondary{TT}, Δt::TT) where {DT,TT,D,ST}
        # get number of stages
        S = tableau.s
        R = tableau.r

        N = 2*D*S + 4*D*R

        if isdefined(tableau, :d) && length(tableau.d) > 0
            N += D
        end

        # create params
        params = ParametersHSPARKsecondary{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, N, params, caches)

        # create initial guess
        iguess = InitialGuessPODE(get_config(:ig_extrapolation), equations[:v̄], equations[:f̄], Δt)

        # create integrator
        IntegratorHSPARKsecondary(params, solver, iguess, caches)
    end

    function IntegratorHSPARKsecondary(equation::HDAE{DT}, tableau::TableauHSPARKsecondary, Δt; kwargs...) where {DT}
        IntegratorHSPARKsecondary{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

GeometricBase.nconstraints(::IntegratorHSPARKsecondary{DT,TT,D}) where {DT,TT,D} = D

function initial_guess!(int::IntegratorHSPARKsecondary{DT}, sol::AtomicSolutionPDAE{DT},
                        cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT}
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄,
                              sol.q, sol.p, sol.v, sol.f,
                              cache.q̃, cache.p̃, cache.ṽ, cache.f̃,
                              tableau(int).q.c[i], tableau(int).p.c[i])

        for k in eachdim(int)
            int.solver.x[2*(ndims(int)*(i-1)+k-1)+1] = (cache.q̃[k] - cache.q[k])/timestep(int)
            int.solver.x[2*(ndims(int)*(i-1)+k-1)+2] = (cache.p̃[k] - cache.p[k])/timestep(int)
        end
    end

    for i in 1:pstages(int)
        evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄,
                              sol.q, sol.p, sol.v, sol.f,
                              cache.q̃, cache.p̃, cache.ṽ, cache.f̃,
                              tableau(int).q̃.c[i], tableau(int).p̃.c[i])

        for k in eachdim(int)
            int.solver.x[2*ndims(int)*nstages(int)+4*(ndims(int)*(i-1)+k-1)+1] = (cache.q̃[k] - cache.q[k])/timestep(int)
            int.solver.x[2*ndims(int)*nstages(int)+4*(ndims(int)*(i-1)+k-1)+2] = (cache.p̃[k] - cache.p[k])/timestep(int)
            int.solver.x[2*ndims(int)*nstages(int)+4*(ndims(int)*(i-1)+k-1)+3] = 0
            int.solver.x[2*ndims(int)*nstages(int)+4*(ndims(int)*(i-1)+k-1)+4] = 0
        end
    end

    if isdefined(tableau(int), :λ) && tableau(int).λ.c[1] == 0
        for k in eachdim(int)
            int.solver.x[2*ndims(int)*nstages(int)+4*ndims(int)*pstages(int)+k] = 0
        end
    end
end

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,D,S,R},
                                        params::ParametersHSPARKsecondary{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local t::TT

    for i in 1:S
        for k in 1:D
            # copy x to Y, Z
            cache.Yi[i][k] = x[2*(D*(i-1)+k-1)+1]
            cache.Zi[i][k] = x[2*(D*(i-1)+k-1)+2]

            # compute Q and P
            cache.Qi[i][k] = params.q[k] + params.Δt * cache.Yi[i][k]
            cache.Pi[i][k] = params.p[k] + params.Δt * cache.Zi[i][k]
        end

        # compute f(X)
        t = params.t + params.Δt * params.tab.p.c[i]
        params.equs[:v](t, cache.Qi[i], cache.Pi[i], cache.Vi[i])
        params.equs[:f](t, cache.Qi[i], cache.Pi[i], cache.Fi[i])
    end

    for i in 1:R
        for k in 1:D
            # copy y to Y, Z and Λ
            cache.Yp[i][k] = x[2*D*S+4*(D*(i-1)+k-1)+1]
            cache.Zp[i][k] = x[2*D*S+4*(D*(i-1)+k-1)+2]
            cache.Up[i][k] = x[2*D*S+4*(D*(i-1)+k-1)+3]
            cache.Λp[i][k] = x[2*D*S+4*(D*(i-1)+k-1)+4]

            # compute Q and V
            cache.Qp[i][k] = params.q[k] + params.Δt * cache.Yp[i][k]
            cache.Pp[i][k] = params.p[k] + params.Δt * cache.Zp[i][k]
        end

        # compute f(X)
        t = params.t + params.Δt * params.tab.p̃.c[i]
        params.equs[:v](t, cache.Qp[i], cache.Pp[i], cache.Vp[i])
        params.equs[:f](t, cache.Qp[i], cache.Pp[i], cache.Fp[i])

        params.equs[:g](t, cache.Qp[i], cache.Pp[i], cache.Up[i], cache.Gp[i])
        params.equs[:g](t, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.G̅p[i])

        params.equs[:ϕ](t, cache.Qp[i], cache.Pp[i], cache.Φp[i])
        params.equs[:ψ](t, cache.Qp[i], cache.Pp[i], cache.Vp[i], cache.Fp[i], cache.Ψp[i])
    end

    # if length(params.tab.d) > 0
    #     for k in 1:D
    #         cache.μ[k] = x[2*D*S+4*D*R+k]
    #     end
    # end

    # compute q and p
    cache.q̃ .= params.q
    cache.p̃ .= params.p
    for i in 1:S
        cache.q̃ .+= params.Δt .* params.tab.q.b[1][i] .* cache.Vi[i]
        cache.p̃ .+= params.Δt .* params.tab.p.b[1][i] .* cache.Fi[i]
    end
    for i in 1:R
        cache.q̃ .+= params.Δt .* params.tab.q.b[2][i] .* cache.Up[i]
        cache.q̃ .+= params.Δt .* params.tab.q.b[3][i] .* cache.Λp[i]
        cache.p̃ .+= params.Δt .* params.tab.p.b[2][i] .* cache.Gp[i]
        cache.p̃ .+= params.Δt .* params.tab.p.b[3][i] .* cache.G̅p[i]
    end

    # compute ϕ(q,p)
    t = params.t + params.Δt
    params.equs[:ϕ](t, cache.q̃, cache.p̃, cache.ϕ̃)
end

"Compute stages of specialised partitioned additive Runge-Kutta methods for variational systems."
function Integrators.function_stages!(y::Vector{ST}, b::Vector{ST}, params::ParametersHSPARKsecondary{DT,TT,D,S,R},
                                      caches::CacheDict) where {ST,DT,TT,D,S,R}

    # get cache for internal stages
    cache = caches[ST]

    compute_stages!(y, cache, params)

    # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ]
    for i in 1:S
        for k in 1:D
            b[2*(D*(i-1)+k-1)+1] = - cache.Yi[i][k]
            b[2*(D*(i-1)+k-1)+2] = - cache.Zi[i][k]
            for j in 1:S
                b[2*(D*(i-1)+k-1)+1] += params.tab.q.a[1][i,j] * cache.Vi[j][k]
                b[2*(D*(i-1)+k-1)+2] += params.tab.p.a[1][i,j] * cache.Fi[j][k]
            end
            for j in 1:R
                b[2*(D*(i-1)+k-1)+1] += params.tab.q.a[2][i,j] * cache.Up[j][k]
                b[2*(D*(i-1)+k-1)+1] += params.tab.q.a[3][i,j] * cache.Λp[j][k]
                b[2*(D*(i-1)+k-1)+2] += params.tab.p.a[2][i,j] * cache.Gp[j][k]
                b[2*(D*(i-1)+k-1)+2] += params.tab.p.a[3][i,j] * cache.G̅p[j][k]
            end
        end
    end

    # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ, ωΨ]
    for i in 1:R
        for k in 1:D
            b[2*D*S+4*(D*(i-1)+k-1)+1] = - cache.Yp[i][k]
            b[2*D*S+4*(D*(i-1)+k-1)+2] = - cache.Zp[i][k]
            b[2*D*S+4*(D*(i-1)+k-1)+3] = - cache.Φp[i][k]
            # b[2*D*S+4*(D*(i-1)+k-1)+4] = - cache.Ψp[i][k]
            b[2*D*S+4*(D*(i-1)+k-1)+4] = 0
            for j in 1:S
                b[2*D*S+4*(D*(i-1)+k-1)+1] += params.tab.q̃.a[1][i,j] * cache.Vi[j][k]
                b[2*D*S+4*(D*(i-1)+k-1)+2] += params.tab.p̃.a[1][i,j] * cache.Fi[j][k]
            end
            for j in 1:R
                b[2*D*S+4*(D*(i-1)+k-1)+1] += params.tab.q̃.a[2][i,j] * cache.Up[j][k]
                b[2*D*S+4*(D*(i-1)+k-1)+1] += params.tab.q̃.a[3][i,j] * cache.Λp[j][k]
                b[2*D*S+4*(D*(i-1)+k-1)+2] += params.tab.p̃.a[2][i,j] * cache.Gp[j][k]
                b[2*D*S+4*(D*(i-1)+k-1)+2] += params.tab.p̃.a[3][i,j] * cache.G̅p[j][k]
            end
            for j in 1:R
                b[2*D*S+4*(D*(i-1)+k-1)+4] -= params.tab.ω[i,j] * cache.Ψp[j][k]
            end
            b[2*D*S+4*(D*(i-1)+k-1)+4] -= params.tab.ω[i,R+1] * cache.ϕ̃[k]
        end
    end

    # if length(params.tab.d) > 0
    #     for i in 1:R
    #         for k in 1:D
    #             b[2*D*S+4*(D*(i-1)+k-1)+2] -= cache.μ[k] * params.tab.d[i] / params.tab.p.b[2][i]
    #         end
    #     end
    #
    #     for k in 1:D
    #         b[2*D*S+4*D*R+k] = 0
    #         for i in 1:R
    #             b[2*D*S+4*D*R+k] -= cache.Vp[i][k] * params.tab.d[i]
    #         end
    #     end
    # end
end

function update_solution!(int::IntegratorHSPARKsecondary{DT,TT}, sol::AtomicSolutionPDAE{DT,TT},
                          cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT,TT}
    # compute final update
    update_solution!(sol.q, sol.q̃, cache.Vi, int.params.tab.q.b[1], timestep(int))
    update_solution!(sol.p, sol.p̃, cache.Fi, int.params.tab.p.b[1], timestep(int))

    # compute projection
    update_solution!(sol.q, sol.q̃, cache.Up, int.params.tab.q.b[2], timestep(int))
    update_solution!(sol.q, sol.q̃, cache.Λp, int.params.tab.q.b[3], timestep(int))
    update_solution!(sol.p, sol.p̃, cache.Gp, int.params.tab.p.b[2], timestep(int))
    update_solution!(sol.p, sol.p̃, cache.G̅p, int.params.tab.p.b[3], timestep(int))
end


text:
const R = R_GAS
const F = FARADAY

struct SimpleLinearModel{T<:Real}
    standardpotential::T
    slope::T
    rsquared::T
    SimpleLinearModel{T}(a,b,c) where {T<:Real} = new(a,b,c)
end

SimpleLinearModel(a::Float64, b::Float64, c::Float64) = SimpleLinearModel{Float64}(a,b,c)
#SimpleLinearModel(a::Real, b::Real)  = SimpleLinearModel{Real}(promote(a, b)...)
function (s::SimpleLinearModel)(x)
    return s.slope .* x .+ s.standardpotential
end

struct ModelResult{S, T}
    themodel::S
    theplot::T
end

function gen_model(n, T=723.0)
  coeff = R*T*log(10.0)/(n*F)
    model(x, E0) = E0[1] .+ coeff .* x
    return model
end

function get_rsquared(model, xdata, ydata)
    ȳ = mean(ydata)
    sstot = sum((ydata .- ȳ).^2)
    ssres = sum((ydata .- model(xdata)).^2)
    R² = 1 - ssres/sstot
    return R²
end

function fitamodel(tofit, xdata, ydata)
    p0 = [1.0]
    fit = curve_fit(tofit, xdata, ydata, p0)
    e0 = fit.param[1]
    slope = tofit(1.0, e0)  - tofit(0.0, e0)
    Rsq = get_rsquared(x->tofit(x,e0), xdata, ydata)
    return SimpleLinearModel(e0, slope, Rsq)
end

function infinitedilution_extrapolate(xdata, vdata, n, T)
    logxs = log10.(xdata)
    fitmodel = gen_model(n, T)
    initialmodel = fitamodel(fitmodel, logxs, vdata)
    Rsq = initialmodel.rsquared
    finalmodel = initialmodel
    modeltwopoints = fitamodel(fitmodel, logxs[1:2], vdata[1:2])
    oldRsq = modeltwopoints.rsquared
    # this is to initialize finalmodel with the right type; we're going
    # to figure out what the REAL final model is in a bit
    for i in 3:length(xdata)
        if Rsq > 0.99
            break # assume if Rsq is already high, model is already good
        elseif i == 3
            if oldRsq > Rsq finalmodel = modeltwopoints end
        end
        xn = logxs[1:i]
        vn = vdata[1:i]
        newmodel = fitamodel(fitmodel, xn, vn)
        Rsq = get_rsquared(newmodel, xn, vn)
        if Rsq < oldRsq
            if Rsq - oldRsq > 0.02 || Rsq < 0.95 || Rsq /oldRsq < 0.95
                break
            end
        end
        oldRsq = Rsq
        finalmodel = newmodel
    end
    return finalmodel
end

function plot_model(model, xdata, vdata, mixname="Mixture", addname="Additive")
    logxs = log10.(xdata)
    myxrange = minimum(logxs) .- 0.1:0.01:0.0
    mynewxs = collect(myxrange)
    mynewys = broadcast(model, mynewxs)
    scene = scatter(logxs, vdata, label="Experimental", legend=:topleft)
    plot!(scene, mynewxs, mynewys, label="Extrapolated")
    title!(scene, "Standard Potential Extrapolation")
    xlabel!(scene, "Log10 of Concetration of $addname in $mixname")
    ylabel!(scene, "Measured Potential (V)")
    return ModelResult(model, scene)
end

function compare_actcoeffs(extrap_model, purecomp_enot, xdata, edata,mixname="Mixture", addname="Additive")
    logxs = log10.(xdata)
    rightslope = extrap_model(1.0) - extrap_model(0.0)
    puremodel(x) = x .* rightslope .+ purecomp_enot
    unrefined_pure = edata .- puremodel(logxs)
    unrefined_infd = edata .- extrap_model(logxs)
    newcoeff = 1.0/rightslope
    pureγ = 10 .^ (unrefined_pure .* newcoeff)
    infdγ = 10 .^ (unrefined_infd .* newcoeff)
    print("Pure Component Activities for $addname in $mixname: ", pureγ)
    print("Infinite Dilution Activities for $addname in $mixname: ", infdγ)
    scene = scatter(log10.(xdata), log10.(infdγ), label="Extrapolated", legend=:topleft, size=(400,400))
    scatter!(scene, log10.(xdata), log10.(pureγ), label="Pure Component")
    #title!(scene, latexstring("\\mathrm{Activity \\: Coefficient \\: Comparison \\: for \\: $addname \top\: in \\: $mixname}"))
    xlabel!(scene, latexstring("\\mathrm{Log}_{10} \\mathrm{ \\: of \\: concentration \\: (mole \\: fraction)}"))
    ylabel!(scene, latexstring("\\mathrm{Log}_{10} \\mathrm{ \\: of \\: activity \\: coefficient}"))
    return scene
end

ext: jl


text:
using Documenter, DocumenterMarkdown, Turing, AdvancedHMC, Bijectors
using LibGit2: clone

# Include the utility functions.
include("make-utils.jl")

# Paths.
source_path = joinpath(@__DIR__, "src")
build_path = joinpath(@__DIR__, "_docs")
tutorial_path = joinpath(@__DIR__, "_tutorials")

with_clean_docs(source_path, build_path) do source, build
    makedocs(
        sitename = "Turing.jl",
        source = source,
        build = build,
        format = Markdown(),
        checkdocs = :all
    )
end

# You can skip this part if you are on a metered
# connection by calling `julia make.jl no-tutorials`
in("no-tutorials", ARGS) || copy_tutorial(tutorial_path)

version_rx = r"v\d.\d.\d"
baseurl = "/dev"
ghref = replace(get(ENV, "GITHUB_REF", ""), "/refs/tags" => "")
if get(ENV, "TRAVIS_TAG", "") != ""
    baseurl = "/" * ENV["TRAVIS_TAG"]
elseif !isnothing(match(version_rx, ghref))
    baseurl = "/" * ghref
end
jekyll_build = joinpath(@__DIR__, "jekyll-build")
with_baseurl(() -> run(`$jekyll_build`), baseurl)

# deploy
devurl = "dev"
repo = "github.com:TuringLang/Turing.jl.git"
deploydocs(
    target = "_site",
    repo = repo,
    branch = "gh-pages",
    devbranch = "master",
    devurl = devurl,
    versions = ["stable" => "v^", "v#.#", devurl => devurl]
)


text:
struct Lock <: Base.AbstractLock
    owner::typeof(Reagents.Ref{Union{Nothing,Task}}())
    notify::typeof(Reagents.channel(Nothing)[1])
    wait::typeof(Reagents.channel(Nothing)[2])
    depth::typeof(Reagents.Ref{Int}())
end

function Lock()
    notify, wait = Reagents.channel(Nothing)
    return Lock(
        Reagents.Ref{Union{Nothing,Task}}(nothing),
        notify,
        wait,
        Reagents.Ref{Int}(0),
    )
end

event(::typeof(lock), l::Lock, requester::Task = current_task()) =
    Computed() do _
        owner = l.owner[]
        @trace(
            label = :locking,
            ownerid = owner === nothing ? UInt(0) : taskid(owner),
            requesterid = taskid(requester),
        )
        if owner === requester
            # Note: Since the increment for `l.depth` can be executed on any
            # task if there is a downstream `Swap` reagent, we need to
            # simultaneously verify that `l.owner[]` is not changed.
            d = l.depth[]
            return CAS(l.owner, owner, owner) ⨟ CAS(l.depth, d, d + 1) ⨟ PostCommit() do _
                @trace(label = :locked_rec, requesterid = taskid(requester))
            end
        end
        cas = CAS(l.owner, nothing, requester)
        if owner === nothing
            return cas
        else
            return (l.wait ⨟ cas | Map(_ -> l.owner[] === nothing ? Retry() : Block()))
        end
    end

function Base.lock(l::Lock)
    apply(lock, l)
    owner = l.owner[]
    @trace(label = :locked_base, ownerid = owner === nothing ? UInt(0) : taskid(owner))
    @assert owner === current_task()
    return
end
Base.trylock(l::Lock) = Reagents.try!(event(lock, l)) !== nothing

function Base.unlock(l::Lock)
    owner = l.owner[]
    @trace(label = :unlock, ownerid = owner === nothing ? UInt(0) : taskid(owner))
    @assert owner === current_task()
    if l.depth[] > 0
        l.depth[] -= 1
        return
    end
    l.owner[] = nothing
    Julio.shield() do
        Reagents.trysync!(l.notify)
    end
    return
end


text:
using Base.Test
using MachO

@testset "library" begin
    # First, load in libfoo.dylib
    h = readmeta("./libfoo.dylib")

    # Assert some basic properties of this image
    @test h.is64
    @test !h.bswapped

    # Inspect the load commands.  We know what this file contains, and so we
    # can assert certain things about the load commands, etc...
    load_cmds = collect(LoadCmds(h))
    @test length(load_cmds) == 12
    @test load_cmds[1].cmd_id == MachO.LC_SEGMENT_64

    # Look at all the dylibs we will try to load
    dylibs = [c.cmd.name for c in load_cmds if c.cmd_id == MachO.LC_LOAD_DYLIB]
    @test length(dylibs) == 1
    @test "/usr/lib/libSystem.B.dylib" in dylibs
end

@testset "executable" begin
    h = readmeta("./fooifier")
    @test h.is64
    @test !h.bswapped

    load_cmds = collect(LoadCmds(h))
    @test length(load_cmds) == 18
    @test load_cmds[1].cmd_id == MachO.LC_SEGMENT_64

    dylibs = [c.cmd.name for c in load_cmds if c.cmd_id == MachO.LC_LOAD_DYLIB]
    @test length(dylibs) == 2
    @test "/usr/lib/libSystem.B.dylib" in dylibs
    @test "@rpath/libfoo.dylib" in dylibs

    # Inspect symbols
    symbols = collect(MachO.Symbols(h))
    symbol_names = [MachO.symname(s) for s in symbols]

    # Check that _main is a symbol and is defined
    main_idx = get(find(symbol_names .==  "_main"), 1, 0)
    @test main_idx > 0
    @test !MachO.isundef(MachO.deref(symbols[main_idx]))

    # Check that _foo is a symbol and is not defined
    foo_idx = get(find(symbol_names .==  "_foo"), 1, 0)
    @test foo_idx > 0
    @test MachO.isundef(MachO.deref(symbols[foo_idx]))
end


text:
using ArrayLisp
using Test

@testset "ArrayLisp.jl" begin
    # Write your tests here.
end


text:
export Model, @model
using MLStyle

using MacroTools: @q, striplines

import SimplePosets
using GeneralizedGenerated

struct Model{A,B,M<:GeneralizedGenerated.TypeLevel} <: AbstractModel{A,B,M,Nothing,Nothing}
    args  :: Vector{Symbol}
    vals  :: NamedTuple
    dists :: NamedTuple
    retn  :: Union{Nothing, Symbol, Expr}
end

function Model(theModule::Module, args, vals, dists, retn)
    M = to_type(theModule)
    A = NamedTuple{Tuple(args)}
    m = Model{A,Any,M}(args, vals, dists, retn)
    B = convert(Expr, m).args[end] |> to_type
    Model{A,B,M}(args, vals, dists, retn)
end

function type2model(::Type{Model{A,B,M}}) where {A,B,M}
    args = [fieldnames(A)...]
    body = from_type(B)
    Model(from_type(M), convert(Vector{Symbol},args), body)
end

function emptyModel(theModule::Module)
    M = to_type(theModule)
    A = NamedTuple{(),Tuple{}}
    B = (@q begin end) |> to_type
    Model{A,B,M}([], NamedTuple(), NamedTuple(), nothing)
end

function Base.merge(m1::Model, m2::Model)
    theModule = getmodule(m1)
    @assert theModule == getmodule(m2)
    vals = merge(m1.vals, m2.vals)
    args = setdiff(union(m1.args, m2.args), keys(vals))
    dists = merge(m1.dists, m2.dists)
    retn = maybesomething(m2.retn, m1.retn) # m2 first so it gets priority

    Model(theModule, args, vals, dists, retn)
end

Base.merge(m::Model, ::Nothing) = m

Model(theModule::Module, arg::Arg) = Model(theModule, Symbol[arg.x], NamedTuple(), NamedTuple(), nothing)
Model(theModule::Module, st::Assign) = Model(theModule, Symbol[], namedtuple(st.x)([st.rhs]), NamedTuple(), nothing)
Model(theModule::Module, st::Sample) = Model(theModule, Symbol[], NamedTuple(), namedtuple(st.x)([st.rhs]), nothing)
Model(theModule::Module, st::Return) = Model(theModule, Symbol[], NamedTuple(), NamedTuple(), st.rhs)
Model(theModule::Module, st::LineNumber) = emptyModel(theModule)

Model(theModule::Module, ::LineNumberNode) = emptyModel(theModule)

function Model(theModule::Module, expr :: Expr)
    nt = NamedTuple()
    @match expr begin
        :($k = $v)   => Model(theModule, Assign(k,v))
        :($k ~ $v)   => Model(theModule, Sample(k,v))
        :($k .~ $v)  => Model(theModule, Sample(k, :(For(identity, $v))))
        Expr(:return, x...) => Model(theModule, Return(x[1]))
        Expr(:block, body...) => foldl(merge, map(body) do line Model(theModule, line) end)
        :(@model $lnn $body) => Model(theModule, body)
        :(@model $lnn $args $body) => Model(theModule, args.args, body)

        x => begin
            @error "Bad argument to Model(::Expr)" expr=x
        end
    end
end

function Model(vs::Expr,expr::Expr)
    @assert vs.head == :tuple
    @assert expr.head == :block
    Model(Vector{Symbol}(vs.args), expr)
end

function Model{A,B,M}(args::Vector{Symbol}, expr::Expr) where {A,B,M}
    m1 = Model{A,B,M}(args, NamedTuple(), NamedTuple(), nothing)
    m2 = Model{A,B,M}(expr)
    merge(m1, m2)
end

function Model(theModule::Module, args::Vector{Symbol}, expr::Expr)
    m1 = Model(theModule, args, NamedTuple(), NamedTuple(), nothing)
    m2 = Model(theModule, expr)
    merge(m1, m2)
end

Model(m::Model) = m

Expr(m::Model,v) = convert(Expr,findStatement(m,v) )

Model(::LineNumberNode) = emptyModel

toargs(vs :: Vector{Symbol}) = Tuple(vs)
toargs(vs :: NTuple{N,Symbol} where {N}) = vs

macro model(vs::Expr,expr::Expr)
    theModule = __module__
    @assert vs.head == :tuple
    @assert expr.head == :block
    Model(theModule,Vector{Symbol}(vs.args), expr)
end

macro model(v::Symbol, expr::Expr)
    theModule = __module__
    Model(theModule,[v], expr)
end

macro model(expr :: Expr)
    theModule = __module__
    Model(theModule,Vector{Symbol}(), expr)
end

function Base.convert(::Type{Expr}, m::Model{T} where T)
    numArgs = length(m.args)
    args = if numArgs == 1
       m.args[1]
    elseif numArgs > 1
        Expr(:tuple, [x for x in m.args]...)
    end

    body = @q begin end

    for v ∈ setdiff(toposort(m), arguments(m))
        push!(body.args, Expr(m,v))
    end

    isnothing(m.retn) || push!(body.args, :(return $(m.retn)))

    q = if numArgs == 0
        @q begin
            @model $body
        end
    else
        @q begin
            @model $(args) $body
        end
    end

    striplines(q).args[1]
end

# function Base.get(m::Model, k::Symbol)
#     result = []

#     if k ∈ keys(m.vals)
#         push!(result, Assign(k,getproperty(m.vals, k)))
#     elseif k ∈ keys(m.dists)
#         if k ∈ keys(m.data)
#             push!(result, Observe(k,getproperty(m.dists, k)))
#         else
#             push!(result, Sample(k,getproperty(m.dists, k)))
#         end
#     end
#     return result
# end

# For pretty-printing in the REPL
Base.show(io::IO, m :: Model) = println(io, convert(Expr, m))

# export observe
# observe(m,v::Symbol) = merge(m, Model(Symbol[], NamedTuple(), NamedTuple(), nothing, Symbol[v]))
# observe(m,vs::Vector{Symbol}) = merge(m, Model(Symbol[], NamedTuple(), NamedTuple(), nothing, vs))

function findStatement(am::AbstractModel, x::Symbol)
    m = Model(am)
    x == :return && return Return(m.retn)
    x ∈ keys(m.vals) && return Assign(x,m.vals[x])
    x ∈ keys(m.dists) && return Sample(x,m.dists[x])
    x ∈ arguments(m) && return Arg(x)
    error("statement not found")
end

function statements(am::AbstractModel)
    m = Model(am)
    s = Statement[Soss.findStatement(m, v) for v in variables(m)]
    return s
end


text:
function varhac(dat,imax,ilag,imodel)

    ## VARHAC ESTIMATOR FOR THE COVARIANCE MATRIX (SPECTRAL DENSITY AT FREQUENCY ZERO)

    ## INPUT:

    ## dat:        input matrix where each row corresponds to a different observation
    ## imax:       maximum lag order considered
    ##             0 = no correction for serial correlation will be made
    ## ilag:       1 = all elements enter with the same lag
    ##             2 = own element can enter with a different lag
    ##             3 = only the own lag enters
    ## imodel:     1 = AIC is used
    ##             2 = BIC is used
    ##             3 = a fixed lag order equal to imax is used
    ## iprint:     1 = the chosen lag length will be printed on the screen

    ## OUTPUT:

    ## ccc:        VARHAC estimator

    #global ex, ex1, ex2, dat2, dep

    nt   = size(dat, 1)
    kdim = size(dat, 2)

    ex   = dat[imax:nt-1, :];
    dat2 = dat[imax:nt-1, :];
    ex1  = dat[imax:nt-1, :]
    ex2  = dat[imax:nt-1, :]
    dep  = dat[imax+1:nt, 1]

    ddd      = dat[imax+1:nt, :]
    minres   = ddd
    aic      = log.(sum(minres.^2, 1))
    minorder = zeros(kdim, 2)

    if imax>0
        minpar   = zeros(kdim, kdim*imax)
        ## ALL ELEMENTS ENTER WITH THE SAME LAG (ILAG = 1)
        if ilag==1
            for k = 1:kdim
                dep = ddd[:, k]
                ex = dat[imax:nt-1, :];

                for iorder = 1:imax
                    if iorder==1
                        ex = dat[imax:nt-1, :]
                    else
                        ex = [ex dat[imax+1-iorder:nt-iorder, :]]
                    end

                    if imodel!=3
                        ## Run the VAR
                        b = (dep\ex)'
                        resid = dep-ex[:,:]*b
                        ## Compute the model selection criterion
                        npar = kdim*iorder
                        if imodel==1
                            aicnew = log.(sum(resid.^2, 1))+2*npar/(nt-imax)
                        else
                            aicnew = log.(sum(resid.^2, 1))+log(nt-imax)*npar/(nt-imax)
                        end
                        if aicnew[1]<aic[k]
                            aic[k] = aicnew[1]
                            minpar[k, 1:kdim*iorder] = b'
                            minres[:, k] = resid
                            minorder[k] = iorder
                        end
                    end

                    if imodel==3
                        b = (dep\ex)'
                        resid = dep-ex*b
                        minpar[k, :] = b'
                        minres[: ,k] = resid
                        minorder[k]  = imax
                    end
                end
            end
                ## ONLY THE OWN LAG ENTERS (ILAG = 3
        elseif ilag==3

            for k = 1:kdim
                dep = ddd[:, k]

                for iorder = 1:imax
                    if iorder == 1
                        ex = dat[imax:nt-1, k]
                    else
                        ex = [ex dat[imax+1-iorder:nt-iorder,k]]
                    end

                    if imodel!=3;
                        ## Run the VAR
                        b = ex\dep
                        resid = dep-ex[:,:]*b
                        ## Compute the model selection criterion
                        npar = iorder
                        if imodel==1
                            aicnew = log.(sum(resid.^2, 1))+2*npar/(nt-imax);
                        else
                            aicnew = log.(sum(resid.^2, 1))+log(nt-imax)*npar/(nt-imax);
                        end

                        if aicnew[1] < aic[k]
                            aic[k] = aicnew[1]
                            for i = 1:iorder
                                minpar[k,k+(i-1)*kdim] = b[i]
                            end
                            minres[:, k] = resid
                            minorder[k]  = iorder
                        end
                    end
                end

                if imodel==3
                    b = (dep\ex)'
                    resid = dep-ex[:,:]*b
                    for i = 1:imax
                        minpar[k, k+(i-1)*kdim] = b[i]
                    end
                    minres[:, k] = resid
                    minorder[k]  = imax
                end
            end
            ## OWN ELEMENT CAN ENTER WITH A DIFFERENT LAG (ILAG = 2)
        else
            begin
            for k = 1:kdim
                dep = ddd[:, k]
                if k==1
                    dat2 = dat[:, 2:kdim];
                elseif k==kdim
                    dat2 = dat[:, 1:kdim-1]
                else
                    dat2 = [dat[:, 1:k-1] dat[:, k+1:kdim]]
                end

                for iorder = 0:imax
                    if iorder==1
                        ex1 = dat[imax:nt-1,k]
                    elseif iorder>1
                        ex1 = [ex1 dat[imax+1-iorder:nt-iorder,k]]
                    end

                    for iorder2 = 0:imax
                        if iorder2 == 1
                            ex2 = dat2[imax:nt-1, :]
                        elseif iorder2 > 1
                            ex2 = [ex2 dat2[imax+1-iorder2:nt-iorder2, :]]
                        end

                        if iorder+iorder2==0
                            break
                        elseif iorder==0
                            ex = ex2
                        elseif iorder2==0
                            ex = ex1
                        else
                            ex = [ex1 ex2]
                        end
                        ## Run the VAR
                        b = ex\dep
                        #println(size(b))
                        #println(size(ex))
                        resid = dep.-ex[:,:]*b
                        ## Compute the model selection criterion

                        npar = iorder+iorder2*(kdim-1)

                        if imodel==1
                            aicnew = log.(sum(resid.^2, 1))+2*npar/(nt-imax);
                        else
                            aicnew = log.(sum(resid.^2, 1))+log(nt-imax)*npar/(nt-imax);
                        end
                        #println(aic[k])
                        if aicnew[1]<aic[k]
                            aic[k] = aicnew[1]
                            minpar[k, :] = zeros(1,kdim*imax)
                            for i = 1:iorder
                                minpar[k,k+(i-1)*kdim] = b[i]
                            end

                            for i = 1:iorder2
                                if k==1
                                    minpar[k,2+(i-1)*kdim:kdim+(i-1)*kdim] =
                                        b[iorder+1+(i-1)*(kdim-1):iorder+i*(kdim-1)]';
                                elseif k==kdim
                                    minpar[k,1+(i-1)*kdim:kdim-1+(i-1)*kdim] =
                                        b[iorder+1+(i-1)*(kdim-1):iorder+i*(kdim-1)]';
                                else
                                    minpar[k,1+(i-1)*kdim:k-1+(i-1)*kdim] =
                                        b[iorder+1+(i-1)*(kdim-1):iorder+k-1+(i-1)*(kdim-1)]';
                                    minpar[k,k+1+(i-1)*kdim:kdim+(i-1)*kdim] =
                                        b[iorder+k+(i-1)*(kdim-1):iorder+i*(kdim-1)]';
                                end
                            end
                            minres[:, k]   = resid
                            minorder[k, :] = [iorder iorder2]
                        end
                    end
                end
            end
            end
        end
    end

    ## COMPUTE THE VARHAC ESTIMATOR
    covar = minres'minres/(nt-imax)

    bbb = eye(kdim)

    if imax>0
        for iorder = 1:imax
            bbb = bbb-minpar[:,(iorder-1)*kdim+1:iorder*kdim]
        end
    end

    inv(bbb)*covar*inv(bbb')

end


text:
abstract type Layoutable end

macro Layoutable(name::Symbol, fields::Expr = Expr(:block))

    if !(fields.head == :block)
        error("Fields need to be within a begin end block")
    end

    structdef = quote
        mutable struct $name <: Layoutable
            parent::Union{Figure, Scene, Nothing}
            layoutobservables::LayoutObservables
            attributes::Attributes
            elements::Dict{Symbol, Any}
        end
    end

    # append user defined fields to struct definition
    # linenumbernode block, struct block, fields block

    allfields = structdef.args[2].args[3].args

    append!(allfields, fields.args)

    fieldnames = map(filter(x -> !(x isa LineNumberNode), allfields)) do field
        if field isa Symbol
            return field
        end
        if field isa Expr && field.head == Symbol("::")
            return field.args[1]
        end
        error("Unexpected field format. Neither Symbol nor x::T")
    end

    # define an inner constructor that calls _on_create so that there can be some default
    # initialization routines for all layoutables right after creating the actual objects
    constructor = quote
        function $name($(fieldnames...))
            layoutable = new($(fieldnames...))
            _on_create(layoutable, $(fieldnames...))
            layoutable
        end
    end

    push!(allfields, constructor)

    structdef
end

"""
Get the scene which layoutables need from their parent to plot stuff into
"""
get_topscene(f::Figure) = f.scene
function get_topscene(s::Scene)
    if !(s.camera_controls[] isa AbstractPlotting.PixelCamera)
        error("Can only use scenes with PixelCamera as topscene")
    end
    s
end

function register_in_figure!(fig::Figure, @nospecialize layoutable::Layoutable)
    if layoutable.parent !== fig
        error("Can't register a layoutable with a different parent in a figure.")
    end
    if !(layoutable in fig.content)
        push!(fig.content, layoutable)
    end
    nothing
end

function _on_create(layoutable, parent, layoutobservables, attributes, elements, args...)
    if parent isa Figure
        register_in_figure!(parent, layoutable)
    end
    nothing
end

# almost like in AbstractPlotting
# make fields type inferrable
# just access attributes directly instead of via indexing detour

@generated Base.hasfield(x::T, ::Val{key}) where {T<:Layoutable, key} = :($(key in fieldnames(T)))

@inline function Base.getproperty(x::T, key::Symbol) where T <: Layoutable
    if hasfield(x, Val(key))
        getfield(x, key)
    else
        x.attributes[key]
    end
end

@inline function Base.setproperty!(x::T, key::Symbol, value) where T <: Layoutable
    if hasfield(x, Val(key))
        setfield!(x, key, value)
    else
        x.attributes[key][] = value
    end
end

# propertynames should list fields and attributes
function Base.propertynames(layoutable::T) where T <: Layoutable
    [fieldnames(T)..., keys(layoutable.attributes)...]
end

# treat all layoutables as scalars when broadcasting
Base.Broadcast.broadcastable(l::Layoutable) = Ref(l)

function Base.show(io::IO, ::T) where T <: Layoutable
    print(io, "$T()")
end

function Base.delete!(layoutable::Layoutable)
    for (key, d) in layoutable.elements
        try
            remove_element(d)
        catch e
            @info "Failed to remove element $key of $(typeof(layoutable))."
            rethrow(e)
        end
    end

    if hasfield(typeof(layoutable), :scene)
        delete_scene!(layoutable.scene)
    end

    GridLayoutBase.remove_from_gridlayout!(GridLayoutBase.gridcontent(layoutable))

    on_delete(layoutable)
    delete_from_parent!(layoutable.parent, layoutable)
    layoutable.parent = nothing

    nothing
end

# do nothing for scene and nothing
function delete_from_parent!(parent, layoutable::Layoutable)
end

function delete_from_parent!(figure::Figure, layoutable::Layoutable)
    filter!(x -> x !== layoutable, figure.content)
    if current_axis(figure) === layoutable
        current_axis!(figure, nothing)
    end
    nothing
end

"""
Overload to execute cleanup actions for specific layoutables that go beyond
deleting elements and removing from gridlayout
"""
function on_delete(layoutable)
end

function remove_element(x)
    delete!(x)
end

function remove_element(x::AbstractPlot)
    delete!(x.parent, x)
end

function remove_element(xs::AbstractArray)
    foreach(remove_element, xs)
end

function remove_element(::Nothing)
end

function delete_scene!(s::Scene)
    for p in copy(s.plots)
        delete!(s, p)
    end
    deleteat!(s.parent.children, findfirst(x -> x === s, s.parent.children))
    nothing
end


text:
config_file = isempty(ARGS) ? joinpath(@__DIR__, "config", "config_precompute_parameters.jl") : ARGS[1]
isfile(config_file) || error("Cannot read '$config_file'")
println("Config supplied: '$config_file'")
include(config_file)

using SVDD, OneClassActiveLearning, Distributed, Gurobi, JuMP, MLKernels
using Serialization

include(joinpath(@__DIR__, "util", "setup_workers.jl"))

@info "Loading packages on all workers."
@everywhere using SVDD, OneClassActiveLearning, OneClassSampling, Gurobi, JuMP, Memento

@everywhere function Memento.debug(logger::Logger, error::ErrorException)
    Memento.warn(logger, "Caught ErrorException, msg='$(error.msg)')")
end

@info "Searching for data sets in '$(data_input_root)'"
data_files = Vector{String}()
for data_set in vcat(values(data_dirs)...)
    data_dir = joinpath(data_input_root, data_set)
    if !isdir(data_dir)
        @info "Could not find dataset $data_dir"
        continue
    end
    local_names = readdir(data_dir)
    full_names = realpath.(joinpath.(data_input_root, data_dir, local_names))
    push!(data_files, full_names...)
end

@info "Initializing models."
all_model_params = @distributed vcat for file in data_files
    data, labels = load_data(file)
    n_obs = length(labels)
    params = Dict{Symbol, Any}(
        :file => file,
        :n_observations => n_obs
    )
    # SVDD
    model = instantiate(VanillaSVDD, data, fill(:U, n_obs), Dict{Symbol, Any}())
    initialize!(model, init_strategy[:svdd])
    model_parameters = get_model_params(model)
    params[:svdd] = Dict{Symbol, Any}(
        :type => Symbol(typeof(model)),
        :C => get_model_params(model)[:C],
        :gamma => MLKernels.getvalue(model.kernel_fct.alpha)
    )
    # PWC
    pwc_gamma = init_strategy[:pwc][:gamma](data)
    pwc_threshold = calculate_threshold(init_strategy[:pwc][:threshold_strat], data, labels, pwc_gamma)
    params[:pwc] = Dict{Symbol, Any}(
        :gamma => pwc_gamma,
        :threshold => pwc_threshold,
        :threshold_strat => Symbol(typeof(init_strategy[:pwc][:threshold_strat]))
    )
    params
end

@info "Saving parameters."
model_params_by_filename = Dict{String, Dict{Symbol, Any}}()
for model_params in all_model_params
    model_params_by_filename[basename(model_params[:file])] = model_params
    delete!(model_params, :file)
end

serialize(data_parameter_file, model_params_by_filename)
@info "Done."

max_stars_repo_name: englhardt/ocs-evaluation

max_issues_repo_name: englhardt/ocs-evaluation

max_forks_repo_name: englhardt/ocs-evaluation


text:
@defcomp SulphateForcing begin
    region = Index()

    se0_sulphateemissionsbase = Parameter(index=[region], unit="TgS/year")
    pse_sulphatevsbase = Parameter(index=[time, region], unit="%")
    se_sulphateemissions = Variable(index=[time, region], unit="TgS/year")
    area = Parameter(index=[region], unit ="km^2")

    sfx_sulphateflux = Variable(index=[time, region], unit="TgS/km^2/yr")

    d_sulphateforcingbase = Parameter(unit="W/m2", default=-0.46666666666666673)
    ind_slopeSEforcing_indirect = Parameter(unit="W/m2", default=-0.4000000000000001)
    nf_naturalsfx = Parameter(index=[region], unit="TgS/km^2/yr")

    fs_sulphateforcing = Variable(index=[time, region], unit="W/m2")

    function run_timestep(p, v, d, tt)
        bigSFX0 = p.se0_sulphateemissionsbase ./ p.area

        for rr in d.region
            v.se_sulphateemissions[tt, rr] = p.se0_sulphateemissionsbase[rr] * p.pse_sulphatevsbase[tt, rr] / 100

            # Eq.17 from Hope (2006) - sulfate flux
            v.sfx_sulphateflux[tt,rr] = v.se_sulphateemissions[tt,rr] / p.area[rr]
            # Update for Eq. 18 from Hope (2009) - sulfate radiative forcing effect
            bigSFD0 = p.d_sulphateforcingbase * bigSFX0[rr] / (sum(bigSFX0 .* p.area) / sum(p.area))
            fsd_term = bigSFD0 * v.sfx_sulphateflux[tt,rr] / bigSFX0[rr]
            fsi_term = p.ind_slopeSEforcing_indirect/log(2) * log((p.nf_naturalsfx[rr] + v.sfx_sulphateflux[tt, rr]) / p.nf_naturalsfx[rr])

            v.fs_sulphateforcing[tt, rr] = fsd_term + fsi_term
        end
    end
end


text:
function calculateToeplitzKernelBasis(img_shape_os, trj::Vector{Matrix{T}}, U::Matrix{Complex{T}}; verbose = false) where {T}

    Ncoeff = size(U, 2)

    FFTW.set_num_threads(Threads.nthreads())
    fftplan = plan_fft(Array{Complex{T}}(undef, img_shape_os); flags = FFTW.MEASURE)
    nfftplan = plan_nfft(trj[1], img_shape_os; precompute = LINEAR, blocking = false, fftflags = FFTW.ESTIMATE)

    λ = Array{Complex{T}}(undef, img_shape_os)
    Λ = Array{Complex{T}}(undef, Ncoeff, Ncoeff, prod(img_shape_os))
    Λ .= 0

    for i ∈ eachindex(trj)
        t_kernel = @elapsed calculateToeplitzKernel!(λ, nfftplan, trj[i], fftplan)

        @views U2 = conj.(U[i, :]) * transpose(U[i, :])
        t_multiplcation = @elapsed begin
            Threads.@threads for j ∈ eachindex(λ)
                @simd for iu ∈ CartesianIndices(U2)
                    @inbounds Λ[iu, j] += U2[iu] * λ[j]
                end
            end
        end

        if verbose
            println("Time frame $i: t_kernel = $t_kernel; t_multiplcation = $t_multiplcation")
            flush(stdout)
        end
    end

    return Λ
end

############################################################################################
# NFFTNormalOpBasisFunc
############################################################################################
struct NFFTNormalOpBasisFunc{S,T,E,F,G}
    shape::S
    Ncoeff::Int
    fftplan::E
    ifftplan::F
    Λ::Array{Complex{T},3}
    kL1::Array{Complex{T}}
    kL2::Array{Complex{T}}
    cmaps::G
end

function NFFTNormalOpBasisFunc(
    img_shape,
    trj::Vector{Matrix{T}},
    U::Matrix{Complex{T}};
    cmaps = (1,),
    verbose = false,
    Λ = calculateToeplitzKernelBasis(2 .* img_shape, trj, U; verbose = verbose),
) where {T}

    FFTW.set_num_threads(1)
    fftplan  = plan_fft!(Array{Complex{T}}(undef, 2 .* img_shape); flags = FFTW.MEASURE)
    ifftplan = plan_ifft!(Array{Complex{T}}(undef, 2 .* img_shape); flags = FFTW.MEASURE)
    Ncoeff = size(U, 2)
    kL1 = Array{Complex{T}}(undef, (2 .* img_shape)..., Ncoeff)
    kL2 = similar(kL1)

    return NFFTNormalOpBasisFunc(img_shape, Ncoeff, fftplan, ifftplan, Λ, kL1, kL2, cmaps)
end

function LinearAlgebra.mul!(x::Vector{T}, S::NFFTNormalOpBasisFunc, b, α, β) where {T}
    idx = CartesianIndices(S.shape)
    idxos = CartesianIndices(2 .* S.shape)
    Ncoils = length(S.cmaps)

    b = reshape(b, S.shape..., S.Ncoeff)
    if β == 0
        fill!(x, zero(T)) # to avoid 0 * NaN == NaN
    else
        x .*= β
    end

    bthreads = BLAS.get_num_threads()
    try
        BLAS.set_num_threads(1)
        for icoil ∈ 1:Ncoils
            Threads.@threads for i = 1:S.Ncoeff
                S.kL1[idxos, i] .= 0
                @views S.kL1[idx, i] .= S.cmaps[icoil] .* b[idx, i]
                @views S.fftplan * S.kL1[idxos, i]
            end

            kL1_rs = reshape(S.kL1, :, S.Ncoeff)
            kL2_rs = reshape(S.kL2, :, S.Ncoeff)
            Threads.@threads for i ∈ eachindex(view(S.Λ, 1, 1, :))
                @views @inbounds mul!(kL2_rs[i, :], S.Λ[:, :, i], kL1_rs[i, :])
            end

            Threads.@threads for i = 1:S.Ncoeff
                @views S.ifftplan * S.kL2[idxos, i]
            end

            @views x .+= α .* vec(conj.(S.cmaps[icoil]) .* S.kL2[idx, :])
        end
    finally
        BLAS.set_num_threads(bthreads)
    end
    return x
end

Base.:*(S::NFFTNormalOpBasisFunc, b::AbstractVector) = mul!(similar(b), S, b)
Base.size(S::NFFTNormalOpBasisFunc) = S.shape
Base.size(S::NFFTNormalOpBasisFunc, dim) = S.shape[dim]
Base.eltype(::Type{NFFTNormalOpBasisFunc{S,D,T}}) where {S,D,T} = T

############################################################################################
# LinearOperator of NFFTNormalOpBasisFunc
############################################################################################
function NFFTNormalOpBasisFuncLO(A::NFFTNormalOpBasisFunc{S,T,E,F,G}) where {S,T,E,F,G}
    return LinearOperator(
        Complex{T},
        prod(A.shape) * A.Ncoeff,
        prod(A.shape) * A.Ncoeff,
        true,
        true,
        (res, x, α, β) -> mul!(res, A, x, α, β),
        nothing,
        (res, x, α, β) -> mul!(res, A, x, α, β),
    )
end

function NFFTNormalOpBasisFuncLO(
    img_shape,
    trj::Vector{Matrix{T}},
    U::Matrix{Complex{T}};
    cmaps = (1,),
    verbose = false,
    Λ = calculateToeplitzKernelBasis(2 .* img_shape, trj, U; verbose = verbose),
    ) where {T}

    S = NFFTNormalOpBasisFunc(img_shape, trj, U; cmaps = cmaps, Λ = Λ)
    return NFFTNormalOpBasisFuncLO(S)
end


text:
function DiffusionDefinition.diffusion_message(::Val{:Sine})
    println("* * * * * * * ")
    println("Sine diffusion")
    println("* * * * * * * ")
    println("Description")
    println("------------")
    println("A one dimensional diffusion Xₜ:")
    println("dXₜ=(a + b×sin(c Xₜ))dt + σdWₜ.")
    println("Xₜ ∈ ℝ.")
    println("The model has four parameters:")
    println("\ta, b, c, σ")
    println()
    println("Constructors")
    println("------------")
    println("You may initialize the process either by passing the parameters")
    println("in the order above, as in:")
    println("a, b, c, σ = ...")
    println("P = Sine(a, b, c, σ)")
    println("or by keyword arguments")
    println("p = (a = ..., b = ..., c = ..., σ = ...)")
    println("P = Sine(;p...)")
    println()
end

@diffusion_process Sine{T} begin
    :parameters
    (a, b, c, σ) --> T

    :additional
    constdiff --> true
end

DiffusionDefinition.b(t, x, P::Sine) = P.a + P.b*sin(P.c * x)
DiffusionDefinition.σ(t, x, P::Sine) = P.σ

DiffusionDefinition.default_type(::Sine) = Float64
DiffusionDefinition.default_wiener_type(::Sine) = Float64


text:
include("../src/components/co2cycle_component.jl")

@testset "co2cycle" begin

    Precision = 1.0e-11
    T = length(MimiDICE2010.model_years)
    f = readxlsx(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))

    m = Model()

    set_dimension!(m, :time, MimiDICE2010.model_years)

    add_comp!(m, co2cycle, :co2cycle)

# Set the parameters that would normally be internal connection from their Excel values
    set_param!(m, :co2cycle, :E, read_params(f, "B109:BI109", T))

# Load the rest of the external parameters
    p = dice2010_excel_parameters(joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx"))
    set_param!(m, :co2cycle, :mat0, p[:mat0])
    set_param!(m, :co2cycle, :mat1, p[:mat1])
    set_param!(m, :co2cycle, :mu0, p[:mu0])
    set_param!(m, :co2cycle, :ml0, p[:ml0])
    set_param!(m, :co2cycle, :b12, p[:b12])
    set_param!(m, :co2cycle, :b23, p[:b23])
    set_param!(m, :co2cycle, :b11, p[:b11])
    set_param!(m, :co2cycle, :b21, p[:b21])
    set_param!(m, :co2cycle, :b22, p[:b22])
    set_param!(m, :co2cycle, :b32, p[:b32])
    set_param!(m, :co2cycle, :b33, p[:b33])

# Run the one-component model
    run(m)

# Extract the generated variables
    MAT     = m[:co2cycle, :MAT]
    MAT_final   = m[:co2cycle, :MAT_final]
    ML      = m[:co2cycle, :ML]
    MU      = m[:co2cycle, :MU]

# Extract the true values
    True_MAT    = read_params(f, "B112:BI112", T)
    True_MAT_final  = read_params(f, "BJ112")
    True_ML     = read_params(f, "B115:BI115", T)
    True_MU     = read_params(f, "B114:BI114", T)

# Test that the values are the same
    @test maximum(abs, MAT .- True_MAT) ≈ 0. atol = Precision
    @test abs(MAT_final - True_MAT_final) ≈ 0. atol = Precision
    @test maximum(abs, ML .- True_ML) ≈ 0. atol = Precision
    @test maximum(abs, MU .- True_MU) ≈ 0. atol = Precision

end


text:
"""
    ArcStandard()

Transition system for for Arc-Standard dependency parsing.

# Transitions

| Transition  | Definition                                         |
|:----------- |:-------------------------------------------------- |
| LeftArc(l)  | (σ\\|s1\\|s0, β, A) → (σ\\|s0, β, A ∪ (s0, l, s1)) |
| RightArc(l) | (σ\\|s, b\\|β, A) → (σ, b\\|β, A ∪ (b, l, s))      |
| Shift	      | (σ,  b\\|β, A) → (σ\\|b, β, A)                     |

# Preconditions

| Transition  | Condition                          |
|:----------- |:---------------------------------- |
| LeftArc(l)  | ¬[s1 = 0], ¬∃k∃l'[(k, l', s1) ϵ A] |
| RightArc(l) | ¬∃k∃l'[(k, l', s0) ϵ A]            |

See [Nivre 2004](https://www.aclweb.org/anthology/W04-0308.pdf).
"""
struct ArcStandard <: AbstractTransitionSystem end

initconfig(s::ArcStandard, graph::DependencyTree) = ArcStandardConfig(graph)
initconfig(s::ArcStandard, words) = ArcStandardConfig(words)

projective_only(::ArcStandard) = true

transition_space(::ArcStandard, labels=[]) =
    isempty(labels) ? [LeftArc(), RightArc(), Shift()] :
    [LeftArc.(labels)..., RightArc.(labels)..., Shift()]

struct ArcStandardConfig <: AbstractParserConfiguration
    stack::Vector{Int}
    buffer::Vector{Int}
    A::Vector{Token}
end

ArcStandardConfig(sentence) = stack_buffer_config(ArcStandardConfig, sentence)

buffer(cfg::ArcStandardConfig) = cfg.buffer
stack(cfg::ArcStandardConfig)  = cfg.stack
tokens(cfg::ArcStandardConfig) = cfg.A

function Base.show(io::IO, cfg::ArcStandardConfig)
    print(io, "$(typeof(cfg))($(stack(cfg)),$(buffer(cfg)))")
    for (i, t) in enumerate(tokens(cfg))
        print("(",join((i, t.form, t.head), ","), ")")
    end
end

function apply_transition(f, cfg::ArcStandardConfig, a...; k...)
    σ, β, A = f(cfg.stack, cfg.buffer, cfg.A, a...; k...)
    return ArcStandardConfig(σ, β, A)
end

leftarc(cfg::ArcStandardConfig, args...; kwargs...) =
    apply_transition(leftarc_reduce2, cfg, args...; kwargs...)

rightarc(cfg::ArcStandardConfig, args...; kwargs...) =
    apply_transition(rightarc_reduce, cfg, args...; kwargs...)

shift(cfg::ArcStandardConfig) =
    apply_transition(shift, cfg)

isfinal(cfg::ArcStandardConfig) =
    length(stack(cfg)) == 1 && stack(cfg)[1] == 0 && length(buffer(cfg)) == 0

"""
    static_oracle(cfg, gold_tree)

Static oracle function for arc-standard dependency parsing.
"""
function static_oracle(cfg::ArcStandardConfig, gold_tree, arc=untyped)
    l = i -> arc(gold_tree[i])
    s = stack(cfg)
    if length(s) >= 2
        s1, s0 = s[end-1], s[end]
        if has_arc(gold_tree, s0, s1)
            return LeftArc(l(s1)...)
        elseif has_arc(gold_tree, s1, s0)
            if !any(k -> (k in s || k in buffer(cfg)), deps(gold_tree, s0))
                return RightArc(l(s0)...)
            end
        end
    end
    return Shift()
end

function is_possible(::LeftArc, cfg::ArcStandardConfig)
    if stacklength(cfg) >= 2
        s, s1, s0 = popstack(cfg, 2)
        return s1 != 0 && !has_head(token(cfg, s1))
    else
        return false
    end
end

is_possible(::RightArc, cfg::ArcStandardConfig) =
    stacklength(cfg) > 1 && !has_head(token(cfg, last(stack(cfg))))

is_possible(::Shift, cfg::ArcStandardConfig) = stacklength(cfg) > 0

function possible_transitions(cfg::ArcStandardConfig, arc=untyped)
    l = i -> arc(token(cfg, i))
    if stacklength(cfg) >= 2
        σ, s1, s0 = popstack(cfg, 2)
        transitions = [LeftArc(l(s1)...), RightArc(l(s0)...), Shift()]
        return filter(t -> is_possible(t, cfg), transitions)
    else
        return filter(t -> is_possible(t, cfg), TransitionOperator[Shift()])
    end
end

possible_transitions(cfg::ArcStandardConfig, gold_tree::DependencyTree, arc=untyped) =
    possible_transitions(cfg, arc)

==(cfg1::ArcStandardConfig, cfg2::ArcStandardConfig) =
    cfg1.stack == cfg2.stack && cfg1.buffer == cfg2.buffer && cfg1.A == cfg2.A


text:
# AUTO GENERATED FILE - DO NOT EDIT

export mediaquery

"""
    mediaquery(;kwargs...)
    mediaquery(children::Any;kwargs...)
    mediaquery(children_maker::Function;kwargs...)

A MediaQuery component.
Apply styles to children if media query matches. For more information, see: https://mantine.dev/core/media-query/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Child that should be shown at given breakpoint, it must accept className prop
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `largerThan` (optional): Larger than given breakpoint or value in px, will be hidden on smaller viewport
- `query` (String; optional): Any other media query
- `smallerThan` (optional): Smaller than given breakpoint or value in px, will be hidden on larger viewport
"""
function mediaquery(; kwargs...)
        available_props = Symbol[:children, :id, :className, :largerThan, :query, :smallerThan]
        wild_props = Symbol[]
        return Component("mediaquery", "MediaQuery", "dash_mantine_components", available_props, wild_props; kwargs...)
end

mediaquery(children::Any; kwargs...) = mediaquery(;kwargs..., children = children)
mediaquery(children_maker::Function; kwargs...) = mediaquery(children_maker(); kwargs...)


text:
#holds the time-indexed parameters and a pre-allocation for the row-indexed versions
#note- we ahve to be careful about using methods from AbstractVolumeParts since we ahve no V₀

abstract type AbstractVolumePartsMCMC{TM, TV, elTΠ} <: AbstractVolumeParts{TM, TV, elTΠ} end

struct VolumePartsMCMC{
    TΨ<:AbstractVector,
    Tcc<:CapacityChains,
    TXv<:AbstractXYIter,
    TA<:AbstractMatrix,
    #Tv<:AbstractVector,
    Tts<:AbstractVector,
    Texpand<:Dict,
    Txsection<:Dict,
    TM<:AbstractMatrix,
    TV<:AbstractVector,
    elTΠ<:Real} <: AbstractVolumePartsMCMC{TM, TV, elTΠ}

  Ψ::TΨ
  cc::Tcc
  Xv::TXv

  Ã::TA

  ts::Tts
  tsM1::Tts
  expand::Texpand
  xsection::Txsection

  #inner constructor needed to identify elTΠ
  function VolumePartsMCMC(
    Ψ::TΨ, cc::Tcc, Xv::TXv, Ã::TA,
    ts::Tts, tsM1::Tts, expand::Texpand, xsection::Txsection,
    ::Type{TM}, ::Type{TV}, ::Type{T}) where {
      TΨ, Tcc, TXv, TA, Tts, Texpand, Txsection, TM, TV, T}

    (eltype(TΨ) == T) || error("Eltype of Ψ must match T in TM,TV,T")
      #replace the below
    #@assert (Π == reduce(vcat, [vec(A), vec(Ã), vec(G)]))

    return new{TΨ, Tcc, TXv, TA, Tts, Texpand, Txsection,TM, TV, T}(
      Ψ, cc, Xv, Ã, ts, tsM1, expand, xsection)
  end
end

#pre-allocates for the vpexpand
function VolumePartsMCMC(Ψ::CCState, cc::CapacityChains,  Xv::AbstractXYIter,
  ts::Tts, ::Type{TM}, ::Type{TV}, ::Type{T};
  tsM1::Tts=ts .- 1) where {TM, TV, T, Tts<:AbstractVector}

  dims=(K=size(Ψ.G,1), T=size(Ψ.G,2)+1, N=size(Xv.ws,1))

  Ã = rand(T, dims.K, dims.T)

  expand = Dict{Union{Symbol, Tuple{Symbol, Int}},AbstractArray{T}}()
  xsection = Dict{Union{Tuple{Symbol, Int}},AbstractArray{T}}()

  for t ∈ 1:(dims.T-1)

    #these expansions include all values equal or greaer than t
    tpart = ts[ts .≥ t+1]
    tpartM1 = tsM1[tsM1 .≥ t]
    expand[:ws,t] = view(Xv.ws, ts .≥ t+1, :)
    expand[:RLws,t] = view(Xv.RLws, ts .≥ t+1, :)
    expand[:v,t] = view(Xv.v, ts .≥ t+1)

    expand[:G,t] = view(Ψ.G, :, tpartM1)
    expand[:σ²,t] = view(Ψ.σ², tpartM1)
    expand[:ω,t] = view(Ψ.ω, tpartM1)

    expand[:Ã,t] = view(Ã, :, tpart)
    expand[:LÃ,t] = view(Ã, :, tpartM1)

    #these expansions are simple xsections for a particular t
    txsection = ts[ts .== t+1]
    txsectionM1 = tsM1[tsM1 .== t]
    xsection[:ws,t] = view(Xv.ws, ts .== t+1, :)
    xsection[:RLws,t] = view(Xv.RLws, ts .== t+1, :)
    xsection[:v,t] = view(Xv.v, ts .== t+1)
    xsection[:X̃,t] = TM(undef, sum(ts .== t+1), dims.K) #cross-sectional allocations
    xsection[:ṽ,t] = TV(undef, sum(ts .== t+1))

    xsection[:G,t] = view(Ψ.G, :, txsectionM1)
    xsection[:σ²,t] = view(Ψ.σ², txsectionM1)
    xsection[:ω,t] = view(Ψ.ω, txsectionM1)
    xsection[:Ã,t] = view(Ã, :, txsection)
    xsection[:LÃ,t] = view(Ã, :, txsectionM1)
  end

  #convenience alias for the complete case
  expand[:ws] = expand[:ws,1]
  expand[:RLws] = expand[:RLws,1]
  expand[:v] = expand[:v,1]
  #expand[:X̃] = expand[:X̃,1]
  #expand[:ṽ] = expand[:ṽ,1]

  expand[:G] = expand[:G,1]
  expand[:σ²] = expand[:σ²,1]
  expand[:ω] = expand[:ω,1]
  expand[:Ã] = expand[:Ã,1]
  expand[:LÃ] = expand[:LÃ,1]

  Θ = VolumePartsMCMC(Ψ, cc, Xv, Ã, ts, tsM1, expand, xsection, TM, TV, T)
  return Θ
end

#basic constructor from dimension arguments
function VolumePartsMCMC(
  ::Type{TM}=PARAM[:mcmcgpu] ? CUDA.CuMatrix{T} : Matrix{T},
  ::Type{TV}=PARAM[:mcmcgpu] ? CUDA.CuVector{T} : Vector{T},
  ::Type{T}=PARAM[:mcmctype];
    cc::CapacityChains=error("cc is required"),
    Xv::AbstractXYIter=error("Xv is required"),
    ts::Vector{Int}=error("ts is required"),
    initial::Union{NamedTuple, Dict}=error("initial ")
    ) where {T<:Real, TV<:AbstractVector{T}, TM<:AbstractMatrix{T}}

  Ψ = CCState(cc, initial)
  #=params=(
    A₀ = TV(rand(T, dims.K)),
    G = TM(ones(T, dims.K, dims.T-1)),
    σ² = rand(T, dims.T-1),
    ω = rand(T, dims.T-1)
    )=#

  return VolumePartsMCMC(Ψ, cc, Xv, ts, TM, TV, T)
end

#deepcopy avoids copying the static objects
function Base.deepcopy(Θ::AbstractVolumePartsMCMC{TM, TV, T}) where {TM,TV,T}

  Ψ=deepcopy(Θ.Ψ)
  cc=deepcopy(Θ.cc)
  Xv=Θ.Xv
  ts=Θ.ts
  tsM1=Θ.tsM1

  return VolumePartsMCMC(Ψ, cc, Xv, ts, TM, TV, T, tsM1=tsM1)
end

#=
#the below funciton is now only for testing purposes
function (Θ::VolumePartsIter)(Xv::XYIterLevel, _Xv::AbstractXYIterAlloc{TM,TV,T},
  ::Val{:level}, RegressionType::Val, absμ::Tabsμ = _Xv.absμ_alloc
  ) where {TM<:Matrix,TV<:Vector,T,Tabsμ<:Vector}
  @unpack G,A = Θ

  #initialize A₁.==1 and all other As to the amount implied by G
  #print("prod(Θ.G): f0=$(prod(Θ.G|>Matrix)) x64=$(prod(Float64.(Θ.G |> Matrix) )) |")
  prevA₁ = A[:,1]
  factored::TM = initializeÃandfactor(Θ,Xv,_Xv)

  #run the regression for A₁
  try
    A[:,1] .= reg(factored, Xv.v, RegressionType)
  catch err
    (err == InterruptException()) && throw(err)
    errmessage = "$err"
    errmessage = length(errmessage > 1000) ? errmessage[1:1000] : errmessage
    all((isfinite.(factored))) || error("Factored includes nonfinite values!
      prevA₁: $(prevA₁)
      ****************\nG: $G
      ****************\nfactored[1:1000,:]: $(factored[1:1000,:])")

    A[:,1] .= reg(factored, Xv.v, Val{:fallbackreg}())
    @warn "regression method $RegressionType failed due to $errmessage. Using fallbackreg fallback val=$(A[:,1])"
  end

  #println("A[:,1]: $(A[:,1])")

  #absμ .= vec(sum(absμₖ,dims=2))
  absμ .= factored * (x->Float64(abs(x))).(A[:,1])

  return absμ
end

#if absμₖ is already formed, generating the demeaned predictions should be cheap
function (Θ::VolumePartsIter)(Xv::AbstractXYIter{TM,TV,T}, _Xv::AbstractXYIterAlloc{TM,TV,T},
  absμₖ::Tabsμₖ, ::Val{:leveldemean},
  ::Type{Tabsμ}=typeof(_Xv.absμ_alloc),
  ::Type{elTabsμₖ}=eltype(Tabsμₖ)) where {TM,TV,T,Tabsμₖ,Tabsμ,elTabsμₖ}

  dims = (K=size(Θ.A,1), T = size(Θ.A,2))
  absμ::Tabsμ = _Xv.absμ_alloc
  LinearAlgebra.BLAS.gemv!('N',elTabsμₖ(1.0),absμₖ,ones(elTabsμₖ,dims.K),zero(elTabsμₖ),absμ)

  return absμ
end

#primary prediction alogrithm for the demeaned version
function (Θ::VolumePartsIter)(Xv::XYIterDemean, _Xv::AbstractXYIterAlloc{TM,TV,T},
    ::Val{:leveldemean}, RegressionType::Val, absμ::Tabsμ = _Xv.absμ_alloc
    ) where {TM,TV,T,Tabsμ}
  @unpack G,A = Θ

  #initialize A₁.==1 and all other As to the amount implied by G
  #print("prod(Θ.G): f0=$(prod(Θ.G|>Matrix)) x64=$(prod(Float64.(Θ.G |> Matrix) )) |")
  prevA₁ = A[:,1]
  factored::TM = initializeÃandfactor(Θ,Xv,_Xv)

  #run the regression for A₁
  try
    A[:,1] .= reg(factored, Xv.ṽ, RegressionType)
  catch err
    (err == InterruptException()) && throw(err)
    errmessage = "$err"
    errmessage = length(errmessage > 1000) ? errmessage[1:1000] : errmessage
    all((isfinite.(factored))) || error("Factored includes nonfinite values!
      prevA₁: $(prevA₁)
      ****************\nG: $G
      ****************\nfactored[1:1000,:]: $(factored[1:1000,:])")

    A[:,1] .= reg(factored, Xv.ṽ, Val{:fallbackreg}())
    @warn "regression method $RegressionType failed due to $errmessage. Using fallback val=$(A[:,1])"
  end

  #println("A[:,1]: $(A[:,1])")

  #absμₖ = foldl(Completing(vcat), Map(M->vec(sum(M,1))), xabsμₖ)
  #absμ .= vec(sum(absμₖ,dims=2))
  absμ .= factored * (x->Float64(abs(x))).(A[:,1])

  return absμ
end

#this is a cheap version which relies on the already calculated absμₖ
(Θ::VolumePartsIter)(Xv, _Xv, absμₖ, ::Val{:leveldemean}) = (Θ::VolumePartsIter)(
  Xv, _Xv, absμₖ, Val{:level}())

(Θ::VolumePartsIter)(Xv::AbstractXYIter{TM,TV,T},
  _Xv::AbstractXYIterAlloc{TM,TV,T},
    ::Val{:leveldemean}) where {TM<:CuMatrix,TV<:CuVector,T} = error(
    "Do not use CUDA for the demeaning version")

function (Θ::VolumePartsIter)(Xv::XYIterLevel,
  _Xv::AbstractXYIterAlloc{TM,TV,T}, ::Val{:level}, RegressionType::Val,
  absμ::Tabsμ = _Xv.absμ_alloc
  ) where {TM<:CuMatrix,TV<:CuVector,T,Tabsμ<:CuArray}

  @unpack G,A = Θ

  #initialize A₁.==1 and all other As to the amount implied by G
  #print("prod(Θ.G): f0=$(prod(Θ.G|>Matrix)) x64=$(prod(Float64.(Θ.G |> Matrix) )) |")
  prevA₁ = A[:,1]
  factored::TM = initializeÃandfactor(Θ,Xv,_Xv)

  #run the regression for A₁
  try
    A[:,1] .= reg(factored, Xv.v, RegressionType)
    #@info "culevel A[:,1]: $(A[:,1])"
    #@info "culevel factored[20,:]: $(factored[20,:])"
    #@info "culevel A[:,20]: $(A[:,20])"
  catch err
    (err == InterruptException()) && throw(err)
    errmessage = "$err"
    errmessage = length(errmessage > 1000) ? errmessage[1:1000] : errmessage
    all((isfinite.(factored))) || error("Factored includes nonfinite values!
      prevA₁: $(prevA₁)
      ****************\nG: $G
      ****************\nfactored[1:1000,:]: $(factored[1:1000,:])")

    A[:,1] .= reg(factored, Xv.v, Val{:fallbackreg}())
    @warn "regression method $RegressionType failed due to $errmessage. Using fallback val=$(A[:,1])"
  end

  #println("A[:,1]: $(A[:,1])")

  absμ .= factored * (x->Float64(abs(x))).(A[:,1])

  return absμ
end

#this is a cheap version which relies on the already calculated absμₖ
function (Θ::VolumePartsIter)(Xv::XYIterLevel,
    _Xv::AbstractXYIterAlloc{TM,TV,T},
    absμₖ::Tabsμₖ, ::Type{elTabsμₖ}=eltype(Tabsμₖ)
    ) where {TM<:CuMatrix,TV<:CuVector,T, Tabsμₖ<:CuMatrix, elTabsμₖ}

  #first map each coefficient to its appropriate
  dims = (K=size(Θ.A,1), T = size(Θ.A,2))
  absμ::CuVector{elTabsμₖ} = _Xv.absμ_alloc
  CUBLAS.gemv!('N',elTabsμₖ(1.0),absμₖ,CUDA.ones(elTabsμₖ,dims.K),zero(elTabsμₖ),absμ)

  return absμ
end

#computs the predicted values for each strategy k and time t and demeans the result
#via the projection matrix
function projectvolumefromA(Θ::VolumePartsIter, Xv::XYIterDemean,
  _Xv::AbstractXYIterAlloc{TM,TV,T}, absμₖ::Tabsμₖ = _Xv.absμₖ_alloc) where {TM,TV,T,Tabsμₖ}

  @unpack A, LA = Θ.expand
  @unpack ws, RLws, M = Xv
  #first map each coefficient to its appropriate

  absμₖ::Tabsμₖ = _Xv.absμₖ_alloc
  absμₖ .= genabsμₖ.(A', LA', ws, RLws)

  multiply!(M, absμₖ)
  return absμₖ
end

#same as the above without demeaning
function projectvolumefromA(Θ::VolumePartsIter, Xv::XYIterLevel,
  _Xv::AbstractXYIterAlloc{TM,TV,T}, absμₖ::Tabsμₖ = _Xv.absμₖ_alloc) where {TM,TV,T,Tabsμₖ}

  @unpack A, LA = Θ.expand
  @unpack ws, RLws= Xv
  #first map each coefficient to its appropriate

  absμₖ::Tabsμₖ = _Xv.absμₖ_alloc
  broadcast!(genabsμₖ, absμₖ, A', LA', ws, RLws)

  return absμₖ
end

function projectvolumefromA(Θ::VolumePartsIter, Xv::XYIterLevel,
  _Xv::AbstractXYIterAlloc{TM,TV,T}, absμₖ::Tabsμₖ = _Xv.absμₖ_alloc
  ) where {TM<:CuMatrix, TV<:CuVector, T, Tabsμₖ}

  @unpack A, LA = Θ.expand
  @unpack ws, RLws = Xv
  #first map each coefficient to its appropriate

  dims = (K=size(Θ.A,1), T = size(Θ.A,2))

  broadcast!(cudagenabsμₖ, absμₖ, A', LA', ws, RLws)
  #absμₖ .= cudagenabsμₖ.(A', LA', ws, RLws)

  return absμₖ
end

#Updates G based on values of A
#of course, this changes the RHS side of teh regressions, so any
#factored cross-sections need to be recalculated
updateGfromA!(args...;kwargs...) = error("I dropped updateGfromA!. If needed, just replace
  with Θ.A[:,2:end] ./ Θ.A[:,1:(end-1)] (though be careful about Ã vs A)")

#uses a base value of A and the G matrix to update all other values of A
function updateAfromGAₜ!(Θ::VolumePartsIter, Xv::AbstractXYIter{TM,TV,T},
  _Xv::AbstractXYIterAlloc{TM,TV,T}, t::Int,
  boundA::Bool = PARAM[:iterboundA]) where {TM,TV,T}

  local A₀::TV
  local Aₜ::TV = Θ.A[:,t]
  dims = (T=size(Θ.A,2),K=size(Θ.A,1))

  #compute the cumulative growth
  ΠG::TM = _Xv.a_alloc
  ΠG[:,1] .= T(1.0)
  ΠG[:,2:end] .= cumprodbyrow(Θ.G, Xv.Ascale)

  #if boundA #maintain consistency- comprodrow may adjust for small values
  #  Θ.G .= ΠG[:,2:end] ./ ΠG[:,1:(end-1)]
  #end

  A₀ = Aₜ ./ ΠG[:,t] #identify the intial value of assets

  Θ.A .= A₀ .* ΠG #use the growth to compute A
  (vec(Θ.A[:,t]) ≈ Aₜ) || error("vec(Θ.A[:,t]) ($(vec(Θ.A[:,t]))) ≠ Aₜ ($(Aₜ))! ") #quick integrity check
  Θ.A[:,t] .= Aₜ #preserve exact values to partly mitigate roundoff error

  return nothing
end

function updateAfromGAₜ!(Θ::VolumePartsIter, Xv::AbstractXYIter{TM,TV,T},
  _Xv::AbstractXYIterAlloc{TM,TV,T}, t::Int,
  boundA::Bool = PARAM[:iterboundA]) where {TM<:CuMatrix,TV<:CuVector,T}

  local Aₜ::TV = Θ.A[:,t]
  dims = (T=size(Θ.A,2),K=size(Θ.A,1))

  #compute the cumulative growth
  ΠG::TM = _Xv.a_alloc
  ΠG[:,1] .= T(1.0)
  ΠG[:,2:end] .= cumprodbyrow(Θ.G, Xv.Ascale)
  #ΠG .= hcat(CUDA.ones(T, dims.K, 1), cumprod(Θ.G,dims=2))

  #A₀ = Aₜ ./ ΠG[:,t] #identify the intial value of assets
  #Θ.A .= A₀ .* ΠG
  A = (Aₜ ./ ΠG[:,t]) .* ΠG
  copyto!(Θ.A, A) #use the growth to compute A
  #@assert vec(Θ.A[:,t]) ≈ Aₜ #quick integrity check
  Θ.A[:,t] .= Aₜ #preserve exact values to partly mitigate roundoff error

  return nothing
end

#NOTE: could be used in some future GPU version
cudaGᵢ(Aᵢ,LAᵢ)= Aᵢ/LAᵢ
cudafactorAᵢ(wsᵢ, RLwsᵢ, Gᵢ) = CUDA.abs(wsᵢ - RLwsᵢ/Gᵢ)
cudafactorLAᵢ(wsᵢ, RLwsᵢ, Gᵢ) = CUDA.abs(Gᵢ*wsᵢ - RLwsᵢ)

#takes a volume part object and writes the results to a df
function formresults(panel::AbstractDataFrame, ms::MeasureSpec,
  Θ::VolumePartsIter,
  grossexposures::AbstractMatrix,
  growthrates::AbstractMatrix)

  Fcoefs::Vector{Symbol} = Fcoefficients0(ms)
  Fgcoefs::Vector{Symbol} = (s->Symbol(:G,s)).(Fcoefs)

  tidx::Dict = gettidx(panel, ms)
  tidxinv::Vector{Date} = sort(collect(keys(tidx)))

  if length(Fcoefs) ≠ size(Θ.A,1)
    error("Coeffcients not in expected form. Fcoefs: $Fcoefs")
  end

  results::DataFrame = DataFrame(date=tidxinv)
  resultsmat = hcat(grossexposures,growthrates)
  results = hcat(results, DataFrame(resultsmat, [Fcoefs; Fgcoefs;]))

  return results

end
=#

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
    abstract type AbstractDice end

	abstract type AbstractCounter end

#Base
push!(gs[:overlay][:dice], gs[:stage][:six_sided_die])

# "sticky" counters
dice_and_counters = [ dc.faces[begin] for dc in [
    values(gs[:overlay][:dice])...,
    values(gs[:overlay][:dice])...
    ] if SDL2.HasIntersection(
        Ref(ib[end].position),
        Ref(dc.faces[begin].position)) && !(ib[end] in [
        values(gs[:overlay][:dice])...,
        values(gs[:overlay][:dice])... ]
        )
    ]

for dc in dice_and_counters
    dc.data[:mouse_offset] = [
        dc.x - gs[:MOUSE_POS][1], dc.y - gs[:MOUSE_POS][2] ]
end

push!(gs[:group][:selected], [ ib[end] ]...) #, dice_and_counters... )

for dc in [ gs[:stage][:six_sided_die], gs[:stage][:glass_counter] ]

    if ib[end].data[:parent_id] == dc.id

        for f in dc.faces
            f.position = ib[end].position
        end

        dc.faces = circshift(dc.faces, 1)
        filter!(x->x!==ib[end], gs[:group][:clickables])
        push!(gs[:group][:clickables], dc.faces[begin])
    end
end

#clickables
[ d.faces[begin] for d in values(gs[:overlay][:dice]) ]...,
[ c.faces[begin] for c in values(gs[:overlay][:counters]) ]...,

#on mouse up
elseif ib[end] in values(gs[:overlay][:dice])

    if g.keyboard.LCTRL || g.keyboard.RCTRL
        push!(gs[:group][:selected], ib[end])

    else
        copy = GIF("die_$(randstring(5))", ib[end].data[:gif],
            x=gs[:MOUSE_POS][1], y=gs[:MOUSE_POS][2])
        push!(gs[:overlay][:dice], copy)
        push!(gs[:group][:selected], copy)
    end

!(a in [ c.faces[begin] for c in values(gs[:overlay][:counters]) ]) &&
!(a in [ d.faces[begin] for d in values(gs[:overlay][:dice]) ])

#overlay
:dice => Dice[],
:counters => Counter[],

#game rules
# ╔═╡ 94d2e13c-55ec-11eb-0c93-91aa4e14ceff
mutable struct Dice <: AbstractDice
	id::String
	faces::Vector{Actor}
	sides::Int
end

# ╔═╡ 3819f224-558c-11eb-2498-95b826027aa2
mutable struct Counter <: AbstractCounter
	id::String
	faces::Vector{Actor}
	on::Union{Card,Player,Nothing}
end

# game stage
dice_faces = [ load("$gd/Base/ui/dice/$fn") for fn in
    readdir("$gd/Base/ui/dice") if occursin("gif", fn) ]

serialize("$gd/Base/ui/dice/dice_faces.jls", dice_faces)

dfs = deserialize("$gd/Base/ui/dice/dice_faces.jls")

function create_die(dfs; id="die_$(randstring(5))", x=0, y=0)
    d = Dice(id,
        [ GIF("dieface_$(randstring(5))", f) for f in dfs ],
        length(dfs),
        )

    for f in d.faces
        f.data[:parent_id] = id
        f.x = x
        f.y = y
    end

    return d
end

function create_glass_counter(id="ctr_$(randstring(5))"; x=0, y=0)
    c = Counter("ctr_$(randstring(5))",
        [ GIF("ctrface_$(randstring(5))", load("$gd/Base/ui/counters/glass.gif")) ],
        nothing,
        )
    c.faces[begin].data[:parent_id] = id
    c.faces[begin].x = x
    c.faces[begin].y = y
    return c
end

#STAGE
:six_sided_die => create_die(dfs,
    x=ceil(Int32, 0.85SCREEN_WIDTH),
    y=ceil(Int32, 0.9SCREEN_HEIGHT)
    ),
:glass_counter => create_glass_counter(
    x=ceil(Int32, 0.8SCREEN_WIDTH),
    y=ceil(Int32, 0.9SCREEN_HEIGHT)
    ),


text:
# Some utility functions to make it more convenient to explore the
# type hierarchy.

using InteractiveUtils

export allsubtypes, showsubtypes, pedigree

"""
    allsubtypes(t::Type)
Return a Vector of `t` and all of its subtypes.
"""
function allsubtypes(t::Type, result=Vector{Type}())
    push!(result, t)
    for st in subtypes(t)
	allsubtypes(st, result)
    end
    return result
end

"""
   showsubtypes(t::Type)
Print a hierarchical list of `t` and all of its subtypes.
"""
function showsubtypes(t::Type, level=0)
    indent1 = "  "
    println("$(repeat(indent1, level))$t")
    for st in subtypes(t)
        showsubtypes(st, level + 1)
    end
end

"""
    pedigree(t::Type)
Return a Vector of `t` and its supertypes.
"""
function pedigree(t::Type, result=Vector{Type}())
    push!(result, t)
    if supertype(t) != t
        pedigree(supertype(t), result)
    end
    result
end


text:
export PendingScheduler, release!

import Base: show, similar

mutable struct PendingActor{L, A} <: Actor{L}
    actor :: A
    last  :: Union{Nothing, L}
end

make_pending_actor(::Type{L}, actor::A) where { L, A } = PendingActor{L, A}(actor, nothing)

getlast(actor::PendingActor)        = actor.last
setlast!(actor::PendingActor, last) = actor.last = last

function release!(actor::PendingActor, reset::Bool)
    last = getlast(actor)
    if last !== nothing
        next!(actor.actor, getlast(actor))
        if reset
            setlast!(actor, nothing)
        end
    end
end

struct PendingScheduler <: AbstractScheduler
    reset   :: Bool
    pending :: Vector{PendingActor}

    PendingScheduler(reset::Bool = true) = new(reset, Vector{PendingActor}())
end

Base.show(io::IO, ::PendingScheduler) = print(io, "PendingScheduler()")

Base.similar(::PendingScheduler) = PendingScheduler()

makeinstance(_, scheduler::PendingScheduler) = scheduler

instancetype(_, ::Type{ <: PendingScheduler }) = PendingScheduler

release!(scheduler::PendingScheduler) = foreach((d) -> release!(d, scheduler.reset), scheduler.pending)

function scheduled_next!(actor::PendingActor{L}, data::L, ::PendingScheduler) where L
    setlast!(actor, data)
    return nothing
end

function scheduled_error!(actor::PendingActor, err, instance::PendingScheduler)
    release!(actor, instance.reset)
    deregister!(instance, actor)
    error!(actor.actor, err)
end

function scheduled_complete!(actor::PendingActor, instance::PendingScheduler)
    release!(actor, instance.reset)
    deregister!(instance, actor)
    complete!(actor.actor)
end

struct PendingSchedulerSubscription <: Teardown
    subscription :: Teardown
    scheduler    :: PendingScheduler
    actor        :: PendingActor
end

as_teardown(::Type{ <: PendingSchedulerSubscription }) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::PendingSchedulerSubscription)
    unsubscribe!(subscription.subscription)
    deregister!(subscription.scheduler, subscription.actor)
    return nothing
end

function scheduled_subscription!(source::S, actor, instance::PendingScheduler) where S
    pending = make_pending_actor(subscribable_extract_type(S), actor)
    register!(instance, pending)
    subscription = on_subscribe!(source, pending, instance)
    return PendingSchedulerSubscription(subscription, instance, pending)
end

function register!(scheduler::PendingScheduler, actor::PendingActor)
    push!(scheduler.pending, actor)
end

function deregister!(scheduler::PendingScheduler, actor::PendingActor)
    filter!(a -> a !== actor, scheduler.pending)
end


text:
#
println("Ed) DR rate coefficients for neon-like Ar: AverageSCA model.")
using JLD
#
if  true
    # Calculation of the Ne^+  2s^2 2p^6; 2s, 2p --> 3s, 3p, 3d, ... 5d electron capture and subsequent autoionization and stabilization
    setDefaults("print summary: open", "Ne-like-DR-rate.sum")
    setDefaults("method: continuum, Galerkin");            setDefaults("method: normalization, pure sine")
    grid = Radial.Grid(Radial.Grid(false), rnt = 1.0e-5, h = 5.0e-2, hp = 2.0e-2, rbox = 15.0)
    name = "DR rate coefficient for neon-like argon and for temperatures < 500 K"
    asfSettings = AsfSettings(AsfSettings(), maxIterationsScf = 36)
    wa   = Cascade.Computation(Cascade.Computation(); name=name, nuclearModel=Nuclear.Model(10.), grid=grid, asfSettings=asfSettings,
                               approach=Cascade.AverageSCA(),
                               scheme=Cascade.DrRateCoefficientScheme([ElecCapture(), Auger(), Radiative()], [E1], 2.5, 0.0, 1, 3, 2,
                                                                      [Shell("2s"), Shell("2p")]),
                               initialConfigs=[Configuration("1s^2 2s^2 2p^6")] )
    println(wa)
    wb = perform(wa; output=true)
    setDefaults("print summary: close", "")
    #
elseif false
    # Calculation of the Ne^+  2s^2 2p^5; photo-ionization cross sections
    setDefaults("print summary: open", "Ne-plus-photoionization.sum")

    name = "Photoionization of Si- "
    grid = Radial.Grid(Radial.Grid(false); rnt = 3.0e-6, h = 2.0e-2, hp = 3.0e-2, NoPoints=1110)
    wa   = Cascade.Computation(Cascade.Computation(); name=name, nuclearModel=Nuclear.Model(10.), grid=grid, approach=Cascade.AverageSCA(),
                               scheme=Cascade.PhotonIonizationScheme([Photo()], 1, [1.0, 3.0, 5.0]),
                               initialConfigs=[Configuration("1s^2 2s^2 2p^5")] )
    println(wa)
    wb = perform(wa; output=true)
    setDefaults("print summary: close", "")

else
    # Simulation of the Ne^+  2s^2 2p^5; photo-ionization cross sections
    setDefaults("print summary: open", "Ne-plus-photoabsorption-simulation.sum")

    data = [JLD.load("zzz-cascade-ionizing-computations-2020-07-10T19.jld"), JLD.load("zzz-cascade-excitation-computations-2020-07-10T19.jld")]
    ##x JLD2.@load "zzz-cascade-ionizing-computations-2020-06-23T09.jld"
    ##x resIon  = results
    ##x JLD2.@load "zzz-cascade-excitation-computations-2020-06-23T09.jld"
    ##x resExc = results
    ##x data = [resExc]

    name = "Simulation of photoabsorption for Ne^+"

    wc   = Cascade.Simulation(Cascade.Simulation(), name=name, property=Cascade.PhotoAbsorption(),
                              settings=Cascade.SimulationSettings(0., 0., 1., 4., 3., [(1, 1.0)]), computationData=data )
    println(wc)
    wd = perform(wc; output=true)
    setDefaults("print summary: close", "")

end


text:
struct NPMLE{T<:Distribution}
    prior_grid::Vector{Float64}
    prior_mixing::Vector{Float64}
    marginal_grid::Vector{Float64}
    marginal::T
    σ::Float64
end

function npmle_dual_hist(prior_grid, marginal_grid, Xs; L=0.0, θ=0.0, maxiter=1_000_000,
     σ_std=1.0)

    h = marginal_grid[2] - marginal_grid[1]
    edges = [marginal_grid-h/2; maximum(marginal_grid)+h/2]
    edges[1] = -Inf
    edges[end] = Inf

    hist = fit(Histogram, Xs, edges; closed=:right)

    idx = findall( hist.weights .> 0)

    marginal_grid_active = marginal_grid[idx]
    #A = NormalConvolutionMatrix(prior_grid,Xs)
    A = convolution_matrix(DiscretizedNormalConvolutionProblem, prior_grid, marginal_grid_active; σ=σ_std);

    if L==0.0
        L = real(tsvd(A)[2][1])
        @show L
    end
    σ = 1/L/1.1
    τ = 1/L/1.1

    N = sum(hist.weights)
    m = length(prior_grid)

    wts = hist.weights[idx]

    x_n = wts./pdf.(Normal(0,std(Xs)),marginal_grid_active)# ν we are interested in
    x_n_tmp = copy(x_n)

    y_n = zeros(prior_grid)
    y_n[ceil(Int,m/2)] = 1.0
    y_n_tmp = copy(y_n)

    bar_x_n = copy(x_n)

    function prox_σFstar(x, σ)
        x - σ*max(0,min(x/σ, N))
    end

    function prox_τG(y, ct, τ )
        (y + sqrt(y^2 + 4τ*ct))/2
    end

    for i=1:maxiter
       At_mul_B!(y_n_tmp, A, bar_x_n)
       y_n .+= σ.*y_n_tmp
       y_n .= prox_σFstar.(y_n, σ)

       bar_x_n .= x_n
       idx = findall(y_n .> 0)

       A_mul_B!(x_n_tmp, view(A, :, idx), view(y_n, idx))
       x_n .= x_n .- τ.*x_n_tmp
       x_n .= prox_τG.(x_n, wts, τ)

       bar_x_n .= (1+θ).*x_n .- θ.*bar_x_n

       if mod(i-1,10000) == 0
           if i>=20_000 && (abs(sum(y_n)-1) <= eps_tol) && norm(x_n_tmp .- wts./x_n, 2) <= 2*eps_tol
                break
           end
       #    @show sum(h.*wts./x_n[1:end])
        end
    end

    #A_full=convolution_matrix(DiscretizedNormalConvolutionProblem, prior_grid, marginal_grid);

    #(y_n, A_full*y_n)
    y_n
end

function StatsBase.fit(::Type{NPMLE}, prior_grid, marginal_grid, Xs; σ=0.0, args...)
    σ_std = sqrt(σ^2 + 1)
   res= npmle_dual_hist(prior_grid, marginal_grid, Xs; σ_std=σ_std, args...)
   # add some checks that everything went well..
   idx = findall(res .> 0)
   prior_grid = prior_grid[idx]
   prior_mixing = res[idx]
   prior_mixing = prior_mixing/sum(prior_mixing)
   m = MixtureModel(Normal, collect(zip(prior_grid, σ_std*ones(prior_grid))), prior_mixing)
   NPMLE(prior_grid, prior_mixing, marginal_grid, m, σ)
end

function MixingNormalConvolutionProblem(a::NPMLE, marginal_grid)
    priors = a.marginal.components
    mixing_coef = a.marginal.prior.p
    MixingNormalConvolutionProblem(priors ,marginal_grid, mixing_coef)
end

function NormalConvolutionProblem(a::NPMLE, marginal_grid)
   m =  MixingNormalConvolutionProblem(a ,marginal_grid)
   NormalConvolutionProblem(m, a.prior_mixing)
end


text:
#Written By Johnathan Bizzano
module Data

    include("DataSet.jl")
    include("MemoryDataSet.jl")
    include("DataSetReader.jl")
    include("LazyDataSet.jl")
    include("RandomBatch.jl")
   # include("DataSetWriter.jl")
   # include("Sorter.jl")

end


text:
"""
    loglikelihood!(m, [needgrad=false], [needhess=false])

Evaluate the log-likelihood and optionally gradient and Hessian of
ordered multinomial model.

# Positional arguments
- `m::PolrModel`: a `PolrModel` type. Log-likelihood is evaluated based on
fields `Y`, `X`, `α`, `β`, and `link` of `m`.
- `needgrad::Bool=false`: evaluate gradient or not.
- `needhess::Bool=false`: evaluate Fisher information matrix (FIM) (negative Hessian) or not.

# Output
- `logl`: log-likelihood. If `needgrad=true`, field `m.∇` is overwritten by the
gradient (score) with respect to `(θ, β)`. If `needhess=true`, field `m.FIM` is
overwritten by the Fisher information matrix (FIM), equivalent to negative Hessian,
with respect to `(θ, β)`.
"""
function loglikelihood!(
    m::OrdinalMultinomialModel,
    needgrad = false,
    needhess = false)

    T = eltype(m.X)
    needgrad && fill!(m.∇, 0)
    needhess && fill!(m.FIM, 0)
    # update η according to X and β
    mul!(m.η, m.X, m.β)
    # update θ according to α
    # α[2] to α[J-2] overwritten by their exponential
    fill!(m.dθdα, 0)
    m.θ[1] = m.α[1]
    m.dθdα[1, :] .= 1
    for j in 2:m.J-1
        m.θ[j] = m.θ[j-1] + exp(m.α[j])
        m.dθdα[j, j:end] .= exp(m.α[j])
    end
    minη, maxη = convert(T, -100), convert(T, 100)
    logl = zero(T)
    for obs in 1:m.n
        wtobs = isempty(m.wts) ? one(T) : m.wts[obs]
        wtobs == 0 && continue
        ### log-likelihood
        yobs = m.Y[obs]
        ηy   = yobs == m.J ? maxη : m.θ[yobs]   - m.η[obs]
        ηym1 = yobs == 1   ? minη : m.θ[yobs-1] - m.η[obs]
        ηy   = clamp(ηy  , minη, maxη) # for numerical stability
        ηym1 = clamp(ηym1, minη, maxη) # ym1 reads y-minus-one
        cumproby   = GLM.linkinv(m.link, ηy)
        cumprobym1 = GLM.linkinv(m.link, ηym1)
        py = cumproby - cumprobym1
        logl += py ≤ 0 ? typemin(T) : wtobs * log(py)
        ### gradient
        needgrad || continue
        if py ≤ 0
            # derivative of θ
            yobs < m.J && (m.∇[yobs  ] += typemax(T))
            yobs > 1   && (m.∇[yobs-1] += typemin(T))
            # derivative of β = wt[obs] * X[:. obs]
            m.wt∂β[obs] = zero(T)
        else
            # derivative of θ
            pyinv    = inv(py)
            derivy   = GLM.mueta(m.link, ηy)
            derivym1 = GLM.mueta(m.link, ηym1)
            derivΔ = derivy - derivym1
            yobs < m.J && (m.∇[yobs  ] += wtobs * pyinv * derivy)
            yobs > 1   && (m.∇[yobs-1] -= wtobs * pyinv * derivym1)
            # derivative of β = wt[obs] * X[:. obs]
            m.wt∂β[obs] = - pyinv * derivΔ
        end
        ### hessian
        needhess || continue
        ηjm1 = ηj = ηjp1 = zero(T)
        cumprobjm1 = cumprobj = cumprobjp1 = zero(T)
        derivjm1 = derivj = derivjp1 = zero(T)
        for j in 1:m.J-1
            # systematic components at j-1, j, j+1
            ηjm1       = j == 1 ? minη : ηj
            ηj         = j == 1 ? m.θ[1] - m.η[obs] : ηjp1
            ηjp1       = j == m.J-1 ? maxη : m.θ[j+1] - m.η[obs]
            # cumulative probabilities at j-1, j, j+1
            cumprobjm1 = j == 1 ? zero(T) : cumprobj
            cumprobj   = j == 1 ? GLM.linkinv(m.link, ηj) : cumprobjp1
            cumprobjp1 = j == m.J-1 ? one(T) : GLM.linkinv(m.link, ηjp1)
            # point mass at j-1, j, j+1
            pj         = cumprobj - cumprobjm1
            pjp1       = cumprobjp1 - cumprobj
            pjinv      =   pj ≤ 0 ? zero(T) : inv(pj)
            pjp1inv    = pjp1 ≤ 0 ? zero(T) : inv(pjp1)
            # mueta at j-1, j, j+1
            derivjm1   = j == 1 ? zero(T) : derivj
            derivj     = j == 1 ? GLM.mueta(m.link, ηj) : derivjp1
            derivjp1   = j == m.J-1 ? zero(T) : GLM.mueta(m.link, ηjp1)
            # FIM for ∂θ^2
            m.FIM[j, j] += wtobs * (pjinv + pjp1inv) * derivj * derivj
            j > 1 && (m.FIM[j, j-1] -= wtobs * pjinv * derivj * derivjm1)
            # FIM for ∂θ∂β, weights only
            m.wt∂θ∂β[obs, j] += (pjinv * derivjm1 - (pjinv + pjp1inv) * derivj +
                pjp1inv * derivjp1) * derivj
            # FIM for ∂β∂β, weights only
            m.wt∂β∂β[obs] += pjinv * (derivj - derivjm1)^2
            j == m.J - 1 &&  (m.wt∂β∂β[obs] += pjp1inv * derivj * derivj) # contribution from p_iJ
        end
    end
    if needgrad
        if isempty(m.wts)
            @views mul!(m.∇[m.J:end], transpose(m.X), m.wt∂β)
        else
            m.wtwk .= m.wts .* m.wt∂β
            @views mul!(m.∇[m.J:end], transpose(m.X), m.wtwk)
        end
    end
    if needhess
        if isempty(m.wts)
            # ∂θ∂β
            @views mul!(m.FIM[m.J:end, 1:m.J-1], transpose(m.X), m.wt∂θ∂β)
            # ∂β∂β
            copyto!(m.scratchm1, m.X)
            lmul!(Diagonal(m.wt∂β∂β), m.scratchm1)
            @views mul!(m.FIM[m.J:end, m.J:end], transpose(m.X), m.scratchm1)
        else
            # ∂θ∂β
            copyto!(m.scratchm1, m.X)
            lmul!(Diagonal(m.wts), m.scratchm1)
            @views mul!(m.FIM[m.J:end, 1:m.J-1], transpose(m.scratchm1), m.wt∂θ∂β)
            # ∂β∂β
            m.wtwk .= m.wts .* m.wt∂β∂β
            copyto!(m.scratchm1, m.X)
            lmul!(Diagonal(m.wtwk), m.scratchm1)
            @views mul!(m.FIM[m.J:end, m.J:end], transpose(m.X), m.scratchm1)
        end
        LinearAlgebra.copytri!(m.FIM, 'L') # copy lower triangular to upper triangular
    end
    logl
end

"""
    polr(formula, df, link, solver=NLoptSolver(algorithm=:LD_SLSQP, maxeval=4000))
    polr(X, y, link, solver=NLoptSolver(algorithm=:LD_SLSQP, maxeval=4000))

Fit ordered multinomial model by maximum likelihood estimation.

# Positional arguments
- `formula::Formula`: a model formula specifying responses and regressors.
- `df::DataFrame`: a dataframe. Response variable should take integer values
    starting from 1.
- `y::Vector`: integer vector taking values in `1,...,J`.
- `X::Matrix`: `n x p` covariate matrix excluding intercept.
- `link::GLM.Link`: `LogitLink()` (default), `ProbitLink()`, `CauchitLink()`, or `CloglogLink()`.
- `solver`: An instance of `Ipopt.Optimizer()` or `NLopt.Optimizer()` (default).

# Keyword arguments
- `wts::AbstractVector=similar(X, 0)`: observation weights.

# Output
- `dd:PolrModel`: a `PolrModel` type.
"""
polr(X::AbstractMatrix, y::AbstractVector, args...; kwargs...) =
    fit(AbstractOrdinalMultinomialModel, X, y, args...; kwargs...)
polr(f::FormulaTerm, df, args...; kwargs...) =
    fit(AbstractOrdinalMultinomialModel, f, df, args...; kwargs...)

"""
    fit(AbstractOrdinalMultinomialModel, X, y, link, solver)

Fit ordered multinomial model by maximum likelihood estimation.

# Input
- `y::Vector`: integer vector taking values in `1,...,J`.
- `X::Matrix`: `n x p` covariate matrix excluding intercept.
- `link::GLM.Link`: `LogitLink()`, `ProbitLink()`, `CauchitLink()`, or `CloglogLink()`.
- `solver`: An instance of `Ipopt.Optimizer()` or `NLopt.Optimizer()` (default).

# Output
- `dd:OrdinalMultinomialModel`: an `OrdinalMultinomialModel` type.
"""
function fit(
    ::Type{M},
    X::AbstractMatrix,
    y::AbstractVecOrMat,
    link::GLM.Link = LogitLink(),
    solver::TSOLVER = NLopt.Optimizer();
    wts::AbstractVector = similar(X, 0)
    ) where {M<:AbstractOrdinalMultinomialModel, TSOLVER<:MOI.AbstractOptimizer}
    T = eltype(X)
    ydata = Vector{Int}(undef, size(y, 1))
    # set up optimization
    MOI.empty!(solver)
    # This is a hack since we can't initialize a solver with parameters already set
    if typeof(solver) == NLopt.Optimizer
        algo = MOI.get(solver, MOI.RawOptimizerAttribute("algorithm"))
        if algo == :none
            MOI.set(solver, MOI.RawOptimizerAttribute("algorithm"), :LD_SLSQP)
            MOI.set(solver, MOI.RawOptimizerAttribute("max_iter"), 4000)
        end
    end

    if size(y, 2) == 1
        ydata = denserank(y)
    else #y is encoded via dummy-encoding
        for i in 1:size(y, 1)
            idx = findfirst(view(y, i, :) .== 1)
            ydata[i] = idx == nothing ? 1 : idx + 1
        end
    end
    #ydata = denserank(y) # dense ranking of y, http://juliastats.github.io/StatsBase.jl/stable/ranking.html#StatsBase.denserank
    dd = OrdinalMultinomialModel(X, ydata, convert(Vector{eltype(X)}, wts), link)
    lb = T[]
    ub = T[]
    NLPBlock = MOI.NLPBlockData(
        MOI.NLPBoundsPair.(lb, ub), dd, true
    )
    # initialize from LS solution
    β0 = [ones(length(ydata)) X] \ ydata
    par0 = [β0[1] - dd.J / 2 + 1; zeros(dd.J - 2); β0[2:end]]

    params = MOI.add_variables(solver, dd.npar)
    for i in 1:dd.npar
        MOI.set(solver, MOI.VariablePrimalStart(), params[i], par0[i])
    end

    MOI.set(solver, MOI.NLPBlock(), NLPBlock)
    MOI.set(solver, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(solver)
    # output
    stat = MOI.get(solver, MOI.TerminationStatus())
    stat == MOI.LOCALLY_SOLVED || @warn("Optimization unsuccessful; got $stat")
    xsol = similar(par0)
    for i in eachindex(xsol)
        xsol[i] = MOI.get(solver, MOI.VariablePrimal(), MOI.VariableIndex(i))
    end
    copyto!(dd.α, 1, xsol, 1, dd.J - 1)
    copyto!(dd.β, 1, xsol, dd.J, dd.p)
    loglikelihood!(dd, true, true)
    FIMbk = bunchkaufman(Symmetric(dd.FIM), check=false)
    if issuccess(FIMbk)
        dd.vcov[:] = inv(FIMbk)
    else # FIM is singular
        dd.vcov .= Inf
    end
    return dd

end

function MathOptInterface.initialize(
    m::OrdinalMultinomialModel,
    requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in MOI.features_available(m))
            error("Unsupported feature $feat")
        end
    end
end

MathOptInterface.features_available(m::OrdinalMultinomialModel) = [:Grad]

function MathOptInterface.eval_objective(
    m::OrdinalMultinomialModel,
    par::Vector)
    copyto!(m.α, 1, par, 1, m.J - 1)
    copyto!(m.β, 1, par, m.J, m.p)
    loglikelihood!(m, false, false)
end

function MathOptInterface.eval_objective_gradient(
    m::OrdinalMultinomialModel,
    grad::Vector,
    par::Vector)
    copyto!(m.α, 1, par, 1, m.J - 1)
    copyto!(m.β, 1, par, m.J, m.p)
    loglikelihood!(m, true, false)
    @views mul!(grad[1:m.J-1], m.dθdα, m.∇[1:m.J-1])
    copyto!(grad, m.J, m.∇, m.J, m.p)
end

function StatsModels.coeftable(mod::StatsModels.TableRegressionModel{T, S}
    where {T <: OrdinalMultinomialModel, S <: Matrix} )
    ct = coeftable(mod.model)
    cfnames = [["intercept$i|$(i+1)" for i in 1:(mod.model.J - 1)]; coefnames(mod)]
    if length(ct.rownms) == length(cfnames)
        ct.rownms = cfnames
    end
    ct
end