text:
using Test
using NCDatasets
using NCDatasets: nc_create, NC_NETCDF4, NC_CLOBBER, NC_NOWRITE, nc_def_dim, nc_def_compound, nc_insert_compound, nc_def_var, nc_put_var, nc_close, NC_INT, nc_unsafe_put_var, libnetcdf, check, ncType, nc_open, nc_inq_vartype, nc_inq_compound_nfields, nc_inq_compound_size, nc_inq_compound_name, nc_inq_compound_fieldoffset,nc_inq_compound_fieldndims,nc_inq_compound_fielddim_sizes, nc_inq_compound_fieldname, nc_inq_compound_fieldindex, nc_inq_compound_fieldtype, nc_inq_compound, nc_inq_varid, nc_get_var!

# mutable struct are not supported
# https://discourse.julialang.org/t/passing-an-array-of-structures-through-ccall/5194

struct s8
   i1::Cint
   i2::Cint
   f1::Cfloat
   d1::Cdouble
end

sz = (2,3)

data = Array{s8,2}(undef,sz)

for j = 1:sz[2]
    for i = 1:sz[1]
        data[i,j] = s8(i,j,1.2,2.3)
    end
end

T = eltype(data)
filename = tempname()
ncid = nc_create(filename, NC_NETCDF4|NC_CLOBBER)

x_dimid = nc_def_dim(ncid, "x", sz[1])
y_dimid = nc_def_dim(ncid, "y", sz[2])

dimids = [x_dimid, y_dimid]

typeid = nc_def_compound(ncid, sizeof(T), "sample_compound_type")

for i = 1:fieldcount(T)
    nctype = ncType[fieldtype(T,i)]
    nc_insert_compound(
        ncid, typeid,
        fieldname(T,i),
        fieldoffset(T,i), nctype)
end

varid = nc_def_var(ncid, "data", typeid, reverse(dimids))

nc_put_var(ncid, varid, data)
nc_close(ncid)

#=
run(`ncdump $filename`)
=#

module NCReconstructedTypes end

ncid = nc_open(filename,NC_NOWRITE)

varid = nc_inq_varid(ncid,"data")
xtype = nc_inq_vartype(ncid,varid)

@test nc_inq_compound_name(ncid,xtype) == "sample_compound_type"
@test nc_inq_compound_size(ncid,xtype) == sizeof(T)
@test nc_inq_compound_nfields(ncid,xtype) == fieldcount(T)

fieldid = 0

@test nc_inq_compound_fieldname(ncid,xtype,fieldid) == String(fieldnames(T)[fieldid+1])

_fieldname = String(fieldnames(T)[1])
@test nc_inq_compound_fieldindex(ncid,xtype,_fieldname) == 0

@test nc_inq_compound_fieldoffset(ncid,xtype,fieldid) == 0

nc_inq_compound_fieldoffset(ncid,xtype,1)

@test nc_inq_compound_fieldtype(ncid,xtype,fieldid) == ncType[fieldtype(T,fieldid+1)]

@test nc_inq_compound_fieldndims(ncid,xtype,fieldid) == 0

dim_sizes = nc_inq_compound_fielddim_sizes(ncid,xtype,fieldid)

type_name,type_size,type_nfields = nc_inq_compound(ncid,xtype)

@test type_name == "sample_compound_type"
@test type_size == sizeof(T)
@test type_nfields == fieldcount(T)

nfields = nc_inq_compound_nfields(ncid,xtype)
names = Symbol.(nc_inq_compound_fieldname.(ncid,xtype,0:(nfields-1)))

types = [NCDatasets.jlType[nc_inq_compound_fieldtype(ncid,xtype,fieldid)] for fieldid = 0:(nfields-1)]

# assume scalars
for fieldid = 0:nfields-1
    @assert nc_inq_compound_fieldndims(ncid,xtype,fieldid) == 0
end

using Random
reconname = Symbol(string(nc_inq_compound_name(ncid,xtype),"_",randstring(12)))

# from JLD2, MIT "Expat" License
# https://github.com/JuliaIO/JLD2.jl/blob/abb9e5920bbe956a4d9fd2f92550cd7ea0a715aa/src/data/reconstructing_datatypes.jl#L493

Core.eval(
    NCReconstructedTypes,
    Expr(:struct, false, reconname,
         Expr(:block,
              Any[ Expr(Symbol("::"), names[i], types[i]) for i = 1:length(types) ]...,
              # suppress default constructors, plus a bogus `new()` call to make sure
              # ninitialized is zero.
              Expr(:if, false, Expr(:call, :new))
              )))

T2 = getfield(NCReconstructedTypes, reconname)

data2 = Array{T2,2}(undef,sz...)

nc_get_var!(ncid, varid, data2)

data2[1,1] == data[1,1]

for fn = fieldnames(eltype(data2))
    @test getproperty.(data,fn) == getproperty.(data2,fn)
end


text:
using Roots
using DifferentialEquations
using Printf
include("leap_frog_definitions.jl")
#equations used for describing motion and conserved quantities

#condition for when to stop integrating if we hit the PSS enough times
condition1(u,t,integrator)= u[5]>max_hit
#affect used when counter hits max
affect1!(integrator) = terminate!(integrator)

#condition for when we hit PSS, X=0
function condition2(u,t,integrator)
u[1]
end

#affect when we hit X=0, update counter
function affect2!(integrator)
integrator.u[5]=integrator.u[5]+1
end

#construct callback for when counter hits and terminates integration
cb1 = DiscreteCallback(condition1,affect1!)

#callback for when we hit PSS

# affect2, nothing  ===> upcrossing    neg->pos
# nothing, affect2  ===> downcrossing  pos ->neg
# affect2           ===> event happens at both
#save_postion in list of points only saves after the even rather than before and/or after
cb2 = ContinuousCallback(condition2,affect2!,nothing)

#creates callback set, put ContinuousCallback first since the update will affect DiscreteCallback
cb=CallbackSet(cb2,cb1)

function PSS_function(Q,P, Energy,  t_end,    max_hit)
    #setting tolerance for determing if vortices move off to infinity
    tol_dist=1e-3

    #converting energy into the hamiltonian used to find the Eq_of_M
    H=(2*Energy)^2

    #for a given Q,P,H with X=0
    Y=Yfind(Q,P,H)

    if ~isempty(Y)
        u0=zeros(5)
        u0[1]=0 #X
        u0[2]=P #P
        u0[3]=Q #Q
        u0[4]=Y #Y
        u0[5]=0 # counter variable

        #constructor for ODE
        prob = ODEProblem(Eq_of_M,u0,(0., t_end))

        #solve ode , save_everystep=false is important to prevent sol to include all points, not just event points
        sol=solve(prob, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb,save_start=true,save_end=true,save_everystep=false)
        # sol=solve(prob, RK4(),maxiters=1e20, reltol=1e-6,abstol=1e-8,callback=cb,save_start=true,save_end=true,save_everystep=false)

        #find final value to check if energy was conserved
        uf=zeros(5)
        uf[1]=sol[1,end] #X
        uf[2]=sol[2,end] #P
        uf[3]=sol[3,end] #Q
        uf[4]=sol[4,end] #Y
        dH=abs(H_test(uf)-H)

        #output Q, P and dH
        return sol[:,2:end-1][3,:],sol[:,2:end-1][2,:], dH
    else
    #need to return 3 values, dH=1 flags that there is no Y
        return 0,0,1
    end
end

ext: jl


text:
# Autogenerated wrapper script for qwt_jll for x86_64-w64-mingw32
export qwt

using Qt_jll
using GLEW_jll
using GLFW_jll
using GLU_jll
JLLWrappers.@generate_wrapper_header("qwt")
JLLWrappers.@declare_library_product(qwt, "qwt.dll")
function __init__()
    JLLWrappers.@generate_init_header(Qt_jll, GLEW_jll, GLFW_jll, GLU_jll)
    JLLWrappers.@init_library_product(
        qwt,
        "bin\\qwt.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    eltypes = (N0f8, N0f16, Float32, Float64)        # eltypes of parametric colors
    pctypes = (Gray, RGB, AGray, GrayA, ARGB, RGBA)  # parametric colors
    cctypes = (Gray24, AGray32, RGB24, ARGB32)       # non-parametric colors
    dims  = (1, 2, 3, 4)

    for T in eltypes
        @assert precompile(clamp01, (T,))
        @assert precompile(clamp01nan, (T,))
        @assert precompile(scaleminmax, (T, T))
        @assert precompile(scalesigned, (T,))
        @assert precompile(scalesigned, (T,T,T))
        for C in pctypes
            @assert precompile(clamp01, (C{T},))
            @assert precompile(clamp01nan, (C{T},))
            @assert precompile(colorsigned, (C{T},C{T}))
        end
    end
    for C in cctypes
        @assert precompile(clamp01, (C,))
        @assert precompile(clamp01nan, (C,))
        @assert precompile(colorsigned, (C,C))
end
    for n in dims
        for T in eltypes
            @assert precompile(rawview, (Array{T,n},))
        end
        @assert precompile(normedview, (Array{UInt8,n},))
        @assert precompile(normedview, (Type{N0f8}, Array{UInt8,n}))
        @assert precompile(normedview, (Type{N0f16}, Array{UInt16,n}))
    end
    for T in eltypes
        for n in dims
            for C in pctypes
                @assert precompile(colorview, (Type{C}, Array{T,n}))
                @assert precompile(colorview, (Type{C{T}}, Array{T,n}))
                if T<:FixedPoint
                    R = FixedPointNumbers.rawtype(T)
                    RA = Base.ReinterpretArray{T,n,R,Array{R,n}}
                    precompile(colorview, (Type{C}, RA))
                    precompile(colorview, (Type{C{T}}, RA))
                    precompile(reinterpretc, (Type{C{T}}, RA))
                end
                @assert precompile(channelview, (Array{C{T},n},))
            end
            for C in cctypes
                @assert precompile(colorview, (Type{C}, Array{T,n}))
                if T<:FixedPoint
                    R = FixedPointNumbers.rawtype(T)
                    RA = Base.ReinterpretArray{T,n,R,Array{R,n}}
                    precompile(colorview, (Type{C}, RA))
                end
                @assert precompile(channelview, (Array{C,n},))
            end
        end
    end
end


text:
export BilinearLens

@doc doc"""

    BilinearLens(ϕ)

BilinearLens is a lensing operator that computes lensing with bilinear
interpolation. The action of the operator, as well as its adjoint, inverse,
inverse-adjoint, and gradient w.r.t. ϕ can all be computed. The log-determinant
of the operation is non-zero and can't be computed.

Internally, BilinearLens forms a sparse matrix with the interpolation weights,
which can be applied and adjoint-ed extremely fast (e.g. at least an order of
magnitude faster than LenseFlow). Inverse and inverse-adjoint lensing is
somewhat slower as it is implemented with several steps of the [preconditioned
generalized minimal residual](https://en.wikipedia.org/wiki/Generalized_minimal_residual_method)
algorithm, taking anti-lensing as the preconditioner.

!!! warning

    Due to [this bug](https://github.com/JuliaLang/PackageCompiler.jl/issues/379)
    in PackageCompiler, currently you have to run `using SparseArrays` by hand
    in your Julia session before `BilinearLens` is available.

"""
mutable struct BilinearLens{Φ,S} <: ImplicitOp{Basis,Spin,Pix}
    ϕ :: Φ
    sparse_repr :: S
    anti_lensing_sparse_repr :: Union{S, Nothing}
end

function BilinearLens(ϕ::FlatS0)

    # if ϕ == 0 then just return identity operator
    if norm(ϕ) == 0
        return BilinearLens(ϕ,I,I)
    end

    @unpack Nx,Ny,Nside,Δx,T = fieldinfo(ϕ)

    # the (i,j)-th pixel is deflected to (ĩs[i],j̃s[j])
    j̃s,ĩs = getindex.((∇*ϕ)./Δx, :Ix)
    ĩs .=  ĩs  .+ (1:Ny)
    j̃s .= (j̃s' .+ (1:Nx))'

    # sub2ind converts a 2D index to 1D index, including wrapping at edges
    indexwrapi(i) = mod(i - 1, Ny) + 1
    indexwrapj(i) = mod(i - 1, Nx) + 1
    sub2ind(i,j) = Base._sub2ind((Ny,Nx),indexwrapi(i),indexwrapj(j))

    # compute the 4 non-zero entries in L[I,:] (ie the Ith row of the sparse
    # lensing representation, L) and add these to the sparse constructor
    # matrices, M, and V, accordingly. this function is split off so it can be
    # called directly or used as a CUDA kernel
    function compute_row!(I, ĩ, j̃, M, V)

        # (i,j) indices of the 4 nearest neighbors
        left,right = floor(Int,ĩ) .+ (0, 1)
        top,bottom = floor(Int,j̃) .+ (0, 1)

        # 1-D indices of the 4 nearest neighbors
        M[4I-3:4I] .= @SVector[sub2ind(left,top), sub2ind(right,top), sub2ind(left,bottom), sub2ind(right,bottom)]

        # weights of these neighbors in the bilinear interpolation
        Δx⁻, Δx⁺ = ((left,right) .- ĩ)
        Δy⁻, Δy⁺ = ((top,bottom) .- j̃)
        A = @SMatrix[
            1 Δx⁻ Δy⁻ Δx⁻*Δy⁻;
            1 Δx⁺ Δy⁻ Δx⁺*Δy⁻;
            1 Δx⁻ Δy⁺ Δx⁻*Δy⁺;
            1 Δx⁺ Δy⁺ Δx⁺*Δy⁺
        ]
        V[4I-3:4I] .= inv(A)[1,:]

    end

    # a surprisingly large fraction of the computation for large Nside, so memoize it:
    @memoize getK(Nx,Ny) = Int32.((4:4*Nx*Ny+3) .÷ 4)

    # CPU
    function compute_sparse_repr(is_gpu_backed::Val{false})
        K = Vector{Int32}(getK(Nx,Ny))
        M = similar(K)
        V = similar(K,T)
        for I in 1:length(ĩs)
            compute_row!(I, ĩs[I], j̃s[I], M, V)
        end
        sparse(K,M,V,Nx*Ny,Nx*Ny)
    end

    # GPU
    function compute_sparse_repr(is_gpu_backed::Val{true})
        K = CuVector{Cint}(getK(Nx,Ny))
        M = similar(K)
        V = similar(K,T)
        cuda(ĩs, j̃s, M, V; threads=256) do ĩs, j̃s, M, V
            index = threadIdx().x
            stride = blockDim().x
            for I in index:stride:length(ĩs)
                compute_row!(I, ĩs[I], j̃s[I], M, V)
            end
        end
        # remove once CuSparseMatrixCOO makes it into official CUDA.jl:
        if !Base.isdefined(CUSPARSE,:CuSparseMatrixCOO)
            error("To use BilinearLens on GPU, run `using Pkg; pkg\"add https://github.com/marius311/CUDA.jl#coo\"` and restart Julia.")
        end
        switch2csr(CUSPARSE.CuSparseMatrixCOO{T}(K,M,V,(Nx*Ny,Nx*Ny)))
    end

    BilinearLens(ϕ, compute_sparse_repr(Val(is_gpu_backed(ϕ))), nothing)

end

# lazily computing the sparse representation for anti-lensing

function get_anti_lensing_sparse_repr!(Lϕ::BilinearLens)
    if Lϕ.anti_lensing_sparse_repr == nothing
        Lϕ.anti_lensing_sparse_repr = BilinearLens(-Lϕ.ϕ).sparse_repr
    end
    Lϕ.anti_lensing_sparse_repr
end

getϕ(Lϕ::BilinearLens) = Lϕ.ϕ
(Lϕ::BilinearLens)(ϕ) = BilinearLens(ϕ)

# applying various forms of the operator

function *(Lϕ::BilinearLens, f::FlatS0{P}) where {N,D,P<:Flat{N,<:Any,<:Any,D}}
    Lϕ.sparse_repr===I && return f
    Łf = Ł(f)
    f̃ = similar(Łf)
    ds = (D == 1 ? ((),) : tuple.(1:D))
    for d in ds
        mul!(@views(f̃.Ix[:,:,d...][:]), Lϕ.sparse_repr, @views(Łf.Ix[:,:,d...][:]))
    end
    f̃
end

function *(Lϕ::Adjoint{<:Any,<:BilinearLens}, f::FlatS0{P}) where {N,D,P<:Flat{N,<:Any,<:Any,D}}
    parent(Lϕ).sparse_repr===I && return f
    Łf = Ł(f)
    f̃ = similar(Łf)
    ds = (D == 1 ? ((),) : tuple.(1:D))
    for d in ds
        mul!(@views(f̃.Ix[:,:,d...][:]), parent(Lϕ).sparse_repr', @views(Łf.Ix[:,:,d...][:]))
    end
    f̃
end

function \(Lϕ::BilinearLens, f̃::FlatS0{P}) where {N,D,P<:Flat{N,<:Any,<:Any,D}}
    Łf̃ = Ł(f̃)
    f = similar(Łf̃)
    ds = (D == 1 ? ((),) : tuple.(1:D))
    for d in ds
        @views(f.Ix[:,:,d...][:]) .= gmres(
            Lϕ.sparse_repr, @views(Łf̃.Ix[:,:,d...][:]),
            Pl = get_anti_lensing_sparse_repr!(Lϕ), maxiter = 5
        )
    end
    f
end

function \(Lϕ::Adjoint{<:Any,<:BilinearLens}, f̃::FlatS0{P}) where {N,D,P<:Flat{N,<:Any,<:Any,D}}
    Łf̃ = Ł(f̃)
    f = similar(Łf̃)
    ds = (D == 1 ? ((),) : tuple.(1:D))
    for d in ds
        @views(f.Ix[:,:,d...][:]) .= gmres(
            parent(Lϕ).sparse_repr', @views(Łf̃.Ix[:,:,d...][:]),
            Pl = get_anti_lensing_sparse_repr!(parent(Lϕ))', maxiter = 5
        )
    end
    f
end

# optimizations for BilinearLens(0ϕ)
\(Lϕ::BilinearLens{<:Any,<:UniformScaling}, f::FlatS0{P}) where {N,P<:Flat{N}} = f
*(Lϕ::BilinearLens{<:Any,<:UniformScaling}, f::FlatS0{P}) where {N,P<:Flat{N}} = f
\(Lϕ::Adjoint{<:Any,<:BilinearLens{<:Any,<:UniformScaling}}, f::FlatS0{P}) where {N,P<:Flat{N}} = f
*(Lϕ::Adjoint{<:Any,<:BilinearLens{<:Any,<:UniformScaling}}, f::FlatS0{P}) where {N,P<:Flat{N}} = f

for op in (:*, :\)
    @eval function ($op)(Lϕ::Union{BilinearLens, Adjoint{<:Any,<:BilinearLens}}, f::FieldTuple)
        Łf = Ł(f)
        F = typeof(Łf)
        F(map(f->($op)(Lϕ,f), Łf.fs))
    end
end

# gradients

@adjoint BilinearLens(ϕ) = BilinearLens(ϕ), Δ -> (Δ,)

@adjoint function *(Lϕ::BilinearLens, f::Field{B}) where {B}
    f̃ = Lϕ * f
    function back(Δ)
        (∇' * (Ref(tuple_adjoint(Ł(Δ))) .* Ł(∇*f̃))), B(Lϕ*Δ)
    end
    f̃, back
end

# gpu

adapt_structure(storage, Lϕ::BilinearLens) = BilinearLens(adapt(storage, fieldvalues(Lϕ))...)


text:
## high level function dealing with NAs and formulas
##--------------------------------------------------

function reduceReturns(da::DataVector, allowedNAs::Float64 = 0.1;
                   log::Bool = true, pct::Bool = false)
    ## reduce series of returns to single number
    ## if more than allowedNAs NAs, return NA
    nObs = size(da, 1)
    nNAs = sum(isna(da))

    if nNAs/nObs > allowedNAs
        return NA
    else # apply reduction to values only
        vals = da[!isna(da)]
        if log
            reduce = reduceLog(vals)
        elseif !log & pct
            reduce = reduceDiscrPercent(vals)
        else
            reduce = reduceDiscr(vals)
        end
    end
    return reduce[1]
end

## different reduction formulas
##-----------------------------

function reduceLog(da::DataVector)
    return sum(da, 1)
end

function reduceDiscr(da::DataVector)
    return prod(da .+ 1, 1) .- 1
end

function reduceDiscrPercent(da::DataVector)
    return (prod(da./100 .+ 1, 1) .- 1).*100
end

##################
## partitioning ##
##################

function aggregate(da::DataVector, freq::Int;
          allowedNAs::Float64 = 0.1,
          log::Bool = true, pct::Bool = false)
    ## aggregate most recent observations, skipping start
    nObs = size(da, 1)
    nAggrRets = div(nObs, freq)
    firstObs = mod(nObs, freq)
    aggrRets = DataArray(eltype(da), nAggrRets)
    for ii=1:nAggrRets
        startSubInterval = (ii - 1)*freq + 1 + firstObs
        endSubInterval = ii*freq + firstObs

        aggrRets[ii] =
            reduceReturns(da[startSubInterval:endSubInterval],
                          allowedNAs;
                          log = log, pct = pct)
    end
    return aggrRets
end

function aggregate(tn::Timenum, freq::Int;
          allowedNAs::Float64 = 0.1,
          log::Bool = true, pct::Bool = false)
    ## apply aggregation to each column and attach reduced time index

end

function sampleAggregationSubsamples(nObs::Int, freq::Int, nSample::Int)
    ## number of aggregates that observations will be reduced to
    nAggrRets = div(nObs, freq) - 1
    if nAggrRets < 1
        error("to less observations to aggregate at this frequency")
    end

    subsampleSize = nAggrRets*freq
    nPossibleStarts = freq + mod(nObs, freq) + 1

    if nPossibleStarts < nSample
        # take each possible start
        startIndices = [1:nPossibleStarts]
    else
        # draw random starting points without replacement
        startIndices = randperm(nPossibleStarts)[1:nSample]
    end
    endIndices = startIndices .+ subsampleSize - 1
    return (startIndices, endIndices)
end


text:
struct Rect{T <: Number}
    lxv::T
    lyv::T
    rxv::T
    ryv::T
    Rect{T}(lxv::T, lyv::T, rxv::T, ryv::T) where T <: Number =
        new(min(lxv, rxv), min(lyv, ryv), max(lxv, rxv), max(lyv, ryv))

end

function Rect{T}(m::Matrix{T}) where T <: Number
    @assert size(m) == (2, 2) "Invalid values"
    return Rect{T}(m[1, 1], m[2, 1], m[1, 2], m[2, 2])
end
Rect{T}(r::Rect) where T <: Number =
    Rect{T}(convert(T, r.lxv), convert(T, r.lyv),
            convert(T, r.rxv), convert(T, r.ryv))
Rect(m::Matrix{T}) where T <: Number = Rect{T}(m)

matrix(r::Rect) = [r.lxv r.rxv; r.lyv r.ryv]

function Rect(lx::Number, ly::Number, rx::Number, ry::Number)
    t = promote(lx, ly, rx, ry)
    return Rect{eltype(t)}(t...)
end

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

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) = (r.lxv, r.lyv)
ru(r) = (r.rxv, r.ryv)
lx(r) = r.lxv
ly(r) = r.lyv
rx(r) = r.rxv
ry(r) = r.ryv

xplot(r::Rect{T}) where T = T[lx(r), rx(r), rx(r), lx(r)]
yplot(r::Rect{T}) where T = T[lx(r), lx(r), rx(r), rx(r)]

coord(r::Rect, axis::Int) = axis == 1 ? (r.lxv, r.rxv) : (r.lyv, r.ryv)
x(r) = coord(r, 1)
y(r) = coord(r, 2)
c_lo(r, axis)   = axis == 1 ? r.lxv : r.lyv
c_hi(r, axis)   = axis == 1 ? r.rxv : r.ryv

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) = append!(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)

Base.union(r1::Rect, r2::Rect) =
    Rect(min(lx(r1), lx(r2)), min(ly(r1), ly(r2)),
         max(rx(r1), rx(r2)), max(ry(r1), ry(r2)))

Base.union(r1::Rect, r2::Rect, y::Rect...) = Base.union(r1, union(r2, y...))
Base.union(r::Rect) = r
function Base.union(rs::AbstractVector{Rect{T}}) where T
    lxv, lyv, rxv, ryv = typemax(T), typemax(T), typemin(T), typemin(T)
    for r in rs
        lxv > lx(r) && (lxv = lx(r))
        lyv > ly(r) && (lyv = ly(r))
        rxv < rx(r) && (rxv = rx(r))
        ryv < ry(r) && (ryv = ry(r))
    end
    return Rect(lxv, lyv, rxv, ryv)
end

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} =
    abs(r1.lxv - r2.lxv) <= pcTol(T) &&
    abs(r1.lyv - r2.lyv) <= pcTol(T) &&
    abs(r1.rxv - r2.rxv) <= pcTol(T) &&
    abs(r1.ryv - r2.ryv) <= pcTol(T)
Base.:(==)(r1::Rect, r2::Rect) = ==(promote(r1, r2)...)

inside(p::Tuple{T, T}, r::Rect{T}) where T <: Number =
    r.lxv <= p[1] <= r.rxv && r.lyv <= p[2] <= r.ryv

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)
    (inside((sx(l), sy(l)), r) || inside((ex(l), ey(l)), r)) &&
        return true
    for tl in lines(r)
        intersects(tl, l) && return true
    end
    return false
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 && c_lo(r1, axis) > c_lo(r2, axis)) ? (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.lxv)
    tol = isopen ? pcTol(T) : zero(T)
    rt = axis == 1 ? Rect(lx(r1) + tol, ly(r1), rx(r1) - tol, ry(r1)) :
                     Rect(lx(r1), ly(r1) + tol, rx(r1), ry(r1) - tol)
    r1, r2 = sortr(rt, r2, reverse=false, axis=axis)
    return c_lo(r1, axis) <= c_lo(r2, axis) <= c_hi(r1, axis)
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)
        l = axis == 1 ? Rect{T}(axl[1], ly(lr), axl[2], ry(lr)) :
                        Rect{T}(lx(lr), axl[1], rx(lr), axl[2])
        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 = matrix(r1), matrix(r2)
    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) ? matrix(r1) : matrix(r2)
        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 = matrix(tr1)
    m[saxis, 1] = c_hi(tr1, saxis)
    m[saxis, 2] + pcTol(T) > c_lo(tr2, saxis) &&
        return Line(lx(tr2), ly(tr2), lx(tr2), ry(tr2))
    m[saxis, 2] = c_lo(tr2, saxis)
    return Rect(m)
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
    minv1, maxv1 = coord(rect, dir)
    if orm.reverseMax != zero(T1)
        minv1, maxv1 = (orm.reverseMax - maxv1), (orm.reverseMax - minv1)
    end
    minv2, maxv2 = coord(rect, odir)
    tol = isopen ? pcTol(T1) : zero(T1)
    imv1 = intersect(orm.data, Interval(minv1 + tol, maxv1 - tol))
    retr = Vector{Rect{T1}}()
    retv = Vector{V}()
    nimv1 = iterate(imv1)
    while nimv1 !== nothing
        v1, iv1 = nimv1
        imv2 = intersect(v1[2], Interval(minv2 + tol, maxv2 - 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
    minv1, maxv1 = coord(rect, dir)
    if orm.reverseMax != zero(T1)
        minv1, maxv1 = (orm.reverseMax - maxv1), (orm.reverseMax - minv1)
    end
    imv = get!(orm.data, Interval(minv1, maxv1), IntervalTree{T1, V}())
    minv2, maxv2 = coord(rect, odir)
    itv2 = Interval(minv2, maxv2)
    tv = get!(imv, itv2, v)
    tv === v && return nothing
    tv = delete!(imv, itv2)
    imv[itv2] = 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
    minv1, maxv1 = coord(rect, dir)
    if orm.reverseMax != zero(T1)
        minv1, maxv1 = (orm.reverseMax - maxv1), (orm.reverseMax - minv1)
    end
    itv1 = Interval(minv1, maxv1)
    imv = get(orm.data, itv1, IntervalTree{T1, V}())
    isempty(imv) && return nothing
    minv2, maxv2 = coord(rect, odir)
    ret = delete!(imv, Interval(minv2, maxv2))
    isempty(imv) && delete!(orm.data, itv1)
    return ret === nothing ? ret : ret[2]
end

function line_xsection(rect::Rect{T1}, ls::Vector{Line{T2}},
                       flines::Function, flt::Function,
                       o::Tuple{Int, Int}) where {T1 <: Number, T2 <: Number}
    r = convert(Rect{T2}, rect)
    l = flines(r)
    xlines = Vector{Int}()
    minv = searchsortedfirst(ls, l[o[1]], lt=flt)
    maxv = searchsortedlast(ls,  l[o[2]], lt=flt)
    for k = minv:maxv
        !intersects(r, ls[k]) && continue
        push!(xlines, k)
    end
    return xlines
end

"""
    Given a Rectangle and a set of already sorted set of vertical lines ordered
    left to right, provides the indices that intersect the rectangle.
"""
vline_xsection(r::Rect{T1}, vls::Vector{Line{T2}}) where {T1 <: Number,
                                                          T2 <: Number} =
    line_xsection(r, vls, vlines, vert_asc, (1, 2))

"""
    Given a Rectangle and a set of a sorted set of horizontal lines ordered
    top to bottom, provides the indices that intersect the rectangle.
"""
hline_xsection(r::Rect{T1}, hls::Vector{Line{T2}}) where {T1 <: Number,
                                                          T2 <: Number} =
    line_xsection(r, hls, hlines, horiz_desc, (2, 1))


text:
@add_problem socp function socp_norm_2_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable(2, 1)
    A = [1 2; 2 1; 3 4]
    b = [2; 3; 4]
    p = minimize(norm2(A * x + b))
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 0.64888 atol=atol rtol=rtol
        @test evaluate(norm2(A * x + b)) ≈ 0.64888 atol=atol rtol=rtol
    end

    x = Variable(2, 1)
    A = [1 2; 2 1; 3 4]
    b = [2; 3; 4]
    lambda = 1
    p = minimize(norm2(A * x + b) + lambda * norm2(x), x >= 1)
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 14.9049 atol=atol rtol=rtol
        @test evaluate(norm2(A * x + b) + lambda * norm2(x)) ≈ 14.9049 atol=atol rtol=rtol
    end

    x = Variable(2)

    p = minimize(norm2([x[1] + 2x[2] + 2; 2x[1] + x[2] + 3; 3x[1]+4x[2] + 4]) + lambda * norm2(x), x >= 1)
    if test
        @test vexity(p) == ConvexVexity()
    end

    handle_problem!(p)
    if test
        @test p.optval ≈ 14.9049 atol=atol rtol=rtol
        @test evaluate(norm2(A * x + b) + lambda * norm2(x)) ≈ 14.9049 atol=atol rtol=rtol
    end

    x = Variable(2, 1)
    A = [1 2; 2 1; 3 4]
    b = [2; 3; 4]
    lambda = 1
    p = minimize(norm2(A * x + b) + lambda * norm_1(x), x >= 1)
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 15.4907 atol=atol rtol=rtol
        @test evaluate(norm2(A * x + b) + lambda * norm_1(x)) ≈ 15.4907 atol=atol rtol=rtol
    end
end

@add_problem socp function socp_frobenius_norm_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    m = Variable(4, 5)
    c = [m[3, 3] == 4, m >= 1]
    p = minimize(norm(vec(m), 2), c)
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ sqrt(35) atol=atol rtol=rtol
        @test evaluate(norm(vec(m), 2)) ≈ sqrt(35) atol=atol rtol=rtol
    end
end

@add_problem socp function socp_quad_over_lin_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable(3, 1)
    A = [2 -3 5; -2 9 -3; 5 -8 3]
    b = [-3; 9; 5]
    c = [3 2 4]
    d = -3
    p = minimize(quadoverlin(A*x + b, c*x + d))
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 17.7831 atol=atol rtol=rtol
        @test (evaluate(quadoverlin(A * x + b, c * x + d)))[1] ≈ 17.7831 atol=atol rtol=rtol
    end
end

@add_problem socp function socp_sum_squares_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable(2, 1)
    A = [1 2; 2 1; 3 4]
    b = [2; 3; 4]
    p = minimize(sumsquares(A*x + b))
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 0.42105 atol=atol rtol=rtol
        @test (evaluate(sumsquares(A * x + b)))[1] ≈ 0.42105 atol=atol rtol=rtol
    end
end

@add_problem socp function socp_square_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable(2, 1)
    A = [1 2; 2 1; 3 4]
    b = [2; 3; 4]
    p = minimize(sum(square(A*x + b)))
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 0.42105 atol=atol rtol=rtol
        @test evaluate(sum(square(A * x + b))) ≈ 0.42105 atol=atol rtol=rtol
    end

    x = Variable(2, 1)
    A = [1 2; 2 1; 3 4]
    b = [2; 3; 4]
    expr = A * x + b
    p = minimize(sum(dot(^)(expr,2))) # elementwise ^
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 0.42105 atol=atol rtol=rtol
        @test evaluate(sum(broadcast(^, expr, 2))) ≈ 0.42105 atol=atol rtol=rtol
    end

    p = minimize(sum(dot(*)(expr, expr))) # elementwise *
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 0.42105 atol=atol rtol=rtol
        @test evaluate(sum((dot(*))(expr, expr))) ≈ 0.42105 atol=atol rtol=rtol
    end
end

@add_problem socp function socp_inv_pos_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable(4)
    p = minimize(sum(invpos(x)), invpos(x) < 2, x > 1, x == 2, 2 == x)
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 2 atol=atol rtol=rtol
        @test evaluate(sum(invpos(x))) ≈ 2 atol=atol rtol=rtol
    end

    x = Variable(3)
    p = minimize(sum(dot(/)([3,6,9], x)), x<=3)
    handle_problem!(p)
    if test
        @test x.value ≈ fill(3.0, (3, 1)) atol=atol rtol=rtol
        @test p.optval ≈ 6 atol=atol rtol=rtol
        @test evaluate(sum((dot(/))([3, 6, 9], x))) ≈ 6 atol=atol rtol=rtol
    end

    x = Variable()
    p = minimize(sum([3,6,9]/x), x<=3)
    handle_problem!(p)
    if test
        @test x.value ≈ 3 atol=atol rtol=rtol
        @test p.optval ≈ 6 atol=atol rtol=rtol
        @test evaluate(sum([3, 6, 9] / x)) ≈ 6 atol=atol rtol=rtol
    end
end

@add_problem socp function socp_geo_mean_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable(2)
    y = Variable(2)
    p = minimize(geomean(x, y), x >= 1, y >= 2)
    # not DCP compliant
    if test
        @test vexity(p) == ConcaveVexity()
    end
    p = maximize(geomean(x, y), 1 < x, x < 2, y < 2)
    # Just gave it a vector as an objective, not okay
    if test
        @test_throws Exception handle_problem!(p)
    end

    p = maximize(sum(geomean(x, y)), 1 < x, x < 2, y < 2)
    handle_problem!(p)
    if test
        @test p.optval ≈ 4 atol=atol rtol=rtol
        @test evaluate(sum(geomean(x, y))) ≈ 4 atol=atol rtol=rtol
    end
end

@add_problem socp function socp_sqrt_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable()
    p = maximize(sqrt(x), 1 >= x)
end

@add_problem socp function socp_quad_form_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable(3, 1)
    A = [0.8608 0.3131 0.5458; 0.3131 0.8584 0.5836; 0.5458 0.5836 1.5422]
    p = minimize(quadform(x, A), [x >= 1])
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 6.1464 atol=atol rtol=rtol
        @test (evaluate(quadform(x, A)))[1] ≈ 6.1464 atol=atol rtol=rtol
    end

    x = Variable(3, 1)
    A = -1.0*[0.8608 0.3131 0.5458; 0.3131 0.8584 0.5836; 0.5458 0.5836 1.5422]
    c = [3 2 4]
    p = maximize(c*x , [quadform(x, A) >= -1])
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 3.7713 atol=atol rtol=rtol
        @test (evaluate(quadform(x, A)))[1] ≈ -1 atol=atol rtol=rtol
    end
end

@add_problem socp function socp_huber_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable(3)
    p = minimize(sum(huber(x, 1)), x >= 2)
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)
    if test
        @test p.optval ≈ 9 atol=atol rtol=rtol
        @test evaluate(sum(huber(x, 1))) ≈ 9 atol=atol rtol=rtol
    end
end

@add_problem socp function socp_rational_norm_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    A = [1 2 3; -1 2 3]
    b = A * ones(3)
    x = Variable(3)
    p = minimize(norm(x, 4.5), [A * x == b])
    if test
        @test vexity(p) == ConvexVexity()
    end
    # Solution is approximately x = [1, .93138, 1.04575]
    handle_problem!(p)
    if test
        @test p.optval ≈ 1.2717 atol=atol rtol=rtol
        @test evaluate(norm(x, 4.5)) ≈ 1.2717 atol=atol rtol=rtol
    end
end

@add_problem socp function socp_rational_norm_dual_norm(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    v = [0.463339, 0.0216084, -2.07914, 0.99581, 0.889391]
    x = Variable(5)
    q = 1.379;  # q norm constraint that generates many inequalities
    qs = q / (q - 1);  # Conjugate to q
    p = minimize(x' * v)
    p.constraints += (norm(x, q) <= 1)
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p) # Solution is -norm(v, q / (q - 1))
    if test
        @test p.optval ≈ -2.144087 atol=atol rtol=rtol
        @test sum(evaluate(x' * v)) ≈ -2.144087 atol=atol rtol=rtol
        @test evaluate(norm(x, q)) ≈ 1 atol=atol rtol=rtol
        @test sum(evaluate(x' * v)) ≈ -(sum(abs.(v) .^ qs) ^ (1 / qs)) atol=atol rtol=rtol
    end
end

@add_problem socp function socp_rational_norm_atom_sum(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    A = [-0.719255  -0.229089
        -1.33632   -1.37121
        0.703447  -1.4482]
    b = [-1.82041, -1.67516, -0.866884]
    q = 1.5
    xvar = Variable(2)
    p = minimize(.5 * sumsquares(xvar) + norm(A * xvar - b, q))
    if test
        @test vexity(p) == ConvexVexity()
    end
    handle_problem!(p)

    if test
        # Compute gradient, check it is zero(ish)
        x_opt = xvar.value
        margins = A * x_opt - b
        qs = q / (q - 1);  # Conjugate
        denom = sum(abs.(margins).^q)^(1/qs)
        g = x_opt + A' * (abs.(margins).^(q-1) .* sign.(margins)) / denom
        @test p.optval ≈ 1.7227 atol=atol rtol=rtol
        @test norm(g, 2) ^ 2 ≈ 0 atol=atol rtol=rtol
    end
end

@add_problem socp function socp_norm_consistent_with_Base_for_matrix_variables(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    A = randn(4, 4)
    x = Variable(4, 4)
    x.value = A
    # Matrix norm
    if test
        @test evaluate(opnorm(x)) ≈ opnorm(A) atol=atol rtol=rtol
        @test evaluate(opnorm(x, 1)) ≈ opnorm(A, 1) atol=atol rtol=rtol
        @test evaluate(opnorm(x, 2)) ≈ opnorm(A, 2) atol=atol rtol=rtol
        @test evaluate(opnorm(x, Inf)) ≈ opnorm(A, Inf) atol=atol rtol=rtol
    end
    # Vector norm
    # TODO: Once the deprecation for norm on matrices is removed, remove the `vec` calls
    if test
        @test evaluate(norm(vec(x), 1)) ≈ norm(vec(A), 1) atol=atol rtol=rtol
        @test evaluate(norm(vec(x), 2)) ≈ norm(vec(A), 2) atol=atol rtol=rtol
        @test evaluate(norm(vec(x), 7)) ≈ norm(vec(A), 7) atol=atol rtol=rtol
        @test evaluate(norm(vec(x), Inf)) ≈ norm(vec(A), Inf) atol=atol rtol=rtol
    end
end

@add_problem socp function socp_fix_and_free_addition(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable()
    y = Variable()

    p = minimize(x+y, x>=0, y>=0)
    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol=atol rtol=rtol
    end

    y.value = 4
    fix!(y)
    handle_problem!(p)
    if test
        @test p.optval ≈ 4 atol=atol rtol=rtol
    end

    free!(y)
    handle_problem!(p)
    if test
        @test p.optval ≈ 0 atol=atol rtol=rtol
    end
end

@add_problem socp function socp_fix_multiplication(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    a = [1,2,3,2,1]
    x = Variable(length(a))
    gamma = Variable(Positive())
    fix!(gamma, 0.7)

    p = minimize(norm(x-a) + gamma*norm(x[1:end-1] - x[2:end]))
    handle_problem!(p)
    if test
        o1 = p.optval
        # x should be very close to a
        @test o1 ≈ 0.7 * norm(a[1:end - 1] - a[2:end]) atol=atol rtol=rtol
    end
    # increase regularization
    fix!(gamma, 1.0)
    handle_problem!(p)

    if test
        o2 = p.optval
        # x should be very close to mean(a)
        @test o2 ≈ norm(a .- mean(a)) atol=atol rtol=rtol
    end

    if test
        @test o1 <= o2
    end
end

@add_problem socp function socp_dual_minimal_norm_solutions(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
    x = Variable(2)
    A = [1 2; 2 4];
    b = [3, 6];
    p = minimize(norm(x, 1), A*x==b)

    if test
        @test vexity(p) == ConvexVexity()
    end

    handle_problem!(p)
    if test
        @test p.optval ≈ 1.5 atol=atol rtol=rtol
        @test evaluate(x) ≈ [0, 1.5] atol=atol rtol=rtol
        @test evaluate(norm(x, 1)) ≈ p.optval atol=atol rtol=rtol
        @test dot(b, p.constraints[1].dual) ≈ p.optval atol=atol rtol=rtol
    end

    x = Variable(2)
    A = [1 2; 2 4];
    b = [3, 6];
    p = minimize(norm(x, 2), A*x==b)

    test && @test vexity(p) == ConvexVexity()

    handle_problem!(p)

    if test
        @test p.optval ≈ 3/sqrt(5) atol=atol rtol=rtol
        @test evaluate(x) ≈ [3/5, 6/5] atol=atol rtol=rtol
        @test evaluate(norm(x, 2)) ≈ p.optval atol=atol rtol=rtol
        @test dot(b, p.constraints[1].dual) ≈ p.optval atol=atol rtol=rtol
    end

    x = Variable(2)
    A = [1 2; 2 4];
    b = [3, 6];
    p = minimize(norm(x, Inf), A*x==b)

    test && @test vexity(p) == ConvexVexity()

    handle_problem!(p)

    if test
        @test p.optval ≈ 1.0 atol=atol rtol=rtol
        @test evaluate(x) ≈ [1, 1] atol=atol rtol=rtol
        @test evaluate(norm(x, Inf)) ≈ p.optval atol=atol rtol=rtol
        @test dot(b, p.constraints[1].dual) ≈ p.optval atol=atol rtol=rtol
    end
end


text:
# Functions for handling effect formulae

"Representation of an effect as additions, deletions, and assignments."
mutable struct Diff
    add::Vector{Term} # List of additions
    del::Vector{Term} # List of deletions
    ops::Dict{Term,Any} # Dictionary of assignment operations
end

Diff() = Diff(Term[], Term[], Dict{Term,Any}())

"Returns true if the diff contains a matching term."
function contains_term(diff::Diff, term::Term)
    return (any(has_subterm(d, term) for d in diff.add) ||
            any(has_subterm(d, term) for d in diff.del) ||
            any(has_subterm(d, term) for d in keys(diff.ops)))
end

"Combine state differences, modifying the first diff in place."
function combine!(d::Diff, diffs::Diff...)
    d1 = diffs[1]
    append!(d.add, d1.add)
    append!(d.del, d1.del)
    merge!(d.ops, d1.ops)
    if (length(diffs) > 1) merge!(d, diffs[2:end]) end
    return d
end

"Combine state differences, returning a fresh diff."
function combine(diffs::Diff...)
    return combine!(Diff(), diffs...)
end

"Convert effect formula to a state difference (additions, deletions, etc.)"
function effect_diff(effect::Term, state::Union{State,Nothing}=nothing,
                     domain::Union{Domain,Nothing}=nothing)
    assign_ops = Dict{Symbol,Function}(
        :assign => (x, y) -> y, :increase => +, :decrease => -,
        Symbol("scale-up") => *, Symbol("scale-down") => /)
    diff = Diff()
    if effect.name == :and
        for eff in get_args(effect)
            combine!(diff, effect_diff(eff, state, domain))
        end
    elseif effect.name == :when
        cond, eff = effect.args[1], effect.args[2]
        if isnothing(state) || satisfy([cond], state, domain)[1] == true
            # Return eff if cond is satisfied
            combine!(diff, effect_diff(eff, state, domain))
        end
    elseif effect.name == :forall
        cond, eff = effect.args[1], effect.args[2]
        if !isnothing(state)
            # Find objects matching cond and apply effects for each
            _, subst = satisfy([cond], state, domain; mode=:all)
            for s in subst
                combine!(diff, effect_diff(substitute(eff, s), state, domain))
            end
        else
            push!(diff.add, effect)
        end
    elseif effect.name == :probabilistic
        n_effs = Int(length(get_args(effect))/2)
        r, cum_prob = rand(), 0.0
        for i in 1:n_effs
            # Sample a random effect
            prob, eff = effect.args[2*i-1].name, effect.args[2*i]
            if cum_prob <= r < (cum_prob + prob)
                diff = combine!(diff, effect_diff(eff, state, domain))
            end
            cum_prob += prob
        end
    elseif effect.name in keys(assign_ops)
        term, val = effect.args[1], effect.args[2]
        val = isnothing(state) ? val : evaluate(val, state, domain)
        diff.ops[term] = (assign_ops[effect.name], val.name)
    elseif effect.name in [:not, :!]
        effect = effect.args[1]
        if !isnothing(state) # Evaluated all nested functions
            effect = eval_term(effect, Subst(), state.fluents)
        end
        push!(diff.del, effect)
    else
        if !isnothing(state) # Evaluated all nested functions
            effect = eval_term(effect, Subst(), state.fluents)
        end
        push!(diff.add, effect)
    end
    return diff
end

"Relax actions"
function compute_relaxed_actions(actions::Dict{Symbol,Action})
    relaxed_effects = Dict{Symbol, Vector{Term}}()
    relaxed_preconditions = Dict{Symbol, Vector{Term}}()
    for (key, act) in actions
        diff_effect = effect_diff(act.effect)
        diff_precond = precond_diff(act.precond)
        relaxed_effects[key] = diff_effect.add
        relaxed_preconditions[key] = diff_precond.add
    end
    return relaxed_effects, relaxed_preconditions
end

"Convert precondition formula to a state difference."
function precond_diff(precond::Term, state::Union{State,Nothing}=nothing)
    diff = Diff()
    # TODO: Handle disjunctions and numeric conditions
    if precond.name == :and
        for eff in get_args(precond)
            combine!(diff, precond_diff(eff, state))
        end
    elseif precond.name == :forall
        cond, eff = precond.args[1], precond.args[2]
        if !isnothing(state)
            # Find objects matching cond and apply effects for each
            _, subst = satisfy([cond], state; mode=:all)
            for s in subst
                combine!(diff, precond_diff(substitute(eff, s), state))
            end
        else
            push!(diff.add, precond)
        end
    elseif precond.name in [:not, :!]
        precond = precond.args[1]
        push!(diff.del, precond)
    else
        push!(diff.add, precond)
    end
    return diff
end

"Update a world state (in-place) with a state difference."
function update!(state::State, diff::Diff)
    filter!(c -> !(c in diff.del), state.facts)
    union!(state.facts, diff.add)
    for (term, (op, val)) in diff.ops
        if isa(term, Const)
            oldval = get(state.fluents, term.name, 0)
            state.fluents[term.name] = op(oldval, val)
        elseif isa(term, Compound)
            valdict = get!(state.fluents, term.name, Dict())
            args = Tuple(a.name for a in term.args)
            oldval = get(valdict, args, 0)
            valdict[args] = op(oldval, val)
        end
    end
    return state
end

"Update a world state with a state difference."
function update(state::State, diff::Diff)
    return update!(copy(state), diff)
end

"A (categorical) distribution over possible state differences."
mutable struct DiffDist
    probs::Vector{Float64}
    diffs::Vector{Diff}
end

DiffDist() = DiffDist([1.0], [Diff()])

"Returns Cartesian product of two distributions over state differences."
function product(d1::DiffDist, d2::DiffDist)
    probs, diffs = Float64[], Diff[]
    for (p1, diff1) in zip(d1.probs, d1.diffs)
        for (p2, diff2) in zip(d2.probs, d2.diffs)
            push!(probs, p1 * p2)
            push!(diffs, combine(diff1, diff2))
        end
    end
    return DiffDist(probs, diffs)
end

"Combine distributions over state differences via Cartesian product."
function combine(dists::DiffDist...)
    return foldl(product, dists; init=DiffDist())
end

"Return distribution over possible effects, given a state."
function effect_dist(effect::Term, state::State,
                     domain::Union{Domain,Nothing}=nothing)
    if effect.name == :and
        sub_dists = [effect_dist(arg, state, domain) for arg in effect.args]
        dist = combine(sub_dists...)
    elseif effect.name == :when
        cond, eff = effect.args[1], effect.args[2]
        if satisfy([cond], state, domain)[1] == true
            dist = effect_dist(eff, state, domain)
        end
    elseif effect.name == :forall
        cond, eff = effect.args[1], effect.args[2]
        # Find objects matching cond and apply effects for each
        _, subst = satisfy([cond], state, domain; mode=:all)
        sub_dists = [effect_dist(substitute(eff, s), state, domain)
                     for s in subst]
        dist = combine(sub_dists...)
    elseif effect.name == :probabilistic
        n_effs = Int(length(effect.args)/2)
        probs, diffs = Float64[], Diff[]
        for i in 1:n_effs
            prob, eff = effect.args[2*i-1].name, effect.args[2*i]
            sub_dist = effect_dist(eff, state, domain)
            for (p, diff) in sub_dist
                push!(probs, prob * p)
                push!(diffs, diff)
            end
        end
        dist = DiffDist(probs, diffs)
    else
        # Bottom out at deterministic additions, deletions, or assignments
        diff = effect_diff(effect, state, domain)
        dist = DiffDist([1.0], [diff])
    end
    return dist
end

"Return the effect of a null operation as either a difference or distribution."
function no_effect(as_dist::Bool=false)
    if as_dist
        return DiffDist()
    else
        return Diff()
    end
end

"Update a state with a distribution over state differences."
function update(state::State, dist::DiffDist)
    return [(p, update(state, d)) for (p, d) in zip(dist.probs, dist.diffs)]
end


text:
#=
depth_first_search:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-10-06

search node as far as possible

Tree definition
     1
    / \
   2   5
 / \  / \
3  4 6  7
=#

module DFS
    # node class
    mutable struct Node
        val
        left
        right
    end

    function search(root, target)
        # define empty stack
        node_stack = []

        # add root node to stack
        push!(node_stack, root)

        # iterate search until empty
        while length(node_stack) > 0
            # get node from tail of stack
            node = node_stack[end]
            pop!(node_stack)

            if node != nothing
                println("Search Node $(node.val)")

                # found target
                if node.val == target
                    return true
                end

                # add child nodes to stack
                push!(node_stack, node.left)
                push!(node_stack, node.right)
            end
        end
    end

    function main()
        # define tree
        root = Node(1, nothing, nothing)
        root.left = Node(2, nothing, nothing)
        root.right = Node(5, nothing, nothing)
        root.left.left = Node(3, nothing, nothing)
        root.left.right = Node(4, nothing, nothing)
        root.right.left = Node(6, nothing, nothing)
        root.right.right = Node(7, nothing, nothing)

        # target node
        target = 4

        # search target node
        result = search(root, target)
        if result == true
            println("Found target node $(target)")
        else
            println("Not found")
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .DFS
    DFS.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using LinearAlgebra
using Test

using SetProg, SetProg.Sets
using Polyhedra
using MultivariatePolynomials

using DynamicPolynomials
import MultivariateBases
const MB = MultivariateBases
const MonoBasis = MB.MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}

using JuMP
const MOIT = MOI.Test

function switched_system_test(optimizer, config::MOIT.Config,
                              variable::SetProg.AbstractVariable, γ,
                              feasible, objective_value, set_test, dual_test)
    model = _model(optimizer)

    @variable(model, ◯, variable)

    # See https://github.com/blegat/SwitchOnSafety.jl/blob/master/examples/LPJ17e43.ipynb
    A1 = [1  0
          1  0] / γ
    A2 = [0  1
          0 -1] / γ

    cref1 = @constraint(model, A1 * ◯ ⊆ ◯)
    cref2 = @constraint(model, A2 * ◯ ⊆ ◯)

	if objective_value !== nothing
		@objective(model, MOI.MAX_SENSE, L1_heuristic(volume(◯), ones(2)))
	end

    SetProg.optimize!(model)
	if objective_value === nothing
		@test JuMP.objective_sense(model) == MOI.FEASIBILITY_SENSE
	else
		@test JuMP.objective_sense(model) == MOI.MIN_SENSE
	end
    if feasible
        @test JuMP.termination_status(model) == MOI.OPTIMAL
        @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
		if objective_value !== nothing
			@test JuMP.objective_value(model) ≈ objective_value atol=config.atol rtol=config.rtol
		end
        return set_test(JuMP.value(◯))
    else
        #@test JuMP.termination_status(model) == MOI.INFEASIBLE # Mosek returns SLOW_PROGRESS
        @test JuMP.dual_status(model) == MOI.INFEASIBILITY_CERTIFICATE
        return dual_test(cref1, cref2)
    end
end

function feasible_switched_system_ell_test(optimizer, config, ε=1e-3)
    Q = [1.0 0.0
         0.0 1.0]
    superset = SetProg.Sets.Ellipsoid(Symmetric(Q))
    switched_system_test(
        optimizer, config,
        Ellipsoid(symmetric=true, superset=superset),
        √2 + ε, true, 4/3,
        ◯ -> begin
            @test ◯ isa Sets.Ellipsoid{Float64}
            @test ◯.Q ≈ Q atol=config.atol rtol=config.rtol
        end,
        (cref1, cref2) -> begin end)
end
function infeasible_switched_system_ell_test(optimizer, config, ε=1e-3)
    Q = [1.0 0.0
         0.0 1.0]
    superset = SetProg.Sets.Ellipsoid(Symmetric(Q))
    switched_system_test(
        optimizer, config,
        Ellipsoid(symmetric=true, superset=superset),
        √2 - ε, false, NaN,
        ◯ -> begin end,
        (cref1, cref2) -> begin
            function ispsd(q)
                Q = [q[1] q[2]; q[2] q[3]]
                return all(eigvals(Q) .≥ -config.atol)
            end
            @test ispsd(JuMP.dual(cref1))
            @test ispsd(JuMP.dual(cref2))
        end)
end

function superset(x, d)
    q = SetProg.SumOfSquares.GramMatrix(SetProg.SumOfSquares.SOSDecomposition(x.^d))
    return SetProg.Sets.PolySet(2d, q)
end

function feasible_switched_system_quad_test(optimizer, config, ε=1e-3)
    @polyvar x[1:2]
    switched_system_test(
        optimizer, config,
        PolySet(symmetric=true, degree=2, superset=superset(x, 1)),
        √2 + ε, true, 8/3,
        ◯ -> begin
            @test ◯ isa Sets.PolySet{Float64, MonoBasis}
            @test polynomial(◯.p) ≈ x[1]^2 + x[2]^2 atol=config.atol rtol=config.rtol
        end,
        (cref1, cref2) -> begin end)
end
function infeasible_switched_system_quad_test(optimizer, config, ε=1e-3)
    @polyvar x[1:2]
    switched_system_test(
        optimizer, config,
        PolySet(symmetric=true, degree=2, superset=superset(x, 1)),
        √2 - ε, false, NaN,
        ◯ -> begin end,
        (cref1, cref2) -> begin
            function ispsd(M)
                return all(eigvals(Matrix(M.Q)) .≥ -config.atol)
            end
            @test ispsd(SetProg.SumOfSquares.moment_matrix(cref1))
            @test ispsd(SetProg.SumOfSquares.moment_matrix(cref2))
        end)
end

function feasible_switched_system_quartic_test(optimizer, config, ε=1e-2)
    @polyvar x[1:2]
    switched_system_test(
        optimizer, config,
        PolySet(symmetric=true, degree=4, superset=superset(x, 2)),
        1.0 + ε, true, 10.001105454190741,
        ◯ -> begin
            @test ◯ isa Sets.PolySet{Float64, MonoBasis}
			α = 11.814054544955727
            @test polynomial(◯.p) ≈ (α+1) * x[1]^4 - 2α * x[1]^2*x[2]^2 + (α+1) * x[2]^4 atol=config.atol rtol=config.rtol
        end,
        (cref1, cref2) -> begin end)
end
function infeasible_switched_system_quartic_test(optimizer, config, ε=2e-1)
    @polyvar x[1:2]
    switched_system_test(
        optimizer, config,
        PolySet(symmetric=true, degree=4, superset=superset(x, 2)),
        1.0 - ε, false, nothing,
        ◯ -> begin end,
        (cref1, cref2) -> begin
            function ispsd(M)
                return all(eigvals(Matrix(M.Q)) .≥ -config.atol)
            end
            @test ispsd(SetProg.SumOfSquares.moment_matrix(cref1))
            @test ispsd(SetProg.SumOfSquares.moment_matrix(cref2))
        end)
end

const switched_system_tests = Dict("feasible_switched_system_ell"     => feasible_switched_system_ell_test,
                                   "infeasible_switched_system_ell"   => infeasible_switched_system_ell_test,
                                   "feasible_switched_system_quad"    => feasible_switched_system_quad_test,
                                   "infeasible_switched_system_quad"  => infeasible_switched_system_quad_test,
                                   "feasible_switched_system_quartic"    => feasible_switched_system_quartic_test,
                                   "infeasible_switched_system_quartic"  => infeasible_switched_system_quartic_test)

@test_suite switched_system


text:
module ArrayInterfaceStaticArraysCore

import StaticArraysCore, ArrayInterfaceCore, Adapt
using LinearAlgebra

ArrayInterfaceCore.ismutable(::Type{<:StaticArraysCore.StaticArray}) = false
ArrayInterfaceCore.ismutable(::Type{<:StaticArraysCore.MArray}) = true
ArrayInterfaceCore.ismutable(::Type{<:StaticArraysCore.SizedArray}) = true

ArrayInterfaceCore.can_setindex(::Type{<:StaticArraysCore.StaticArray}) = false
ArrayInterfaceCore.buffer(A::Union{StaticArraysCore.SArray,StaticArraysCore.MArray}) = getfield(A, :data)

function ArrayInterfaceCore.lu_instance(_A::StaticArraysCore.StaticMatrix{N,N}) where {N}
    lu(one(_A))
end

function ArrayInterfaceCore.restructure(x::StaticArraysCore.SArray{S,T,N}, y::StaticArraysCore.SArray) where {S,T,N}
    StaticArraysCore.SArray{S,T,N}(y)
end
ArrayInterfaceCore.restructure(x::StaticArraysCore.SArray{S}, y) where {S} = StaticArraysCore.SArray{S}(y)

Adapt.adapt_storage(::Type{<:StaticArraysCore.SArray{S}}, xs::Array) where {S} = SArray{S}(xs)

end