text:
#--------------------------------------------------------------------
# DNSS.jl
# Soham 03-2022
# Simulate the Minkowski spacetime
#--------------------------------------------------------------------
using NLsolve, Printf

function minkowski(PS::ProductSpace{S1, S2})::NTuple{2, Field{ProductSpace{S1, S2}}} where {S1, S2}
    a = Field(PS, (u,v)->2)
    η = Field(PS, (u,v)->v-u)
    return (a, η)
end

function computeUboundary(PS::ProductSpace{S1, S2})::NTuple{2, Field{S2}} where {S1, S2}
    a0 = extractUboundary(Field(PS, (u,v)->2), :incoming)
    η0 = extractUboundary(Field(PS, (u,v)->(v-u)), :incoming)
    return (a0, η0)
end

function computeVboundary(PS::ProductSpace{S1, S2})::NTuple{2, Field{S1}} where {S1, S2}
    a0 = extractVboundary(Field(PS, (u,v)->2), :incoming)
    η0 = extractVboundary(Field(PS, (u,v)->(v-u)), :incoming)
    return (a0, η0)
end

function initialguess(PS::ProductSpace{S1, S2})::NTuple{2, Field{ProductSpace{S1, S2}}} where {S1, S2}
    (a0, η0) = minkowski(PS)
    n1       = Field(PS, (u,v)-> 1e-3 * rand())
    n2       = Field(PS, (u,v)-> 1e-3 * rand())
    return (a0 + n1, η0 + n2)
end

function computePatch(PS::ProductSpace{S1, S2}, ubnd::NTuple{2, Field{S2}}, vbnd::NTuple{2, Field{S1}})::NTuple{2, Field{ProductSpace{S1, S2}}} where {S1, S2}
    B = incomingboundary(PS)
    I = identity(PS)
    DU, DV = derivative(PS)
    ∂U = combineUVboundary(ubnd, vbnd, :incoming)

    function constr(U::NTuple{2, Field{S}})::NTuple{2, Field{S}} where {S}
        (a, η) = U
        C1 = DU*(DU*η) - (2/a)*(DU*a)*(DU*η)
        C2 = DV*(DV*η) - (2/a)*(DV*a)*(DV*η)
        C  = sqrt(C1*C1 + C2*C2)
        return (C, C)
    end

    function offaxis(U::NTuple{2, Field{S}})::NTuple{2, Field{S}} where {S}
        (a, η) = U
        F1 = DU*(DV*a) - (1/a)*(DU*a)*(DV*a) + (a/η)*(DU*(DV*η))
        F2 = DU*(DV*η) + (1/η)*(DU*η)*(DV*η) + (1/4)*(1/η)*(a^2)
        return (F1, F2)
    end

    function onaxis(U::NTuple{2, Field{S}})::NTuple{2, Field{S}} where {S}
        (a, η) = U
        A1 = DU*(DV*a) - (1/a)*(DU*a)*(DV*a)
        A2 = DU*(DV*η)
        return (A1, A2)
    end

    function boundary(U::NTuple{2, Field{S}})::NTuple{2, Field{S}} where {S}
        return U - ∂U
    end

    function residual(U::NTuple{2, Field{S}})::NTuple{2, Field{S}} where {S}
        F = offaxis(U) + constr(U)
        enforceregularity!(F, onaxis(U))
        enforcebc!(F, boundary(U))
        return F
    end

    function f!(fvec::Array{T,1}, x::Array{T,1}) where {T}
        fvec[:] = reshapeFromTuple(residual(reshapeToTuple(PS, 2, x)))
    end

    r = initialguess(PS)
    s = nlsolve(f!, reshapeFromTuple(r); method=solver, autodiff=:forward, show_trace=debug >= 3, ftol=1e-8, iterations=80)
    k = reshapeToTuple(PS, 2, s.zero)

    if debug >= 2
        println(" [Before solve] norm(constraints) = ", norm(constr(r)))
        println(" Converged?                       = ", converged(s))
        println(" [After solve] norm(constraints)  = ", norm(constr(k)))
    end

    return k
end

function constraints(U::NTuple{2, Field{S}})::NTuple{1, Field{S}} where {S}
    (a, η) = U
    DU, DV = derivative(a.space)
    C1 = DU*(DU*η) - (2/a)*(DU*a)*(DU*η)
    C2 = DV*(DV*η) - (2/a)*(DV*a)*(DV*η)
    C  = sqrt(C1*C1 + C2*C2)
    return (C, )
end

function pconv(range::StepRange, nh::Int=1)
    println("Testing p convergence")
    n_ = collect(range)
    l_ = zeros(size(n_))
    for index in CartesianIndices(n_)
        np = n_[index]
        params = Parameters((np, np), (nh,nh), urange, vrange, 2)
        q = distribute(params, computePatch, computeUboundary, computeVboundary)
        l_[index] = rmse(extract(map(constraints, q), 1))
        @printf("  n = %i, rmse = %e\n", n_[index], l_[index])
    end
    return (n_, l_)
end

function hconv(range::StepRange, np::Int=4)
    println("Testing h convergence")
    n_ = collect(range)
    l_ = zeros(size(n_))
    for index in CartesianIndices(n_)
        nh = n_[index] = 2^n_[index]
        params = Parameters((np, np), (nh, nh), urange, vrange, 2)
        l_[index] = rmse(extract(map(constraints, distribute(params, computePatch, computeUboundary, computeVboundary)), 1))
        l0 = index[1] > 1 ? l_[index[1] - 1] : 1
        @printf("  n = %3i, rmse[2h] / rmse[h] = %e\n", n_[index], l0 / l_[index])
    end
    return (n_, l_)
end

#-----------------------------------------
# Setup simulation grid parameters,
#-----------------------------------------
npoints  = (8, 8)
npatches = (4, 4)
urange   = (0.0, 1.0)
vrange   = (2.0, 3.0)
nfields  = 2
params   = Parameters(npoints, npatches, urange, vrange, nfields)
solver   = :trust_region
debug    = 1

#-----------------------------------------
# Check solution solution
#-----------------------------------------
# NOTE: We found the solution to be okayshily stable
# with the modified regularity conditions and after including the constraints.
# TODO: Check for convergence
# S = distribute(params, computePatch, computeUboundary, computeVboundary, debug)
# save("output/data/minkowski/constraints/minkowski_guv", extract(S, 1))
# save("output/data/minkowski/constraints/minkowski_grr", extract(S, 2))

#-----------------------------------------
# Now test h-p convergence of the error
#-----------------------------------------
np, ep = pconv(2:2:12, 3)
nh, eh = hconv(0:1:7, 4)
save("output/data/minkowski/minkowski_constraints_pconv", np, ep)
save("output/data/minkowski/minkowski_constraints_hconv", nh, eh)

# Notes
# [Adding noise]
# 1/ Adding noise to guv results in the non-linear solver failing to converge.
# 2/ Adding noise to psi allows the non-linear solver to converge.
# 3/ However the higher modes saturate around 1e-5, which explains why the constraint
#    violations are so large (i.e., once you take the second derivatives of the basis functions
#    their values at the boundary are orders of magnitude higher than in the interior.)
#    Filtering, at least in the Minkowski case is tricky, since most modes saturate around 1e-4. Filtering
#    the top-half is already quite aggresive but that doesn't work either. If we set all the modes below
#    say 1e-3 to zero, we get the constraints to the satisfied well.


text:
struct ExpressionKey{T <: ExpressionType, U <: Union{PSY.Component, PSY.System}} <:
       OptimizationContainerKey
    meta::String
end

function ExpressionKey(
    ::Type{T},
    ::Type{U},
    meta=CONTAINER_KEY_EMPTY_META,
) where {T <: ExpressionType, U <: Union{PSY.Component, PSY.System}}
    if isabstracttype(U)
        error("Type $U can't be abstract")
    end
    check_meta_chars(meta)
    return ExpressionKey{T, U}(meta)
end

get_entry_type(
    ::ExpressionKey{T, U},
) where {T <: ExpressionType, U <: Union{PSY.Component, PSY.System}} = T

get_component_type(
    ::ExpressionKey{T, U},
) where {T <: ExpressionType, U <: Union{PSY.Component, PSY.System}} = U

function encode_key(key::ExpressionKey)
    return encode_symbol(get_component_type(key), get_entry_type(key), key.meta)
end

Base.convert(::Type{ExpressionKey}, name::Symbol) = ExpressionKey(decode_symbol(name)...)

abstract type SystemBalanceExpressions <: ExpressionType end
abstract type RangeConstraintLBExpressions <: ExpressionType end
abstract type RangeConstraintUBExpressions <: ExpressionType end
abstract type CostExpressions <: ExpressionType end
struct ActivePowerBalance <: SystemBalanceExpressions end
struct ReactivePowerBalance <: SystemBalanceExpressions end
struct EmergencyUp <: ExpressionType end
struct EmergencyDown <: ExpressionType end
struct RawACE <: ExpressionType end
struct ProductionCostExpression <: CostExpressions end
struct ActivePowerRangeExpressionLB <: RangeConstraintLBExpressions end
struct ComponentActivePowerRangeExpressionLB <: RangeConstraintLBExpressions end
struct ReserveRangeExpressionLB <: RangeConstraintLBExpressions end
struct ActivePowerRangeExpressionUB <: RangeConstraintUBExpressions end
struct ReserveRangeExpressionUB <: RangeConstraintUBExpressions end
struct ComponentActivePowerRangeExpressionUB <: RangeConstraintUBExpressions end
struct ComponentReserveUpBalanceExpression <: ExpressionType end
struct ComponentReserveDownBalanceExpression <: ExpressionType end

should_write_resulting_value(::Type{<:ExpressionType}) = false
should_write_resulting_value(::Type{<:CostExpressions}) = true
should_write_resulting_value(::Type{RawACE}) = true

convert_result_to_natural_units(::Type{<:ExpressionType}) = false

ext: jl


text:
using EngEconomics

# Given
B = 5000
i = 0.015
A = B * i
ieq = 0.08 / 4
C = A + A * seriesPresentAmountFactor(ieq, 25) + B * presentWorthFactor(ieq, 25)


text:
# define and test barycentric velocity interpolation

using eFEMpart: squareMesh
using BenchmarkTools: @btime

include("../../src/part/SeepageForce.jl")  # loads functions to be tested

function interp_timing(N=5)
    # generate mesh
    Dom = [-1.0,1.0,-1.0,1.0]
    #N    = 5
    mesh = squareMesh(Dom,N,2)

    # generate u,v data
    u = zeros(length(mesh.xy))
    v = zeros(length(mesh.xy))
    for ti=1:length(mesh.xy)
        u[ti] = 1.0 - mesh.xy[ti].y^2
    end

    # generate particle position
    xpos=0.999; ypos=0.999
    p = Point2D(xpos,ypos)

    # print out FEM stats
    #FEMstats(mesh)

    # initialize memory so that allocation isn't an issue
    polygon    = [Point2D(0.0,0.0),Point2D(0.0,0.0),Point2D(0.0,0.0),Point2D(0.0,0.0)]
    extremeArr = Point2D(100_000.0,0.0)
    w          = zeros(4)
    uEl        = zeros(4)
    vEl        = zeros(4)
    a = zeros(3)
    b = zeros(3)
    c = zeros(3)
    d = zeros(3)

    # call
    #@code_warntype BarycentricVelocityInterp(mesh,u,v,p,polygon,extremeArr)
    #@btime BarycentricVelocityInterp($mesh,$u,$v,$p,$polygon,$extremeArr,$w,$uEl,$vEl,$a,$b,$c,$d)
    sVx,sVy = BarycentricVelocityInterp(mesh,u,v,p,polygon,extremeArr,w,uEl,vEl,a,b,c,d)

    if true
        println()
        println("expected u:    ",1-ypos^2)
        println("inteprolate u: ",sVx)
        println()
        println("expected v:    ",0.0)
        println("inteprolate v: ",sVy)
        println()
    end

    nothing
end

# timing of everything the same but using a cell list
function interp_timing_withCL(N=5,L=0.1)
    # generate mesh
    Dom = [-1.0,1.0,-1.0,1.0]
    #N    = 5
    mesh = squareMesh(Dom,N,2)

    # generate u,v data
    u = zeros(length(mesh.xy))
    v = zeros(length(mesh.xy))
    for ti=1:length(mesh.xy)
        u[ti] = 1.0 - mesh.xy[ti].y^2
    end

    # generate particle position
    xpos=1.0; ypos=1.0
    p = Point2D(xpos,ypos)

    # generate node cell list with p
    nodeList = [p]
    totalBounds = [-1.0,1.0,-1.0,1.0]
    particleCL = generateCellList(nodeList,totalBounds,L)

    # generate mesh cell list map
    meshCLmap = femGenerateMap(mesh,totalBounds,L)

    # initialize memory so that allocation isn't an issue
    polygon    = [Point2D(0.0,0.0),Point2D(0.0,0.0),Point2D(0.0,0.0),Point2D(0.0,0.0)]
    extremeArr = Point2D(100_000.0,0.0)
    w          = zeros(4)
    uEl        = zeros(4)
    vEl        = zeros(4)
    a          = zeros(3)
    b          = zeros(3)
    c          = zeros(3)
    d          = zeros(3)
    uInterp    = zeros(1)
    vInterp    = zeros(1)

    # call
    #@code_warntype BarycentricVelocityInterp_CL!(mesh,u,v,nodeList,polygon,particleCL,meshCLmap,extremeArr,w,uEl,vEl,a,b,c,d,uInterp,vInterp)
    @btime BarycentricVelocityInterp_CL!($uInterp,$vInterp,$mesh,$u,$v,$nodeList,$polygon,$particleCL,$meshCLmap,$extremeArr,$w,$uEl,$vEl,$a,$b,$c,$d)
    BarycentricVelocityInterp_CL!(uInterp,vInterp,mesh,u,v,nodeList,polygon,particleCL,meshCLmap,extremeArr,w,uEl,vEl,a,b,c,d)

    if false
        println()
        println("expected u:    ",1-ypos^2)
        println("inteprolate u: ",uInterp[1])
        println()
        println("expected v:    ",0.0)
        println("inteprolate v: ",vInterp[1])
        println()
    end

    # testing
    return mesh
end

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
struct Rect{T <: Number}
    m::Matrix{T}
    function Rect{T}(m::Matrix{T}) where T <: Number
        @assert size(m) == (2, 2) && all(lb(m) .< ru(m)) "Invalid values."
        new(m)
    end
    Rect{T}(lx::T, ly::T, rx::T, ry::T) where T <: Number =
        new(Matrix([min(lx, rx) max(lx, rx); min(ly, ry) max(ly, ry)]))
    Rect{T}(r::Rect) where T <: Number = new(Matrix{T}(r.m))
end

Rect(m::Matrix{T}) where T <: Number = Rect{T}(m)

function Rect(lx::Number, ly::Number, rx::Number, ry::Number)
    t = promote(lx, ly, rx, ry)
    T = typeof(t[1])
    return Rect{T}(t...)
end

Base.convert(::Type{Rect{T}}, r::Rect{S}) where {T <: Number, S <: Number} =
    Rect{T}(Matrix{T}(r.m))

Base.promote_rule(::Type{Rect{T}},
                  ::Type{Rect{S}}) where {T <: Number, S <: Number} =
    Rect{promote_type(T, S)}

Base.show(io::IO, r::Rect) =
    print(io, "Rect:[$(lx(r)) $(ly(r)) $(rx(r)) $(ry(r))]")

lb(r) = lb(r.m)
ru(r) = ru(r.m)
lb(m::Matrix) = m[:, 1]
ru(m::Matrix) = m[:, 2]
lx(r) = r.m[1, 1]
ly(r) = r.m[2, 1]
rx(r) = r.m[1, 2]
ry(r) = r.m[2, 2]

coord(r, axis) = r.m[axis, :]
x(r) = coord(r, 1)
y(r) = coord(r, 2)
c_lo(r, axis)   = r.m[axis, 1]
c_hi(r, axis)   = r.m[axis, 2]

hlines(r::Rect) =
    [Line(lx(r), ly(r), rx(r), ly(r)), Line(lx(r), ry(r), rx(r), ry(r))]
vlines(r::Rect) =
    [Line(lx(r), ly(r), lx(r), ry(r)), Line(rx(r), ly(r), rx(r), ry(r))]
lines(r::Rect)  = [hlines(r)..., vlines(r)...]
olines(r::Rect) =
    (ls = lines(r); [ls[1], ls[4], reverse(ls[2]), reverse(ls[3])])
diags(r::Rect)  =
    [Line(lx(r), ly(r), rx(r), ry(r)), Line(rx(r), ly(r), lx(r), ry(r))]
cg(r::Rect) = div(diags(r)[1], 1//2)

function Base.union(r1::Rect, r2::Rect)
    l = min.(lb(r1), lb(r2))
    r = max.(ru(r1), ru(r2))
    return Rect(l[1], l[2], r[1], r[2])
end

Base.union(r1::Rect, r2::Rect, y...) = Base.union(r1, union(r2, y...))
Base.union(r::Rect) = r

Base.intersect(r1::Rect, r2::Rect) = intersect(promote(r1, r2)...)

function Base.intersect(r1::Rect{T}, r2::Rect{T}) where T <: Number
    l = max.(lb(r1), lb(r2))
    r = min.(ru(r1), ru(r2))
    l1 = l .+ pcTol(T)
    any(l1 .>= r) && return nothing
    return Rect(l[1], l[2], r[1], r[2])
end

Base.:(==)(r1::Rect{T}, r2::Rect{T}) where {T <: Number} =
    all(abs.(r1.m - r2.m) .<= pcTol(T))
Base.:(==)(r1::Rect, r2::Rect) = ==(promote(r1, r2)...)

inside(p::Tuple{T, T}, r::Rect{T}) where T <: Number =
    all(r.m[:, 1] .<= p .<= r.m[:, 2])

function inside(p::Tuple{T1, T2}, r::Rect{S}) where {T1 <: Number,
                                                     T2 <: Number,
                                                     S  <: Number}
    ST = promote_type(S, T1, T2)
    inside(convert(Tuple{ST, ST}, p), convert(Rect{ST}, r))
end

inside(ri::Rect, ro::Rect) = intersect(ri, ro) == ri

intersects(r1::Rect, r2::Rect) = intersect(r1, r2) != nothing

function intersects(r::Rect, l::Line)
    ml = l.m
    (inside((ml[1, 1], ml[2, 1]), r) || inside((ml[1, 2], ml[2, 2]), r)) &&
        return true
    return  any(intersects.(lines(r), Ref(l)))
end

h(r::Rect) = ry(r) - ly(r)
w(r::Rect) = rx(r) - lx(r)
area(r::Rect) = h(r) * w(r)

perimeter(r::Rect) =  (s = h(r) + w(r); s + s)
to_plot_shape(r::Rect) =
    ([lx(r), rx(r), rx(r), lx(r)], [ly(r), ly(r), ry(r), ry(r)])

xsort(r1::Rect, r2::Rect, reverse=false) =
    sortr(r1, r2, reverse=reverse, axis=1)
ysort(r1::Rect, r2::Rect, reverse=false) =
    sortr(r1, r2, reverse=reverse, axis=2)

sortr(r1::Rect, r2::Rect; reverse=false, axis=1) =
    (!reverse && r1.m[axis, 1] > r2.m[axis, 1]) ? (r2, r1) : (r1, r2)

has_x_overlap(r1::Rect, r2::Rect; isopen=false) =
    has_overlap(r1::Rect, r2::Rect, axis=1, isopen=isopen)
has_y_overlap(r1::Rect, r2::Rect; isopen=false) =
    has_overlap(r1::Rect, r2::Rect, axis=2, isopen=isopen)

function has_overlap(ri1::Rect{T1}, ri2::Rect{T2};
                     axis=1, isopen=false) where {T1, T2}
    r1, r2 = promote(ri1, ri2)
    T = typeof(r1.m[1, 1])
    tol = isopen ? pcTol(T) : zero(T)
    m = copy(r1.m)
    m[axis, 1] += tol
    m[axis, 2] -= tol
    rt = Rect(m)
    r1, r2 = sortr(rt, r2, reverse=false, axis=axis)
    return r1.m[axis, 1] <= r2.m[axis, 1] <= r1.m[axis, 2]
end

"""
```
    projectX(r1::Rect, r2::Rect) -> (left, overlap, right)
```
Projects the rectangles along the X-axis and returns three parts of rectangles.

`left`: The left segment of the projection
`overlap`: If there is any overlap between the rectangles
`right`: The right segment of the projection

Each portion is returned as a tuple.

If the rectangle is part of the first rectangle, it's returned as the first
element of the tuple. `nothing` is returned for a part when a portion is not
available.
"""
projectX(r1::Rect, r2::Rect) = project(r1, r2, axis=1)

"""
```
    projectY(r1::Rect, r2::Rect) -> (bottom, overlap, top)
```
Projects the rectangles along the Y-axis and returns three parts of rectangles.

`bottom`: The bottom segment of the projection
`overlap`: If there is any overlap between the rectangles
`top`: The top segment of the projection

Each portion is returned as a tuple.

If the rectangle is part of the first rectangle, it's returned as the first
element of the tuple. `nothing` is returned for a part when a portion is not
available.
"""
projectY(r1::Rect, r2::Rect) = project(r1, r2, axis=2)

project(r1::Rect, r2::Rect; axis::Int=1) = project(promote(r1, r2)...; axis=axis)
function project(r1::Rect{T}, r2::Rect{T}; axis::Int=1) where T <: Number
    # l=low and h=high
    lr, hr = sortr(r1, r2, axis=axis)
    flip = lr != r1

    cl = coord(lr, axis)
    ch = coord(hr, axis)
    axl = sort([cl[1], cl[2], ch[1], ch[2]])

    # Low segment
    if c_lo(lr, axis) + pcTol(T) < c_lo(hr, axis)
        lm = copy(lr.m)
        lm[axis, 1] = axl[1]
        lm[axis, 2] = axl[2]
        l = Rect{T}(lm)
        low = !flip ? (l, nothing) : (nothing, l)
        # Non-overlap
        l == lr &&
            return low, (nothing, nothing), !flip ? (nothing, hr) : (hr, nothing)
    else
        low = nothing, nothing
    end

    #Overlap segment
    om1, om2 = copy(r1.m), copy(r2.m)
    om1[axis, 1] = om2[axis, 1] = axl[2]
    om1[axis, 2] = om2[axis, 2] = axl[3]

    overlap = Rect{T}(om1), Rect{T}(om2)

    #High segment
    if -pcTol(T) <= c_hi(r1, axis) - c_hi(r2, axis) <= pcTol(T)
        high = nothing, nothing
    else
        hm = c_hi(r1, axis) > c_hi(r2, axis) ? copy(r1.m) : copy(r2.m)
        hm[axis, 1] = axl[3]
        hm[axis, 2] = axl[4]
        high = c_hi(r1, axis) > c_hi(r2, axis) ? (Rect{T}(hm), nothing) :
                                                 (nothing, Rect{T}(hm))
    end
    return low, overlap, high
end

"""
```
    visibleX(r1::Rect, r2::Rect) -> Rect
    visibleY(r1::Rect, r2::Rect) -> Rect
```
Projects the rectangles along the X-axis (Y-axis) and returns a rectangle area
which is completely visible from both rectangles.

`nothing` is returned when there is no overlap along the X-axis.
"""
visibleX(r1::Rect, r2::Rect) = visible(r1::Rect, r2::Rect, axis=1)
visibleY(r1::Rect, r2::Rect) = visible(r1::Rect, r2::Rect, axis=2)

visible(r1::Rect, r2::Rect; axis::Int=1) = visible(promote(r1, r2)...; axis=axis)

function visible(r1::Rect{T}, r2::Rect{T}; axis::Int=1) where T <: Number
    l, ox, r = project(r1, r2, axis=axis)
    ox == (nothing, nothing) && return nothing
    saxis = (axis == 1 ? 2 : 1)
    tr1, tr2 = sortr(ox[1], ox[2], axis=saxis)
    m = copy(tr1.m)
    m[saxis, 1] = tr1.m[saxis, 2]
    m[saxis, 2] + pcTol(T) > tr2.m[saxis, 1] &&
        return Line(lx(tr2), ly(tr2), lx(tr2), ry(tr2))
    m[saxis, 2] = tr2.m[saxis, 1]
    return Rect(m[1, 1], m[2, 1], m[1, 2], m[2, 2])
end

function get_overlapped_dist(l, h, f, a1, a2)
    # l = low and h = high
    lr, adl = l[1] === nothing ? (l[2], a2) : (notvoid(l[1]), a1)
    wlal    = lr   === nothing ? zero(a1)   : f(notvoid(lr))*area(notvoid(lr))

    hr, adh = h[1] === nothing ? (h[2], a2) : (notvoid(h[1]), a1)
    whah    = hr   === nothing ? zero(a1)   : f(notvoid(hr))*area(notvoid(hr))

    d = Float64(wlal)/Float64(adl) + Float64(whah)/Float64(adh)
    return d/2.0
end

"""
```
    avg_min_dist(r1::Rect, r2::Rect) -> dx::Float64, dy::Float64
```
Rectangles are essentially point sets. Hence, one can perceive existence of a
minimum distance of one point in `r1` from `r2`. Similar, distance would also
exist for every point in `r2` from `r1`. While, technically Euclidean distance
metric can exist, the computation is fairly cumborsome. Here, we use the city
block distance or L1-metric.

`dx`: The distance in the x-direction
`dy`: The distance in the y-direction

The minimum distance will be `zero` when the rectangles are intersecting.
The distance also will be lower in a specific direction if there is an overlap of
 the rectangles in that direction
"""
function avg_min_dist(r1::Rect, r2::Rect)
    intersect(r1, r2) != nothing && return 0.0, 0.0
    h1, h2 = Float64(h(r1)),    Float64(h(r2))
    w1, w2 = Float64(w(r1)),    Float64(w(r2))
    a1, a2 = Float64(area(r1)), Float64(area(r2))

    hg = Float64(min(abs(ry(r2) - ly(r1)), abs(ry(r1) - ly(r2))))
    wg = Float64(min(abs(rx(r2) - lx(r1)), abs(rx(r1) - lx(r2))))

    dy = hg + (h1*a1 + h2*a2)/(2.0*(a1 + a2))
    dx = wg + (w1*a1 + w2*a2)/(2.0*(a1 + a2))

    l, ox, r = projectX(r1, r2)
    ox != (nothing, nothing) &&
        return get_overlapped_dist(l, r, w, a1, a2), dy
    b, oy, t = projectY(r1, r2)
    oy != (nothing, nothing) && return dx, get_overlapped_dist(b, t, h, a1, a2)
    return dx, dy
end

"""
```
    min_dist(r1::Rect, r2::Rect) -> dx, dy
```
Minimum distance or gap between two rectangles.

`dx`: The distance in the x-direction
`dy`: The distance in the y-direction

The minimum distance will be `zero` when the rectangles are overlapping in a
direction.
"""
function min_dist(r1::Rect{T}, r2::Rect{T}) where T <: Number
    r1, r2 = xsort(r1, r2)
    dx = has_x_overlap(r1, r2) ? zero(T) : lx(r2) - rx(r1)
    r1, r2 = ysort(r1, r2)
    dy = has_y_overlap(r1, r2) ? zero(T) : ly(r2) - ry(r1)
    return dx, dy
end
min_dist(r1::Rect, r2::Rect) = min_dist(promote(r1, r2)...)

mutable struct OrderedRectMap{T <: Number, V, D}
    data::IntervalTree{T, IntervalTree{T, V}}
    reverseMax::T
    function OrderedRectMap{T, V, D}(
        ;reverseMax::T1=zero(T)) where {T <: Number, V, D, T1 <: Number}
        @assert reverseMax >= zero(T) "Invalid max value for reverse ordering"
        new(IntervalTree{T, IntervalTree{T, V}}(), convert(T, reverseMax))
    end
end

const OrderedRectMapX{T, V} = OrderedRectMap{T, V, 1}
const OrderedRectMapY{T, V} = OrderedRectMap{T, V, 2}

function create_ordered_map(rects::AbstractVector{Rect{T}},
                            values::AbstractVector{V};
                            dir::Int=1,
                            reverseMax::T1=zero(T)) where {T <: Number,
                                                           V,
                                                           T1 <: Number}
    map = OrderedRectMap{T, V, dir}(reverseMax=reverseMax)
    nr = iterate(rects)
    nv = iterate(values)
    odir = dir == 1 ? 2 : 1
    while nr !== nothing
        rect, itr = nr
        v2,   itv = nv
        nr = iterate(rects, itr)
        nv = iterate(values, itv)
        insert_rect!(map, rect, v2)
    end
    return map
end

function Base.intersect(orm::OrderedRectMap{T1, V, D},
                        rect::Rect{T2}, dX::T1, dY::T1;
                        dirX=0, dirY=0,
                        fstop::Function=(x) -> false) where {T1 <: Number,
                                                             T2 <: Number, V, D}
    dl = dirX > 0  ? zero(T1) : -dX
    dr = dirX < 0  ? zero(T1) :  dX
    dd = dirY > 0  ? zero(T1) : -dY
    du = dirY < 0  ? zero(T1) :  dY

    r = convert(Rect{T1}, rect)
    ors, ovs = Vector{Rect{T1}}(), Vector{V}()
    while true
        tr = Rect(lx(r) + dl, ly(r) + dd, rx(r) + dr, ry(r) + du)
        rs, vs = intersect(orm, tr)
        length(rs) == 0 && return rs, vs
        ttr = union(rs...)
        fstop(ttr) && return ors, ovs
        r == ttr && return rs, vs
        r, ors, ovs = ttr, rs, vs
    end
end

function Base.intersect(orm::OrderedRectMap{T1, V, D},
                        r::Rect{T2};
                        isopen = true) where {T1 <: Number, T2 <: Number, V, D}
    rect = convert(Rect{T1}, r)
    dir = D
    odir = dir == 1 ? 2 : 1
    r1 = coord(rect, dir)
    if orm.reverseMax != zero(T1)
        r1[1], r1[2] = (orm.reverseMax - r1[2]), (orm.reverseMax - r1[1])
    end
    r2 = coord(rect, odir)
    tol = isopen ? pcTol(T1) : zero(T1)
    imv1 = intersect(orm.data, Interval(r1[1] + tol, r1[2] - tol))
    retr = Vector{Rect{T1}}()
    retv = Vector{V}()
    nimv1 = iterate(imv1)
    while nimv1 !== nothing
        v1, iv1 = nimv1
        imv2 = intersect(v1[2], Interval(r2[1] + tol, r2[2] - tol))
        nimv2 = iterate(imv2)
        while nimv2 !== nothing
            v2, iv2 = nimv2
            m = zeros(T1, (2,2))
            m[ dir, 1], m[ dir, 2] =
                (orm.reverseMax == zero(T1)) ? (v1[1].lo, v1[1].hi) :
                (orm.reverseMax - v1[1].hi, orm.reverseMax - v1[1].lo)
            m[odir, 1], m[odir, 2] = v2[1].lo, v2[1].hi
            push!(retr, Rect{T1}(m))
            push!(retv, v2[2])
            nimv2 = iterate(imv2, iv2)
        end
        nimv1 = iterate(imv1, iv1)
    end
    return retr, retv
end

function insert_rect!(orm::OrderedRectMap{T1, V, D},
    r::Rect{T2}, v::V) where {T1 <: Number, T2 <: Number, V, D}
    rect = convert(Rect{T1}, r)
    dir = D
    odir = dir == 1 ? 2 : 1
    r1 = coord(rect, dir)
    if orm.reverseMax != zero(T1)
        r1[1], r1[2] = (orm.reverseMax - r1[2]), (orm.reverseMax - r1[1])
    end
    imv = get!(orm.data, Interval(r1[1], r1[2]), IntervalTree{T1, V}())
    r2 = coord(rect, odir)
    tv = get!(imv, Interval(r2[1], r2[2]), v)
    tv === v && return nothing
    tv = delete!(imv, Interval(r2[1], r2[2]))
    imv[Interval(r2[1], r2[2])] = v
    return tv
end

function delete_rect!(orm::OrderedRectMap{T1, V, D},
    r::Rect{T2}) where {T1 <: Number, T2 <: Number, V, D}
    rect = convert(Rect{T1}, r)
    dir = D
    odir = dir == 1 ? 2 : 1
    r1 = coord(rect, dir)
    if orm.reverseMax != zero(T1)
        r1[1], r1[2] = (orm.reverseMax - r1[2]), (orm.reverseMax - r1[1])
    end
    r2 = coord(rect, odir)
    imv = get(orm.data, Interval(r1[1], r1[2]), IntervalTree{T1, V}())
    isempty(imv) && return nothing
    ret = delete!(imv, Interval(r2[1], r2[2]))
    isempty(imv) && delete!(orm.data, Interval(r1[1], r1[2]))
    return ret === nothing ? ret : ret[2]
end


text:
module MarketEquilibria

using Random
using JuMP
using Ipopt

export homogeneousfisher, tatonnement, exchange, production

"""
    homogeneousfisher(endowments, A, supplies, form, ρ)

§6.2: Fisher Model with Homogeneous Consumers. `form` specifies which utility
function consumers use, among `:linear`, `:CES`, `:cobb_douglas`, and `:leontief`. `ρ` is the parameter for
`:CES`.
"""
function homogeneousfisher(endowments::Array{Float64,1},
                           A::Array{Float64,2},
                           supplies::Array{Float64,1},
                           form=:linear,
                           ρ::Union{Float64,Nothing}=0.5)::Tuple{Array{Float64,2},
                                                                 Array{Float64,1}}
    @assert form in [:linear, :CES, :cobb_douglas, :leontief] "Unknown form"
    if form == :cobb_douglas
        ρ = 1e-8          # Ideally 0
    end

    (n, m) = size(A)
    @assert size(endowments) == (n, ) "Dim mismatch between A and endowments"
    @assert size(supplies) == (m, ) "Dim mismatch between A and supplies"

    model = Model(Ipopt.Optimizer)
    set_optimizer_attribute(model, "print_level", 0)

    @variable(model, X[1:n, 1:m], lower_bound=0)
    @constraint(model, Supply[j in 1:m], sum(X[:, j]) ≤ supplies[j])

    if form == :linear || ρ > .97
        if form != :linear
            @warn "ρ > .97 tends to fail, so assuming linear"
        end
        # Very bad results if you allow the default scaling 🐸☕
        set_optimizer_attribute(model, "nlp_scaling_method", "none")
        @NLobjective(model, Max, sum(endowments[i] * log(sum(A[i, j] * X[i, j] for j in 1:m)) for i in 1:n))
    elseif form == :leontief
        @variable(model, t)
        @objective(model, Max, t)
        @NLconstraint(model, tLessThan[j in 1:m],
                    t ≤ sum(endowments[i] * log(A[i, j] * X[i, j]) for i in 1:n))
    else
        @NLobjective(model, Max, sum(endowments[i] * log(sum(A[i, j] * X[i, j] ^ ρ for j in 1:m) ^ (1 / ρ)) for i in 1:n))
    end

        # Explicit Cobb–Douglas; very poorly conditioned. I couldn't get it to converge.
        # @NLobjective(model, Max, sum(endowments[i] * log(prod(X[i, j] ^ A[i, j] for j in 1:m)) for i in 1:n))

        # Explicit Leontief is as follows, but minimum() not allowed by JuMP.
        # Got around this using the linearization trick.
        # @NLobjective(model, Max, sum(endowments[i] * log(minimum(A[i, j] * X[i, j] for j in 1:m)) for i in 1:n))

        # Attempt at "softmin" differential substitute; returns errors.
        # @NLobjective(model, Max, sum(endowments[i] * log(-log(sum(exp(-A[i, j] * X[i, j]) for j in 1:m))) for i in 1:n))

    optimize!(model)
    return value.(X), dual.(Supply)
end

"""
    tatonnement(endowments, A, ρ, π0;
                η=1e-4, δ=0.5, α=0.5, verbose=false, maxit=1000, tol=1e-5)

§6.3.1.B: The Discrete Tâtonnement Process. Simulate a bidding process in search
of equilibrium prices for CES exchange.

Price update step is quite different from the book: I *multiply* prices
entrywise by the excess demand: `π .*= exp.(δ * Z * k ^ -α)`. This works pretty well.

`η` is a lower bound on the prices before normalization.
"""
function tatonnement(endowments ::AbstractArray{<:AbstractFloat,2},
                     A          ::AbstractArray{<:AbstractFloat,2},
                     ρ          ::AbstractArray{<:AbstractFloat,1},
                     π0         ::Union{AbstractArray{<:AbstractFloat,1}, Nothing};
                     η          ::AbstractFloat=1e-4,
                     δ          ::AbstractFloat=0.5,
                     α          ::AbstractFloat=0.4,
                     verbose    ::Bool=false,
                     maxit      ::Int=2000,
                     tol        ::AbstractFloat=1e-5)::Tuple{Array{Float64,1},
                                                             Array{Float64,2}}

    (n, m) = size(A)
    @assert (n, m) == size(endowments)
    W = endowments      # sorry

    @assert (n, )  == size(ρ)
    @assert all(ρ .< 1)
    @assert all(0 .< (tol, η, δ, α) .< 1)

    if π0 == nothing
        π0 = rand(m)
    end

    @assert (m, )  == size(π0)

    total_supply = sum(W, dims=1)
    # Express utility as fractions of endowment
    A .*= total_supply
    # So we can normalize W
    W ./= total_supply

    σ = 1 ./ (1 .- ρ)
    π_LB = η / (4 * n)
    π_UB = 1 + π_LB
    π_LB_tight = η / (2 * n)

    π = copy(π0)
    D = zeros(Float64, n, m)
    Z = zeros(Float64, m)

    for k in 1:maxit
        # Oracle that computes player demand under CES (eq. 6.7, p. 149)
        for i in 1:n, j in 1:m
            D[i, j] = A[i, j] ^ σ[i] * sum(π[k] * W[i, k] for k in 1:m) /
                       (π[j] ^ σ[i] * sum((A[i, k] * π[k] ^ -ρ[i]) ^ σ[i] for k in 1:m))
        end

        for j in 1:m
            # Total supply, after normalization, is ones
            Z[j] = sum(D[i, j] for i in 1:n) - 1
        end

        if NaN in Z || all(-tol .< Z .< tol)
            break
        end

        if verbose
            println("Iteration $k")
            println("  Current prices: ", round.(π, digits=5))
            println("  Excess demand:  ", round.(Z, digits=5))
        end

        # This equilibrium step process is completely different from the book.
        # Just something that makes intuitive sense. Book is unclear on how to
        # calculate step size.
        if all(π_LB .≤ π .≤ π_UB)
            verbose ? println("  Step toward equilibrium") : nothing
            π .*= exp.(δ * Z * k ^ -α)
        else                             # Bottom of p. 144
            verbose ? println("  Regularization step") : nothing
            for j in 1:m
                if π[j] > 1
                    π[j] = 1
                elseif π[j] < π_LB_tight
                    π[j] = π_LB_tight
                end
            end
        end
    end

    return π ./ sum(π), D
end

"""
    exchange(endowments, A, ρ=:linear)

§6.4: Convex programs for linear and CES exchange economies. Note that unlike `homogeneousfisher()`,
here `ρ` is a vector of CES parameters used by each player. Or pass `ρ=:linear` for linear model.
"""
function exchange(endowments::Array{Float64,2},
                  A::Array{Float64,2},
                  ρ::Union{Array{Float64,1},Symbol}=:linear)::Tuple{Array{Float64,1},
                                                                    Array{Float64,2}}
    (n, m) = size(A)
    @assert size(endowments) == (n, m)         "Dim mismatch between A and endowments"

    model = Model(Ipopt.Optimizer)
    set_optimizer_attribute(model, "print_level", 0)

    if ρ == :linear
        # set_optimizer_attribute(model, "nlp_scaling_method", "none")

        @variable(model, ψ[1:m])
        @variable(model, X[1:n, 1:m], lower_bound=0)
        @NLconstraint(model, IndividualRationality[i in 1:n, j in 1:m],
                        sum(A[i, k] * X[i, k] for k in 1:m) ≥
                        A[i, j] * sum(endowments[i, k] * exp(ψ[k] - ψ[j]) for k in 1:m))
        @constraint(model, Supply[j in 1:m], sum(X[i, j] for i in 1:n) ==
                                             sum(endowments[i, j] for i in 1:n))

        optimize!(model)

        prices = exp.(value.(ψ))
        demands = value.(X)

    else
        @assert size(ρ) == (n, )                   "Dim mismatch between A and ρ"
        @assert all(-1 .≤ ρ .< 0)                  "Need ρ ∈ [-1, 0)"

        @variable(model, σ[1:m], lower_bound=0)
        @NLconstraint(model, Supply[j in 1:m], sum(A[i, j] ^ (1 / (1 - ρ[i])) *
                                               sum(σ[k] ^ 2 * endowments[i, k] for k in 1:m) /
                                               (
                                                   σ[j] ^ ((ρ[i] - 2) / (1 - ρ[i])) *
                                                   sum(A[i, k] ^ (1 / (1 - ρ[i])) * σ[k] ^ ((-2 * ρ[i]) / (1 - ρ[i]))
                                                       for k in 1:m)
                                               )
                                               for i in 1:n) ≤ σ[j] * sum(endowments[i, j] for i in 1:n))

        optimize!(model)

        prices = value.(σ) .^ 2

        demands = zeros(n, m)
        for i in 1:n, j in 1:m
        demands[i, j] = A[i, j] ^ (1 / (1 - ρ[i])) * sum(prices[k] * endowments[i, k] for k in 1:m) /
                        (
                            prices[j] ^ (1 / 1 - ρ[i]) *
                            sum(A[i, k] ^ (1 / (1 - ρ[i])) * prices[k] ^ (-ρ[i] / (1 - ρ[i])) for k in 1:m)
                        )
        end
    end

    return prices, demands
end

"""
    exchange(endowments, A_consumers, A_producers,
             ρ_consumers, ρ_producers, o_producers)

§6.6: Models with production. `ρ_consumers` describes
each consumer's CES utility function, while `ρ_producers` describes the CES
production function. Convergence is very iffy, but best when `A_consumers` is
tall and entries of `ρ` vectors are not too close to 0 or 1.
"""
function production(endowments::Array{Float64,2},
                    A_consumers::Array{Float64,2},
                    A_producers::Array{Float64,2},
                    ρ_consumers::Array{Float64,1},
                    ρ_producers::Array{Float64,1},
                    o_producers::Array{Int,1})::Tuple{Array{Float64,1},
                                                      Array{Float64,1},
                                                      Array{Float64,2}}

    (n, m) = size(A_consumers)
    (l, ) = size(ρ_producers)

    @assert size(A_producers) == (l, m)   "Dim mismatch between A_consumers and A_producers"
    @assert size(endowments) == (n, m)    "Dim mismatch between A_consumers and endowments"
    @assert size(ρ_consumers) == (n, )    "Dim mismatch between A_consumers and ρ_consumers"
    @assert size(o_producers) == (l, )    "Dim mismatch between o and ρ_producers"
    @assert all(0 .< ρ_consumers .< 1)    "Need ρ ∈ (0, 1)"
    @assert all(0 .< ρ_producers .< 1)    "Need ρ ∈ (0, 1)"

    A = vcat(A_producers, A_consumers)
    w = hcat(endowments, zeros(n, l))
    ρ = vcat(ρ_consumers, ρ_producers)
    σ = 1 ./ (1 .- ρ)
    o = vcat(o_producers, m .+ (1:n))

    # Comparison idx since comparison forbidden in lincon
    O = [o[k] == j for k in 1:l, j in 1:m]

    model = Model(Ipopt.Optimizer)
    set_optimizer_attribute(model, "print_level", 0)
    set_optimizer_attribute(model, "nlp_scaling_method", "none")

    @variable(model, ψ[1:m + n], lower_bound=0)
    @variable(model, x[1:m + n], lower_bound=0)
    @variable(model, Z[1:l + n, 1:m], lower_bound=0)
    @variable(model, q[1:l + n], lower_bound=0)

    @NLconstraint(model, ConsumerRationality[i in 1:n],                   # 6.10
                    exp(ψ[m + i]) * x[m + i] ≥ # or ==
                    sum(exp(ψ[j]) * w[i, j] for j in 1:m))
    @NLconstraint(model, ProducerSolvency[k in 1:l + n],                  # 6.11
                    q[k] ≤
                    sum(A[k, j] * x[j] ^ ρ[k] for j in 1:m) ^ (1 / ρ[k]))
    @NLconstraint(model, ProducerRationality[k in 1:l + n],               # 6.12
                    exp(ψ[o[k]] * (1 - σ[k])) ≥
                    sum(A[k, j] ^ σ[k] * exp(ψ[j] * (1 - σ[k])) for j in 1:m))
    @constraint(model, MarketClearing[j in 1:m],                          # 6.13
                    sum(Z[k, j] for k in 1:l + n) ≤
                    sum(w[i, j] for i in 1:n) + sum(O[k, j] * q[k] for k in 1:l))
                                            # Equiv
                                            # + sum((o[k] == j) * q[k] for k in 1:l + n))
    @constraint(model, Supply[i in 1:n], x[m + i] ≤ q[l + i])             # 6.14

    #=  Final constraint given in final paragraph of 154. This should be redundant, but
        aids computation.                   =#
    @NLconstraint(model, ZeroProfit[k in 1:l + n],
                    sum(exp(ψ[j]) * Z[k, j] for j in 1:m) == exp(ψ[o[k]]) * q[k])

    optimize!(model)
    prices = exp.(value.(ψ)[1:m])

    # Player utility, specifically, demand of each player for her own utility item.
    demands = value.(x)[m + 1:end]
    inputs = value.(Z)[l + 1:end, :]

    return prices, demands, inputs
end

end

max_stars_repo_name: maxkapur/MarketEquilibria

max_issues_repo_name: maxkapur/MarketEquilibria

max_forks_repo_name: maxkapur/MarketEquilibria


text:
############################################################
# joLinearOperator - subtype operators #####################
############################################################

# joKron
include("joLinearOperatorConstructors/joKron.jl")

# joCoreBlock
include("joLinearOperatorConstructors/joCoreBlock.jl")
# joCoreBlock derivatives: joBlockDiag joDict joStack joBlock
include("joLinearOperatorConstructors/joCoreBlockConstructors.jl")

############################################################
# joLinearOperator - constructor-only operators ############
############################################################


text:
function bash_jobstate(jobdir::String)
    job = load_job(jobdir)
    n = now()
    u = username()
    i = last_running_calculation(jobdir)
    i === nothing && return false
    l = job[i]
    codeexec = l.exec.exec
    try
        pids = parse.(Int, split(read(`pgrep $codeexec`, String)))
        if isempty(pids)
            return (-1, Jobs.Completed)
        end
        pwds = map(x -> abspath(split(strip(read(`pwdx $x`, String)))[end]), pids)
        id = findfirst(isequal(abspath(job)), pwds)
        if id === nothing
            return (-1, Jobs.Completed)
        else
            return (pids[id], Jobs.Running)
        end
    catch
        return (-1, Jobs.Completed)
    end
end

function bash_queue!(q, init)
    if init
        jobdirs = filter(exists_job, readlines(RUNNING_JOBS_FILE))
    else
        jobdirs = keys(q)
    end
    to_running = String[]
    for j in jobdirs
        info = bash_jobstate(j)
        q[j] = info
        if info[2] == Jobs.Running
            push!(to_running, j)
        end
    end
    write(RUNNING_JOBS_FILE, join(to_running, "\n"))
    return q
end

function bash_submit(j::String)
    cd(j)
    run(Cmd(`bash job.tt`, detach=true), wait=false)
    open(RUNNING_JOBS_FILE, "a") do f
        write(f, j * "\n")
    end
    return (-1, Jobs.Submitted)
end


text:
using Query
using Test
using QuerySQLite
using SQLite: DB, drop!, execute!, Stmt
using QueryTables

@testset "SQLite tutorial" begin

filename = joinpath(@__DIR__, "Chinook_Sqlite.sqlite")
database = Database(filename)
database2 = Database(DB(filename))

@test database.Track |>
    @map({_.TrackId, _.Name, _.Composer, _.UnitPrice}) |>
    collect |>
    first |>
    propertynames == (:TrackId, :Name, :Composer, :UnitPrice)

@test (database2.Customer |>
    @map({_.City, _.Country}) |>
    @orderby(_.Country) |>
    DataTable).Country[1] == "Argentina"

@test database.Customer |>
    @map({_.City}) |>
    @unique() |>
    collect |>
    length == 53

@test database.Track |>
    @map({_.TrackId, _.Name}) |>
    @take(10) |>
    collect |>
    length == 10

@test first((database.Track |>
    @map({_.TrackId, _.Name}) |>
    @take(10) |>
    @drop(10) |>
    DataTable).Name)  == "C.O.D."

@test first((
    @from i in database.Track begin
        @orderby descending(i.Bytes), i.Name
        @select {i.TrackId, i.Name, i.Bytes}
        @collect DataTable
    end
    ).Bytes) == 1059546140

@test database.Artist |>
    @join(database.Album, _.ArtistId, _.ArtistId, {_.ArtistId, __.AlbumId}) |>
    DataTable |>
    length == 347

@test (database.Track |>
    @groupby(_.AlbumId) |>
    @map({AlbumId = key(_), Count = length(_.AlbumId)}) |>
    collect |>
    first).Count == 10

end

@testset "Systematic tests" begin

filename = joinpath(@__DIR__, "test.sqlite")
connection = DB(filename)
execute!(Stmt(connection, """DROP TABLE IF EXISTS test"""))
execute!(Stmt(connection, """
    CREATE TABLE test (
        a Int,
        b Int,
        c Int
    )"""))
execute!(Stmt(connection, """
    INSERT INTO test VALUES(0, 1, -1)
"""))
database = Database(connection)
result =
    database.test |>
    @map({
        c = _.a == _.b,
        d = _.a != _.b,
        e = !(_.a),
        f = _.a & _.b,
        g = _.a | _.b,
        h = _.a * _.b,
        i = _.a + _.b,
        j = _.a % _.b,
        k = abs(_.c),
        l = _.a in (0, 1)
    }) |>
    collect |>
    first

@test result.c == 0
@test result.d == 1
@test result.e == 1
@test result.f == 0
@test result.g == 1
@test result.h == 0
@test result.i == 1
@test result.j == 0
@test result.k == 1
@test result.l == 1

drop!(connection, "test")

@test_throws ArgumentError Database("file.not_sqlite")

end


text: println("hello julia")


text:
using JuMP, EAGO, Gurobi

m = Model(with_optimizer(EAGO.Optimizer, relaxed_optimizer = Gurobi.Optimizer(OutputFlag=0)))

# ----- Variables ----- #
x_Idx = Any[2, 3, 4, 5, 6, 7]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[2], 1.0e-7)
JuMP.set_upper_bound(x[2], 0.2)
JuMP.set_lower_bound(x[3], 1.0e-7)
JuMP.set_upper_bound(x[3], 0.2)
JuMP.set_lower_bound(x[4], 1.0e-7)
JuMP.set_upper_bound(x[4], 0.4)
JuMP.set_lower_bound(x[5], 1.0e-7)
JuMP.set_upper_bound(x[5], 0.4)
JuMP.set_lower_bound(x[6], 1.0e-7)
JuMP.set_upper_bound(x[6], 0.4)
JuMP.set_lower_bound(x[7], 1.0e-7)
JuMP.set_upper_bound(x[7], 0.4)

# ----- Constraints ----- #
@constraint(m, e2, x[2]+x[3] == 0.2)
@constraint(m, e3, x[4]+x[5] == 0.4)
@constraint(m, e4, x[6]+x[7] == 0.4)

# ----- Objective ----- #
@NLobjective(m, Min, ((15.3261663216011*x[2]+23.2043471859416*x[4]+6.69678129464404*x[6])*log(2.1055*x[2]+3.1878*x[4]+0.92*x[6])-2.46348749603266*x[2]-4.33155441248417*x[4]-0.626542690017204*x[6]+6.4661663216011*x[2]*log(x[2]/(2.1055*x[2]+3.1878*x[4]+0.92*x[6]))+12.2043471859416*x[4]*log(x[4]/(2.1055*x[2]+3.1878*x[4]+0.92*x[6]))+0.696781294644034*x[6]*log(x[6]/(2.1055*x[2]+3.1878*x[4]+0.92*x[6]))+9.86*x[2]*log(x[2]/(1.972*x[2]+2.4*x[4]+1.4*x[6]))+12*x[4]*log(x[4]/(1.972*x[2]+2.4*x[4]+1.4*x[6]))+7*x[6]*log(x[6]/(1.972*x[2]+2.4*x[4]+1.4*x[6]))+(1.972*x[2]+2.4*x[4]+1.4*x[6])*log(1.972*x[2]+2.4*x[4]+1.4*x[6])+1.972*x[2]*log(x[2]/(1.972*x[2]+0.283910843616504*x[4]+3.02002220174195*x[6]))+2.4*x[4]*log(x[4]/(1.45991339466884*x[2]+2.4*x[4]+0.415073537580851*x[6]))+1.4*x[6]*log(x[6]/(0.602183324335333*x[2]+0.115623371371275*x[4]+1.4*x[6]))+(15.3261663216011*x[3]+23.2043471859416*x[5]+6.69678129464404*x[7])*log(2.1055*x[3]+3.1878*x[5]+0.92*x[7])-2.46348749603266*x[3]-4.33155441248417*x[5]-0.626542690017204*x[7]+6.4661663216011*x[3]*log(x[3]/(2.1055*x[3]+3.1878*x[5]+0.92*x[7]))+12.2043471859416*x[5]*log(x[5]/(2.1055*x[3]+3.1878*x[5]+0.92*x[7]))+0.696781294644034*x[7]*log(x[7]/(2.1055*x[3]+3.1878*x[5]+0.92*x[7]))+9.86*x[3]*log(x[3]/(1.972*x[3]+2.4*x[5]+1.4*x[7]))+12*x[5]*log(x[5]/(1.972*x[3]+2.4*x[5]+1.4*x[7]))+7*x[7]*log(x[7]/(1.972*x[3]+2.4*x[5]+1.4*x[7]))+(1.972*x[3]+2.4*x[5]+1.4*x[7])*log(1.972*x[3]+2.4*x[5]+1.4*x[7])+1.972*x[3]*log(x[3]/(1.972*x[3]+0.283910843616504*x[5]+3.02002220174195*x[7]))+2.4*x[5]*log(x[5]/(1.45991339466884*x[3]+2.4*x[5]+0.415073537580851*x[7]))+1.4*x[7]*log(x[7]/(0.602183324335333*x[3]+0.115623371371275*x[5]+1.4*x[7]))-17.2981663216011*x[2]*log(x[2])-25.6043471859416*x[4]*log(x[4])-8.09678129464404*x[6]*log(x[6])-17.2981663216011*x[3]*log(x[3])-25.6043471859416*x[5]*log(x[5])-8.09678129464404*x[7]*log(x[7])))

JuMP.optimize!(m)
run_time = backend(m).optimizer.model.optimizer._run_time
println("run time: $run_time")


text:
# Pkg.add("Spirograph")
using Spirograph

spirograph(40,-11,18,linecolor=:red)
spirograph!(40,-11,16,linecolor=:blue)
spirograph!(40,-11,14,linecolor=:green)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


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

mutable struct IoK8sApiCoreV1NodeSystemInfo <: SwaggerModel
    architecture::Any # spec type: Union{ Nothing, String } # spec name: architecture
    bootID::Any # spec type: Union{ Nothing, String } # spec name: bootID
    containerRuntimeVersion::Any # spec type: Union{ Nothing, String } # spec name: containerRuntimeVersion
    kernelVersion::Any # spec type: Union{ Nothing, String } # spec name: kernelVersion
    kubeProxyVersion::Any # spec type: Union{ Nothing, String } # spec name: kubeProxyVersion
    kubeletVersion::Any # spec type: Union{ Nothing, String } # spec name: kubeletVersion
    machineID::Any # spec type: Union{ Nothing, String } # spec name: machineID
    operatingSystem::Any # spec type: Union{ Nothing, String } # spec name: operatingSystem
    osImage::Any # spec type: Union{ Nothing, String } # spec name: osImage
    systemUUID::Any # spec type: Union{ Nothing, String } # spec name: systemUUID

    function IoK8sApiCoreV1NodeSystemInfo(;architecture=nothing, bootID=nothing, containerRuntimeVersion=nothing, kernelVersion=nothing, kubeProxyVersion=nothing, kubeletVersion=nothing, machineID=nothing, operatingSystem=nothing, osImage=nothing, systemUUID=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("architecture"), architecture)
        setfield!(o, Symbol("architecture"), architecture)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("bootID"), bootID)
        setfield!(o, Symbol("bootID"), bootID)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("containerRuntimeVersion"), containerRuntimeVersion)
        setfield!(o, Symbol("containerRuntimeVersion"), containerRuntimeVersion)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("kernelVersion"), kernelVersion)
        setfield!(o, Symbol("kernelVersion"), kernelVersion)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("kubeProxyVersion"), kubeProxyVersion)
        setfield!(o, Symbol("kubeProxyVersion"), kubeProxyVersion)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("kubeletVersion"), kubeletVersion)
        setfield!(o, Symbol("kubeletVersion"), kubeletVersion)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("machineID"), machineID)
        setfield!(o, Symbol("machineID"), machineID)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("operatingSystem"), operatingSystem)
        setfield!(o, Symbol("operatingSystem"), operatingSystem)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("osImage"), osImage)
        setfield!(o, Symbol("osImage"), osImage)
        validate_property(IoK8sApiCoreV1NodeSystemInfo, Symbol("systemUUID"), systemUUID)
        setfield!(o, Symbol("systemUUID"), systemUUID)
        o
    end
end # type IoK8sApiCoreV1NodeSystemInfo

const _property_map_IoK8sApiCoreV1NodeSystemInfo = Dict{Symbol,Symbol}(Symbol("architecture")=>Symbol("architecture"), Symbol("bootID")=>Symbol("bootID"), Symbol("containerRuntimeVersion")=>Symbol("containerRuntimeVersion"), Symbol("kernelVersion")=>Symbol("kernelVersion"), Symbol("kubeProxyVersion")=>Symbol("kubeProxyVersion"), Symbol("kubeletVersion")=>Symbol("kubeletVersion"), Symbol("machineID")=>Symbol("machineID"), Symbol("operatingSystem")=>Symbol("operatingSystem"), Symbol("osImage")=>Symbol("osImage"), Symbol("systemUUID")=>Symbol("systemUUID"))
const _property_types_IoK8sApiCoreV1NodeSystemInfo = Dict{Symbol,String}(Symbol("architecture")=>"String", Symbol("bootID")=>"String", Symbol("containerRuntimeVersion")=>"String", Symbol("kernelVersion")=>"String", Symbol("kubeProxyVersion")=>"String", Symbol("kubeletVersion")=>"String", Symbol("machineID")=>"String", Symbol("operatingSystem")=>"String", Symbol("osImage")=>"String", Symbol("systemUUID")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1NodeSystemInfo }) = collect(keys(_property_map_IoK8sApiCoreV1NodeSystemInfo))
Swagger.property_type(::Type{ IoK8sApiCoreV1NodeSystemInfo }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1NodeSystemInfo[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1NodeSystemInfo }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1NodeSystemInfo[property_name]

function check_required(o::IoK8sApiCoreV1NodeSystemInfo)
    (getproperty(o, Symbol("architecture")) === nothing) && (return false)
    (getproperty(o, Symbol("bootID")) === nothing) && (return false)
    (getproperty(o, Symbol("containerRuntimeVersion")) === nothing) && (return false)
    (getproperty(o, Symbol("kernelVersion")) === nothing) && (return false)
    (getproperty(o, Symbol("kubeProxyVersion")) === nothing) && (return false)
    (getproperty(o, Symbol("kubeletVersion")) === nothing) && (return false)
    (getproperty(o, Symbol("machineID")) === nothing) && (return false)
    (getproperty(o, Symbol("operatingSystem")) === nothing) && (return false)
    (getproperty(o, Symbol("osImage")) === nothing) && (return false)
    (getproperty(o, Symbol("systemUUID")) === nothing) && (return false)
    true
end

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


text:
module MyExample

include("file1.jl")
export f

end


text:
######## GVar, a bundle that records gradient
"""
    GVar{T,GT} <: IWrapper{T}
    GVar(x)

Attach a gradient field to `x`.
The gradient type `GT` can be a `GVar` for storing second order gradients, the `AutoBcast` type for broadcasting or the same type of `x`.
"""
struct GVar{T,GT} <: IWrapper{T}
    x::T
    g::GT
    function GVar{T,GT}(x::T, g::GT) where {T,GT}
        new{T,GT}(x, g)
    end
    function GVar(x::T, g::T) where T<:Real
        new{T,T}(x, g)
    end
    function GVar{T,GT}(x::T) where {T, GT}
        new{T,GT}(x, zero(GT))
    end
    function GVar(x::T, g::GT) where {T,GT}
        new{T,GT}(x, g)
    end
end

# `GVar` and `~GVar` on composite types
@generated function GVar(x::Type{T}) where T
    :($(getfield(T.name.module, nameof(T))){$(GVar.(T.parameters)...)})
end
@generated function GVar(x::Type{T}, y::Type{T}) where T
    :($(getfield(T.name.module, nameof(T))){$(GVar.(T.parameters, T.parameters)...)})
end
@generated function (_::Type{Inv{GVar}})(x::Type{T}) where T
    :($(getfield(T.name.module, nameof(T))){$((~GVar).(T.parameters)...)})
end
# `GVar` and `~GVar` on composite vars
@generated function GVar(x::T) where T
    Expr(:new, GVar(T), [:(GVar(x.$NAME)) for NAME in fieldnames(T)]...)
end
@generated function GVar(x::T, g::T) where T
    Expr(:new, GVar(T, T), [:(GVar(x.$NAME, g.$NAME)) for NAME in fieldnames(T)]...)
end
@generated function (_::Type{Inv{GVar}})(x::T) where T
    Expr(:new, (~GVar)(T), [:((~GVar)(x.$NAME)) for NAME in fieldnames(T)]...)
end

for T in [:Real]
    ## differentiable elementary types
    @eval GVar(::Type{ET}) where ET<:$T = GVar{ET,ET}
    @eval GVar(::Type{ET}, ::Type{ET}) where ET<:$T = GVar{ET,ET}
    @eval (_::Type{Inv{GVar}})(::Type{GVar{ET,GT}}) where {ET<:$T,GT} = ET

    ## differentiable elementary vars
    @eval GVar(x::$T) = GVar(x, zero(x))
    @eval @inline function (_::Type{Inv{GVar}})(x::GVar{<:$T})
        @invcheck x.g zero(x.x)
        x.x
    end
end

for T in [:Integer, :Bool, :Function, :String, :Char, :Nothing]
    ## non-differentiable elementary types
    @eval GVar(::Type{ET}) where ET<:$T = ET
    @eval GVar(::Type{ET}, ::Type{ET}) where ET<:$T = GVar{ET,ET}
    @eval (_::Type{Inv{GVar}})(::Type{ET}) where ET<:$T = ET

    ## non-differentiable elementary vars
    @eval GVar(x::$T) = x
    @eval (_::Type{Inv{GVar}})(x::$T) = x
end

for T in [:Tuple, :AbstractArray]
    ## broadcastable elementary types
    @eval GVar(x::$T) = GVar.(x)
    @eval GVar(x::$T, y::$T) = GVar.(x, y)
    @eval (_::Type{Inv{GVar}})(x::$T) = (~GVar).(x)
end

# no gradient wrapper
GVar(x::NoGrad) = (~NoGrad)(x)

# define on complex numbers to fix ambiguity errors
GVar(x::Complex) = Complex(GVar(x.re), GVar(x.im))
GVar(x::Complex, y::Complex) = Complex(GVar(x.re, y.re), GVar(x.im, y.im))
(_::Type{Inv{GVar}})(x::Complex) = Complex((~GVar)(x.re), (~GVar)(x.im))

Base.copy(b::GVar) = GVar(b.x, copy(b.g))
Base.zero(x::GVar) = GVar(Base.zero(x.x), Base.zero(x.g))
Base.zero(::Type{<:GVar{T,GT}}) where {T,GT} = GVar(zero(T), zero(GT))
Base.one(x::GVar) = GVar(Base.one(x.x), Base.zero(x.g))
Base.one(::Type{<:GVar{T}}) where T = GVar(one(T))
Base.adjoint(b::GVar) = GVar(b.x', b.g')
Base.isapprox(x::GVar, y::GVar; kwargs...) = isapprox(x.x, y.x; kwargs...) && isapprox(x.g, y.g; kwargs...)

# define kernel and field views
"""
    grad(var)

Get the gradient field of `var`.
"""
@fieldview grad(gv::GVar) = gv.g
@fieldview value(gv::GVar) = gv.x
# TODO: fix the problem causing this patch, the field type can not change?!
chfield(x::GVar, ::typeof(value), xval::GVar) = GVar(xval, x.g)

@generated function grad(x::T) where T
    isprimitivetype(T) && throw("not supported type to obtain gradients: $T.")
    Expr(:new, typegrad(T), [:(grad(x.$NAME)) for NAME in fieldnames(T)]...)
end
typegrad(x) = x
@generated function typegrad(x::Type{T}) where T
    if isprimitivetype(T)
        T
    else
        :($(getfield(T.name.module, nameof(T))){$(typegrad.(T.parameters)...)})
    end
end
typegrad(::Type{GVar{ET,GT}}) where {ET,GT} = ET
grad(gv::T) where T<:Real = zero(T)
grad(gv::AbstractArray{T}) where T = grad.(gv)
grad(gv::Function) = 0
grad(gv::String) = ""
grad(t::Tuple) = grad.(t)
chfield(x::T, ::typeof(grad), g::T) where T = (@invcheck iszero(g) || g≈0; x)
chfield(x::GVar, ::typeof(grad), g::GVar) = GVar(x.x, g)
chfield(x::Complex{<:GVar}, ::typeof(grad), g::Complex) = Complex(GVar(value(x.re), g.re), GVar(value(x.im), g.im))

# NOTE: superwarning: check value only to make ancilla gradient descardable.
NiLangCore.deanc(x::GVar, val::GVar) = NiLangCore.deanc(value(x), value(val))
function deanc(x::T, val::T) where {T<:AbstractArray}
   x === val || deanc.(x, val)
end

# constructors and deconstructors
Base.iszero(x::GVar) = iszero(x.x)

## variable mapping
function (_::Type{Inv{GVar}})(x::GVar{<:GVar,<:GVar})
    Partial{:x}(x)
end

Base.show(io::IO, gv::GVar) = print(io, "GVar($(gv.x), $(gv.g))")
Base.show(io::IO, ::MIME"plain/text", gv::GVar) = Base.show(io, gv)

# used in log number iszero function.
Base.isfinite(x::GVar) = isfinite(x.x)
# interfaces

"""
    @nograd f(args...)

Mark `f(args...)` as having no gradients.
"""
macro nograd(ex)
    @smatch ex begin
        :($f($(args...))) => begin
            newargs = []
            for arg in args
                push!(newargs, @smatch arg begin
                    :($x::GVar) => :($x.x)
                    :($x::VecGVar) => :($x.x)
                    :($x::GVar{$tp}) => :($x.x)
                    _ => NiLangCore.get_argname(arg)
                end
                )
            end
            esc(quote
                @i function $f($(args...))
                    $f($(newargs...))
                end
            end)
        end
        _ => error("expect `f(args...)`, got $ex")
    end
end

# load data from stack
function loaddata(::Type{TG}, x::T) where {T,TG<:GVar{T}}
    TG(x)
end

function loaddata(::Type{T}, x::T) where T <: GVar
    x
end

function loaddata(::Type{AGT}, x::AT) where {T, GT, AT<:AbstractArray{T}, AGT<:AbstractArray{GVar{T,T}}}
    map(x->GVar(x, zero(x)), x)
end

# ULogarithmic
_content(x::ULogarithmic) = x.log
NiLang.AD.GVar(x::ULogarithmic) = exp(ULogarithmic, GVar(_content(x), zero(_content(x))))
(_::Type{Inv{GVar}})(x::ULogarithmic{GVar{TE}}) where TE = exp(ULogarithmic{TE}, (~GVar)(_content(x)))

Base.one(x::ULogarithmic{GVar{T,GT}}) where {T, GT} = one(ULogarithmic{GVar{T,GT}})
Base.one(::Type{ULogarithmic{GVar{T,GT}}}) where {T,GT} = exp(ULogarithmic, GVar(zero(T), zero(GT)))
Base.zero(x::ULogarithmic{GVar{T,GT}}) where {T,GT} =zero(ULogarithmic{GVar{T,GT}})
Base.zero(::Type{ULogarithmic{GVar{T,T}}}) where T = exp(ULogarithmic, GVar(zero(T), zero(T)))

function NiLang.loaddata(::Type{Array{<:ULogarithmic{GVar{T,T}}}}, data::Array{<:ULogarithmic{T}}) where {T}
    GVar.(data)
end

# the patch for dicts
function GVar(d::Dict)
    Dict([(k=>GVar(v)) for (k, v) in d])
end

function (_::Type{Inv{GVar}})(d::Dict)
    Dict([(k=>(~GVar)(v)) for (k, v) in d])
end

function grad(d::Dict)
    Dict([(k=>grad(v)) for (k, v) in d])
end


text:
using MathOptInterface, Test

const CleverDicts = MathOptInterface.Utilities.CleverDicts

# Note: `MyKey` is just for testing. You wouldn't want to use it in practice
# because the key type of the dictionary isn't a concrete type.
struct MyKey{X} end
CleverDicts.key_to_index(::MyKey{X}) where {X} = X
CleverDicts.index_to_key(::Type{MyKey}, index::Int) = MyKey{index}()

@testset "CleverDict" begin
    @testset "MyKey type" begin
        d = CleverDicts.CleverDict{MyKey, String}()
        key = CleverDicts.add_item(d, "first")
        @test key == MyKey{1}()
        @test d[MyKey{1}()] == "first"
    end

    @testset "get/set" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        key = CleverDicts.add_item(d, "first")
        @test key == MathOptInterface.VariableIndex(1)
        @test d[key] == "first"
        @test haskey(d, key) == true
        @test_throws KeyError d[MathOptInterface.VariableIndex(2)]
        delete!(d, key)
        @test_throws KeyError d[key]
        @test_throws KeyError d[key] = "key"
        @test haskey(d, key) == false
        key2 = CleverDicts.add_item(d, "second")
        @test key2 == MathOptInterface.VariableIndex(2)
        @test d[key2] == "second"
        @test d.vector === nothing
        @test d.dict !== nothing
        d[key2] = "third"
        @test d[key2] == "third"

        empty!(d)

        key = CleverDicts.add_item(d, "first")
        @test key == MathOptInterface.VariableIndex(1)
        @test d[key] == "first"
        d[key] = "zeroth"
        @test d[key] == "zeroth"
        @test haskey(d, key) == true
        @test_throws KeyError d[MathOptInterface.VariableIndex(2)]
        delete!(d, key)
        @test_throws KeyError d[key]
        @test_throws KeyError d[key] = "key"
        @test haskey(d, key) == false
        key2 = CleverDicts.add_item(d, "second")
        @test key2 == MathOptInterface.VariableIndex(2)
        @test d[key2] == "second"
        @test d.vector === nothing
        @test d.dict !== nothing
    end

    @testset "LinearIndex" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        key = CleverDicts.add_item(d, "first")
        @test d[CleverDicts.LinearIndex(1)] == "first"
        key2 = CleverDicts.add_item(d, "second")
        @test d[CleverDicts.LinearIndex(2)] == "second"
        @test length(d) == 2
        delete!(d, key)
        @test d.vector === nothing
        @test d[CleverDicts.LinearIndex(1)] == "second"
        @test_throws KeyError d[CleverDicts.LinearIndex(2)]
        @test length(d) == 1
        @test d.vector !== nothing
    end

    @testset "keys/values" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        key = CleverDicts.add_item(d, "first")
        key2 = CleverDicts.add_item(d, "second")
        @test collect(keys(d)) == [MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(2)]
        @test collect(values(d)) == ["first", "second"]
        delete!(d, key)
        key3 = CleverDicts.add_item(d, "third")
        @test collect(keys(d)) == [MathOptInterface.VariableIndex(2), MathOptInterface.VariableIndex(3)]
        @test collect(values(d)) == ["second", "third"]
    end

    @testset "iterate" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        key = CleverDicts.add_item(d, "first")
        key2 = CleverDicts.add_item(d, "second")
        my_keys = MathOptInterface.VariableIndex[]
        my_values = String[]
        for (k, v) in d
           push!(my_keys, k)
           push!(my_values, v)
        end
        @test my_keys == [MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(2)]
        @test my_values == ["first", "second"]
        delete!(d, key)
        key3 = CleverDicts.add_item(d, "third")
        my_keys = MathOptInterface.VariableIndex[]
        my_values = String[]
        for (k, v) in d
           push!(my_keys, k)
           push!(my_values, v)
        end
        @test my_keys == [MathOptInterface.VariableIndex(2), MathOptInterface.VariableIndex(3)]
        @test my_values == ["second", "third"]
    end

    @testset "iterate ii" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        key = CleverDicts.add_item(d, "first")
        key2 = CleverDicts.add_item(d, "second")
        my_keys = MathOptInterface.VariableIndex[]
        my_values = String[]
        for (k, v) in d
            push!(my_keys, k)
            push!(my_values, v)
        end
        @test my_keys == [MathOptInterface.VariableIndex(1), MathOptInterface.VariableIndex(2)]
        @test my_values == ["first", "second"]
        delete!(d, key)
        @test d[CleverDicts.LinearIndex(1)] == "second"
        key3 = CleverDicts.add_item(d, "third")
        my_keys = MathOptInterface.VariableIndex[]
        my_values = String[]
        for (k, v) in d
            push!(my_keys, k)
            push!(my_values, v)
        end
        @test my_keys == [MathOptInterface.VariableIndex(2), MathOptInterface.VariableIndex(3)]
        @test my_values == ["second", "third"]
    end

    @testset "iterate iii" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        y = 0
        for (k, v) in d
            y += 1
        end
        @test y == 0
    end

    @testset "haskey" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        @test !haskey(d, 1)
        k = CleverDicts.add_item(d, "a")
        @test haskey(d, k)
        j = CleverDicts.add_item(d, "b")
        @test haskey(d, j)
        delete!(d, k)
        @test !haskey(d, k)
        @test haskey(d, j)
    end

    @testset "haskey" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        @test isempty(d) == true
        k = CleverDicts.add_item(d, "a")
        @test isempty(d) == false
        delete!(d, k)
        @test isempty(d) == true
        j = CleverDicts.add_item(d, "b")
        @test isempty(d) == false
    end

    @testset "delete!" begin
        d = CleverDicts.CleverDict{MathOptInterface.VariableIndex, String}()
        @test length(d) == 0
        @test delete!(d, MathOptInterface.VariableIndex(0)) == nothing
        k1 = CleverDicts.add_item(d, "a")
        k2 = CleverDicts.add_item(d, "b")
        d[CleverDicts.LinearIndex(2)] == "b"
        delete!(d, k1)
        d[CleverDicts.LinearIndex(1)] == "b"
        k3 = CleverDicts.add_item(d, "c")
        @test d[k3] == "c"
        @test d[CleverDicts.LinearIndex(1)] == "b"
        @test d[CleverDicts.LinearIndex(2)] == "c"
    end
end


text:
# Julia wrapper for header: highs_c_api.h
# Automatically generated using Clang.jl

function Highs_call(numcol, numrow, numnz, colcost, collower, colupper, rowlower, rowupper, astart, aindex, avalue, colvalue, coldual, rowvalue, rowdual, colbasisstatus, rowbasisstatus, modelstatus)
    ccall((:Highs_call, libhighs), Cint, (Cint, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), numcol, numrow, numnz, colcost, collower, colupper, rowlower, rowupper, astart, aindex, avalue, colvalue, coldual, rowvalue, rowdual, colbasisstatus, rowbasisstatus, modelstatus)
end

function Highs_create()
    ccall((:Highs_create, libhighs), Ptr{Cvoid}, ())
end

function Highs_destroy(highs)
    ccall((:Highs_destroy, libhighs), Cvoid, (Ptr{Cvoid},), highs)
end

function Highs_readModel(highs, filename)
    ccall((:Highs_readModel, libhighs), Cint, (Ptr{Cvoid}, Cstring), highs, filename)
end

function Highs_writeModel(highs, filename)
    ccall((:Highs_writeModel, libhighs), Cint, (Ptr{Cvoid}, Cstring), highs, filename)
end

function Highs_clearModel(highs)
    ccall((:Highs_clearModel, libhighs), Cint, (Ptr{Cvoid},), highs)
end

function Highs_runQuiet(highs)
    ccall((:Highs_runQuiet, libhighs), Cint, (Ptr{Cvoid},), highs)
end

function Highs_run(highs)
    ccall((:Highs_run, libhighs), Cint, (Ptr{Cvoid},), highs)
end

function Highs_writeSolution(highs, filename)
    ccall((:Highs_writeSolution, libhighs), Cint, (Ptr{Cvoid}, Cstring), highs, filename)
end

function Highs_passLp(highs, numcol, numrow, numnz, colcost, collower, colupper, rowlower, rowupper, astart, aindex, avalue)
    ccall((:Highs_passLp, libhighs), Cint, (Ptr{Cvoid}, Cint, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}), highs, numcol, numrow, numnz, colcost, collower, colupper, rowlower, rowupper, astart, aindex, avalue)
end

function Highs_setHighsBoolOptionValue(highs, option, value)
    ccall((:Highs_setHighsBoolOptionValue, libhighs), Cint, (Ptr{Cvoid}, Cstring, Cint), highs, option, value)
end

function Highs_setHighsIntOptionValue(highs, option, value)
    ccall((:Highs_setHighsIntOptionValue, libhighs), Cint, (Ptr{Cvoid}, Cstring, Cint), highs, option, value)
end

function Highs_setHighsDoubleOptionValue(highs, option, value)
    ccall((:Highs_setHighsDoubleOptionValue, libhighs), Cint, (Ptr{Cvoid}, Cstring, Cdouble), highs, option, value)
end

function Highs_setHighsStringOptionValue(highs, option, value)
    ccall((:Highs_setHighsStringOptionValue, libhighs), Cint, (Ptr{Cvoid}, Cstring, Cstring), highs, option, value)
end

function Highs_setHighsOptionValue(highs, option, value)
    ccall((:Highs_setHighsOptionValue, libhighs), Cint, (Ptr{Cvoid}, Cstring, Cstring), highs, option, value)
end

function Highs_getHighsBoolOptionValue(highs, option, value)
    ccall((:Highs_getHighsBoolOptionValue, libhighs), Cint, (Ptr{Cvoid}, Cstring, Ptr{Cint}), highs, option, value)
end

function Highs_getHighsIntOptionValue(highs, option, value)
    ccall((:Highs_getHighsIntOptionValue, libhighs), Cint, (Ptr{Cvoid}, Cstring, Ptr{Cint}), highs, option, value)
end

function Highs_getHighsDoubleOptionValue(highs, option, value)
    ccall((:Highs_getHighsDoubleOptionValue, libhighs), Cint, (Ptr{Cvoid}, Cstring, Ptr{Cdouble}), highs, option, value)
end

function Highs_getHighsStringOptionValue(highs, option, value)
    ccall((:Highs_getHighsStringOptionValue, libhighs), Cint, (Ptr{Cvoid}, Cstring, Cstring), highs, option, value)
end

function Highs_resetHighsOptions(highs)
    ccall((:Highs_resetHighsOptions, libhighs), Cint, (Ptr{Cvoid},), highs)
end

function Highs_getHighsIntInfoValue(highs, info, value)
    ccall((:Highs_getHighsIntInfoValue, libhighs), Cint, (Ptr{Cvoid}, Cstring, Ptr{Cint}), highs, info, value)
end

function Highs_getHighsDoubleInfoValue(highs, info, value)
    ccall((:Highs_getHighsDoubleInfoValue, libhighs), Cint, (Ptr{Cvoid}, Cstring, Ptr{Cdouble}), highs, info, value)
end

function Highs_getSolution(highs, colvalue, coldual, rowvalue, rowdual)
    ccall((:Highs_getSolution, libhighs), Cvoid, (Ptr{Cvoid}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), highs, colvalue, coldual, rowvalue, rowdual)
end

function Highs_getBasis(highs, colstatus, rowstatus)
    ccall((:Highs_getBasis, libhighs), Cvoid, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}), highs, colstatus, rowstatus)
end

function Highs_getModelStatus(highs, scaled_model)
    ccall((:Highs_getModelStatus, libhighs), Cint, (Ptr{Cvoid}, Cint), highs, scaled_model)
end

function Highs_getObjectiveValue(highs)
    ccall((:Highs_getObjectiveValue, libhighs), Cdouble, (Ptr{Cvoid},), highs)
end

function Highs_getIterationCount(highs)
    ccall((:Highs_getIterationCount, libhighs), Cint, (Ptr{Cvoid},), highs)
end

function Highs_getSimplexIterationCount(highs)
    ccall((:Highs_getSimplexIterationCount, libhighs), Cint, (Ptr{Cvoid},), highs)
end

function Highs_getBasicVariables(highs, basic_variables)
    ccall((:Highs_getBasicVariables, libhighs), Cint, (Ptr{Cvoid}, Ptr{Cint}), highs, basic_variables)
end

function Highs_getBasisInverseRow(highs, row, row_vector, row_num_nz, row_indices)
    ccall((:Highs_getBasisInverseRow, libhighs), Cint, (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}), highs, row, row_vector, row_num_nz, row_indices)
end

function Highs_getBasisInverseCol(highs, col, col_vector, col_num_nz, col_indices)
    ccall((:Highs_getBasisInverseCol, libhighs), Cint, (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}), highs, col, col_vector, col_num_nz, col_indices)
end

function Highs_getBasisSolve(highs, rhs, solution_vector, solution_num_nz, solution_indices)
    ccall((:Highs_getBasisSolve, libhighs), Cint, (Ptr{Cvoid}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}), highs, rhs, solution_vector, solution_num_nz, solution_indices)
end

function Highs_getBasisTransposeSolve(highs, rhs, solution_vector, solution_nz, solution_indices)
    ccall((:Highs_getBasisTransposeSolve, libhighs), Cint, (Ptr{Cvoid}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}), highs, rhs, solution_vector, solution_nz, solution_indices)
end

function Highs_getReducedRow(highs, row, row_vector, row_num_nz, row_indices)
    ccall((:Highs_getReducedRow, libhighs), Cint, (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}), highs, row, row_vector, row_num_nz, row_indices)
end

function Highs_getReducedColumn(highs, col, col_vector, col_num_nz, col_indices)
    ccall((:Highs_getReducedColumn, libhighs), Cint, (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}), highs, col, col_vector, col_num_nz, col_indices)
end

function Highs_setBasis(highs, colstatus, rowstatus)
    ccall((:Highs_setBasis, libhighs), Cint, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}), highs, colstatus, rowstatus)
end

function Highs_setLogicalBasis(highs)
    ccall((:Highs_setLogicalBasis, libhighs), Cint, (Ptr{Cvoid},), highs)
end

function Highs_getHighsRunTime(highs)
    ccall((:Highs_getHighsRunTime, libhighs), Cdouble, (Ptr{Cvoid},), highs)
end

function Highs_addRow(highs, lower, upper, num_new_nz, indices, values)
    ccall((:Highs_addRow, libhighs), Cint, (Ptr{Cvoid}, Cdouble, Cdouble, Cint, Ptr{Cint}, Ptr{Cdouble}), highs, lower, upper, num_new_nz, indices, values)
end

function Highs_addRows(highs, num_new_row, lower, upper, num_new_nz, starts, indices, values)
    ccall((:Highs_addRows, libhighs), Cint, (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}), highs, num_new_row, lower, upper, num_new_nz, starts, indices, values)
end

function Highs_addCol(highs, cost, lower, upper, num_new_nz, indices, values)
    ccall((:Highs_addCol, libhighs), Cint, (Ptr{Cvoid}, Cdouble, Cdouble, Cdouble, Cint, Ptr{Cint}, Ptr{Cdouble}), highs, cost, lower, upper, num_new_nz, indices, values)
end

function Highs_addCols(highs, num_new_col, costs, lower, upper, num_new_nz, starts, indices, values)
    ccall((:Highs_addCols, libhighs), Cint, (Ptr{Cvoid}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}), highs, num_new_col, costs, lower, upper, num_new_nz, starts, indices, values)
end

function Highs_changeObjectiveSense(highs, sense)
    ccall((:Highs_changeObjectiveSense, libhighs), Cint, (Ptr{Cvoid}, Cint), highs, sense)
end

function Highs_changeColCost(highs, col, cost)
    ccall((:Highs_changeColCost, libhighs), Cint, (Ptr{Cvoid}, Cint, Cdouble), highs, col, cost)
end

function Highs_changeColsCostBySet(highs, num_set_entries, set, cost)
    ccall((:Highs_changeColsCostBySet, libhighs), Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}), highs, num_set_entries, set, cost)
end

function Highs_changeColsCostByMask(highs, mask, cost)
    ccall((:Highs_changeColsCostByMask, libhighs), Cint, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}), highs, mask, cost)
end

function Highs_changeColBounds(highs, col, lower, upper)
    ccall((:Highs_changeColBounds, libhighs), Cint, (Ptr{Cvoid}, Cint, Cdouble, Cdouble), highs, col, lower, upper)
end

function Highs_changeColsBoundsByRange(highs, from_col, to_col, lower, upper)
    ccall((:Highs_changeColsBoundsByRange, libhighs), Cint, (Ptr{Cvoid}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}), highs, from_col, to_col, lower, upper)
end

function Highs_changeColsBoundsBySet(highs, num_set_entries, set, lower, upper)
    ccall((:Highs_changeColsBoundsBySet, libhighs), Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}), highs, num_set_entries, set, lower, upper)
end

function Highs_changeColsBoundsByMask(highs, mask, lower, upper)
    ccall((:Highs_changeColsBoundsByMask, libhighs), Cint, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}), highs, mask, lower, upper)
end

function Highs_changeRowBounds(highs, row, lower, upper)
    ccall((:Highs_changeRowBounds, libhighs), Cint, (Ptr{Cvoid}, Cint, Cdouble, Cdouble), highs, row, lower, upper)
end

function Highs_changeRowsBoundsBySet(highs, num_set_entries, set, lower, upper)
    ccall((:Highs_changeRowsBoundsBySet, libhighs), Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}), highs, num_set_entries, set, lower, upper)
end

function Highs_changeRowsBoundsByMask(highs, mask, lower, upper)
    ccall((:Highs_changeRowsBoundsByMask, libhighs), Cint, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}), highs, mask, lower, upper)
end

function Highs_getObjectiveSense(highs, sense)
    ccall((:Highs_getObjectiveSense, libhighs), Cint, (Ptr{Cvoid}, Ptr{Cint}), highs, sense)
end

function Highs_getColsByRange(highs, from_col, to_col, num_col, costs, lower, upper, num_nz, matrix_start, matrix_index, matrix_value)
    ccall((:Highs_getColsByRange, libhighs), Cint, (Ptr{Cvoid}, Cint, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}), highs, from_col, to_col, num_col, costs, lower, upper, num_nz, matrix_start, matrix_index, matrix_value)
end

function Highs_getColsBySet(highs, num_set_entries, set, num_col, costs, lower, upper, num_nz, matrix_start, matrix_index, matrix_value)
    ccall((:Highs_getColsBySet, libhighs), Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}), highs, num_set_entries, set, num_col, costs, lower, upper, num_nz, matrix_start, matrix_index, matrix_value)
end

function Highs_getColsByMask(highs, mask, num_col, costs, lower, upper, num_nz, matrix_start, matrix_index, matrix_value)
    ccall((:Highs_getColsByMask, libhighs), Cint, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}), highs, mask, num_col, costs, lower, upper, num_nz, matrix_start, matrix_index, matrix_value)
end

function Highs_getRowsByRange(highs, from_row, to_row, num_row, lower, upper, num_nz, matrix_start, matrix_index, matrix_value)
    ccall((:Highs_getRowsByRange, libhighs), Cint, (Ptr{Cvoid}, Cint, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}), highs, from_row, to_row, num_row, lower, upper, num_nz, matrix_start, matrix_index, matrix_value)
end

function Highs_getRowsBySet(highs, num_set_entries, set, num_row, lower, upper, num_nz, matrix_start, matrix_index, matrix_value)
    ccall((:Highs_getRowsBySet, libhighs), Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}), highs, num_set_entries, set, num_row, lower, upper, num_nz, matrix_start, matrix_index, matrix_value)
end

function Highs_getRowsByMask(highs, mask, num_row, lower, upper, num_nz, matrix_start, matrix_index, matrix_value)
    ccall((:Highs_getRowsByMask, libhighs), Cint, (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}), highs, mask, num_row, lower, upper, num_nz, matrix_start, matrix_index, matrix_value)
end

function Highs_deleteColsByRange(highs, from_col, to_col)
    ccall((:Highs_deleteColsByRange, libhighs), Cint, (Ptr{Cvoid}, Cint, Cint), highs, from_col, to_col)
end

function Highs_deleteColsBySet(highs, num_set_entries, set)
    ccall((:Highs_deleteColsBySet, libhighs), Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}), highs, num_set_entries, set)
end

function Highs_deleteColsByMask(highs, mask)
    ccall((:Highs_deleteColsByMask, libhighs), Cint, (Ptr{Cvoid}, Ptr{Cint}), highs, mask)
end

function Highs_deleteRowsByRange(highs, from_row, to_row)
    ccall((:Highs_deleteRowsByRange, libhighs), Cint, (Ptr{Cvoid}, Cint, Cint), highs, from_row, to_row)
end

function Highs_deleteRowsBySet(highs, num_set_entries, set)
    ccall((:Highs_deleteRowsBySet, libhighs), Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}), highs, num_set_entries, set)
end

function Highs_deleteRowsByMask(highs, mask)
    ccall((:Highs_deleteRowsByMask, libhighs), Cint, (Ptr{Cvoid}, Ptr{Cint}), highs, mask)
end

function Highs_getHighsInfinity(highs)
    ccall((:Highs_getHighsInfinity, libhighs), Cdouble, (Ptr{Cvoid},), highs)
end

function Highs_getNumCols(highs)
    ccall((:Highs_getNumCols, libhighs), Cint, (Ptr{Cvoid},), highs)
end

function Highs_getNumRows(highs)
    ccall((:Highs_getNumRows, libhighs), Cint, (Ptr{Cvoid},), highs)
end

function Highs_getNumNz(highs)
    ccall((:Highs_getNumNz, libhighs), Cint, (Ptr{Cvoid},), highs)
end

function Highs_highsModelStatusToChar(highs, int_highs_model_status)
    ccall((:Highs_highsModelStatusToChar, libhighs), Cstring, (Ptr{Cvoid}, Cint), highs, int_highs_model_status)
end

function Highs_primalDualStatusToChar(highs, int_primal_dual_status)
    ccall((:Highs_primalDualStatusToChar, libhighs), Cstring, (Ptr{Cvoid}, Cint), highs, int_primal_dual_status)
end


text:
using CUDA, GLMakie

scene = scatter(rand(Point2f0, 10_000), show_axis=false)
screen = display(scene)

buffer = screen.renderlist[1][3][:position]
resource =
cuGraphicsGLRegisterBuffer(&resource, pbo, cudaGraphicsRegisterFlagsReadOnly)

CUDA.cuGraphicsMapResources


text:
api_path = joinpath(dirname(@__FILE__), "api")

files = [
    "keys.jl",
    "leaders.jl",
    "machines.jl",
    "members.jl",
    "stats.jl",
    "watch.jl",
]

for file in files
    include(joinpath(api_path, file))
end


text:
using Test
import CryptoGroups.Fields: F2PB, F2GNB, print_poly, red!, mul, mul_gnb, construct_integer_order_prime, FP, order, reducer, tobits
import CryptoGroups: @bin_str

f = bin"10011"
c = reverse(bin"1100101")
@test red!(c, reverse(f)) == bin"1111" # I already changed rem and order functions and that this is the result

### Additional test for multiplication

@test mul(bin"0100", bin"0100") == bin"00100000"

a = reverse(bin"1101")
b = reverse(bin"1001")

c = mul(a, b)

@test red!(c, reverse(f)) == bin"1111"

@test mul(reverse(bin"1101"), reverse(bin"1001"), reverse(bin"10011")) == bin"1111"
@test mul(reverse(bin"0010"), reverse(bin"0010"), reverse(bin"10011")) == reverse(bin"0100")
@test mul(reverse(bin"0010"), reverse(bin"0100"), reverse(bin"10011")) == reverse(bin"1000")

a = F2PB(f)(bin"1101")
b = F2PB(f)(bin"1001")

c = F2PB(f)(bin"1111")

@test a*b == c

### Testing generation of basis elements

α = F2PB(f)(bin"0010")

@test α^(order(α) + 1) == α

@test α^8 * α^8 == α^16
@test α * α * α * α * α == α^5
@test α^14 * α == α^15

@test α^3 * α^5 == α^8

@test inv(α) * α == one(α)

@test zero(F2GNB{5, 1}) == F2GNB{5, 1}(bin"00000")
@test one(F2GNB{5, 1}) == F2GNB{5, 1}(bin"11111")

a = F2GNB{4, 1}(bin"1101")
b = F2GNB{4, 1}(bin"1011")

@test a + b == F2GNB{4, 1}(bin"0110")

@test construct_integer_order_prime(3, 3*4 + 1) in [3, 9, 1]

### Now the last bit of getting multiplication right

@test mul_gnb(bin"1000", bin"1101", 3) == bin"0010"

###

a = F2GNB{4, 3}(bin"1000")
b = F2GNB{4, 3}(bin"1101")

c = F2GNB{4, 3}(bin"0010")

@test a * b == c

@test c^(order(c) + 1) == c
@test inv(c) * c == one(c)

a = FP{23}(17)
b = FP{23}(3)

@test a * b == FP{23}(5)

@test b^11 == FP{23}(1)

#### More diligent testing for polynomial basis

# function field2poly(z::BitVector)
#     poly = Int[]
#     for (i, b) in enumerate(z)
#         if b==true
#             push!(poly, i - 1)
#         end
#     end
#     return poly
# end

#field2poly(z) = field2poly(z.x)

# Arguments are written with little endianess in mind
# let

#     F = F2PB([163, 7, 6, 3, 0])
#     f = reducer(F)

#     x19 = F(BitVector(i==20 for i in 1:163))
#     @test findfirst(x -> x== 1, mul(x19.x, x19.x)) == 39

#     x38 = F(BitVector(i==39 for i in 1:163))
#     @test findfirst(x -> x==1, red!(x38.x, f)) == 39

#     @test field2poly(x19 * x19) == [38]

#     x144 = F(BitVector(i==145 for i in 1:163))
#     x145 = F(BitVector(i==146 for i in 1:163))
#     @test field2poly(x145 * x144) == [126, 129, 132, 133]

#     x160 = F(BitVector(i==161 for i in 1:163))
#     x159 = F(BitVector(i==160 for i in 1:163))

#     @test field2poly(x160 * x159) == [0, 3, 6, 7, 156, 159, 162]

#     _x319 = BitVector(i==320 for i in 1:163*2)
#     @test field2poly(red!(_x319, f)) == [0, 3, 6, 7, 156, 159, 162]

# end


text:
# Test Linear System Optimizaiton
using SwitchTimeOpt
using MathOptInterface
const MOI = MathOptInterface
using Ipopt

# Define test options
objtol = 1e-04
primaltol = 1e-03

# Define Solver options
maxiter = 50
maxtime = 100.0
verbose = 0
tolerance = 1e-06

#Define solver
solver = Ipopt.Optimizer()
MOI.set(solver, MOI.RawOptimizerAttribute("tol"), tolerance)
MOI.set(solver, MOI.RawOptimizerAttribute("print_level"), verbose)
MOI.set(solver, MOI.RawOptimizerAttribute("max_cpu_time"), maxtime)
MOI.set(solver, MOI.RawOptimizerAttribute("max_iter"), maxiter)

# Define Time Interval
t0 = 0.0; tf = 12.0

uvec = [repeat([0.0; 1.0], 4, 1); 0.0]  # Input Vector
uvec = copy(uvec')

# Cost Funcction Matrix
C = [1.0 0.0 -1.0;
      0.0 1.0 -1.0]
Q = C'*C

# Define Initial State
x0 = [0.5; 0.7; 1]

### Define System Dynamics
function nldyn(x::Array{Float64,1}, u::Array{Float64,1})
  n = length(x)
  f = zeros(n)
  f[1] = x[1] - x[1]*x[2] - 0.4*x[1]*u[1]
  f[2] = -x[2] + x[1]*x[2] - 0.2*x[2]*u[1]
  return f

end

function nldyn_deriv(x::Array{Float64,1}, u::Array{Float64,1})
  df = [1.0-x[2]-0.4*u[1]       -x[1]                   0;
        x[2]                     -1+x[1]-0.2*u[1]       0;
        0                        0                      0]
end

m = stoproblem(x0, nldyn, nldyn_deriv, uvec, ngrid = 300, t0=t0, tf=tf, Q=Q, solver=solver)

solve!(m)

# Test Optimal Solution
@testset "Test optimal switching times" begin
  tauopt = gettau(m)

  @test isapprox(tauopt[1], 2.4435555103155107, atol=primaltol)
  @test isapprox(tauopt[2], 4.122622390312436, atol=primaltol)
  @test isapprox(tauopt[3], 4.433725972515174, atol=primaltol)
  @test isapprox(tauopt[4], 4.682336152264894,  atol=primaltol)
  @test isapprox(tauopt[5], 5.199622121419033, atol=primaltol)
  @test isapprox(tauopt[6], 5.367438990179279, atol=primaltol)
  @test isapprox(tauopt[7], 6.370950931911677, atol=primaltol)
  @test isapprox(tauopt[8], 6.4660389494043295, atol=primaltol)
end

@testset "Test status and optimal objective value" begin
  # Test Optimal Value
  @test isapprox(getobjval(m), 1.3454355602054182, atol=objtol)
  @test string(getstat(m)) == "LOCALLY_SOLVED"
end


text:
module QuadraticFormsMGHyp

using LinearAlgebra
using Roots
using SpecialFunctions: besselk #, lgamma
using QuadGK: quadgk
using StatsFuns: normcdf, normpdf, logtwo
# work around https://github.com/JuliaMath/SpecialFunctions.jl/issues/186
# until https://github.com/JuliaDiff/ForwardDiff.jl/pull/419/ is merged
using Base.Math: libm
using ForwardDiff: Dual, value, partials, derivative
@inline lgamma(x::Float64) = ccall((:lgamma, libm), Float64, (Float64,), x)
@inline lgamma(x::Float32) = ccall((:lgammaf, libm), Float32, (Float32,), x)
@inline lgamma(d::Dual{T}) where {T} =
    Dual{T}(lgamma(value(d)), digamma(value(d)) * partials(d))
Base.@irrational rp 0.3183098861837906715 1 / big(pi)

export qfmgh

"""
    qfmgh(x::Union{AbstractVector{<:Real}, Real}, a0, a, A, C, mu, gam, lam, chi, psi; do_spa=false, order=2)

Survivor function P(L>x) = 1-F(x) and
tail conditional mean, E(L|L>x), of

L = a0 + a' * X + X' * A * X, where:

   X = mu + W * gam + sqrt(W) * C * Z,
   Z ~ N(0, I), and W ~ GIG(lam, chi, psi), i.e.,
   X is distributed as multivariate GHyp.

Keyword arguments:

    `do_spa`: whether to return the exact result or a saddlepoint approximation
    `order`: order of the saddlepoint approximation

 (c) 2020 S.A. Broda
"""
function qfmgh end

qfmgh(x::Real, args...; kwargs...) = getindex.(qfmgh([x], args...; kwargs...), 1)

function qfmgh(
    x::AbstractVector{<:Real},
    a0,
    a,
    A,
    C,
    mu,
    gam,
    lam,
    chi,
    psi;
    do_spa::Bool = false,
    order::Int = 2,
)
    if (do_spa && (lam>=0 || chi<= 0 ||  any(gam .!= 0 )))
        @warn "Saddlepoint approximation is inaccurate with these parameters."
    end

    CAC = Symmetric(C' * A * C)
    E = eigen(CAC)
    omega = E.values
    P = E.vectors
    muA = mu' * A
    CP = C * P
    gA = gam' * A
    c = a' * gam + 2 * muA * gam
    d = a' * CP + 2 * muA * CP
    e = 2 * gA * CP
    de = d .* e
    d2 = d .* d
    e2 = e .* e
    k = gA * gam
    kk = a0 + a' * mu + muA * mu
    LK2 = lklam(lam, chi, psi)
    qq = x .- kk
    ccdf = similar(float(x))
    pm = similar(float(x))
    lM, alpha2p, ldM0da1, alpha1p, lM0, lrhop = get_funcs(omega, de[:], e2[:], d2[:], c, k, LK2, lam, chi, psi)
    shat0 = zeros(Threads.nthreads())
    shat2 = zeros(Threads.nthreads())
    shat3 = zeros(Threads.nthreads())
    Threads.@threads for i = 1:length(qq)
        q = qq[i]
        if ~do_spa
            ccdf[i], _ = quadgk(s -> imag(exp(lM(1im * s, -1im * s * q)) / s), 0.0, Inf)
            ccdf[i] = exp(lM(0, 0)) / 2 + rp * ccdf[i]
            M2(s, t) = exp(lM(s, -q * s)) * alpha2p(s) + exp(ldM0da1(s, -q * s)) * alpha1p(s) + exp(lM0(s, -q * s)) * lrhop(s)
            pm[i], _ = quadgk(s -> imag(M2(1im * s, -1im * s * q) / s), 0.0, Inf)
            pm[i] = (M2(0, 0) / 2 + rp * pm[i]) / ccdf[i] + kk
        else
            ccdf[i], shat0[Threads.threadid()] = compute_spa(s -> 1, s -> lM(s, -q * s), order, shat0[Threads.threadid()])
            I1 = all(d.==0) ? 0. : compute_spa(alpha2p, s -> lM(s, -q * s), order, shat0[Threads.threadid()], false)[1]
            I2, shat2[Threads.threadid()] = all(gam.==0) ? (0., shat2[Threads.threadid()]) : compute_spa(alpha1p, s -> ldM0da1(s, -q * s), order, shat2[Threads.threadid()])
            I3, shat3[Threads.threadid()] = compute_spa(lrhop, s -> lM0(s, -q * s), order, shat3[Threads.threadid()])
            pm[i] = (I1 + I2 + I3) / ccdf[i] + kk
        end
    end
    return ccdf, pm
end

@inline function lklam(lam, chi, psi)
    if chi == 0
        if real(psi) < 0
            return Inf
        else
            return -lam * log(psi / 2) + lgamma(lam)
        end
    elseif psi == 0
        if real(chi) < 0
            return Inf
        else
            return lam * log(chi / 2) + lgamma(-lam)
        end
    elseif real(chi * psi) < 0
        return Inf
    else
        scp = sqrt(chi * psi)
        return logtwo + lam * log(chi / psi) / 2 + log(besselk(lam, scp))
    end
end

function compute_spa(g, h, order, shat=0., solve=true)
    h0 = h(0.0)
    g0 = g(0.0)
    hp = s -> derivative(h, s)
    hpp = s -> derivative(hp, s)

    if solve # otherwise, keep starting value
        shat = find_zero(hp, shat)
    end

    if abs(h(shat) - h0) < 1e-5
        @warn("Saddlepoint approximation is inaccurate; returning NaN.")
        spa, shat = NaN, 0.
    else
        what = sign(shat) * sqrt(-2 * (h(shat) - h0))
        H = hpp(shat)
        uhat = shat * sqrt(H)
        ghat = g(shat)
        spa = exp(h0) * (g0 * (1 - normcdf(what)) + normpdf(what) * (ghat / uhat - g0 / what))
        if order == 2
            gp = s -> derivative(g, s)
            gpp = s -> derivative(gp, s)
            hppp = s -> derivative(hpp, s)
            hpppp = s -> derivative(hppp, s)
            k3 = hppp(shat) / H^(3 // 2)
            k4 = hpppp(shat) / H^2
            T1 = ghat / uhat * ((k4 / 8 - 5 * k3^2 / 24) - uhat^-2 - k3 / (2 * uhat))
            T2 = shat * gp(shat) / uhat * (1 / uhat^2 + k3 / (2 * uhat))
            T3 = -gpp(shat) * shat^2 / (2 * uhat^3) + g0 * what^-3
            spa = spa + exp(h0) * normpdf(what) * (T1 + T2 + T3)
        end
    end
    return spa, shat
end

function get_funcs(omega, de, e2, d2, c, k, LK2, lam, chi, psi)
    @inline logTheta(s, t, j) = @fastmath (t1 = 0.; t2 = 0.; t3 = 0.; t4 = 0.;
                              @inbounds @simd ivdep for i = 1:length(omega)
                                  nu = 1 / (1 - 2 * omega[i] * s)
                                  t1 += d2[i] * nu
                                  t2 += e2[i] * nu
                                  t3 += de[i] * nu
                                  t4 += log(nu)
                              end;
                              lklam(lam + j,
                                        chi - 2 * (0.5 * s^2 * t1 + t),
                                        psi - 2 * (k * s + 0.5 * s^2 * t2),
                                        ) - LK2 + s * c + s^2 * t3 + 0.5 * t4,
                                    )
    lM(s, t) = logTheta(s, t, 0)
    lM0(s, t) = logTheta(s, t, 1)
    ldM0da1(s, t) = logTheta(s, t, 2)
    lrhop(s) = @fastmath (t=0.; @inbounds @simd ivdep for i = 1:length(omega)
                        nu = 1 / (1 - 2 * omega[i] * s)
                        t += 2 * s * de[i] * nu + 2 * s^2 * de[i] * omega[i] * nu +  omega[i] * nu
                      end;
                      t + c
                )
    alpha1p(s) = @fastmath (t=0.; @inbounds @simd ivdep for i = 1:length(omega)
                          nu = 1 / (1 - 2 * omega[i] * s)
                          t += s * e2[i] * nu + s^2 * e2[i] * omega[i] * nu^2
                        end;
                        t + k
                )

    alpha2p(s) = @fastmath (t=0.; @inbounds @simd ivdep for i = 1:length(omega)
                          nu = 1 / (1 - 2 * omega[i] * s)
                          t += s * d2[i] * nu + s^2 * d2[i] * omega[i] * nu^2
                        end;
                        t
                 )
    return lM, alpha2p, ldM0da1, alpha1p, lM0, lrhop
end
end


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

include("model_Access.jl")
include("model_ApplicationSecurityGroup.jl")
include("model_ApplicationSecurityGroupPropertiesFormat.jl")
include("model_AvailableProvidersList.jl")
include("model_AvailableProvidersListCity.jl")
include("model_AvailableProvidersListCountry.jl")
include("model_AvailableProvidersListParameters.jl")
include("model_AvailableProvidersListState.jl")
include("model_AzureReachabilityReport.jl")
include("model_AzureReachabilityReportItem.jl")
include("model_AzureReachabilityReportLatencyInfo.jl")
include("model_AzureReachabilityReportLocation.jl")
include("model_AzureReachabilityReportParameters.jl")
include("model_ConnectionMonitor.jl")
include("model_ConnectionMonitorDestination.jl")
include("model_ConnectionMonitorEndpoint.jl")
include("model_ConnectionMonitorEndpointFilter.jl")
include("model_ConnectionMonitorEndpointFilterItem.jl")
include("model_ConnectionMonitorHttpConfiguration.jl")
include("model_ConnectionMonitorIcmpConfiguration.jl")
include("model_ConnectionMonitorListResult.jl")
include("model_ConnectionMonitorOutput.jl")
include("model_ConnectionMonitorParameters.jl")
include("model_ConnectionMonitorQueryResult.jl")
include("model_ConnectionMonitorResult.jl")
include("model_ConnectionMonitorResultProperties.jl")
include("model_ConnectionMonitorSource.jl")
include("model_ConnectionMonitorSuccessThreshold.jl")
include("model_ConnectionMonitorTcpConfiguration.jl")
include("model_ConnectionMonitorTestConfiguration.jl")
include("model_ConnectionMonitorTestGroup.jl")
include("model_ConnectionMonitorWorkspaceSettings.jl")
include("model_ConnectionStateSnapshot.jl")
include("model_ConnectivityDestination.jl")
include("model_ConnectivityHop.jl")
include("model_ConnectivityHop2.jl")
include("model_ConnectivityInformation.jl")
include("model_ConnectivityIssue.jl")
include("model_ConnectivityIssue2.jl")
include("model_ConnectivityParameters.jl")
include("model_ConnectivitySource.jl")
include("model_Direction.jl")
include("model_EffectiveNetworkSecurityRule.jl")
include("model_ErrorDetails.jl")
include("model_ErrorResponse.jl")
include("model_ErrorResponse2.jl")
include("model_EvaluatedNetworkSecurityGroup.jl")
include("model_FlowLog.jl")
include("model_FlowLogFormatParameters.jl")
include("model_FlowLogInformation.jl")
include("model_FlowLogListResult.jl")
include("model_FlowLogProperties.jl")
include("model_FlowLogPropertiesFormat.jl")
include("model_FlowLogStatusParameters.jl")
include("model_HTTPConfiguration.jl")
include("model_HTTPHeader.jl")
include("model_HopLink.jl")
include("model_HopLink2.jl")
include("model_HopLinkProperties.jl")
include("model_HopLinkProperties2.jl")
include("model_IPVersion.jl")
include("model_IssueContext.jl")
include("model_IssueContext2.jl")
include("model_MatchedRule.jl")
include("model_NetworkConfigurationDiagnosticParameters.jl")
include("model_NetworkConfigurationDiagnosticProfile.jl")
include("model_NetworkConfigurationDiagnosticResponse.jl")
include("model_NetworkConfigurationDiagnosticResult.jl")
include("model_NetworkInterfaceAssociation.jl")
include("model_NetworkSecurityGroupResult.jl")
include("model_NetworkSecurityRulesEvaluationResult.jl")
include("model_NetworkWatcher.jl")
include("model_NetworkWatcherListResult.jl")
include("model_NetworkWatcherPropertiesFormat.jl")
include("model_NextHopParameters.jl")
include("model_NextHopResult.jl")
include("model_PacketCapture.jl")
include("model_PacketCaptureFilter.jl")
include("model_PacketCaptureListResult.jl")
include("model_PacketCaptureParameters.jl")
include("model_PacketCaptureQueryStatusResult.jl")
include("model_PacketCaptureResult.jl")
include("model_PacketCaptureResultProperties.jl")
include("model_PacketCaptureStorageLocation.jl")
include("model_ProtocolConfiguration.jl")
include("model_ProvisioningState.jl")
include("model_QueryTroubleshootingParameters.jl")
include("model_Resource.jl")
include("model_RetentionPolicyParameters.jl")
include("model_SecurityGroupNetworkInterface.jl")
include("model_SecurityGroupViewParameters.jl")
include("model_SecurityGroupViewResult.jl")
include("model_SecurityRule.jl")
include("model_SecurityRuleAccess.jl")
include("model_SecurityRuleAssociations.jl")
include("model_SecurityRuleDirection.jl")
include("model_SecurityRulePropertiesFormat.jl")
include("model_SubResource.jl")
include("model_SubnetAssociation.jl")
include("model_TagsObject.jl")
include("model_Topology.jl")
include("model_TopologyAssociation.jl")
include("model_TopologyParameters.jl")
include("model_TopologyResource.jl")
include("model_TrafficAnalyticsConfigurationProperties.jl")
include("model_TrafficAnalyticsProperties.jl")
include("model_TroubleshootingDetails.jl")
include("model_TroubleshootingParameters.jl")
include("model_TroubleshootingProperties.jl")
include("model_TroubleshootingRecommendedActions.jl")
include("model_TroubleshootingResult.jl")
include("model_VerificationIPFlowParameters.jl")
include("model_VerificationIPFlowResult.jl")