text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell
include("useful_routines.jl")

function bulk_modulus(l, m)
    return l + 2m / 3
end

function lame_lambda(k, m)
    return k - 2m / 3
end

function analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
    a = zeros(3, 3)
    a[1, 1] = inradius
    a[1, 2] = -inradius
    a[1, 3] = -1.0 / inradius
    a[2, 1] = 2 * (lc + mc)
    a[2, 2] = -2 * (ls + ms)
    a[2, 3] = 2ms / inradius^2
    a[3, 2] = 2(ls + ms)
    a[3, 3] = -2ms / outradius^2
    return a
end

function analytical_coefficient_rhs(ls, ms, theta0)
    r = zeros(3)
    Ks = bulk_modulus(ls, ms)
    r[2] = -Ks * theta0
    r[3] = Ks * theta0
    return r
end

struct AnalyticalSolution
    inradius::Any
    outradius::Any
    center::Any
    A1c::Any
    A1s::Any
    A2s::Any
    ls::Any
    ms::Any
    lc::Any
    mc::Any
    theta0::Any
    function AnalyticalSolution(
        inradius,
        outradius,
        center,
        ls,
        ms,
        lc,
        mc,
        theta0,
    )
        a = analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
        r = analytical_coefficient_rhs(ls, ms, theta0)
        coeffs = a \ r
        new(
            inradius,
            outradius,
            center,
            coeffs[1],
            coeffs[2],
            coeffs[3],
            ls,
            ms,
            lc,
            mc,
            theta0,
        )
    end
end

function radial_displacement(A::AnalyticalSolution, r)
    if r <= A.inradius
        return A.A1c * r
    else
        return A.A1s * r + A.A2s / r
    end
end

function (A::AnalyticalSolution)(x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    ur = radial_displacement(A, r)
    if ur ≈ 0.0
        [0.0, 0.0]
    else
        costheta = (x[1] - A.center[1]) / r
        sintheta = (x[2] - A.center[2]) / r
        u1 = ur * costheta
        u2 = ur * sintheta
        return [u1, u2]
    end
end

function shell_radial_stress(ls, ms, theta0, A1, A2, r)
    return (ls + 2ms) * (A1 - A2 / r^2) + ls * (A1 + A2 / r^2) -
           (ls + 2ms / 3) * theta0
end

function shell_circumferential_stress(ls, ms, theta0, A1, A2, r)
    return ls * (A1 - A2 / r^2) + (ls + 2ms) * (A1 + A2 / r^2) -
           (ls + 2ms / 3) * theta0
end

function shell_out_of_plane_stress(ls, ms, A1, theta0)
    return 2 * ls * A1 - (ls + 2ms / 3) * theta0
end

function core_in_plane_stress(lc, mc, A1)
    return (lc + 2mc) * A1 + lc * A1
end

function core_out_of_plane_stress(lc, A1)
    return 2 * lc * A1
end

function rotation_matrix(x, r)
    costheta = x[1] / r
    sintheta = x[2] / r
    Q = [
        costheta -sintheta
        sintheta costheta
    ]
    return Q
end

function shell_stress(A::AnalyticalSolution, x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    Q = rotation_matrix(relpos, r)

    srr = shell_radial_stress(A.ls, A.ms, A.theta0, A.A1s, A.A2s, r)
    stt = shell_circumferential_stress(A.ls, A.ms, A.theta0, A.A1s, A.A2s, r)

    cylstress = [
        srr 0.0
        0.0 stt
    ]

    cartstress = Q * cylstress * Q'
    s11 = cartstress[1, 1]
    s22 = cartstress[2, 2]
    s12 = cartstress[1, 2]
    s33 = shell_out_of_plane_stress(A.ls, A.ms, A.A1s, A.theta0)

    return [s11, s22, s12, s33]
end

function core_stress(A::AnalyticalSolution)
    s11 = core_in_plane_stress(A.lc, A.mc, A.A1c)
    s33 = core_out_of_plane_stress(A.lc, A.A1c)
    return [s11, s11, 0.0, s33]
end

function exact_stress(A::AnalyticalSolution, x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    if r < A.inradius
        return core_stress(A)
    else
        return shell_stress(A, x)
    end
end

function onboundary(x, L, W)
    return x[2] ≈ 0.0 || x[1] ≈ L || x[2] ≈ W || x[1] ≈ 0.0
end

function update_core_stress_error!(
    err,
    basis,
    stiffness,
    celldisp,
    quad,
    jac,
    detjac,
    vectosymmconverter,
    analyticalsolution,
)
    dim = CutCell.dimension(basis)
    lambda, mu = CutCell.lame_coefficients(stiffness, -1)
    for (p, w) in quad
        numericalstress = CutCell.parent_stress(
            p,
            basis,
            stiffness,
            celldisp,
            jac,
            vectosymmconverter,
        )
        exactstress = core_stress(analyticalsolution)

        err .+= (numericalstress - exactstress) .^ 2 * detjac * w
    end
end

function update_interface_core_stress_error!(
    err,
    basis,
    stiffness,
    celldisp,
    quad,
    jac,
    facescale,
    vectosymmconverter,
    analyticalsolution,
)

    @assert length(facescale) == length(quad)
    dim = CutCell.dimension(basis)
    lambda, mu = CutCell.lame_coefficients(stiffness, -1)
    for (idx, (p, w)) in enumerate(quad)
        numericalstress = CutCell.parent_stress(
            p,
            basis,
            stiffness,
            celldisp,
            jac,
            vectosymmconverter,
        )
        exactstress = core_stress(analyticalsolution)

        err .+= (numericalstress - exactstress) .^ 2 * facescale[idx] * w
    end
end

function update_shell_stress_error!(
    err,
    basis,
    stiffness,
    transfstress,
    celldisp,
    quad,
    jac,
    detjac,
    cellmap,
    vectosymmconverter,
    analyticalsolution,
)

    dim = CutCell.dimension(basis)
    lambda, mu = CutCell.lame_coefficients(stiffness, +1)
    theta0 = analyticalsolution.theta0
    for (p, w) in quad
        numericalstress = CutCell.product_stress(
            p,
            basis,
            stiffness,
            transfstress,
            theta0,
            celldisp,
            jac,
            vectosymmconverter,
        )
        exactstress = shell_stress(analyticalsolution, cellmap(p))

        err .+= (numericalstress - exactstress) .^ 2 * detjac * w
    end
end

function update_interface_shell_stress_error!(
    err,
    basis,
    stiffness,
    transfstress,
    celldisp,
    quad,
    jac,
    facescale,
    cellmap,
    vectosymmconverter,
    analyticalsolution,
)
    @assert length(facescale) == length(quad)
    dim = CutCell.dimension(basis)
    lambda, mu = CutCell.lame_coefficients(stiffness, +1)
    theta0 = analyticalsolution.theta0
    for (idx, (p, w)) in enumerate(quad)
        numericalstress = CutCell.product_stress(
            p,
            basis,
            stiffness,
            transfstress,
            theta0,
            celldisp,
            jac,
            vectosymmconverter,
        )
        exactstress = shell_stress(analyticalsolution, cellmap(p))

        err .+= (numericalstress - exactstress) .^ 2 * facescale[idx] * w
    end
end

function compute_stress_error(
    displacement,
    basis,
    cellquads,
    stiffness,
    transfstress,
    cutmesh,
    analyticalsolution,
)
    err = zeros(4)
    numcells = CutCell.number_of_cells(cutmesh)

    dim = CutCell.dimension(basis)
    jac = CutCell.jacobian(cutmesh)
    detjac = CutCell.determinant_jacobian(cutmesh)
    vectosymmconverter = CutCell.vector_to_symmetric_matrix_converter()

    for cellid = 1:numcells
        s = CutCell.cell_sign(cutmesh, cellid)

        @assert s == -1 || s == 0 || s == +1
        if s == -1 || s == 0
            nodeids = CutCell.nodal_connectivity(cutmesh, -1, cellid)
            celldofs = CutCell.element_dofs(nodeids, dim)
            celldisp = displacement[celldofs]
            quad = cellquads[-1, cellid]

            update_core_stress_error!(
                err,
                basis,
                stiffness,
                celldisp,
                quad,
                jac,
                detjac,
                vectosymmconverter,
                analyticalsolution,
            )
        end

        if s == +1 || s == 0
            nodeids = CutCell.nodal_connectivity(cutmesh, +1, cellid)
            celldofs = CutCell.element_dofs(nodeids, dim)
            celldisp = displacement[celldofs]
            quad = cellquads[+1, cellid]
            cellmap = CutCell.cell_map(cutmesh, cellid)

            update_shell_stress_error!(
                err,
                basis,
                stiffness,
                transfstress,
                celldisp,
                quad,
                jac,
                detjac,
                cellmap,
                vectosymmconverter,
                analyticalsolution,
            )
        end
    end

    den = integral_norm_on_cut_mesh(
        x -> exact_stress(analyticalsolution, x),
        cellquads,
        cutmesh,
        4,
    )

    return sqrt.(err) ./ den
end

function compute_interface_stress_error(
    displacement,
    basis,
    interfacequads,
    stiffness,
    transfstress,
    cutmesh,
    analyticalsolution,
)

    perr = zeros(4)
    nerr = zeros(4)

    cellsign = CutCell.cell_sign(cutmesh)
    cellids = findall(cellsign .== 0)
    dim = CutCell.dimension(cutmesh)
    jac = CutCell.jacobian(cutmesh)
    vectosymmconverter = CutCell.vector_to_symmetric_matrix_converter()

    for cellid in cellids
        cellmap = CutCell.cell_map(cutmesh, cellid)
        normals = CutCell.interface_normals(interfacequads, cellid)
        facescale = CutCell.scale_area(cellmap, normals)

        nnodeids = CutCell.nodal_connectivity(cutmesh, -1, cellid)
        celldofs = CutCell.element_dofs(nnodeids, dim)
        celldisp = displacement[celldofs]
        nquad = interfacequads[-1, cellid]

        update_interface_core_stress_error!(
            nerr,
            basis,
            stiffness,
            celldisp,
            nquad,
            jac,
            facescale,
            vectosymmconverter,
            analyticalsolution,
        )

        pnodeids = CutCell.nodal_connectivity(cutmesh, +1, cellid)
        celldofs = CutCell.element_dofs(pnodeids, dim)
        celldisp = displacement[celldofs]
        pquad = interfacequads[+1, cellid]

        update_interface_shell_stress_error!(
            perr,
            basis,
            stiffness,
            transfstress,
            celldisp,
            pquad,
            jac,
            facescale,
            cellmap,
            vectosymmconverter,
            analyticalsolution,
        )
    end

    nden = integral_norm_on_interface(
        x -> core_stress(analyticalsolution),
        interfacequads,
        -1,
        cutmesh,
        4,
    )
    pden = integral_norm_on_interface(
        x -> shell_stress(analyticalsolution, x),
        interfacequads,
        +1,
        cutmesh,
        4,
    )

    return perr ./ pden, nerr ./ nden

end

function solve_and_compute_stress_error(
    width,
    center,
    inradius,
    outradius,
    stiffness,
    theta0,
    nelmts,
    polyorder,
    numqp,
    penaltyfactor,
)

    lambda1, mu1 = CutCell.lame_coefficients(stiffness, +1)
    lambda2, mu2 = CutCell.lame_coefficients(stiffness, -1)
    transfstress =
        CutCell.plane_strain_transformation_stress(lambda1, mu1, theta0)

    analyticalsolution = AnalyticalSolution(
        inradius,
        outradius,
        center,
        lambda1,
        mu1,
        lambda2,
        mu2,
        theta0,
    )

    dx = width / nelmts
    meanmoduli = 0.5 * (lambda1 + lambda2 + mu1 + mu2)
    penalty = penaltyfactor / dx * meanmoduli

    basis = TensorProductBasis(2, polyorder)
    mesh = CutCell.Mesh([0.0, 0.0], [width, width], [nelmts, nelmts], basis)
    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs = CutCell.levelset_coefficients(
        x -> -circle_distance_function(x, center, inradius),
        mesh,
    )

    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    cellquads =
        CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    interfacequads =
        CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    facequads =
        CutCell.FaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition = CutCell.coherent_interface_condition(
        basis,
        interfacequads,
        stiffness,
        cutmesh,
        penalty,
    )

    displacementbc = CutCell.DisplacementCondition(
        analyticalsolution,
        basis,
        facequads,
        stiffness,
        cutmesh,
        x -> onboundary(x, width, width),
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(
        sysmatrix,
        interfacecondition,
        cutmesh,
    )
    CutCell.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        cutmesh,
    )
    CutCell.assemble_coherent_interface_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        displacementbc,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        cutmesh,
        x -> onboundary(x, width, width),
    )

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    sol = matrix \ rhs

    stresserr = compute_stress_error(
        sol,
        basis,
        cellquads,
        stiffness,
        transfstress,
        cutmesh,
        analyticalsolution,
    )

    return stresserr
end

function solve_and_compute_interface_stress_error(
    width,
    center,
    inradius,
    outradius,
    stiffness,
    theta0,
    nelmts,
    polyorder,
    numqp,
    penaltyfactor,
)

    lambda1, mu1 = CutCell.lame_coefficients(stiffness, +1)
    lambda2, mu2 = CutCell.lame_coefficients(stiffness, -1)
    transfstress =
        CutCell.plane_strain_transformation_stress(lambda1, mu1, theta0)

    analyticalsolution = AnalyticalSolution(
        inradius,
        outradius,
        center,
        lambda1,
        mu1,
        lambda2,
        mu2,
        theta0,
    )

    dx = width / nelmts
    meanmoduli = 0.5 * (lambda1 + lambda2 + mu1 + mu2)
    penalty = penaltyfactor / dx * meanmoduli

    basis = TensorProductBasis(2, polyorder)
    mesh = CutCell.Mesh([0.0, 0.0], [width, width], [nelmts, nelmts], basis)
    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs = CutCell.levelset_coefficients(
        x -> -circle_distance_function(x, center, inradius),
        mesh,
    )

    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    cellquads =
        CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    interfacequads =
        CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    facequads =
        CutCell.FaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition = CutCell.coherent_interface_condition(
        basis,
        interfacequads,
        stiffness,
        cutmesh,
        penalty,
    )

    displacementbc = CutCell.DisplacementCondition(
        analyticalsolution,
        basis,
        facequads,
        stiffness,
        cutmesh,
        x -> onboundary(x, width, width),
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(
        sysmatrix,
        interfacecondition,
        cutmesh,
    )
    CutCell.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        cutmesh,
    )
    CutCell.assemble_coherent_interface_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        displacementbc,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        cutmesh,
        x -> onboundary(x, width, width),
    )

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    sol = matrix \ rhs

    perr, nerr = compute_interface_stress_error(
        sol,
        basis,
        interfacequads,
        stiffness,
        transfstress,
        cutmesh,
        analyticalsolution,
    )

    return perr, nerr
end

function mean(v)
    return sum(v) / length(v)
end

function convergence_rate(dx, err)
    return diff(log.(err)) ./ diff(log.(dx))
end

K1, K2 = 247.0, 192.0
mu1, mu2 = 126.0, 87.0
lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)

stiffness = CutCell.HookeStiffness(lambda1, mu1, lambda2, mu2)

theta0 = -0.067

width = 1.0

center = [width / 2, width / 2]
inradius = width / 4
outradius = width

polyorder = 2
numqp = required_quadrature_order(polyorder) + 2
penaltyfactor = 1e2

powers = [3, 4, 5]
nelmts = [2^p + 1 for p in powers]
dx = 1.0 ./ nelmts

stresserr = [
    solve_and_compute_stress_error(
        width,
        center,
        inradius,
        outradius,
        stiffness,
        theta0,
        ne,
        polyorder,
        numqp,
        penaltyfactor,
    ) for ne in nelmts
]

s11err = [st[1] for st in stresserr]
s22err = [st[2] for st in stresserr]
s12err = [st[3] for st in stresserr]
s33err = [st[4] for st in stresserr]

s11rate = convergence_rate(dx, s11err)
s22rate = convergence_rate(dx, s22err)
s12rate = convergence_rate(dx, s12err)
s33rate = convergence_rate(dx, s33err)

@test allapprox(s11rate, repeat([2.0], length(s11rate)), 0.1)
@test allapprox(s22rate, repeat([2.0], length(s11rate)), 0.1)
@test allapprox(s12rate, repeat([2.0], length(s11rate)), 0.1)
@test allapprox(s33rate, repeat([2.0], length(s11rate)), 0.1)

# K1, K2 = 247.0, 192.0
# mu1, mu2 = 126.0, 87.0
# lambda1 = lame_lambda(K1, mu1)
# lambda2 = lame_lambda(K2, mu2)
#
# stiffness = CutCell.HookeStiffness(lambda1, mu1, lambda2, mu2)
#
# theta0 = -0.067
#
# width = 1.0
#
# center = [width / 2, width / 2]
# inradius = width / 4
# outradius = width
#
# polyorder = 2
# numqp = required_quadrature_order(polyorder) + 2
# penaltyfactor = 1e2
#
# powers = [2, 3, 4, 5]
# nelmts = [2^p + 1 for p in powers]
# dx = 1.0 ./ nelmts
#
#
# stresserr = [
#     solve_and_compute_interface_stress_error(
#         width,
#         center,
#         inradius,
#         outradius,
#         stiffness,
#         theta0,
#         ne,
#         polyorder,
#         numqp,
#         penaltyfactor,
#     ) for ne in nelmts
# ]
#
# perr = [st[1] for st in stresserr]
# nerr = [st[2] for st in stresserr]
#
# ps11err = [st[1] for st in perr]
# ps22err = [st[2] for st in perr]
# ps12err = [st[3] for st in perr]
# ps33err = [st[4] for st in perr]
#
# ns11err = [st[1] for st in nerr]
# ns22err = [st[2] for st in nerr]
# ns12err = [st[3] for st in nerr]
# ns33err = [st[4] for st in nerr]
#
# ps11rate = convergence_rate(dx, ps11err)
# ps22rate = convergence_rate(dx, ps22err)
# ps12rate = convergence_rate(dx, ps12err)
# ps33rate = convergence_rate(dx, ps33err)
#
# ns11rate = convergence_rate(dx, ns11err)
# ns22rate = convergence_rate(dx, ns22err)
# ns12rate = convergence_rate(dx, ns12err)
# ns33rate = convergence_rate(dx, ns33err)

# @test allapprox(s11rate, repeat([2.0], length(s11rate)), 0.1)
# @test allapprox(s22rate, repeat([2.0], length(s11rate)), 0.1)
# @test allapprox(s12rate, repeat([2.0], length(s11rate)), 0.1)
# @test allapprox(s33rate, repeat([2.0], length(s11rate)), 0.1)


text:
_TEMPLATE_KWARGS = [:network, :devices, :services]

function _filter_kwargs(kwargs)
    template_kwargs = Dict(kwargs)
    build_kwargs = Dict()
    for kw in setdiff(keys(template_kwargs), _TEMPLATE_KWARGS)
        build_kwargs[kw] = pop!(template_kwargs, kw)
    end
    return template_kwargs, build_kwargs
end

"""
    EconomicDispatchProblem(system::PSY.System; kwargs...)

Creates a `ProblemTemplate` with default DeviceModels for an EconomicDispatch
problem. Uses the template to create an `DecisionProblem`.

# Example

ed_problem = EconomicDispatchProblem(system)

```

# Accepted Key Words
- `network::Type{<:PM.AbstractPowerModel}` : override default network model settings
- `devices::Dict{String, DeviceModel}` : override default `DeviceModel` settings
- `services::Dict{String, ServiceModel}` : override default `ServiceModel` settings
- Key word arguments supported by `DecisionProblem`
```
"""
function EconomicDispatchProblem(system::PSY.System; kwargs...)
    kwargs, problem_kwargs = _filter_kwargs(kwargs)
    output_dir = pop!(problem_kwargs, :output_dir)
    template = template_economic_dispatch(; kwargs...)
    model = DecisionModel(EconomicDispatchProblem, template, system; problem_kwargs...)
    res = build!(model; output_dir=output_dir)
    if res != BuildStatus.BUILT
        error("The EconomicDispatch problem didn't build successfully")
    end
    return model
end

"""
    UnitCommitmentProblem(system::PSY.System; kwargs...)

Creates a `ProblemTemplate` with default DeviceModels for a Unit Commitment
problem. Uses the template to create an `DecisionProblem`.

# Example

uc_problem = UnitCommitmentProblem(system)

```

# Accepted Key Words
- `network::Type{<:PM.AbstractPowerModel}` : override default network model settings
- `devices::Dict{String, DeviceModel}` : override default `DeviceModel` settings
- `services::Dict{String, ServiceModel}` : override default `ServiceModel` settings
- Key word arguments supported by `DecisionProblem`
```
"""
function UnitCommitmentProblem(system::PSY.System; kwargs...)
    kwargs, problem_kwargs = _filter_kwargs(kwargs)
    output_dir = pop!(problem_kwargs, :output_dir)
    template = template_unit_commitment(; kwargs...)
    model = DecisionModel(UnitCommitmentProblem, template, system; problem_kwargs...)
    res = build!(model; output_dir=output_dir)
    if res != BuildStatus.BUILT
        error("The EconomicDispatch problem didn't build successfully")
    end
    return model
end

"""
    AGCReserveDeployment(system::PSY.System; kwargs...)

Creates a `ProblemTemplate` with default DeviceModels for an AGC Reserve Deplyoment Problem.
Uses the template to create an `DecisionProblem`.

# Example

agc_problem = AGCReserveDeployment(system)

```

# Accepted Key Words
- Key word arguments supported by `DecisionProblem`
```
"""
function AGCReserveDeployment(system::PSY.System; kwargs...)
    kwargs, problem_kwargs = _filter_kwargs(kwargs)
    output_dir = pop!(problem_kwargs, :output_dir)
    template = template_agc_reserve_deployment(; kwargs...)
    model = DecisionModel(UnitCommitmentProblem, template, system; problem_kwargs...)
    res = build!(model; output_dir=output_dir)
    if res != BuildStatus.BUILT
        error("The EconomicDispatch problem didn't build successfully")
    end
    return model
end

"""
    run_unit_commitment(system::PSY.System; kwargs...)

Creates a `ProblemTemplate` with default DeviceModels for a Unit Commitment
problem. Uses the template to create an `DecisionProblem`. Solves the created operations problem.

# Example
results = run_unit_commitment(system; optimizer = optimizer)
```

# Accepted Key Words
- `network::Type{<:PM.AbstractPowerModel}` : override default network model settings
- `devices::Dict{String, DeviceModel}` : override default `DeviceModel` settings
- `services::Dict{String, ServiceModel}` : override default `ServiceModel` settings
- `optimizer::JuMP Optimizer` : An optimizer is a required key word
- `output_dir::AbstractString`  : Path to save outputs
- Key word arguments supported by `DecisionProblem`
"""

function run_unit_commitment(sys::PSY.System; kwargs...)
    model = UnitCommitmentProblem(sys; kwargs...)
    solve_status = solve!(model)
    return solve_status
end

"""
    run_economic_dispatch(system::PSY.System; kwargs...)

Creates a `ProblemTemplate` with default DeviceModels for an EconomicDispatch
problem. Uses the template to create an `DecisionProblem`.

# Example

results = run_economic_dispatch(system; optimizer = optimizer)

```

# Accepted Key Words
- `network::Type{<:PM.AbstractPowerModel}` : override default network model settings
- `devices::Dict{String, DeviceModel}` : override default `DeviceModel` settings
- `services::Dict{String, ServiceModel}` : override default `ServiceModel` settings
- `optimizer::JuMP optimizer` : a JuMP optimizer is a required key word
- `output_dir::AbstractString`  : Path to save outputs
- Key word arguments supported by `DecisionProblem`
```
"""
function run_economic_dispatch(sys::PSY.System; kwargs...)
    model = EconomicDispatchProblem(sys; kwargs...)
    solve_status = solve!(model)
    return solve_status
end

ext: jl


text:
using DiffEqBayes, CmdStan, DynamicHMC, DataFrames
using Distributions, BenchmarkTools
using OrdinaryDiffEq, RecursiveArrayTools, ParameterizedFunctions
using ModelingToolkit
using StatsPlots, CSV
gr(size=(600,900))

ProjDir = @__DIR__
cd(ProjDir)
isdir("tmp") && rm("tmp", recursive=true)

include("../stan_inference.jl")

Random.seed!(123)

df = CSV.read("$(ProjDir)/../lynx_hare.csv", DataFrame; delim=",")

f = @ode_def LotkaVolterraTest begin
    dx = a*x - b*x*y
    dy = -c*y + d*x*y
end a b c d

u0 = [30.0, 4.0]
tspan = (0.0, 21.0)
p = [0.55, 0.028, 0.84, 0.026]

prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Tsit5())

t = collect(range(1,stop=20,length=20))

sig = 0.3
data = convert(Array, VectorOfArray([(sol(t[i]) + sig*randn(2)) for i in 1:length(t)]))

#data = transpose(hcat(df[2:end, :Hare], df[2:end, :Lynx]))

scatter(t, data[1,:], lab="#prey (data)")
scatter!(t, data[2,:], lab="#predator (data)")
plot!(sol)
savefig("$(ProjDir)/fig_01.png")

priors = [Truncated(Normal(1,0.5),0,3),Truncated(Normal(0.05,0.05),0,2),
  Truncated(Normal(1,0.5),0,3),Truncated(Normal(0.05,0.05),0,2)]

# Stan.jl backend

#The solution converges for tolerance values lower than 1e-3, lower tolerance
#leads to better accuracy in result but is accompanied by longer warmup and
#sampling time, truncated normal priors are used for preventing Stan from stepping into negative values.

nchains = 4
num_samples = 800
@time bayesian_result_stan = stan_inference(prob,t,data,priors,
  num_samples=num_samples, printsummary=false, nchains=nchains);
sdf1  = CmdStan.read_summary(bayesian_result_stan.model)
println()

@time bayesian_result_stan = stan_inference(prob,t,data,
  priors, bayesian_result_stan.model,
  num_samples=num_samples, printsummary=false, nchains=nchains);
sdf2  = CmdStan.read_summary(bayesian_result_stan.model)
println()

cnames = ["theta1", "theta2", "theta3", "theta4", "sigma1", "sigma2", ]

df = CmdStan.convert_a3d(bayesian_result_stan.chains, bayesian_result_stan.cnames,
  Val(:dataframes));
df_stan = append!(append!(df[1], df[2]), append!(df[3], df[4]));
df_stan = df_stan[:, [10, 11, 12, 13, 8, 9]]
rename!(df_stan, Symbol.(cnames))

sdf1 |> display
sdf2 |> display

# DynamicHMC.jl backend

@time bayesian_result_dynamichmc = dynamichmc_inference(prob,Tsit5(),t,data,priors,
  num_samples=Int(4*num_samples));

N = length(bayesian_result_dynamichmc.posterior)

res_array = zeros(N, 6);
for i in 1:length(bayesian_result_dynamichmc.posterior)
  res_array[i,:] = vcat(
      bayesian_result_dynamichmc.posterior[i].parameters,
      bayesian_result_dynamichmc.posterior[i].σ
    )
end

df_dhmc = DataFrame(
  :theta1 => res_array[:, 1],
  :theta2 => res_array[:, 2],
  :theta3 => res_array[:, 3],
  :theta4 => res_array[:, 4],
  :sigma1 => res_array[:, 5],
  :sigma2 => res_array[:, 6]
)

println("\nDhmc results\n")
mean(res_array, dims=1) |> display
std(res_array, dims=1) |> display
println()

function plot_theta_results(
    dfs::DataFrame,
    dfd::DataFrame,
    varn::AbstractString,
    idx::Int,
    priors=priors
  )
  prior_theta = Symbol("prior_"*varn*string(idx))
  theta = Symbol(varn*string(idx))
  prior_theta_array = rand(priors[idx], size(dfs, 1))
  p1 = plot(prior_theta_array, xlab="prior_theta_$idx", leg=false)
  p2 = density(prior_theta_array, xlab="prior_theta_$idx", leg=false)
  p3 = plot(dfs[:, theta], xlab="theta_$idx", lab="stan")
  plot!(dfs[:, theta], lab="dhmc")
  p4 = density(dfs[:, theta], xlab="theta_$idx", lab="stan")
  density!(dfd[:, theta], lab="dhmc")
  p5 = density(dfs[:, :sigma1], xlab="sigma1", lab="stan")
  density!(dfd[:, :sigma1], lab="dhmc")
  p6 = density(dfs[:, :sigma2], xlab="sigma2", lab="stan")
  density!(dfd[:, :sigma2], lab="dhmc")

  plot(p1, p2, p3, p4, p5, p6, layout=(3,2))
  savefig("$(ProjDir)/$(varn)_$idx")

end

for i in 1:4
  plot_theta_results(df_stan, df_dhmc, "theta", i, priors)
end

stan_results = "
                mean se_mean    sd    10%    50%    90% n_eff Rhat
    theta[1]   0.549   0.002 0.065  0.469  0.545  0.636  1163    1
    theta[2]   0.028   0.000 0.004  0.023  0.028  0.034  1281    1
    theta[3]   0.797   0.003 0.091  0.684  0.791  0.918  1125    1
    theta[4]   0.024   0.000 0.004  0.020  0.024  0.029  1170    1
    sigma[1]   0.248   0.001 0.045  0.198  0.241  0.306  2625    1
    sigma[2]   0.252   0.001 0.044  0.201  0.246  0.310  2808    1
    z_init[1] 33.960   0.056 2.909 30.363 33.871 37.649  2684    1
    z_init[2]  5.949   0.011 0.533  5.294  5.926  6.644  2235    1
  ";


text:
lazy(f; kw...) = delayed(f; kw...)

# If any input is lazy, make the output lazy
maybe_lazy(f, x) = any(x->x isa Union{Thunk, Chunk}, x) ? lazy(f)(x...) : f(x...)
maybe_lazy(f) = (x...) -> maybe_lazy(f, x)

"""
    compute(tree::FileTree; cache=true)

Compute any lazy values (Thunks) in `tree` and return a new tree where the values refer to
the computed values (maybe on remote processes). The tree still behaves as a Lazy tree. `exec` on it will fetch the values from remote processes.
"""
compute(d::FileTree; cache=true, kw...) = compute(Dagger.Context(), d; cache=cache, kw...)
function compute(ctx, d::FileTree; cache=true, kw...)
    thunks = []
    mapvalues(d; lazy=false) do x
        if x isa Thunk
            if cache
                x.cache = true
            end
            push!(thunks, x)
        end
    end

    vals = compute(delayed((xs...)->[xs...]; meta=true)(thunks...); kw...)

    i = 0
    mapvalues(d; lazy=false) do x
        i += 1
        vals[i]
    end
end

exec(d::FileTree) = mapvalues(exec, compute(d, cache=false); lazy=false)
exec(d::Union{Thunk, Chunk}) = collect(compute(d))
"""
    exec(x)

If `x` is a FileTree, computes any uncomputed `Thunk`s stored as values in it. Returns a new tree with the computed values.
If `x` is a `Thunk` (such as the result of a `reducevalues`), then exec will compute the result.
If `x` is anything else, `exec` just returns the same value.
"""
exec(x) = x


text:
module NamedTuples

"""
    `function Core.NamedTuple(kwargs::Dict{Symbol,T})::NamedTuple where {T}`

Converts the `Dict` `kwargs` with keys of type `Symbol` to a `NamedTuple`.

### Example

```julia
julia> NamedTuple(Dict(:a => "a", :b => "b"))
(a = "a", b = "b")
```
"""
function Core.NamedTuple(kwargs::Dict{Symbol,T})::NamedTuple where {T}
  NamedTuple{Tuple(keys(kwargs))}(collect(values(kwargs)))
end

"""
    `function Core.NamedTuple(kwargs::Dict{Symbol,T}, property::Symbol, value::String)::NamedTuple where {T}`

Prepends `value` to `kwargs[property]` if defined or adds a new `kwargs[property] = value` and then converts the
resulting `kwargs` dict to a `NamedTuple`.

### Example

```julia
julia> NamedTuple(Dict(:a => "a", :b => "b"), :d, "h")
(a = "a", b = "b", d = "h")

julia> NamedTuple(Dict(:a => "a", :b => "b"), :a, "h")
(a = "h a", b = "b")
```
"""
function Core.NamedTuple(kwargs::Dict{Symbol,T}, property::Symbol, value::String)::NamedTuple where {T}
  value = "$value $(get!(kwargs, property, ""))" |> strip
  kwargs = delete!(kwargs, property)
  kwargs[property] = value

  NamedTuple(kwargs)
end

end


text:
using FBD
using Test

include("util.jl")

include("dchsbm/model.jl")

include("hyper_pa.jl")

include("kronecker/util.jl")
include("kronecker/model.jl")

include("typing.jl")

include("ER.jl")

include("io.jl")

include("examples.jl")

max_stars_repo_name: LilithHafner/FBD

max_issues_repo_name: LilithHafner/FBD

max_forks_repo_name: LilithHafner/FBD


text:
"""
    @hprofile

A macro which
- resets the default `TimerOutputs.get_defaulttimer` to zero
- execute the code block
- print the profiling details

This is useful as a coarse-grained profiling strategy in `HMatrices`
to get a rough idea of where time is spent. Note that this relies on
`TimerOutputs` annotations manually inserted in the code.
"""
macro hprofile(block)
    return quote
        TimerOutputs.enable_debug_timings(HMatrices)
        reset_timer!()
        $(esc(block))
        print_timer()
    end
end

"""
    PermutedMatrix{K,T} <: AbstractMatrix{T}

Structured used to reprensent the permutation of a matrix-like object. The
original matrix is stored in the `data::K` field, and the permutations are
stored in `rowperm` and `colperm`.
"""
struct PermutedMatrix{K,T} <: AbstractMatrix{T}
    data::K # original matrix
    rowperm::Vector{Int}
    colperm::Vector{Int}
    function PermutedMatrix(orig,rowperm,colperm)
        K = typeof(orig)
        T = eltype(orig)
        new{K,T}(orig,rowperm,colperm)
    end
end
Base.size(M::PermutedMatrix) = size(M.data)

function Base.getindex(M::PermutedMatrix,i,j)
    ip = M.rowperm[i]
    jp = M.colperm[j]
    M.data[ip,jp]
end

"""
    hilbert_cartesian_to_linear(n,x,y)

Convert the cartesian indices `x,y` into a linear index `d` using a hilbert
curve of order `n`. The coordinates `x,y` range from `0` to `n-1`, and the
output `d` ranges from `0` to `n^2-1`.

See [https://en.wikipedia.org/wiki/Hilbert_curve](https://en.wikipedia.org/wiki/Hilbert_curve).
"""
function hilbert_cartesian_to_linear(n::Integer,x,y)
    @assert ispow2(n)
    @assert 0 ≤ x ≤ n-1
    @assert 0 ≤ y ≤ n-1
    d = 0
    s = n >> 1
    while s>0
        rx = (x & s) > 0
        ry = (y & s) > 0
        d += s^2*((3*rx)⊻ry)
        x,y = _rot(n,x,y,rx,ry)
        s = s >> 1
    end
    @assert 0 ≤ d ≤ n^2-1
    return d
end

"""
    hilbert_linear_to_cartesian(n,d)

Convert the linear index `0 ≤ d ≤ n^2-1` into the cartesian coordinates `0 ≤ x <
n-1` and `0 ≤ y ≤ n-1` on the Hilbert curve of order `n`.

See [https://en.wikipedia.org/wiki/Hilbert_curve](https://en.wikipedia.org/wiki/Hilbert_curve).
"""
function hilbert_linear_to_cartesian(n::Integer,d)
    @assert ispow2(n)
    @assert 0 ≤ d ≤ n^2-1
    x,y = 0,0
    s = 1
    while s<n
        rx = 1 & (d >> 1)
        ry = 1 & (d ⊻ rx)
        x,y = _rot(s,x,y,rx,ry)
        x +=  s*rx
        y +=  s*ry
        d = d >> 2
        s = s << 1
    end
    @assert 0 ≤ x ≤ n-1
    @assert 0 ≤ y ≤ n-1
    return x,y
end

# auxiliary function using in hilbert curve. Rotates the points x,y
function _rot(n,x,y,rx,ry)
    if ry == 0
        if rx == 1
            x = n-1 - x
            y = n-1 - y
        end
        x,y = y,x
    end
    return x,y
end

function hilbert_points(n::Integer)
    @assert ispow2(n)
    xx = Int[]
    yy = Int[]
    for d in 0:n^2-1
        x,y = hilbert_linear_to_cartesian(n,d)
        push!(xx,x)
        push!(yy,y)
    end
    xx,yy
end

"""
    build_sequence_partition(seq,nq,cost,nmax)

Partition the sequence `seq` into `nq` contiguous subsequences with a maximum of
cost of `nmax` per set. Note that if `nmax` is too small, this may not be
possible (see [`has_partition`](@ref)).
"""
function build_sequence_partition(seq,np,cost,cmax)
    acc = 0
    partition = [empty(seq) for _ in 1:np]
    k = 1
    for el in seq
        c = cost(el)
        acc += c
        if acc > cmax
            k += 1
            push!(partition[k],el)
            acc = c
            @assert k <= np "unable to build sequence partition. Value of  `cmax` may be too small."
        else
            push!(partition[k],el)
        end
    end
    return partition
end

"""
    find_optimal_cost(seq,nq,cost,tol)

Find an approximation to the cost of an optimal partitioning of `seq` into `nq`
contiguous segments. The optimal cost is the smallest number `cmax` such that
`has_partition(seq,nq,cost,cmax)` returns `true`.
"""
function find_optimal_cost(seq,np,cost=identity,tol=1)
    lbound = maximum(cost,seq) |> Float64
    ubound = sum(cost,seq)     |> Float64
    guess  = (lbound+ubound)/2
    while ubound - lbound ≥ tol
        if has_partition(seq,np,guess,cost)
            ubound = guess
        else
            lbound = guess
        end
        guess  = (lbound+ubound)/2
    end
    return ubound
end

"""
    find_optimal_partition(seq,nq,cost,tol=1)

Find an approximation to the optimal partition `seq` into `nq` contiguous
segments according to the `cost` function. The optimal partition is the one
which minimizes the maximum `cost` over all possible partitions of `seq` into
`nq` segments.

The generated partition is optimal up to a tolerance `tol`; for integer valued
`cost`, setting `tol=1` means the partition is optimal.
"""
function find_optimal_partition(seq,np,cost=(x)->1,tol=1)
    cmax = find_optimal_cost(seq,np,cost,tol)
    p = build_sequence_partition(seq,np,cost,cmax)
    return p
end

"""
    has_partition(v,np,cost,cmax)

Given a vector `v`, determine whether or not a partition into `np` segments is
possible where the `cost` of each partition does not exceed `cmax`.
"""
function has_partition(seq,np,cmax,cost=identity)
    acc = 0
    k   = 1
    for el in seq
        c = cost(el)
        acc += c
        if acc > cmax
            k += 1
            acc = c
            k > np && (return false)
        end
    end
    return true
end

"""
    disable_getindex()

Call this function to disable the `getindex` method on `AbstractHMatrix`. This
is useful to avoid performance pitfalls associated with linear algebra methods
falling back to a generic implementation which uses the `getindex` method.
Calling `getindex(H,i,j)` will error after calling this function.
"""
disable_getindex() = (ALLOW_GETINDEX[] = false)

"""
    enable_getindex()

The opposite of [`disable_getindex`](@ref).
"""
enable_getindex()  = (ALLOW_GETINDEX[] = true)

max_stars_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
using Mill, JSON, Flux, JsonGrinder, Test

using JsonGrinder: DictEntry, suggestextractor, schema
using Mill: reflectinmodel

@testset "testing pipeline mixing different types of arrays" begin
	j1 = JSON.parse("""{"a": 1, "b": "hello works", "c":{ "a":1 ,"b": "hello world"}}""")
	j2 = JSON.parse("""{"a": 2, "b": "hello world", "c":{ "a":2 ,"b": "hello"}}""")

	sch = schema([j1,j2])
	extractor = suggestextractor(sch)
	ds = map(s-> extractor(s), [j1,j2])
	dss = reduce(catobs, ds)

	m = reflectinmodel(dss, k -> Dense(k,10, relu));
	o = m(dss).data
	for i in 1:length(ds)
		@test o[:,i] ≈ m(ds[i]).data
	end
end

@testset "testing pipeline with simple arrays and missing values" begin
	j1 = JSON.parse("""{"a": 4, "b": {"a":[1,2,3],"b": 1},"c": { "a": {"a":[1,2,3],"b":[4,5,6]}}}""",inttype=Float64)
	j2 = JSON.parse("""{"a": 4, "c": { "a": {"a":[2,3],"b":[5,6]}}}""")
	j3 = JSON.parse("""{"a": 4, "b": {"a":[1,2,3],"b": 1}}""")
	j4 = JSON.parse("""{"a": 4, "b": {}}""")
	j5 = JSON.parse("""{"b": {}}""")
	j6 = JSON.parse("""{}""")

	sch = JsonGrinder.schema([j1,j2,j3])
	extractor = suggestextractor(sch)
	dss = map(s-> extractor(s), [j1,j2,j3,j4,j5,j6])
	ds = reduce(catobs, dss)
	m = reflectinmodel(ds, k -> Dense(k,10, relu));
	o = m(ds).data

	for i in 1:length(dss)
		@test o[:,i] ≈ m(dss[i]).data
	end
end

@testset "testing pipeline with arrays of dicts and missing values" begin
	j1 = JSON.parse("""{"a": [{"a":1},{"b":2}]}""")
	j2 = JSON.parse("""{"a": [{"a":1,"b":3},{"b":2,"a" : 1}]}""")
	j3 = JSON.parse("""{"a": [{"a":2,"b":3}]}""")
	j4 = JSON.parse("""{"a": []}""")
	j5 = JSON.parse("""{}""")

	sch = JsonGrinder.schema([j1,j2,j3])
	extractor = suggestextractor(sch)
	dss = map(s-> extractor(s), [j1,j2,j3,j4,j5])
	ds = reduce(catobs, dss)
	m = reflectinmodel(ds, k -> Dense(k,10, relu));
	o = m(ds).data
	for i in 1:length(dss)
		@test o[:,i] ≈ m(dss[i]).data
	end
end


text:
export DiGraph, add_edge!, remove_edge!, nodes, edges, query

struct DiGraph
    edges::Set{Pair{}}

    DiGraph() = new(Set{Pair}())
end

function DiGraph(graph::DiGraph)
    new_graph = DiGraph()
    for e in edges(graph)
        add_edge!(new_graph, e.first => e.second)
    end
    new_graph
end

"""
    add_edge!(graph, from, to)
    add_edge!(::DiGraph, ::Pair)
Add an edge to the DiGraph going between the specified nodes.
The graph is returned.
"""
function add_edge! end

function add_edge!(graph, from, to)
    add_edge!(graph, from => to)
    graph
end

function add_edge!(graph::DiGraph, edge::Pair)
    push!(graph.edges, edge)
    graph
end

"""
    remove_edge!(graph, from, to)
    remove_edge!(::DiGraph, ::Pair)
Remove the edge of the DiGraph going between the specified nodes.
The graph is returned.
"""
function remove_edge! end

function remove_edge!(graph::DiGraph, edge::Pair)::DiGraph
    delete!(graph.edges, edge)
    graph
end

function remove_edge!(graph, from, to)
    remove_edge!(graph, from => to)
    graph
end

# Dict interface:
Base.keys(graph::DiGraph) = Set([p.first for p in graph.edges])
Base.values(graph::DiGraph) = Set([p.second for p in graph.edges])

function Base.haskey(graph::DiGraph, key)::Bool
    !isempty(filter(edge -> edge.first == key, graph.edges))
end

function Base.getindex(graph::DiGraph, key)::Set
    Set((p -> p.second).(filter(p -> p.first == key, graph.edges)))
end

nodes(graph::DiGraph) = union(keys(graph), values(graph))

edges(graph::DiGraph) = graph.edges

"""
    query(::DiGraph, from, to)
Return the `Set` of edges that match the query.
`from` and `to` can be objects or types, so, for example,
`query(graph, Number, :a)` returns all of the edges that
go from any `Number` to the `Symbol` `:a`.
"""
function query(graph::DiGraph, from, to)::Set
    function querytest(elt, val)
        if val isa Type
            elt isa val
        else
            elt == val
        end
    end
    filter(p -> querytest(p.first, from) && querytest(p.second, to),
           graph.edges)
end


text:
using Pkg
pkg"activate Burner"
using Revise
using Burner
using Images
using VideoIO

width = 512
height = 288

demodisp = Burner.Display(ones(3, height, width), 0.1)
imageburn = imresize(load("17498.jpg"), (height, width))
imageshow = Float64.(Images.channelview(imresize(load("descargar-images-full-hd-1080p.jpg"), (height, width))))

for i = 1:50
    rotimg = RGB.(Burner.permuteImage(imageburn, 0))
    Burner.evolve(demodisp, channelview(rotimg), 0.1)
end
imagepresented = colorview(RGB, Burner.present(demodisp, imageshow))

max_stars_repo_name: thekineticeffect/burninsim

max_issues_repo_name: thekineticeffect/burninsim

max_forks_repo_name: thekineticeffect/burninsim


text:
#=
timeseries decompositions, i.e. separating some data into components
=#
export seasonal_decomposition
"""
    seasonal_decomposition(A::ClimArray, fs = [1, 2]) → seasonal, residual
Decompose `A` into a seasonal and residual components, where the seasonal contains the
periodic parts of `A`, with frequencies given in `fs`, and residual contains what's left.

`fs` is measured in 1/year. This function works even for non-equispaced time axis (e.g.
monthly averages) and uses LPVSpectral.jl and SignalDecomposition.jl.
"""
seasonal_decomposition(A::AbDimArray, b) = seasonal_decomposition(dims(A, Time), A, b)

function seasonal_decomposition(t, A::AbDimArray, fs::Vector)
    @assert hasdim(A, Time)
    E = _numbertype(T)
    method = Sinusoidal(E.(fs ./ DAYS_IN_YEAR))
    seasonal = DimensionalData.basetypeof(A)(copy(Array(A)), dims(A); name = A.name*"seasonal")
    residual = DimensionalData.basetypeof(A)(copy(Array(A)), dims(A); name = A.name*"residual")

    # TODO: Fix this to use "real time"
    truetime = Float32.(cumsum(daysinmonth.(t)))
    for i in otheridxs(T, Time)
        y = Array(A[i...])
        sea, res = SignalDecomposition.decompose(truetime, y, method)
        seasonal[i...] .= sea
        residual[i...] .= res
    end
    return seasonal, residual
end


text:
# Maximum likelihood estimation from a given initial point
using EstimateSDE
using Random
using Printf
using Distributions

## Settings
proc = testprocess(:C)
N = 200
Tav = 0.5
Tkld = 0.5
engine = SdeRoots

data = rand(MersenneTwister(12202),proc,N,Tav)
problem = ProblemLike(data;p=EstimateSDE.order(proc))
tor = MlEstimator()
procstart = transform_ss(problem,proc)
proc_ml,opt = ml_estimate_init(:dynamichmc,tor,problem,procstart)

D = kldiv(proc,proc_ml,(0:Tkld:N*Tav))

## Report
printstyled("Actual:\n",color=:blue)
describe_model(proc)
printstyled("Estimate:\n",color=:blue)
describe_model(proc_ml)
@printf "KL divergence: %0.1f\n" D


text:
using Test
using QuakeML

@testset "Comparison" begin
    # Two identical EventParameters
    datafile = joinpath(@__DIR__, "data", "2004-12-26_mag5+.qml")
    qml = QuakeML.read(datafile)
    qml′ = QuakeML.read(datafile)
    # A different EventParameters: one field is missing in one and set in another
    qml″ = deepcopy(qml)
    qml″.event[end].magnitude[1].mag.uncertainty = 0.1

    @testset "==" begin
        @test qml == qml′
        @test qml != qml″
    end

    @testset "hash" begin
        @test hash(qml) == hash(qml′)
        @test hash(qml) != hash(qml″)
    end

    # These should follow from the above
    @testset "Others" begin
        @test unique([qml, qml′]) == [qml]
        @test unique([qml, qml″]) == [qml, qml″]
        @test Dict(qml=>1) == Dict(qml′=>1) != Dict(qml″=>1)
    end
end

max_stars_repo_name: anowacki/QuakeML


text:
using Documenter, DataFramesMeta

makedocs(
	modules = [DataFramesMeta],
	sitename = "DataFramesMeta Documentation",
	format = Documenter.HTML(
		canonical = "https://juliadata.github.io/DataFramesMeta.jl/stable/"
	),
	pages = Any[
		"Introduction" => "index.md",
		"Tutorial for coming from dplyr" => "dplyr.md",
		"API" => "api/api.md"])

deploydocs(
    repo = "github.com/JuliaData/DataFramesMeta.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
using Test;

active = '#'
inactive = '.'

struct Point
    x::Int
    y::Int
    z::Int
    w::Int
end

function get_neighbors(p::Point)
    neighbors = []
    for x_offset in [1,0,-1]
        for y_offset in [1,0,-1]
            for z_offset in [1,0,-1]
                for w_offset in [1,0,-1]
                    if x_offset == 0 && y_offset == 0 && z_offset == 0 && w_offset == 0
                        # skip returning the same point
                        continue
                    end
                    push!(neighbors, Point(p.x + x_offset, p.y + y_offset, p.z + z_offset, p.w + w_offset))
                end
            end
        end
    end
    return neighbors
end

function is_active(values, pt)
    return in(pt, keys(values)) && values[pt] == active
end

function copy(p::Point)
    return Point(p.x, p.y, p.z, p.w)
end

function do_iter(values, low::Point, high::Point)
    new_values = Dict()
    # new_low = copy(low)
    # new_high = copy(high)

    # We need not just the current points, but also all things adjacent to them
    all_points = Set()
    for p in keys(values)
        for n in get_neighbors(p)
            push!(all_points, n)
        end
    end

    for p in all_points
        active_neighbor_count = sum(map((x) -> is_active(values, x), get_neighbors(p)))
        # If a cube is active
        if is_active(values, p)
            # and exactly 2 or 3 of its neighbors are also active, the cube remains active.
            # NOTE: This was my bug :( that slowed me down a lot.. I had swapped
            # arguments and written in([2,3], active_neighbor_count)
            if in(active_neighbor_count, [2,3])
                new_values[p] = active
            # Otherwise, the cube becomes inactive.
            else
                new_values[p] = inactive
            end
        # If a cube is inactive
        else
            # but exactly 3 of its neighbors are active, the cube becomes active.
            if active_neighbor_count == 3
                new_values[p] = active
            # Otherwise, the cube remains inactive.
            else
                new_values[p] = inactive
            end
        end

        # TODO: Try to have better boundaries
        # if new_values[p] == active
        #     # update boundary values
        #     new_low.x = min(p.x, new_low.x)
        #     new_low.y = min(p.y, new_low.y)
        #     new_low.z = min(p.z, new_low.z)
        #     new_high.x = max(p.x, new_high.x)
        #     new_high.y = max(p.y, new_high.y)
        #     new_high.z = max(p.z, new_high.z)
        # end
    end

    new_low = Point(low.x - 1, low.y - 1, low.z - 1, low.w - 1)
    new_high = Point(high.x + 1, high.y + 1, high.z + 1, high.w + 1)
    return new_values, new_low, new_high
end

DEBUG = false
function print_values(values, low, high)
    if ! DEBUG
        return
    end

    for z in low.z:high.z
        println("z=", z)
        for y in low.y:high.y
            for x in low.x:high.x
                print(values[Point(x, y, z, w)])
            end
            println("")
        end
        println("")
    end
end

function count_active(values)
    result = 0
    for k in keys(values)
        v = values[k]
        if v == active
            result += 1
        end
    end
    return result
end

function run(fname, cycles=6)
    # parse input
    f = open(fname, "r")
    lines = readlines(f)
    close(f)

    # initialize
    values = Dict{Point,Char}()
    for y in 1:length(lines)
        for (x, val) in enumerate(lines[y])
            values[Point(x, y, 0, 0)] = val
        end
    end

    low_point = Point(1, 1, 0, 0)
    high_point = Point(3, 3, 0, 0)

    # print base state
    println("After ", 0, " cycle...")
    print_values(values, low_point, high_point)

    # do cycles...
    for c in 1:cycles
        # println(length(values), low_point, high_point)
        values, low_point, high_point = do_iter(values, low_point, high_point)
        println("After ", c, " cycle...")
        print_values(values, low_point, high_point)
        println("")
    end

    # Count number of active (`#`)
    result = count_active(values)
    return result
end

# @test run("17ex.txt", 0) == 5
# @test run("17ex.txt", 1) == 11
# @test run("17ex.txt", 2) == 21

# @test run("17ex.txt") == 112

# print(run("17.txt"))

@test run("17ex.txt") == 848

println(run("17.txt"))


text:
# MathOptInterface results
MOI.get(model::Optimizer, ::MOI.RawStatusString) = string(get_status(model.inner))

KN_TO_MOI_RETURN_STATUS = Dict{Int,MOI.TerminationStatusCode}(
    0 => MOI.LOCALLY_SOLVED,
    -100 => MOI.ALMOST_OPTIMAL,
    -101 => MOI.SLOW_PROGRESS,
    -102 => MOI.SLOW_PROGRESS,
    -103 => MOI.SLOW_PROGRESS,
    -200 => MOI.LOCALLY_INFEASIBLE,
    -201 => MOI.LOCALLY_INFEASIBLE,
    -202 => MOI.LOCALLY_INFEASIBLE,
    -203 => MOI.LOCALLY_INFEASIBLE,
    -204 => MOI.LOCALLY_INFEASIBLE,
    -205 => MOI.LOCALLY_INFEASIBLE,
    -300 => MOI.DUAL_INFEASIBLE,
    -301 => MOI.DUAL_INFEASIBLE,
    -400 => MOI.ITERATION_LIMIT,
    -401 => MOI.TIME_LIMIT,
    -402 => MOI.OTHER_LIMIT,
    -403 => MOI.OTHER_LIMIT,
    -404 => MOI.OTHER_LIMIT,
    -405 => MOI.OTHER_LIMIT,
    -406 => MOI.NODE_LIMIT,
    -410 => MOI.ITERATION_LIMIT,
    -411 => MOI.TIME_LIMIT,
    -412 => MOI.INFEASIBLE,
    -413 => MOI.INFEASIBLE,
    -414 => MOI.OTHER_LIMIT,
    -415 => MOI.OTHER_LIMIT,
    -416 => MOI.NODE_LIMIT,
    -500 => MOI.INVALID_MODEL,
    -501 => MOI.NUMERICAL_ERROR,
    -502 => MOI.INVALID_MODEL,
    -503 => MOI.MEMORY_LIMIT,
    -504 => MOI.INTERRUPTED,
    -505 => MOI.OTHER_ERROR,
    -506 => MOI.OTHER_ERROR,
    -507 => MOI.OTHER_ERROR,
    -508 => MOI.OTHER_ERROR,
    -509 => MOI.OTHER_ERROR,
    -510 => MOI.OTHER_ERROR,
    -511 => MOI.OTHER_ERROR,
    -512 => MOI.OTHER_ERROR,
    -513 => MOI.OTHER_ERROR,
    -514 => MOI.OTHER_ERROR,
    -515 => MOI.OTHER_ERROR,
    -516 => MOI.OTHER_ERROR,
    -517 => MOI.OTHER_ERROR,
    -518 => MOI.OTHER_ERROR,
    -519 => MOI.OTHER_ERROR,
    -519 => MOI.OTHER_ERROR,
    -520 => MOI.OTHER_ERROR,
    -521 => MOI.OTHER_ERROR,
    -522 => MOI.OTHER_ERROR,
    -523 => MOI.OTHER_ERROR,
    -524 => MOI.OTHER_ERROR,
    -525 => MOI.OTHER_ERROR,
    -526 => MOI.OTHER_ERROR,
    -527 => MOI.OTHER_ERROR,
    -528 => MOI.OTHER_ERROR,
    -529 => MOI.OTHER_ERROR,
    -530 => MOI.OTHER_ERROR,
    -531 => MOI.OTHER_ERROR,
    -532 => MOI.OTHER_ERROR,
    -600 => MOI.OTHER_ERROR,
)

# Refer to KNITRO manual for solver status:
# https://www.artelys.com/tools/knitro_doc/3_referenceManual/returnCodes.html#returncodes
function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    if model.number_solved == 0
        return MOI.OPTIMIZE_NOT_CALLED
    end
    status = get_status(model.inner)
    if haskey(KN_TO_MOI_RETURN_STATUS, status)
        return KN_TO_MOI_RETURN_STATUS[status]
    else
        error("Unrecognized KNITRO status $status")
    end
end

# TODO
function MOI.get(model::Optimizer, ::MOI.ResultCount)
    return (model.number_solved >= 1) ? 1 : 0
end

function MOI.get(model::Optimizer, status::MOI.PrimalStatus)
    if (
        (model.number_solved == 0) ||   # no solution available if the model is unsolved
        (status.result_index > 1)       # KNITRO stores only a single solution
    )
        return MOI.NO_SOLUTION
    end
    status = get_status(model.inner)
    if status == 0
        return MOI.FEASIBLE_POINT
    elseif -109 <= status <= -100
        return MOI.FEASIBLE_POINT
    elseif -209 <= status <= -200
        return MOI.INFEASIBLE_POINT
    elseif status == -300
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif -409 <= status <= -400
        return MOI.FEASIBLE_POINT
    elseif -419 <= status <= -410
        return MOI.INFEASIBLE_POINT
    elseif -599 <= status <= -500
        return MOI.UNKNOWN_RESULT_STATUS
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

function MOI.get(model::Optimizer, status::MOI.DualStatus)
    if (
        (model.number_solved == 0) ||   # no solution available if the model is unsolved
        (status.result_index > 1)       # KNITRO stores only a single solution
    )
        return MOI.NO_SOLUTION
    end
    status = get_status(model.inner)
    if status == 0
        return MOI.FEASIBLE_POINT
    elseif -109 <= status <= -100
        return MOI.FEASIBLE_POINT
    elseif -209 <= status <= -200
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif status == -300
        return MOI.NO_SOLUTION
    elseif -409 <= status <= -400
        return MOI.FEASIBLE_POINT
    elseif -419 <= status <= -410
        return MOI.INFEASIBLE_POINT
    elseif -599 <= status <= -500
        return MOI.UNKNOWN_RESULT_STATUS
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

function MOI.get(model::Optimizer, obj::MOI.ObjectiveValue)
    if model.number_solved == 0
        error("ObjectiveValue not available.")
    elseif obj.result_index > 1
        throw(MOI.ResultIndexBoundsError{MOI.ObjectiveValue}(obj, 1))
    end
    return get_objective(model.inner)
end

function MOI.get(model::Optimizer, v::MOI.VariablePrimal, vi::MOI.VariableIndex)
    if model.number_solved == 0
        error("VariablePrimal not available.")
    elseif v.result_index > 1
        throw(MOI.ResultIndexBoundsError{MOI.VariablePrimal}(v, 1))
    end
    check_inbounds(model, vi)
    return get_solution(model.inner, vi.value)
end
function MOI.get(model::Optimizer, v::MOI.VariablePrimal, vi::Vector{MOI.VariableIndex})
    if model.number_solved == 0
        error("VariablePrimal not available.")
    elseif v.result_index > 1
        throw(MOI.ResultIndexBoundsError{MOI.VariablePrimal}(v, 1))
    end
    x = get_solution(model.inner)
    return [x[v.value] for v in vi]
end

macro checkcons(model, ci, cp)
    quote
        if $(esc(model)).number_solved == 0
            error("Solve problem before accessing solution.")
        elseif $(esc(cp)).result_index > 1
            throw(MOI.ResultIndexBoundsError{typeof($(esc(cp)))}($(esc(cp)), 1))
        end
        if !(0 <= $(esc(ci)).value <= number_constraints($(esc(model))) - 1)
            error("Invalid constraint index ", $(esc(ci)).value)
        end
    end
end

##################################################
## ConstraintPrimal
function MOI.get(
    model::Optimizer,
    cp::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{S,T},
) where {S<:SF,T<:SS}
    @checkcons(model, ci, cp)
    g = KN_get_con_values(model.inner)
    index = model.constraint_mapping[ci] .+ 1
    return g[index]
end

function MOI.get(
    model::Optimizer,
    cp::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{S,T},
) where {S<:VAF,T<:Union{MOI.Nonnegatives,MOI.Nonpositives}}
    @checkcons(model, ci, cp)
    g = KN_get_con_values(model.inner)
    index = model.constraint_mapping[ci] .+ 1
    return g[index]
end

function MOI.get(
    model::Optimizer,
    cp::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{S,T},
) where {S<:VOV,T<:Union{MOI.Nonnegatives,MOI.Nonpositives}}
    @checkcons(model, ci, cp)
    x = get_solution(model.inner)
    index = model.constraint_mapping[ci] .+ 1
    return x[index]
end

function MOI.get(
    model::Optimizer,
    cp::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{S,T},
) where {S<:Union{VAF,VOV},T<:MOI.Zeros}
    @checkcons(model, ci, cp)
    ncons = length(model.constraint_mapping[ci])
    return zeros(ncons)
end

function MOI.get(
    model::Optimizer,
    cp::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{S,T},
) where {S<:Union{VAF,VOV},T<:MOI.SecondOrderCone}
    @checkcons(model, ci, cp)
    x = get_solution(model.inner)
    index = model.constraint_mapping[ci] .+ 1
    return x[index]
end

function MOI.get(
    model::Optimizer,
    cp::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,<:LS},
)
    if model.number_solved == 0
        error("ConstraintPrimal not available.")
    elseif cp.result_index > 1
        throw(MOI.ResultIndexBoundsError{MOI.ConstraintPrimal}(cp, 1))
    end
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)
    return get_solution(model.inner, vi.value)
end

function MOI.get(
    model::Optimizer,
    cp::MOI.ConstraintPrimal,
    ci::Vector{MOI.ConstraintIndex{MOI.VariableIndex,<:LS}},
)
    if model.number_solved == 0
        error("ConstraintPrimal not available.")
    elseif cp.result_index > 1
        throw(MOI.ResultIndexBoundsError{MOI.ConstraintPrimal}(cp, 1))
    end
    x = get_solution(model.inner)
    return [x[c.value] for c in ci]
end

##################################################
## ConstraintDual
#
# KNITRO's dual sign depends on optimization sense.
sense_dual(model::Optimizer) = (model.sense == MOI.MAX_SENSE) ? 1.0 : -1.0

function MOI.get(
    model::Optimizer,
    cd::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{S,T},
) where {S<:SF,T<:SS}
    @checkcons(model, ci, cd)

    index = model.constraint_mapping[ci] + 1
    lambda = get_dual(model.inner)
    return sense_dual(model) * lambda[index]
end

function MOI.get(
    model::Optimizer,
    cd::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{S,T},
) where {S<:VAF,T<:VLS}
    @checkcons(model, ci, cd)
    index = model.constraint_mapping[ci] .+ 1
    lambda = get_dual(model.inner)
    return sense_dual(model) * lambda[index]
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{S,T},
) where {S<:VOV,T<:VLS}
    offset = number_constraints(model)
    index = model.constraint_mapping[ci] .+ 1 .+ offset
    lambda = get_dual(model.inner)
    return sense_dual(model) * lambda[index]
end

###
# Get constraint of a SOC constraint.
#
# Use the following mathematical property.  Let
#
#   ||u_i || <= t_i      with dual constraint    || z_i || <= w_i
#
# At optimality, we have
#
#   w_i * u_i  = - t_i z_i
#
###
function MOI.get(
    model::Optimizer,
    cd::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{S,T},
) where {S<:Union{VAF,VOV},T<:MOI.SecondOrderCone}
    @checkcons(model, ci, cd)
    index_var = model.constraint_mapping[ci] .+ 1
    index = model.constraint_mapping[ci] .+ 1
    index_con = ci.value + 1
    x = get_solution(model.inner)[index_var]
    # By construction.
    t_i = x[1]
    u_i = x[2:end]
    w_i = get_dual(model.inner)[index_con]

    dual = [-w_i; 1 / t_i * w_i * u_i]

    return dual
end

## Reduced costs.
function MOI.get(
    model::Optimizer,
    cd::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}},
)
    if model.number_solved == 0
        error("ConstraintDual not available.")
    elseif cd.result_index > 1
        throw(MOI.ResultIndexBoundsError{MOI.ConstraintDual}(cd, 1))
    end
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)

    # Constraints' duals are before reduced costs in KNITRO.
    offset = number_constraints(model)
    lambda = sense_dual(model) * get_dual(model.inner, vi.value + offset)
    if lambda < 0
        return lambda
    else
        return 0
    end
end

function MOI.get(
    model::Optimizer,
    cd::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}},
)
    if model.number_solved == 0
        error("ConstraintDual not available.")
    elseif cd.result_index > 1
        throw(MOI.ResultIndexBoundsError{MOI.ConstraintDual}(cd, 1))
    end
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)

    # Constraints' duals are before reduced costs in KNITRO.
    offset = number_constraints(model)
    lambda = sense_dual(model) * get_dual(model.inner, vi.value + offset)
    if lambda > 0
        return lambda
    else
        return 0
    end
end

function MOI.get(
    model::Optimizer,
    cd::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}},
)
    if model.number_solved == 0
        error("ConstraintDual not available.")
    elseif cd.result_index > 1
        throw(MOI.ResultIndexBoundsError{MOI.ConstraintDual}(cd, 1))
    end
    vi = MOI.VariableIndex(ci.value)
    check_inbounds(model, vi)

    # Constraints' duals are before reduced costs in KNITRO.
    offset = number_constraints(model)
    lambda = get_dual(model.inner, vi.value + offset)
    return sense_dual(model) * lambda
end

function MOI.get(model::Optimizer, ::MOI.NLPBlockDual)
    if model.number_solved == 0
        error("NLPBlockDual not available.")
    end
    # Get first index corresponding to a non-linear constraint:
    lambda = get_dual(model.inner)
    # FIXME: Assume that lambda has same sense as for linear
    # and quadratic constraint, but this is not tested inside MOI.
    return sense_dual(model) .* [lambda[i+1] for i in model.nlp_index_cons]
end

###
if KNITRO_VERSION >= v"12.0"
    MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = KN_get_solve_time_cpu(model.inner)
end
# Additional getters
MOI.get(model::Optimizer, ::MOI.NodeCount) = KN_get_mip_number_nodes(model.inner)
MOI.get(model::Optimizer, ::MOI.BarrierIterations) = KN_get_number_iters(model.inner)
MOI.get(model::Optimizer, ::MOI.RelativeGap) = KN_get_mip_rel_gap(model.inner)
MOI.get(model::Optimizer, ::MOI.ObjectiveBound) = KN_get_mip_relaxation_bnd(model.inner)


text:
__precompile__()
module FINUFFT

## Export
export nufft1d1, nufft1d2, nufft1d3
export nufft2d1, nufft2d2, nufft2d3
export nufft3d1, nufft3d2, nufft3d3

export nufft1d1!, nufft1d2!, nufft1d3!
export nufft2d1!, nufft2d2!, nufft2d3!
export nufft3d1!, nufft3d2!, nufft3d3!

export finufft_default_opts
export nufft_opts
export nufft_c_opts # backward-compability

## External dependencies
using Libdl

const depsfile = joinpath(dirname(@__DIR__), "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("FINUFFT is not properly installed. Please build it first.")
end

function __init__()
    Libdl.dlopen(fftw, Libdl.RTLD_GLOBAL)
    if !Sys.iswindows()
        Libdl.dlopen(fftw_threads, Libdl.RTLD_GLOBAL)
    end
end

const BIGINT = Int64 # defined in src/finufft.h

## FINUFFT opts struct from src/finufft.h
"""
    mutable struct nufft_opts
        debug              :: Cint
        spread_debug       :: Cint
        spread_sort        :: Cint
        spread_kerevalmeth :: Cint
        spread_kerpad      :: Cint
        chkbnds            :: Cint
        fftw               :: Cint
        modeord            :: Cint
        upsampfac          :: Cdouble
    end

Options struct passed to the FINUFFT library.

# Fields

    debug :: Cint
0: silent, 1: text basic timing output

    spread_debug :: Cint
passed to spread_opts, 0 (no text) 1 (some) or 2 (lots)

    spread_sort :: Cint
passed to spread_opts, 0 (don't sort) 1 (do) or 2 (heuristic)

    spread_kerevalmeth :: Cint
passed to spread_opts, 0: exp(sqrt()), 1: Horner ppval (faster)

    spread_kerpad :: Cint
passed to spread_opts, 0: don't pad to mult of 4, 1: do

    chkbnds :: Cint
0: don't check if input NU pts in [-3pi,3pi], 1: do

    fftw :: Cint
0:`FFTW_ESTIMATE`, or 1:`FFTW_MEASURE` (slow plan but faster)

    modeord :: Cint
0: CMCL-style increasing mode ordering (neg to pos), or\\
1: FFT-style mode ordering (affects type-1,2 only)

    upsampfac::Cdouble
upsampling ratio sigma, either 2.0 (standard) or 1.25 (small FFT)
"""
mutable struct nufft_opts
    debug              :: Cint
    spread_debug       :: Cint
    spread_sort        :: Cint
    spread_kerevalmeth :: Cint
    spread_kerpad      :: Cint
    chkbnds            :: Cint
    fftw               :: Cint
    modeord            :: Cint
    upsampfac          :: Cdouble
end

const nufft_c_opts = nufft_opts # backward compability

"""
    finufft_default_opts()

Return a [`nufft_opts`](@ref) struct with the default FINUFFT settings.\\
See: <https://finufft.readthedocs.io/en/latest/usage.html#options>
"""
function finufft_default_opts()
    opts = nufft_opts(0,0,0,0,0,0,0,0,0)
    ccall( (:finufft_default_opts, libfinufft),
           Nothing,
           (Ref{nufft_opts},),
           opts
           )
    return opts
end

### Error handling
const ERR_EPS_TOO_SMALL        = 1
const ERR_MAXNALLOC            = 2
const ERR_SPREAD_BOX_SMALL     = 3
const ERR_SPREAD_PTS_OUT_RANGE = 4
const ERR_SPREAD_ALLOC         = 5
const ERR_SPREAD_DIR           = 6
const ERR_UPSAMPFAC_TOO_SMALL  = 7
const HORNER_WRONG_BETA        = 8
const ERR_NDATA_NOTVALID       = 9

struct FINUFFTError <: Exception
    errno::Cint
    msg::String
end
Base.showerror(io::IO, e::FINUFFTError) = print(io, "FINUFFT Error ($(e.errno)): ", e.msg)

function check_ret(ret)
    # Check return value and output error messages
    if ret==0
        return
    elseif ret==ERR_EPS_TOO_SMALL
        msg = "requested tolerance epsilon too small"
    elseif ret==ERR_MAXNALLOC
        msg = "attemped to allocate internal arrays larger than MAX_NF (defined in common.h)"
    elseif ret==ERR_SPREAD_BOX_SMALL
        msg = "spreader: fine grid too small"
    elseif ret==ERR_SPREAD_PTS_OUT_RANGE
        msg = "spreader: if chkbnds=1, a nonuniform point out of input range [-3pi,3pi]^d"
    elseif ret==ERR_SPREAD_ALLOC
        msg = "spreader: array allocation error"
    elseif ret==ERR_SPREAD_DIR
        msg = "spreader: illegal direction (should be 1 or 2)"
    elseif ret==ERR_UPSAMPFAC_TOO_SMALL
        msg = "upsampfac too small (should be >1)"
    elseif ret==HORNER_WRONG_BETA
        msg = "upsampfac not a value with known Horner eval: currently 2.0 or 1.25 only"
    elseif ret==ERR_NDATA_NOTVALID
        msg = "ndata not valid (should be >= 1)"
    else
        msg = "unknown error"
    end
    throw(FINUFFTError(ret, msg))
end

### Simple Interfaces (allocate output)

## Type-1

"""
    nufft1d1(xj      :: Array{Float64},
             cj      :: Array{ComplexF64},
             iflag   :: Integer,
             eps     :: Float64,
             ms      :: Integer
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-1 1D complex nonuniform FFT.
"""
function nufft1d1(xj::Array{Float64},
                  cj::Array{ComplexF64},
                  iflag::Integer,
                  eps::Float64,
                  ms::Integer,
                  opts::nufft_opts=finufft_default_opts())
    fk = Array{ComplexF64}(undef, ms)
    nufft1d1!(xj, cj, iflag, eps, fk, opts)
    return fk
end

"""
    nufft2d1(xj      :: Array{Float64},
             yj      :: Array{Float64},
             cj      :: Array{ComplexF64},
             iflag   :: Integer,
             eps     :: Float64,
             ms      :: Integer,
             mt      :: Integer,
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-1 2D complex nonuniform FFT.
"""
function nufft2d1(xj      :: Array{Float64},
                  yj      :: Array{Float64},
                  cj      :: Array{ComplexF64},
                  iflag   :: Integer,
                  eps     :: Float64,
                  ms      :: Integer,
                  mt      :: Integer,
                  opts    :: nufft_opts = finufft_default_opts())
    fk = Array{ComplexF64}(undef, ms, mt)
    nufft2d1!(xj, yj, cj, iflag, eps, fk, opts)
    return fk
end

"""
    nufft3d1(xj      :: Array{Float64},
             yj      :: Array{Float64},
             zj      :: Array{Float64},
             cj      :: Array{ComplexF64},
             iflag   :: Integer,
             eps     :: Float64,
             ms      :: Integer,
             mt      :: Integer,
             mu      :: Integer,
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-1 3D complex nonuniform FFT.
"""
function nufft3d1(xj      :: Array{Float64},
                  yj      :: Array{Float64},
                  zj      :: Array{Float64},
                  cj      :: Array{ComplexF64},
                  iflag   :: Integer,
                  eps     :: Float64,
                  ms      :: Integer,
                  mt      :: Integer,
                  mu      :: Integer,
                  opts    :: nufft_opts = finufft_default_opts())
    fk = Array{ComplexF64}(undef, ms, mt, mu)
    nufft3d1!(xj, yj, zj, cj, iflag, eps, fk, opts)
    return fk
end

## Type-2

"""
    nufft1d2(xj      :: Array{Float64},
             iflag   :: Integer,
             eps     :: Float64,
             fk      :: Array{ComplexF64}
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-2 1D complex nonuniform FFT.
"""
function nufft1d2(xj      :: Array{Float64},
                  iflag   :: Integer,
                  eps     :: Float64,
                  fk      :: Array{ComplexF64},
                  opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    cj = Array{ComplexF64}(undef, nj)
    nufft1d2!(xj, cj, iflag, eps, fk, opts)
    return cj
end

"""
    nufft2d2(xj      :: Array{Float64},
             yj      :: Array{Float64},
             iflag   :: Integer,
             eps     :: Float64,
             fk      :: Array{ComplexF64}
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-2 2D complex nonuniform FFT.
"""
function nufft2d2(xj      :: Array{Float64},
                  yj      :: Array{Float64},
                  iflag   :: Integer,
                  eps     :: Float64,
                  fk      :: Array{ComplexF64},
                  opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    cj = Array{ComplexF64}(undef, nj)
    nufft2d2!(xj, yj, cj, iflag, eps, fk, opts)
    return cj
end

"""
    nufft3d2(xj      :: Array{Float64},
             yj      :: Array{Float64},
             zj      :: Array{Float64},
             iflag   :: Integer,
             eps     :: Float64,
             fk      :: Array{ComplexF64}
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-2 3D complex nonuniform FFT.
"""
function nufft3d2(xj      :: Array{Float64},
                  yj      :: Array{Float64},
                  zj      :: Array{Float64},
                  iflag   :: Integer,
                  eps     :: Float64,
                  fk      :: Array{ComplexF64},
                  opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    cj = Array{ComplexF64}(undef, nj)
    nufft3d2!(xj, yj, zj, cj, iflag, eps, fk, opts)
    return cj
end

## Type-3

"""
    nufft1d3(xj      :: Array{Float64},
             cj      :: Array{ComplexF64},
             iflag   :: Integer,
             eps     :: Float64,
             sk      :: Array{Float64},
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-3 1D complex nonuniform FFT.
"""
function nufft1d3(xj      :: Array{Float64},
                  cj      :: Array{ComplexF64},
                  iflag   :: Integer,
                  eps     :: Float64,
                  sk      :: Array{Float64},
                  opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(cj)==nj
    nk = length(sk)
    fk = Array{ComplexF64}(undef, nk)
    nufft1d3!(xj, cj, iflag, eps, sk, fk, opts);
    return fk
end

"""
    nufft2d3(xj      :: Array{Float64},
             yj      :: Array{Float64},
             cj      :: Array{ComplexF64},
             iflag   :: Integer,
             eps     :: Float64,
             sk      :: Array{Float64},
             tk      :: Array{Float64}
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-3 2D complex nonuniform FFT.
"""
function nufft2d3(xj      :: Array{Float64},
                  yj      :: Array{Float64},
                  cj      :: Array{ComplexF64},
                  iflag   :: Integer,
                  eps     :: Float64,
                  sk      :: Array{Float64},
                  tk      :: Array{Float64},
                  opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(cj)==nj
    nk = length(sk)
    fk = Array{ComplexF64}(undef, nk)
    nufft2d3!(xj, yj, cj, iflag, eps, sk, tk, fk, opts);
    return fk
end

"""
    nufft3d3(xj      :: Array{Float64},
             yj      :: Array{Float64},
             zj      :: Array{Float64},
             cj      :: Array{ComplexF64},
             iflag   :: Integer,
             eps     :: Float64,
             sk      :: Array{Float64},
             tk      :: Array{Float64}
             uk      :: Array{Float64}
             [, opts :: nufft_opts]
            ) -> Array{ComplexF64}

Compute type-3 3D complex nonuniform FFT.
"""
function nufft3d3(xj      :: Array{Float64},
                  yj      :: Array{Float64},
                  zj      :: Array{Float64},
                  cj      :: Array{ComplexF64},
                  iflag   :: Integer,
                  eps     :: Float64,
                  sk      :: Array{Float64},
                  tk      :: Array{Float64},
                  uk      :: Array{Float64},
                  opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(cj)==nj
    nk = length(sk)
    fk = Array{ComplexF64}(undef, nk)
    nufft3d3!(xj, yj, zj, cj, iflag, eps, sk, tk, uk, fk, opts);
    return fk
end

### Direct interfaces (No allocation)

## 1D

"""
    nufft1d1!(xj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
            )

Compute type-1 1D complex nonuniform FFT. Output stored in fk.
"""
function nufft1d1!(xj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(cj)==nj
    ms = length(fk)
    # Calling interface
    # int finufft1d1(BIGINT nj,FLT* xj,CPX* cj,int iflag,FLT eps,BIGINT ms,
    # 	       CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft1d1, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Cdouble},
                  Ref{ComplexF64},
                  Cint,
                  Cdouble,
                  BIGINT,
                  Ref{ComplexF64},
                  nufft_opts),
                 nj, xj, cj, iflag, eps, ms, fk, opts
                 )
    check_ret(ret)
end

"""
    nufft1d2!(xj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
            )

Compute type-2 1D complex nonuniform FFT. Output stored in cj.
"""
function nufft1d2!(xj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(cj)==nj
    ms = length(fk)
    # Calling interface
    # int finufft1d2(BIGINT nj,FLT* xj,CPX* cj,int iflag,FLT eps,BIGINT ms,
    #                CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft1d2, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Cdouble},
                  Ref{ComplexF64},
                  Cint,
                  Cdouble,
                  BIGINT,
                  Ref{ComplexF64},
                  nufft_opts),
                 nj, xj, cj, iflag, eps, ms, fk, opts
                 )
    check_ret(ret)
end

"""
    nufft1d3!(xj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              sk      :: Array{Float64},
              fk      :: Array{ComplexF64},
              [, opts :: nufft_opts]
             )

Compute type-3 1D complex nonuniform FFT. Output stored in fk.
"""
function nufft1d3!(xj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   sk      :: Array{Float64},
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(cj)==nj
    nk = length(sk)
    @assert length(fk)==nk
    # Calling interface
    # int finufft1d3(BIGINT nj,FLT* x,CPX* c,int iflag,FLT eps,BIGINT nk, FLT* s, CPX* f, nufft_opts opts);
    ret = ccall( (:finufft1d3, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Cdouble},
                  Ref{ComplexF64},
                  Cint,
                  Cdouble,
                  BIGINT,
                  Ref{Cdouble},
                  Ref{ComplexF64},
                  nufft_opts),
                 nj, xj, cj, iflag, eps, nk, sk, fk, opts
                 )
    check_ret(ret)
end

## 2D

"""
    nufft2d1!(xj      :: Array{Float64},
              yj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
            )

Compute type-1 2D complex nonuniform FFT. Output stored in fk.
"""
function nufft2d1!(xj      :: Array{Float64},
                   yj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(cj)==nj
    ms, mt = size(fk)
    # Calling interface
    # int finufft2d1(BIGINT nj,FLT* xj,FLT *yj,CPX* cj,int iflag,FLT eps,
    #                BIGINT ms, BIGINT mt, CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft2d1, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Cdouble},
                  Ref{Cdouble},
                  Ref{ComplexF64},
                  Cint,
                  Cdouble,
                  BIGINT,
                  BIGINT,
                  Ref{ComplexF64},
                  nufft_opts),
                 nj, xj, yj, cj, iflag, eps, ms, mt, fk, opts
                 )
    check_ret(ret)
end

"""
    nufft2d2!(xj      :: Array{Float64},
              yj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
            )

Compute type-2 2D complex nonuniform FFT. Output stored in cj.
"""
function nufft2d2!(xj      :: Array{Float64},
                   yj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(cj)==nj
    ms, mt = size(fk)
    # Calling interface
    # int finufft2d2(BIGINT nj,FLT* xj,FLT *yj,CPX* cj,int iflag,FLT eps,
    #                BIGINT ms, BIGINT mt, CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft2d2, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Cdouble},
                  Ref{Cdouble},
                  Ref{ComplexF64},
                  Cint,
                  Cdouble,
                  BIGINT,
                  BIGINT,
                  Ref{ComplexF64},
                  nufft_opts),
                 nj, xj, yj, cj, iflag, eps, ms, mt, fk, opts
                 )
    check_ret(ret)
end

"""
    nufft2d3!(xj      :: Array{Float64},
              yj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              sk      :: Array{Float64},
              tk      :: Array{Float64},
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
             )

Compute type-3 2D complex nonuniform FFT. Output stored in fk.
"""
function nufft2d3!(xj      :: Array{Float64},
                   yj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   sk      :: Array{Float64},
                   tk      :: Array{Float64},
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(cj)==nj
    nk = length(sk)
    @assert length(tk)==nk
    @assert length(fk)==nk
    # Calling interface
    # int finufft2d3(BIGINT nj,FLT* x,FLT *y,CPX* cj,int iflag,FLT eps,BIGINT nk, FLT* s, FLT* t, CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft2d3, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Cdouble},
                  Ref{Cdouble},
                  Ref{ComplexF64},
                  Cint,
                  Cdouble,
                  BIGINT,
                  Ref{Cdouble},
                  Ref{Cdouble},
                  Ref{ComplexF64},
                  nufft_opts),
                 nj, xj, yj, cj, iflag, eps, nk, sk, tk, fk, opts
                 )
    check_ret(ret)
end

## 3D

"""
    nufft3d1!(xj      :: Array{Float64},
              yj      :: Array{Float64},
              zj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
            )

Compute type-1 3D complex nonuniform FFT. Output stored in fk.
"""
function nufft3d1!(xj      :: Array{Float64},
                   yj      :: Array{Float64},
                   zj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(zj)==nj
    @assert length(cj)==nj
    ms, mt, mu = size(fk)
    # Calling interface
    # int finufft3d1(BIGINT nj,FLT* xj,FLT *yj,FLT *zj,CPX* cj,int iflag,FLT eps,
    # 	       BIGINT ms, BIGINT mt, BIGINT mu, CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft3d1, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Cdouble},
                  Ref{Cdouble},
                  Ref{Cdouble},
                  Ref{ComplexF64},
                  Cint,
                  Cdouble,
                  BIGINT,
                  BIGINT,
                  BIGINT,
                  Ref{ComplexF64},
                  nufft_opts),
                 nj, xj, yj, zj, cj, iflag, eps, ms, mt, mu, fk, opts
                 )
    check_ret(ret)
end

"""
    nufft3d2!(xj      :: Array{Float64},
              yj      :: Array{Float64},
              zj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
            )

Compute type-2 3D complex nonuniform FFT. Output stored in cj.
"""
function nufft3d2!(xj      :: Array{Float64},
                   yj      :: Array{Float64},
                   zj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(zj)==nj
    @assert length(cj)==nj
    ms, mt, mu = size(fk)
    # Calling interface
    # int finufft3d2(BIGINT nj,FLT* xj,FLT *yj,FLT *zj,CPX* cj,int iflag,FLT eps,
    #                BIGINT ms, BIGINT mt, BIGINT mu, CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft3d2, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Cdouble},
                  Ref{Cdouble},
                  Ref{Cdouble},
                  Ref{ComplexF64},
                  Cint,
                  Cdouble,
                  BIGINT,
                  BIGINT,
                  BIGINT,
                  Ref{ComplexF64},
                  nufft_opts),
                 nj, xj, yj, zj, cj, iflag, eps, ms, mt, mu, fk, opts
                 )
    check_ret(ret)
end

"""
    nufft3d3!(xj      :: Array{Float64},
              yj      :: Array{Float64},
              zj      :: Array{Float64},
              cj      :: Array{ComplexF64},
              iflag   :: Integer,
              eps     :: Float64,
              sk      :: Array{Float64},
              tk      :: Array{Float64},
              uk      :: Array{Float64},
              fk      :: Array{ComplexF64}
              [, opts :: nufft_opts]
             )

Compute type-3 3D complex nonuniform FFT. Output stored in fk.
"""
function nufft3d3!(xj      :: Array{Float64},
                   yj      :: Array{Float64},
                   zj      :: Array{Float64},
                   cj      :: Array{ComplexF64},
                   iflag   :: Integer,
                   eps     :: Float64,
                   sk      :: Array{Float64},
                   tk      :: Array{Float64},
                   uk      :: Array{Float64},
                   fk      :: Array{ComplexF64},
                   opts    :: nufft_opts = finufft_default_opts())
    nj = length(xj)
    @assert length(yj)==nj
    @assert length(zj)==nj
    @assert length(cj)==nj
    nk = length(sk)
    @assert length(tk)==nk
    @assert length(uk)==nk
    @assert length(fk)==nk
    # Calling interface
    # int finufft3d3(BIGINT nj,FLT* x,FLT *y,FLT *z, CPX* cj,int iflag,
    #                FLT eps,BIGINT nk,FLT* s, FLT* t, FLT *u,
    #                CPX* fk, nufft_opts opts);
    ret = ccall( (:finufft3d3, libfinufft),
                 Cint,
                 (BIGINT,
                  Ref{Cdouble},
                  Ref{Cdouble},
                  Ref{Cdouble},
                  Ref{ComplexF64},
                  Cint,
                  Cdouble,
                  BIGINT,
                  Ref{Cdouble},
                  Ref{Cdouble},
                  Ref{Cdouble},
                  Ref{ComplexF64},
                  nufft_opts),
                 nj, xj, yj, zj, cj, iflag, eps, nk, sk, tk, uk, fk, opts
                 )
    check_ret(ret)
end

end # module


text:
#=
  From https://github.polettix.it/ETOOBUSY/2021/07/19/brute-force-puzzle/
  """
  Consider a display with N digits in base b

      A random value id shown on the display; each possible value is equiprobable.
      What is the expected number of missing digit values?

  As an example, suppose that we consider base-4 digits (i.e. 0 to 3) and a display that
  has 6 slots. An example value shown on the display would be 232133, which includes digits 1, 2,
  and 3 but not digit 0. In this arrangement, digit 0 is considered missing. Similarly, in
  arrangement 111111 there are three missing digits, namely 0, 2, and 3.

  ....

  [By hand and Brute force:]
  2916 / 4096 ≅ 0.711914
  """

  This Turing.jl model:
  """
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
       Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64
  ....
  numMissing    0.7188    0.6277     0.0063    0.0079   6612.2186    0.9999      943.1206

  ...
  Distributions of variable numMissing (num:0)
  1.00000 =>    5291  (0.529100)
  0.00000 =>    3764  (0.376400)
  2.00000 =>     938  (0.093800)
  3.00000 =>       7  (0.000700)
  """

  Cf the WebPPL model probability_of_missing_values.wppl which give the exact probability
  0.7119140625000009

=#

using Turing
include("jl_utils.jl")

@model function probability_of_missing_values(base=4,n=6)
    # digit = tzeros(n)
    # for i in 1:n
    #     digit[i] ~ DiscreteUniform(0,base-1)
    # end
    # Faster:
    digit ~ filldist(DiscreteUniform(0,base-1), n)
    missingValues = tzeros(base)
    for i in 0:base-1
        missingValues[i+1] ~ Dirac(sum([digit[j] == i for j in 1:n]) == 0)
    end
    numMissing ~ Dirac(sum(missingValues))
end

base = 4
n = 6
model = probability_of_missing_values(base,n)

num_chns = 4
# chns = sample(model, Prior(), 1000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, IS(), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, SMC(), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, SGLD(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns,:numMissing)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
export FieldType, BC  # types
export EH, nE, nH, nEH  # instances
export ft2gt  # functions

# Field types
const nE, nH = 1, 2  # E-, H-fields
const nEH = SVec(nE, nH)
@enum FieldType EE=nE HH
const EH = SVec(EE, HH)
for ins in instances(FieldType); @eval export $(Symbol(ins)); end  # export all instances
Base.string(ins::FieldType) = ins==EE ? "E" : "H"
StaggeredGridCalculus.alter(ins::FieldType) = ins==EE ? HH : EE

# Given boundary field types, determine whether the grid planes specified by the given field
# type ft are primal or dual grid planes in the Cartesian directions.
# Note that the boundary field type is the type of the fields tangential to the boundary.
# The fields defined on but normal to the boundary have the complementary type.
ft2gt(ft::FieldType, boundft::FieldType) = PD[2 - (boundft==ft)]

# Boundary conditions
@enum BC PERIODIC=1 CONDUCTING  # periodic boundary condition, perfect conductor boundary condition
for ins in instances(BC); @eval export $(Symbol(ins)); end  # export all instances
Base.string(ins::BC) = ins==PERIODIC ? "periodic" : "perfect conductor"


text:
using PyPlot
using JLD

path = "/Users/milan/cams"
allvars = ["no2","go3","so2","aermr04","aermr05","aermr06","ch4","co","co2","no"]
vari = "so2"

# volume and error norms (columns):
# vol, mean, median, 90% and max  of decimal error
# mean, median, 90% and max of abs error
# total 9

# compression methods (rows):
# LinQuant24, LogQuant16, RoundNearest16/24 + Blosc/Deflate/Zstd/LZ4
# total 10

D = load("cams/error/$vari.jld")
EV = D["EV"]
lz4hcc = D["lz4"]

# theoretical
theosize = [32/i for i in 32:-1:5]
x0 = 1f0
theoerror = [log10(nextfloat(x0,2^i)) for i in 0:27]

## plotting
pygui(true)
n = size(EV)[1]

cm1 = ColorMap("viridis")
cm2 = ColorMap("inferno")

function RGB(cv::Tuple)
    r = UInt8(round(cv[1]*255))
    g = UInt8(round(cv[2]*255))
    b = UInt8(round(cv[3]*255))
    return "#"*repr(r)[3:4]*repr(g)[3:4]*repr(b)[3:4]
end

colours = ["0.5","0.3",RGB(cm2(0.9)),RGB(cm2(0.8)),RGB(cm2(0.6)),RGB(cm2(0.4)),
                    RGB(cm1(0.8)),RGB(cm1(0.6)),RGB(cm1(0.4)),RGB(cm1(0.2))]

names = ["LinQuant24","LogQuant16","RoundNearest16+Blosc","RoundNearest16+LZ4HC",
            "RoundNearest16+Zstd","RoundNearest16+Deflate","RoundNearest24+Blosc",
            "RoundNearest24+LZ4HC","RoundNearest24+Zstd","RoundNearest24+Deflate"]

markers = ["s","s","d","d","d","d","o","o","o","o"]

fig,ax = subplots(1,1,figsize=(8,6))
ax.set_xscale("log")

ax.plot(theoerror[1],1,"w",markersize=15,marker="*",lw=0,markeredgewidth=1.4,
            markeredgecolor="k",alpha=1,label="Float32",zorder=10)
ax.plot(theoerror[1],lz4hcc,"w",markersize=10,marker="v",lw=0,markeredgewidth=1.4,
            markeredgecolor="k",alpha=0.7,label="Float32+LZ4HC",zorder=10)

for i in [1,2,3,4,5,6,7,8,9,10]
    # ax.plot(EV[i,[3,5]],[EV[i,1],EV[i,1]],"k",lw=5.5,zorder=-1)
    ax.plot(EV[i,3:4],[EV[i,1],EV[i,1]],colours[i],
            marker=markers[i],markevery=2,
            markeredgecolor="k",markersize=10,
            lw=5,label=names[i])
    ax.plot(EV[i,4:5],[EV[i,1],EV[i,1]],colours[i],lw=2)
end

ax.plot(theoerror,theosize,"k",label="Theoretical rounding",zorder=-1)
ax.set_xlim(theoerror[1]/3,100)
ax.set_ylim(0.3,1+ceil(maximum(EV[:,1])))
# ax.set_ylim(0,110)
ax.set_xlabel("decimal error")
ax.set_ylabel("compression factor")
ax.legend(loc=1,ncol=1,fontsize=9)
ax.set_title("NO compression",loc="left")

ax.text(EV[1,3],EV[1,1],"median  ",rotation=90,ha="center",va="top",fontsize=9)
ax.text(EV[1,4],EV[1,1],"90%  ",rotation=90,ha="center",va="top",fontsize=9)
ax.text(EV[1,5],EV[1,1],"max  ",rotation=90,ha="center",va="top",fontsize=9)

tight_layout()


text:
using Convex
using Convex.ProblemDepot: run_tests
using Test
using SCS, ECOS, GLPK

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

# Seed random number stream to improve test reliability
using Random
Random.seed!(2)

@testset "ProblemDepot" begin
    @testset "Problems can run without `solve!`ing if `test==false`; T=$T" for T in (Float64, BigFloat)
        Convex.ProblemDepot.foreach_problem() do name, func
            @testset "$name" begin
                # We want to check to make sure this does not throw
                func(Val(false), 0.0, 0.0, T) do problem
                    @test problem isa Convex.Problem{T} # check numeric type
                    model = MOIU.MockOptimizer(MOIU.Model{T}())
                    Convex.load_MOI_model!(model, problem) # make sure it loads without throwing
                end
            end
        end
    end
end

@testset "Convex" begin
    include("test_utilities.jl")

    @testset "SCS with warmstarts" begin
        # We exclude `sdp_matrix_frac_atom` due to the bug https://github.com/JuliaOpt/SCS.jl/issues/153
        # Likewise, `sdp_sdp_constraints` hits a bug https://github.com/JuliaOpt/SCS.jl/issues/167
        run_tests(; exclude=[r"mip", r"sdp_matrix_frac_atom", r"sdp_sdp_constraints"]) do p
            solve!(p, SCS.Optimizer(verbose=0, eps=1e-6); warmstart = true)
        end
    end

    @testset "SCS" begin
        # Exclusions same as for "SCS with warmstarts"
        run_tests(; exclude=[r"mip", r"sdp_matrix_frac_atom", r"sdp_sdp_constraints"]) do p
            solve!(p, SCS.Optimizer(verbose=0, eps=1e-6))
        end
    end

    @testset "ECOS" begin
        run_tests(; exclude=[r"mip", r"sdp"]) do p
            solve!(p, ECOS.Optimizer(verbose=0))
        end
    end

    @testset "GLPK" begin
        run_tests(; exclude=[r"exp", r"sdp", r"socp"]) do p
            solve!(p, GLPK.Optimizer(msg_lev = GLPK.OFF))
        end
    end
end


text:
#=
run MOI tests
=#

using Test
import MathOptInterface
const MOI = MathOptInterface
import Hypatia
include(joinpath(@__DIR__, "moicones.jl"))

@testset "MathOptInterface wrapper tests" begin

@testset "MOI cone tests" begin
    println("starting MOI wrapper cone tests")
    real_types = [
        Float64,
        # Float32,
        BigFloat,
        ]
    for T in real_types
        @testset "$T" begin
            println(T, " ...")
            test_moi_cones(T)
        end
    end
end

# real types, tolerances, and tests to include for MOI.Test tests
test_T = [
    (Float64, 2 * sqrt(sqrt(eps())), 4, String[]),
    # TODO add test_linear after MOI 0.10.7 is tagged:
    (BigFloat, 2 * eps(BigFloat)^0.2, 1, String["test_conic", "test_modification"]),
]

@testset "MOI.Test tests: $T" for (T, tol_test, tol_relax, includes) in test_T
    println("\nstarting MOI.Test tests: $T")
    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{T}()),
            Hypatia.Optimizer{T}(; default_tol_relax = tol_relax),
        ),
        T,
    )
    MOI.set(model, MOI.Silent(), true)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            T,
            atol = tol_test,
            rtol = tol_test,
            exclude = Any[
                MOI.ConstraintBasisStatus,
                MOI.VariableBasisStatus,
                MOI.ObjectiveBound,
                MOI.SolverVersion,
            ],
        ),
        include = includes,
        exclude = String[
            # TODO(odow): unexpected failure, probably in the bridge layer
            "test_model_UpperBoundAlreadySet",
            "test_model_LowerBoundAlreadySet",
        ],
    )
end

end
;