text: const SYMATA_VERSION = v"0.4.5"


text:
function update_and_get_max_abs_diff!(  pssrb::PotentialSimulationSetupRB{T, S, DAT, N1, N2},
                                        depletion_handling::Val{depletion_handling_enabled},
                                        only2d::Val{only_2d} = Val{false}(),
                                        is_weighting_potential::Val{_is_weighting_potential} = Val{false}(),
                                        use_nthreads::Int = Base.Threads.nthreads()
                                        )::T where {T, S, DAT, N1, N2, depletion_handling_enabled, only_2d, _is_weighting_potential}
    tmp_potential::Array{T, N2} = copy(pssrb.potential)
    if depletion_handling_enabled
        update!(pssrb, use_nthreads = use_nthreads, depletion_handling = depletion_handling, only2d = only2d, is_weighting_potential = is_weighting_potential)
        slopes::Array{T, N2} = tmp_potential - pssrb.potential
        @inbounds for i in 1:19
            tmp_potential[:] = pssrb.potential[:]
            update!(pssrb, use_nthreads = use_nthreads, depletion_handling = depletion_handling, only2d = only2d, is_weighting_potential = is_weighting_potential)
            slopes += tmp_potential - pssrb.potential
        end
        @inbounds slopes /= 20
        return maximum(abs.(slopes))
    else
        for i in 1:10
            update!(pssrb, use_nthreads = use_nthreads, depletion_handling = depletion_handling, only2d = only2d, is_weighting_potential = is_weighting_potential)
        end
        max_diff::T = maximum(abs.(tmp_potential - pssrb.potential))
        return max_diff
    end
end

function _update_till_convergence!( pssrb::PotentialSimulationSetupRB{T, S, 3},
                                    convergence_limit::T,
                                    ::Type{Array};
                                    n_iterations_between_checks = 500,
                                    depletion_handling::Val{depletion_handling_enabled} = Val{false}(),
                                    only2d::Val{only_2d} = Val{false}(),
                                    is_weighting_potential::Val{_is_weighting_potential} = Val{false}(),
                                    use_nthreads::Int = Base.Threads.nthreads(),
                                    max_n_iterations::Int = -1,
                                    verbose::Bool = true
                                    )::T where {T, S, depletion_handling_enabled, only_2d, _is_weighting_potential}
    n_iterations::Int = 0
    cf::T = Inf
    cfs::Vector{T} = fill(cf, 10)
    cl::T = _is_weighting_potential ? convergence_limit : abs(convergence_limit * pssrb.bias_voltage) # to get relative change in respect to bias voltage
    # To disable automatically ProgressMeters in CI builds:
    is_logging(io) = isa(io, Base.TTY) == false || (get(ENV, "CI", nothing) == "true")
    if verbose prog = ProgressThresh(cl; dt = 0.1, desc = "Convergence: ", output = stderr, enabled = !is_logging(stderr)) end
    while cf > cl
        for i in 1:n_iterations_between_checks
            update!(pssrb, use_nthreads = use_nthreads, depletion_handling = depletion_handling, only2d = only2d, is_weighting_potential = is_weighting_potential)
        end
        cf = update_and_get_max_abs_diff!(pssrb, depletion_handling, only2d, is_weighting_potential, use_nthreads)
        @inbounds cfs[1:end-1] = cfs[2:end]
        @inbounds cfs[end] = cf
        slope::T = abs(mean(diff(cfs)))
        if verbose ProgressMeter.update!(prog, cf) end
        n_iterations += n_iterations_between_checks
        if slope < cl
            # @info "Slope is basically 0 -> Converged: $slope"
            cf = slope
        end
        if max_n_iterations > 0 && n_iterations > max_n_iterations
            # @show n_iterations_between_checks
            @info "Maximum number of iterations reached. (`n_iterations = $(n_iterations)`)"
            break
        end
    end
    if depletion_handling_enabled
        tmp_point_types = pssrb.point_types .& undepleted_bit
        @showprogress "Checking undepleted regions " for i in 1:10
            update!(pssrb, use_nthreads = use_nthreads, depletion_handling = depletion_handling, only2d = only2d, is_weighting_potential = is_weighting_potential)
            @inbounds for i in eachindex(pssrb.point_types)
                if (pssrb.point_types[i] & undepleted_bit == 0) && (tmp_point_types[i] > 0)
                    pssrb.point_types[i] += undepleted_bit
                elseif (pssrb.point_types[i] & undepleted_bit > 0) && (tmp_point_types[i] == 0)
                    tmp_point_types[i] += undepleted_bit
                end
            end
        end
    end

    if verbose ProgressMeter.finish!(prog) end
    return cf
end

# """
#     refine_scalar_potential(p::ScalarPotential{T}, max_diffs::NTuple{3, T}, minimum_distances::NTuple{3, T};
#         only2d::Val{only_2d} = Val(size(p.data, 2)==1)) where {T, only_2d}
#
# Refine any scalar potential `p`.
#
# 1. Extent the grid to be a closed grid in all dimensions.
# 2. Refine the axis of the grid based on `max_diffs` and `minimum_applied_potential`:
#    Insert N new ticks between to existing ticks such that the potential difference between each tick becomes
#    smaller than `max_diff[i]` (i -> dimension) but that the distances between the ticks stays larger than `minimum_distances[i]`.
# 3. Create the new data array for the refined grid and fill it by interpolation of the the initial (coarse) grid.
# """
function refine_scalar_potential(p::ScalarPotential{T}, max_diffs::NTuple{3, T}, minimum_distances::NTuple{3, T};
        only2d::Val{only_2d} = Val(size(p.data, 2)==1)) where {T, only_2d}
    closed_potential = _get_closed_potential(p)
    new_grid = _create_refined_grid(closed_potential, T.(max_diffs), T.(minimum_distances))
    new_data = Array{T, 3}(undef, size(new_grid))
    if only_2d
        int = interpolate_closed_potential(closed_potential, only2d)
        for i3 in axes(new_data, 3)
            x3 = new_grid.axes[3].ticks[i3]
            for i1 in axes(new_data, 1)
                x1 = new_grid.axes[1].ticks[i1]
                new_data[i1, 1, i3] = int(x1, x3)
            end
        end
    else
        int = interpolate_closed_potential(closed_potential, only2d)
        for i3 in axes(new_data, 3)
            x3 = new_grid.axes[3].ticks[i3]
            for i2 in axes(new_data, 2)
                x2 = new_grid.axes[2].ticks[i2]
                for i1 in axes(new_data, 1)
                    x1 = new_grid.axes[1].ticks[i1]
                    new_data[i1, i2, i3] = int(x1, x2, x3)
                end
            end
        end
    end
    return _convert_to_original_potential(p, new_data, new_grid)
end

function interpolate_closed_potential(p::ScalarPotential, ::Val{true}) where {T}
    interpolate((p.grid.axes[1], p.grid.axes[3]), p.data[:,1,:], Gridded(Linear()))
end
function interpolate_closed_potential(p::ScalarPotential, ::Val{false}) where {T}
    interpolate(p.grid.axes, p.data, Gridded(Linear()))
end

_get_closed_ticks(ticks::Vector{T}, int::ClosedInterval{T}) where {T} = ticks
_get_closed_ticks(ticks::Vector{T}, int::Interval{:closed, :open, T}) where {T} = vcat(ticks, int.right)
_get_closed_ticks(ticks::Vector{T}, int::Interval{:open, :closed, T}) where {T} = vcat(int.left, ticks)

_get_closed_values(values::Array{T,3}, dim::Int, int::ClosedInterval{T}) where {T} = values

function _get_closed_values(values::Array{T,3}, dim::Int, int::Interval{:closed, :open, T})::Array{T,3} where {T}
    if dim == 1
        cat(values, values[1:1,:,:], dims=dim)
    elseif dim == 2
        cat(values, values[:,1:1,:], dims=dim)
    else # dim == 3
        cat(values, values[:,:,1:1], dims=dim)
    end
end
function _get_closed_values(values::Array{T,3}, dim::Int, int::Interval{:open, :closed, T})::Array{T,3} where {T}
    if dim == 1
        cat(values[end:end,:,:], values, dims=dim)
    elseif dim == 2
        cat(values[:,end:end,:], values, dims=dim)
    else # dim == 3
        cat(values[:,:,end:end], values, dims=dim)
    end
end

function _get_closed_axis(ax::DiscreteAxis{T, BL, BR}) where {T, BL, BR}
    ticks = _get_closed_ticks(ax.ticks, ax.interval)
    return DiscreteAxis{T, BL, BR, ClosedInterval{T}}(ClosedInterval(ax.interval.left, ax.interval.right), ticks)
end

# """
#     _get_closed_potential(p::ScalarPotential{T,3,CS}) where {T, CS}
#
# Returns an closed Grid & Potential:
# E.g. if one of the axis is {:closed,:open} it will turn this into {:closed,:closed}
# and also extend the `data` field of the potential in the respective dimension and fill
# it with the respective values.
# """
function _get_closed_potential(p::ScalarPotential{T,3,CS}) where {T, CS}
    grid = p.grid
    closed_values = _get_closed_values(p.data, 1, grid.axes[1].interval)
    closed_values = _get_closed_values(closed_values, 2, grid.axes[2].interval)
    closed_values = _get_closed_values(closed_values, 3, grid.axes[3].interval)
    closed_axes = broadcast(i -> _get_closed_axis(grid.axes[i]), (1, 2, 3))
    AT = typeof(closed_axes)
    closed_grid = Grid{T, 3, CS, AT}(closed_axes)
    ScalarPotential(p, closed_values, closed_grid)
end

# """
#     _convert_to_original_potential(::Type{P}, data, grid) where {P, T, CS}
#
# Basically the counterpart to `_get_closed_potential`.
# """
function _convert_to_original_potential(p::ScalarPotential{T,3,CS,ATO}, data::Array{T, 3}, grid::Grid{T, 3, CS, AT}) where {T, CS, AT, ATO}
    ATs = get_axes_type(ATO)
    axs = broadcast(i -> _convert_closed_axis(ATs[i], grid.axes[i]), (1, 2, 3))
    new_data = _reduce_closed_values(ATs, data)
    new_grid = Grid{T,3,CS,ATO}(axs)
    ScalarPotential(p, new_data, new_grid)
end

_convert_closed_axis(::Type{DiscreteAxis{T, BL, BR, I}}, ax::DiscreteAxis{T, BL, BR, I}) where {T, BL, BR, I} = ax

function _convert_closed_axis(::Type{DiscreteAxis{T, BL, BR, Interval{:closed, :open, T}}}, ax::DiscreteAxis{T, BL, BR, ClosedInterval{T}}) where {T, BL, BR}
    DiscreteAxis{T, BL, BR, Interval{:closed, :open, T}}(Interval{:closed, :open, T}(ax.interval.left, ax.interval.right), ax.ticks[1:end-1])
end
function _convert_closed_axis(::Type{DiscreteAxis{T, BL, BR, Interval{:open, :closed, T}}}, ax::DiscreteAxis{T, BL, BR, ClosedInterval{T}}) where {T, BL, BR}
    DiscreteAxis{T, BL, BR, Interval{:open, :closed, T}}(Interval{:open, :closed, T}(ax.interval.left, ax.interval.right), ax.ticks[2:end])
end

function _reduce_closed_values(axTypes::Tuple, a::Array{T, 3}) where {T}
   inds = broadcast(i -> _get_ind_range_of_closed_axis(size(a, i), axTypes[i]), (1, 2, 3))
   a[inds...]
end

_get_ind_range_of_closed_axis(l::Int, ::Type{DiscreteAxis{T, BL, BR, ClosedInterval{T}}}) where {T, BL, BR} = 1:l
_get_ind_range_of_closed_axis(l::Int, ::Type{DiscreteAxis{T, BL, BR, Interval{:closed, :open, T}}}) where {T, BL, BR} = 1:(l-1)
_get_ind_range_of_closed_axis(l::Int, ::Type{DiscreteAxis{T, BL, BR, Interval{:open, :closed, T}}}) where {T, BL, BR} = 2:l

function _create_refined_grid(p::ScalarPotential{T,3}, max_diffs::NTuple{3, T}, minimum_distances::NTuple{3, T}) where {T}
    max_diffs = broadcast(md -> iszero(md) ? Inf : md, max_diffs)
    n_1 = floor.(Int, [maximum(abs.(p.data[i+1,:,:] .- p.data[i,:,:])) for i in 1:size(p.data, 1)-1] ./ max_diffs[1])
    n_2 = floor.(Int, [maximum(abs.(p.data[:,i+1,:] .- p.data[:,i,:])) for i in 1:size(p.data, 2)-1] ./ max_diffs[2])
    n_3 = floor.(Int, [maximum(abs.(p.data[:,:,i+1] .- p.data[:,:,i])) for i in 1:size(p.data, 3)-1] ./ max_diffs[3])
    ns = (n_1, n_2, n_3)
    widths = diff.((p.grid.axes[1].ticks, p.grid.axes[2].ticks, p.grid.axes[3].ticks))
    sub_widths = broadcast(ia -> [widths[ia][i] / (ns[ia][i]+1) for i in eachindex(ns[ia])], (1,2,3))
    for ia in 1:3
        for i in eachindex(ns[ia])
            while sub_widths[ia][i] < minimum_distances[ia] && ns[ia][i] > 0
                ns[ia][i] -= 1
                sub_widths[ia][i] = widths[ia][i] / (ns[ia][i]+1)
            end
        end
    end
    for i in 1:3 # always add an even number of ticks
        if isodd(sum(ns[i]))
            i_max_width = findmax(sub_widths[i])[2]
            ns[i][i_max_width] += 1
            sub_widths[i][i_max_width] = widths[i][i_max_width] / (ns[i][i_max_width]+1)
        end
    end
    new_axes = broadcast(i -> _refine_axis(p.grid.axes[i], ns[i], sub_widths[i]), (1, 2, 3))
    return typeof(p.grid)(new_axes)
end

function _refine_axis(ax::DiscreteAxis{T, <:Any, <:Any, ClosedInterval{T}}, ns::Vector{Int}, sub_widths::Vector{T}) where {T, I}
    @assert length(ns) == length(ax.ticks)-1 # for ClosedInterval axis
    ticks = Vector{T}(undef, length(ax.ticks) + sum(ns))
    i = 1
    for j in eachindex(ns)
        ticks[i] = ax.ticks[j]
        for k in 1:ns[j]
            i += 1
            ticks[i] = ax.ticks[j] + k * sub_widths[j]
        end
        i += 1
    end
    ticks[end] = ax.ticks[end]
    typeof(ax)(ax.interval, ticks)
end

_extend_refinement_limits(rl::Real) = (rl, rl, rl )
_extend_refinement_limits(rl::Tuple{<:Real,<:Real,<:Real}) = rl

ext: jl


text:
"""
Computes the path of consumption and debt for the previously described
complete markets model where exogenous income follows a linear
state space
"""
function complete_ss(β::AbstractFloat,
                     b0::Union{AbstractFloat, Array},
                     x0::Union{AbstractFloat, Array},
                     A::Union{AbstractFloat, Array},
                     C::Union{AbstractFloat, Array},
                     S_y::Union{AbstractFloat, Array},
                     T::Integer=12)

    # Create a linear state space for simulation purposes
    # This adds "b" as a state to the linear state space system
    # so that setting the seed places shocks in same place for
    # both the complete and incomplete markets economy
    # Atilde = vcat(hcat(A, zeros(size(A,1), 1)),
    #               zeros(1, size(A,2) + 1))
    # Ctilde = vcat(C, zeros(1, 1))
    # S_ytilde = hcat(S_y, zeros(1, 1))

    lss = LSS(A, C, S_y, mu_0=x0)

    # Add extra state to initial condition
    # x0 = hcat(x0, 0)

    # Compute the (I - β*A)^{-1}
    rm = inv(eye(size(A, 1)) - β * A)

    # Constant level of consumption
    cbar = (1 - β) * (S_y * rm * x0 - b0)
    c_hist = ones(T) * cbar[1]

    # Debt
    x_hist, y_hist = simulate(lss, T)
    b_hist = (S_y * rm * x_hist - cbar[1] / (1.0 - β))

    return c_hist, vec(b_hist), vec(y_hist), x_hist
end

N_simul = 150

# Define parameters
α, ρ1, ρ2 = 10.0, 0.9, 0.0
σ = 1.0
# N_simul = 1
# T = N_simul
A = [1.0 0.0 0.0;
     α    ρ1  ρ2;
     0.0 1.0 0.0]
C = [0.0, σ, 0.0]
S_y = [1.0 1.0 0.0]
β, b0 = 0.95, -10.0
x0 = [1.0, α / (1 - ρ1), α / (1 - ρ1)]

# Do simulation for complete markets
s = rand(1:10000)
srand(s)  # Seeds get set the same for both economies
out = complete_ss(β, b0, x0, A, C, S_y, 150)
c_hist_com, b_hist_com, y_hist_com, x_hist_com = out

fig, ax = subplots(1, 2, figsize=(15, 5))

# Consumption plots
ax[1][:set_title]("Cons and income")
ax[1][:plot](1:N_simul, c_hist_com, label="consumption")
ax[1][:plot](1:N_simul, y_hist_com, label="income",
            lw=2, alpha=0.6, ls="--")
ax[1][:legend]()
ax[1][:set_xlabel]("Periods")
ax[1][:set_ylim]([-5.0, 110])

# Debt plots
ax[2][:set_title]("Debt and income")
ax[2][:plot](1:N_simul, b_hist_com, label="debt")
ax[2][:plot](1:N_simul, y_hist_com, label="Income",
            lw=2, alpha=0.6, ls="--")
ax[2][:legend]()
ax[2][:axhline](0, color="k", ls="--")
ax[2][:set_xlabel]("Periods")


text:
module GESVfunctions

const dlmach_val = dlmach()

# source (LAPACK is a software package provided by Univ. of Tennessee):
# this is a modification for Julia language, for square matrices!
# http://www.netlib.org/lapack/explore-html/d8/d72/dgesv_8f_source.html
function dgesv!(n, A, B, ipiv)
    # n : range of the square matrix
    # A: dimension matrix (n,n), lds
    # IPIV is INTEGER array, dimension (n)
    # B: dimension matrix (n,n)
    # INFO is INTEGER see http://www.netlib.org/lapack/explore-html/d8/d72/dgesv_8f_source.html

    #* Compute the LU factorization of A.
    # CALL dgetrf( n, n, a, lda, ipiv, info )
    dgetrf!(A, n, ipiv, info)
    # A is mutated according to:
    # http://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga0019443faea08275ca60a734d0593e60.html#ga0019443faea08275ca60a734d0593e60
    # "A is DOUBLE PRECISION array, dimension (LDA,N)
    # On entry, the M-by-N matrix to be factored.
    # On exit, the factors L and U from the factorization
    # A = P*L*U; the unit diagonal elements of L are not stored.    "

    if info == 0
        #* Solve the system A*X = B, overwriting B with X.
        # * A and B are mutated!
        # CALL dgetrs( 'No transpose', n, nrhs, a, lda, ipiv, b, ldb, info )
        dgetrs!("No transpose", n, nrhs, A, lda, ipiv, B, ldb, info)
    end
end
# needs to define dgetrf, dgetrs for dgesv

# needs to define dgemm, dgetrf2, dlaswp, dtrsm, xerbla, ilaenv for dgetrf
# A is mutated
function dgetrf!(A, n, ipiv, info)
    one = 1.0

    # Determine the block size for this environment.
    # nb = ilaenv(1, "DGETRF", " ", n, n, -1, -1)
    # for the above case, it resulted nb=64:
    nb = 64 #
    if (nb <= 1) || (nb >= n)
        # Use unblocked code.
        dgetrf2(n, n, A, n, ipiv, info)
    else
        #  Use blocked code.
        for j in 1:nb:n # 1, 1+nb, 1+(2*nb), ..., n
            jb = min(n - j + 1, nb)

            # Factor diagonal and subdiagonal blocks and test for exact singularity.
            dgetrf2( n - j + 1, jb, A[j,j], n, ipiv[j], iinfo )
            # TODO: it returns iinfo, rearrange it

            #  Adjust INFO and the pivot indices.
            if (info == 0) && (iinfo > 0)
                info = iinfo + j - 1
            end

            for i in j: min(n, j + jb - 1)
                ipiv[i] += j - 1
            end

            # Apply interchanges to columns 1:J-1.
            dlaswp( j - 1, A, n, j, j + jb - 1, ipiv, 1 )

            if j + jb <= n
                # Apply interchanges to columns J+JB:N
                dlaswp(n - j - jb + 1, A[1, j + jb], n, j, j + jb - 1, ipiv, 1)

                # Compute block row of U.
                dtrsm("Left", "Lower", "No transpose", "Unit", jb, n - j - jb + 1, one, A[j,j], n, A[j, j + jb], n )

                if j + jb <= m
                    # Update trailing submatrix.
                    dgemm("No transpose", "No transpose", n - j - jb + 1, n - j - jb + 1, jb, -one, A[j + jb, j], lda, A[j, j + jb], n, one, A[j + jb, j + jb], n )
                end
            end
end

#TODO: define ilaenv function
#* I found that ilaenv just returns nb=64 for `ilaenv(1, "DGETRF", " ", n, n, -1, -1)` in dgetrf!(...)
# nb =   ilaenv(1,     "DGETRF", " ", n, n, -1, -1)
# function ilaenv( ISPEC, NAME, OPTS,  N1, N2, N3, N4 )
# end

# dgetrf2(n, m, X, p, ipiv, info)
function dgetrf2( m, n, A, lda, IPIV, info )
    one  = 1.0
    zero = 0.0
    #
    if m == 1
        #* Use unblocked code for one row case
        #* Just need to handle IPIV and INFO
        IPIV[iStart] = 1
        if A[iStart, iStart] == zero
            info = 1
        end
    elseif n == 1
        #* Use unblocked code for one column case

        #* Compute machine safe minimum
        # sfmin = dlamch('S')
        sfmin = dlmach_val # defined as global variable in the current module

        #* Find pivot and test for singularity
        # i = idamax( m, A[1, 1], 1 )
        i = finIndxMax(A, iStart, iStart + m - 1) # see Utils.jl

        IPIV[iStart] = i

        if  A[i, iStart] != zero
            #* Apply the interchange
            swap1(X, iStart, i) # see Utils.jl

            #* Compute elements 2:M of the column
            if abs( A[1, 1] ) >= sfmin
                dscal( m - 1, one / A[1, 1], A[2, 1], 1 )
            else
                for i in 1:(m - 1)
                    A[1 + i, 1] = A[1 + i, 1] / A[1, 1]
                end
            end
        else
            info = 1
        end
    else
        min_mn = min(m, n)
        #* Use recursive code
        n1 = min_mn / 2
        n2 = n - n1

        # *               [ A11 ]
        # *        Factor [ --- ]
        # *               [ A21 ]
        dgetrf2( m, n1, A, lda, IPIV, iinfo )

        if (info == 0) && (iinfo > 0)
            info = iinfo
        end

        # *                              [ A12 ]
        # *        Apply interchanges to [ --- ]
        # *                              [ A22 ]
        n1plus1  = n1 + 1
        dlaswp( n2, A[1, n1plus1], lda, 1, n1, IPIV, 1 )

        # *        Solve A12
        dtrsm( 'L', 'L', 'N', 'U', n1, n2, one, A, lda, A[1, n1plus1], lda )

        # *        Update A22
        dgetrf2( m - n1, n2, A[n1plus1, n1plus1], lda, IPIV[n1plus1], iinfo )

        # *        Factor A22
        dgetrf2( m - n1, n2, A[n1plus1, n1plus1], lda, IPIV[n1plus1], iinfo )

        #* Adjust INFO and the pivot indices
        if (info == 0) && (iinfo > 0)
            info = iinfo + n1
        end

        for i in n1plus1:min_mn
            IPIV[i] += n1
        end

        #* Apply interchanges to A21
        dlaswp( n1, A[1, 1], lda, n1plus1, min_mn, IPIV, 1 )
    end
end

# http://www.netlib.org/lapack/explore-html/d5/dd4/dlamch_8f_source.html#l00068
function dlmach() # for dlamch('S')
    # https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/
    #
    # The distance between 1.0 and the next larger representable floating-point value of Float64.
    # eps() : same as eps(Float64) : 2.220446049250313e-16
    # eps(0.0) : 5.0e-324
    # prevfloat(Inf)  : 1.7976931348623157e308
    # nextfloat(-Inf) : -1.7976931348623157e308

    zero = 0.0
    one  = 1.0
    #
    # eps_ = eps(zero) * 0.5 # it will return 0.0 in REPL julia
    eps_ = eps(zero) # it returns 5.0e-324
    #
    small = one / prevfloat(Inf)
    sfmin = eps(zero) # tiny(zero)
    if small >= sfmin
        #* Use SMALL plus a bit, to avoid the possibility of rounding
        #* causing overflow when computing  1/sfmin.
        sfmin = small * (one + eps_)
    end
    return sfmin
end

end # module

max_stars_repo_name: chinchay/gesvLapack

max_issues_repo_name: chinchay/gesvLapack

max_forks_repo_name: chinchay/gesvLapack


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct StaticReserve{T <: ReserveDirection} <: Reserve{T}
        name::String
        timeframe::Float64
        requirement::Float64
        ext::Dict{String, Any}
        internal::InfrastructureSystemsInternal
    end

Data Structure for a proportional reserve product for system simulations.

# Arguments
- `name::String`
- `timeframe::Float64`: the relative saturation timeframe, validation range: (0, nothing), action if invalid: error
- `requirement::Float64`: the static value of required reserves, validation range: (0, nothing), action if invalid: error
- `ext::Dict{String, Any}`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct StaticReserve{T <: ReserveDirection} <: Reserve{T}
    name::String
    "the relative saturation timeframe"
    timeframe::Float64
    "the static value of required reserves"
    requirement::Float64
    ext::Dict{String, Any}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function StaticReserve{T}(name, timeframe, requirement, ext=Dict{String, Any}(), ) where T <: ReserveDirection
    StaticReserve{T}(name, timeframe, requirement, ext, InfrastructureSystemsInternal(), )
end

function StaticReserve{T}(; name, timeframe, requirement, ext=Dict{String, Any}(), ) where T <: ReserveDirection
    StaticReserve{T}(name, timeframe, requirement, ext, )
end

# Constructor for demo purposes; non-functional.
function StaticReserve{T}(::Nothing) where T <: ReserveDirection
    StaticReserve{T}(;
        name="init",
        timeframe=0.0,
        requirement=0.0,
        ext=Dict{String, Any}(),
    )
end

InfrastructureSystems.get_name(value::StaticReserve) = value.name
"""Get StaticReserve timeframe."""
get_timeframe(value::StaticReserve) = value.timeframe
"""Get StaticReserve requirement."""
get_requirement(value::StaticReserve) = value.requirement
"""Get StaticReserve ext."""
get_ext(value::StaticReserve) = value.ext
"""Get StaticReserve internal."""
get_internal(value::StaticReserve) = value.internal


text:
z = zeros(length(x))
z[1] ~ D1()
x[1] ~ D2(z[1])

%1 = length(x)
z = zeros(%1)
%2 = D1()
%3 = tilde(... @varname(z[1]) ... %2)
%4 = setindex!(z, 1, %3)
%5 = getindex(z, 1)
%6 = D2(%5)
%7 = tilde(... @varname(x[1] ... %6)
%8 = setindex!(x, 1, %7)


text:
module PDFEstimators

    using
        AverageShiftedHistograms,
        CUDA,
        ComputationalResources,
        DataFrames,
        DiffEqFlux,
        DiffEqSensitivity,
        DifferentialEquations,
        Distances,
        Distributions,
        Flux,
        ForwardDiff,
        GalacticFlux,
        GalacticOptim,
        GalacticOptimJL,
        KernelDensity,
        LineSearches,
        MLJBase,
        MLJFlux,
        MLJModelInterface,
        ModelingToolkit,
        Optim,
        Parameters,
        ReverseDiff,
        SciMLBase,
        ScientificTypes,
        Zygote,
        LinearAlgebra,
        Statistics

    include("core.jl")
    include("actual.jl")
    include("kde.jl")
    include("ash.jl")
    include("mvn.jl")
    include("ffjord.jl")
    include("utils.jl")

    MLJBase.metadata_pkg.(
        [ActualModel, KDEModel, FFJORDModel, MvnModel],
        package_name="PDFEstimators",
        package_uuid="4d826980-ec0a-4f65-a989-f1052e211ebd",
        package_url="https://github.com/prbzrg/PDFEstimators.jl",
        is_pure_julia=true,
        package_license="MIT",
        is_wrapper=false,
    )

    MLJBase.metadata_model(
        ActualModel,
        input_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
        target_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
        output_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
        supports_weights=false,
        docstring="ActualModel",
        load_path="PDFEstimators.ActualModel",
    )
    MLJBase.metadata_model(
        KDEModel,
        input_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
        target_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
        output_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
        supports_weights=false,
        docstring="KDEModel",
        load_path="PDFEstimators.KDEModel",
    )
    MLJBase.metadata_model(
        FFJORDModel,
        input_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
        target_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
        output_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
        supports_weights=false,
        docstring="FFJORDModel",
        load_path="PDFEstimators.FFJORDModel",
    )
    MLJBase.metadata_model(
        MvnModel,
        input_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
        target_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
        output_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
        supports_weights=false,
        docstring="MvnModel",
        load_path="PDFEstimators.MvnModel",
    )
end


text:
using Distributed

# Now add 2 procs that can exec in parallel (obviously it depends on your CPU
# what you actually gain from this though)
addprocs(2)

# Ensure BlackBoxOptim loaded on all workers
@everywhere using BlackBoxOptim

# define the function to optimize on all workers. Parallel eval only gives a gain
# if function to optimize is slow. For this example we introduce a fake sleep
# to make it slow since the function is actually very quick to eval...
@everywhere function slow_rosenbrock(x)
  sleep(0.001) # Fake a slower func to be optimized...
  return BlackBoxOptim.rosenbrock(x)
end

# First run without any parallel procs used in eval
opt1 = bbsetup(slow_rosenbrock; Method=:xnes, SearchRange = (-5.0, 5.0),
               NumDimensions = 50, MaxFuncEvals = 5000)
el1 = @elapsed res1 = bboptimize(opt1)
t1 = round(el1, digits=3)

# When Workers= option is given, BlackBoxOptim enables parallel
# evaluation of fitness using the specified worker processes
opt2 = bbsetup(slow_rosenbrock; Method=:xnes, SearchRange = (-5.0, 5.0),
               NumDimensions = 50, MaxFuncEvals = 5000, Workers = workers())
el2 = @elapsed res2 = bboptimize(opt2)
t2 = round(el2, digits=3)

println("Time: serial = $(t1)s, parallel = $(t2)s")
if t2 < t1
  println("Speedup is $(round(t1/t2, digits=1))x")
else
  println("Slowdown is $(round(t2/t1, digits=1))x")
end


text:
using Documenter, FranklinContent
FranklinContent.franklincontent_hfuncs()

# generate index.md
#
@inline docitem(item) =  "\n```@docs\n$(item)\n```\n"

c = IOBuffer()
exps = Set(names(FranklinContent)[2:end])
# first exported
for sym in exps
    write(c, docitem(sym))
end

for sym in names(FranklinContent; all=true)
    sym ∈ exps && continue
    syms = string(sym)
    if startswith(syms, "hfun_") ||
        startswith(syms, "lx_")
        write(c, docitem(syms))
    end
end
write("src/index.md", String(take!(c)))
close(c)

# [sym for sym in names(@__MODULE__; all=true) if startswith(string(sym), "hfun_")]

makedocs(sitename="Documentation for the FranklinContent Plugin", pages=["API" => "index.md" ])


text:
module NetLogo

using JavaCall

const JDouble  = JavaObject{Symbol("java.lang.Double")}
const LogoList = JavaObject{Symbol("org.nlogo.core.LogoList")}
const Agent    = JavaObject{Symbol("org.nlogo.api.Agent")}

export JDouble, LogoList, Agent

function init(nlpath)
  cd(nlpath)
  if JavaCall.isloaded()
    warn("JVM already initialized.")
  else
    foreach(JavaCall.addOpts, ["-Xmx1024m", "-Dfile.encoding=UTF-8"])
    JavaCall.addClassPath(joinpath(nlpath, "*"))
    JavaCall.init()
  end
end

function companion{T}(t::Type{JavaCall.JavaObject{T}})
  const O = JavaObject{Symbol(string(T) * "\$")}
  jfield(O, "MODULE\$", O)
end

HeadlessWorkspace = @jimport org.nlogo.headless.HeadlessWorkspace

newWorkspace() = jcall(companion(HeadlessWorkspace), "newInstance", HeadlessWorkspace, ())
command(w, source) = jcall(w, "command", Void, (JString,), source)
report(w, source) = unpack(jcall(w, "report", JObject, (JString,), source))
open(w, modelname) = jcall(w, "open", Void, (JString,), modelname)
dispose(w) = jcall(w, "dispose", Void, ())

function withworkspace(f)
  w = newWorkspace()
  try
    f(w)
  finally
    dispose(w)
  end
end

function withmodel(f, modelname)
  withworkspace() do w
    open(w, modelname)
    f(w)
  end
end

function unpack(o::JavaObject)
  realclass = ccall(JavaCall.jnifunc.GetObjectClass, Ptr{Void},
    (Ptr{JavaCall.JNIEnv}, Ptr{Void}), JavaCall.penv, o.ptr)
  is(t) = JavaCall.isConvertible(t, realclass)
  to(t) = convert(t, o)
  if is(JDouble)
    jcall(to(JDouble), "doubleValue", jdouble, ())
  elseif is(JString)
    JavaCall.unsafe_string(to(JString))
  elseif is(LogoList)
    to(LogoList)
  elseif is(Agent)
    println("agent!")
    to(Agent)
  else
    o
  end
end

Base.length(lst::LogoList) = jcall(lst, "length", jint, ())
Base.getindex(lst::LogoList, i::Int) = unpack(jcall(lst, "get", JObject, (jint,), i - 1))
Base.getindex(lst::LogoList, range::Range) = [getindex(lst, i) for i = range]
Base.endof(lst::LogoList) = Base.length(lst)
Base.start(::LogoList) = 1
Base.next(lst::LogoList, state) = (getindex(lst, state), state + 1)
Base.done(lst::LogoList, state) = state > endof(lst)

type UndefinedVariableError <: Exception
  variable::String
end
UndefinedVariableError(jvariable::JString) = UndefinedVariableError(JavaCall.unsafe_string(jvariable))
Base.showerror(io::IO, e::UndefinedVariableError) = print(io, "Nothing named $(e.variable) has been defined.");

function indexofvariable(a::Agent, variable::String)
  const W = JavaObject{Symbol("org.nlogo.agent.World")}
  const A = JavaObject{Symbol("org.nlogo.agent.Agent")}
  agent = convert(A, a)
  world = jcall(agent, "world", W, ())
  jvariable = JString(uppercase(variable)) # TODO: uppercase conversion should be done on JVM side to be safe
  zerobasedindex = jcall(world, "indexOfVariable", jint, (A, JString), agent, jvariable)
  zerobasedindex == -1 ? throw(UndefinedVariableError(jvariable)) : zerobasedindex + 1
end

Base.length(a::Agent) = length(jcall(a, "variables", Array{JObject, 1}, ()))
Base.getindex(a::Agent, i::Int) = unpack(jcall(a, "getVariable", JObject, (jint,), i - 1))
Base.getindex(a::Agent, range::Range) = [getindex(a, i) for i = range]
Base.getindex(a::Agent, variable::String) = getindex(a, indexofvariable(a, variable))
Base.endof(a::Agent) = Base.length(a)
Base.start(::Agent) = 1
Base.next(a::Agent, state) = (getindex(a, state), state + 1)
Base.done(a::Agent, state) = state > endof(a)

end


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

mutable struct _ConstraintInfo{F,S}
    func::F
    set::S
    dual_start::Union{Nothing,Float64}
end

_ConstraintInfo(func, set) = _ConstraintInfo(func, set, nothing)

"""
    Optimizer()

Create a new Ipopt optimizer.
"""
mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::Union{Nothing,IpoptProblem}
    invalid_model::Bool
    variables::MOI.Utilities.VariablesContainer{Float64}
    variable_primal_start::Vector{Union{Nothing,Float64}}
    variable_lower_start::Vector{Union{Nothing,Float64}}
    variable_upper_start::Vector{Union{Nothing,Float64}}
    nlp_data::MOI.NLPBlockData
    sense::MOI.OptimizationSense
    objective::Union{
        Nothing,
        MOI.VariableIndex,
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    }
    linear_le_constraints::Vector{
        _ConstraintInfo{
            MOI.ScalarAffineFunction{Float64},
            MOI.LessThan{Float64},
        },
    }
    linear_ge_constraints::Vector{
        _ConstraintInfo{
            MOI.ScalarAffineFunction{Float64},
            MOI.GreaterThan{Float64},
        },
    }
    linear_eq_constraints::Vector{
        _ConstraintInfo{MOI.ScalarAffineFunction{Float64},MOI.EqualTo{Float64}},
    }
    quadratic_le_constraints::Vector{
        _ConstraintInfo{
            MOI.ScalarQuadraticFunction{Float64},
            MOI.LessThan{Float64},
        },
    }
    quadratic_ge_constraints::Vector{
        _ConstraintInfo{
            MOI.ScalarQuadraticFunction{Float64},
            MOI.GreaterThan{Float64},
        },
    }
    quadratic_eq_constraints::Vector{
        _ConstraintInfo{
            MOI.ScalarQuadraticFunction{Float64},
            MOI.EqualTo{Float64},
        },
    }
    nlp_dual_start::Union{Nothing,Vector{Float64}}
    silent::Bool
    options::Dict{String,Any}
    solve_time::Float64
    callback::Union{Nothing,Function}

    function Optimizer()
        return new(
            nothing,
            false,
            MOI.Utilities.VariablesContainer{Float64}(),
            Union{Nothing,Float64}[],
            Union{Nothing,Float64}[],
            Union{Nothing,Float64}[],
            MOI.NLPBlockData([], _EmptyNLPEvaluator(), false),
            MOI.FEASIBILITY_SENSE,
            nothing,
            _ConstraintInfo{
                MOI.ScalarAffineFunction{Float64},
                MOI.LessThan{Float64},
            }[],
            _ConstraintInfo{
                MOI.ScalarAffineFunction{Float64},
                MOI.GreaterThan{Float64},
            }[],
            _ConstraintInfo{
                MOI.ScalarAffineFunction{Float64},
                MOI.EqualTo{Float64},
            }[],
            _ConstraintInfo{
                MOI.ScalarQuadraticFunction{Float64},
                MOI.LessThan{Float64},
            }[],
            _ConstraintInfo{
                MOI.ScalarQuadraticFunction{Float64},
                MOI.GreaterThan{Float64},
            }[],
            _ConstraintInfo{
                MOI.ScalarQuadraticFunction{Float64},
                MOI.EqualTo{Float64},
            }[],
            nothing,
            false,
            Dict{String,Any}(),
            NaN,
            nothing,
        )
    end
end

MOI.get(::Optimizer, ::MOI.SolverVersion) = "3.14.4"

### _EmptyNLPEvaluator

struct _EmptyNLPEvaluator <: MOI.AbstractNLPEvaluator end

MOI.features_available(::_EmptyNLPEvaluator) = [:Grad, :Jac, :Hess]
MOI.initialize(::_EmptyNLPEvaluator, ::Any) = nothing
MOI.eval_constraint(::_EmptyNLPEvaluator, g, x) = nothing
MOI.jacobian_structure(::_EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
MOI.hessian_lagrangian_structure(::_EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
MOI.eval_constraint_jacobian(::_EmptyNLPEvaluator, J, x) = nothing
MOI.eval_hessian_lagrangian(::_EmptyNLPEvaluator, H, x, σ, μ) = nothing

function MOI.empty!(model::Optimizer)
    model.inner = nothing
    model.invalid_model = false
    MOI.empty!(model.variables)
    empty!(model.variable_primal_start)
    empty!(model.variable_lower_start)
    empty!(model.variable_upper_start)
    model.nlp_data = MOI.NLPBlockData([], _EmptyNLPEvaluator(), false)
    model.sense = MOI.FEASIBILITY_SENSE
    model.objective = nothing
    empty!(model.linear_le_constraints)
    empty!(model.linear_ge_constraints)
    empty!(model.linear_eq_constraints)
    empty!(model.quadratic_le_constraints)
    empty!(model.quadratic_ge_constraints)
    empty!(model.quadratic_eq_constraints)
    model.nlp_dual_start = nothing
    return
end

function MOI.is_empty(model::Optimizer)
    return MOI.is_empty(model.variables) &&
           isempty(model.variable_primal_start) &&
           isempty(model.variable_lower_start) &&
           isempty(model.variable_upper_start) &&
           model.nlp_data.evaluator isa _EmptyNLPEvaluator &&
           model.sense == MOI.FEASIBILITY_SENSE &&
           isempty(model.linear_le_constraints) &&
           isempty(model.linear_ge_constraints) &&
           isempty(model.linear_eq_constraints) &&
           isempty(model.quadratic_le_constraints) &&
           isempty(model.quadratic_ge_constraints) &&
           isempty(model.quadratic_eq_constraints)
end

MOI.supports_incremental_interface(::Optimizer) = true

function MOI.copy_to(model::Optimizer, src::MOI.ModelLike)
    return MOI.Utilities.default_copy_to(model, src)
end

MOI.get(::Optimizer, ::MOI.SolverName) = "Ipopt"

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{
        <:Union{
            MOI.VariableIndex,
            MOI.ScalarAffineFunction{Float64},
            MOI.ScalarQuadraticFunction{Float64},
        },
    },
    ::Type{
        <:Union{
            MOI.LessThan{Float64},
            MOI.GreaterThan{Float64},
            MOI.EqualTo{Float64},
        },
    },
)
    return true
end

function MOI.get(model::Optimizer, ::MOI.ListOfConstraintTypesPresent)
    ret = MOI.get(model.variables, MOI.ListOfConstraintTypesPresent())
    constraints = Set{Tuple{Type,Type}}()
    for F in (
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    )
        for S in (
            MOI.LessThan{Float64},
            MOI.GreaterThan{Float64},
            MOI.EqualTo{Float64},
        )
            if !isempty(_constraints(model, F, S))
                push!(constraints, (F, S))
            end
        end
    end
    return append!(ret, collect(constraints))
end

### MOI.Silent

MOI.supports(::Optimizer, ::MOI.Silent) = true

function MOI.set(model::Optimizer, ::MOI.Silent, value)
    model.silent = value
    return
end

MOI.get(model::Optimizer, ::MOI.Silent) = model.silent

### MOI.TimeLimitSec

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true

function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value::Real)
    MOI.set(model, MOI.RawOptimizerAttribute("max_cpu_time"), Float64(value))
    return
end

function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, ::Nothing)
    delete!(model.options, "max_cpu_time")
    return
end

function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return get(model.options, "max_cpu_time", nothing)
end

### MOI.RawOptimizerAttribute

MOI.supports(::Optimizer, ::MOI.RawOptimizerAttribute) = true

function MOI.set(model::Optimizer, p::MOI.RawOptimizerAttribute, value)
    model.options[p.name] = value
    return
end

function MOI.get(model::Optimizer, p::MOI.RawOptimizerAttribute)
    if !haskey(model.options, p.name)
        error("RawParameter with name $(p.name) is not set.")
    end
    return model.options[p.name]
end

### Variables

"""
    column(x::MOI.VariableIndex)

Return the column associated with a variable.
"""
column(x::MOI.VariableIndex) = x.value

function MOI.add_variable(model::Optimizer)
    push!(model.variable_primal_start, nothing)
    push!(model.variable_lower_start, nothing)
    push!(model.variable_upper_start, nothing)
    return MOI.add_variable(model.variables)
end

function MOI.is_valid(model::Optimizer, x::MOI.VariableIndex)
    return MOI.is_valid(model.variables, x)
end

function MOI.get(
    model::Optimizer,
    attr::Union{MOI.NumberOfVariables,MOI.ListOfVariableIndices},
)
    return MOI.get(model.variables, attr)
end

function MOI.is_valid(
    model::Optimizer,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {S<:Union{MOI.LessThan,MOI.GreaterThan,MOI.EqualTo}}
    return MOI.is_valid(model.variables, ci)
end

function MOI.get(
    model::Optimizer,
    attr::Union{
        MOI.NumberOfConstraints{MOI.VariableIndex,S},
        MOI.ListOfConstraintIndices{MOI.VariableIndex,S},
    },
) where {S<:Union{MOI.LessThan,MOI.GreaterThan,MOI.EqualTo}}
    return MOI.get(model.variables, attr)
end

function MOI.get(
    model::Optimizer,
    attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet},
    c::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {S<:Union{MOI.LessThan,MOI.GreaterThan,MOI.EqualTo}}
    return MOI.get(model.variables, attr, c)
end

function MOI.add_constraint(
    model::Optimizer,
    x::MOI.VariableIndex,
    set::Union{
        MOI.LessThan{Float64},
        MOI.GreaterThan{Float64},
        MOI.EqualTo{Float64},
    },
)
    return MOI.add_constraint(model.variables, x, set)
end

function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintSet,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
    set::S,
) where {S<:Union{MOI.LessThan,MOI.GreaterThan,MOI.EqualTo}}
    MOI.set(model.variables, MOI.ConstraintSet(), ci, set)
    return
end

function MOI.delete(
    model::Optimizer,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,S},
) where {S<:Union{MOI.LessThan,MOI.GreaterThan,MOI.EqualTo}}
    MOI.delete(model.variables, ci)
    return
end

### ScalarAffineFunction and ScalarQuadraticFunction constraints

function MOI.is_valid(
    model::Optimizer,
    ci::MOI.ConstraintIndex{F,S},
) where {
    F<:Union{
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    },
    S<:Union{MOI.LessThan,MOI.GreaterThan,MOI.EqualTo},
}
    return 1 <= ci.value <= length(_constraints(model, F, S))
end

function _constraints(
    model::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{MOI.LessThan{Float64}},
)
    return model.linear_le_constraints
end

function _constraints(
    model::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{MOI.GreaterThan{Float64}},
)
    return model.linear_ge_constraints
end

function _constraints(
    model::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{MOI.EqualTo{Float64}},
)
    return model.linear_eq_constraints
end

function _constraints(
    model::Optimizer,
    ::Type{MOI.ScalarQuadraticFunction{Float64}},
    ::Type{MOI.LessThan{Float64}},
)
    return model.quadratic_le_constraints
end

function _constraints(
    model::Optimizer,
    ::Type{MOI.ScalarQuadraticFunction{Float64}},
    ::Type{MOI.GreaterThan{Float64}},
)
    return model.quadratic_ge_constraints
end

function _constraints(
    model::Optimizer,
    ::Type{MOI.ScalarQuadraticFunction{Float64}},
    ::Type{MOI.EqualTo{Float64}},
)
    return model.quadratic_eq_constraints
end

function _check_inbounds(model::Optimizer, var::MOI.VariableIndex)
    MOI.throw_if_not_valid(model, var)
    return
end

function _check_inbounds(model::Optimizer, aff::MOI.ScalarAffineFunction)
    for term in aff.terms
        MOI.throw_if_not_valid(model, term.variable)
    end
    return
end

function _check_inbounds(model::Optimizer, quad::MOI.ScalarQuadraticFunction)
    for term in quad.affine_terms
        MOI.throw_if_not_valid(model, term.variable)
    end
    for term in quad.quadratic_terms
        MOI.throw_if_not_valid(model, term.variable_1)
        MOI.throw_if_not_valid(model, term.variable_2)
    end
    return
end

function MOI.add_constraint(
    model::Optimizer,
    func::F,
    set::S,
) where {
    F<:Union{
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    },
    S<:MOI.AbstractScalarSet,
}
    _check_inbounds(model, func)
    constraints = _constraints(model, F, S)
    push!(constraints, _ConstraintInfo(func, set))
    return MOI.ConstraintIndex{F,S}(length(constraints))
end

function MOI.get(
    model::Optimizer,
    ::MOI.NumberOfConstraints{F,S},
) where {
    F<:Union{
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    },
    S,
}
    return length(_constraints(model, F, S))
end

function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{F,S},
) where {
    F<:Union{
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    },
    S,
}
    return MOI.ConstraintIndex{F,S}[
        MOI.ConstraintIndex{F,S}(i) for
        i in eachindex(_constraints(model, F, S))
    ]
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{F,S},
) where {
    F<:Union{
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    },
    S,
}
    return _constraints(model, F, S)[c.value].func
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{F,S},
) where {
    F<:Union{
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    },
    S,
}
    return _constraints(model, F, S)[c.value].set
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ConstraintDualStart,
    ::Type{MOI.ConstraintIndex{F,S}},
) where {
    F<:Union{
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    },
    S,
}
    return true
end

function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintDualStart,
    ci::MOI.ConstraintIndex{F,S},
    value::Union{Real,Nothing},
) where {
    F<:Union{
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    },
    S,
}
    MOI.throw_if_not_valid(model, ci)
    constraints = _constraints(model, F, S)
    constraints[ci.value].dual_start = value
    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintDualStart,
    ci::MOI.ConstraintIndex{F,S},
) where {
    F<:Union{
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    },
    S,
}
    MOI.throw_if_not_valid(model, ci)
    constraints = _constraints(model, F, S)
    return constraints[ci.value].dual_start
end

### MOI.VariablePrimalStart

function MOI.supports(
    ::Optimizer,
    ::MOI.VariablePrimalStart,
    ::Type{MOI.VariableIndex},
)
    return true
end

function MOI.set(
    model::Optimizer,
    ::MOI.VariablePrimalStart,
    vi::MOI.VariableIndex,
    value::Union{Real,Nothing},
)
    MOI.throw_if_not_valid(model, vi)
    model.variable_primal_start[column(vi)] = value
    return
end

### MOI.ConstraintDualStart

_dual_start(::Optimizer, ::Nothing, ::Int = 1) = 0.0

function _dual_start(model::Optimizer, value::Real, scale::Int = 1)
    return _dual_multiplier(model) * value * scale
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ConstraintDualStart,
    ::Type{
        MOI.ConstraintIndex{
            MOI.VariableIndex,
            <:Union{MOI.GreaterThan,MOI.LessThan,MOI.EqualTo},
        },
    },
)
    return true
end

function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintDualStart,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}},
    value::Union{Real,Nothing},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_lower_start[ci.value] = value
    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintDualStart,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}},
)
    MOI.throw_if_not_valid(model, ci)
    return model.variable_lower_start[ci.value]
end

function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintDualStart,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}},
    value::Union{Real,Nothing},
)
    MOI.throw_if_not_valid(model, ci)
    model.variable_upper_start[ci.value] = value
    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintDualStart,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}},
)
    MOI.throw_if_not_valid(model, ci)
    return model.variable_upper_start[ci.value]
end

function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintDualStart,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}},
    value::Union{Real,Nothing},
)
    MOI.throw_if_not_valid(model, ci)
    if value === nothing
        model.variable_lower_start[ci.value] = nothing
        model.variable_upper_start[ci.value] = nothing
    elseif value >= 0.0
        model.variable_lower_start[ci.value] = value
        model.variable_upper_start[ci.value] = 0.0
    else
        model.variable_lower_start[ci.value] = 0.0
        model.variable_upper_start[ci.value] = value
    end
    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintDualStart,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}},
)
    MOI.throw_if_not_valid(model, ci)
    l = model.variable_lower_start[ci.value]
    u = model.variable_upper_start[ci.value]
    return (l === u === nothing) ? nothing : (l + u)
end

### MOI.NLPBlockDualStart

MOI.supports(::Optimizer, ::MOI.NLPBlockDualStart) = true

function MOI.set(
    model::Optimizer,
    ::MOI.NLPBlockDualStart,
    values::Union{Nothing,Vector},
)
    model.nlp_dual_start = values
    return
end

MOI.get(model::Optimizer, ::MOI.NLPBlockDualStart) = model.nlp_dual_start

### MOI.NLPBlock

MOI.supports(::Optimizer, ::MOI.NLPBlock) = true

function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    model.nlp_data = nlp_data
    return
end

### ObjectiveSense

MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

function MOI.set(
    model::Optimizer,
    ::MOI.ObjectiveSense,
    sense::MOI.OptimizationSense,
)
    model.sense = sense
    return
end

MOI.get(model::Optimizer, ::MOI.ObjectiveSense) = model.sense

### ObjectiveFunction

MOI.get(model::Optimizer, ::MOI.ObjectiveFunctionType) = typeof(model.objective)

function MOI.get(model::Optimizer, ::MOI.ObjectiveFunction{F}) where {F}
    return convert(F, model.objective)::F
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ObjectiveFunction{
        <:Union{
            MOI.VariableIndex,
            MOI.ScalarAffineFunction{Float64},
            MOI.ScalarQuadraticFunction{Float64},
        },
    },
)
    return true
end

function MOI.set(
    model::Optimizer,
    ::MOI.ObjectiveFunction{F},
    func::F,
) where {
    F<:Union{
        MOI.VariableIndex,
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    },
}
    _check_inbounds(model, func)
    model.objective = func
    return
end

### Ipopt callback functions
### In setting up the data for Ipopt, we order the constraints as follows:
### - linear_le_constraints
### - linear_ge_constraints
### - linear_eq_constraints
### - quadratic_le_constraints
### - quadratic_ge_constraints
### - quadratic_eq_constraints
### - nonlinear constraints from nlp_data

const _CONSTRAINT_ORDERING = (
    :linear_le_constraints,
    :linear_ge_constraints,
    :linear_eq_constraints,
    :quadratic_le_constraints,
    :quadratic_ge_constraints,
    :quadratic_eq_constraints,
)

function _offset(
    ::Optimizer,
    ::Type{<:MOI.ScalarAffineFunction},
    ::Type{<:MOI.LessThan},
)
    return 0
end

function _offset(
    model::Optimizer,
    ::Type{<:MOI.ScalarAffineFunction},
    ::Type{<:MOI.GreaterThan},
)
    return length(model.linear_le_constraints)
end

function _offset(
    model::Optimizer,
    F::Type{<:MOI.ScalarAffineFunction},
    ::Type{<:MOI.EqualTo},
)
    return _offset(model, F, MOI.GreaterThan{Float64}) +
           length(model.linear_ge_constraints)
end

function _offset(
    model::Optimizer,
    ::Type{<:MOI.ScalarQuadraticFunction},
    ::Type{<:MOI.LessThan},
)
    x = _offset(model, MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})
    return x + length(model.linear_eq_constraints)
end

function _offset(
    model::Optimizer,
    F::Type{<:MOI.ScalarQuadraticFunction},
    ::Type{<:MOI.GreaterThan},
)
    return _offset(model, F, MOI.LessThan{Float64}) +
           length(model.quadratic_le_constraints)
end

function _offset(
    model::Optimizer,
    F::Type{<:MOI.ScalarQuadraticFunction},
    ::Type{<:MOI.EqualTo},
)
    return _offset(model, F, MOI.GreaterThan{Float64}) +
           length(model.quadratic_ge_constraints)
end

function _nlp_constraint_offset(model::Optimizer)
    x = _offset(
        model,
        MOI.ScalarQuadraticFunction{Float64},
        MOI.EqualTo{Float64},
    )
    return x + length(model.quadratic_eq_constraints)
end

_eval_function(::Nothing, ::Any) = 0.0

_eval_function(f, x) = MOI.Utilities.eval_variables(xi -> x[xi.value], f)

### Eval_F_CB

function _eval_objective(model::Optimizer, x)
    if model.nlp_data.has_objective
        return MOI.eval_objective(model.nlp_data.evaluator, x)
    end
    return _eval_function(model.objective, x)
end

### Eval_Grad_F_CB

_fill_gradient(::Any, ::Any, ::Nothing) = nothing

function _fill_gradient(grad, ::Vector, f::MOI.VariableIndex)
    grad[f.value] = 1.0
    return
end

function _fill_gradient(grad, ::Vector, f::MOI.ScalarAffineFunction{Float64})
    for term in f.terms
        grad[term.variable.value] += term.coefficient
    end
    return
end

function _fill_gradient(
    grad,
    x::Vector,
    quad::MOI.ScalarQuadraticFunction{Float64},
)
    for term in quad.affine_terms
        grad[term.variable.value] += term.coefficient
    end
    for term in quad.quadratic_terms
        row_idx = term.variable_1
        col_idx = term.variable_2
        if row_idx == col_idx
            grad[row_idx.value] += term.coefficient * x[row_idx.value]
        else
            grad[row_idx.value] += term.coefficient * x[col_idx.value]
            grad[col_idx.value] += term.coefficient * x[row_idx.value]
        end
    end
    return
end

function _eval_objective_gradient(model::Optimizer, grad, x)
    if model.nlp_data.has_objective
        MOI.eval_objective_gradient(model.nlp_data.evaluator, grad, x)
    else
        fill!(grad, 0.0)
        _fill_gradient(grad, x, model.objective)
    end
    return
end

### Eval_G_CB

function _eval_constraint(model::Optimizer, g, x)
    row = 1
    for key in _CONSTRAINT_ORDERING
        for info in getfield(model, key)
            g[row] = _eval_function(info.func, x)
            row += 1
        end
    end
    nlp_g = view(g, row:length(g))
    MOI.eval_constraint(model.nlp_data.evaluator, nlp_g, x)
    return
end

### Eval_Jac_G_CB

function _append_to_jacobian_sparsity(J, f::MOI.ScalarAffineFunction, row)
    for term in f.terms
        push!(J, (row, term.variable.value))
    end
    return
end

function _append_to_jacobian_sparsity(J, f::MOI.ScalarQuadraticFunction, row)
    for term in f.affine_terms
        push!(J, (row, term.variable.value))
    end
    for term in f.quadratic_terms
        row_idx = term.variable_1
        col_idx = term.variable_2
        if row_idx == col_idx
            push!(J, (row, row_idx.value))
        else
            push!(J, (row, row_idx.value))
            push!(J, (row, col_idx.value))
        end
    end
    return
end

function _jacobian_structure(model::Optimizer)
    J = Tuple{Int64,Int64}[]
    row = 1
    for key in _CONSTRAINT_ORDERING
        for info in getfield(model, key)
            _append_to_jacobian_sparsity(J, info.func, row)
            row += 1
        end
    end
    if length(model.nlp_data.constraint_bounds) > 0
        for (nlp_row, col) in MOI.jacobian_structure(model.nlp_data.evaluator)
            push!(J, (nlp_row + row - 1, col))
        end
    end
    return J
end

function _fill_constraint_jacobian(
    values,
    offset,
    ::Vector,
    f::MOI.ScalarAffineFunction,
)
    num_coefficients = length(f.terms)
    for i in 1:num_coefficients
        values[offset+i] = f.terms[i].coefficient
    end
    return num_coefficients
end

function _fill_constraint_jacobian(
    values,
    offset,
    x,
    f::MOI.ScalarQuadraticFunction,
)
    nterms = 0
    for term in f.affine_terms
        nterms += 1
        values[offset+nterms] = term.coefficient
    end
    for term in f.quadratic_terms
        row_idx = term.variable_1
        col_idx = term.variable_2
        if row_idx == col_idx
            nterms += 1
            values[offset+nterms] = term.coefficient * x[col_idx.value]
        else
            # Note that the order matches the Jacobian sparsity pattern.
            nterms += 2
            values[offset+nterms-1] = term.coefficient * x[col_idx.value]
            values[offset+nterms] = term.coefficient * x[row_idx.value]
        end
    end
    return nterms
end

function _eval_constraint_jacobian(model::Optimizer, values, x)
    offset = 0
    for key in _CONSTRAINT_ORDERING
        for info in getfield(model, key)
            offset += _fill_constraint_jacobian(values, offset, x, info.func)
        end
    end
    nlp_values = view(values, (1+offset):length(values))
    MOI.eval_constraint_jacobian(model.nlp_data.evaluator, nlp_values, x)
    return
end

### Eval_H_CB

_append_to_hessian_sparsity(::Any, ::Any) = nothing

function _append_to_hessian_sparsity(H, f::MOI.ScalarQuadraticFunction)
    for term in f.quadratic_terms
        push!(H, (term.variable_1.value, term.variable_2.value))
    end
    return
end

function _append_hessian_lagrangian_structure(H, model::Optimizer)
    if !model.nlp_data.has_objective
        _append_to_hessian_sparsity(H, model.objective)
    end
    for info in model.quadratic_le_constraints
        _append_to_hessian_sparsity(H, info.func)
    end
    for info in model.quadratic_ge_constraints
        _append_to_hessian_sparsity(H, info.func)
    end
    for info in model.quadratic_eq_constraints
        _append_to_hessian_sparsity(H, info.func)
    end
    append!(H, MOI.hessian_lagrangian_structure(model.nlp_data.evaluator))
    return
end

_fill_hessian_lagrangian(::Any, ::Any, ::Any, ::Any) = 0

function _fill_hessian_lagrangian(H, offset, λ, f::MOI.ScalarQuadraticFunction)
    for term in f.quadratic_terms
        H[offset+1] = λ * term.coefficient
        offset += 1
    end
    return length(f.quadratic_terms)
end

function _eval_hessian_lagrangian(
    ::Type{S},
    model::Optimizer,
    H,
    μ,
    offset,
) where {S}
    F = MOI.ScalarQuadraticFunction{Float64}
    offset_start = _offset(model, F, S)
    for (i, info) in enumerate(_constraints(model, F, S))
        offset +=
            _fill_hessian_lagrangian(H, offset, μ[offset_start+i], info.func)
    end
    return offset
end

function _eval_hessian_lagrangian(model::Optimizer, H, x, σ, μ)
    offset = 0
    if !model.nlp_data.has_objective
        offset += _fill_hessian_lagrangian(H, 0, σ, model.objective)
    end
    # Handles any quadratic constraints that are present. The order matters.
    offset =
        _eval_hessian_lagrangian(MOI.LessThan{Float64}, model, H, μ, offset)
    offset =
        _eval_hessian_lagrangian(MOI.GreaterThan{Float64}, model, H, μ, offset)
    offset = _eval_hessian_lagrangian(MOI.EqualTo{Float64}, model, H, μ, offset)
    # Handles the Hessian in the nonlinear block
    MOI.eval_hessian_lagrangian(
        model.nlp_data.evaluator,
        view(H, 1+offset:length(H)),
        x,
        σ,
        view(μ, 1+_nlp_constraint_offset(model):length(μ)),
    )
    return
end

### MOI.optimize!

_bounds(s::MOI.LessThan) = (-Inf, s.upper)
_bounds(s::MOI.GreaterThan) = (s.lower, Inf)
_bounds(s::MOI.EqualTo) = (s.value, s.value)

function MOI.optimize!(model::Optimizer)
    # TODO: Reuse model.inner for incremental solves if possible.
    num_quadratic_constraints =
        length(model.quadratic_le_constraints) +
        length(model.quadratic_ge_constraints) +
        length(model.quadratic_eq_constraints)
    num_nlp_constraints = length(model.nlp_data.constraint_bounds)
    has_hessian = :Hess in MOI.features_available(model.nlp_data.evaluator)
    init_feat = [:Grad]
    if has_hessian
        push!(init_feat, :Hess)
    end
    if num_nlp_constraints > 0
        push!(init_feat, :Jac)
    end
    MOI.initialize(model.nlp_data.evaluator, init_feat)
    jacobian_sparsity = _jacobian_structure(model)
    hessian_sparsity = Tuple{Int,Int}[]
    if has_hessian
        _append_hessian_lagrangian_structure(hessian_sparsity, model)
    end
    function eval_f_cb(x)
        # TODO(odow): FEASIBILITY_SENSE could produce confusing solver output if
        # a nonzero objective is set.
        if model.sense == MOI.FEASIBILITY_SENSE
            return 0.0
        end
        return _eval_objective(model, x)
    end
    function eval_grad_f_cb(x, grad_f)
        if model.sense == MOI.FEASIBILITY_SENSE
            grad_f .= zero(eltype(grad_f))
        else
            _eval_objective_gradient(model, grad_f, x)
        end
        return
    end
    eval_g_cb(x, g) = _eval_constraint(model, g, x)
    function eval_jac_g_cb(x, rows, cols, values)
        if values === nothing
            for i in 1:length(jacobian_sparsity)
                rows[i], cols[i] = jacobian_sparsity[i]
            end
        else
            _eval_constraint_jacobian(model, values, x)
        end
        return
    end
    function eval_h_cb(x, rows, cols, obj_factor, lambda, values)
        if values === nothing
            for i in 1:length(hessian_sparsity)
                rows[i], cols[i] = hessian_sparsity[i]
            end
        else
            _eval_hessian_lagrangian(model, values, x, obj_factor, lambda)
        end
        return
    end
    g_L, g_U = Float64[], Float64[]
    for key in _CONSTRAINT_ORDERING
        for info in getfield(model, key)
            l, u = _bounds(info.set)
            push!(g_L, l)
            push!(g_U, u)
        end
    end
    for bound in model.nlp_data.constraint_bounds
        push!(g_L, bound.lower)
        push!(g_U, bound.upper)
    end
    start_time = time()
    if length(model.variables.lower) == 0
        # Don't attempt to create a problem because Ipopt will error.
        model.invalid_model = true
        return
    end
    model.inner = CreateIpoptProblem(
        length(model.variables.lower),
        model.variables.lower,
        model.variables.upper,
        length(g_L),
        g_L,
        g_U,
        length(jacobian_sparsity),
        length(hessian_sparsity),
        eval_f_cb,
        eval_g_cb,
        eval_grad_f_cb,
        eval_jac_g_cb,
        has_hessian ? eval_h_cb : nothing,
    )
    if model.sense == MOI.MIN_SENSE
        AddIpoptNumOption(model.inner, "obj_scaling_factor", 1.0)
    elseif model.sense == MOI.MAX_SENSE
        AddIpoptNumOption(model.inner, "obj_scaling_factor", -1.0)
    end
    # Ipopt crashes by default if NaN/Inf values are returned from the
    # evaluation callbacks. This option tells Ipopt to explicitly check for them
    # and return Invalid_Number_Detected instead. This setting may result in a
    # minor performance loss and can be overwritten by specifying
    # check_derivatives_for_naninf="no".
    AddIpoptStrOption(model.inner, "check_derivatives_for_naninf", "yes")
    if !has_hessian
        AddIpoptStrOption(
            model.inner,
            "hessian_approximation",
            "limited-memory",
        )
    end
    if num_nlp_constraints == 0 && num_quadratic_constraints == 0
        AddIpoptStrOption(model.inner, "jac_c_constant", "yes")
        AddIpoptStrOption(model.inner, "jac_d_constant", "yes")
        if !model.nlp_data.has_objective
            # We turn on this option if all constraints are linear and the
            # objective is linear or quadratic. From the documentation, it's
            # unclear if it may also apply if the constraints are at most
            # quadratic.
            AddIpoptStrOption(model.inner, "hessian_constant", "yes")
        end
    end
    if model.silent
        AddIpoptIntOption(model.inner, "print_level", 0)
    end
    # Other misc options that over-ride the ones set above.
    for (name, value) in model.options
        if value isa String
            AddIpoptStrOption(model.inner, name, value)
        elseif value isa Integer
            AddIpoptIntOption(model.inner, name, value)
        else
            @assert value isa Float64
            AddIpoptNumOption(model.inner, name, value)
        end
    end
    # Initialize the starting point, projecting variables from 0 onto their
    # bounds if VariablePrimalStart  is not provided.
    for (i, v) in enumerate(model.variable_primal_start)
        if v !== nothing
            model.inner.x[i] = v
        else
            model.inner.x[i] = max(0.0, model.variables.lower[i])
            model.inner.x[i] = min(model.inner.x[i], model.variables.upper[i])
        end
    end
    # Initialize the dual start to 0.0 if NLPBlockDualStart is not provided.
    if model.nlp_dual_start === nothing
        model.nlp_dual_start = zeros(Float64, num_nlp_constraints)
    end
    # ConstraintDualStart
    row = 1
    for key in _CONSTRAINT_ORDERING
        for info in getfield(model, key)
            model.inner.mult_g[row] = _dual_start(model, info.dual_start, -1)
            row += 1
        end
    end
    for dual_start in model.nlp_dual_start
        model.inner.mult_g[row] = _dual_start(model, dual_start, -1)
        row += 1
    end
    # ConstraintDualStart for variable bounds
    for i in 1:length(model.inner.n)
        model.inner.mult_x_L[i] =
            _dual_start(model, model.variable_lower_start[i])
        model.inner.mult_x_U[i] =
            _dual_start(model, model.variable_upper_start[i], -1)
    end
    # Initialize CallbackFunction.
    if model.callback !== nothing
        SetIntermediateCallback(model.inner, model.callback)
    end
    IpoptSolve(model.inner)
    # Store SolveTimeSec.
    model.solve_time = time() - start_time
    return
end

# From Ipopt/src/Interfaces/IpReturnCodes_inc.h
const _STATUS_CODES = Dict(
    0 => :Solve_Succeeded,
    1 => :Solved_To_Acceptable_Level,
    2 => :Infeasible_Problem_Detected,
    3 => :Search_Direction_Becomes_Too_Small,
    4 => :Diverging_Iterates,
    5 => :User_Requested_Stop,
    6 => :Feasible_Point_Found,
    -1 => :Maximum_Iterations_Exceeded,
    -2 => :Restoration_Failed,
    -3 => :Error_In_Step_Computation,
    -4 => :Maximum_CpuTime_Exceeded,
    -5 => :Maximum_WallTime_Exceeded,
    -10 => :Not_Enough_Degrees_Of_Freedom,
    -11 => :Invalid_Problem_Definition,
    -12 => :Invalid_Option,
    -13 => :Invalid_Number_Detected,
    -100 => :Unrecoverable_Exception,
    -101 => :NonIpopt_Exception_Thrown,
    -102 => :Insufficient_Memory,
    -199 => :Internal_Error,
)

### MOI.ResultCount

# Ipopt always has an iterate available.
function MOI.get(model::Optimizer, ::MOI.ResultCount)
    return (model.inner !== nothing) ? 1 : 0
end

### MOI.TerminationStatus

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    if model.invalid_model
        return MOI.INVALID_MODEL
    elseif model.inner === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    end
    status = _STATUS_CODES[model.inner.status]
    if status == :Solve_Succeeded || status == :Feasible_Point_Found
        return MOI.LOCALLY_SOLVED
    elseif status == :Infeasible_Problem_Detected
        return MOI.LOCALLY_INFEASIBLE
    elseif status == :Solved_To_Acceptable_Level
        return MOI.ALMOST_LOCALLY_SOLVED
    elseif status == :Search_Direction_Becomes_Too_Small
        return MOI.NUMERICAL_ERROR
    elseif status == :Diverging_Iterates
        return MOI.NORM_LIMIT
    elseif status == :User_Requested_Stop
        return MOI.INTERRUPTED
    elseif status == :Maximum_Iterations_Exceeded
        return MOI.ITERATION_LIMIT
    elseif status == :Maximum_CpuTime_Exceeded
        return MOI.TIME_LIMIT
    elseif status == :Maximum_WallTime_Exceeded
        return MOI.TIME_LIMIT
    elseif status == :Restoration_Failed
        return MOI.NUMERICAL_ERROR
    elseif status == :Error_In_Step_Computation
        return MOI.NUMERICAL_ERROR
    elseif status == :Invalid_Option
        return MOI.INVALID_OPTION
    elseif status == :Not_Enough_Degrees_Of_Freedom
        return MOI.INVALID_MODEL
    elseif status == :Invalid_Problem_Definition
        return MOI.INVALID_MODEL
    elseif status == :Invalid_Number_Detected
        return MOI.INVALID_MODEL
    elseif status == :Unrecoverable_Exception
        return MOI.OTHER_ERROR
    elseif status == :NonIpopt_Exception_Thrown
        return MOI.OTHER_ERROR
    else
        @assert status == :Insufficient_Memory
        return MOI.MEMORY_LIMIT
    end
end

### MOI.RawStatusString

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    if model.invalid_model
        return "The model has no variable"
    elseif model.inner === nothing
        return "Optimize not called"
    else
        return string(_STATUS_CODES[model.inner.status])
    end
end

### MOI.PrimalStatus

function MOI.get(model::Optimizer, attr::MOI.PrimalStatus)
    if !(1 <= attr.result_index <= MOI.get(model, MOI.ResultCount()))
        return MOI.NO_SOLUTION
    end
    status = _STATUS_CODES[model.inner.status]
    if status == :Solve_Succeeded
        return MOI.FEASIBLE_POINT
    elseif status == :Feasible_Point_Found
        return MOI.FEASIBLE_POINT
    elseif status == :Solved_To_Acceptable_Level
        # Solutions are only guaranteed to satisfy the "acceptable" convergence
        # tolerances.
        return MOI.NEARLY_FEASIBLE_POINT
    elseif status == :Infeasible_Problem_Detected
        return MOI.INFEASIBLE_POINT
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

### MOI.DualStatus

function MOI.get(model::Optimizer, attr::MOI.DualStatus)
    if !(1 <= attr.result_index <= MOI.get(model, MOI.ResultCount()))
        return MOI.NO_SOLUTION
    end
    status = _STATUS_CODES[model.inner.status]
    if status == :Solve_Succeeded
        return MOI.FEASIBLE_POINT
    elseif status == :Feasible_Point_Found
        return MOI.FEASIBLE_POINT
    elseif status == :Solved_To_Acceptable_Level
        # Solutions are only guaranteed to satisfy the "acceptable" convergence
        # tolerances.
        return MOI.NEARLY_FEASIBLE_POINT
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

### MOI.SolveTimeSec

MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = model.solve_time

### MOI.ObjectiveValue

function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    return model.inner.obj_val
end

### MOI.VariablePrimal

function MOI.get(
    model::Optimizer,
    attr::MOI.VariablePrimal,
    vi::MOI.VariableIndex,
)
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, vi)
    return model.inner.x[column(vi)]
end

### MOI.ConstraintPrimal

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{F,S},
) where {
    F<:Union{
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    },
    S,
}
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, ci)
    return model.inner.g[_offset(model, F, S)+ci.value]
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{
        MOI.VariableIndex,
        <:Union{
            MOI.LessThan{Float64},
            MOI.GreaterThan{Float64},
            MOI.EqualTo{Float64},
        },
    },
)
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, ci)
    return model.inner.x[ci.value]
end

### MOI.ConstraintDual

_dual_multiplier(model::Optimizer) = model.sense == MOI.MIN_SENSE ? 1.0 : -1.0

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{F,S},
) where {
    F<:Union{
        MOI.ScalarAffineFunction{Float64},
        MOI.ScalarQuadraticFunction{Float64},
    },
    S,
}
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, ci)
    s = -_dual_multiplier(model)
    return s * model.inner.mult_g[_offset(model, F, S)+ci.value]
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}},
)
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, ci)
    rc = model.inner.mult_x_L[ci.value] - model.inner.mult_x_U[ci.value]
    return min(0.0, _dual_multiplier(model) * rc)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}},
)
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, ci)
    rc = model.inner.mult_x_L[ci.value] - model.inner.mult_x_U[ci.value]
    return max(0.0, _dual_multiplier(model) * rc)
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}},
)
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, ci)
    rc = model.inner.mult_x_L[ci.value] - model.inner.mult_x_U[ci.value]
    return _dual_multiplier(model) * rc
end

### MOI.NLPBlockDual

function MOI.get(model::Optimizer, attr::MOI.NLPBlockDual)
    MOI.check_result_index_bounds(model, attr)
    s = -_dual_multiplier(model)
    return s .* model.inner.mult_g[(1+_nlp_constraint_offset(model)):end]
end

### Ipopt.CallbackFunction

"""
    CallbackFunction()

A solver-dependent callback for Ipopt's IntermediateCallback.

The callback should be a function like the following:
```julia
function my_intermediate_callback(
    alg_mod::Cint,
    iter_count::Cint,
    obj_value::Float64,
    inf_pr::Float64,
    inf_du::Float64,
    mu::Float64,
    d_norm::Float64,
    regularization_size::Float64,
    alpha_du::Float64,
    alpha_pr::Float64,
    ls_trials::Cint,
)
    # ... user code ...
    return true # or `return false` to terminate the solve.
end
```

The arguments are defined in the Ipopt documentation:
https://coin-or.github.io/Ipopt/OUTPUT.html

Note: Calling `SetIntermediateCallback` will over-write this callback! Don't
call both.
"""
struct CallbackFunction <: MOI.AbstractCallback end

function MOI.set(model::Optimizer, ::CallbackFunction, f::Function)
    model.callback = f
    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.CallbackVariablePrimal,
    x::MOI.VariableIndex,
)
    return model.inner.x[column(x)]
end


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

mutable struct WhatIfOperationResult <: SwaggerModel
    status::Any # spec type: Union{ Nothing, String } # spec name: status
    properties::Any # spec type: Union{ Nothing, WhatIfOperationProperties } # spec name: properties
    error::Any # spec type: Union{ Nothing, ErrorResponse } # spec name: error

    function WhatIfOperationResult(;status=nothing, properties=nothing, error=nothing)
        o = new()
        validate_property(WhatIfOperationResult, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        validate_property(WhatIfOperationResult, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        validate_property(WhatIfOperationResult, Symbol("error"), error)
        setfield!(o, Symbol("error"), error)
        o
    end
end # type WhatIfOperationResult

const _property_map_WhatIfOperationResult = Dict{Symbol,Symbol}(Symbol("status")=>Symbol("status"), Symbol("properties")=>Symbol("properties"), Symbol("error")=>Symbol("error"))
const _property_types_WhatIfOperationResult = Dict{Symbol,String}(Symbol("status")=>"String", Symbol("properties")=>"WhatIfOperationProperties", Symbol("error")=>"ErrorResponse")
Base.propertynames(::Type{ WhatIfOperationResult }) = collect(keys(_property_map_WhatIfOperationResult))
Swagger.property_type(::Type{ WhatIfOperationResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_WhatIfOperationResult[name]))}
Swagger.field_name(::Type{ WhatIfOperationResult }, property_name::Symbol) =  _property_map_WhatIfOperationResult[property_name]

function check_required(o::WhatIfOperationResult)
    true
end

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


text:
using Pkg

Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

using NBInclude

cd(@__DIR__) do
    @nbinclude("Base_Tables.ipynb")
    @nbinclude("GlobalWgt_Tables.ipynb")
    @nbinclude("IncomeElas0_Tables.ipynb")
    @nbinclude("Poverty_Tables.ipynb")
    @nbinclude("VSL0_Tables.ipynb")
end


text:
"""
    NaNChecker{F} <: AbstractDiagnostic

A diagnostic that checks for `NaN` values and aborts the simulation if any are found.
"""
struct NaNChecker{T, F} <: AbstractDiagnostic
    schedule :: T
      fields :: F
end

"""
    NaNChecker(; schedule, fields)

Returns a `NaNChecker` that checks for `NaN` anywhere within `fields`
when `schedule` actuates.
"""
NaNChecker(model=nothing; schedule, fields) = NaNChecker(schedule, fields)

function run_diagnostic!(nc::NaNChecker, model)
    for (name, field) in nc.fields
        if any(isnan, field.data.parent)
            t, i = model.clock.time, model.clock.iteration
            error("time = $t, iteration = $i: NaN found in $name. Aborting simulation.")
        end
    end
end


text:
using SafeTestsets

@safetestset "Memory Tests" begin

    @safetestset "baselevel" begin
        using ACTRModels, Test
        import ACTRModels: baselevel, baselevel!
        include("memory.jl")

        lags = [1]
        d = .5
        @test baselevel(d, lags) == 0

        lags = [1]
        d = .1
        @test baselevel(d, lags) == 0

        lags = [2,10]
        d = .5
        @test baselevel(d, lags) ≈ 0.023066 atol=1e-5

        lags = [.5,10]
        d = .4
        @test baselevel(d, lags) ≈ 0.540936 atol=1e-5
    end

    @safetestset "baselevel!" begin
        using ACTRModels, Test
        import ACTRModels: baselevel, baselevel!
        include("memory.jl")

        actr,memory,chunks = initializeACTR(;bll=true)
        c = get_chunks(actr; animal=:dog,name=:Sigma)[1]
        update_lags!(c, 2.0)
        baselevel!(actr, c)
        @test c.act_bll ≈ -0.346573 atol=1e-5
        c.act_bll = 0.0
        update_lags!(c, 5.0)
        baselevel!(actr, c)
        @test c.act_bll ≈ -0.80471 atol=1e-5
        c.act_bll = 0.0
        update_recent!(c, 10.0)
        update_chunk!(c, 12)
        update_lags!(c, 14)
        baselevel!(actr, c)
        @test c.act_bll ≈ 0.09076 atol=1e-5
        c.act_bll = 0.0
        actr.parms.d = .8
        baselevel!(actr, c)
        @test c.act_bll ≈ -0.220564 atol=1e-5

        actr,memory,chunks = initializeACTR(; bll=true)
        c = get_chunks(actr; animal=:dog, name=:Sigma)[1]
        c.k = 2
        update_chunk!(c, 1)
        update_chunk!(c, 2)
        update_chunk!(c, 4)
        update_lags!(c, 10)
        @test c.N == 4
        @test c.L == 10
        @test all(c.lags .== [6,8])
        @test all(c.recent .== [4,2])
    end

    @safetestset "Chunk" begin
        using ACTRModels, Test
        import ACTRModels: baselevel, baselevel!
        include("memory.jl")

        chunk = Chunk(;N=20, animal=:dog, name=:Sigma)
        @test chunk.slots.animal == :dog
        @test chunk.slots.name == :Sigma
        @test chunk.N == 20
    end

    @safetestset "add Chunk" begin
        using ACTRModels, Test
        import ACTRModels: baselevel, baselevel!
        include("memory.jl")

        actr,memory,chunks = initializeACTR()
        c = get_chunks(actr; animal=:human, name=:Wilford, last_name=:Brimley)
        @test isempty(c)
        add_chunk!(actr, 10.0; animal=:human, name=:Wilford, last_name=:Brimley)
        c = get_chunks(actr; animal=:human, name=:Wilford, last_name=:Brimley)
        @test !isempty(c)
        @test c[1].N == 1
        @test c[1].time_created == 10.0
        @test c[1].recent[1] == 10.0

        add_chunk!(actr, 20.0; animal=:human, name=:Wilford, last_name=:Brimley)
        @test c[1].N == 2
        @test c[1].time_created == 10.0
        @test c[1].recent[1] == 20.0
    end

    @safetestset "retrieval_prob" begin
        using ACTRModels, Test
        import ACTRModels: baselevel, baselevel!
        include("memory.jl")

        actr,memory,chunks = initializeACTR()
        @test retrieval_prob(actr, chunks[2]; animal=:cat, name=:Butters) == (.5,.5)
        @test retrieval_prob(actr, chunks[1]; isa=:mammal) == (0,1)

        actr,memory,chunks = initializeACTR(;τ=.5)
        c = get_chunks(actr; animal=:rat)
        @test retrieval_prob(actr, c; animal=:rat)[1] ≈ 0.38098 atol=1e-5

        actr,memory,chunks = initializeACTR(;τ=.5, mmp=true, δ=1.0)
        c = get_chunks(actr; animal=:rat, name=:Joy)
        p1,_ = retrieval_prob(actr, c)
        @test p1 ≈ 0.137939 atol=1e-5

        p2,_ = retrieval_prob(actr, c; animal=:rat)
        @test p2 ≈ 0.183859 atol=1e-5

        p3,_ = retrieval_prob(actr, c; animal=:rat, name=:Joy)
        @test p3 ≈ 0.229243 atol=1e-5

        p3,_ = retrieval_prob(actr, c; isa=:rock, animal=:rat, name=:Joy)
        @test p3 == 0

        actr,memory,chunks = initializeACTR(;τ=.5, bll=true)
        c = get_chunks(actr; animal=:rat, name=:Joy)
        update_lags!(actr, 1.0)
        p1,_ = retrieval_prob(actr, c, 1.0)
        update_lags!(actr, 3.0)
        p2,_ = retrieval_prob(actr, c, 3.0)
        actr.parms.d = .8
        p3,_ = retrieval_prob(actr, c, 3.0)
        @test p1 > p2 > p3

        actr,memory,chunks = initializeACTR(;d=.5, bll=true)
        c = get_chunks(actr; animal=:rat, name=:Joy)
        update_recent!.(c, 5.0)
        update_lags!(actr, 7)
        p1,_ = retrieval_prob(actr, c, 7.0)
        c[1].N = 2
        p2,_=retrieval_prob(actr, c, 7.0)
        c[1].N = 5
        p3,_ = retrieval_prob(actr, c, 7.0)
        @test p1 < p2 < p3

        actr,memory,chunks = initializeACTR(;sa=true, γ=1.0)
        c = get_chunks(actr; animal=:rat, name=:Joy)
        actr.imaginal.buffer[1] = Chunk(;animal=:rat)
        p1,_ = retrieval_prob(actr, c)
        actr.imaginal.buffer[1] = Chunk(;name=:Joy)
        p2,_ = retrieval_prob(actr, c)
        @test p2 > p1
    end

    @safetestset "modify!" begin
        using ACTRModels, Test
        chunks = Chunk[Chunk(;isa=:bafoon, animal=:dog,name=:Sigma, retrieved=[false]),
        Chunk(;isa=:mammal,animal=:cat,name=:Butters, retrieved=[false])]
        memory = Declarative(;memory=chunks)
        actr = ACTR(;declarative=memory)
        modify!(chunks[2].slots, retrieved=true)
        @test chunks[2].slots.retrieved[1] == true
    end

    @safetestset "filter" begin
        using ACTRModels, Test
        chunks = Chunk[Chunk(;isa=:bafoon, animal=:dog,name=:Sigma, retrieved=[false]),
        Chunk(;isa=:mammal,animal=:cat,name=:Butters, retrieved=[false])]
        memory = Declarative(;memory=chunks)
        actr = ACTR(;declarative=memory, mmp=true)
        request = retrieval_request(actr; isa=:mammal)
        @test request[1].slots.name == :Butters
    end

    @safetestset "count_values" begin
        using ACTRModels, Test
        import ACTRModels: count_values
        chunk = Chunk(a=:a, b=:a, c=:v)
        @test count_values(chunk, :a) == 2
        @test count_values(chunk, :v) == 1
        @test count_values(chunk, :c) == 0
    end

    @safetestset "update_recent!" begin
        using ACTRModels, Test
        chunk = Chunk(a=:a, b=:a, c=:v)
        update_recent!(chunk, .10)
        @test chunk.recent[1] == .1
        update_recent!(chunk, .20)
        @test chunk.recent[1] == .2
        @test length(chunk.recent) == 1
        chunk.k = 2
        update_recent!(chunk, .10)
        update_recent!(chunk, .20)
        update_recent!(chunk, .30)
        @test chunk.recent[1] == .3
        @test chunk.recent[2] == .2
        @test length(chunk.recent) == 2
    end

    @safetestset "update_lags!" begin
        using ACTRModels, Test
        chunk = Chunk(a=:a, b=:a, c=:v, k=2)
        update_recent!(chunk, .10)
        update_recent!(chunk, .20)
        update_lags!(chunk, 1.0)
        @test chunk.lags[1] == 1 - .2
        @test chunk.lags[2] == 1 - .1
    end

    @safetestset "update_chunk!" begin
        using ACTRModels, Test
        chunk = Chunk(a=:a, b=:a, c=:v)
        update_chunk!(chunk, .10)
        @test chunk.recent[1] == .1
        @test chunk.N == 2
        update_chunk!(chunk, .20)
        @test chunk.recent[1] == .2
        @test chunk.N == 3
        @test length(chunk.recent) == 1
    end

    @safetestset "spreading activation" begin
        using ACTRModels, Test
        chunks = [Chunk(a=:a, b=:b, c=:c), Chunk(a=:a, b=:b, c=:a)]
        memory = Declarative(memory=chunks)
        imaginal = Imaginal(buffer=Chunk(a=:a, b=:b))
        actr = ACTR(declarative=memory, imaginal=imaginal, γ=1.6)
        compute_activation!(actr)
        @test chunks[1].act == 0.0
        @test chunks[2].act == 0.0
        actr.parms.sa = true
        compute_activation!(actr)
        @test chunks[1].act ≈ 0.3575467 atol = 1e-5
        @test chunks[2].act ≈ 0.7041203 atol = 1e-5
    end

    @safetestset "spreading activation zero" begin
        using ACTRModels, Test
        chunks = [Chunk(;isa=:mammal, animal=:dog,name=:Sigma),
        Chunk(;isa=:mammal,animal=:cat,name=:Butters)]
        memory = Declarative(;memory=chunks)
        actr = ACTR(;declarative=memory, sa=true, γ=1.0)
        p,_ = retrieval_probs(actr)
        sa = map(x->x.act_sa, chunks)
        @test sa[1] ≈ 0.0 atol = 1e-10
        @test sa[2] ≈ 0.0 atol = 1e-10

        memory = Declarative(;memory=chunks)
        T = typeof(chunks)
        imaginal = Imaginal(buffer=T)
        empty!(imaginal.buffer)
        actr = ACTR(;declarative=memory, imaginal=imaginal, sa=true, γ=1.0)
        p,_ = retrieval_probs(actr)
        sa = map(x->x.act_sa, chunks)
        @test sa[1] ≈ 0.0 atol = 1e-10
        @test sa[2] ≈ 0.0 atol = 1e-10

        chunk = Chunk(;isa=:mammal,animal=:cat,name=:Butters)
        push!(imaginal.buffer, chunk)
        p,_ = retrieval_probs(actr)
        sa = map(x->x.act_sa, chunks)
        @test sa[1] != 0.0
        @test sa[2] != 0.0

        empty!(imaginal.buffer)
        p,_ = retrieval_probs(actr)
        sa = map(x->x.act_sa, chunks)
        @test sa[1] ≈ 0.0 atol = 1e-10
        @test sa[2] ≈ 0.0 atol = 1e-10

        chunk = Chunk(;isa=:arachnid,animal=:spider,name=:Soen)
        push!(imaginal.buffer, chunk)
        p,_ = retrieval_probs(actr)
        sa = map(x->x.act_sa, chunks)
        @test sa[1] ≈ 0.0 atol = 1e-10
        @test sa[2] ≈ 0.0 atol = 1e-10
    end

    @safetestset "partial matching" begin
        using ACTRModels, Test
        chunk = Chunk(a=:a, b=:b)
        memory = Declarative(memory=[chunk])
        actr = ACTR(declarative=memory, mmp=true, δ=1.0)
        compute_activation!(actr)
        @test chunk.act == 0.0
        compute_activation!(actr; a=:a)
        @test chunk.act == 0.0
        compute_activation!(actr; a=:a, b=:b)
        @test chunk.act == 0.0
        compute_activation!(actr; a=:b, b=:b)
        @test chunk.act == -1.0
        compute_activation!(actr; a=:a, b=:a)
        @test chunk.act == -1.0
        compute_activation!(actr; a=:c, b=:c)
        @test chunk.act == -2.0
    end

    @safetestset "get_chunks" begin
        using ACTRModels, Test
        chunks = [Chunk(a=:a, b=:c),Chunk(a=:a, b=:b)]
        memory = Declarative(memory=chunks)
        actr = ACTR(declarative=memory)
        result = get_chunks(actr; a=:b)
        @test isempty(result)

        result = get_chunks(actr; a=:a)
        @test !isempty(result)
        @test length(result) == 2
        @test result[1].slots.a == :a
        @test result[1].slots.b == :c
        @test result[2].slots.a == :a
        @test result[2].slots.b == :b

       result = get_chunks(actr; a=:a, b=:b)
       @test !isempty(result)
       @test length(result) == 1
       @test result[1].slots.a == :a
       @test result[1].slots.b == :b
    end

    @safetestset "first_chunk" begin
        using ACTRModels, Test
        chunks = [Chunk(a=:a, b=:c),Chunk(a=:a, b=:b)]
        memory = Declarative(memory=chunks)
        actr = ACTR(declarative=memory)
        result = get_chunks(actr; a=:b)
        @test isempty(result)

        result = first_chunk(actr; a=:a)
        @test !isempty(result)
        @test length(result) == 1
        @test result[1].slots.a == :a
        @test result[1].slots.b == :c
    end

    @safetestset "reset_activation!" begin
        using ACTRModels, Test
        import ACTRModels: reset_activation!
        chunks = [Chunk(a=:a, b=:b, c=:c), Chunk(a=:a, b=:b, c=:a)]
        chunk = chunks[1]
        memory = Declarative(memory=chunks)
        imaginal = Imaginal(buffer=Chunk(a=:a, b=:b))
        actr = ACTR(declarative=memory, imaginal=imaginal, mmp=true, noise=true, bll=true,
            sa=true, γ=1.6, δ=1.0)
        compute_activation!(actr, 3.0; a=:b)
        @test chunk.act_bll != 0
        @test chunk.act_pm != 0
        @test chunk.act_sa != 0
        @test chunk.act_noise != 0
        reset_activation!(chunk)
        @test chunk.act_bll == 0
        @test chunk.act_pm == 0
        @test chunk.act_sa == 0
        @test chunk.act_noise == 0
    end

    @safetestset "total_activation!" begin
        using ACTRModels, Test
        import ACTRModels: total_activation!
        chunk = Chunk()
        chunk.act_blc = .5
        chunk.act_bll = 1.0
        chunk.act_pm = 1.0
        chunk.act_sa = 1.0
        chunk.act_noise = 1.0
        total_activation!(chunk)
        @test chunk.act ≈ 2.5 atol = 1e-5
    end

    @safetestset "get_parm" begin
        using ACTRModels, Test
        chunks = [Chunk(a=:a, b=:b, c=:c), Chunk(a=:a, b=:b, c=:a)]
        chunk = chunks[1]
        memory = Declarative(memory=chunks)
        actr = ACTR(declarative=memory, mmp=true, noise=true, bll=true, sa=true, γ=1.6,
        δ=1.0, Σ=.3)
        @test get_parm(actr, :δ) ≈ 1.0 atol = 1e-5
        @test get_parm(actr, :Σ) ≈ 0.3 atol = 1e-5
        @test get_parm(actr, :noise)
    end

    @safetestset "compute_RT" begin
        using ACTRModels, Test, Distributions, Random
        Random.seed!(41140)
        chunks = [Chunk(a=:a, b=:b, c=:c), Chunk(a=:a, b=:b, c=:a)]
        chunk = chunks[1]
        memory = Declarative(memory=chunks)
        actr = ACTR(declarative=memory, blc=1.5)
        retrieved = retrieve(actr; a=:a, b=:b, c=:c)
        rt = compute_RT(actr, retrieved)
        @test rt ≈ exp(-1.5) atol = 1e-5

        retrieved = retrieve(actr; a=:z)
        rt = compute_RT(actr, retrieved)
        @test rt ≈ exp(0.0) atol = 1e-5

        chunks = [Chunk(a=:a, b=:b, c=:c), Chunk(a=:a, b=:b, c=:a)]
        chunk = chunks[1]
        s = .2
        memory = Declarative(memory=chunks)
        actr = ACTR(declarative=memory, blc=1.5, noise=true, s=s)
        retrieved = retrieve(actr; a=:a, b=:b, c=:c)
        function sim(actr, chunk)
            compute_activation!(actr, chunk)
            return compute_RT(actr, retrieved)
        end
        σ = s * pi / sqrt(3)
        rts = map(_->sim(actr, retrieved), 1:20_000)
        @test mean(rts) ≈ mean(LogNormal(-1.5, σ)) atol = 1e-3
        @test std(rts) ≈ std(LogNormal(-1.5, σ)) atol = 1e-3
    end

    @safetestset "match" begin
        using ACTRModels, Test
        chunk = Chunk(a=:a, b=:b, c=:c)
        @test match(chunk, a=:a)
        @test !match(chunk, a=:b)
        @test !match(chunk, d=:b)
        @test match(chunk, !=, ==, a=:b, b=:b)
        @test !match(chunk, !=, ==, a=:a, b=:b)

        @test match(chunk, a=:a; check_value=false)
        @test match(chunk, a=:b; check_value=false)
        @test !match(chunk, d=:a; check_value=false)
    end

    @safetestset "Threshold" begin
        using ACTRModels, Test, Distributions, Random
        Random.seed!(41140)
        chunks = [Chunk(a=:a, b=:b, c=:c, bl=-1.0), Chunk(a=:a, b=:b, c=:a, bl=-1.0)]
        memory = Declarative(memory=chunks)
        τ = 0.0
        actr = ACTR(;declarative=memory, τ, noise=false)
        retrieved = retrieve(actr; a=:a)
        @test isempty(retrieved)
        @test actr.parms.τ′ == τ

        τ = -2.0
        actr = ACTR(;declarative=memory, τ, noise=false)
        retrieved = retrieve(actr; a=:a)
        @test !isempty(retrieved)
        @test actr.parms.τ′ == τ

        τ = -20.0
        actr = ACTR(;declarative=memory, τ, noise=true)
        retrieved = retrieve(actr; a=:a)
        @test !isempty(retrieved)
        @test actr.parms.τ′ != τ
    end

    @safetestset "base level learning example" begin
        using ACTRModels, Test, Distributions, Random
        memory = Declarative(memory=Chunk[])
        bll = true
        d = .5
        noise = false
        actr = ACTR(;declarative=memory, bll, d, noise)
        cur_time = 45.185
        for _ in 1:10
            add_chunk!(actr, cur_time; a=:a)
            cur_time += rand(Uniform(0, 100))
        end
        cur_time = 1265.185
        add_chunk!(actr, cur_time; a=:a)
        cur_time = 1590.185
        compute_activation!(actr, cur_time)
        chunk = actr.declarative.memory[1]
        @test chunk.act_bll ≈ -.905594 rtol = .0001
        # based on lisp act-r example:
        # 1590.185   DECLARATIVE            START-RETRIEVAL
        # Computing base-level from 11 references (1265.185)
        #   creation time: 45.185 decay: 0.5  Optimized-learning: 1
        # base-level value: -0.905594
    end

    @safetestset "blend_chunks" begin
        using ACTRModels, Test, Random
        Random.seed!(598)
        chunks = [Chunk(;a=1, b=0), Chunk(;a=1, b=3)]
        parms = (mmp = true, δ=1.0, noise=true, s=.2)
        declarative = Declarative(;memory=chunks)
        actr = ACTR(;declarative, parms...)

        request = (a=2,)
        blended_slots = :b
        n_sim = 10_000
        mean_value1 = map(_->blend_chunks(actr, blended_slots; request...), 1:n_sim) |> mean
        @test mean_value1 ≈ 1.5 atol = .01

        chunks = [Chunk(;a=1, b=0), Chunk(;a=2, b=3)]
        parms = (mmp = true, δ=1.0, noise=true, s=.2)
        declarative = Declarative(;memory=chunks)
        actr = ACTR(;declarative, parms...)
        request = (a=2,)
        blended_slots = :b
        n_sim = 10_000
        mean_value2 = map(_->blend_chunks(actr, blended_slots; request...), 1:n_sim) |> mean
        @test mean_value1 < mean_value2

        chunks = [Chunk(;a=2, b=0), Chunk(;a=1, b=3)]
        parms = (mmp = true, δ=1.0, noise=true, s=.2)
        declarative = Declarative(;memory=chunks)
        actr = ACTR(;declarative, parms...)
        request = (a=2,)
        blended_slots = :b
        n_sim = 10_000
        mean_value3 = map(_->blend_chunks(actr, blended_slots; request...), 1:n_sim) |> mean
        @test mean_value1 > mean_value3

        chunks = [Chunk(;a=2, b=0), Chunk(;a=1, b=3)]
        parms = (mmp = true, δ=.5, noise=true, s=.2)
        declarative = Declarative(;memory=chunks)
        actr = ACTR(;declarative, parms...)
        request = (a=2,)
        blended_slots = :b
        n_sim = 10_000
        mean_value4 = map(_->blend_chunks(actr, blended_slots; request...), 1:n_sim) |> mean
        @test mean_value4 > mean_value3
    end

    @safetestset "blend_slots" begin
        using ACTRModels, Test, Random
        import ACTRModels: blend_slots
        Random.seed!(598)
        chunks = [Chunk(;a=1, b=0), Chunk(;a=1, b=3)]
        parms = (mmp = true, δ=1.0, noise=true, s=.2)
        declarative = Declarative(;memory=chunks)
        actr = ACTR(;declarative, parms...)

        blended_slots = :b
        probs = [.3,.7]
        v = blend_slots(chunks, probs, blended_slots)
        @test v ≈ 2.1 atol = 1e-4
    end

    @safetestset "get_chunks_exact" begin
        using ACTRModels, Test
        using ACTRModels: get_chunks_exact
        chunks = [Chunk(;a=1, b=0)]
        declarative = Declarative(;memory=chunks)

        result = get_chunks_exact(declarative; a=1, b=0)
        @test !isempty(result)

        result = get_chunks_exact(declarative; a=1, b=3)
        @test isempty(result)

        result = get_chunks_exact(declarative; a=1)
        @test isempty(result)
    end
end


text:
# Append (UNION ALL) node.

mutable struct AppendNode <: TabularNode
    over::Union{SQLNode, Nothing}
    args::Vector{SQLNode}

    AppendNode(; over = nothing, args) =
        new(over, args)
end

AppendNode(args...; over = nothing) =
    AppendNode(over = over, args = SQLNode[args...])

"""
    Append(; over = nothing, args)
    Append(args...; over = nothing)

`Append` concatenates input datasets.

```sql
SELECT ...
FROM \$over
UNION ALL
SELECT ...
FROM \$(args[1])
UNION ALL
...
```

# Examples

```jldoctest
julia> measurement = SQLTable(:measurement, columns = [:measurement_id, :person_id, :measurement_date]);

julia> observation = SQLTable(:observation, columns = [:observation_id, :person_id, :observation_date]);

julia> q = From(measurement) |>
           Define(:date => Get.measurement_date) |>
           Append(From(observation) |>
                  Define(:date => Get.observation_date)) |>
           Select(Get.person_id, Get.date);

julia> print(render(q))
SELECT
  "union_1"."person_id",
  "union_1"."date"
FROM (
  SELECT
    "measurement_1"."person_id",
    "measurement_1"."measurement_date" AS "date"
  FROM "measurement" AS "measurement_1"
  UNION ALL
  SELECT
    "observation_1"."person_id",
    "observation_1"."observation_date" AS "date"
  FROM "observation" AS "observation_1"
) AS "union_1"
```
"""
Append(args...; kws...) =
    AppendNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(Append), pats::Vector{Any}) =
    dissect(scr, AppendNode, pats)

function PrettyPrinting.quoteof(n::AppendNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(Append))
    if isempty(n.args)
        push!(ex.args, Expr(:kw, :args, Expr(:vect)))
    else
        append!(ex.args, quoteof(n.args, ctx))
    end
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

function label(n::AppendNode)
    lbl = label(n.over)
    for l in n.args
        label(l) === lbl || return :union
    end
    lbl
end

rebase(n::AppendNode, n′) =
    AppendNode(over = rebase(n.over, n′), args = n.args)