text:
module Handlers

export handle, gethandler, Handler, HandlerFunction, Router, register!

import ..Nothing, ..Cvoid, ..Val

using HTTP

"""
handle(handler::Handler, request) => Response

Function used to dispatch to a Handler. Called from the core HTTP.serve method with the
initial Handler passed to `HTTP.serve(handler=handler)`.
"""
function handle end

"""
Abstract type representing an object that knows how to "handle" a server request.

Types of handlers include `HandlerFunction` (a julia function of the form `f(request)` and
`Router` (which pattern matches request url paths to other specific `Handler` types).
"""
abstract type Handler end

"""
HandlerFunction(f::Function)

A Function-wrapper type that is a subtype of `Handler`. Takes a single Function as an argument.
The provided argument should be of the form `f(request) => Response`, i.e. it accepts a `Request` returns a `Response`.
"""
struct HandlerFunction{F <: Function} <: Handler
    func::F # func(req)
end

handle(h::HandlerFunction, req) = h.func(req)

"A default 404 Handler"
const FourOhFour = HandlerFunction(req -> HTTP.Response(404))

"""
Router(h::Handler)
Router(f::Function)
Router()

An `HTTP.Handler` type that supports mapping request url paths to other `HTTP.Handler` types.
Can accept a default `Handler` or `Function` that will be used in case no other handlers match; by
default, a 404 response handler is used.
Paths can be mapped to a handler via `HTTP.register!(r::Router, path, handler)`, see `?HTTP.register!` for more details.
"""
struct Router <: Handler
    segments::Dict{String, Val}
    sym::Symbol
    func::Function
    function Router(ff::Union{Handler, Function, Nothing}=nothing)
        sym = gensym()
        if ff == nothing
            f = @eval $sym(args...) = FourOhFour
        else
            f = ff isa Function ? HandlerFunction(ff) : ff
        end
        r = new(Dict{String, Val}(), sym, f)
        return r
    end
end

const SCHEMES = Dict{String, Val}("http" => Val(:http), "https" => Val(:https))
const EMPTYVAL = Val(())

"""
HTTP.register!(r::Router, url, handler)
HTTP.register!(r::Router, m::String, url, handler)

Function to map request urls matching `url` and an optional method `m` to another `handler::HTTP.Handler`.
URLs are registered one at a time, and multiple urls can map to the same handler.
The URL can be passed as a String or `HTTP.URI` object directly. Requests can be routed based on: method, scheme,
hostname, or path.
The following examples show how various urls will direct how a request is routed by a server:

- `"http://*"`: match all HTTP requests, regardless of path
- `"https://*"`: match all HTTPS requests, regardless of path
- `"google"`: regardless of scheme, match requests to the hostname "google"
- `"google/gmail"`: match requests to hostname "google", and path starting with "gmail"
- `"/gmail"`: regardless of scheme or host, match any request with a path starting with "gmail"
- `"/gmail/userId/*/inbox`: match any request matching the path pattern, "*" is used as a wildcard that matches any value between the two "/"
"""
register!(r::Router, url, handler) = register!(r, "", url, handler)

function register!(r::Router, method::String, url, handler)
    m = isempty(method) ? Any : typeof(Val(Symbol(method)))
    # get scheme, host, split path into strings & vals
    uri = url isa String ? HTTP.URI(url) : url
    s = uri.scheme
    sch = !isempty(s) ? typeof(get!(SCHEMES, s, Val(s))) : Any
    h = !isempty(uri.host) ? Val{Symbol(uri.host)} : Any
    hand = handler isa Function ? HandlerFunction(handler) : handler
    register!(r, m, sch, h, uri.path, hand)
end

function splitsegments(r::Router, h::Handler, segments)
    vals = Expr[]
    for s in segments
        if s == "*" #TODO: or variable, keep track of variable types and store in handler
            T = Any
        else
            v = Val(Symbol(s))
            r.segments[s] = v
            T = typeof(v)
        end
        push!(vals, Expr(:(::), T))
    end
    return vals
end

function register!(r::Router, method::DataType, scheme, host, path, handler)
    # save string => Val mappings in r.segments
    segments = map(String, split(path, '/'; keep=false))
    vals = splitsegments(r, handler, segments)
    # return a method to get dispatched to
    #TODO: detect whether defining this method will create ambiguity?
    @eval $(r.sym)(::$method, ::$scheme, ::$host, $(vals...), args...) = $handler
    return
end

function gethandler(r::Router, req)
    # get the url/path of the request
    m = Val(Symbol(req.method))
    # get scheme, host, split path into strings and get Vals
    uri = HTTP.URI(req.target)
    s = get(SCHEMES, uri.scheme, EMPTYVAL)
    h = Val(Symbol(uri.host))
    p = uri.path
    segments = split(p, '/'; keep=false)
    # dispatch to the most specific handler, given the path
    vals = (get(r.segments, s, EMPTYVAL) for s in segments)
    return r.func(m, s, h, vals...)
end

function handle(r::Router, req)
    handler = gethandler(r,req)
    # pass the request to the handler and return
    return handle(handler, req)
end

handle(hf::HandlerFunction, http::HTTP.Stream) = HTTP.handle_request(hf.func,http)

end # module


text:
using Base

# referred base types or modules
import Base:    Indices, CommonReduceResult, SmallSigned, SmallUnsigned,
                OneTo, Bottom, Generator, BitInteger_types, Broadcast,
                ReshapedArray, SparseArrays

# related to traits
import Base:    HasLength, HasShape, SizeUnknown, HasEltype, EltypeUnknown,
                IteratorSize

# they are used only
import Base:    scalarmin, scalarmax, uniontypes, arraylen, sizeof, safe_tail, tail,
                indices1, getindex, setindex!, setindex_shape_check, one, zero

# those are re-defined and overwritten
import Base:    sum, sum!, prod, prod!, minimum, minimum!, maximum, maximum!,
                all, all!, any, any!,
                convert, copy, copy!, size, length, _length,
                reducedim, mapreduce_impl,
                findmin, findmin!, findmax, findmax!

import Base.==

# macros
import Base:    @_inline_meta, @_noinline_meta, @_propagate_inbounds_meta

ext: jl


text:
module PhaseSlopeIndex

using Statistics: mean, std
using FFTW: fft
using Einsum

# Exports
#---
export data2psi

"""
    int(x) = floor(Int64, x)
"""
int(x) = floor(Int64, x)

"""
    dropmean(X, d) = dropdims(mean(X, dims=d), dims=d)
"""
dropmean(X, d) = dropdims(mean(X; dims=d); dims=d)

"""
    squeeze(X::AbstractArray)

removing singleton dimensions
"""
function squeeze(X::AbstractArray)
    keepd = ()
    for (i, d) in enumerate(size(X))
        if d != 1
            keepd = (keepd..., d)
        end
    end
    return reshape(X, keepd)
end

"""
    detrend!(data, n)

(in place) Linear detrend of signals along first axis

### Arguments

  - `data::AbstractArray`: N-dim array where signal is in column-major order
  - `n::Integer`: n=0 subtracts mean from data, n=1 removes linear trend

**Note**: shape of data must be (signal length, ...)
"""
function detrend!(data::AbstractArray, n::Integer)
    original_shape = size(data)
    nsamp = size(data, 1)  # number of samples

    A = Array{Float64}([ones(nsamp) Array(1:nsamp)])

    data = reshape(data, (nsamp, :))  # reshaping data
    if n == 0
        data .-= mean(data; dims=1)
    elseif n == 1
        data .-= A * (A \ data)
    end
    return reshape(data, original_shape)
end

"""
    window = hanning_fun(N)

Hanning window similar to MATLAB implementation
"""
function hanning_fun(N::Integer)
    x = [range(0.0, 1.0; length=N + 2);]
    window = map(x -> 0.5 * (1 - cospi(2 * x)), x)
    return window[2:(end - 1)]
end

"""
    data2para(data, seglen, segshift, eplen, freqlist, method, nboot)

Extracts and builds a named tuple of parameters.

### Arguments

  - `data::AbstractArray`: NxM array for N data points in M channels.
  - `seglen::Integer`: segment length (determines the frequency resolution).
  - `segshift::Integer`: number of bins by which neighboring segments are shifted.
    e.g. segshift=seglen/2 makes overlapping segments
  - `eplen::Integer`: length of epochs
  - `freqlist::AbstractArray`: 2D Array where each column is a frequency band
  - `method::String`: standard deviation (error) estimation method
  - `subave::Bool`: if true, subtract average from CS segments (for continuous data, subave = false)
  - `nboot::Integer`: number of bootstrap resamplings

### Returns

  - `parameters::NamedTuple`: a named tuple of parameters
"""
function data2para(
    data::AbstractArray,
    seglen::Integer,
    segshift::Integer,
    eplen::Integer,
    freqlist::AbstractArray,
    method::String,
    subave::Bool,
    nboot::Integer,
)

    # We would like to avoid transpose and copying the data!
    if size(data, 1) < size(data, 2)
        @info "data is transposed to (#samples, #channels)"
        data = reshape(data, size(data, 2), size(data, 1))
    end

    # number of samples per channel and number of channels
    nsamples, nchannels = size(data)

    # jackknife check
    if eplen == 0 && method == "jackknife"
        @warn "Epoch length = 0 => Bootstrap method will be used."
        method = "bootstrap"
    end

    # if eplen = nsamples: continuous recording
    if eplen == 0
        eplen = nsamples
    end
    nep = int(nsamples / eplen)  # number of epochs

    if nep == 1 && subave == true
        subave = false
        @warn "subave is set to false for continuous data"
    end

    if segshift == 0
        segshift = int(seglen / 2)
    end
    nseg = int((eplen - seglen) / segshift) + 1

    # size(freqlist) = (freqs, nfbands)
    if length(freqlist) == 0
        freqlist = reshape(Array(1:(int(seglen / 2) + 1)), (:, 1))
    elseif ndims(freqlist) == 1
        freqlist = reshape(freqlist, :, 1)
    end
    if size(freqlist, 1) < size(freqlist, 2)
        @info "freqlist is transposed to (#freq, #nfbands)"
        freqlist = freqlist'
    end
    maxfreq = maximum(freqlist)  # TODO: maximum(freqlist, dims=1)
    nfbands = size(freqlist, 2)  # multiple frequency bands

    if method == "jackknife"
        n_resample = nep
    elseif method == "bootstrap"
        n_resample = nboot
    else
        throw(
            method *
            " is not a supported error estimation method!" *
            " Please choose from [jackknife, bootstrap] methods.",
        )
    end

    # we use named tuples to book the parameters
    parameters = (
        data=data,
        nsamples=nsamples,
        nchan=nchannels,
        eplen=eplen,
        nep=nep,
        method=method,
        subave=subave,
        segshift=segshift,
        nseg=nseg,
        freqlist=freqlist,
        maxfreq=maxfreq,
        nfbands=nfbands,
        n_resample=n_resample,
    )

    return parameters
end

"""
    make_eposeg(data, seglen, nep, nseg, nchan, segshift)

Partitioning data into epochs and segments

### Arguments

  - `data::AbstractArray`: NxM array for N data points in M channels.
  - `seglen::Integer`: segment length
  - `nep::Integer`: number of epochs
  - `nseg::Integer`: number of segments
  - `nchan::Integer`: number of channels
  - `segshift::Integer`: number of bins by which neighboring segments are shifted.

### Returns

  - `epseg::AbstractArray`: partitioned data into shape (seglen, nep, nseg, nchan)

**Note**: returned `epseg` may have more data entries than input data.
"""
function make_eposeg(
    data::AbstractArray,
    seglen::Integer,
    eplen::Integer,
    nep::Integer,
    nseg::Integer,
    nchan::Integer,
    segshift::Integer,
)::AbstractArray

    # preallocation
    epseg = Array{Float64}(undef, seglen, nep, nseg, nchan)

    for (i, e) in zip(1:nep, 1:eplen:(nep * eplen))
        for (j, s) in zip(1:nseg, 1:segshift:(nseg * seglen))
            @views epseg[:, i, j, :] = data[e:(e + eplen - 1), :][s:(s + seglen - 1), :]
        end
    end

    return epseg
end

"""
    cs2ps(cs)

Cross Spectra to Phase Slope

### Arguments

  - `cs::AbstractArray`: Cross Spectral array with size (seglen, :, nchan, nchan)

### Returns

  - phase slope index (AbstractArray) as Eq. 3 of the reference paper

**Note**: frequency resolution is assumed to be the resolution of freq. band!
"""
function cs2ps(cs::AbstractArray)

    # size(cs) = (seglen, nchan, nchan)
    # we don't use df but we slice the fband

    # complex coherency
    @einsum coh[f, i, j] := cs[f, i, j] / sqrt(cs[f, i, i] * cs[f, j, j])

    # phase slope (Eq. 3)
    @views imag.(sum(conj(coh[1:(end - 1), :, :]) .* coh[2:end, :, :]; dims=1))
end

"""
    data2ps(data)

Epoched segmented data to Cross Spectra

### Arguments

  - `data::AbstractArray`: Segmented data of shape (maxfreq, nep, nseg, nchan)

### Return

  - `cs::AbstractArray{Complex}`: Cross Spectral as eq. 2 of reference paper.
"""
function data2cs(data::AbstractArray)
    # cs: cross-spectral matrix

    # Eq. 2 size(csepseg) = (maxfreq, nep, nseg, nchan, nchan)
    @einsum cs[f, e, s, i, j] := data[f, e, s, i] * conj(data[f, e, s, j])
end

"""
    _cs_ = cs2cs_(data, cs, fband, nep, segave, subave, method)

preparing Cross Spectra for Phase Slope by segment averaging and subtraction

### Arguments

  - `data::AbstractArray`: Fourier-transformed detrended epoched segmented data.
  - `cs::AbstractArray`: Cross Spectra of data
  - `fband::AbstractArray`: 1D array of frequency range for PSI calculation
  - `nep::Integer`: number of epochs
  - `segave::Bool`: if true, averages across CS segments
  - `subave::Bool`: if true, subtract average across CS segments
  - `method::String`: standard deviation estimation method

### Returns

  - `_cs_::AbstractArray`: segment averaged and subtracted Cross Spectra
"""
function cs2cs_(
    data::AbstractArray,
    cs::AbstractArray,
    fband::AbstractArray,
    nep::Integer,
    segave::Bool,
    subave::Bool,
    method::String,
)
    if segave
        if method == "bootstrap"
            randboot = rand(1:nep, nep)
            cs_ = dropmean(view(cs, :, randboot, :, :, :), (2, 3))
            av_ = dropmean(view(data, fband, randboot, :, :), (2, 3))
        elseif method == "psi"
            cs_ = dropmean(cs, (2, 3))
            av_ = dropmean(view(data, fband, :, :, :), (2, 3))
        elseif method == "jackknife"
            cs_ = dropmean(cs, 3)
            av_ = dropmean(view(data, fband, :, :, :), 3)
        end
    else
        if method == "bootstrap"
            randboot = rand(1:nep, nep)
            cs_ = dropmean(view(cs, :, randboot, 1, :, :), 2)
            av_ = dropmean(view(data, fband, randboot, 1, :), 2)
        elseif method == "psi"
            cs_ = dropmean(view(cs, :, :, 1, :, :), 2)
            av_ = dropmean(view(data, fband, :, 1, :), 2)
        elseif method == "jackknife"
            cs_ = view(cs, :, :, 1, :, :)
            av_ = view(data, fband, :, 1, :)
        end
    end

    if subave
        if method == "bootstrap" || method == "psi"
            @einsum av_ij[f, i, j] := av_[f, i] * conj(av_[f, j])
        elseif method == "jackknife"
            @einsum av_ij[f, e, i, j] := av_[f, e, i] * conj(av_[f, e, j])
        end
        _cs_ = squeeze(cs_ - av_ij)
    else
        if method == "bootstrap" || method == "psi" || method == "jackknife"
            _cs_ = squeeze(cs_)
        end
    end

    return _cs_
end

"""
    data2psi(data, seglen [, segshift, eplen, freqlist, method,
                             nboot, segave, subave, detrend])

calculates phase slope index (PSI)

### Arguments

  - `data::AbstractArray`: NxM array for N data points in M channels
  - `seglen::Integer`: segment length (determines the frequency resolution)

*optional arguments*

  - `segshift::Integer`: number of bins by which neighboring segments are shifted (default=seglen/2)
  - `eplen::Integer`: length of epochs (if eplen=0, eplen is defaulted to number of samples)
  - `freqlist::AbstractArray`: 2D Array where each column is a frequency band (default is full range)
  - `method::String`: standard deviation estimation method (default is "jackknife")
  - `nboot::Integer`: number of bootstrap resamplings (default is 100)
  - `segave::Bool`: if true, average across CS segments (default is false)
  - `subave::Bool`: if true, subtract average across CS segments (default is false)
  - `detrend::Bool`: if true, performs a 0th-order detrend across raw segments (default is false)
  - `window::Function`: window function with interval length as sole necessary argument (default is Hanning)

### Returns

  - `psi::AbstractArray`: Phase Slope Index with shape (channel, channel, frequency bands)
  - `psi_se::AbstractArray`: PSI estimated standard error with shape (channel, channel, frequency bands)
"""
function data2psi(
    data::AbstractArray,
    seglen::Integer;
    segshift::Integer=0,
    eplen::Integer=0,
    freqlist::AbstractArray=Int64[],
    method::String="jackknife",
    nboot::Integer=100,
    segave::Bool=false,
    subave::Bool=false,
    detrend::Bool=false,
    window::Function=hanning_fun,
)
    (data, nsamples, nchan, eplen, nep, method, subave, segshift, nseg, freqlist, maxfreq, nfbands, n_resample) = data2para(
        data, seglen, segshift, eplen, freqlist, method, subave, nboot
    )

    eposeg = make_eposeg(data, seglen, eplen, nep, nseg, nchan, segshift)

    if detrend
        detrend!(eposeg, 0)
    end

    if isa(window, Function)
        window = window(seglen)
    else
        throw(window * " is not a Function!")
    end
    eposeg .*= window

    eposeg = view(fft(eposeg, 1), 1:maxfreq, :, :, :)

    # preallocation
    psi = Array{Float64}(undef, nchan, nchan, nfbands)
    psi_est = Array{Float64}(undef, nchan, nchan, nfbands, n_resample)
    for (f, fband) in enumerate(eachrow(freqlist'))
        cs_full = data2cs(view(eposeg, fband, :, :, :))

        cs_psi = cs2cs_(eposeg, cs_full, fband, nep, segave, subave, "psi")
        psi[:, :, f] = cs2ps(cs_psi)

        if method == "jackknife"
            cs_jack = cs2cs_(eposeg, cs_full, fband, nep, segave, subave, "jackknife")
            for e in 1:nep
                cs_jack_se = (nep * cs_psi - view(cs_jack, :, e, :, :)) / (nep + 1)
                psi_est[:, :, f, e] = cs2ps(cs_jack_se)
            end

        elseif method == "bootstrap"
            for e in 1:nboot
                cs_boot = cs2cs_(eposeg, cs_full, fband, nep, segave, subave, "bootstrap")
                psi_est[:, :, f, e] = cs2ps(cs_boot)
            end
        end
    end

    psi = squeeze(psi)
    if method == "jackknife"
        psi_se = sqrt(nep) * squeeze(std(psi_est; corrected=true, dims=4))
    elseif method == "bootstrap"
        psi_se = squeeze(std(psi_est; corrected=true, dims=4))
    end

    return psi, psi_se
end

end  # module


text:
export solve!

# Convert from sets used within Convex to MOI sets
function get_MOI_set(cone, length_inds)
    if cone == :SDP
        set = MOI.PositiveSemidefiniteConeTriangle(Int(sqrt(.25 + 2 * length_inds) - .5))
    elseif cone == :Zero
        set = MOI.Zeros(length_inds)
    elseif cone == :Free
        set = MOI.Reals(length_inds)
    elseif cone == :NonNeg
        set = MOI.Nonnegatives(length_inds)
    elseif cone == :NonPos
        set = MOI.Nonpositives(length_inds)
    elseif cone == :SOC
        set = MOI.SecondOrderCone(length_inds)
    elseif cone == :SOCRotated
        set = MOI.RotatedSecondOrderCone(length_inds)
    elseif cone == :ExpPrimal
        set = MOI.ExponentialCone()
    elseif cone == :ExpDual
        set = MOI.DualExponentialCone()
    else
        error("Cone $cone not found somehow; please file an issue.")
    end
    return set
end

# add_terms! turns a matrix or vector `X` into a collection of `VectorAffineTerm`'s, which it adds to a list (the first argument).
# This replicates the behavior of the relation `A[row_indices, col_indices] = X`, where `X` could be a vector or matrix.
# Here, `row_indices` corresponds to the output index of a `VectorAffineTerm`
# Instead of `col_indices`, here we have their wrapped counterparts, `var_indices` (a vector of `MOI.VariableIndex`'s).
function add_terms!(terms::Vector{MOI.VectorAffineTerm{T}}, matrix::SparseMatrixCSC, row_indices, var_indices) where T
    CIs = CartesianIndices((row_indices, 1:length(var_indices)))
    I, J, V = findnz(matrix)
    for  k = eachindex(I, J, V)
        term = MOI.VectorAffineTerm{T}(CIs[I[k],J[k]][1], MOI.ScalarAffineTerm{T}(V[k], var_indices[CIs[I[k],J[k]][2]]))
        push!(terms, term)
    end
end

function add_terms!(terms::Vector{MOI.VectorAffineTerm{T}}, vector::SparseVector, row_indices, var_indices) where T
    CIs = CartesianIndices((row_indices, 1:length(var_indices)))
    I, V = findnz(vector)
    for k = eachindex(I, V)
        term =  MOI.VectorAffineTerm{T}(CIs[I[k]][1], MOI.ScalarAffineTerm{T}(V[k], var_indices[CIs[I[k]][2]]))
        push!(terms, term)
    end
end

function add_terms!(terms::Vector{MOI.VectorAffineTerm{T}}, matrix::Matrix, row_indices, var_indices) where T
    CIs = CartesianIndices((row_indices, 1:length(var_indices)))
    for ci = CartesianIndices(size(matrix))
        term = MOI.VectorAffineTerm{T}(CIs[ci][1], MOI.ScalarAffineTerm{T}(matrix[ci], var_indices[CIs[ci][2]]))
        push!(terms, term)
    end
end

function add_terms!(terms::Vector{MOI.VectorAffineTerm{T}}, vector::Vector, row_indices, var_indices) where T
    CIs = CartesianIndices((row_indices, 1:length(var_indices)))
    for k = 1:length(vector)
        term = MOI.VectorAffineTerm{T}(CIs[k][1], MOI.ScalarAffineTerm{T}(vector[k], var_indices[CIs[k][2]]))
        push!(terms, term)
    end
end

# This function generates a MOI set and VectorAffineFunction corresponding to a constraint that
# Convex represents with a `ConicConstr` object.
# The constraint is of the form Ax + b ∈ set
# where `A` is a matrix, `x` a vector of variables, and `b` a constant vector.
# MOI represents this type of constraint as a VectorAffineFunction, where `A` is represented
# by a collection of `VectorAffineTerm`'s.
function make_MOI_constr(conic_constr::ConicConstr, var_to_indices, id_to_variables, T)
    total_constraint_size = sum(conic_constr.sizes)
    constr_index = 0
    # We create a vector of `VectorAffineTerm`'s to build the constraint.
    terms =  MOI.VectorAffineTerm{T}[]
    b = spzeros(T, total_constraint_size)

    for i = 1:length(conic_constr.objs)
        sz = conic_constr.sizes[i]
        for (id, val) in conic_constr.objs[i]
            if id == objectid(:constant)
                for l in 1:sz
                    b[constr_index + l] = ifelse(val[1][l] == 0, val[2][l],  val[1][l])
                end
            else
                var_indices = var_to_indices[id]
                if id_to_variables[id].sign == ComplexSign()
                    l = length(var_indices) ÷ 2
                    # We create MOI terms and add them to `terms`.
                    # Real part:
                    add_terms!(terms, val[1], constr_index + 1 : constr_index + sz, var_indices[1:l])

                    # Imaginary part:
                    add_terms!(terms, val[2], constr_index + 1 : constr_index + sz, var_indices[l+1 : end])

                else
                    # Only a real part:
                    add_terms!(terms, val[1], constr_index + 1 : constr_index + sz, var_indices)
                end
            end
        end
        constr_index += sz
    end
    set = get_MOI_set(conic_constr.cone, total_constraint_size)
    constr_fn =  MOI.VectorAffineFunction{T}(terms, b)
    return set, constr_fn
end

# This function gets MOI variable indices for each variable used as a key
# inside a `ConicObj` used in the problem.
function get_variable_indices!(model, conic_constraints, id_to_variables)
    var_to_indices = Dict{UInt64,Vector{MOI.VariableIndex}}()
    for conic_constr in conic_constraints
        for i = 1:length(conic_constr.objs)
            for (id, val) in conic_constr.objs[i]
                if !haskey(var_to_indices, id) && id != objectid(:constant)
                    var = id_to_variables[id]
                    if var.sign == ComplexSign()
                        var_size = 2 * length(var)
                    else
                        var_size = length(var)
                    end
                    var_to_indices[id] = MOI.add_variables(model, var_size)
                end
            end
        end
    end
    return var_to_indices
end

# This function traverses the problem tree and loads
# an MOI model with the corresponding problem instance.
# The problem is written in a simple epigraph form.
function load_MOI_model!(model, problem::Problem{T}) where {T}
    if length(problem.objective) != 1
        error("Objective must be a scalar")
    end

    unique_conic_forms = UniqueConicForms()
    objective, objective_var_id = conic_form!(problem, unique_conic_forms)
    conic_constraints = unique_conic_forms.constr_list
    conic_constr_to_constr = unique_conic_forms.conic_constr_to_constr
    id_to_variables = unique_conic_forms.id_to_variables

    # `var_to_indices` maps from variable id to the MOI `VariableIndex`'s corresponding to the variable
    var_to_indices = get_variable_indices!(model, conic_constraints, id_to_variables)

    # the objective: maximize or minimize a scalar variable
    objective_index = var_to_indices[objective_var_id][] # get the `MOI.VariableIndex` corresponding to the objective
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(objective_index))
    MOI.set(model, MOI.ObjectiveSense(), problem.head == :maximize ? MOI.MAX_SENSE : MOI.MIN_SENSE)

    # Constraints: Generate a MOI function and a MOI sets for each `ConicConstr` object in the problem
    MOI_constr_fn = Union{MOI.VectorAffineFunction{T},MOI.SingleVariable}[]
    MOI_sets = Any[]
    for conic_constr in conic_constraints
        set, constr_fn = make_MOI_constr(conic_constr, var_to_indices, id_to_variables, T)
        push!(MOI_sets, set)
        push!(MOI_constr_fn, constr_fn)
    end

    # Add integral and boolean constraints
    for var_id in keys(var_to_indices)
        variable = id_to_variables[var_id]
        if :Int in variable.sets
            var_indices = var_to_indices[var_id]
            for idx = eachindex(var_indices)
                push!(MOI_constr_fn, MOI.SingleVariable(var_indices[idx]))
                push!(MOI_sets, MOI.Integer())
            end
        end
        if :Bin in variable.sets
            var_indices = var_to_indices[var_id]
            for idx in eachindex(var_indices)
                push!(MOI_constr_fn, MOI.SingleVariable(var_indices[idx]))
                push!(MOI_sets, MOI.ZeroOne())
            end
        end
    end

    # Add all the constraints to the model and collect the corresponding MOI indices
    constraint_indices = MOI.add_constraints(model, MOI_constr_fn, MOI_sets)

    return id_to_variables, conic_constr_to_constr, conic_constraints, var_to_indices, constraint_indices
end

function solve!(problem::Problem{T}, optimizer::MOI.ModelLike;
    check_vexity = true,
    verbose = true,
    warmstart = false) where {T}

    if check_vexity
        vex = vexity(problem)
    end

    model = MOIB.full_bridge_optimizer(
        MOIU.CachingOptimizer(
            MOIU.UniversalFallback(MOIU.Model{T}()),
            optimizer
        ),
        T
    )

    id_to_variables, conic_constr_to_constr, conic_constraints, var_to_indices, constraint_indices = load_MOI_model!(model, problem)

    if warmstart
        warmstart_variables!(model, var_to_indices, id_to_variables, T, verbose)
    end

    MOI.optimize!(model)
    problem.model = model

    # populate the status, primal variables, and dual variables (when possible)
    moi_populate_solution!(model, problem, id_to_variables, conic_constr_to_constr, conic_constraints, var_to_indices, constraint_indices)

    if problem.status != MOI.OPTIMAL && verbose
        @warn "Problem status $(problem.status); solution may be inaccurate."
    end
end

function warmstart_variables!(model, var_to_indices, id_to_variables, T, verbose)
    if MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex)
        for (id, var_inds) in pairs(var_to_indices)
            x = id_to_variables[id]
            value = x.value
            value === nothing && continue
            value_vec = packvec(value, sign(x) == ComplexSign())
            value_vec = convert(Vector{T}, value_vec)
            MOI.set(model, MOI.VariablePrimalStart(), var_inds, value_vec)
        end
    elseif verbose
        @warn "Skipping variable warmstart; the solver does not support it."
    end
end

function moi_populate_solution!(model::MOI.ModelLike, problem, id_to_variables, conic_constr_to_constr, conic_constraints, var_to_indices, constraint_indices)
    status = MOI.get(model, MOI.TerminationStatus())
    problem.status = status

    ## This can throw if the problem was not solved correctly
    # objective = MOI.get(model, MOI.ObjectiveValue())
    # problem.optval = objective
    ## Instead, we use a `getproperty` override to access the objective value only when the user asks.
    ## That way `moi_populate_solution!` will complete successfully.

    dual_status = MOI.get(model, MOI.DualStatus())
    primal_status = MOI.get(model, MOI.PrimalStatus())

    if primal_status != MOI.NO_SOLUTION
        for (id, var_indices) in var_to_indices
            var = id_to_variables[id]
            vectorized_value =  MOI.get(model, MOI.VariablePrimal(), var_indices)
            var.value = unpackvec(vectorized_value, size(var), var.sign == ComplexSign())
        end
    else
        for (id, var_indices) in var_to_indices
            var = id_to_variables[id]
            var.value = nothing
        end
    end

    if dual_status != MOI.NO_SOLUTION
        for (idx, conic_constr) in enumerate(conic_constraints)
            haskey(conic_constr_to_constr, conic_constr) || continue
            constr = conic_constr_to_constr[conic_constr]
            MOI_constr_indices = constraint_indices[idx]
            dual_value_vectorized = MOI.get(model, MOI.ConstraintDual(), MOI_constr_indices)
            iscomplex = sign(constr.lhs) == ComplexSign() || sign(constr.rhs) == ComplexSign()
            constr.dual = unpackvec(dual_value_vectorized, constr.size, iscomplex)
        end
    else
        for (idx, conic_constr) in enumerate(conic_constraints)
            haskey(conic_constr_to_constr, conic_constr) || continue
            constr = conic_constr_to_constr[conic_constr]
            constr.dual = nothing
        end
    end

end

# This is type unstable!
function unpackvec(v::AbstractVector, size::Tuple{Int,Int}, iscomplex::Bool)
    if iscomplex && length(v) == 2
        return v[1] + im * v[2]
    elseif iscomplex
        l = length(v) ÷ 2
        # use views?
        return reshape(v[1:l], size) + im * reshape(v[l+1 : end], size)
    elseif length(v) == 1
        return v[]
    else
        return reshape(v, size)
    end
end

# In `packvec`, we use `real` even in the `iscomplex == false` branch
# for type stability. Note here `iscomplex` refers to the `Sign` of the variable
# associated to the values here, not whether or not the variable actually holds
# a complex number at this point. Hence, one could have `iscomplex==true` and
# `isreal(value)==true` or even `value isa Real` could hold.
function packvec(value::Number, iscomplex::Bool)
    iscomplex ? [real(value), imag(value)] : [real(value)]
end

function packvec(value::AbstractArray, iscomplex::Bool)
    value = reshape(value, length(value))
    if iscomplex
        return [real(value); imag(value)]
    else
        return real(value)
    end
end


text:
@enum BFInstruction inc dec left right accept output debug

struct Loop
  body::Vector{Union{Loop,BFInstruction}}
end

function parsechar(s, i)
  ch = s[i]
  ch == '+' ? (i+1, inc) :
  ch == '-' ? (i+1, dec) :
  ch == '>' ? (i+1, right) :
  ch == '<' ? (i+1, left) :
  ch == '.' ? (i+1, output) :
  ch == ',' ? (i+1, accept) :
  ch == '#' ? (i+1, debug) :
  ch == '[' ? parseloop(s, i) :
  (i+1, nothing)
end

function parseloop(s, i)
  body = Union{Loop,BFInstruction}[]
  i += 1
  while true
    s[i] == ']' && break
    i, op = parsechar(s, i)
    op === nothing || push!(body, op)
  end
  return i+1, Loop(body)
end

function bfparse(s, i = 1)
  is = Union{Loop,BFInstruction}[]
  while i <= length(s)
    i, op = parsechar(s, i)
    op === nothing || push!(is, op)
  end
  return is
end

function debugtape(tape, ptr)
  for i = 1:length(tape)
    print(tape[i])
    i == ptr && print(*)
    print(" ")
  end
  println()
end

function interploop(tape, loop, ptr)
  while tape[ptr] != 0
    ptr = interp(tape, loop.body, ptr)
  end
  return ptr
end

function interp(tape, ops, ptr)
  for op in ops
    if op === inc
      tape[ptr] += 1
    elseif op === dec
      tape[ptr] -= 1
    elseif op === left
      ptr -= 1
    elseif op === right
      ptr += 1
    elseif op === debug
      debugtape(tape, ptr)
    elseif op isa Loop
      ptr = interploop(tape, op, ptr)
    end
  end
  return ptr
end

function interpret(ops, tape = zeros(UInt32, 2^16))
  interp(tape, ops, 1)
  return tape
end

interpret(ops::String, args...) = interpret(bfparse(ops, 1), args...)

mutable struct NVector{T,S} <: AbstractVector{T}
  data::NTuple{S,T}
end

Base.size(xs::NVector) = (length(xs.data),)
Base.getindex(xs::NVector, i::Integer) = xs.data[i]
Base.setindex!(xs::NVector, v, i::Integer)  =
  xs.data = ntuple(j -> j == i ? v : xs.data[j], length(xs.data))


text:
using DiffEqBase: AbstractODEProblem, remake
using StaticArrays: SVector
using Setfield: get, set

using ..Continuations: init, solve!, residual_jacobian, ContinuationSolver,
    ZeroNotFoundError, FindZeroInputError
using ..BifurcationsBase: AbstractSpecialPoint, SpecialPoint, special_intervals
using ..Codim1
using ..Codim1: resolving_points, resolve_point
using ..Codim2
using ..Codim2: Codim2Solver, DiffEqCodim2Problem, first_lyapunov_coefficient
using ..Codim1LimitCycle: Codim1LCSolver, diameter

# See also:
# * [[../codim2/problem.jl::function BifurcationProblem]]
# * [[../codim2/switch.jl::function BifurcationProblem]]
# * [[../codim1lc/factories.jl::function LimitCycleProblem]]
function FoldLimitCycleProblem(
        point::AbstractSpecialPoint,
        solver::Codim2Solver;
        via = :linear,
        lc_solver_opts = [],
        kwargs...)
    @assert point.point_type == Codim2.PointTypes.bautin
    @assert timekind(point) isa Continuous

    solver.prob :: DiffEqCodim2Problem
    de_prob = solver.prob.de_prob :: AbstractODEProblem

    if via == :fold_lc
        # TODO: Stop solving sub-problems in the constructor:
        lc_prob = bautin_to_subcritical_lc(point, solver; kwargs...)
        u_init = nearest_fold_lc(lc_prob; lc_solver_opts...)
    elseif via == :linear
        lc_prob = lc_at_bautin(point, solver; kwargs...)
        u_init = get_u0(lc_prob)
    else
        error("Unsupported option: via = $via (∉ {:fold_lc, :linear})")
    end

    prob_cache = get_prob_cache(lc_prob)
    _H, J_init = residual_jacobian(u_init, prob_cache)

    n = length(lc_prob.xs0) + 1  # +1 to include period
    vals, vecs = eigen(@view J_init[1:n, 1:n])
    _, idx0 = findmin(abs.(vals))  # closest to zero
    opts = as(solver, ContinuationSolver).opts
    if any(x -> abs(imag(x)) > opts.atol, @view vecs[:, idx0])
        @warn "Nearest-to-null vector is complex (eigval: $(vals[idx0]))"
    end
    v0 = normalize(real(@view vecs[:, idx0]))

    xs0 = reshape(u_init[1:length(lc_prob.xs0)], size(lc_prob.xs0))
    l0 = u_init[length(lc_prob.xs0) + 1]
    vs0 = reshape(v0[1:end-1], size(lc_prob.xs0))
    dl0 = v0[end]
    t0 = SVector(u_init[end],
                 get(lc_prob.de_prob.p, solver.prob.param_axis2))

    return FoldLimitCycleProblem(
        lc_prob;
        xs0 = xs0,
        l0 = l0,
        vs0 = vs0,
        dl0 = dl0,
        t0 = t0,
        param_axis1 = solver.prob.param_axis1,
        param_axis2 = solver.prob.param_axis2,
        t_domain = solver.prob.t_domain,
    )
end

function BifurcationProblem(point::AbstractSpecialPoint,
                            solver::Codim1LCSolver,
                            param_axis2::Lens,
                            t2_domain::Tuple;
                            kwargs...)

    cd1_prob = solver.prob :: LimitCycleProblem
    de_prob = cd1_prob.de_prob :: AbstractODEProblem

    xs0 = reshape(point.u[1:end-2], size(cd1_prob.xs0))

    # Initialize perturbation direction
    vs0 = xs0 .- mean(xs0, dims=2)
    normalize!(view(vs0, :))

    l0 = point.u[end - 1]
    t0 = SVector(point.u[end], get(de_prob.p, param_axis2))
    t_domain = (
        SVector(cd1_prob.t_domain[1], t2_domain[1]),
        SVector(cd1_prob.t_domain[2], t2_domain[2]),
    )
    param_axis1 = cd1_prob.param_axis

    return FoldLimitCycleProblem(
        cd1_prob;
        xs0 = xs0,
        l0 = l0,
        vs0 = vs0,
        dl0 = 0.0,  # ???
        t0 = t0,
        param_axis1 = param_axis1,
        param_axis2 = param_axis2,
        t_domain = t_domain,
        kwargs...)
end

function nearest_fold_lc(lc_prob; lc_solver_opts...)
    lc_solver = init(
        lc_prob;
        start_from_nearest_root = true,
        max_branches = 0,
        nominal_angle_rad = 2π * (1 / 360),  # TODO: should it be the default?
        lc_solver_opts...
    )
    solve!(lc_solver)

    PT = Codim1LimitCycle.PointTypes.saddle_node
    point_candidates = special_intervals(lc_solver, PT)
    if isempty(point_candidates)
        error("Cannot find fold bifurcation of the limit cycle.")
    end
    sw1_points = special_intervals(lc_solver.sol.sweeps[1], PT)
    sw2_points = special_intervals(lc_solver.sol.sweeps[2], PT)
    sweep_index = -1
    if length(sw2_points) == 0
        @assert length(sw1_points) > 0
        flc_point, = sw1_points
        sweep_index = 1
    elseif length(sw1_points) == 0
        @assert length(sw2_points) > 0
        flc_point, = sw2_points
        sweep_index = 2
    else
        # Choose increasing sweep
        cont_sweeps = as(lc_solver, ContinuationSolver).sol.sweeps
        at_most = (n, u) -> u[1:max(length(u), n)]
        ds1 = [diameter(u, lc_prob) for u in at_most(10, cont_sweeps[1].u)]
        ds2 = [diameter(u, lc_prob) for u in at_most(10, cont_sweeps[2].u)]
        mean_dd1 = mean(diff(ds1))
        mean_dd2 = mean(diff(ds2))
        if mean_dd1 < 0 && mean_dd2 < 0
            @warn "Both sweeps have decreasing diameters."
        end
        if mean_dd1 > mean_dd2
            flc_point, = sw1_points
            sweep_index = 1
        else
            flc_point, = sw2_points
            sweep_index = 2
        end
    end

    # TODO?: resolve flc_point
    u_init = (flc_point.u0 .+ flc_point.u1) ./ 2

    if length(point_candidates) > 1
        @warn """
        $(length(point_candidates)) fold bifurcations are found.
        Using the first point found in sweep $sweep_index:
            $(lc_prob.param_axis) : $(u_init[end])
            diameter : $(diameter(u_init, lc_prob))
        """
    end

    return u_init
end

function bautin_to_subcritical_lc(
        point::AbstractSpecialPoint,
        solver::Codim2Solver;
        kwargs...)

    @assert point.point_type == Codim2.PointTypes.bautin
    @assert timekind(point) isa Continuous

    solver.prob :: DiffEqCodim2Problem
    de_prob = solver.prob.de_prob :: AbstractODEProblem
    prob_cache = as(solver.cache, ContinuationCache).prob_cache

    if point isa SpecialPoint
        point = let
            if point.sweep.value == nothing
                error("""
                `point` is not associated with a sweep.
                point = $(point)
                 """)
            end
            # See: [[../base/solver.jl::SpecialPointInterval]]
            sweep = as(point.sweep.value, ContinuationSweep)
            u0 = sweep.u[point.point_index - 1]
            u1 = sweep.u[point.point_index]
            _H, J1 = residual_jacobian(u1, prob_cache)
            SpecialPointInterval(
                timekind(point),
                point.point_type,
                point.point_index,
                u0,
                u1,
                J1,
                WeakRef(sweep),
            )
        end
    else
        point :: SpecialPointInterval
    end

    # TODO: Specify target l1 (> 0) value and step into it using zero
    # finder (so that hopefully it's not too close to/far away from
    # the Bautin point).
    coeff_1 = first_lyapunov_coefficient(prob_cache, point.u1, point.J1)
    if coeff_1 > 0  # u1 is subcritical
        u_sc = point.u1
        J_sc = point.J1
    else
        u_sc = point.u0
        _H, J_sc = residual_jacobian(point.u0, prob_cache)
    end
    coeff_0 = first_lyapunov_coefficient(
        prob_cache, point.u0, residual_jacobian(point.u0, prob_cache)[2])
    @assert coeff_0 * coeff_1 < 0

    # [[../codim2/switch.jl::Manually cast]]
    xtype = typeof(de_prob.u0)

    N = length(de_prob.u0)
    x0 = xtype(u_sc[1:N])
    t0 = u_sc[end-1]

    hopf_point = SpecialPoint(
        timekind(point),
        Codim1.PointTypes.hopf,
        -1,  # dummy index
        vcat(x0, t0[1]),
        J_sc[1:N, 1:N+1],
        WeakRef(),  # dummy
    )
    param = set(de_prob.p, solver.prob.param_axis2, u_sc[end])

    return LimitCycleProblem(
        hopf_point,
        solver.opts;
        de_prob = remake(de_prob; p = param),
        param_axis = solver.prob.param_axis1,
        t_domain = (solver.prob.t_domain[1][1],
                    solver.prob.t_domain[2][1]),
        # phase_space = solver.prob.phase_space,
        kwargs...)
end

function lc_at_bautin(
        point::AbstractSpecialPoint,
        solver::Codim2Solver;
        kwargs...)

    @assert point.point_type == Codim2.PointTypes.bautin
    @assert timekind(point) isa Continuous

    solver.prob :: DiffEqCodim2Problem
    de_prob = solver.prob.de_prob :: AbstractODEProblem

    resolved = resolve_point(point, solver)

    # [[../codim2/switch.jl::Manually cast]]
    xtype = typeof(de_prob.u0)

    N = length(de_prob.u0)
    x0 = xtype(resolved.u[1:N])
    t0 = resolved.u[end-1]

    hopf_point = SpecialPoint(
        timekind(point),
        Codim1.PointTypes.hopf,
        -1,  # dummy index
        vcat(x0, t0[1]),
        resolved.J[1:N, 1:N+1],
        WeakRef(),  # dummy
    )
    param = set(de_prob.p, solver.prob.param_axis2, resolved.u[end])

    return LimitCycleProblem(
        hopf_point,
        solver.opts;
        de_prob = remake(de_prob; p = param),
        param_axis = solver.prob.param_axis1,
        t_domain = (solver.prob.t_domain[1][1],
                    solver.prob.t_domain[2][1]),
        # phase_space = solver.prob.phase_space,
        kwargs...)
end


text:
using Statistics

function fuelcost(steps)
    cost = 0
    stepcost = 1
    for i in 1:steps
        cost += stepcost
        stepcost += 1
    end
    return cost
end

function solve(filename, timesteps = 80)
    #load txt file
    values = readlines(filename*".txt")
    values = split(values[1], ',')
    positions = parse.(Int, values)
    goalpos = median(positions)
    avgpos = round(mean(positions))
    cost = 0
    costavg = 0
    costavgsmall = 0
    costavgbig = 0
    for p in positions
        cost += abs(p - goalpos)
        costavg += fuelcost(abs(p - avgpos))
        costavgsmall += fuelcost(abs(p-(avgpos-1)))
        costavgbig += fuelcost(abs(p-(avgpos+1)))
    end
    costavg = min(costavg, costavgsmall, costavgbig)
    println("Median: "*string(goalpos))
    println("Accumulated fuel cost: "*string(cost))
    println("Average: "*string(avgpos))
    println("Avg fuel cost accumulated: "*string(costavg))
end

solve("smalltest")
solve("input")


text:
include("perm.jl")

using DrakeVisualizer

DrakeVisualizer.new_window()

model = Visualizer(poly2)


text:
using Manopt, Manifolds, ManifoldsBase, Test

@testset "Test primal dual plan" begin
    #
    # Perform an really easy test, just compute a mid point
    #
    pixelM = Sphere(2)
    M = PowerManifold(pixelM, NestedPowerRepresentation(), 2)
    data = [[1.0, 0.0, 0.0], 1 / sqrt(2) .* [1.0, 1.0, 0.0]]
    α = 1
    # known minimizer
    δ = min(α / distance(pixelM, data[1], data[2]), 0.5)
    x_hat = shortest_geodesic(M, data, reverse(data), δ)
    N = TangentBundle(M)
    fidelity(M, x) = 1 / 2 * distance(M, x, f)^2
    Λ(M, x) = ProductRepr(x, forward_logs(M, x))
    function Λ!(M, Y, x)
        N = TangentBundle(M)
        copyto!(M, Y[N, :point], x)
        zero_tangent_vector!(N.manifold, Y[N, :vector], Y[N, :point])
        forward_logs!(M, Y[N, :vector], x)
        return Y
    end
    prior(M, x) = norm(norm.(Ref(pixelM), x, (Λ(M, x))[N, :vector]), 1)
    cost(M, x) = (1 / α) * fidelity(M, x) + prior(M, x)
    prox_F(M, λ, x) = prox_distance(M, λ / α, data, x, 2)
    prox_F!(M, y, λ, x) = prox_distance!(M, y, λ / α, data, x, 2)
    function prox_G_dual(N, n, λ, ξ)
        return ProductRepr(
            ξ[N, :point],
            project_collaborative_TV(
                base_manifold(N), λ, n[N, :point], ξ[N, :vector], Inf, Inf, 1.0
            ),
        )
    end
    function prox_G_dual!(N, η, n, λ, ξ)
        η[N, :point] .= ξ[N, :point]
        project_collaborative_TV!(
            base_manifold(N), η[N, :vector], λ, n[N, :point], ξ[N, :vector], Inf, Inf, 1.0
        )
        return η
    end
    DΛ(M, m, X) = ProductRepr(zero_tangent_vector(M, m), differential_forward_logs(M, m, X))
    function DΛ!(M, Y, m, X)
        N = TangentBundle(M)
        zero_tangent_vector!(M, Y[N, :point], m)
        differential_forward_logs!(M, Y[N, :vector], m, X)
        return Y
    end
    function adjoint_DΛ(N, m, n, Y)
        return adjoint_differential_forward_logs(N.manifold, m, Y[N, :vector])
    end
    function adjoint_DΛ!(N, X, m, n, Y)
        return adjoint_differential_forward_logs!(N.manifold, X, m, Y[N, :vector])
    end

    m = fill(mid_point(pixelM, data[1], data[2]), 2)
    n = Λ(M, m)
    x0 = deepcopy(data)
    ξ0 = ProductRepr(zero_tangent_vector(M, m), zero_tangent_vector(M, m))

    p_exact = PrimalDualProblem(M, N, cost, prox_F, prox_G_dual, adjoint_DΛ; Λ=Λ)
    p_linearized = PrimalDualProblem(
        M, N, cost, prox_F, prox_G_dual, adjoint_DΛ; linearized_forward_operator=DΛ
    )
    o_exact = ChambollePockOptions(m, n, x0, ξ0; variant=:exact)
    o_linearized = ChambollePockOptions(m, n, x0, ξ0; variant=:linearized)
    n_old = ProductRepr(n[N, :point], n[N, :vector])
    x_old = copy(x0)
    ξ_old = ProductRepr(ξ0[N, :point], ξ0[N, :vector])

    @testset "test Mutating/Allocation Problem Variants" begin
        p1 = PrimalDualProblem(
            M, N, cost, prox_F, prox_G_dual, adjoint_DΛ; linearized_forward_operator=DΛ, Λ=Λ
        )
        p2 = PrimalDualProblem(
            M,
            N,
            cost,
            prox_F!,
            prox_G_dual!,
            adjoint_DΛ!;
            linearized_forward_operator=DΛ!,
            Λ=Λ!,
            evaluation=MutatingEvaluation(),
        )
        x1 = get_primal_prox(p1, 1.0, x0)
        x2 = get_primal_prox(p2, 1.0, x0)
        @test x1 == x2
        get_primal_prox!(p1, x1, 0.8, x0)
        get_primal_prox!(p2, x2, 0.8, x0)
        @test x1 == x2

        ξ1 = get_dual_prox(p1, n, 1.0, ξ0)
        ξ2 = get_dual_prox(p2, n, 1.0, ξ0)
        @test ξ1[N, :point] == ξ2[N, :point]
        @test ξ1[N, :vector] == ξ2[N, :vector]
        get_dual_prox!(p1, ξ1, n, 1.0, ξ0)
        get_dual_prox!(p2, ξ2, n, 1.0, ξ0)
        @test ξ1[N, :point] == ξ2[N, :point]
        @test ξ1[N, :vector] == ξ2[N, :vector]

        y1 = forward_operator(p1, x0)
        y2 = forward_operator(p2, x0)
        @test y1[N, :point][1] == y2[N, :point][1]
        @test y1[N, :point][2] == y2[N, :point][2]
        @test y1[N, :vector][1] == y2[N, :vector][1]
        @test y1[N, :vector][2] == y2[N, :vector][2]
        forward_operator!(p1, y1, x0)
        forward_operator!(p2, y2, x0)
        @test y1[N, :point][1] == y2[N, :point][1]
        @test y1[N, :point][2] == y2[N, :point][2]
        @test y1[N, :vector][1] == y2[N, :vector][1]
        @test y1[N, :vector][2] == y2[N, :vector][2]

        X = log(M, m, x0)
        Y1 = linearized_forward_operator(p1, m, X, n)
        Y2 = linearized_forward_operator(p2, m, X, n)
        @test Y1[N, :point] == Y2[N, :point]
        @test Y1[N, :vector] == Y2[N, :vector]
        linearized_forward_operator!(p1, Y1, m, X, n)
        linearized_forward_operator!(p2, Y2, m, X, n)
        @test Y1[N, :point] == Y2[N, :point]
        @test Y1[N, :vector] == Y2[N, :vector]

        Z1 = adjoint_linearized_operator(p1, m, n, ξ0)
        Z2 = adjoint_linearized_operator(p2, m, n, ξ0)
        @test Z1 == Z2
        adjoint_linearized_operator!(p1, Z1, m, n, ξ0)
        adjoint_linearized_operator!(p2, Z2, m, n, ξ0)
        @test Z1 == Z2
    end
    @testset "Primal/Dual residual" begin
        p_exact = PrimalDualProblem(M, N, cost, prox_F, prox_G_dual, adjoint_DΛ; Λ=Λ)
        p_linearized = PrimalDualProblem(
            M, N, cost, prox_F, prox_G_dual, adjoint_DΛ; linearized_forward_operator=DΛ
        )
        o_exact = ChambollePockOptions(m, n, x0, ξ0; variant=:exact)
        o_linearized = ChambollePockOptions(m, n, x0, ξ0; variant=:linearized)
        @test primal_residual(p_exact, o_exact, x_old, ξ_old, n_old) ≈ 0 atol = 1e-16
        @test primal_residual(p_linearized, o_linearized, x_old, ξ_old, n_old) ≈ 0 atol =
            1e-16
        @test dual_residual(p_exact, o_exact, x_old, ξ_old, n_old) ≈ 4.0 atol = 1e-16
        @test dual_residual(p_linearized, o_linearized, x_old, ξ_old, n_old) ≈ 0.0 atol =
            1e-16

        step_solver!(p_exact, o_exact, 1)
        step_solver!(p_linearized, o_linearized, 1)
        @test primal_residual(p_exact, o_exact, x_old, ξ_old, n_old) > 0
        @test primal_residual(p_linearized, o_linearized, x_old, ξ_old, n_old) > 0
        @test dual_residual(p_exact, o_exact, x_old, ξ_old, n_old) > 4.0
        @test dual_residual(p_linearized, o_linearized, x_old, ξ_old, n_old) > 0

        o_err = ChambollePockOptions(m, n, x0, ξ0; variant=:err)
        @test_throws DomainError dual_residual(p_exact, o_err, x_old, ξ_old, n_old)
    end
    @testset "Debug prints" begin
        a = StoreOptionsAction((:x, :ξ, :n, :m))
        update_storage!(a, Dict(:x => x_old, :ξ => ξ_old, :n => n_old, :m => copy(m)))
        io = IOBuffer()

        d1 = DebugDualResidual(a, io)
        d1(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "Dual Residual:")

        d2 = DebugPrimalResidual(a, io)
        d2(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "Primal Residual: ")

        d3 = DebugPrimalDualResidual(a, io)
        d3(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "PD Residual: ")

        d4 = DebugPrimalChange(a, "Primal Change: ", io)
        d4(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "Primal Change: ")

        d5 = DebugPrimalIterate(io)
        d5(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "x:")

        d6 = DebugDualIterate(io)
        d6(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "ξ:")

        d7 = DebugDualChange(a, io)
        d7(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "Dual Change:")

        d7a = DebugDualChange((ξ0, n), a, io)
        d7a(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "Dual Change:")

        d8 = DebugDualBaseIterate(io)
        d8(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "n:")

        d9 = DebugDualBaseChange(a, io)
        d9(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "Dual Base Change:")

        d10 = DebugPrimalBaseIterate(io)
        d10(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "m:")

        d11 = DebugPrimalBaseChange(a, io)
        d11(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "Primal Base Change:")

        d12 = DebugDualResidual((x0, ξ0, n), a, io)
        d12(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "Dual Residual:")

        d13 = DebugPrimalDualResidual((x0, ξ0, n), a, io)
        d13(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "PD Residual:")

        d14 = DebugPrimalResidual((x0, ξ0, n), a, io)
        d14(p_exact, o_exact, 1)
        s = String(take!(io))
        @test startswith(s, "Primal Residual:")
    end
    @testset "Records" begin
        a = StoreOptionsAction((:x, :ξ, :n, :m))
        update_storage!(a, Dict(:x => x_old, :ξ => ξ_old, :n => n_old, :m => copy(m)))
        io = IOBuffer()

        for r in [
            RecordPrimalChange(),
            RecordPrimalIterate(x0),
            RecordDualIterate(ξ0),
            RecordDualChange(),
            RecordDualBaseIterate(n),
            RecordDualBaseChange(),
            RecordPrimalBaseIterate(x0),
            RecordPrimalBaseChange(),
        ]
            r(p_exact, o_exact, 1)
            @test length(get_record(r)) == 1
        end
    end
end


text:
module ExteriorCalculus
export Ob, Hom, dom, codom, compose, ⋅, id,
  otimes, ⊗, munit, braid, oplus, ⊕, mzero, swap, coproduct, ⊔,
  mcopy, Δ, delete, ◊, plus, +, zero, antipode,
  MetricFreeExtCalc1D, MetricFreeExtCalc2D, ExtCalc1D, ExtCalc2D, FreeExtCalc2D,
  Space, Chain0, Chain1, Chain2, Form0, Form1, Form2,
  ∂₁, ∂₂, d₀, d₁, ∧₀₀, ∧₁₀, ∧₀₁, ∧₁₁, ∧₂₀, ∧₀₂, ι₁, ι₂, ℒ₀, ℒ₁, ℒ₂,
  DualForm0, DualForm1, DualForm2, ⋆₀, ⋆₁, ⋆₂, ⋆₀⁻¹, ⋆₁⁻¹, ⋆₂⁻¹,
  dual_d₀, dual_d₁, δ₁, δ₂, ∇²₀, ∇²₁, Δ₀, Δ₁, Δ₂

using Catlab, Catlab.Theories
import Catlab.Theories: Ob, Hom, dom, codom, compose, ⋅, id,
  otimes, ⊗, munit, braid, oplus, ⊕, mzero, swap, coproduct,
  mcopy, Δ, delete, ◊, plus, +, zero, antipode

""" Theory of additive (symmetric) monoidal categories.

Additive monoidal categories are additive categories with a symmetric monoidal
product that is bi-additive, i.e., additive in each variable separately.
Sometimes these are called "tensor categories" but that term is ambiguous, as
noted at the nLab. Not to be confused with an SMC that happens to be written
additively (`SymmetricMonoidalCategoryAdditive` in Catlab).

TODO: Migrate to `Catlab.Theories.MonoidalMultiple`.
"""
@theory AdditiveMonoidalCategory{Ob,Hom} <: AdditiveCategory{Ob,Hom} begin
  # Tensor product.
  otimes(A::Ob, B::Ob)::Ob
  otimes(f::(A → B), g::(C → D))::((A ⊗ C) → (B ⊗ D)) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob)
  @op (⊗) := otimes
  munit()::Ob
  braid(A::Ob, B::Ob)::((A ⊗ B) → (B ⊗ A))
  @op (σ) := braid
  # TODO: Tensor product axioms should be inherited along with SMC operators.

  # Distributivity of tensor products over direct sums.
  #
  # In the standard categories like R-Mod, these would be distributors
  # (distributivity natural isomorphisms) but we replace them with equalities
  # following our usual approach to coherence maps like associators and unitors.
  (A ⊕ B) ⊗ C == (A ⊗ C) ⊕ (B ⊗ C) ⊣ (A::Ob, B::Ob, C::Ob)
  mzero() ⊗ A == mzero() ⊣ (A::Ob)
  A ⊗ (B ⊕ C) == (A ⊗ B) ⊕ (A ⊗ C) ⊣ (A::Ob, B::Ob, C::Ob)
  A ⊗ mzero() == mzero() ⊣ (A::Ob)

  # Bi-additivity.
  (f + g) ⊗ h == (f ⊗ h) + (g ⊗ h) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob, f::(A → B), g::(A → B), h::(C → D))
  (delete(A) ⋅ zero(B)) ⊗ f == delete(A⊗C) ⋅ zero(B⊗D) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob, f::(C → D))
  f ⊗ (g + h) == (f ⊗ g) + (f ⊗ h) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob, f::(A → B), g::(C → D), h::(C → D))
  f ⊗ (delete(C) ⋅ zero(D)) == delete(A⊗C) ⋅ zero(B⊗D) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob, f::(A → B))
end

""" Base theory for calculus on manifold-like spaces.

This non-standard theory is for manifold-like spaces (`Space`) equipped with
objects (`Ob`) and morphisms (`Hom`) belonging to an additive category, say of
real vector spaces. The objects are spaces of things like vector fields, chains,
forms, and twisted forms. Elements of these spaces, e.g. particular forms, are
all assumed to be smoothly time-dependent and so are equipped with time
derivative operators. Unlike in the space-time exterior calculus, we reserve the
exterior calculus for the spatial dimensions and handle time separately.
"""
@theory ManifoldCalculus{Ob,Hom,Space} <: AdditiveMonoidalCategory{Ob,Hom} begin
  Space::TYPE

  # Coproduct of spaces. TODO: Is it worth fully axiomatizing?
  coproduct(X::Space, Y::Space)::Space
  @op (⊔) := coproduct

  """ Partial derivative with respect to time, a linear operator.
  """
  ∂ₜ(A::Ob)::Hom(A,A)
end

# Metric-free exterior calculus
###############################

""" Theory of exterior calculus on 1-or-higher-dimensional manifold-like spaces.
"""
@theory MetricFreeExtCalc1D₊{Ob,Hom,Space} <: ManifoldCalculus{Ob,Hom,Space} begin
  Chain0(X::Space)::Ob
  Chain1(X::Space)::Ob
  ∂₁(X::Space)::Hom(Chain1(X), Chain0(X))

  Form0(X::Space)::Ob
  Form1(X::Space)::Ob
  Form0(X⊔Y) == Form0(X)⊕Form0(Y) ⊣ (X::Space, Y::Space)
  Form1(X⊔Y) == Form1(X)⊕Form1(Y) ⊣ (X::Space, Y::Space)
  d₀(X::Space)::Hom(Form0(X), Form1(X))

  ∧₀₀(X::Space)::Hom(Form0(X)⊗Form0(X), Form0(X))
  ∧₁₀(X::Space)::Hom(Form1(X)⊗Form0(X), Form1(X))
  ∧₀₁(X::Space)::Hom(Form0(X)⊗Form1(X), Form1(X))
  σ(Form0(X),Form0(X)) ⋅ ∧₀₀(X) == ∧₀₀(X) ⊣ (X::Space)
  σ(Form0(X),Form1(X)) ⋅ ∧₁₀(X) == ∧₀₁(X) ⊣ (X::Space)
  σ(Form1(X),Form0(X)) ⋅ ∧₀₁(X) == ∧₁₀(X) ⊣ (X::Space)

  ι₁(X::Space)::Hom(Form1(X)⊗Form1(X), Form0(X))
  ℒ₀(X::Space)::Hom(Form1(X)⊗Form0(X), Form0(X))
  ℒ₀(X) == (id(Form1(X))⊗d₀(X)) ⋅ ι₁(X) ⊣ (X::Space)
  ℒ₁(X::Space)::Hom(Form1(X)⊗Form1(X), Form1(X))
end

""" Theory of exterior caclulus on 1D manifold-like spaces.
"""
@theory MetricFreeExtCalc1D{Ob,Hom,Space} <: MetricFreeExtCalc1D₊{Ob,Hom,Space} begin
  ℒ₁(X) == ι₁(X) ⋅ d₀(X) ⊣ (X::Space)
end

""" Theory of exterior calculus on 2-or-higher-dimensional manifold-like spaces.
"""
@theory MetricFreeExtCalc2D₊{Ob,Hom,Space} <: MetricFreeExtCalc1D₊{Ob,Hom,Space} begin
  Chain2(X::Space)::Ob
  ∂₂(X::Space)::Hom(Chain2(X), Chain1(X))
  ∂₂(X) ⋅ ∂₁(X) == zero(Chain2(X), Chain0(X)) ⊣ (X::Space)

  Form2(X::Space)::Ob
  Form2(X⊔Y) == Form2(X)⊕Form2(Y) ⊣ (X::Space, Y::Space)
  d₁(X::Space)::Hom(Form1(X), Form2(X))
  d₀(X) ⋅ d₁(X) == zero(Form0(X), Form2(X)) ⊣ (X::Space)

  ∧₁₁(X::Space)::Hom(Form1(X)⊗Form1(X), Form2(X))
  ∧₂₀(X::Space)::Hom(Form2(X)⊗Form0(X), Form2(X))
  ∧₀₂(X::Space)::Hom(Form0(X)⊗Form2(X), Form2(X))
  σ(Form1(X),Form1(X)) ⋅ ∧₁₁(X) == ∧₁₁(X) ⋅ antipode(Form2(X)) ⊣ (X::Space)
  σ(Form0(X),Form2(X)) ⋅ ∧₂₀(X) == ∧₀₂(X) ⊣ (X::Space)
  σ(Form2(X),Form0(X)) ⋅ ∧₀₂(X) == ∧₂₀(X) ⊣ (X::Space)

  ι₂(X::Space)::Hom(Form1(X)⊗Form2(X), Form1(X))
  ℒ₁(X) == (ι₁(X) ⋅ d₀(X)) + ((id(Form1(X))⊗d₁(X)) ⋅ ι₂(X)) ⊣ (X::Space)
  ℒ₂(X::Space)::Hom(Form1(X)⊗Form2(X), Form2(X))
end

""" Theory of exterior calculus on 2D manifold-like spaces.
"""
@theory MetricFreeExtCalc2D{Ob,Hom,Space} <: MetricFreeExtCalc2D₊{Ob,Hom,Space} begin
  ℒ₂(X) == ι₂(X) ⋅ d₁(X) ⊣ (X::Space)
end

# Exterior calculus (with metric)
#################################

""" Theory of exterior calculus on 1D Riemannian manifold-like space.
"""
@theory ExtCalc1D{Ob,Hom,Space} <: MetricFreeExtCalc1D{Ob,Hom,Space} begin
  DualForm0(X::Space)::Ob
  DualForm1(X::Space)::Ob
  DualForm0(X⊔Y) == DualForm0(X)⊕DualForm0(Y) ⊣ (X::Space, Y::Space)
  DualForm1(X⊔Y) == DualForm1(X)⊕DualForm1(Y) ⊣ (X::Space, Y::Space)
  dual_d₀(X::Space)::Hom(DualForm0(X), DualForm1(X))

  ⋆₀(X::Space)::Hom(Form0(X), DualForm1(X))
  ⋆₁(X::Space)::Hom(Form1(X), DualForm0(X))
  ⋆₀⁻¹(X::Space)::Hom(DualForm1(X), Form0(X))
  ⋆₁⁻¹(X::Space)::Hom(DualForm0(X), Form1(X))
  ⋆₀(X) ⋅ ⋆₀⁻¹(X) == id(X) ⊣ (X::Ob)
  ⋆₁(X) ⋅ ⋆₁⁻¹(X) == id(X) ⊣ (X::Ob)

  δ₁(X::Space)::Hom(Form1(X), Form0(X))
  δ₁(X) == ⋆₁(X) ⋅ dual_d₀(X) ⋅ ⋆₀⁻¹(X) ⋅ antipode(Form0(X))

  ∇²₀(X::Space)::Hom(Form0(X),Form0(X))
  ∇²₀(X) == d₀(X) ⋅ δ₁(X) ⋅ antipode(Form0(X)) ⊣ (X::Space)

  Δ₀(X::Space)::Hom(Form0(X),Form0(X))
  Δ₁(X::Space)::Hom(Form1(X),Form1(X))
  Δ₀(X) == d₀(X) ⋅ δ₁(X) ⊣ (X::Space)
  Δ₁(X) == δ₁(X) ⋅ d₀(X) ⊣ (X::Space)
end

""" Theory of exterior calculus on 2D Riemannian manifold-like space.
"""
@theory ExtCalc2D{Ob,Hom,Space} <: MetricFreeExtCalc2D{Ob,Hom,Space} begin
  DualForm0(X::Space)::Ob
  DualForm1(X::Space)::Ob
  DualForm2(X::Space)::Ob
  DualForm0(X⊔Y) == DualForm0(X)⊕DualForm0(Y) ⊣ (X::Space, Y::Space)
  DualForm1(X⊔Y) == DualForm1(X)⊕DualForm1(Y) ⊣ (X::Space, Y::Space)
  DualForm2(X⊔Y) == DualForm2(X)⊕DualForm2(Y) ⊣ (X::Space, Y::Space)
  dual_d₀(X::Space)::Hom(DualForm0(X), DualForm1(X))
  dual_d₁(X::Space)::Hom(DualForm1(X), DualForm2(X))
  dual_d₀(X) ⋅ dual_d₁(X) == zero(DualForm0(X), DualForm2(X)) ⊣ (X::Space)

  ⋆₀(X::Space)::Hom(Form0(X), DualForm2(X))
  ⋆₁(X::Space)::Hom(Form1(X), DualForm1(X))
  ⋆₂(X::Space)::Hom(Form2(X), DualForm0(X))
  ⋆₀⁻¹(X::Space)::Hom(DualForm2(X), Form0(X))
  ⋆₁⁻¹(X::Space)::Hom(DualForm1(X), Form1(X))
  ⋆₂⁻¹(X::Space)::Hom(DualForm0(X), Form2(X))
  ⋆₀(X) ⋅ ⋆₀⁻¹(X) == id(Form0(X)) ⊣ (X::Ob)
  ⋆₁(X) ⋅ ⋆₁⁻¹(X) == antipode(Form1(X)) ⊣ (X::Ob)
  ⋆₂(X) ⋅ ⋆₂⁻¹(X) == id(Form1(X)) ⊣ (X::Ob)

  δ₁(X::Space)::Hom(Form1(X), Form0(X))
  δ₂(X::Space)::Hom(Form2(X), Form1(X))
  δ₁(X) == ⋆₁(X) ⋅ dual_d₁(X) ⋅ ⋆₀⁻¹(X) ⋅ antipode(Form0(X))
  δ₂(X) == ⋆₂(X) ⋅ dual_d₀(X) ⋅ ⋆₁⁻¹(X)

  ∇²₀(X::Space)::Hom(Form0(X),Form0(X))
  ∇²₁(X::Space)::Hom(Form1(X),Form1(X))
  ∇²₀(X) == d₀(X) ⋅ δ₁(X) ⋅ antipode(Form0(X)) ⊣ (X::Space)
  ∇²₁(X) == d₁(X) ⋅ δ₂(X) ⋅ antipode(Form1(X)) ⊣ (X::Space)

  Δ₀(X::Space)::Hom(Form0(X),Form0(X))
  Δ₁(X::Space)::Hom(Form1(X),Form1(X))
  Δ₂(X::Space)::Hom(Form2(X),Form2(X))
  Δ₀(X) == d₀(X) ⋅ δ₁(X) ⊣ (X::Space)
  Δ₁(X) == (d₁(X) ⋅ δ₂(X)) + (δ₁(X) ⋅ d₀(X)) ⊣ (X::Space)
  Δ₂(X) == δ₂(X) ⋅ d₁(X) ⊣ (X::Space)
end

@syntax FreeExtCalc2D{ObExpr,HomExpr,GATExpr} ExtCalc2D begin
  compose(f::Hom, g::Hom) = associate_unit(new(f,g; strict=true), id)
  oplus(A::Ob, B::Ob) = associate_unit(new(A,B), mzero)
  oplus(f::Hom, g::Hom) = associate(new(f,g))
  otimes(A::Ob, B::Ob) = associate_unit(new(A,B), munit)
  otimes(f::Hom, g::Hom) = associate(new(f,g))
end

end


text:
# This file was generated, do not modify it. # hide
pred_rfr = MLJ.predict(rfr_m, rows=test);
rms(pred_rfr, y[test])


text:
using Test
import LinearAlgebra.norm
include("../src/SSSTensor.jl")

 n = 3
 ord = 5
 tol = 1e-9
 nnz = 10

@testset "NumericalRoutines Tests" begin

    @testset "Low Mem SSHOPM Tests" begin
        #make tensors
        x = rand(n)
        indices = rand(1:n,nnz,ord)
        for i =1:nnz
            indices[i,:] = indices[i,sortperm(indices[i, :])]
        end
        vals = rand(nnz)
        dictA = ssten.SSSTensor([(indices[i,:],vals[i]) for i in 1:nnz])
        COOA = ssten.COOTen(indices,vals)

        max_iter = 500
        x_0 = randn(n)

        #no shift test
        shift =0.0
        COOHOPM_x, COOHOPM_λ, COOHOPM_iter =
          ssten.SSHOPM(COOA,x_0,shift,max_iter,tol)

        dictHOPM_x, dictHOPM_λ, dictHOPM_iter =
           ssten.SSHOPM(dictA, x_0,shift,max_iter,tol)

        @test dictHOPM_iter == COOHOPM_iter
        @test norm(COOHOPM_x - dictHOPM_x)/norm(COOHOPM_x) < tol
        @test abs(COOHOPM_λ - dictHOPM_λ) < tol

        #postive shift test
        shift = rand()

        COOHOPM_x, COOHOPM_λ, COOHOPM_iter =
          ssten.SSHOPM(COOA,x_0,shift,max_iter,tol)

        dictHOPM_x, dictHOPM_λ, dictHOPM_iter =
           ssten.SSHOPM(dictA, x_0,shift,max_iter,tol)

        @test dictHOPM_iter == COOHOPM_iter
        @test norm(COOHOPM_x - dictHOPM_x)/norm(COOHOPM_x) < tol
        @test abs(COOHOPM_λ - dictHOPM_λ) < tol

        #negative shift test
        shift = -rand()

        COOHOPM_x, COOHOPM_λ, COOHOPM_iter =
          ssten.SSHOPM(COOA,x_0,shift,max_iter,tol)

        dictHOPM_x, dictHOPM_λ, dictHOPM_iter =
           ssten.SSHOPM(dictA, x_0,shift,max_iter,tol)

        @test dictHOPM_iter == COOHOPM_iter
        @test norm(COOHOPM_x - dictHOPM_x)/norm(COOHOPM_x) < tol
        @test abs(COOHOPM_λ - dictHOPM_λ) < tol

    end

end


text:
using Project
using Base.Test

test_project = "test_project"

function test_new_project()
  Project.new(test_project,"0.0.1")
  return readall(open(string(pwd(),"/",test_project,"/project.json"))) != ""
end

function test_add_deps()
  cd(test_project)
  Project.add("MNIST")
  return get(Pkg.installed(),"MNIST",0) != 0
end

# write your own tests here
@test 1 == 1
@test test_new_project()
@test test_add_deps()

cd("..")
rm(test_project,recursive=true)


text:
abstract type FluidModel end

# NOTE: Should add a test for matching u0 with grid
mutable struct DiffusionModel{T<:Union{Function,Nothing}} <: FluidModel
    u0::Array{Float64,2}
    u_analytic::T

    grid::Grid
    parameters::Parameters

    Δt::Float64
    t::Array{Float64,1}

    function DiffusionModel(u0, u_analytic::T, grid, parameters) where T<:Union{Function,Nothing}
        Δt = parameters.σ*grid.Δx*grid.Δy/parameters.ν
        t = parameters.tspan[1]:Δt:parameters.tspan[2]
        return new{T}(u0, u_analytic, grid, parameters, Δt ,t)
    end
end
DiffusionModel(grid, parameters) = DiffusionModel(zeros(grid.ny, grid.nx), nothing, grid, parameters)
DiffusionModel(u0, grid, parameters) = DiffusionModel(u0, nothing, grid, parameters);


text:
"""
    forward_star(n_nodes::Int64, edge_list::Array{Int64, 2})

Creates the forward star representation of the input edge list
"""
function forward_star(n_nodes::Int64, edge_list::Array{Int64, 2})::Tuple{Array{Int64, 1}, Array{Int64, 1}}
    first_out = ones(Int64, n_nodes + 1)
    n_edges = size(edge_list, 1)
    end_node = zeros(Int64, n_edges)
    ind = 1
    for node in 1:n_nodes
        if ind > n_edges
            first_out[node] = ind
            continue
        end
        for i in ind:size(edge_list, 1)
            source = edge_list[i, 1]
            target = edge_list[i, 2]
            if source == node
                end_node[ind] = target
                ind += 1
            else
                first_out[node + 1] = ind
                break
            end
        end
    end
    first_out[end] = ind
    return first_out, end_node
end

# TODO: Implement reverse_star and integrate that with forward_star


text:
using NNlib: batched_mul

include("edgeftleayer.jl")
include("graphconv.jl")


text:
# Common functions for dynamics and rendering of pendulums
function pendulum_dynamics(env, s, a, rng::AbstractRNG = Random.GLOBAL_RNG)
    # Deal with terminal states
    # println("failur thresh: ", env.failure_thresh, " val: ", abs(s[1]))
    if (isnothing(env.failure_thresh) ?  false : abs(s[1]) > env.failure_thresh)
        # println("here")
        return fill(-100, size(s)), 0
    elseif env.include_time_in_state && (s[3] > env.maxT || s[3] ≈ env.maxT)
        return fill(100, size(s)), 0
    end

    θ, ω = s[1], s[2]
    dt, g, m, l = env.dt, env.g, env.m, env.l

    a = a[1]
    a = clamp(a, -env.max_torque, env.max_torque)
    costs = angle_normalize(θ)^2 + 0.1f0 * ω^2 + 0.001f0 * a^2

    ω = ω + (-3. * g / (2 * l) * sin(θ + π) + 3. * a / (m * l^2)) * dt
    θ = angle_normalize(θ + ω * dt)
    ω = clamp(ω, -env.max_speed, env.max_speed)

    if env.include_time_in_state
        sp = Float32.([θ, ω, s[3] + dt])
    else
        sp = Float32.([θ, ω])
    end
    r = Float32(env.Rstep - env.λcost*costs)
    return sp, r
end

angle_normalize(x) = mod((x+π), (2*π)) - π

function render_pendulum(env, s, a)
    θ = s[1] + π/2.
    point_array = [(0.5,0.5), (0.5 + 0.3*cos(θ), 0.5 - 0.3*sin(θ))]

    a_rad = abs(a)/10.
    if a < 0
        θstart = -3π/4
        θend = -θstart
        θarr = θend
    else
        θstart = π/4
        θend = -θstart
        θarr = θstart
    end

    # Draw the arrow
    endpt = (0.5, 0.5) .+ a_rad.*(cos(θarr), sin(θarr))
    uparr = endpt .+ 0.1*a_rad.*(cos(θarr)-sign(a)*sin(θarr), sign(a)*cos(θarr)+sin(θarr))
    dwnarr = endpt .+ 0.1*a_rad.*(-cos(θarr)-sign(a)sin(θarr), sign(a)*cos(θarr)-sin(θarr))
    arrow_array = [[endpt, uparr], [endpt, dwnarr]]

    img = compose(context(),
        (context(), line(arrow_array), arc(0.5, 0.5, a_rad, θstart, θend),  linewidth(0.5mm), fillopacity(0.), stroke("red")),
        (context(), circle(0.5, 0.5, 0.01), fill("blue"), stroke("black")),
        (context(), line(point_array), stroke("black"), linewidth(1mm)),
        (context(), rectangle(), fill("white"))
    )
    tmpfilename = tempname()
    img |> PNG(tmpfilename, 10cm, 10cm)
    load(tmpfilename)
end

## POMDP implementation -- Allows for pixel observation
@with_kw struct PendulumPOMDP <: POMDP{Array{Float32}, Float64, Array{Float32}}
    failure_thresh::Union{Nothing, Float64} = nothing # if set, defines the operating range fo the pendulum. Episode terminates if abs(θ) is larger than this.
    θ0 = Distributions.Uniform(-π, π) # Distribution to sample initial angular position
    ω0 = Distributions.Uniform(-1., 1.) # Distribution to sample initial angular velocity
    Rstep = 0 # Reward earned on each step of the simulation
    λcost = 1 # Coefficient to the traditional OpenAIGym Reward
    max_speed::Float64 = 8.
    max_torque::Float64 = 2.
    dt::Float64 = .05
    g::Float64 = 10.
    m::Float64 = 1.
    l::Float64 = 1.
    γ::Float32 = 0.99
    actions::Vector{Float64} = [-1., 1.]
    pixel_observations::Bool = false
    render_fun::Union{Nothing, Function} = nothing
    include_time_in_state = false
    maxT = 99*dt
end

InvertedPendulumPOMDP(failure_thresh = deg2rad(20),
                    θ0 = Distributions.Uniform(-failure_thresh/2., failure_thresh/2.),
                    ω0 = Distributions.Uniform(-.1, .1),
                    Rstep = 1,
                    λcost = 1;
                    kwargs...) = PendulumPOMDP(failure_thresh = failure_thresh, θ0 = θ0, ω0 = ω0, Rstep = Rstep, λcost = λcost; kwargs...)

function POMDPs.gen(mdp::PendulumPOMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG)
    sp, r = pendulum_dynamics(mdp, s, a, rng)
    (sp = sp, o=rand(rng, observation(mdp, sp)), r = r)
end

function POMDPs.observation(mdp::PendulumPOMDP, s)
    o = mdp.pixel_observations ? mdp.render_fun(s) : [angle_normalize(s[1]), s[2]] #[cos(s[1]), sin(s[1]), s[2]]
    Deterministic(Float32.(o))
end

function POMDPs.initialstate(mdp::PendulumPOMDP)
    ImplicitDistribution((rng) -> Float32.([rand(rng, mdp.θ0), rand(rng, mdp.ω0)]))
end

POMDPs.initialobs(mdp::PendulumPOMDP, s) = observation(mdp, s)

POMDPs.actions(mdp::PendulumPOMDP) = mdp.actions
isfailure(mdp::PendulumPOMDP, s) = s[1] == -100
# Terminate if there is a failure or if we exceed the maximum time (in the case that we are tracking such things)
POMDPs.isterminal(mdp::PendulumPOMDP, s) = isfailure(mdp, s) || s[1] == 100
POMDPs.discount(mdp::PendulumPOMDP) = mdp.γ

render(mdp::PendulumPOMDP, s, a::AbstractArray) = render(mdp, s, a...)
render(mdp::PendulumPOMDP, s, a = 0) = render_pendulum(mdp, s, a)

## MDP formulation
@with_kw struct PendulumMDP <: MDP{Array{Float32}, Array{Float32}}
    failure_thresh::Union{Nothing, Float64} = nothing # if set, defines the operating range fo the pendulum. Episode terminates if abs(θ) is larger than this.
    θ0 = Distributions.Uniform(-π, π) # Distribution to sample initial angular position
    ω0 = Distributions.Uniform(-1., 1.) # Distribution to sample initial angular velocity
    Rstep = 0 # Reward earned on each step of the simulation
    λcost = 1 # Coefficient to the traditional OpenAIGym Reward
    max_speed::Float64 = 8.
    max_torque::Float64 = 2.
    dt::Float64 = .05
    g::Float64 = 10.
    m::Float64 = 1.
    l::Float64 = 1.
    γ::Float32 = 0.99
    actions::Vector{Float64} = [-1., 1.]
    include_time_in_state = false
    maxT = 99*dt
end

InvertedPendulumMDP(failure_thresh = deg2rad(20),
                    θ0 = Distributions.Uniform(-failure_thresh/2., failure_thresh/2.),
                    ω0 = Distributions.Uniform(-.1, .1),
                    Rstep = 1,
                    λcost = 1;
                    kwargs...) = PendulumMDP(failure_thresh = failure_thresh, θ0 = θ0, ω0 = ω0, Rstep = Rstep, λcost = λcost; kwargs...)

function POMDPs.gen(mdp::PendulumMDP, s, a, rng::AbstractRNG = Random.GLOBAL_RNG)
    sp, r = pendulum_dynamics(mdp, s, a, rng)
    (sp = sp, r = r)
end

function POMDPs.initialstate(mdp::PendulumMDP)
    if mdp.include_time_in_state
        ImplicitDistribution((rng) -> Float32.([rand(rng, mdp.θ0), rand(rng, mdp.ω0), 0f0]))
    else
        ImplicitDistribution((rng) -> Float32.([rand(rng, mdp.θ0), rand(rng, mdp.ω0)]))
    end
end

POMDPs.actions(mdp::PendulumMDP) = mdp.actions

isfailure(mdp::PendulumMDP, s) = s[1] == -100
# Terminate if there is a failure or if we exceed the maximum time (in the case that we are tracking such things)
POMDPs.isterminal(mdp::PendulumMDP, s) = isfailure(mdp, s) || s[1] == 100
POMDPs.discount(mdp::PendulumMDP) = mdp.γ

render(mdp::PendulumMDP, s, a::AbstractArray) = render(mdp, s, a...)
render(mdp::PendulumMDP, s, a = 0) = render_pendulum(mdp, s, a)

max_forks_repo_name: ancorso/POMDPGym


text:
abstract AbstractCmd

type Cmd <: AbstractCmd
    exec::Vector{ByteString}
    ignorestatus::Bool
    detach::Bool
    env::Union(Array{ByteString},Nothing)
    dir::UTF8String
    Cmd(exec::Vector{ByteString}) = new(exec,false,false,nothing,"")
end

type OrCmds <: AbstractCmd
    a::AbstractCmd
    b::AbstractCmd
    OrCmds(a::AbstractCmd, b::AbstractCmd) = new(a,b)
end

type ErrOrCmds <: AbstractCmd
    a::AbstractCmd
    b::AbstractCmd
    ErrOrCmds(a::AbstractCmd, b::AbstractCmd) = new(a,b)
end

type AndCmds <: AbstractCmd
    a::AbstractCmd
    b::AbstractCmd
    AndCmds(a::AbstractCmd, b::AbstractCmd) = new(a,b)
end

shell_escape(cmd::Cmd) = shell_escape(cmd.exec...)

function show(io::IO, cmd::Cmd)
    print_env = cmd.env !== nothing
    print_dir = !isempty(cmd.dir)
    (print_env || print_dir) && print(io,"setenv(")
    esc = shell_escape(cmd)
    print(io,'`')
    for c in esc
        if c == '`'
            print(io,'\\')
        end
        print(io,c)
    end
    print(io,'`')
    print_env && (print(io,","); show(io,cmd.env))
    print_dir && (print(io,"; dir="); show(io,cmd.dir))
    (print_dir || print_env) && print(io,")")
end

function show(io::IO, cmds::OrCmds)
    if isa(cmds.a, AndCmds) || isa(cmds.a, CmdRedirect)
        print(io,"(")
        show(io, cmds.a)
        print(io,")")
    else
        show(io, cmds.a)
    end
    print(io, " |> ")
    if isa(cmds.b, AndCmds) || isa(cmds.b, CmdRedirect)
        print(io,"(")
        show(io, cmds.b)
        print(io,")")
    else
        show(io, cmds.b)
    end
end

function show(io::IO, cmds::AndCmds)
    if isa(cmds.a, OrCmds) || isa(cmds.a, CmdRedirect)
        print(io,"(")
        show(io, cmds.a)
        print(io,")")
    else
        show(io, cmds.a)
    end
    print(io," & ")
    if isa(cmds.b, OrCmds) || isa(cmds.b, CmdRedirect)
        print(io,"(")
        show(io, cmds.b)
        print(io,")")
    else
        show(io, cmds.b)
    end
end

const STDIN_NO  = 0
const STDOUT_NO = 1
const STDERR_NO = 2

immutable FileRedirect
    filename::String
    append::Bool
    function FileRedirect(filename,append)
        if lowercase(filename) == (@unix? "/dev/null" : "nul")
            warn_once("for portability use DevNull instead of a file redirect")
        end
        new(filename,append)
    end
end

immutable DevNullStream <: AsyncStream end
const DevNull = DevNullStream()
copy(::DevNullStream) = DevNull
uvhandle(::DevNullStream) = C_NULL
uvhandle(x::Ptr) = x
uvtype(::Ptr) = UV_STREAM
uvtype(::DevNullStream) = UV_STREAM

# Not actually a pointer, but that's how we pass it through the C API
# so it's fine
uvhandle(x::RawFD) = convert(Ptr{Void},x.fd)
uvtype(x::RawFD) = UV_RAW_FD

typealias Redirectable Union(UVStream,FS.File,FileRedirect,DevNullStream,IOStream,RawFD)

type CmdRedirect <: AbstractCmd
    cmd::AbstractCmd
    handle::Redirectable
    stream_no::Int
end

function show(io::IO, cr::CmdRedirect)
    if cr.stream_no == STDOUT_NO
        show(io,cr.cmd)
        print(io," |> ")
        show(io,cr.handle)
    elseif cr.stream_no == STDERR_NO
        show(io,cr.cmd)
        print(io," .> ")
        show(io,cr.handle)
    elseif cr.stream_no == STDIN_NO
        show(io,cr.handle)
        print(io," |> ")
        show(io,cr.cmd)
    end
end

ignorestatus(cmd::Cmd) = (cmd.ignorestatus=true; cmd)
ignorestatus(cmd::Union(OrCmds,AndCmds)) = (ignorestatus(cmd.a); ignorestatus(cmd.b); cmd)
detach(cmd::Cmd) = (cmd.detach=true; cmd)

setenv{S<:ByteString}(cmd::Cmd, env::Array{S}; dir="") = (cmd.env = ByteString[x for x in env]; setenv(cmd,dir=dir); cmd)
setenv(cmd::Cmd, env::Associative; dir="") = (cmd.env = ByteString[string(k)*"="*string(v) for (k,v) in env]; setenv(cmd,dir=dir); cmd)
setenv(cmd::Cmd; dir="") = (cmd.dir = dir; cmd)

(&)(left::AbstractCmd,right::AbstractCmd) = AndCmds(left,right)
(|>)(src::AbstractCmd,dest::AbstractCmd) = OrCmds(src,dest)
(.>)(src::AbstractCmd,dest::AbstractCmd) = ErrOrCmds(src,dest)

# Stream Redirects
(|>)(dest::Redirectable,src::AbstractCmd) = CmdRedirect(src,dest,STDIN_NO)
(|>)(src::AbstractCmd,dest::Redirectable) = CmdRedirect(src,dest,STDOUT_NO)
(.>)(src::AbstractCmd,dest::Redirectable) = CmdRedirect(src,dest,STDERR_NO)

# File redirects
(|>)(src::AbstractCmd,dest::String) = CmdRedirect(src,FileRedirect(dest,false),STDOUT_NO)
(|>)(src::String,dest::AbstractCmd) = CmdRedirect(dest,FileRedirect(src,false),STDIN_NO)
(.>)(src::AbstractCmd,dest::String) = CmdRedirect(src,FileRedirect(dest,false),STDERR_NO)
(>>)(src::AbstractCmd,dest::String) = CmdRedirect(src,FileRedirect(dest,true),STDOUT_NO)
(.>>)(src::AbstractCmd,dest::String) = CmdRedirect(src,FileRedirect(dest,true),STDERR_NO)

typealias RawOrBoxedHandle Union(UVHandle,UVStream,Redirectable,IOStream)
typealias StdIOSet NTuple{3,RawOrBoxedHandle}

type Process
    cmd::Cmd
    handle::Ptr{Void}
    in::AsyncStream
    out::AsyncStream
    err::AsyncStream
    exitcode::Int32
    termsignal::Int32
    exitcb::Callback
    exitnotify::Condition
    closecb::Callback
    closenotify::Condition
    function Process(cmd::Cmd,handle::Ptr{Void},in::RawOrBoxedHandle,out::RawOrBoxedHandle,err::RawOrBoxedHandle)
        if !isa(in,AsyncStream) || in === DevNull
            in=DevNull
        end
        if !isa(out,AsyncStream) || out === DevNull
            out=DevNull
        end
        if !isa(err,AsyncStream) || err === DevNull
            err=DevNull
        end
        this = new(cmd,handle,in,out,err,typemin(Int32),typemin(Int32),false,Condition(),false,Condition())
        finalizer(this, uvfinalize)
        this
    end
end

type ProcessChain
    processes::Vector{Process}
    in::Redirectable
    out::Redirectable
    err::Redirectable
    ProcessChain(stdios::StdIOSet) = new(Process[],stdios[1],stdios[2],stdios[3])
end
typealias ProcessChainOrNot Union(Bool,ProcessChain)

function _jl_spawn(cmd::Ptr{Uint8}, argv::Ptr{Ptr{Uint8}}, loop::Ptr{Void}, pp::Process,
                   in, out, err)
    proc = c_malloc(_sizeof_uv_process)
    error = ccall(:jl_spawn, Int32,
        (Ptr{Uint8}, Ptr{Ptr{Uint8}}, Ptr{Void}, Ptr{Void}, Any, Int32,
         Ptr{Void},    Int32,       Ptr{Void},     Int32,       Ptr{Void},
         Int32, Ptr{Ptr{Uint8}}, Ptr{Uint8}),
         cmd,        argv,            loop,      proc,      pp,  uvtype(in),
         uvhandle(in), uvtype(out), uvhandle(out), uvtype(err), uvhandle(err),
         pp.cmd.detach, pp.cmd.env === nothing ? C_NULL : pp.cmd.env, isempty(pp.cmd.dir) ? C_NULL : pp.cmd.dir)
    if error != 0
        disassociate_julia_struct(proc)
        ccall(:jl_forceclose_uv,Void,(Ptr{Void},),proc)
        throw(UVError("could not spawn "*string(pp.cmd), error))
    end
    associate_julia_struct(proc,pp)
    return proc
end

function uvfinalize(proc::Process)
    proc.handle != C_NULL && ccall(:jl_close_uv,Void,(Ptr{Void},),proc.handle)
    disassociate_julia_struct(proc)
    proc.handle = 0
end

function _uv_hook_return_spawn(proc::Process, exit_status::Int64, termsignal::Int32)
    proc.exitcode = exit_status
    proc.termsignal = termsignal
    if isa(proc.exitcb, Function) proc.exitcb(proc, exit_status, termsignal) end
    ccall(:jl_close_uv,Void,(Ptr{Void},),proc.handle)
    notify(proc.exitnotify)
end

function _uv_hook_close(proc::Process)
    proc.handle = 0
    if isa(proc.closecb, Function) proc.closecb(proc) end
    notify(proc.closenotify)
end

function spawn(pc::ProcessChainOrNot,redirect::CmdRedirect,stdios::StdIOSet,exitcb::Callback,closecb::Callback)
        spawn(pc,redirect.cmd,(redirect.stream_no==STDIN_NO ?redirect.handle:stdios[1],
                                                   redirect.stream_no==STDOUT_NO?redirect.handle:stdios[2],
                                                   redirect.stream_no==STDERR_NO?redirect.handle:stdios[3]),exitcb,closecb)
end

function spawn(pc::ProcessChainOrNot,cmds::OrCmds,stdios::StdIOSet,exitcb::Callback,closecb::Callback)
    out_pipe = box(Ptr{Void},Intrinsics.jl_alloca(_sizeof_uv_named_pipe))
    in_pipe = box(Ptr{Void},Intrinsics.jl_alloca(_sizeof_uv_named_pipe))
    #out_pipe = c_malloc(_sizeof_uv_named_pipe)
    #in_pipe = c_malloc(_sizeof_uv_named_pipe)
    link_pipe(in_pipe,false,out_pipe,false)
    if pc == false
        pc = ProcessChain(stdios)
    end
    try
        spawn(pc, cmds.a, (stdios[1], out_pipe, stdios[3]), exitcb, closecb)
        spawn(pc, cmds.b, (in_pipe, stdios[2], stdios[3]), exitcb, closecb)
    catch err
        close_pipe_sync(out_pipe)
        close_pipe_sync(in_pipe)
        rethrow(err)
    end
    close_pipe_sync(out_pipe)
    close_pipe_sync(in_pipe)
    pc
end

function spawn(pc::ProcessChainOrNot,cmds::ErrOrCmds,stdios::StdIOSet,exitcb::Callback,closecb::Callback)
    out_pipe = box(Ptr{Void},Intrinsics.jl_alloca(_sizeof_uv_named_pipe))
    in_pipe = box(Ptr{Void},Intrinsics.jl_alloca(_sizeof_uv_named_pipe))
    #out_pipe = c_malloc(_sizeof_uv_named_pipe)
    #in_pipe = c_malloc(_sizeof_uv_named_pipe)
    link_pipe(in_pipe,false,out_pipe,false)
    if pc == false
        pc = ProcessChain(stdios)
    end
    try
        spawn(pc, cmds.a, (stdios[1], stdios[2], out_pipe), exitcb, closecb)
        spawn(pc, cmds.b, (in_pipe, stdios[2], stdios[3]), exitcb, closecb)
    catch err
        close_pipe_sync(out_pipe)
        close_pipe_sync(in_pipe)
        rethrow(err)
    end
    close_pipe_sync(out_pipe)
    close_pipe_sync(in_pipe)
    pc
end

macro setup_stdio()
    esc(
    quote
        close_in,close_out,close_err = false,false,false
        in,out,err = stdios
        if isa(stdios[1],Pipe)
            if stdios[1].handle==C_NULL
                error("pipes passed to spawn must be initialized")
            end
        elseif isa(stdios[1],FileRedirect)
            in = FS.open(stdios[1].filename,JL_O_RDONLY)
            close_in = true
        elseif isa(stdios[1],IOStream)
            in = FS.File(RawFD(fd(stdios[1])))
        end
        if isa(stdios[2],Pipe)
            if stdios[2].handle==C_NULL
                error("pipes passed to spawn must be initialized")
            end
        elseif isa(stdios[2],FileRedirect)
            out = FS.open(stdios[2].filename,JL_O_WRONLY|JL_O_CREAT|(stdios[2].append?JL_O_APPEND:JL_O_TRUNC),S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH)
            close_out = true
        elseif isa(stdios[2],IOStream)
            out = FS.File(RawFD(fd(stdios[2])))
        end
        if isa(stdios[3],Pipe)
            if stdios[3].handle==C_NULL
                error("pipes passed to spawn must be initialized")
            end
        elseif isa(stdios[3],FileRedirect)
            err = FS.open(stdios[3].filename,JL_O_WRONLY|JL_O_CREAT|(stdios[3].append?JL_O_APPEND:JL_O_TRUNC),S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH)
            close_err = true
        elseif isa(stdios[3],IOStream)
            err = FS.File(RawFD(fd(stdios[3])))
        end
    end)
end

macro cleanup_stdio()
    esc(
    quote
        close_in && close(in)
        close_out && close(out)
        close_err && close(err)
    end)
end

function spawn(pc::ProcessChainOrNot,cmd::Cmd,stdios::StdIOSet,exitcb::Callback,closecb::Callback)
    loop = eventloop()
    pp = Process(cmd,C_NULL,stdios[1],stdios[2],stdios[3]);
    @setup_stdio
    ptrs = _jl_pre_exec(cmd.exec)
    pp.exitcb = exitcb
    pp.closecb = closecb
    pp.handle = _jl_spawn(ptrs[1], convert(Ptr{Ptr{Uint8}}, ptrs), loop, pp,
                          in,out,err)
    @cleanup_stdio
    if isa(pc, ProcessChain)
        push!(pc.processes,pp)
    end
    pp
end

function spawn(pc::ProcessChainOrNot,cmds::AndCmds,stdios::StdIOSet,exitcb::Callback,closecb::Callback)
    if pc == false
        pc = ProcessChain(stdios)
    end
    @setup_stdio
    spawn(pc, cmds.a, (in,out,err), exitcb, closecb)
    spawn(pc, cmds.b, (in,out,err), exitcb, closecb)
    @cleanup_stdio
    pc
end

# INTERNAL
# returns a tuple of function arguments to spawn:
# (stdios, exitcb, closecb)
# |       |        \ The function to be called once the uv handle is closed
# |       \ The function to be called once the process exits
# \ A set of up to 256 stdio instructions, where each entry can be either:
#   | - An AsyncStream to be passed to the child
#   | - DevNull to pass /dev/null
#   | - An FS.File object to redirect the output to
#   \ - An ASCIIString specifying a filename to be opened

spawn_opts_swallow(stdios::StdIOSet,exitcb::Callback=false,closecb::Callback=false) =
    (stdios,exitcb,closecb)
spawn_opts_swallow(in::Redirectable=DevNull,out::Redirectable=DevNull,err::Redirectable=DevNull,args...) =
    (tuple(in,out,err,args...),false,false)
spawn_opts_inherit(stdios::StdIOSet,exitcb::Callback=false,closecb::Callback=false) =
    (stdios,exitcb,closecb)
spawn_opts_inherit(in::Redirectable=STDIN,out::Redirectable=STDOUT,err::Redirectable=STDERR,args...) =
    (tuple(in,out,err,args...),false,false)

spawn(pc::ProcessChainOrNot,cmds::AbstractCmd,args...) = spawn(pc,cmds,spawn_opts_swallow(args...)...)
spawn(cmds::AbstractCmd,args...) = spawn(false,cmds,spawn_opts_swallow(args...)...)

macro tmp_rpipe(pipe,tmppipe,code,args...)
    esc(quote
        $pipe = Pipe(C_NULL)
        $tmppipe = Pipe(C_NULL)
        link_pipe($pipe,true,$tmppipe,false)
        r = begin
            $code
        end
        close_pipe_sync($tmppipe)
        r
    end)
end

macro tmp_wpipe(tmppipe,pipe,code)
    esc(quote
        $pipe = Pipe(C_NULL)
        $tmppipe = Pipe(C_NULL)
        link_pipe($tmppipe,false,$pipe,true)
        r = begin
            $code
        end
        close_pipe_sync($tmppipe)
        r
    end)
end

function eachline(cmd::AbstractCmd,stdin)
    @tmp_rpipe out tmp begin
        processes = spawn(false, cmd, (stdin,tmp,STDERR))
        # implicitly close after reading lines, since we opened
        EachLine(out, ()->close(out))
    end
end
eachline(cmd::AbstractCmd) = eachline(cmd,DevNull)

# return a (Pipe,Process) pair to write/read to/from the pipeline
function open(cmds::AbstractCmd, mode::String="r", stdio::AsyncStream=DevNull)
    if mode == "r"
        processes = @tmp_rpipe out tmp spawn(false, cmds, (stdio,tmp,STDERR))
        start_reading(out)
        (out, processes)
    elseif mode == "w"
        processes = @tmp_wpipe tmp inpipe spawn(false,cmds, (tmp,stdio,STDERR))
        (inpipe, processes)
    else
        throw(ArgumentError("mode must be \"r\" or \"w\", not \"$mode\""))
    end
end

function open(f::Function, cmds::AbstractCmd, args...)
    io, P = open(cmds, args...)
    ret = try
        f(io)
    catch
        kill(P)
        rethrow()
    finally
        close(io)
    end
    success(P) || pipeline_error(P)
    return ret
end

# TODO: convert this to use open(cmd, "r+"), with a single read/write pipe
function readandwrite(cmds::AbstractCmd)
    (out, processes) = @tmp_wpipe tmp inpipe open(cmds, "r", tmp)
    (out, inpipe, processes)
end

function readbytes(cmd::AbstractCmd, stdin::AsyncStream=DevNull)
    (out,pc) = open(cmd, "r", stdin)
    !success(pc) && pipeline_error(pc)
    wait_close(out)
    return takebuf_array(out.buffer)
end

function readall(cmd::AbstractCmd, stdin::AsyncStream=DevNull)
    return bytestring(readbytes(cmd, stdin))
end

function writeall(cmd::AbstractCmd, stdin::String, stdout::AsyncStream=DevNull)
    open(cmd, "w", stdout) do io
        write(io, stdin)
    end
end

function run(cmds::AbstractCmd,args...)
    ps = spawn(cmds,spawn_opts_inherit(args...)...)
    success(ps) ? nothing : pipeline_error(ps)
end

const SIGPIPE = 13
function test_success(proc::Process)
    assert(process_exited(proc))
    if proc.exitcode < 0
        error(UVError("could not start process "*string(proc.cmd),proc.exitcode))
    end
    proc.exitcode==0 && (proc.termsignal == 0 || proc.termsignal == SIGPIPE)
end

function success(x::Process)
    wait(x)
    kill(x)
    test_success(x)
end
success(procs::Vector{Process}) = mapreduce(success, &, procs)
success(procs::ProcessChain) = success(procs.processes)
success(cmd::AbstractCmd) = success(spawn(cmd))

function pipeline_error(proc::Process)
    if !proc.cmd.ignorestatus
        error("failed process: ",proc," [",proc.exitcode,"]")
    end
    nothing
end

function pipeline_error(procs::ProcessChain)
    failed = Process[]
    for p = procs.processes
        if !test_success(p) && !p.cmd.ignorestatus
            push!(failed, p)
        end
    end
    if length(failed)==0 return nothing end
    if length(failed)==1 pipeline_error(failed[1]) end
    msg = "failed processes:"
    for proc in failed
        msg = string(msg,"\n  ",proc," [",proc.exitcode,"]")
    end
    error(msg)
end

_jl_kill(p::Process,signum::Integer) = ccall(:uv_process_kill,Int32,(Ptr{Void},Int32),p.handle,signum)
function kill(p::Process,signum::Integer)
    if process_running(p)
        @assert p.handle != C_NULL
        _jl_kill(p, signum)
    else
        int32(-1)
    end
end
kill(ps::Vector{Process}) = map(kill, ps)
kill(ps::ProcessChain) = map(kill, ps.processes)
kill(p::Process) = kill(p,15) #SIGTERM

function _contains_newline(bufptr::Ptr{Void},len::Int32)
    return (ccall(:memchr,Ptr{Void},(Ptr{Void},Int32,Csize_t),bufptr,'\n',len)!=C_NULL)
end

## process status ##
process_running(s::Process) = s.exitcode == typemin(Int32)
process_running(s::Vector{Process}) = any(process_running,s)
process_running(s::ProcessChain) = process_running(s.processes)

process_exited(s::Process) = !process_running(s)
process_exited(s::Vector{Process}) = all(process_exited,s)
process_exited(s::ProcessChain) = process_exited(s.processes)

process_signaled(s::Process) = (s.termsignal > 0)

#process_stopped (s::Process) = false #not supported by libuv. Do we need this?
#process_stop_signal(s::Process) = false #not supported by libuv. Do we need this?

function process_status(s::Process)
    process_running (s) ? "ProcessRunning" :
    process_signaled(s) ? "ProcessSignaled("*string(s.termsignal)*")" :
    #process_stopped (s) ? "ProcessStopped("*string(process_stop_signal(s))*")" :
    process_exited  (s) ? "ProcessExited("*string(s.exitcode)*")" :
    error("process status error")
end

# WARNING: do not call this and keep the returned array of pointers
# around longer than the args vector and then use array of pointers.
# this could cause a segfault. this is really just for use by the
# spawn function below so that we can exec more efficiently.
#
function _jl_pre_exec(args::Vector{ByteString})
    if length(args) < 1
        error("too few arguments to exec")
    end
    ptrs = Array(Ptr{Uint8}, length(args)+1)
    for i = 1:length(args)
        ptrs[i] = args[i].data
    end
    ptrs[length(args)+1] = C_NULL
    return ptrs
end

## implementation of `cmd` syntax ##

arg_gen()          = ByteString[]
arg_gen(x::String) = ByteString[x]
arg_gen(cmd::Cmd)  = cmd.exec

function arg_gen(head)
    if applicable(start,head)
        vals = ByteString[]
        for x in head
            push!(vals,string(x))
        end
        return vals
    else
        return ByteString[string(head)]
    end
end

function arg_gen(head, tail...)
    head = arg_gen(head)
    tail = arg_gen(tail...)
    vals = ByteString[]
    for h = head, t = tail
        push!(vals,bytestring(h,t))
    end
    vals
end

function cmd_gen(parsed)
    args = ByteString[]
    for arg in parsed
        append!(args,arg_gen(arg...))
    end
    Cmd(args)
end

macro cmd(str)
    :(cmd_gen($(shell_parse(str)[1])))
end

wait(x::Process)      = if !process_exited(x); stream_wait(x,x.exitnotify); end
wait(x::ProcessChain) = for p in x.processes; wait(p); end

show(io::IO, p::Process) = print(io, "Process(", p.cmd, ", ", process_status(p), ")")

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
include("layout.jl")                    # layout
include(joinpath("api", "endpoint.jl")) # Endpoint
include("wasm.jl")                      # WASM
include("front.jl")                     # Front

if PROGRAM_FILE == basename(@__FILE__)

using Bukdu
import Bukdu.Actions: index, show, new, edit, create, delete, update
using .Front # WelcomeController
using .Endpoint # CustomerController
import .WASM: WasmController, hello_js, hello_wast
using Sockets

pipeline(:api) do conn::Conn
end

routes(:api) do
    resources("/customers", CustomerController)
end

routes(:wasm) do
    get("/wasm", WasmController, index)
    get("/hello.js", WasmController, hello_js)
    get("/hello.wast", WasmController, hello_wast)
    plug(Plug.Static, at="/", from=normpath(@__DIR__, "public"))
end

routes(:front) do
    get("/", WelcomeController, index)
end

if haskey(ENV, "ON_HEROKU")
    Bukdu.start(parse(Int, ENV["PORT"]); host=Sockets.IPAddr(0,0,0,0))
else
    Bukdu.start(8080)
end

Router.call(get, "/") #
# CLI.routes()

Base.JLOptions().isinteractive==0 && wait()

# Bukdu.stop()

end # if

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu