text:
function transformation_matrices(Nx::Int32, Ny::Int32)
    # calculates the interpolation and restriction matrices given a certain resolution in each direction
    # Nx = grid.Nx
    # Ny = grid.Ny
    coarsexres = Int32(Nx/2)
    coarseyres = Int32(Ny/2)

    xcoarse = zeros(Float64, coarsexres)
    ycoarse = zeros(Float64, coarseyres)  # using the reduction in number of points from the paper
    # way to define these with each new iteration? Another overarching loop with the number of iterations?
    # could do a while loop potentially, if trying to get under a threshold

    # define the interpolation matrix, then the restriction matrix as its transpose

    # only working in 1D so far
    Interpx = zeros(Nx, coarsexres)
    # Restrx1D = zeros(coarsexres, Nx)

    for i=1:Nx
        for j=1:coarsexres
            # sort of brute force based on the form of the interpolation matrix given in the MIT site
            Interpx[2(j-1)+1, j] = .5
            Interpx[2(j-1)+2, j] = 1
            if 2(j-1)+3 <= i
                Interpx[2(j-1)+3, j] = .5
            end
            # may not need outer loop if indexing only depends on j
        end
    end

    Restrx1D = .5*transpose(Interpx) # definition of restriction matrix from the paper

    # moving to 2D?
    # defining values based on 2D interpolation scheme in the paper
    """for i=1:Nx
        for j=1:coarsexres
            if isodd(i)
                if isodd(j)
                    Interpx[i,j] = .25
                elseif iseven(j)
                    Interpx[i,j] = .5
                end
            elseif iseven(i)
                if isodd(j)
                    Interpx[i,j] = .5
                elseif iseven(j)
                    Interpx[i,j] = 1
                end
            end
        end
    end"""

    # Restrx = .25*transpose(Interpx) # definition of restriction matrix
    Restr2d = kron(Restrx1D, Restrx1D) # assuming we're using the 1D matrix
    Interp2d = 4*transpose(Restr2d)
    # same Interpolation/Restriction matrices in x and y?

    return Restr2d, Interp2d, coarsexres, coarseyres
end

# not sure if fine_to_coarse and coarse_to_fine are needed, but they could be useful so I'll keep them in
function fine_to_coarse(Nx::Int32, Ny::Int32, x::Vector{Float64}, y::Vector{Float64})
    # restricts values from fine grid to fit them onto coarse grid

    Restriction, Interpolation, newxres, newyres = transformation_matrices(Nx, Ny)

    # calculating the coarse resolution for reducing the number of grid points by 2 in each dimension, assuming odd number of points in each dimension
    coarsexres = Int32(Nx/2)
    coarseyres = Int32(Ny/2)

    # initializing coarse grid values in each dimension
    xcoarse = zeros(Float64, coarsexres)
    ycoarse = zeros(Float64, coarseyres)

    # assuming dimensionality of both are the same?
    xcoarse[:] = Restriction * x
    ycoarse[:] = Restriction * y

    return xcoarse, ycoarse
end

function fine_to_coarse(Nx::Int32, Ny::Int32, x::Vector{Float64}, y::Vector{Float64})
    # projects vector values from the fine grid to the coarse grid using interpolation

    Restriction, Interpolation, newxres, newyres = transformation_matrices(Nx, Ny)

    # calculating the coarse resolution for reducing the number of grid points by 2 in each dimension, assuming odd number of points in each dimension
    coarsexres = ceil((Nx-1)/2)
    coarseyres = ceil((Ny-1)/2)

    # initializing coarse grid values in each dimension
    xcoarse = zeros(Float64, coarsexres)
    ycoarse = zeros(Float64, coarseyres)

    # assuming dimensionality of both are the same?
    xcoarse[:] = Interpolation * x
    ycoarse[:] = Interpolation * y

    return xcoarse, ycoarse
end

function project_to_coarse(A::SparseMatrixCSC, Nx::Int32, Ny::Int32, R::Matrix{Float64}, I::Matrix{Float64})
    # projects the matrix from the matrix equation Ax = b onto the coarse grid

    # acquiring restriction and interpolation matrices
    # R, I, xres, yres = transformation_matrices(Nx, Ny)

    # using formula for projection of matrix onto coarse grid
    Acoarse = R * A * I

    return Acoarse
end

function correction_matrix(A::SparseMatrixCSC, Nx::Int32, Ny::Int32)
    # calculates the error correction matrix s

    # acquiring restriction and interpolation matrices
    R, I, xres, yres = transformation_matrices(Nx, Ny)

    # carrying out explicit calculation of s
    S = I * inv(R * A * I) * R * A

    return S
end

# NOTE THAT THE x I AM USING HERE IS THE ACTUAL SOLUTION AND NOT THE GUESSED/KNOWN SOLUTION
# Also the x above is x as in the coordinate, below the x is the solution to the matrix equation Ax = b

function residual(A::SparseMatrixCSC, x::Vector{Float64}, b::Vector{Float64})
    # calculates the residual associated with the vector equation we're trying to solve

    # definition of residual on fine grid
    rfine = b - A * x

    return rfine
end

function residual(A::Matrix{Float64}, x::Vector{Float64}, b::Vector{Float64})
    # calculates the residual associated with the vector equation we're trying to solve

    # definition of residual on fine grid
    rfine = b - A * x

    return rfine
end

function restricted_residual(A::SparseMatrixCSC, x::Vector{Float64}, b::Vector{Float64}, Nx::Int32, Ny::Int32)
    # calculates the residual on the coarse grid

    # getting restriction/interpolation matrices on the given grid
    R, I, xres, yres = transformation_matrices(Nx::Int32, Ny::Int32)

    # obtaining residual on fine grid
    rfine = residual(A, x, b)

    # calculation of restriction residual
    rcoarse = R * rfine

    return rcoarse
end

function multigrid_error(e::Vector{Float64}, A::SparseMatrixCSC, Nx::Int32, Ny::Int32)
    # given the error e = x(guessed) - x(calculated), this gives the multigrid error correction

    # calculation of correction matrix S that maps from coarse to fine grid
    S = correction_matrix(A, Nx, Ny)

    # definition of multigrid correction given true fine grid error
    E = S * e

    return E
end

# Need to solve Acoarse = Ecoarse * rcoarse, then interpolate back and add Efine to x(calculated)

# may not need y vector - check later if issue doesn't get resolved
function multigrid_step(grid::Grid, A::SparseMatrixCSC, xcalc::Vector{Float64}, xguess::Vector{Float64}, b::Vector{Float64})
    # steps through one step of multigrid from fine to coarse and back according to the steps in the MIT page

    Nx = grid.Nx
    Ny = grid.Ny

    R1, I1, Nx1, Ny1 = transformation_matrices(Nx, Ny)
    R2, I2, Nx2, Ny2 = transformation_matrices(Nx1, Ny1)
    R3, I3, Nx3, Ny3 = transformation_matrices(Nx2, Ny2)

    rcoarse1 = restricted_residual(A, xcalc, b, Nx, Ny)
    Acoarse1 = project_to_coarse(A, Nx, Ny, R1, I1)
    xcalc1 = R1 * xcalc
    b1 = R1 * b

    rcoarse2 = restricted_residual(Acoarse1, xcalc1, b1, Nx1, Ny1)
    Acoarse2 = project_to_coarse(Acoarse1, Nx1, Ny1, R2, I2)
    xcalc2 = R2 * xcalc
    b2 = R2 * b

    rcoarse3 = restricted_residual(Acoarse2, xcalc2, b2, Nx2, Ny2)
    Acoarse3 = project_to_coarse(Acoarse2, Nx2, Ny2, R3, I3)
    xcalc3 = R3 * xcalc
    b3 = R3 * b

    # actual calculation - use a vcycle for example?
    # iterate at each step, carry over E into next grid and use it as guess for iteration on that grid
    # potentially start with exact solution at lowest grid level
    # for example
    Ecoarse3 = Acoarse3 \ rcoarse3
    Ecoarse2 = I3 * Ecoarse3
    Ecoarse2_1 = p_jacobi(Acoarse2, Ecoarse2, rcoarse2, 1, 1, 1, .1)
    Ecoarse1 = I2 * Ecoarse2_1
    Ecoarse1_1 = p_jacobi(Acoarse1, Ecoarse1, rcoarse1, 1, 1, 1, .1)
    Efine = I1 * Ecoarse1_1

    xcalc = xcalc + Efine

    return xcalc # may want to return more depending on the information we want
end

function v_cycle(grid::Grid, A::SparseMatrixCSC, xcalc::Vector{Float64}, b::Vector{Float64})
    # steps through one step of multigrid from fine to coarse and back according to the steps in the MIT page

    Nx = grid.Nx
    Ny = grid.Ny

    R1, I1, Nx1, Ny1 = transformation_matrices(Nx, Ny)
    R2, I2, Nx2, Ny2 = transformation_matrices(Nx1, Ny1)
    R3, I3, Nx3, Ny3 = transformation_matrices(Nx2, Ny2)

    rcoarse1 = restricted_residual(A, xcalc, b, Nx, Ny)
    Acoarse1 = project_to_coarse(A, Nx, Ny, R1, I1)
    xcalc1 = R1 * xcalc
    b1 = R1 * b

    rcoarse2 = restricted_residual(Acoarse1, xcalc1, b1, Nx1, Ny1)
    Acoarse2 = project_to_coarse(Acoarse1, Nx1, Ny1, R2, I2)
    xcalc2 = R2 * xcalc1
    b2 = R2 * b1

    rcoarse3 = restricted_residual(Acoarse2, xcalc2, b2, Nx2, Ny2)
    Acoarse3 = project_to_coarse(Acoarse2, Nx2, Ny2, R3, I3)
    xcalc3 = R3 * xcalc2
    b3 = R3 * b2

    Einitial = zeros(Float64, size(xcalc))
    Ecoarse1 = R1 * Einitial
    Ecoarse1_1 = jacobi_smooth(Acoarse1, Ecoarse1, rcoarse1, 1, 1, 1, 1, .1)
    Ecoarse2 = R2 * Ecoarse1_1
    Ecoarse2_1 = jacobi_smooth(Acoarse2, Ecoarse2, rcoarse2, 1, 1, 1, 1, .1)
    Ecoarse3 = R3 * Ecoarse2_1
    Ecoarse3_1 = jacobi_smooth(Acoarse3, Ecoarse3, rcoarse3, 1, 1, 1, 1, .1)
    Ecoarse2_2 = I3 * Ecoarse3_1
    Ecoarse2_3 = jacobi_smooth(Acoarse2, Ecoarse2_2, rcoarse2, 1, 1, 1, 1, .1)
    Ecoarse1_2 = I2 * Ecoarse2_3
    Ecoarse1_3 = jacobi_smooth(Acoarse1, Ecoarse1_2, rcoarse1, 1, 1, 1, 1, .1)
    Efine = I1 * Ecoarse1_3

    xcalc = xcalc + Efine

    return xcalc # may want to return more depending on the information we want
end

function v_cycle_general(grids::Diffresgrids, A::SparseMatrixCSC, x::Vector{Float64}, b::Vector{Float64})
    # doing a v-cycle in the formatting of the rest of the code (no restricted grids)

    # loading the restriction/interpolation matrices into concise variables
    R1 = grids.R1
    R2 = grids.R2
    R3 = grids.R3
    I1 = grids.I1
    I2 = grids.I2
    I3 = grids.I3
    Nx = grids.Nx
    Ny = grids.Ny
    Nx1 = grids.Nx1
    Ny1 = grids.Ny1
    Nx2 = grids.Nx2
    Ny2 = grids.Ny2
    Nx3 = grids.Nx3
    Ny3 = grids.Ny3

    # calculating the form of the matrix/residual at different resolutions
    """- error is calculated through smoothing (using given error and residual)
    - residual is calculated using error and previous residual as well as spacing
    - residual is restricted
    - error is set to 0
    IF NOT AT LOWEST RESOLUTION
    - run cycle again using the new error/residual to get error
    IF AT LOWEST RESOLUTION
    - error is given by smoothing function at lowest resolution"""

    Acoarse1 = sparse(project_to_coarse(A, Nx, Ny, R1, I1))
    x1 = R1 * x
    b1 = R1 * b
    r = residual(Acoarse1, x1, b1)
    Einitial1 = zeros(Float64, size(x1))
    Ecoarse1 = jacobi_smooth(Acoarse1, Einitial1, r, 1., 1, 1, 1)
    rcoarse1 = Acoarse1 * Ecoarse1
    rcoarse2 = R2 * rcoarse1

    Acoarse2 = sparse(project_to_coarse(Acoarse1, Nx1, Ny1, R2, I2))
    x2 = R2 * x1
    b2 = R2 * b1
    Einitial2 = zeros(Float64, size(x2))
    Ecoarse2 = jacobi_smooth(Acoarse2, Einitial2, rcoarse2, 1., 1, 1, 1)
    rcoarse2_1 = Acoarse2 * Ecoarse2
    rcoarse3 = R3 * rcoarse2_1

    Acoarse3 = sparse(project_to_coarse(Acoarse2, Nx2, Ny2, R3, I3))
    x3 = R3 * x2
    b3 = R3 * b2
    Einitial3 = zeros(Float64, size(x3))
    Ecoarse3 = jacobi_smooth(Acoarse3, Einitial3, rcoarse3, 1., 1, 1, 1) # maybe want to do exact solution at lowest res?
    rcoarse3_1 = Acoarse3 * Ecoarse3
    rcoarse2_2 = I3 * rcoarse3_1

    Ecoarse2_1 = jacobi_smooth(Acoarse2, Einitial2, rcoarse2_2, 1., 1, 1, 1)
    rcoarse2_3 = Acoarse2 * Ecoarse2_1
    rcoarse1_2 = I2 * rcoarse2_3

    Ecoarse1_1 = jacobi_smooth(Acoarse1, Einitial1, rcoarse1_2, 1., 1, 1, 1)
    rcoarse1_3 = Acoarse1 * Ecoarse1_1 # these two calculations may not be necessary
    rfine = I1 * rcoarse1_3
    Efine = Ecoarse1_1

    # actual v-cycle steps
    """Einitial = zeros({Float64}, size(x))
    Ecoarse1 = R1 * Einitial
    Ecoarse1_1 = jacobi_smooth(Acoarse1, Ecoarse1, rcoarse1, 1, 1, 1, 1, .1)
    Ecoarse2 = R2 * Ecoarse1_1
    Ecoarse2_1 = jacobi_smooth(Acoarse2, Ecoarse2, rcoarse2, 1, 1, 1, 1, .1)
    Ecoarse3 = R3 * Ecoarse2_1
    Ecoarse3_1 = jacobi_smooth(Acoarse3, Ecoarse3, rcoarse3, 1, 1, 1, 1, .1)
    Ecoarse2_2 = I3 * Ecoarse3_1
    Ecoarse2_3 = jacobi_smooth(Acoarse2, Ecoarse2_2, rcoarse2, 1, 1, 1, 1, .1)
    Ecoarse1_2 = I2 * Ecoarse2_3
    Ecoarse1_3 = jacobi_smooth(Acoarse1, Ecoarse1_2, rcoarse1, 1, 1, 1, 1, .1)
    Efine = I1 * Ecoarse1_3"""

    x = x + Efine

    return x, Efine
end


text:
using Test
using JuMP
import DiffOpt
import MathOptInterface
const MOI = MathOptInterface
import LinearAlgebra: dot, ⋅, Diagonal
import SparseArrays: sparse

macro _test(computed, expected::Symbol)
    exp = esc(expected)
    com = esc(computed)
    ato = esc(:atol)
    rto = esc(:rtol)
    return :(
        if $exp === nothing
            @test ($exp = $com) isa Any
        else
            @test $com ≈ $exp atol = $ato rtol = $rto
        end
    )
end

"""
    qp_test(solver; kws...)

Check our results for QPs using the notations of [AK17].

[AK17] Amos, Brandon, and J. Zico Kolter. "Optnet: Differentiable optimization as a layer in neural networks." International Conference on Machine Learning. PMLR, 2017. https://arxiv.org/pdf/1703.00443.pdf
"""
function qp_test(
    solver,
    lt::Bool,
    set_zero::Bool,
    canonicalize::Bool;
    dzb = nothing,
    n = length(dzb),
    q = nothing,
    dqf = zeros(n),
    dqb = nothing,
    Q = zeros(n, n),
    dQf = zeros(n, n),
    dQb = nothing,
    fix_indices = Int[],
    fix_values = Float64[],
    nfix = length(fix_values),
    ub_indices = Int[],
    ub_values = Float64[],
    nub = length(ub_values),
    lb_indices = Int[],
    lb_values = Float64[],
    nlb = length(lb_values),
    h = zeros(0),
    nle = length(h),
    dhf = zeros(nle + nub + nlb),
    dhb = nothing,
    G = zeros(0, n),
    dGf = zeros(nle + nub + nlb, n),
    dGb = nothing,
    b = zeros(0),
    neq = length(b),
    dbf = zeros(neq + nfix),
    dbb = nothing,
    A = zeros(0, n),
    dAf = zeros(neq + nfix, n),
    dAb = nothing,
    z = nothing,
    dzf = nothing,
    ∇zf = nothing,
    ∇zb = nothing,
    λ = nothing,
    dλf = zeros(nle + nub + nlb),
    dλb = zeros(nle + nub + nlb),
    ∇λf = nothing,
    ∇λb = nothing,
    ν = nothing,
    dνf = nothing,
    dνb = zeros(neq + nfix),
    ∇νf = nothing,
    ∇νb = nothing,
    atol = ATOL,
    rtol = RTOL,
)
    n = length(q)
    @assert n == LinearAlgebra.checksquare(Q)
    @assert n == size(A, 2)
    @assert n == size(G, 2)
    @assert length(fix_values) == length(fix_indices)
    model = DiffOpt.diff_optimizer(solver)
    MOI.set(model, MOI.Silent(), true)

    v = MOI.add_variables(model, n)

    _sign(x, a) = a == lt ? -x : x

    if lt
        cle = MOI.add_constraint.(model, G * v, MOI.LessThan.(h))
    else
        cle = MOI.add_constraint.(model, -G * v, MOI.GreaterThan.(-h))
    end
    if !iszero(nub)
        cub = MOI.add_constraint.(model, v[ub_indices], MOI.LessThan.(ub_values))
        G = vcat(G, sparse(1:nub, ub_indices, ones(nub), nub, n))
        h = vcat(h, ub_values)
    end
    if !iszero(nlb)
        clb = MOI.add_constraint.(model, v[lb_indices], MOI.GreaterThan.(lb_values))
        G = vcat(G, sparse(1:nlb, lb_indices, -ones(nlb), nlb, n))
        h = vcat(h, -lb_values)
    end

    ceq = MOI.add_constraint.(model, A * v, MOI.EqualTo.(b))
    if !iszero(nfix)
        cfix = MOI.add_constraint.(model, v[fix_indices], MOI.EqualTo.(fix_values))
        A = vcat(A, sparse(1:nfix, fix_indices, ones(nfix), nfix, n))
        b = vcat(b, fix_values)
    end

    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    if iszero(Q)
        obj = q ⋅ v
    else
        obj = dot(v, Q/2, v) + dot(q, v)
    end
    MOI.set(model, MOI.ObjectiveFunction{typeof(obj)}(), obj)

    MOI.optimize!(model)

    @_test(MOI.get(model, MOI.VariablePrimal(), v), z)

    # The sign as reversed as [AK17]
    # use a different convention for the dual
    _ν = -MOI.get.(model, MOI.ConstraintDual(), ceq)
    if !iszero(nfix)
        _ν = vcat(_ν, -MOI.get.(model, MOI.ConstraintDual(), cfix))
    end
    @_test(convert(Vector{Float64}, _ν), ν)
    _λ = _sign(MOI.get.(model, MOI.ConstraintDual(), cle), true)
    if !iszero(nub)
        _λ = vcat(_λ, -MOI.get.(model, MOI.ConstraintDual(), cub))
    end
    if !iszero(nlb)
        _λ = vcat(_λ, MOI.get.(model, MOI.ConstraintDual(), clb))
    end
    @_test(convert(Vector{Float64}, _λ), λ)

    #dobjb = v' * (dQb / 2.0) * v + dqb' * v
    # TODO, it should .-
    #dleb = dGb * v .+ dhb
    #deqb = dAb * v .+ dbb
    @assert dzb !== nothing
    @testset "Backward pass" begin
        MOI.set.(model, DiffOpt.ReverseVariablePrimal(), v, dzb)

        DiffOpt.reverse_differentiate!(model)

        dobjb = MOI.get(model, DiffOpt.ReverseObjectiveFunction())
        spb = DiffOpt.sparse_array_representation(
            DiffOpt.standard_form(dobjb),
            n,
            dobjb.index_map,
        )
        @_test(spb.quadratic_terms, dQb)
        @_test(spb.affine_terms, dqb)

        # FIXME should multiply by -1 if lt is false
        funcs = MOI.get.(model, DiffOpt.ReverseConstraintFunction(), cle)
        if !iszero(nub)
            funcs = vcat(funcs, MOI.get.(model, DiffOpt.ReverseConstraintFunction(), cub))
        end
        if !iszero(nlb)
            # FIXME should multiply by -1
            funcs = vcat(funcs, MOI.get.(model, DiffOpt.ReverseConstraintFunction(), clb))
        end
        @_test(convert(Vector{Float64}, _sign(MOI.constant.(funcs), true)), dhb)
        @_test(Float64[_sign(JuMP.coefficient(funcs[i], vi), false) for i in eachindex(funcs), vi in v], dGb)

        funcs = MOI.get.(model, DiffOpt.ReverseConstraintFunction(), ceq)
        if !iszero(nfix)
            funcs = vcat(funcs, MOI.get.(model, DiffOpt.ReverseConstraintFunction(), cfix))
        end
        @_test(convert(Vector{Float64}, -MOI.constant.(funcs)), dbb)
        @_test(Float64[JuMP.coefficient(funcs[i], vi) for i in eachindex(funcs), vi in v], dAb)
    end

    # Test against [AK17, eq. (8)]
    @_test(dqb, ∇zb)
    @_test((∇zb * z' + z * ∇zb') / 2, dQb)
    @_test(-dbb, ∇νb)
    @_test(∇νb * z' + ν * ∇zb', dAb)
    if all(i -> abs(λ[i]) > ATOL, 1:nle)
        @_test(-dhb ./ λ, ∇λb)
    end
    if ∇λb !== nothing
        @_test(Diagonal(λ) * ∇λb * z' + λ * ∇zb', dGb)
    end

    # Test against [AK17, eq. (7)]
    if ∇λb !== nothing
        @_test(-(Q * ∇zb + G' * (λ .* ∇λb) + A' * ∇νb), dzb)
    end
    if ∇λb !== nothing
        @_test(-(G * ∇zb + (G * z - h) .* ∇λb), dλb)
    end
    @_test(-A * ∇zb, dνb)

    dobjf = v' * (dQf / 2.0) * v + dqf' * v
    dlef = dGf * v .- dhf
    deqf = dAf * v .- dbf

    @testset "Forward pass" begin
        MOI.set(model, DiffOpt.ForwardObjectiveFunction(), dobjf)
        for (j, jc) in enumerate(cle)
            func = dlef[j]
            canonicalize && MOI.Utilities.canonicalize!(func)
            if set_zero || !MOI.iszero(dlef[j])
                MOI.set(model, DiffOpt.ForwardConstraintFunction(), jc, _sign(func, false))
            end
        end
        for (j, jc) in enumerate(ceq)
            func = deqf[j]
            canonicalize && MOI.Utilities.canonicalize!(func)
            if set_zero || !MOI.iszero(func)
                MOI.set(model, DiffOpt.ForwardConstraintFunction(), jc, func)
            end
        end
        if !iszero(nfix)
            for (j, jc) in enumerate(cfix)
                func = deqf[length(ceq)+j]
                canonicalize && MOI.Utilities.canonicalize!(func)
                if set_zero || !MOI.iszero(func)
                    # TODO FIXME should work if we drop support for `VariableIndex` and we let the Functionize bridge do the work
                    @test_throws MOI.UnsupportedAttribute MOI.set(model, DiffOpt.ForwardConstraintFunction(), jc, func)
                end
            end
        end

        DiffOpt.forward_differentiate!(model)

        @_test(MOI.get.(model, DiffOpt.ForwardVariablePrimal(), v), dzf)
    end

    # Test against [AK17, eq. (6)]
    @_test(dQf * z + dqf + dAf' * ν + dGf' * λ, ∇zf)
    if dλf !== nothing && dνf !== nothing
        @test Q * dzf + G' * dλf + A' * dνf ≈ ∇zf atol = ATOL rtol = RTOL
    end
    @_test(λ .* (dGf * z - dhf), ∇λf)
    @test (G * z - h) .* dλf + λ .* (G * dzf) ≈ -∇λf atol = ATOL rtol = RTOL
    @_test(dAf * z - dbf, ∇νf)
    @test A * dzf ≈ -∇νf atol = ATOL rtol = RTOL

    # As a kind of integration test, we check that the scalar product is the same whether it is don at the level of
    # 1) (dz, dλ, dν) (dλb and dνb are zero so we ignore their product (appropriate since we have not yet
    #    implemented the getter for dνf))
    dprod = dzf ⋅ dzb # ignored as it is zero : + dλf ⋅ dλb + dνf ⋅ dνb
    # 2) (∇z, ∇λ, ∇ν) which are the LHS of (6) and (7) (which are differentiation
    #    of the gradient of the laplacian with respect to z, λ and ∇ν hence the variable names)
    if ∇λb !== nothing
        ∇prod = ∇zf ⋅ ∇zb + ∇λf ⋅ ∇λb + ∇νf ⋅ ∇νb
        @test dprod ≈ ∇prod atol = ATOL rtol = RTOL
    end
    # 3) the problem data (here we made it so that they are the same for the forward
    #    and backward pass but we could have picked any other dQ, dq, ... for the forward pass
    #    and we would still have pprod = ∇prod = dprod
    pprod = dQf ⋅ dQb + dqf ⋅ dqb + dGf ⋅ dGb + dhf ⋅ dhb + dAf ⋅ dAb + dbf ⋅ dbb
    @test pprod ≈ pprod atol = ATOL rtol = RTOL
end

function qp_test(solver; kws...)
    @testset "With $(lt ? "LessThan" : "GreaterThan") constraints" for lt in [true, false]
        @testset "With$(set_zero ? "" : "out") setting zero tangents" for set_zero in [true, false]
            @testset "With$(canonicalize ? "" : "out") canonicalization" for canonicalize in [true, false]
                qp_test(solver, lt, set_zero, canonicalize; kws...)
            end
        end
    end
end

function qp_test_with_solutions(solver;
    dzb = nothing,
    n = length(dzb),
    q = nothing,
    dqf = zeros(n),
    Q = zeros(n, n),
    dQf = zeros(n, n),
    fix_indices = Int[],
    fix_values = Float64[],
    nfix = length(fix_values),
    ub_indices = Int[],
    ub_values = Float64[],
    nub = length(ub_values),
    lb_indices = Int[],
    lb_values = Float64[],
    nlb = length(lb_values),
    h = zeros(0),
    nle = length(h),
    dhf = zeros(nle + nub + nlb),
    G = zeros(0, n),
    dGf = zeros(nle + nub + nlb, n),
    b = zeros(0),
    neq = length(b),
    dbf = zeros(neq + nfix),
    A = zeros(0, n),
    dAf = zeros(neq + nfix, n),
    kws...
)
    @testset "Without known solutions" begin
        qp_test(solver;
            dzb = dzb,
            q = q,
            dqf = dqf,
            Q = Q,
            dQf = dQf,
            h = h,
            dhf = dhf,
            G = G,
            dGf = dGf,
            b = b,
            dbf = dbf,
            A = A,
            dAf = dAf,
            fix_indices = fix_indices,
            fix_values = fix_values,
            ub_indices = ub_indices,
            ub_values = ub_values,
            lb_indices = lb_indices,
            lb_values = lb_values,
        )
    end
    @testset "With known solutions" begin
        qp_test(solver;
            dzb = dzb,
            q = q,
            dqf = dqf,
            Q = Q,
            dQf = dQf,
            h = h,
            dhf = dhf,
            G = G,
            dGf = dGf,
            b = b,
            dbf = dbf,
            A = A,
            dAf = dAf,
            fix_indices = fix_indices,
            fix_values = fix_values,
            ub_indices = ub_indices,
            ub_values = ub_values,
            lb_indices = lb_indices,
            lb_values = lb_values,
            kws...,
        )
    end
end

ext: jl


text:
"""
# Compressed sensing JLab project code

This package provides three ways to form "single-pixel" images.

- `GhostImage`: Computational ghost imaging.
- `ExplicitCS`: Explicit compressive sensing (BPDN) by convex optimization.
- `ImplicitCS`: Implicit compressive sensing (BPDN) via the GenSPGL package.
"""
module CSJLab

using LinearAlgebra
using LinearAlgebra.BLAS
using DelimitedFiles
using Images.ImageCore
using FFTW
using JOLI
using GenSPGL

include("Utils.jl")
include("GhostImage.jl")
include("ImplicitCS.jl")

end # module


text:
module ConditionalComponent

println("\nRectifier: Demonstrating conditional components")

using Modia
using Modia.Electric
using ModiaMath: plot

@model ConditionalLoad begin
    R=Resistor(R=1)
    r=Resistor(R=0.01)
    C=Capacitor(C=1,start=0)
    D=IdealDiode()
    V=SineVoltage(V=5,freqHz=1.5, offset=0, startTime=0)
    extraLoad = true
    R2=if extraLoad; Resistor(R=1) end
@equations begin
    connect(V.p, D.p)
    connect(D.n, R.p)
    connect(R.n, r.p)
    connect(r.n, V.n)
    connect(C.n, R.n)
    connect(C.p, R.p)
    connect(C.n, R2.n)
    connect(C.p, R2.p)
    end
end

result = simulate(ConditionalLoad, 2, logTranslation=true)
plot(result, ("C.v", "V.v"), heading="Conditional", figure=12)

@model NoExtraLoad begin
    @extends ConditionalLoad(extraLoad = false)
end

result = simulate(NoExtraLoad, 2, logTranslation=true)
plot(result, ("C.v", "V.v"), heading="Conditional", figure=12)

end


text:
@testset "Folding" begin
  @testset "Uniform" begin
    d = CartesianGrid(100,100)
    f = folds(d, UniformFolding(10))
    for (source, target) in f
      @test length(source) == 9000
      @test length(target) == 1000
    end
    @test length(collect(f)) == 10
  end

  @testset "One" begin
    d = CartesianGrid(10,10)
    f = folds(d, OneFolding())
    for (source, target) in f
      @test length(source) == 99
      @test length(target) == 1
    end
    @test length(collect(f)) == 100
  end

  @testset "Block" begin
    d = CartesianGrid(100,100)
    f = folds(d, BlockFolding((10.,10.)))
    for (source, target) in f
      @test length(source) ∈ [9100,9400,9600]
      @test length(target) == 100
    end
    @test length(collect(f)) == 100
  end

  @testset "Ball" begin
    d = CartesianGrid(50,50)
    f = folds(d, BallFolding(MetricBall(10.)))
    @test length(collect(f)) == 2500
    ps = collect(f)
    ms = length.(first.(ps))
    ns = length.(last.(ps))
    @test all(2183 .≤ ms .≤ 2410)
    @test all(ns .== 1)
  end
end


text:
###### This file is part of the MendelImpute.jl package.
###### It contains relevant code to compress a reference haplotype panel into
###### a `.jlso` or `.jld2` compressed panel.

"""
Highest level of abstraction of a compressed object
"""
abstract type Compressed end

"""
Data structure for keeping track of unique haplotypes in a window.

# Every window contains:
- `uniqueindex`: the unique haplotype indices
- `to_unique`: vertor that maps every haplotype index to the corresponding
    unique haplotype index
- `uniqueH`: A BitMatrix storing unique haplotypes in columns
- `hapmap`: Dictionary that maps every unique haplotype to its equivalent
    haplotypes. Singleton haplotypes will not be recorded to conserve memory.

# Example:

Suppose in window `w` there are 9 haplotypes (represented by different letters):

    [a b b c b d a d a]

then

    uniqueindex = [1, 2, 4, 6]
    to_unique = [1, 2, 2, 3, 2, 4, 1, 4, 1]
    uniqueH = [a b c d]
    hapmap = `Dict{Int32,Array{Int32,1}} with 3 entry:
        1 => [1, 7, 9]
        2 => [2, 3, 5]
        6 => [6, 8]`

where `4 => [4]` in `hapmap` has been skipped since it is a singleton.
"""
struct CompressedWindow
    uniqueindex::Vector{Int32}
    hapmap::Dict{Int32, Vector{Int32}}
    to_unique::Vector{Int32}
    uniqueH::BitMatrix
end

"""
Keeps a vector of `CompressedWindow`. Indexing off instances of
`CompressedHaplotypes` means indexing off `CompressedHaplotypes.CW`

- `CW`: Vector of `CompressedWindow`. `CW[i]` stores unique haplotypes filtered
    with respect to all SNPs from `Hstart[i]` to `Hstart[i + 1]`
- `CW_typed`: Vector of `CompressedWindow`. `CW_typed[i]` stores unique
    haplotypes filtered with respect to typed SNPs from `Hstart[i]` to
    `Hstart[i + 1]`
- `Hstart`: `Hstart[i]` to `Hstart[i + 1]` is the range of H's SNPs (typed or
    untyped) that are in window `i`. It includes all SNP until the first
    typed snp of window `i + 1`. Used to directly index into `H` like
    `H[Hstart[i]:Hstart[i+1], :]`
- `X_window_range`: `idx = X_window_range[w]` is the range of `X`'s (typed
    SNP's) index in window `w`. Used to directly index into `X` like `X[idx, :]`.
    Does not include overlapping SNPs.
- `sampleID`: Sample names as listed in the VCF file
- `altfreq`: Alternate allele frequency (frequency of "1" in VCF file)
- `max_unique_haplotypes`: Maximum number of unique haplotypes in each window.
- `overlap`: Boolean indicating whether adjacent windows of typed SNPs overlap.
"""
struct CompressedHaplotypes
    CW::Vector{CompressedWindow}
    CW_typed::Vector{CompressedWindow}
    Hstart::Vector{Int} # we need the starting point more often than the range
    X_window_range::Vector{UnitRange{Int}}
    sampleID::Vector{String}
    chr::Vector{String}
    pos::Vector{Int}
    SNPid::Vector{Vector{String}}
    refallele::Vector{String}
    altallele::Vector{Vector{String}}
    altfreq::Vector{Float32}
    max_unique_haplotypes::Int
    overlap::Bool
end
CompressedHaplotypes(windows::Int, sampleID, chr, pos, SNPid, ref, alt, altfreq,
    max_unique_haplotypes, overlap) = CompressedHaplotypes(
    Vector{CompressedWindow}( undef, windows), Vector{CompressedWindow}(undef,
    windows), zeros(windows), [1:0 for _ in 1:windows], sampleID, chr, pos,
    SNPid, ref, alt, altfreq, max_unique_haplotypes, overlap)

nhaplotypes(x::CompressedHaplotypes) = 2length(x.sampleID)
nwindows(x::CompressedHaplotypes) = length(x.CW)

"""
    max_dim(reffile::String)

Searches each compressed window of the typed SNPs and outputs the maximum window
width (i.e. number of SNPs) and maximum unique haplotypes per window. If windows
are overlapping, extra padding will be automatically inserted
"""
function max_dim(reffile::String)
    endswith(reffile, ".jlso") || error("max_width can only" *
        " be called on a `CompressedHaplotypes` or `.jlso` files.")
    loaded = JLSO.load(reffile)
    compressed_Hunique = loaded[:compressed_Hunique]
    return max_dim(compressed_Hunique)
end
function max_dim(Hunique::CompressedHaplotypes)
    win = nwindows(Hunique)
    maxwidth = 0
    max_d = 0
    for w in 1:win
        n = size(Hunique.CW_typed[w].uniqueH, 1)
        n > maxwidth && (maxwidth = n)

        d = size(Hunique.CW_typed[w].uniqueH, 2)
        d > max_d && (max_d = d)
    end
    return maxwidth, max_d
end

"""
    get_window_widths(reffile::String)

Searches each compressed window of the typed SNPs and output their window widths
(i.e. number of SNPs per window)
"""
function get_window_widths(reffile::String)
    endswith(reffile, ".jlso") || error("max_width can only" *
        " be called on a `CompressedHaplotypes` or `.jlso` files.")
    loaded = JLSO.load(reffile)
    compressed_Hunique = loaded[:compressed_Hunique]
    return get_window_widths(compressed_Hunique)
end
function get_window_widths(Hunique::CompressedHaplotypes)
    win = nwindows(Hunique)
    widths = zeros(Int, win)
    for w in 1:win
        widths[w] = length(Hunique.X_window_range[w])
    end
    return widths
end

"""
    count_haplotypes_per_window(Hunique::CompressedHaplotypes)

Searches each compressed window of the typed SNPs and output number of unique
haplotypes per window
"""
function count_haplotypes_per_window(Hunique::CompressedHaplotypes)
    win = nwindows(Hunique)
    unique_haplotype_counts = zeros(Int, win)
    for w in 1:win
        unique_haplotype_counts[w] = length(Hunique.CW_typed[w].uniqueindex)
    end
    return unique_haplotype_counts
end
function count_haplotypes_per_window(reffile::String)
    endswith(reffile, ".jlso") || error("count_haplotypes_per_window can only" *
        " be called on a `CompressedHaplotypes` or `.jlso` files.")
    loaded = JLSO.load(reffile)
    compressed_Hunique = loaded[:compressed_Hunique]
    return count_haplotypes_per_window(compressed_Hunique)
end

"""
    avg_haplotypes_per_window(Hunique::CompressedHaplotypes)

Searches each compressed window of the typed SNPs and output average number of
unique haplotypes per window
"""
avg_haplotypes_per_window(Hunique::CompressedHaplotypes) =
    mean(count_haplotypes_per_window(Hunique))
avg_haplotypes_per_window(reffile::String) =
    mean(count_haplotypes_per_window(reffile))

"""
For an index in unique haplotype (of typed snps), finds the first occurance of
that haplotype in the complete reference pool for the specified window.
"""
function unique_idx_to_complete_idx(unique_idx, window,
    Hunique::CompressedHaplotypes)
    return Hunique.CW_typed[window].uniqueindex[unique_idx]
end

"""
For an index in unique haplotype (of all snps), finds the first occurance of
that haplotype in the complete reference pool for the specified window.
"""
function unique_all_idx_to_complete_idx(unique_idx, window,
    Hunique::CompressedHaplotypes)
    return Hunique.CW[window].uniqueindex[unique_idx]
end

"""
For an index in the complete haplotype pool, find its index in the unique
haplotype pool of all SNPs (typed + untyped) in specified window.
"""
function complete_idx_to_unique_all_idx(complete_idx, window,
    Hunique::CompressedHaplotypes)
    return Hunique.CW[window].to_unique[complete_idx]
end

"""
For an index in the complete haplotype pool, find its index in the unique
haplotype pool of just the typed SNPs in specified window.
"""
function complete_idx_to_unique_typed_idx(complete_idx, window,
    Hunique::CompressedHaplotypes)
    return Hunique.CW_typed[window].to_unique[complete_idx]
end

"""
    extend_to_overlap_range(Hunique::CompressedHaplotypes, window::Int, overlap::Bool)

If genotype window overlap, the unique haplotype matrix will be larger than the
genotype matrix in window `w`. This function computes the range of genotype
matrix so that the ranges match up.
"""
function extend_to_overlap_range(
    Hunique::CompressedHaplotypes,
    window::Int,
    overlap::Bool
    )
    winranges = Hunique.X_window_range
    Hw = Hunique.CW_typed[window].uniqueH
    Xrange = winranges[window]
    if overlap
        extra_width = size(Hw, 1) - length(winranges[window])
        if window == 1 # first window
            Xstart = first(winranges[window])
            Xend = last(winranges[window]) + extra_width
        elseif window == length(winranges) # last window
            Xstart = first(winranges[window]) - extra_width
            Xend = last(winranges[window])
        else
            Xstart = first(winranges[window]) - (extra_width >> 1)
            Xend = last(winranges[window]) + (extra_width >> 1)
        end
        Xrange = Xstart:Xend
    end
    return Xrange
end

"""
    nonoverlap_range(Hunique::CompressedHaplotypes, window::Int, overlap::Bool)

If windows overlap, the unique haplotype matrix will include extra SNPs on
both ends as paddings. This function calculates the range of SNPs of the unique
haplotype matrix that are not in the overlapping regions.
"""
function nonoverlap_range(
    Hunique::CompressedHaplotypes,
    window::Int
    )
    winranges = Hunique.X_window_range
    Hw = Hunique.CW_typed[window].uniqueH
    Hw_width = size(Hw, 1)
    extra_width = Hw_width - length(winranges[window])

    if window == 1 # first window
        Hstart = 1
        Hend = Hw_width - extra_width
    elseif window == length(winranges) # last window
        Hstart = 1 + extra_width
        Hend = Hw_width
    else
        Hstart = 1 + (extra_width >> 1)
        Hend = Hw_width - (extra_width >> 1)
    end

    return Hstart:Hend
end

"""
    compress_haplotypes(reffile::String, tgtfile::String, outfile::String,
        [d::Int], [minwidth::Int], [overlap::Float64])

Cuts a haplotype matrix `reffile` into windows of variable width so that each
window has less than `d` unique haplotypes. Saves result to `outfile` as
a compressed binary format. All SNPs in `tgtfile` must be present in `reffile`.

# Why is `tgtfile` required?
The unique haplotypes in each window is computed on the typed SNPs only.
A genotype matrix `tgtfile` is used to identify the typed SNPs. In the future,
hopefully we can pre-compute compressed haplotype panels for all genotyping
platforms and provide them as downloadable files. But currently, users must
run this function by themselves.

# Inputs
* `reffile`: reference haplotype file name (ends in `.vcf` or `.vcf.gz`)
* `tgtfile`: target genotype file name (ends in `.vcf` or `.vcf.gz`)
* `outfile`: Output file name (ends in `.jlso`)

# Optional Inputs
* `d`: Max number of unique haplotypes per genotype window (default `d = 1000`).
* `minwidth`: Minimum number of typed SNPs per window (default 0)
* `overlap`: How much overlap between adjacent genotype windows in percentage of
    each window's width (default 0.0)
"""
function compress_haplotypes(
    reffile::AbstractString,
    tgtfile::AbstractString,
    outfile::AbstractString,
    d::Int=1000,
    minwidth::Int=0,
    overlap::Float64=0.0
    )
    endswith(outfile, ".jld2") || endswith(outfile, ".jlso") ||
        error("Unrecognized compression format: `outfile` can only end in " *
        "`.jlso` or `.jld2`")
    0.0 ≤ overlap ≤ 1.0 || error("overlap must be a percentage")

    # import reference haplotypes
    H, H_sampleID, H_chr, H_pos, H_ids, H_ref, H_alt = convert_ht(Bool,
        reffile, trans=true, save_snp_info=true,
        msg="importing reference data...")
    X, X_sampleID, X_chr, X_pos, X_ids, X_ref, X_alt =
        VCFTools.convert_gt(UInt8, tgtfile, trans=true, save_snp_info=true,
        msg = "Importing genotype file...")
    any(isnothing, indexin(X_pos, H_pos)) && error("Found SNPs in target " *
        "file that are not in reference file! Please filter them out first!")

    # compress routine
    compress_haplotypes(H, X, outfile, X_pos, H_sampleID, H_chr, H_pos, H_ids,
        H_ref, H_alt, d, minwidth, overlap)

    return nothing
end

function compress_haplotypes(H::AbstractMatrix, X::AbstractMatrix,
    outfile::AbstractString, X_pos::AbstractVector, H_sampleID::AbstractVector,
    H_chr::AbstractVector, H_pos::AbstractVector, H_ids::AbstractVector,
    H_ref::AbstractVector, H_alt::AbstractVector, d::Int, minwidth::Int,
    overlap::Float64)

    endswith(outfile, ".jld2") || endswith(outfile, ".jlso") ||
        error("Unrecognized compression format: `outfile` can only end in " *
        "`.jlso` or `.jld2`")
    0.0 ≤ overlap ≤ 1.0 || error("overlap must be a percentage")

    # some constants
    ref_snps = size(H, 1)
    tgt_snps = size(X, 1)

    # compute window intervals based on typed SNPs
    XtoH_idx = indexin(X_pos, H_pos) # assumes all SNPs in X are in H
    any(isnothing, XtoH_idx) && error("Found SNPs in target " *
        "file that are not in reference file! Please filter them out first!")
    Hw_typed = H[XtoH_idx, :]        # H with only typed snps
    window_ranges = get_window_intervals(Hw_typed, d, minwidth)
    wins = length(window_ranges)

    # initialize compressed haplotype object
    alt_freq = reshape(sum(H, dims=2), size(H, 1)) ./ size(H, 2)
    compressed_Hunique = CompressedHaplotypes(wins, H_sampleID,
        H_chr, H_pos, H_ids, H_ref, H_alt, convert(Vector{Float32},alt_freq),
        d, overlap > 0.0)
    compressed_Hunique.X_window_range .= window_ranges

    # record unique haplotypes and mappings window by window
    Hw_idx_start = 1
    for w in 1:length(window_ranges)
        # genotype matrix's current window ranges
        Xw_idx_start = first(window_ranges[w])
        Xw_idx_end = last(window_ranges[w])
        if overlap > 0.0 # overlaps between adjacent windows
            flankwidth = round(Int, overlap*(Xw_idx_end-Xw_idx_start), RoundUp)
            w == 1 || (Xw_idx_start -= flankwidth)
            w == wins || (Xw_idx_end += flankwidth)
        end

        # haplotype reference panel's window ranges
        Xw_pos_next_start = (w == wins ? X_pos[end] :
            X_pos[first(window_ranges[w + 1])])
        Hw_idx_end = (w == wins ? length(H_pos) :
            something(findnext(x -> x == Xw_pos_next_start, H_pos,
            Hw_idx_start)) - 1)
        compressed_Hunique.Hstart[w] = Hw_idx_start

        # get current window of H
        Xw_pos = @view(X_pos[Xw_idx_start:Xw_idx_end])
        XwtoH_idx = indexin(Xw_pos, H_pos) # assumes all SNPs in X are in H
        Hw = @view(H[Hw_idx_start:Hw_idx_end, :]) # including all snps
        Hw_typed = @view(H[XwtoH_idx, :])         # including only typed snps

        # find unique haplotypes on all SNPs
        mapping = groupslices(Hw, dims = 2)
        unique_idx = unique(mapping)
        hapmap = Dict{Int32, Vector{Int32}}()
        for idx in unique_idx
            redundant_haplotypes = findall(x -> x == idx, mapping)
            if length(redundant_haplotypes) == 1
                continue # skip singletons
            else
                hapmap[idx] = redundant_haplotypes
            end
        end
        complete_to_unique = indexin(mapping, unique_idx)
        uniqueH = Hw[:, unique_idx]
        compressed_Hunique.CW[w] = MendelImpute.CompressedWindow(unique_idx,
            hapmap, complete_to_unique, uniqueH)

        # find unique haplotypes on typed SNPs
        mapping_typed = groupslices(Hw_typed, dims = 2)
        unique_idx_typed = unique(mapping_typed)
        hapmap_typed = Dict{Int32, Vector{Int32}}()
        for idx in unique_idx_typed
            redundant_haplotypes = findall(x -> x == idx, mapping_typed)
            if length(redundant_haplotypes) == 1
                continue # skip singletons
            else
                hapmap_typed[idx] = redundant_haplotypes
            end
        end
        complete_to_unique_typed = indexin(mapping_typed, unique_idx_typed)
        uniqueH_typed = Hw_typed[:, unique_idx_typed]
        compressed_Hunique.CW_typed[w] =
            MendelImpute.CompressedWindow(unique_idx_typed, hapmap_typed,
            complete_to_unique_typed, uniqueH_typed)

        # update Hw_idx_start
        Hw_idx_start = Hw_idx_end + 1
    end

    # save using JLSO or JLD2
    if endswith(outfile, ".jld2")
        jldopen(outfile, "w", compress=true) do file
            file["compressed_Hunique"] = compressed_Hunique
        end
    else #jlso
        JLSO.save(outfile, :compressed_Hunique =>
            compressed_Hunique, format=:julia_serialize, compression=:gzip)
    end

    return compressed_Hunique
end

"""
    get_window_intervals(H, d, minwidth, low, high, intervals, seen)

Find window intervals so that every window has `d` or less unique haplotypes
by recursively dividing windows into halves.

# Inputs
- `H`: The full haplotype matrix (on typed SNPs). Each column is a haplotype.
- `d`: Number of unique haplotypes in each window
- `minwidth`: Minimum number of SNPs per window (default 20)
- `low`: start of current window
- `high`: end of current window
- `intervals`: Vector of window ranges. This is also the return vector
- `seen`: storage container

# Output
- `intervals`: Vector of window ranges. Each range have ≤ `d` unique haplotypes
"""
function get_window_intervals(
    H::AbstractMatrix,
    d::Int,
    minwidth::Int=20,
    low::Int=1,
    high::Int=size(H, 1),
    intervals = UnitRange[],
    seen=BitSet(),
    )

    mid = (low + high) >> 1
    unique_columns_maps = groupslices(view(H, low:high, :), dims = 2)
    k = count_unique(unique_columns_maps, seen)
    if k ≤ d || length(low:mid) ≤ minwidth
        push!(intervals, low:high)
    else
        get_window_intervals(H, d, minwidth, low,     mid,  intervals, seen)
        get_window_intervals(H, d, minwidth, mid + 1, high, intervals, seen)
    end

    return sort!(intervals)
end

"""
    count_unique(v, seen)

Count the number of unique elements in `v`, using `seen` as storage.
"""
function count_unique(v::AbstractVector{<:Integer}, seen::AbstractSet=BitSet())
    empty!(seen)
    s = 0
    for i in v
        if i ∉ seen
            s += 1
            push!(seen, i)
        end
    end
    return s
end


text:
using CSV
using DataFrames
using DelimitedFiles
using YAML
using Query

"""
    share_utd!(d::Dict, set::Dict)
`utd`: Share of total trade by region.
"""
function share_utd!(d::Dict, set::Dict)
    println("  Calculating share of total trade by region")
    d[:utd] = fill_zero(d[:utd])
    df = d[:utd] / transform_over(d[:utd], :r)

    df_yr = transform_over(d[:utd], :yr) / transform_over(d[:utd], [:yr,:r])
    df[isnan.(df[:,:value]), :value] .= df_yr[isnan.(df[:,:value]),:value]

    # Check import and export shares.
    verify_over(filter_with(df, (t = "imports",)), :r) !== true && @error("Import shares don't sum to 1.")
    verify_over(filter_with(df, (t = "exports",)), :r) !== true && @error("Export shares don't sum to 1.")

    d[:utd] = dropnan(dropzero(df))
    set[:notrd] = setdiff(set[:s], d[:utd][:,:s])
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
module FileHelpers
import Pluto: is_pluto_notebook

flatmap(args...) = vcat(map(args...)...)

list_files_recursive(dir=".") = let
    paths = flatmap(walkdir(dir)) do (root, dirs, files)
        joinpath.([root], files)
    end
    relpath.(paths, [dir])
end

"""
Search recursively for Pluto notebook files.

Return paths relative to the search directory.
"""
function find_notebook_files_recursive(start_dir)
    notebook_files = filter(is_pluto_notebook, list_files_recursive(start_dir))

    not_interal_notebook_files = filter(notebook_files) do f
        !occursin(".julia", f) || occursin(".julia", start_dir)
    end

    # reverse alphabetical order so that week5 becomes available before week4 :)
    reverse(not_interal_notebook_files)
end
end


text:
using Plots,ApproxFun,SingularIntegralEquations;  gadfly()

##
#  Ideal fluid flow consists of level sets of the imagainary part of a function
# that is asymptotic to c*z and whose imaginary part vanishes on Γ
#
#
# On the unit interval, -2*hilbert(ui) gives the imaginary part of cauchy(ui,z)
#  So if we want to find ui defined on Γ so that hilbert(ui) = imag(c*z)
#  then c*z + 2cauchy(u,z) vanishes on Γ
##

u(x,y)=α*(x+im*y)+2cauchy(ui,x+im*y)

m=80;x = linspace(-2.,2.,m);y = linspace(-1.,1.,m+1)
    xx,yy = x.+0.*y',0.*x.+y'

c

ui|>space
Number(c)+hilbert(ui)(0.25+0.125im)

imag(α*(0.25+0.125im))

Hilbert()*ui
x=Fun()
f=Fun(exp)*sqrt(1-x^2)
cauchy(f,.1+.000001im)
2cauchy(f,.1,+)
f(.1)-im*hilbert(f,.1)

d=Circle()
z=Fun(d)
f=exp(real(z))
f(exp(im*0.1))

2imag(cauchy(f,exp(im*0.1),+))+real(hilbert(f,exp(im*0.1)))

im*sum(f/z)
linesum(f)

∪

k=50
    Γ=Interval(0.,1+0.5im)
    z=Fun(Γ)
    α=exp(-π*k/50im)
    c,ui=[1 Hilbert()]\imag(α*z)
    plot(Γ)
    contour!(x,y,imag(u(xx,yy))).o

##
# On an arc, the Hilbert transform no longer gives the imaginary part
#  However, pseudohilbert gives the imaginary part of pseudocauchy
#  So if we want to find ui defined on Γ so that pseudohilbert(ui) = imag(c*z)
#  then c*z + 2pseudocauchy(u,z) vanishes on Γ
##

u(x,y)=α*(x+im*y)+2pseudocauchy(ui,x+im*y)

m=80;x = linspace(-2.,2.,m);y = linspace(-2.,2.,m+1)
    xx,yy = x.+0.*y',0.*x.+y'

k=227;
    Γ=0.5+exp(im*Interval(0.1,-42))
    z=Fun(Γ)
    α=exp(-k/50im)
    c,ui=[1 PseudoHilbert()]\imag(α*z)
    plot(Γ)
    contour!(x,y,imag(u(xx,yy))).o

##
#  Circle
##

Γ=Circle()
z=Fun(Fourier(Γ))

u(x,y)=α*(x+im*y)+2cauchy(ui,x+im*y)

m=80;x = linspace(-2.,2.,m);y = linspace(-2.,2.,m+1)
    xx,yy = x.+0.*y',0.*x.+y'

k=239;
    α=exp(-k/45im)
    c,ui=[0 BasisFunctional(1);
          1 real(Hilbert())]\[0.,imag(α*z)]

    plot(Γ)
    contour!(x,y,imag(u(xx,yy))).o

##
# On a curve, the Hilbert transform may be complex, so we
# take the real part
##

u(x,y)=α*(x+im*y)+2cauchy(ui,x+im*y)

m=80;x = linspace(-2.,2.,m);y = linspace(-2.,2.,m+1)
    xx,yy = x.+0.*y',0.*x.+y'

Γ=Curve(Fun(x->exp(0.8im)*(x+x^2-1+im*(x-4x^3+x^4)/6)))
    z=Fun(Γ)
    α=im
    c,ui=[1 real(Hilbert())]\imag(α*z)
    plot(Γ)
    contour!(x,y,imag(u(xx,yy))).o

##
#  Two intervals requires explicitely stating the space (for now)
##

Γ=Interval(-1.,-0.5)∪Interval(-0.3,1.)
z=Fun(Γ)

ds=PiecewiseSpace(map(d->JacobiWeight(0.5,0.5,Ultraspherical(1,d)),Γ.domains))

u(x,y)=α*(x+im*y)+2cauchy(ui,x+im*y)

m=80;x = linspace(-2.,2.,m);y = linspace(-1.,1.,m+1)
    xx,yy = x.+0.*y',0.*x.+y'

k=114;
    α=exp(k/50*im)
    a,b,ui=[ones(Γ[1]) ones(Γ[2]) Hilbert(ds)]\imag(α*z)
    plot(Γ)
    contour!(x,y,imag(u(xx,yy))).o

Γ=Interval(-1.,0.)∪Interval(0.5im,1.)
z=Fun(Γ)
ds=PiecewiseSpace(map(d->JacobiWeight(0.5,0.5,Ultraspherical(1,d)),Γ.domains))

u(x,y)=α*(x+im*y)+2cauchy(ui,x+im*y)

m=80;x = linspace(-2.,2.,m);y = linspace(-1.,1.,m+1)
    xx,yy = x.+0.*y',0.*x.+y'

k=114;
    α=exp(k/50*im)
    a,b,ui=[ones(Γ[1]) ones(Γ[2]) real(Hilbert(ds))]\imag(α*z)
    plot(Γ)
    contour!(x,y,imag(u(xx,yy)))

Γ=Interval(-im,1.0-im)∪Curve(Fun(x->exp(0.8im)*(x+x^2-1+im*(x-4x^3+x^4)/6)))
z=Fun(Γ)

ds=PiecewiseSpace([JacobiWeight(0.5,0.5,Ultraspherical(1,Γ[1])),
                       MappedSpace(Γ[2],JacobiWeight(0.5,0.5,Ultraspherical(1)))])

m=80;x = linspace(-2.,2.,m);y = linspace(-3.,2.,m+1)
    xx,yy = x.+0.*y',0.*x.+y'

k=114;
    α=exp(k/50*im)
    a,b,ui=[ones(Γ[1]) ones(Γ[2]) real(Hilbert(ds))]\imag(α*z)
Gadfly.plot(ApproxFun.layer(Γ),
                layer(x=x,y=y,z=imag(u(xx,yy)),Main.Gadfly.Geom.contour(levels=[-3.:.05:4.])))

## Channel

# TODO: rename
using ApproxFun,SingularIntegralEquations,SO,Plots
    import ApproxFun:RealUnivariateSpace,ComplexBasis,canonicalspace,bandinds,addentries!,conversion_rule,
            spacescompatible,toeplitz_addentries!,domainscompatible,rangespace, ConcreteConversion,
            coefficients
    import SingularIntegralEquations:cauchy,ConcreteHilbert, stieltjes
immutable FourierDirichlet{DD} <: RealUnivariateSpace{DD}
    domain::DD
end

FourierDirichlet()=FourierDirichlet(PeriodicInterval())

    spacescompatible(a::FourierDirichlet,b::FourierDirichlet)=domainscompatible(a,b)

    canonicalspace(S::FourierDirichlet)=Fourier(domain(S))

    bandinds{DD}(::ConcreteConversion{FourierDirichlet{DD},Fourier{DD}})=-1,1
    function addentries!{DD}(C::ConcreteConversion{FourierDirichlet{DD},Fourier{DD}},A,kr::Range,::Colon)
        toeplitz_addentries!([1.],[0.,1.],A,kr)
    end

    conversion_rule{DD}(b::FourierDirichlet{DD},a::Fourier{DD})=b

    setdomain(F::FourierDirichlet,d::Domain)=FourierDirichlet(d)

    function coefficients(v::Vector,a::FourierDirichlet,b::Fourier)
        ret=zeros(typeof(v),length(v)+1)
        ret[2]=v[1]
        for k=2:length(v)
            ret[k-1]+=v[k]
            ret[k+1]+=v[k]
        end
        ret
    end

bandinds{L<:PeriodicLine}(H::ConcreteHilbert{FourierDirichlet{L}})=-2,2
    rangespace{L<:PeriodicLine}(H::ConcreteHilbert{FourierDirichlet{L}})=Fourier(domain(H))

    function addentries!{L<:PeriodicLine}(H::ConcreteHilbert{FourierDirichlet{L}},A,kr::Range,::Colon)
        if 1 in kr
            A[1,1]+=1
        end
        addentries!(real(Hilbert(FourierDirichlet(Circle()))),A,kr,:)
        A
    end

stieltjes{DD}(sp::FourierDirichlet{DD},cfs,z)=stieltjes(Fun(Fun(f,sp),Fourier),z)
    stieltjes{PL}(f::Fun{Fourier{PL}},z::Vector)=Complex128[stieltjes(f,zk) for zk in z]
    stieltjes{PL}(f::Fun{Fourier{PL}},z::Matrix)=reshape(stieltjes(f,vec(z)),size(z,1),size(z,2))::Matrix{Complex128}
    stieltjes{PL}(f::Fun{FourierDirichlet{PL}},z)=stieltjes(Fun(f,Fourier),z)

S=FourierDirichlet(PeriodicLine()-1.0im)
PS=PiecewiseSpace([S,JacobiWeight(0.5,0.5,Ultraspherical(1,Interval(-0.5im,1.+1.1im))),
                    JacobiWeight(0.5,0.5,Ultraspherical(1,Interval(-0.5,-0.5+.5im)))])

H=Hilbert(PS)

f=Fun(z->imag(z+1.0im),rangespace(H))
Γ=domain(PS)

a,b,ui=linsolve([ones(Γ[2]) ones(Γ[3]) real(H)],f;tolerance=1E-10)
m=80;x = linspace(-2.,3.,m);y = linspace(-1.,3.,m+1)
    xx,yy = x.+0.*y',0.*x.+y'

u(x,y)=(x+im*y)+2cauchy(ui,x+im*y)
    vals=imag(u(xx,yy))

pyplot()
plot(Γ[2];color=:Blue,legend=false)
    plot!(Γ[3];color=:Blue)
    plot!(Interval(-2.-im,3.0-im);color=:Blue)
    contour!(x,y,vals;ylims=(-2.,3.),levels=-1.:.1:3.)

Main.PyPlot.savefig("halfg.png";dpi=200,bbox_inches="tight")

png("halfplane")

gadfly()
pdf("halfplane")

Gadfly.plot(ApproxFun.layer(Γ[2]),

            ApproxFun.layer(Interval(-2.-im,3.0-im)),
                layer(x=x,y=y,z=imag(u(xx,yy)),Main.Gadfly.Geom.contour(levels=[-1.:0.05:2.])))

pieces(ui)[2]|>ApproxFun.plot
g=pieces(ui)[1]
    g=Fun(g.coefficients,space(g).space)
    g=Fun(g,Fourier)

g|>ApproxFun.plot

S1=MappedSpace(PeriodicLine()-1.0im,FourierDirichlet())
S2=MappedSpace(PeriodicLine()+1.0im,FourierDirichlet())

PS=PiecewiseSpace([S1,S2])
H=Hilbert(PS)

using SO

ds=PS[1]
rs=MappedSpace(domain(PS[2]),Fourier())
b=Fun([zeros(18);1.],ds);
    Fun(x->-stieltjes(b,x)/π,rs).coefficients|>chopm

g=Fun(b,MappedSpace(domain(PS[1]),Fourier()))
@which cauchy(b,5.0-1.0im)
cauchy(g,5.0-1.0im)
Fun(x->-stieltjes(b,x)/π,rs,20)

v1=chop(Fun(x->-stieltjes(b,x)/π,rs).coefficients,tol)

f=Fun(z->imag(z+1.0im),rangespace(H))
Γ=domain(PS)

a,ui=linsolve([ones(Γ[2]) real(H)],f;tolerance=1E-5)
m=80;x = linspace(-2.,3.,m);y = linspace(-1.5,2.,m+1)
    xx,yy = x.+0.*y',0.*x.+y'

u(x,y)=(x+im*y)+2cauchy(ui,x+im*y)
Gadfly.plot(ApproxFun.layer(Γ[2]),
            ApproxFun.layer(Interval(-2.-im,3.0-im)),
                layer(x=x,y=y,z=imag(u(xx,yy)),Main.Gadfly.Geom.contour(levels=[-1.:0.05:2.])))

d=Circle()
f=Fun(exp,d)
g=Fun(cos,d)

s=z->imag(conj(z)*cauchy(f,z)+cauchy(g,z))
pyplot()
x=y=linspace(-1.1,1.1,40)
z=Float64[s(x+im*y) for x in x, y in y]

s(0.999999999exp(im*1.))-imag(exp(-im)*f(exp(im))+g(exp(im)))
s(0.999999999exp(im*1.))
exp(dual(0.99999999999exp(im*1.),exp(im*1.)))
z=dual(0.99999999999exp(im*1.),exp(im*1.))
conj(z(&exp(

exp(im*1.)*exp(exp(im*1.))
z=Fun(d)
Multiplication(Fun([0.,0.,1.],Fourier(d)),Fourier(d))
Multiplication(conj(z),Laurent(d))
imag(f(z) +  conj(z)*g(z))

g(z)==r(z)/2
f(z)==q(z)-z/2*r(z)

imag(q(z) +  (z+conj(z))/2*r(z))

imag(exp(z))
imag(conj(z)(1/abs(z)*f(z)      +z*f'(z)))+z*g(z)) == 0

realpart(imag(z*exp(z))           )
realpart(imag(conj(z)(1/abs(z)*f(z)      +z*f'(z)))+z*g(z))        )
dualpart(imag(conj(z)*exp(z))          )-(realpart(imag(conj(z)/abs(z)*exp(z)           + conj(z)*z*exp(z))))

conj(z)

s(dual(0.99999999999exp(im*1.),exp(im*1.)))
dualpart(cos((10E5)*dual(1,2)))+2*10E5*sin(10E5)
2exp(1.)
expm([1 1 1/2;
      0 1 1;
      0 0 1])
cauchy(f,exp(im),+)
cauchy(g,exp(im),+)
gui()
surface(x,y,z)
s(.1+.2im)

real(Hilbert())  # Im g(z)

max_stars_repo_path: examples/Ideal Fluid Flow.jl

max_issues_repo_path: examples/Ideal Fluid Flow.jl

max_forks_repo_path: examples/Ideal Fluid Flow.jl


text:
using LinearAlgebra
using LoopVectorization
using Octavian

function RCCSDpT(rhf::RHF, Alg::ijk2)
    aoints = IntegralHelper()
    moints = IntegralHelper(orbitals=rhf.orbitals)
    ccsd = RCCSD(moints, aoints)
    return RCCSDpT(ccsd, moints, Alg)
end

function RCCSDpT(Alg::ijk2)
    aoints = IntegralHelper{Float64}()
    rhf = RHF(aoints)
    moints = IntegralHelper(orbitals=rhf.orbitals)
    ccsd = RCCSD(moints, aoints)
    return RCCSDpT(ccsd, moints, Alg)
end

function RCCSDpT(ccsd::RCCSD, moints::IntegralHelper{T,E,O}, Alg::ijk2) where {T<:AbstractFloat, E<:AbstractERI, O<:AbstractRestrictedOrbitals}

    output("\n   • Perturbative Triples Started\n")
    output("   - Contraction Engine: Octavian")

    T1 = permutedims(ccsd.T1.data, (2,1))
    T2 = permutedims(ccsd.T2.data, (4,3,2,1))
    v,o = size(T1)

    # Switch to physicist notation for better memory layout
    ovoo = permutedims(moints["OOOV"].data, (1,4,2,3))
    vvoo = permutedims(moints["OVOV"].data, (2,4,1,3))
    vvvo = reshape(permutedims(moints["OVVV"].data, (3,2,4,1)), v*v, v, o)

    fo = moints["Fii"].data
    fv = moints["Faa"].data

    # Pre-allocate Intermediate arrays
    Ws  = [Array{T}(undef, v,v,v) for _ = 1:Threads.nthreads()]
    Vs  = [Array{T}(undef, v,v,v) for _ = 1:Threads.nthreads()]
    tVs  = [Array{T}(undef, v,v,v) for _ = 1:Threads.nthreads()]
    Evals = zeros(T, Threads.nthreads())

    output("Computing energy contribution from occupied orbitals:")
    BLAS_THREADS = BLAS.get_num_threads()
    BLAS.set_num_threads(1)
    TBLIS_THREADS = TBLIS.set_num_threads(1)

    t = @elapsed begin
    Threads.@threads for i in 1:o
    @inbounds begin
        id = Threads.threadid()
        W = Ws[id]
        V = Vs[id]
        tV = tVs[id]

        # T1 views → V
        @views T1_i   = T1[:,i]

        # T2 views → W
        @views T2_i   = reshape(T2[:,:,:,i], v*v, o)

        @views vvvo_i = vvvo[:,:,i]

        for j in 1:i
            Dij = fo[i] + fo[j]

            @views vvvo_j = vvvo[:,:,j]

            # T1 views → V
            @views T1_j  = T1[:,j]

            # T2 views → W
            @views T2_j   = reshape(T2[:,:,:,j], v*v, o)
            ## Twice
            @views T2_ji = T2[:,:,j,i]
            @views T2_ij = T2[:,:,i,j]

            @views vvoo_ij = vvoo[:,:,i,j]

            @views ovoo_ji = ovoo[:,:,j,i]
            @views ovoo_ij = ovoo[:,:,i,j]

            δij = i == j
            for k in 1:j
                Dijk = Dij + fo[k]

                # T1 views → V
                @views T1_k = T1[:,k]
                @views vvvo_k = vvvo[:,:,k]

                # T2 views → W
                @views T2_k   = reshape(T2[:,:,:,k], v*v, o)
                ## Twice
                @views T2_jk = T2[:,:,j,k]
                @views T2_ik = T2[:,:,i,k]
                @views T2_ki = T2[:,:,k,i]
                @views T2_kj = T2[:,:,k,j]

                @views vvoo_jk = vvoo[:,:,j,k]
                @views vvoo_ik = vvoo[:,:,i,k]
                @views ovoo_ik = ovoo[:,:,i,k]
                @views ovoo_jk = ovoo[:,:,j,k]
                @views ovoo_ki = ovoo[:,:,k,i]
                @views ovoo_kj = ovoo[:,:,k,j]

                # Build intermediates W and V for the set i j k

                # abc
                W = reshape(W, v^2, v)
                matmul_serial!(W, vvvo_j, T2_ik)
                matmul_serial!(W, T2_j, ovoo_ik, -1.0, 1.0)
                W = reshape(W, v,v,v)

                # acb
                V = reshape(V, v^2, v)
                matmul_serial!(V, vvvo_k, T2_ij)
                matmul_serial!(V, T2_k, ovoo_ij, -1.0, 1.0)
                V = reshape(V, v,v,v)
                permutedims!(tV, V, (1,3,2))
                W .+= tV

                # bac
                V = reshape(V, v^2, v)
                matmul_serial!(V, vvvo_i, T2_jk)
                matmul_serial!(V, T2_i, ovoo_jk, -1.0, 1.0)
                V = reshape(V, v,v,v)
                permutedims!(tV, V, (2,1,3))
                W .+= tV

                # bca
                V = reshape(V, v^2, v)
                matmul_serial!(V, vvvo_k, T2_ji)
                matmul_serial!(V, T2_k, ovoo_ji, -1.0, 1.0)
                V = reshape(V, v,v,v)
                permutedims!(tV, V, (3,1,2))
                W .+= tV

                # cab
                V = reshape(V, v^2, v)
                matmul_serial!(V, vvvo_i, T2_kj)
                matmul_serial!(V, T2_i, ovoo_kj, -1.0, 1.0)
                V = reshape(V, v,v,v)
                permutedims!(tV, V, (2,3,1))
                W .+= tV

                # cba
                V = reshape(V, v^2, v)
                matmul_serial!(V, vvvo_j, T2_ki)
                matmul_serial!(V, T2_j, ovoo_ki, -1.0, 1.0)
                V = reshape(V, v,v,v)
                permutedims!(tV, V, (3,2,1))
                W .+= tV

                @tensor begin
                    V[a,b,c]  = W[a,b,c] + T1_i[a]*vvoo_jk[b,c] + vvoo_ik[a,c]*T1_j[b] + vvoo_ij[a,b]*T1_k[c]
                end

                # Compute Energy contribution
                δjk = j === k
                for a in 1:v
                    Dijka = Dijk - fv[a]
                    for b in 1:a
                        Dijkab = Dijka - fv[b]
                        δab = a === b
                        @fastmath for c in 1:b
                            Dd = Dijkab - fv[c]
                            δbc = b === c
                            X = (W[a,b,c]*V[a,b,c] + W[a,c,b]*V[a,c,b] + W[b,a,c]*V[b,a,c] + W[b,c,a]*V[b,c,a] + W[c,a,b]*V[c,a,b] + W[c,b,a]*V[c,b,a])
                            Y = (V[a,b,c] + V[b,c,a] + V[c,a,b])
                            Z = (V[a,c,b] + V[b,a,c] + V[c,b,a])
                            Ef = (Y - 2*Z)*(W[a,b,c] + W[b,c,a] + W[c,a,b]) + (Z - 2*Y)*(W[a,c,b]+W[b,a,c]+W[c,b,a]) + 3*X
                            Evals[id] += Ef*(2 - δij - δjk) / (Dd*(1 + δab + δbc))
                        end
                    end
                end
            end
        end
        output("  Orbital {} ✔", i)
    end # Threads (i loop)
    end # inbounds
    end # time

    BLAS.set_num_threads(BLAS_THREADS)
    TBLIS.set_num_threads(TBLIS_THREADS)
    Et = sum(Evals)
    output("Finished in {:5.5f} s", t)
    output("Final (T) contribution: {:15.10f}", Et)
    output("CCSD(T) energy:         {:15.10f}", Et+ccsd.energy)
    return RCCSDpT{T}(ccsd, Et+ccsd.energy, Et)
end


text:
using Documenter, DiDySim

makedocs(;
    modules=[DiDySim],
    format=:html,
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/dpriego87/DiDySim.jl/blob/{commit}{path}#L{line}",
    sitename="DiDySim.jl",
    authors="Daniel Priego-Espinosa",
    assets=[],
)

deploydocs(;
    repo="github.com/dpriego87/DiDySim.jl",
    target="build",
    julia="0.6",
    deps=nothing,
    make=nothing,
)


text:
using Sandboxes: baremod, sandboxify, sandboxed_eval_expr
using Test

"Test `ex` evaluates without throwing"
macro test_nothrow(ex)
    :(@test ($(esc(ex)); true))
end

@testset "Sandboxes.jl" begin
    s = baremod()

    seval(expr) = eval(sandboxed_eval_expr(s, expr))

    # We can inject functionality into the sandbox
    @test seval(:( nand(a, b) = $(~)($(&)(a, b)) ) ) == s.nand
    @test seval(:( not(a) = nand(a, true) ) ) == s.not
    @test seval(:( not(false) ) ) == true == s.not(false)

    # We can access properties and create symbols with names that are exported from Core
    @test_nothrow sandboxify( :( X.Int ) )
    @test_nothrow sandboxify( :( :Int ) )
    @test_nothrow sandboxify( :( Int = "foo" ) )

    # Some expression types are denied
    @test_throws Exception sandboxify( :( import X ) )
    @test_throws Exception sandboxify( :( import X: y ) )
    @test_throws Exception sandboxify( :( using X ) )
    @test_throws Exception sandboxify( :( using X: y ) )
    @test_throws Exception sandboxify( :( baremodule X end ) )
    @test_throws Exception sandboxify( :( baremodule X end ) )
    @test_throws Exception sandboxify( :( module X end; X.eval ) )

    # Access to values of `Core` is denied
    @test_throws UndefVarError seval(:( Core ) )
    @test s.eval == Core.eval
    @test_throws UndefVarError seval(:( eval ) )

    # Can't use a macro to escape the sandbox
    @test_throws Exception sandboxify(:(macro x() :($(:Core)) end))

    # Eventually it would be nice to allow macro definitions within the
    # sandbox, so here's a broken test to remind us.
    # seval(:(macro x() :($(:Core)) end)))
    @test_broken seval(:(@x().eval)) != Core.eval

    # DataTypes let you get a reference to their defining module, which we
    # should probably block :)
    @test_broken seval(:(f(x::T) = T; f(1).name.module)) != Core

    # Until you define something with that name
    @test seval(:( eval() = "hello" ) ) != s.eval
    @test seval(:( eval() ) ) == "hello"
end

max_stars_repo_name: cmcaine/Sandboxes

max_issues_repo_name: cmcaine/Sandboxes

max_forks_repo_name: cmcaine/Sandboxes


text:
module Interventions
include("intervene.jl")
include("higherintervene.jl")
include("intervenepass.jl")
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using Test
using LinearAlgebra
using LinearAlgebra: norm
using Ferrite
using StaticArrays

using TopOpt
using TrussTopOpt.TrussTopOptProblems
using TrussTopOpt.TrussTopOptProblems: buckling
using Crayons.Box
using GLMakie

problem_json = ["mgz_geom_stiff_ex9.1.json", "buckling_2d_debug.json"]
ins_dir = joinpath(@__DIR__, "instances", "fea_examples");

# @testset "Buckling problem solve - $(problem_json[i])" for i in 1:length(problem_json)
    i = 1
    file_name = problem_json[i]
    problem_file = joinpath(ins_dir, file_name)

    node_points, elements, mats, crosssecs, fixities, _ = parse_truss_json(problem_file);

    P_val = 400 # kN
    P = Dict(2 => SVector{2}(P_val * [0.05, -1.0]))
    nnode = length(node_points)

    # Newton iteration setup
    max_newton_itr = 50
    NEWTON_TOL = 1e-8

    function dict2vec(dict, veclen)
        v = zeros(veclen)
        for (nid, val) in dict
            v[2*nid-1:2*nid] = val
        end
        return v
    end

    # Δu = Dict(nid => SVector{2,Float64}(0,0) for nid in keys(node_points))

    updated_nodes = copy(node_points)
    problem = TrussProblem(Val{:Linear}, node_points, elements, P, fixities, mats, crosssecs);
    solver = FEASolver(Displacement, Direct, problem)
    solver()
    K, Kσ = buckling(problem, solver.globalinfo, solver.elementinfo; u=zeros(2*nnode))
    Kt = K + Kσ
    P_vec = solver.globalinfo.f
    # * initialize displacements as zero
    ut = zeros(nnode*2)

    # tracking iteration history for plotting
    ub_t = []
    p_t = []

    newton_itr = -1

    # Newton-Raphson Iterations
    # https://people.duke.edu/~hpgavin/cee421/truss-finite-def.pdf
    # https://kristofferc.github.io/JuAFEM.jl/dev/examples/hyperelasticity/
    while true
        global newton_itr, Kt, P_vec, ut, updated_nodes
        newton_itr += 1
        reaction = Kt * ut
        @show res_norm = norm(reaction - P_vec)
        push!(ub_t, ut[2*2-1])
        push!(p_t,-reaction[2*2])

        if res_norm < NEWTON_TOL
            break
        elseif newton_itr > max_newton_itr
            break
            # error("Reached maximum Newton iterations, aborting")
        end

        # construct new stiffness matrix
        ut = Kt \ P_vec
        updated_nodes = Dict(nid => pt + ut[2*nid-1:2*nid] for (nid, pt) in updated_nodes)
        problem = TrussProblem(Val{:Linear}, updated_nodes, elements, P, fixities, mats, crosssecs);
        solver = FEASolver(Displacement, Direct, problem)
        solver()
        K, Kσ = buckling(problem, solver.globalinfo, solver.elementinfo; u=ut)
        Kt = K + Kσ
    end
    scatter(ub_t, p_t)
    lines!(ub_t, p_t, color = :blue, linewidth = 3)
    current_figure()

    # scene, layout = draw_truss_problem(problem; u=v, default_exagg=1.0, exagg_range=10.0,
    #     default_load_scale=0.2, default_support_scale=0.2, default_arrow_size=0.03)

# end # end test set

# =#


text:
@testset "Variable type" begin
    d = Degree(2,3)

    z = MPC.Variable("z", Complex)
    x = MPC.Variable("x", Real)
    b = MPC.Variable("b", Bool)

    @test iscomplex(z) == true
    @test isreal(x) == true
    @test isbool(b) == true

    @test hash(Degree(1,3)) == hash(Degree(1,3))
    @test hash(MPC.Variable("a",Complex)) == hash(MPC.Variable("a",Complex))
end

@testset "Exponent type" begin
    z = MPC.Variable("z", Complex)
    x = MPC.Variable("x", Real)

    one = Exponent()
    zexpo = Exponent(z)
    expo = Exponent(SortedDict(z=>Degree(3,2),
                               x=>Degree(2,0)))

    @test zexpo == Exponent(SortedDict(z=>Degree(1,0)))
    @test Exponent(SortedDict(z=>Degree(3,2),
                              x=>Degree(0,0))) == Exponent(SortedDict(z=>Degree(3,2)))
end

@testset "Polynomial type" begin
    emptypoly = Polynomial()

    a = MPC.Variable("a", Complex)
    b = MPC.Variable("b", Real)

    one = Exponent()
    expoa2 = Exponent(SortedDict(a=>Degree(1,1)))
    expob2 = Exponent(SortedDict(b=>Degree(1,0)))
    expoa = Exponent(SortedDict(a=>Degree(1,0)))
    expob = Exponent(SortedDict(b=>Degree(1,0)))

    p_obj = Polynomial(SortedDict{Exponent, Number}(expoa2=>1.0,
                                                    expob2=>3.5,
                                                    expoa=>5+3im,
                                                    expob=>0,
                                                    one=>-1))

    io = IOBuffer()
    print(io, p_obj)
    @test String(take!(io)) == "-1 + (5 + 3im)*a + conj(a) * a"
end

@testset "Point type" begin

    pt = Point()

    z = MPC.Variable("z", Complex)
    x1 = MPC.Variable("x1", Real)
    x2 = MPC.Variable("x2", Real)
    b = MPC.Variable("b", Bool)
    b1 = MPC.Variable("b1", Bool)

    pointdict = SortedDict(z =>1+2im,
                          x1 => 3.5,
                          x2 => 0,
                          b => 1,
                          b1 => 3.5)

    pt = Point(pointdict)
    @test pt == Point([z x1 x2 b b1], [1+2im 3.5 0 1 3.5])
    @test length(pt) == 4

    io = IOBuffer()
    print(io, pt)
    @test String(take!(io)) == "b  1\nb1 1\nx1 3.5\nz  1 + 2im\n"
end

@testset "Constraint type" begin
    a = MPC.Variable("a", Complex)
    b = MPC.Variable("b", Real)

    one = Exponent()
    expoa2 = Exponent(SortedDict(a=>Degree(1,1)))
    expob2 = Exponent(SortedDict(b=>Degree(1,0)))
    expoa = Exponent(SortedDict(a=>Degree(1,0)))
    expob = Exponent(SortedDict(b=>Degree(1,0)))

    p1 = Polynomial(SortedDict{Exponent, Number}(expoa2=>1.0,
                                                 expob2=>3.5,
                                                 expoa=>5+3im,
                                                 expob=>0,
                                                 one=>-1))

    ub_c = 3+5im
    lb_c = -3-2im

    C1 = p1 << ub_c
    C2 = ub_c >> p1
    @test C1 == Constraint(p1, -Inf-Inf*im, ub_c)
    @test C2 == Constraint(p1, -Inf-Inf*im, ub_c)

    @test (lb_c << C1) == Constraint(p1, lb_c, ub_c)
    @test (C2 >> lb_c) == Constraint(p1, lb_c, ub_c)

    C1 = lb_c << p1
    C2 = p1 >> lb_c
    @test C1 == Constraint(p1, lb_c, +Inf+Inf*im)
    @test C2 == Constraint(p1, lb_c, +Inf+Inf*im)

    @test (C1 << ub_c) == Constraint(p1, lb_c, ub_c)
    @test (ub_c >> C2) == Constraint(p1, lb_c, ub_c)

    @test (p1 == ub_c) == Constraint(p1, ub_c, ub_c)
    @test (p1 == ub_c) != (p1 << ub_c)

    io = IOBuffer()
    print(io, Constraint(p1, lb_c, ub_c))
    @test String(take!(io)) == "-3 - 2im < -1 + (5 + 3im)*a + conj(a) * a < 3 + 5im"

    print(io, p1 << ub_c)
    @test String(take!(io)) == "-1 + (5 + 3im)*a + conj(a) * a < 3 + 5im"
    print(io, lb_c << p1)
    @test String(take!(io)) == "-3 - 2im < -1 + (5 + 3im)*a + conj(a) * a"

    print(io, Constraint(p1, ub_c, ub_c))
    @test String(take!(io)) == "-1 + (5 + 3im)*a + conj(a) * a = 3 + 5im"
end

@testset "Problem type" begin
    a = MPC.Variable("a", Complex)
    b = MPC.Variable("b", Real)

    one = Exponent()
    expoa2 = Exponent(SortedDict(a=>Degree(1,1)))
    expob2 = Exponent(SortedDict(b=>Degree(1,0)))
    expoa = Exponent(SortedDict(a=>Degree(1,0)))
    expob = Exponent(SortedDict(b=>Degree(1,0)))

    p1 = Polynomial(SortedDict{Exponent, Number}(expoa2=>1.0,
                                                expob2=>3.5,
                                                expoa=>5+3im,
                                                expob=>0,
                                                one=>-1))

    vars = SortedDict{String, Type}("a" => Complex, "b" => Real)
    ctrs = SortedDict{String, Constraint}("ctr1" => ((-3+5im) << p1),
                                          "ctr2" => (p1 == 2.5))
    pb = Problem(p1, ctrs, vars)
    Problem()

    io = IOBuffer()
    print(io, pb)
    @test String(take!(io)) == "▶ variables: a b \n▶ objective: -1 + (5 + 3im)*a + conj(a) * a\n▶ constraints: \n →       ctr1: -3 + 5im < -1 + (5 + 3im)*a + conj(a) * a\n →       ctr2: -1 + (5 + 3im)*a + conj(a) * a = 2.5\n"

end


text:
using Microscopy
using Documenter

DocMeta.setdocmeta!(Microscopy, :DocTestSetup, :(using Microscopy); recursive=true)

makedocs(;
    modules=[Microscopy],
    authors="roflmaostc <fxw+git@mailbox.org> and contributors",
    repo="https://github.com/JuliaMicroscopy/Microscopy.jl/blob/{commit}{path}#{line}",
    sitename="Microscopy.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaMicroscopy.github.io/Microscopy.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaMicroscopy/Microscopy.jl",
    devbranch = "main"
)

max_stars_repo_name: JuliaMicroscopy/Microscopy

max_forks_repo_name: JuliaMicroscopy/Microscopy


text:
using Compat.SparseArrays

linear_input_type(x::T, nbatch) where T = linear_input_type(T, nbatch)

linear_input_type(::Type{Matrix{T}}, nbatch) where T = nbatch > 1 ? Matrix{T} : Vector{T}
linear_input_type(::Type{SparseMatrixCSC{Tv, Ti}}, nbatch) where {Tv, Ti} =
    nbatch > 1 ? SparseMatrixCSC{Tv, Ti} : SparseVector{Tv, Ti}

function linear_input_type(::Type{SArray{Tuple{out, in}, T, 2, L}}, nbatch) where {out, in, T, L}
    if nbatch > 1
        return SArray{Tuple{in, nbatch}, T, 2, in}
    else
        return SArray{Tuple{in}, T, 1, in}
    end
end

function linear_input_type(::Type{MArray{Tuple{out, in}, T, 2, L}}, nbatch) where {out, in, T, L}
    if nbatch > 1
        return MArray{Tuple{in, nbatch}, T, 2, in}
    else
        return MArray{Tuple{in}, T, 1, in}
    end
end

# TODO: GPU array
function linear_input_type end


text:
module Jib

using DataFrames

include("client.jl")
include("enums.jl")
include("types.jl")
include("types_condition.jl")
include("types_mutable.jl")
include("types_private.jl")
include("requests.jl")        ; using .Requests
include("wrapper.jl")
include("reader.jl")          ; using .Reader: check_all, start_reader
include("utils.jl")
include("connect.jl")         ; using .Connect: connect, disconnect

end


text:
import Plots: plot, xlabel!, ylabel!, title!, plot!
using MCMCChains

tpsA549 = @SVector Float64[60, 240]
gasA549 = @SVector Float64[64, 16, 4, 1, 0.25, 0]
pYA549 = @SMatrix [10.8 8.3; 7.4 7.1; 7.1 7.7; 4.6 8.2; 6.1 7.2; 7.5 7.5]
totA549 = @SMatrix [3443.1 3219.7; 3143.4 3353.8; 3018.9 3611.8; 2608.9 3448.2; 2690.2 3168.1; 2672.0 2672.0]
surfA549 = @SMatrix [0.206 0.239; 0.274 0.316; 0.281 0.251; 0.220 0.302; 0.256 0.281; 0.257 0.337]

" Run the model calculations for the A549 measurements. "
function dataModelCalc(tps, g6conc, params, scale, scaleSurf)
    pYresids = Matrix{typeof(scale)}(undef, length(tps), length(g6conc))
    totalresids = Matrix{typeof(scale)}(undef, length(tps), length(g6conc))
    surfresids = Matrix{typeof(scale)}(undef, length(tps), length(g6conc))

    # Setup reductions
    pYAXL = pY .* recpSpecific[1]
    surfAXL = surface .* recpSpecific[1]
    totAXL = total .* recpSpecific[1]

    paramsStart = param(params)
    solInit = getAutocrine(paramsStart)

    Threads.@threads for ii = 1:length(g6conc)
        params = deepcopy(paramsStart)

        params.gasCur += g6conc[ii]
        data = runTAMinit(tps, params, solInit)

        pYresids[:, ii] = (data * pYAXL) * scale
        surfresids[:, ii] = (data * surfAXL) * scaleSurf
        totalresids[:, ii] = (data * totAXL)
    end

    return pYresids, totalresids, surfresids
end

@model AXLfit(pYDataExp, surfDataExp, totDataExp, sqResid, tps, g6conc, ::Type{TV} = Vector{Float64}) where {TV} = begin
    internalize ~ LogNormal(log(0.1), 0.1)
    pYinternalize ~ LogNormal(log(1.0), 0.1)
    sortF ~ Beta(2.0, 20.0)
    kRec ~ LogNormal(log(0.1), 0.1)
    kDeg ~ LogNormal(log(0.01), 0.1)
    xFwd ~ LogNormal(-3.0, 1.0)
    gasCur ~ LogNormal(log(0.1), 0.1)
    AXLexpr ~ LogNormal(log(1.0), 1.0)
    Ig2rev ~ LogNormal(0.0, 1.0)
    scale ~ LogNormal(0.0, 1.0)
    scaleSurf ~ LogNormal(0.0, 1.0)

    params = vcat(internalize, pYinternalize, sortF, kRec, kDeg, xFwd, gasCur, AXLexpr, zeros(2), Ig2rev, ones(4))
    pYresids, totalresids, surfresids = dataModelCalc(tps, g6conc, params, scale, scaleSurf)

    # The scaling is based on the average stderr, to make these proportional to the std MvNorm
    # TODO: Identify exact values
    sqResid = sqeuclidean(pYresids, transpose(pYDataExp))
    sqResid += sqeuclidean(totalresids, transpose(totDataExp)) / 100.0
    sqResid += sqeuclidean(surfresids, transpose(surfDataExp)) / 0.038

    sqResid ~ Chisq(length(pYDataExp) + length(totDataExp) + length(surfDataExp))
end

function plot_overlay(chn)
    Ig2rev = get(chn, :Ig2rev)[1]
    scale = get(chn, :scale)[1]
    scaleSurf = get(chn, :scaleSurf)[1]

    x = get(chn, [:internalize, :pYinternalize, :sortF, :kRec, :kDeg, :xFwd, :gasCur, :AXLexpr])
    samp_params = hcat(x.internalize, x.pYinternalize, x.sortF, x.kRec, x.kDeg, x.xFwd, x.gasCur, x.AXLexpr)

    pY = Array{Float64}(undef, size(samp_params, 1), length(tpsA549), length(gasA549))
    tot = Array{Float64}(undef, size(samp_params, 1), length(tpsA549), length(gasA549))
    surf = Array{Float64}(undef, size(samp_params, 1), length(tpsA549), length(gasA549))

    for iter = 1:size(samp_params, 1)
        params = vcat(samp_params[iter, :], zeros(2), Ig2rev[iter], ones(4))
        pY[iter, :, :], tot[iter, :, :], surf[iter, :, :] = dataModelCalc(tpsA549, gasA549, params, scale[iter], scaleSurf[iter])
    end

    medpY = Statistics.median(pY, dims = 1)
    medtot = Statistics.median(tot, dims = 1)
    medsurf = Statistics.median(surf, dims = 1)

    tp1_calcmed = hcat(transpose(medpY[:, 1, :]), transpose(medsurf[:, 1, :]), transpose(medtot[:, 1, :]))
    tp2_calcmed = hcat(transpose(medpY[:, 2, :]), transpose(medsurf[:, 2, :]), transpose(medtot[:, 2, :]))
    tp1_exp = hcat(pYA549[:, 1], surfA549[:, 1], totA549[:, 1])
    tp2_exp = hcat(pYA549[:, 2], surfA549[:, 2], totA549[:, 2])

    plot(
        gasA549,
        [tp1_calcmed, tp2_calcmed],
        label = ["1 hr, calc" "4 hr, calc"],
        title = ["Phosphorylated receptor" "Surface receptor" "Total receptor"],
        lw = 3,
        layout = (1, 3),
        size = (1200, 400),
    )
    plot!(gasA549, [tp1_exp, tp2_exp], label = ["1 hr, exp" "4 hr, exp"], lw = 3, layout = (1, 3))
    xlabel!("Gas6 Concentration (nM)")
end

A549model = AXLfit(pYA549, surfA549, totA549, [0.0], tpsA549, gasA549)


text:
using Test

using Algorithm101
const ALGO = Algorithm101

using DataStructure101
const DS = DataStructure101

@testset "Tree Traverse" begin

    bt = DS.BinaryTree()
    @test bt isa DS.BinaryTree

    for i=1:10
        bt[i] = i
    end
    @test ALGO.pre_order(DS.root(bt)) == [1,2,4,8,9,5,10,3,6,7]
    @test ALGO.in_order(DS.root(bt)) == [8,4,9,2,10,5,1,6,3,7]
    @test ALGO.pos_order(DS.root(bt)) == [8,9,4,10,5,2,6,7,3,1]

end


text:
function assess_examinee_mmle!(
    examinee::AbstractExaminee,
    items::Vector{<:AbstractItem},
    responses::Vector{Response};
    clean = false,
    rescale_latent = false,
    dist::Distributions.DiscreteUnivariateDistribution = Distributions.DiscreteNonParametric([0.0,1.0], [0.5, 0.5]),
    max_time::Int64 = 100,
    max_iter::Int64 = 10,
    x_tol_rel::Float64 = 0.001,
    f_tol_rel::Float64 = 0.00001,
    kwargs...
    )
    if !examinee.latent.assessed
        if !clean
            responses = sort(filter(r -> r.examinee_idx == examinee.idx, responses), by = r -> r.item_idx)
            items_idx = map(r -> r.item_idx, responses)
            items = sort(filter( i -> i.idx in items_idx, items), by = i -> i.idx)
        end

        parameters = get_parameters(items)
        latent = examinee.latent
        latent.prior = dist

        I = size(parameters, 1)
        responses_per_examinee = map( r -> r.val ,responses)
        iter = 1

        # starting values
        before_time = time()
        stop = false
        old_likelihood = 0.0
        old_latents = hcat(_get_latents.([latent])...)
        start_time = time()

        latent = _update_posterior(latent, parameters, responses_per_examinee),
        while !stop
            #rescale dist
            if rescale_latent
                dist = _rescale(
                    dist,
                    [latent];
                    metric = metric
                )
                latent.prior = dist
            end
            latent = _update_posterior(latent, parameters[i_index], responses_per_examinee),

            if any([
                check_iter(iter; max_iter = max_iter),
                check_time(start_time; max_time = max_time),
                check_f_tol_rel!(
                    [latent],
                    old_likelihood;
                    f_tol_rel = f_tol_rel
                    ),
                check_x_tol_rel!(
                    [latent],
                    old_latents;
                    x_tol_rel = x_tol_rel
                )]
                )
                stop = true
            end
            iter += 1
        end
        latent.posterior = Distributions.DiscreteNonParametric(dist.support, latent.posterior.p)
        latent.val = latent.posterior.p' * latent.posterior.support
        examinee = Examinee(examinee.idx, examinee.id, latent)
    end
    return nothing
end

function assess_examinee_mmle!(
    examinees::Vector{<:AbstractExaminee},
    items::Vector{<:AbstractItem},
    responses::Vector{Response};
    dist::Distributions.DiscreteUnivariateDistribution = Distributions.DiscreteNonParametric([0.0,1.0], [0.5, 0.5]),
    clean = false,
    max_time::Int64 = 100,
    max_iter::Int64 = 10,
    x_tol_rel::Float64 = 0.001,
    f_tol_rel::Float64 = 0.00001,
    kwargs...
    )
    pmap( e -> assess_examinee_mmle!(
            e,
            items,
            responses;
            clean = false,
            max_time = max_time,
            max_iter = max_iter,
            x_tol_rel = x_tol_rel,
            f_tol_rel = f_tol_rel
            ),
            examinees,
            batch_size = Int(ceil(size(examinees, 1) / (nprocs()-1)))
            )
    return nothing
end


text:
immutable BilinearSpline <: AbstractInterpolation
    BilinearSpline() = BilinearSpline
end

immutable BilinearSplineInterpolation{T}
    x::Vector{T}
    y::Vector{T}
    z::Array{T,2}
end

function default_settings{D,T}(::Type{BilinearSpline},axes::NTuple{D,T})
    nx = length(axes[1])
    ny = length(axes[2])
    return Irregular{2,Val{(nx,ny)}}, Nearest
end

function BilinearSpline(x::AbstractVector, y::AbstractVector, z::AbstractArray)
  issorted(x) || throw(ArgumentError("x points must be in ascending order"))
  issorted(y) || throw(ArgumentError("y points must be in ascending order"))
  x = convert(Array{Float64,1},x)
  y = convert(Array{Float64,1},y)
  z = convert(Array{Float64,2},z)

  return BilinearSplineInterpolation{eltype(x)}(x,y,z)
end

@inline function preprocess(::Type{BilinearSpline}, axes, f)
    return BilinearSpline(axes...,f)
end

function apply_bc(::Type{BilinearSpline},BC::Type{Nearest}, I::Symbol, x::Symbol)
    xbc = apply_nearest(I,1, x, 1)
    ybc = apply_nearest(I,2, x, 2)
    quote
        $(xbc)
        $(ybc)
    end
end

function apply_bc(::Type{BilinearSpline},BC::Type{Reflect}, I::Symbol, x::Symbol)
    xbc = apply_reflect(I,1,x,1)
    ybc = apply_reflect(I,2,x,2)
    quote
        $(xbc)
        $(ybc)
    end
end

function apply_bc{V}(::Type{BilinearSpline},BC::Type{Value{V}}, I::Symbol, x::Symbol)
    xbc = apply_value(BC,I,1,x,1)
    ybc = apply_value(BC,I,2,x,2)
    quote
        $(xbc)
        $(ybc)
    end
end

function apply_bc(::Type{BilinearSpline},BC::Type{None}, I::Symbol, x::Symbol)
    apply_none()
end

function apply_bc(::Type{BilinearSpline},BC::Type{Error}, I::Symbol, x::Symbol)
    xbc = apply_error(I,1,x,1)
    ybc = apply_error(I,2,x,2)
    quote
        $(xbc)
        $(ybc)
    end
end

@inline function get_indices{S}(::Type{BilinearSpline}, G::Type{Uniform{2,S}}, I, x)
    dx = I.x[2]-I.x[1]
    ix = clamp(Int(floor((x[1] - I.x[1])/dx)) + 1, 1, S.parameters[1][1]-1)
    dy = I.y[2]-I.y[1]
    iy = clamp(Int(floor((x[2] - I.y[1])/dy)) + 1, 1, S.parameters[1][2]-1)
    return (ix, iy)
end

@inline function get_indices{S}(::Type{BilinearSpline}, G::Type{Irregular{2,S}}, I, x)
    xr = searchsortedlast(I.x,x[1],Base.Order.Forward)
    ix = clamp(xr, 1, S.parameters[1][1]-1)
    yr = searchsortedlast(I.y,x[2],Base.Order.Forward)
    iy = clamp(yr, 1, S.parameters[1][2]-1)
    return (ix, iy)
end

@inline function get_value(::Type{BilinearSpline}, I, x, inds)
    ix, iy = inds
    xp = x[1]
    yp = x[2]
    x1, x2 = I.x[ix], I.x[ix+1]
    y1, y2 = I.y[iy], I.y[iy+1]
    dxdy = (x2-x1)*(y2-y1)

    b11 = (x2 - xp)*(y2 - yp)
    b21 = (xp - x1)*(y2 - yp)
    b12 = (x2 - xp)*(yp - y1)
    b22 = (xp - x1)*(yp - y1)
    v = b11*I.z[ix,iy] + b12*I.z[ix,iy+1] + b21*I.z[ix+1,iy] + b22*I.z[ix+1,iy+1]
    return v/dxdy
end


text:
function Plane(arg0::Plane)
    return Plane((Plane,), arg0)
end

function Plane(arg0::Vector3D, arg1::Vector3D, arg2::Vector3D, arg3::jdouble)
    return Plane((Vector3D, Vector3D, Vector3D, jdouble), arg0, arg1, arg2, arg3)
end

function Plane(arg0::Vector3D, arg1::Vector3D, arg2::jdouble)
    return Plane((Vector3D, Vector3D, jdouble), arg0, arg1, arg2)
end

function Plane(arg0::Vector3D, arg1::jdouble)
    return Plane((Vector3D, jdouble), arg0, arg1)
end

function contains(obj::Plane, arg0::Vector3D)
    return jcall(obj, "contains", jboolean, (Vector3D,), arg0)
end

function copy_self(obj::Plane)
    return jcall(obj, "copySelf", Plane, ())
end

function empty_hyperplane(obj::Plane)
    return jcall(obj, "emptyHyperplane", SubPlane, ())
end

function get_normal(obj::Plane)
    return jcall(obj, "getNormal", Vector3D, ())
end

function get_offset(obj::Plane, arg0::JVector)
    return jcall(obj, "getOffset", jdouble, (JVector,), arg0)
end

function get_offset(obj::Plane, arg0::Plane)
    return jcall(obj, "getOffset", jdouble, (Plane,), arg0)
end

function get_offset(obj::Plane, arg0::Point)
    return jcall(obj, "getOffset", jdouble, (Point,), arg0)
end

function get_origin(obj::Plane)
    return jcall(obj, "getOrigin", Vector3D, ())
end

function get_point_at(obj::Plane, arg0::Vector2D, arg1::jdouble)
    return jcall(obj, "getPointAt", Vector3D, (Vector2D, jdouble), arg0, arg1)
end

function get_tolerance(obj::Plane)
    return jcall(obj, "getTolerance", jdouble, ())
end

function get_u(obj::Plane)
    return jcall(obj, "getU", Vector3D, ())
end

function get_v(obj::Plane)
    return jcall(obj, "getV", Vector3D, ())
end

function intersection(::Type{Plane}, arg0::Plane, arg1::Plane, arg2::Plane)
    return jcall(Plane, "intersection", Vector3D, (Plane, Plane, Plane), arg0, arg1, arg2)
end

function intersection(obj::Plane, arg0::Line)
    return jcall(obj, "intersection", Vector3D, (Line,), arg0)
end

function intersection(obj::Plane, arg0::Plane)
    return jcall(obj, "intersection", Line, (Plane,), arg0)
end

function is_similar_to(obj::Plane, arg0::Plane)
    return jcall(obj, "isSimilarTo", jboolean, (Plane,), arg0)
end

function project(obj::Plane, arg0::Point)
    return jcall(obj, "project", Point, (Point,), arg0)
end

function reset(obj::Plane, arg0::Plane)
    return jcall(obj, "reset", void, (Plane,), arg0)
end

function reset(obj::Plane, arg0::Vector3D, arg1::Vector3D)
    return jcall(obj, "reset", void, (Vector3D, Vector3D), arg0, arg1)
end

function revert_self(obj::Plane)
    return jcall(obj, "revertSelf", void, ())
end

function rotate(obj::Plane, arg0::Vector3D, arg1::Rotation)
    return jcall(obj, "rotate", Plane, (Vector3D, Rotation), arg0, arg1)
end

function same_orientation_as(obj::Plane, arg0::Hyperplane)
    return jcall(obj, "sameOrientationAs", jboolean, (Hyperplane,), arg0)
end

function to_space(obj::Plane, arg0::JVector)
    return jcall(obj, "toSpace", Vector3D, (JVector,), arg0)
end

function to_space(obj::Plane, arg0::Point)
    return jcall(obj, "toSpace", Vector3D, (Point,), arg0)
end

function to_sub_space(obj::Plane, arg0::JVector)
    return jcall(obj, "toSubSpace", Vector2D, (JVector,), arg0)
end

function to_sub_space(obj::Plane, arg0::Point)
    return jcall(obj, "toSubSpace", Vector2D, (Point,), arg0)
end

function translate(obj::Plane, arg0::Vector3D)
    return jcall(obj, "translate", Plane, (Vector3D,), arg0)
end

function whole_hyperplane(obj::Plane)
    return jcall(obj, "wholeHyperplane", SubPlane, ())
end

function whole_space(obj::Plane)
    return jcall(obj, "wholeSpace", PolyhedronsSet, ())
end


text:
include("header.jl")
# @time include("core.jl")
@time include("primitive.jl")
@time include("nonlinear.jl")
@time include("cnn.jl")
@time include("rnn.jl")
@time include("loss.jl")
@time include("special.jl")
@time include("data.jl")


text:
################################################################################
# PWSet.jl
#
# The Persistent weighted collection PWSet which is built on top of the
# persistent heap type in PHeap.jl transient heap type.
#
# @author Noah C. Benson
#
# MIT License
# Copyright (c) 2020 Noah C. Benson

# ==============================================================================
# AbstractPWSet
"""
    AbstractPWSet{T,W}

AbstractPWSet is an abstract type implemented by all persistent weighted set
types.

See also: [`PWSet`](@ref), [`AbstractPSet`](@ref).
"""
abstract type AbstractPWSet{T,W<:Number} <: AbstractPSet{T} end

# ==============================================================================
# PWSet
macro _pwset_code(name::Symbol, dicttype::Symbol)
    return quote
        struct $name{T,W} <: AbstractPWSet{T,W}
            heap::PHeap{T,W,typeof(>),$dicttype{T,Int}}
        end
        function $name{T,W}(itr::AbstractArray) where {T,W<:Number}
            return $name{T,W}(itr...)
        end
        function $name{T,W}(itr::AbstractSet) where {T,W<:Number}
            return $name{T,W}(itr...)
        end
        function $name{T,W}(itr::AbstractDict) where {T,W<:Number}
            return $name{T,W}(itr...)
        end
        function $name{T,W}() where {T,W<:Number}
            return $name{T,W}(PHeap{T,W,typeof(>),$dicttype{T,Int}}(>))
        end
        function $name{T,W}(d::$name{T,W}) where {T,W}
            return d
        end
        function $name{T,W}(ps::Union{Tuple,Pair}...) where {T,W<:Number}
            return reduce(push, ps, init=$name{T,W}())
        end
        function $name{T}(itr) where {T}
            return $name{T}(itr...)
        end
        function $name{T}(::Tuple{}) where {T,W<:Number}
            return $name{T,Float64}(PHeap{T,Float64,typeof(>),$dicttype{T,Int}}(>))
        end
        function $name{T}() where {T,W<:Number}
            return $name{T,Float64}(PHeap{T,Float64,typeof(>),$dicttype{T,Int}}(>))
        end
        function $name{T}(d::$name{T,W}) where {T,W}
            return d
        end
        function $name{T}(ps::Union{Tuple,Pair}...) where {T}
            return $name{T,Float64}(ps...)
        end
        function $name(itr)
            return $name(itr...)
        end
        function $name(d::$name{T,W}) where {T,W}
            return d
        end
        function $name(ps::Union{Tuple,Pair}...)
            return $name{Any,Float64}(ps...)
        end
        function $name(::Tuple{})
            return $name{Any,Float64}()
        end
        function $name()
            return $name{Any,Float64}(PHeap{Any,Float64,typeof(>),$dicttype{T,Int}}(>))
        end
        # Document the equal/hash types.
        equalfn(::Type{T}) where {T <: $name} = equalfn($(dicttype))
        hashfn(::Type{T}) where {T <: $name} = hashfn($(dicttype))
        equalfn(::$name) = equalfn($(dicttype))
        hashfn(::$name) = hashfn($(dicttype))
        # Generic construction.
        # Base methods.
        Base.length(s::$name) = length(s.heap)
        Base.iterate(u::$name{T,W}) where {T,W} = iterate(u.heap)
        Base.iterate(u::$name{T,W}, x) where {T,W} = iterate(u.heap, x)
        Base.in(s::S, u::$name{T,W}, eqfn::Function) where {S,T,W} = in(s, u.heap, eqfn)
        Base.in(s::S, u::$name{T,W}) where {S,T,W} = in(s, u.heap)
        # And the Air methods.
        Air.push(u::$name{T,W}, x::Tuple{S,X}) where {T,W,S,X<:Number} = begin
            heap = push(u.heap, x)
            (heap === u.heap) && return u
            return $name{T,W}(heap)
        end
        Air.push(u::$name{T,W}, x::Pair{S,X}) where {T,W,S,X<:Number} = push(u, (x[1],x[2]))
        Air.pop(u::$name{T,W}) where {T,W} = begin
            (length(u) == 0) && throw(ArgumentError("$name must be non-empty"))
            return $name{T,W}(pop(u.heap))
        end
        Base.first(u::$name{T,W}) where {T,W} = begin
            (length(u) == 0) && throw(ArgumentError("$name must be non-empty"))
            return first(u.heap)
        end
        Air.delete(u::$name{T,W}, s::S) where {T,W,S} = begin
            heap = delete(u.heap, s)
            (heap === u.heap) && return u
            return $name{T,W}(heap)
        end
        Air.getweight(u::$name{T,W}, x::J) where {T,W,J} = getweight(u.heap, x)
        Air.setweight(u::$name{T,W}, x::J, w::X) where {T,W,J,X<:Number} = begin
            heap = setweight(u.heap, x, w)
            (heap === u.heap) && return u
            return $name{T,W}(heap)
        end
        Random.rand(u::$name{T,W}) where {T,W} = begin
            (length(u) > 0) || throw(ArgumentError("PWSet must be non-empty"))
            return rand(u.heap)
        end
    end |> esc
end
# Declare the types.
@_pwset_code PWSet PDict
@_pwset_code PWIdSet PIdDict
# Document the types.
@doc """
    PWSet{K,V}

A persistent set with weighted elements. As such, a `PWSet` supports the
typical operations of a `PSet` as well as the following:
* The `first` function yields the element with the highest weight.
* The `pop` function yields a copy of the set without the element that has
  the highest weight.
* Iteration occurs in the order of greatest to least weight.
* The weights can be changed with the `getweight` and `setweight` functions;
  `setweight` yields a duplicate dictionary with updated weights.

See also: [`PSet`](@ref), [`getweight`](@ref), [`setweight`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest; filter=r"PWSet{Symbol, ?Float64} with 3 elements:"
julia> PWSet{Symbol}(:a => 0.1, :b => 0.2, :c => 0.3)
PWSet{Symbol,Float64} with 3 elements:
  :c
  :b
  :a
```
""" PWSet
@doc """
    PWIdSet{T}

A persistent set with weighted elements. As such, a `PWIdSet` supports the
typical operations of a `PIdSet` as well as the following:
* The `first` function yields the element with the highest weight.
* The `pop` function yields a copy of the set without the element that has
  the highest weight.
* Iteration occurs in the order of greatest to least weight.
* The weights can be changed with the `getweight` and `setweight` functions;
  `setweight` yields a duplicate dictionary with updated weights.

See also: [`PIdSet`](@ref), [`PWSet`](@ref), [`getweight`](@ref),
[`setweight`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest; filter=r"PWIdSet{Symbol, ?Float64} with 3 elements:"
julia> PWIdSet{Symbol}(:a => 0.1, :b => 0.2, :c => 0.3)
PWIdSet{Symbol,Float64} with 3 elements:
  :c
  :b
  :a
```
""" PWIdSet
# A few functions needed here instead of inside the macro.
Base.empty(s::PWSet{T,W}, ::Type{S}=T, ::Type{X}=W) where {T,W,S,X} = PWSet{S,X}()
Base.empty(s::PWIdSet{T,W}, ::Type{S}=T, ::Type{X}=W) where {T,W,S,X} = PWIdSet{S,X}()
Base.isequal(s::ST, t) where {ST<:AbstractPWSet} = false
Base.isequal(t, s::ST) where {ST<:AbstractPWSet} = false
Base.isequal(s::SS, t::ST) where {
    S,T,
    SS <: AbstractPWSet{S},
    ST <: AbstractPWSet{T}} = begin
    (length(t) == length(s)) || return false
    for x in s
        (x in t) || return false
        (getweight(s, x) == getweight(t, x)) || return false
    end
    (equalfn(SS) === equalfn(ST)) && return true
    for x in t
        (x in s) || return false
    end
    return true
end
Base.hash(u::PWS) where {T,W,PWS<:AbstractPWSet{T,W}} = begin
    h = length(u)
    for s in u.heap
        w = getweight(u.heap, s)
        h += hash(s) ⊻ hash(w)
    end
    return h
end

# Export the relevant symbols.
export PWSet, PWIdSet