text:
struct SegmentedMax{T} <: AggregationFunction
    C::T
end

Flux.@functor SegmentedMax

SegmentedMax(d::Int) = SegmentedMax(zeros(Float32, d))

Base.show(io::IO, sm::SegmentedMax) = print(io, "SegmentedMax($(length(sm.C)))\n")
modelprint(io::IO, sm::SegmentedMax; pad=[]) = paddedprint(io, "SegmentedMax($(length(sm.C)))")

(m::SegmentedMax)(x::MaybeMatrix, bags::AbstractBags, w=nothing) = segmented_max_forw(x, m.C, bags)
function (m::SegmentedMax)(x::AbstractMatrix, bags::AbstractBags, w::AggregationWeights, mask::AbstractVector)
    segmented_max_forw(x .+ typemin(T) * mask', m.C, bags)
end

segmented_max_forw(::Missing, C::AbstractVector, bags::AbstractBags) = repeat(C, 1, length(bags))
function segmented_max_forw(x::AbstractMatrix, C::AbstractVector, bags::AbstractBags)
    y = fill(typemin(eltype(x)), size(x, 1), length(bags))
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(C)
                y[i, bi] = C[i]
            end
        else
            for j in b
                for i in 1:size(x, 1)
                    y[i, bi] = max(y[i, bi], x[i, j])
                end
            end
        end
    end
    y
end

function segmented_max_back(Δ, y, x, C, bags)
    dx = zero(x)
    dC = zero(C)
    v = similar(x, size(x, 1))
    idxs = zeros(Int, size(x, 1))
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(C)
                dC[i] += Δ[i, bi]
            end
        else
            fi = first(b)
            v .= x[:,fi]
            idxs .= fi
            for j in b
                for i in 1:size(x,1)
                    if v[i] < x[i, j]
                        idxs[i] = j
                        v[i] = x[i, j]
                    end
                end
            end
            for i in 1:size(x, 1)
                dx[i, idxs[i]] = Δ[i, bi]
            end
        end
    end
    dx, dC, nothing, nothing
end

function segmented_max_back(Δ, y, x::Missing, C, bags)
    dC = zero(C)
    @inbounds for (bi, b) in enumerate(bags)
        for i in eachindex(C)
            dC[i] += Δ[i, bi]
        end
    end
    nothing, dC, nothing, nothing
end

@adjoint function segmented_max_forw(args...)
    y = segmented_max_forw(args...)
    grad = Δ -> segmented_max_back(Δ, y, args...)
    y, grad
end


text:
export reconstructionPeriodicMotion

function reconstructionPeriodicMotion(bSF::MPIFile, bMeas::MPIFile;
  minFreq=0, maxFreq=1.25e6, SNRThresh=-1,maxMixingOrder=-1, numUsedFreqs=-1, sortBySNR=false, recChannels=1:numReceivers(bMeas), kargs...)

  freq = filterFrequencies(bSF,minFreq=minFreq, maxFreq=maxFreq,recChannels=recChannels, SNRThresh=SNRThresh, numUsedFreqs=numUsedFreqs, sortBySNR=sortBySNR)

  @debug "selecting $(length(freq)) frequencies"

  return reconstructionPeriodicMotion(bSF, bMeas, freq; kargs...)
end

"""
    reconstructionPeriodicMotion(bSF::MPIFile, bMeas::MPIFile, freq::Array{Int64,1};
				bEmpty=nothing, bgFrames=nothing,
				alpha::Float64=3.0, choosePeak::Int64=1, frames::UnitRange=1:acqNumFrames(bMeas),
				samplingPrecision::Bool=true, windowType::Int64=1,
				bSFFrequencyAnalysis::MPIFile=bSF,higherHarmonic::Int64=1,
				kargs...)

	Performs multi-patch reconstruction of raw data from an object with periodic motion

- bSF:			System functions for reconstruction
- bMeas:		Raw data of the measurement
- freq:                 Selected frequencies for reconstruction
- bEmpty:		Background measurement
- bgFrames:		Background frames
- choosePeak:		Number of chosen peak for motion frequency
- alpha:        	Window width relative to DF cycle
- frames: 		Selected frame
- lambda:		Regularization parameter for reconstruction
- iterations:		Number of iterations
- samplingPrecision:    true: rounding motion period to sampling precision, false: rounding to DF period precision
- windowFunction:       1: Hann window, 2:FT1A05, 3: Rectangle
- bSFFrequencyAnalysis: System function for frequency analysis

"""
function reconstructionPeriodicMotion(bSF::MPIFile, bMeas::MPIFile, freq::Array{Int64,1};
				bEmpty=nothing, bgFrames=nothing,
				alpha::Float64=3.0, choosePeak::Int64=1,
		                frames::UnitRange=1:acqNumFrames(bMeas),
				samplingPrecision::Bool=true, windowType::Int64=1,
				bSFFrequencyAnalysis::MPIFile=bSF, higherHarmonic::Int64=1,
				kargs...)

  FFP = squeeze(acqOffsetFieldShift(bMeas))

  motFreq = getMotionFreq(bSFFrequencyAnalysis, bMeas, choosePeak) ./ higherHarmonic
  tmot = getRepetitionsOfSameState(bMeas, motFreq, frames)

  # sort measured data in virtual frames
  uReco = getMeasurementsMotionCompFD(bMeas, motFreq, tmot, freq, frames, alpha,
                                      samplingPrecision, windowType)

  # subtract background measurement
  if bEmpty != nothing
    uEmpty = getMeasurementsFD(bEmpty, frequencies=freq, numAverages=1, spectralLeakageCorrection=true)
    if bgFrames == nothing
      numFrames = acqNumPeriodsPerPatch(bEmpty)
      bgFrames = [1+(i-1)*numFrames:i*numFrames for i=1:acqNumPatches(bEmpty)]
    end
    for i=1:acqNumPatches(bMeas)
      uReco[:,i,:] = uReco[:,i,:] .- mean(uEmpty[:,bgFrames[i],:],dims=(2,3))
    end
  end

  P = numDFPeriodsInMotionCycle(motFreq, frames, dfCycle(bMeas))

  mapping = collect(1:acqNumPatches(bMeas))
  resortedInd = zeros(Int64, acqNumPatches(bMeas), P)
  for i=1:acqNumPatches(bMeas)
    resortedInd[i,:] = unflattenOffsetFieldShift(FFP)[i][1:P]
  end
  FFOp = MultiPatchOperator(bSF, bMeas, freq, false,
			    indFFPos=resortedInd[:,1],
			    FFPos=FFP[:,resortedInd[:,1]], mapping=mapping,
			    FFPosSF=FFP[:,resortedInd[:,1]])

  image = initImage(bSF[1],bMeas,size(uReco,3),acqNumAverages(bMeas),FFOp.grid,false)
  c =  reconstruction(FFOp, uReco; kargs...)
  writePartToImage!(image, c, 1, 1:size(uReco,3), acqNumAverages(bMeas))

  return image
end

ext: jl


text:
LinearDiscretizer([0.0,0.5,1.0], UInt8)
@test_throws ErrorException LinearDiscretizer([0.0])
@test_throws ErrorException LinearDiscretizer([0.5,0.0])
@test_throws ErrorException LinearDiscretizer([0.0,1.0,1.0])

ld = LinearDiscretizer([0.0,0.5,1.0])
@test encode(ld, -1.0) == 1
@test encode(ld,  0.0) == 1
@test encode(ld,  0.2) == 1
@test encode(ld,  0.5) == 2
@test encode(ld,  0.7) == 2
@test encode(ld,  1.0) == 2
@test encode(ld,  1.2) == 2
@test encode(ld,  Inf) == 2
@test_throws ErrorException encode(ld,  NaN)
@test encode(ld, Float32(0.2)) == 1
@test encode(ld, [-1.0, 0.0, 0.2, 0.5, 0.7, 1.0, 1.2, Inf]) ==
        [1,1,1,2,2,2,2,2]
@test encode(ld, convert(Vector{Float32}, [0.2,0.6])) ==
        [1,2]
@test encode(ld, [0.0 1.0; 0.5 0.2]) == [1 2; 2 1]

@test 0.0 ≤ decode(ld, 1) ≤ 0.5
@test 0.5 ≤ decode(ld, 2) ≤ 1.0
@test 0.0 ≤ decode(ld, 1, SAMPLE_UNIFORM) ≤ 0.5
@test 0.5 ≤ decode(ld, 2, SAMPLE_UNIFORM) ≤ 1.0
@test isapprox(decode(ld, 1, SAMPLE_BIN_CENTER), 0.25)
@test isapprox(decode(ld, 2, SAMPLE_BIN_CENTER), 0.75)

@test_throws BoundsError decode(ld,  0)
@test_throws BoundsError decode(ld,  3)
@test 0.5 ≤ decode(ld, 2 % UInt8) ≤ 1.0
mat = decode(ld, [2 1; 1 2])
@test 0.5 ≤ mat[1,1] ≤ 1.0
@test 0.5 ≤ mat[2,2] ≤ 1.0
@test 0.0 ≤ mat[1,2] ≤ 0.5
@test 0.0 ≤ mat[2,1] ≤ 0.5

@test array_matches([extrema(ld)...], [0.0, 1.0], 1e-8)
@test array_matches([extrema(ld, 1)...], [0.0, 0.5], 1e-8)
@test array_matches([extrema(ld, 2)...], [0.5, 1.0], 1e-8)

@test nlabels(ld) == 2
@test array_matches(bincenters(ld), [0.25,0.75], 1e-8)
@test isapprox([bincenter(ld, i) for i in 1:nlabels(ld)], bincenters(ld))

@test encoded_type(ld) == Int
@test decoded_type(ld) == Float64
@test encoded_type(ld) != Float64
@test decoded_type(ld) != Int

@test supports_encoding(ld, 0.0)
@test supports_encoding(ld, 0.5)
@test supports_encoding(ld, 1.0)
@test supports_encoding(ld, 0.25)
@test supports_encoding(ld, 0.75)
@test supports_encoding(ld, -0.1)
@test supports_encoding(ld, 1.1)

@test supports_decoding(ld, 1)
@test supports_decoding(ld, 2)
@test !supports_decoding(ld, 0)
@test !supports_decoding(ld, 3)

ld = LinearDiscretizer([0.0,0.5,1.0], Int, force_outliers_to_closest=false)
@test supports_encoding(ld, 0.0)
@test supports_encoding(ld, 0.5)
@test supports_encoding(ld, 1.0)
@test supports_encoding(ld, 0.25)
@test supports_encoding(ld, 0.75)
@test !supports_encoding(ld, -0.1)
@test !supports_encoding(ld, 1.1)
@test !supports_encoding(ld, Inf)

############################################################

ld = LinearDiscretizer([0,10,20])
@test encode(ld,  -1) == 1
@test encode(ld,   0) == 1
@test encode(ld,   4) == 1
@test encode(ld,  10) == 2
@test encode(ld,  14) == 2
@test encode(ld,  20) == 2
@test encode(ld,  25) == 2
@test_throws InexactError encode(ld,  NaN)
@test encode(ld, 4 % UInt8 ) == 1
@test encode(ld, [-1,0,4,10,14,20,25]) == [1,1,1,2,2,2,2]
@test encode(ld, convert(Vector{UInt8}, [4,14])) == [1,2]
@test encode(ld, [-1 4; 14 25]) == [1 1; 2 2]

@test  0 ≤ decode(ld, 1) < 10
@test 10 ≤ decode(ld, 2) ≤ 20
@test  0 ≤ decode(ld, 1, SAMPLE_UNIFORM) < 10
@test 10 ≤ decode(ld, 2, SAMPLE_UNIFORM) ≤ 20
@test decode(ld, 1, SAMPLE_BIN_CENTER) == 5
@test decode(ld, 2, SAMPLE_BIN_CENTER) == 15

@test_throws BoundsError decode(ld,  0)
@test_throws BoundsError decode(ld,  3)
@test 10 ≤ decode(ld, 2 % UInt8) ≤ 20
mat = decode(ld, [2 1; 1 2])
@test 10 ≤ mat[1,1] ≤ 20
@test 10 ≤ mat[2,2] ≤ 20
@test  0 ≤ mat[1,2] < 10
@test  0 ≤ mat[2,1] < 10

@test extrema(ld) == (0,20)
@test extrema(ld, 1) == (0,9)
@test extrema(ld, 2) == (10,20)

@test nlabels(ld) == 2
@test array_matches(bincenters(ld), [4.5,15], 1e-8)
@test isapprox([bincenter(ld, i) for i in 1:nlabels(ld)], bincenters(ld))

@test encoded_type(ld) == Int
@test decoded_type(ld) == Int
@test encoded_type(ld) != Float64
@test decoded_type(ld) != Float64

@test supports_encoding(ld, 0)
@test supports_encoding(ld, 10)
@test supports_encoding(ld, 20)
@test supports_encoding(ld, 4)
@test supports_encoding(ld, 16)
@test supports_encoding(ld, -1)
@test supports_encoding(ld, 21)

@test supports_decoding(ld, 1)
@test supports_decoding(ld, 2)
@test !supports_decoding(ld, 0)
@test !supports_decoding(ld, 3)

ld = LinearDiscretizer([0,10,20], Int, force_outliers_to_closest=false)
@test supports_encoding(ld, 0)
@test supports_encoding(ld, 10)
@test supports_encoding(ld, 20)
@test supports_encoding(ld, 4)
@test supports_encoding(ld, 16)
@test !supports_encoding(ld, -1)
@test !supports_encoding(ld, 21)


text:
@doc raw"""
    NelderMead(M, F [, p])
perform a nelder mead minimization problem for the cost funciton `F` on the
manifold `M`. If the initial population `p` is not given, a random set of
points is chosen.

This algorithm is adapted from the Euclidean Nelder-Mead method, see
[https://en.wikipedia.org/wiki/Nelder–Mead_method](https://en.wikipedia.org/wiki/Nelder–Mead_method)
and
[http://www.optimization-online.org/DB_FILE/2007/08/1742.pdf](http://www.optimization-online.org/DB_FILE/2007/08/1742.pdf).

# Input

* `M` – a manifold $\mathcal M$
* `F` – a cost function $F\colon\mathcal M\to\mathbb R$ to minimize
* `population` – (n+1 `random_point(M)`) an initial population of $n+1$ points, where $n$
  is the dimension of the manifold `M`.

# Optional

* `stopping_criterion` – ([`StopAfterIteration`](@ref)`(2000)`) a [`StoppingCriterion`](@ref)
* `α` – (`1.`) reflection parameter ($\alpha > 0$)
* `γ` – (`2.`) expansion parameter ($\gamma$)
* `ρ` – (`1/2`) contraction parameter, $0 < \rho \leq \frac{1}{2}$,
* `σ` – (`1/2`) shrink coefficient, $0 < \sigma \leq 1$
* `retraction_method` – (`ExponentialRetraction`) the rectraction to use
* `inverse_retraction_method` - (`LogarithmicInverseRetraction`) an inverse retraction to use.

and the ones that are passed to [`decorate_options`](@ref) for decorators.

# Output
* either `x` the last iterate or the complete options depending on the optional
  keyword `return_options`, which is false by default (hence then only `x` is
  returned).
"""
function NelderMead(
    M::Manifold,
    F::TF,
    population=[random_point(M) for i in 1:(manifold_dimension(M) + 1)];
    kwargs...,
) where {TF}
    res_population = allocate.(population)
    copyto!.(res_population, population)
    return NelderMead!(M, F, res_population; kwargs...)
end
@doc raw"""
    NelderMead(M, F [, p])
perform a nelder mead minimization problem for the cost funciton `F` on the
manifold `M`. If the initial population `p` is not given, a random set of
points is chosen. If it is given, the computation is done in place of `p`.

For more options see [`NelderMead`](@ref).
"""
function NelderMead!(
    M::Manifold,
    F::TF,
    population=[random_point(M) for i in 1:(manifold_dimension(M) + 1)];
    stopping_criterion::StoppingCriterion=StopAfterIteration(200000),
    α=1.0,
    γ=2.0,
    ρ=1 / 2,
    σ=1 / 2,
    retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
    inverse_retraction_method::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(),
    return_options=false,
    kwargs..., #collect rest
) where {TF}
    p = CostProblem(M, F)
    o = NelderMeadOptions(
        population,
        stopping_criterion;
        α=α,
        γ=γ,
        ρ=ρ,
        σ=σ,
        retraction_method=retraction_method,
        inverse_retraction_method=inverse_retraction_method,
    )
    o = decorate_options(o; kwargs...)
    resultO = solve(p, o)
    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end
#
# Solver functions
#
function initialize_solver!(p::CostProblem, o::NelderMeadOptions)
    # init cost and x
    o.costs = get_cost.(Ref(p), o.population)
    return o.x = o.population[argmin(o.costs)] # select min
end
function step_solver!(p::CostProblem, o::NelderMeadOptions, iter)
    m = mean(p.M, o.population)
    ind = sortperm(o.costs) # reordering for cost and p, i.e. minimizer is at ind[1]
    ξ = inverse_retract(p.M, m, o.population[last(ind)], o.inverse_retraction_method)
    # reflect last
    xr = retract(p.M, m, -o.α * ξ, o.retraction_method)
    Costr = get_cost(p, xr)
    # is it better than the worst but not better than the best?
    if Costr >= o.costs[first(ind)] && Costr < o.costs[last(ind)]
        # store as last
        o.population[last(ind)] = xr
        o.costs[last(ind)] = Costr
    end
    # --- Expansion ---
    if Costr < o.costs[first(ind)] # reflected is better than fist -> expand
        xe = retract(p.M, m, -o.γ * o.α * ξ, o.retraction_method)
        Coste = get_cost(p, xe)
        # successful? use the expanded, otherwise still use xr
        o.population[last(ind)] .= Coste < Costr ? xe : xr
        o.costs[last(ind)] = min(Coste, Costr)
    end
    # --- Contraction ---
    if Costr > o.costs[ind[end - 1]] # even worse than second worst
        s = (Costr < o.costs[last(ind)] ? -o.ρ : o.ρ)
        xc = retract(p.M, m, s * ξ, o.retraction_method)
        Costc = get_cost(p, xc)
        if Costc < o.costs[last(ind)] # better than last ? -> store
            o.population[last(ind)] = xc
            o.costs[last(ind)] = Costc
        end
    end
    # --- Shrink ---
    for i in 2:length(ind)
        o.population[ind[i]] = retract(
            p.M,
            o.population[ind[1]],
            inverse_retract(
                p.M, o.population[ind[1]], o.population[ind[i]], o.inverse_retraction_method
            ),
            o.σ,
            o.retraction_method,
        )
        # update cost
        o.costs[ind[i]] = get_cost(p, o.population[ind[i]])
    end
    # store best
    return o.x = o.population[argmin(o.costs)]
end


text:
# example for dead-reckoning branch on factor graph

# using Revise

using DistributedFactorGraphs
using IncrementalInference
using RoME

using Test

@testset "testing against solvable..." begin

# start with an empty factor graph object
fg = initfg()

# Add the first pose :x0
addVariable!(fg, :x0, Pose2)

# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
addFactor!(fg, [:x0], PriorPose2( MvNormal([10; 10; pi/6.0], Matrix(Diagonal([0.1;0.1;0.05].^2))) ) )

# Drive around in a hexagon
for i in 0:5
    psym = Symbol("x$i")
    nsym = Symbol("x$(i+1)")
    addVariable!(fg, nsym, Pose2)
    pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    addFactor!(fg, [psym;nsym], pp )
end

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, Point2, tags=[:LANDMARK])
p2br = Pose2Point2BearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1], p2br )

## async solving with dead-reckon branch

addVariable!(fg, :deadreckon_x0, Pose2, solvable=0)

drec = MutablePose2Pose2Gaussian(MvNormal(zeros(3), Matrix{Float64}(LinearAlgebra.I, 3,3)))

addFactor!(fg, [:x0; :deadreckon_x0], drec, solvable=0)

#
@test length(map( x->x.label, getVariables(fg, solvable=1))) == 8
@test length(map( x->x.label, getVariables(fg, solvable=0))) == 9
#
# # make sure
@test length(getEliminationOrder(fg, solvable=1)) == 8
# check default
@test length(getEliminationOrder(fg)) == 8

# default check
vo = getEliminationOrder(fg)
@test length(vo) == 8

tree = buildTreeFromOrdering!(fg,vo)

# syms = map( x->x.label, getVariables(fg, solvable=1))

# union(map(x->getCliqAllVarIds(x), map(n->getClique(tree, n), syms))...)

getSolverParams(fg).dbg = true
# getSolverParams(fg).drawtree = true
# getSolverParams(fg).showtree = true

# breaks because internal saveDFG somewhere is failing on bad converter for MutablePose2Pose2Gaussian
tree2, smt, hist = solveTree!(fg, recordcliqs=ls(fg));

@test !isInitialized(fg, :deadreckon_x0)

val = accumulateFactorMeans(fg, [:x0deadreckon_x0f1])

@test norm(val - calcVariablePPE(fg, :x0).suggested) < 1e-8

end

#


text:
# by JMW
macro delegate(source, targets)
    typename = esc(source.args[1])
    fieldname = source.args[2].args[1]
    funcnames = targets.args
    n = length(funcnames)
    fdefs = Vector{Any}(n)
    for i in 1:n
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typename), args...) =
                       ($funcname)(a.$fieldname, args...)
                   end
    end
    return Expr(:block, fdefs...)
end

macro delegate_return_parent(source, targets)
    typename = esc(source.args[1])
    fieldname = source.args[2].args[1]
    funcnames = targets.args
    n = length(funcnames)
    fdefs = Vector{Any}(n)
    for i in 1:n
        funcname = esc(funcnames[i])
        fdefs[i] = quote
                     ($funcname)(a::($typename), args...) =
                       (($funcname)(a.$fieldname, args...); a)
                   end
    end
    return Expr(:block, fdefs...)
end


text:
type ArbResamplerState
    rate::Float64
    N𝜙::Int
    Δ::Float64
    counter::Float64
    𝜙Accumulator::Float64
    𝜙Idx::Int
    δ::Float64
    𝜙IdxVirtual::Float64
    yLower::Number
    function ArbResamplerState( rate::Real, N𝜙::Integer = 32 )
        rate         = rate
        N𝜙           = N𝜙
        Δ            = 1.0/rate
        counter      = 0
        𝜙Accumulator = 0
        𝜙IdxVirtual  = 1
        𝜙Idx         = 1
        δ            = 0
        yLower       = 0
        self = new( rate, N𝜙, Δ, counter, 𝜙Accumulator, 𝜙Idx, δ, 𝜙IdxVirtual, yLower )
        # update( self )
        self
    end
end

function update( self::ArbResamplerState )
    self.counter += self.Δ
    self.𝜙Accumulator += self.Δ
    if self.𝜙Accumulator >= 1
        self.𝜙Accumulator = mod( self.𝜙Accumulator, 1 )
        self.counter += 1
    end

    self.𝜙IdxVirtual = self.𝜙Accumulator * ( self.N𝜙 - 1 ) + 1
    self.𝜙Idx        = floor( Int, self.𝜙IdxVirtual )
    self.δ           = self.𝜙IdxVirtual - self.𝜙Idx
    self
end

resamp = 2
N𝜙     = 10
yCount = 0
xCount = 0
self   = ArbResamplerState( resamp, N𝜙 )

while xCount < 11
    xCount += 1
    @printf( "%d:\tcounter = %f, 𝜙Accumulator = %f, 𝜙IdxVirtual = %f, 𝜙Idx = %d, δ = %f\n", xCount, self.counter, self.𝜙Accumulator, self.𝜙IdxVirtual, self.𝜙Idx, self.δ)
    update( self )
end
#
# resamp       = 0.9
# N𝜙           = 32
# yCount       = 0
# xCount       = 0
# 𝜙Idx         = 0
# δ            = 0.0
# 𝜙IdxVirtual  = 0.0
# 𝜙Accumulator = 0.0
# Δ = Int(resamp)
#
# while xCount < 10
#     xCount += 1
#     while 𝜙Idx <= N𝜙
#         yCount       += 1
#
#         println( "$yCount: 𝜙Idx = $𝜙Idx, δ = $δ, 𝜙IdxVirtual = $𝜙IdxVirtual, 𝜙Accumulator = $𝜙Accumulator")
#
#         𝜙Accumulator += Δ
#         𝜙IdxVirtual   = 𝜙Accumulator * N𝜙
#         𝜙Idx          = floor( Int, 𝜙IdxVirtual ) + 1
#         δ             = mod( 𝜙Accumulator, 1 )
#
#     end
#
#     𝜙Accumulator -= 1
#     𝜙IdxVirtual  -= N𝜙
#     𝜙Idx          = floor( Int, 𝜙IdxVirtual )
#
# end

max_stars_repo_path: test/Δ Test.jl

max_issues_repo_path: test/Δ Test.jl

max_forks_repo_path: test/Δ Test.jl


text:
#=
This file has some quick speed checks of AxisKeys functionality,
and of other packages of similar or overlapping concerns.
Plus generally a place to collect their various syntax, for comparison.
=#

using AxisKeys, BenchmarkTools # Julia 1.4 + macbook escape

#==============================#
#===== getkey vs getindex =====#
#==============================#

mat = wrapdims(rand(3,4), 11:13, 21:24)
bothmat = wrapdims(mat.data, x=11:13, y=21:24)
bothmat2 = wrapdims(mat.data, x=collect(11:13), y=collect(21:24))

@btime $mat[3, 4]    # 1.699 ns
@btime $mat(13, 24)  # 5.312 ns

@btime $bothmat[3,4]        # 1.700 ns
@btime $bothmat[x=3, y=4]   # 1.701 ns
@btime $bothmat(13, 24)     # 5.874 ns -- fast range lookup
@btime $bothmat(x=13, y=24) # 14.063 ns
@btime $bothmat2(13, 24)    # 16.719 ns -- collected vector

ind_collect(A) = [@inbounds(A[ijk...]) for ijk in Iterators.ProductIterator(axes(A))]
key_collect(A) = [@inbounds(A(vals...)) for vals in Iterators.ProductIterator(axiskeys(A))]

bigmat = wrapdims(rand(100,100), 1:100, 1:100);
bigmat2 = wrapdims(rand(100,100), collect(1:100), collect(1:100));

@btime ind_collect($(bigmat.data)); #  9.117 μs (4 allocations: 78.25 KiB)
@btime ind_collect($bigmat);        # 11.530 μs (4 allocations: 78.25 KiB)
@btime key_collect($bigmat);        # 64.064 μs (4 allocations: 78.27 KiB) -- fast range lookup, was 20μs!
@btime key_collect($bigmat2);      # 718.804 μs (5 allocations: 78.27 KiB) -- findfirst(..., vector) lookup

twomat = wrapdims(mat.data, x=[:a, :b, :c], y=21:24)
@btime $twomat(x=:a, y=24)  # 36.734 ns (2 allocations: 64 bytes)

@btime $twomat(24.0)        # 26.686 ns (4 allocations: 112 bytes)
@btime $twomat(y=24.0)      # 33.860 ns (4 allocations: 112 bytes)
@btime view($twomat, :,3)   # 24.951 ns (4 allocations: 112 bytes)

#==========================#
#===== other packages =====#
#==========================#

using OffsetArrays                      #===== OffsetArrays =====#

of1 = OffsetArray(rand(3,4), 11:13, 21:24)

@btime $of1[13,24] # 1.700 ns

bigoff = OffsetArray(bigmat.data, 1:100, 1:100);
@btime ind_collect($bigoff);        # 15.372 μs (5 allocations: 78.30 KiB)

using AxisArrays                        #===== AxisArrays =====#

ax1 = AxisArray(rand(3,4), 11:13, 21:24)
ax2 = AxisArray(rand(3,4), :x, :y)

@btime $ax1[3,4] # 1.421 ns
@btime $ax2[3,4] # 1.421 ns
@btime $ax2[Axis{:x}(3), Axis{:y}(4)] # 1.696 ns

@btime $ax1[atvalue(13), atvalue(24)] # 21.348 ns

bigax = AxisArray(bigmat.data, 1:100, 1:100);
bigax2 = AxisArray(bigmat.data, collect(1:100), collect(1:100));
ax_collect(A) = [@inbounds(A[atvalue(vals[1]), atvalue(vals[2])]) for vals in Iterators.ProductIterator(AxisArrays.axes(A))]
@btime ind_collect($bigax);        # 9.946 μs (4 allocations: 78.25 KiB) -- indexing
@btime ax_collect($bigax);         # 212.160 μs (6 allocations: 78.34 KiB) -- range lookup?
@btime ax_collect($bigax2);        # 511.157 μs (5 allocations: 78.27 KiB) -- compare to findfirst(...) time

using NamedArrays                       #===== NamedArrays =====#

na1 = NamedArray(rand(3,4))

na2 = copy(na1);
setnames!(na2, ["a", "b", "c"], 1); na2
setnames!(na2, string.(1:4) .* "th", 2); na2 # must be strings, and all different

@btime $na1[3,4]     # 4.467 ns
@btime $na2["c",4]   # 22.592 ns
@btime $na2["c","4th"] # 40.851 ns
# This lookup is via OrderedDict{String,Int64}
@btime $na2.dicts[2]["4th"] # 22.225 ns

na1.dimnames == (:A, :B)
# @btime na1[:A => 3, :B => 4] # error
@btime $na2[:A => "c", :B => "4th"] # 1.077 μs (29 allocations: 2.02 KiB)

using NamedDims                         #===== NamedDims =====#

nd1 = NamedDimsArray(rand(3,4), (:x, :y))

@btime $nd1[3,4]       #  1.421 ns
@btime $nd1[x=3, y=4]  #  1.696 ns # fixed by https://github.com/invenia/NamedDims.jl/pull/84

using DimensionalData                   #===== DimensionalData =====#
using DimensionalData: X, Y, @dim, Forward

dd1 = DimensionalArray(rand(3,4), (Dim{:x}(11:13), Dim{:y}(21:24)))
@dim dd_x; @dim dd_y
dd2 = DimensionalArray(rand(3,4), (dd_x(11:13), dd_y(21:24)))
dd3 = DimensionalArray(rand(3,4), (X(11:13), Y(21:24)))

@btime $dd1[3,4] # 1.421 ns
@btime $dd1[At(13), At(24)] # 23.149 ns
# StackOverflowError? no, was missing a (,) see https://github.com/rafaqz/DimensionalData.jl/issues/4
@btime $dd1[Dim{:x}(1), Dim{:y}(2)] # 1.421 ns

@btime $dd2[At(13), At(24)] # 23.394 ns
@btime $dd2[dd_x <| At(13), dd_y <| At(24)] # 23.172 ns

@btime $dd3[X <| At(13), Y <| At(24)] # 23.256 ns
@btime (m -> m[X <| At(13), Y <| At(24)])($dd3) # 23.396 ns

using AbstractIndices                   #===== AbstractIndices =====#
# https://github.com/Tokazama/AbstractIndices.jl

ai1 = IndicesArray(rand(3,4), 11:13, 21:24)
ai2 = IndicesArray(rand(3,4); x=11:13, y=21:24)
ai3 = IndicesArray(rand(3,4), 11:13.0, 21:24.0)

# ai1[3,4] # all errors?
# ai1[13,24]
# ai2[x=13,y=24]
@btime $ai3[3,4] # 2.530 ns
@btime $ai3[13.0, 24.0] # 31.437 ns

bigai = IndicesArray(bigmat.data, 1:100.0, 1:100.0);
ai_collect(A) = [@inbounds(A[vals...]) for vals in Iterators.ProductIterator(map(collect, axes(A)))]
# I think axes(A) in ind_collect() is resulting in lookup not indexing:
@btime ind_collect($bigai);        # 645.921 μs (30205 allocations: 1.00 MiB) -- compare to findfirst(...) time?
@btime ai_collect($bigai);         #  24.656 μs (8 allocations: 80.05 KiB) -- compare to indexing

using AxisIndices                       #===== AxisIndices =====#
# https://github.com/Tokazama/AxisIndices.jl

ai4 = AxisIndicesArray(rand(3,4), (11:13, 21:24))
ai5 = AxisIndicesArray(rand(3,4), (collect(11:13), collect(21:24)))

ax41 = axes(ai4, 1) # an Axis type
collect(ax41) == 1:3 == ax41._values
ax41._keys
zeros(ax41) isa OffsetArray # like zeros(1:3), unlike zeros(Base.OneTo(3))

@btime $ai4[3,4] # 1.696 ns
@btime $ai4[==(13),==(24)] #  7.257 ns -- fast lookup
@btime $ai5[==(13),==(24)] # 20.676 ns -- search

ai4[:, x -> x<23] # arbitrary functions allowed

ai4' .+ rand(4) # ok
ai4 .+ ai5 # error
ai4 .+ AxisIndicesArray(rand(3,4), (13:15, 21:24)) # has 13:15

using IndexedDims                       #===== IndexedDims =====#
# https://github.com/invenia/IndexedDims.jl

id1 = IndexedDimsArray(rand(3,4), 11:13, 21:24)
id2 = IndexedDimsArray(rand(3,4), 11.0:13.0, 21.0:24.0)

@btime $id1[3,4] # 1.695 ns
@btime $id2[11.0, 24.0] # 544.219 ns

bigid = IndexedDimsArray(bigmat.data, 1:100.0, 1:100.0);
idx_collect(A) = [@inbounds(A[vals...]) for vals in Iterators.ProductIterator(A.indexes)]

@btime ind_collect($bigid);         # 12.167 μs (4 allocations: 78.25 KiB)
@btime idx_collect($bigid);     # 26.015 ms (150004 allocations: 14.42 MiB) -- compare to findfirst(...) time?

#=============================#
#===== fast range lookup =====#
#=============================#

const A = AxisKeys
const B = Base
using Base: OneTo

@btime B.findfirst(isequal(300), OneTo(1000)) #  173.225 ns
@btime A.findfirst(isequal(300), OneTo(1000)) #  0.029 ns

@btime B.findfirst(isequal(300), 0:1000)      #  90.410 ns
@btime A.findfirst(isequal(300), 0:1000)      #   0.029 ns

@btime B.findall( <(10), OneTo(1000));        # 742.162 ns
@btime A.findall( <(10), OneTo(1000));        #   0.029 ns
@btime collect(A.findall( <(10), OneTo(1000)));# 28.855 ns

@btime B.findall( <(10), 1:1000);             # 890.300 ns
@btime A.findall( <(10), 1:1000);             #   0.029 ns
@btime collect(A.findall( <(10), 1:1000));    #  28.888 ns

B.findall( <(10), 1:1000) isa Vector
A.findall( <(10), 1:1000) isa OneTo

#==============================#
#===== accelerated lookup =====#
#==============================#

using AcceleratedArrays

s1 = sort(rand(1:25, 100));
s2 = accelerate(s1, SortIndex);
@btime findall(isequal(10), $s1)  # 176.629 ns (6 allocations: 288 bytes)
@btime findall(isequal(10), $s2)  #  52.173 ns (3 allocations: 176 bytes)

w1 = wrapdims(1:100.0, x=s1);
w2 = wrapdims(1:100.0, x=s2);
@btime $w1(==(10))       # 538.672 ns (21 allocations: 976 bytes)
@btime $w2(==(10))       # 397.184 ns (18 allocations: 864 bytes)

@btime $w1(10) # 98.399 ns (4 allocations: 256 bytes)
@btime $w2(10) # no help for findfirst

r1 = collect(1:1000);
r2 = accelerate(r1, SortIndex);
@btime findall(<(99), $r1);   # 1.426 μs (11 allocations: 2.36 KiB)
@btime findall(<(99), $r2);   # 1.688 μs (11 allocations: 2.36 KiB)
r0 = 1:1000
@btime A.findall(<(99), $r0); # 0.029 ns (0 allocations: 0 bytes)

bigmat3 = wrapdims(rand(100,100), accelerate(collect(1:100),SortIndex), accelerate(collect(1:100),SortIndex));
@btime key_collect($bigmat3);      # 765.999 μs (5 allocations: 78.27 KiB) -- slower, findfirst.

using DataStructures

v1 = (1:100)[sortperm(rand(100))];
@btime findfirst(isequal(10), $v1) # 54.311 ns

v2 = let
    d = OrderedDict{Int, Int}()
    for (i,k) in enumerate(v1)
        d[k] = i
    end
    d
end
@btime $v2[10] # 7.818 ns

v3 = accelerate(v1, UniqueHashIndex);
@btime findfirst(isequal(10), v3) # 144.299 ns

using CategoricalArrays

ca1 = wrapdims(rand(4), x = CategoricalArray(["a", "b", "a", "b"]))
ca1("a")
ca1(==("a"))

v4 = CategoricalArray(v1)
v4[1] == v1[1]
@btime findfirst(isequal(10), $v4) # 187.328 ns

#===================#
#===== The End =====#
#===================#


text:
#
# Structure to store auxiliary arrays
#
struct VAux
  g :: Vector{Float64}
  xn :: Vector{Float64}
  gn :: Vector{Float64}
  fprev :: Vector{Float64}
end
VAux(n,m) = VAux( Vector{Float64}(undef,n),
                  Vector{Float64}(undef,n),
                  Vector{Float64}(undef,n),
                  Vector{Float64}(undef,m) )
# The default value for m
VAux(n) = VAux(n,10)


text:
# Gaussian Stein Distribution
#
# Represents a multivariate or univariate Gaussian distribution
using Distributions: MvNormalCanon, Normal
import Distributions

type SteinGaussian{T<:Number} <: SteinDistribution
    mu::Array{Float64}
    precision::AbstractArray{T,2}
    c1::Float64
    c2::Float64
    c3::Float64
end

# Constructor with Wasserstein-bounding non-uniform setting of
# (c1,c2,c3) = (1,1,1)
SteinGaussian(mu::Array{Float64}, precision::AbstractArray{Float64}) =
SteinGaussian(mu, precision, 1.0, 1.0, 1.0)

# Constructor for p independent standard Gaussian variables
SteinGaussian(p::Int64, c1::Float64, c2::Float64, c3::Float64) =
SteinGaussian(zeros(p), eye(p), c1, c2, c3);

# Constructor for p independent standard Gaussian variables with canonical
# constants c1, c2, c3
SteinGaussian(p::Int64) = SteinGaussian(zeros(p), eye(p))

# Constructor for Univariate standard Gaussian with canonical
# constants c1, c2, c3
SteinGaussian() = SteinGaussian(1)

# Draw n independent samples from distribution
function rand(d::SteinGaussian, n::Int64)
    Distributions.rand(MvNormalCanon(d.precision*d.mu, d.precision), n).';
end

function supportlowerbound(d::SteinGaussian, j::Int64)
    -Inf
end

function supportupperbound(d::SteinGaussian, j::Int64)
    Inf
end

function gradlogdensity(d::SteinGaussian, x::Array{Float64,1})
    d.precision * (d.mu - x)
end

# Gradient of log density of distribution evaluated at each row of
# n x p matrix X
function gradlogdensity(d::SteinGaussian, X::Array{Float64,2})
    (d.mu.' .- X) * d.precision;
end

# Cumulative distribution function (only valid when X is univariate)
function cdf(d::SteinGaussian, t)
    Distributions.cdf(Normal(d.mu[1],1/sqrt(d.precision[1])),t);
end

function numdimensions(d::SteinGaussian)
    length(d.mu)
end

# diaghessian(x) = (grad_1 b1(x), grad_2 b2(x), ....) where b = grad log p.
function diaghessian(d::SteinGaussian, x::Array{Float64,1})
    -diag(d.precision)
end


text:
"""
A set of callback functions to be used with an `AbstractODESolver`
"""
module GenericCallbacks
using MPI

using ..ODESolvers

"""
    EveryXWallTimeSeconds(f, time, mpicomm)

This callback will run the function 'f()' every `time` wallclock time seconds.
The `mpicomm` is used to syncronize runtime across MPI ranks.
"""
struct EveryXWallTimeSeconds
    "time of the last callback"
    lastcbtime_ns::Array{UInt64}
    "time between callbacks"
    Δtime::Real
    "MPI communicator"
    mpicomm
    "function to execute for callback"
    func::Function
    function EveryXWallTimeSeconds(func, Δtime, mpicomm)
        lastcbtime_ns = [time_ns()]
        new(lastcbtime_ns, Δtime, mpicomm, func)
    end
end
function (cb::EveryXWallTimeSeconds)(initialize::Bool = false)
    # Is this an initialization call? If so, start the timers
    if initialize
        cb.lastcbtime_ns[1] = time_ns()
        # If this is initialization init the callback too
        try
            cb.func(true)
        catch
        end
        return 0
    end

    # Check whether we should do a callback
    currtime_ns = time_ns()
    runtime = (currtime_ns - cb.lastcbtime_ns[1]) * 1e-9
    runtime = MPI.Allreduce(runtime, MPI.MAX, cb.mpicomm)
    if runtime < cb.Δtime
        return 0
    else
        # Compute the next time to do a callback
        cb.lastcbtime_ns[1] = currtime_ns
        return cb.func()
    end
end

"""
   EveryXSimulationTime(f, time, state)

This callback will run the function 'f()' every `time` wallclock time seconds.
The `state` is used to query for the simulation time.
"""
struct EveryXSimulationTime
    "time of the last callback"
    lastcbtime::Array{Real}
    "time between callbacks"
    Δtime::Real
    "function to execute for callback"
    func::Function
    "timestepper, used to query for time"
    solver
    function EveryXSimulationTime(func, Δtime, solver)
        lastcbtime = [ODESolvers.gettime(solver)]
        new(lastcbtime, Δtime, func, solver)
    end
end
function (cb::EveryXSimulationTime)(initialize::Bool = false)
    # Is this an initialization call? If so, start the timers
    if initialize
        cb.lastcbtime[1] = ODESolvers.gettime(cb.solver)
        # If this is initialization init the callback too
        try
            cb.func(true)
        catch
        end
        return 0
    end

    # Check whether we should do a callback
    currtime = ODESolvers.gettime(cb.solver)
    if (currtime - cb.lastcbtime[1]) < cb.Δtime
        return 0
    else
        # Compute the next time to do a callback
        cb.lastcbtime[1] = currtime
        return cb.func()
    end
end

"""
   EveryXSimulationSteps(f, steps)

This callback will run the function 'f()' every `steps` of the time stepper
"""
struct EveryXSimulationSteps
    "number of steps since last callback"
    steps::Array{Int}
    "number of steps between callbacks"
    Δsteps::Integer
    "function to execute for callback"
    func::Function
    function EveryXSimulationSteps(func, Δsteps)
        steps = [Int(0)]
        new(steps, Δsteps, func)
    end
end
function (cb::EveryXSimulationSteps)(initialize::Bool = false)
    # Is this an initialization call? If so, start the timers
    if initialize
        cb.steps[1] = 0
        # If this is initialization init the callback too
        try
            cb.func(true)
        catch
        end
        return 0
    end

    # Check whether we should do a callback
    cb.steps[1] += 1
    if cb.steps[1] < cb.Δsteps
        return 0
    else
        cb.steps[1] = 0
        return cb.func()
    end
end

end


text:
# This file is part of Kpax3. License is MIT.

include("functions.jl")

include("io.jl")

include("merge.jl")
include("split.jl")
include("gibbs.jl")
include("biased_random_walk.jl")

include("algorithm.jl")

include("diagnostics.jl")


text:
using Test

@testset "VSCodeServer" begin
end


text:
#=
# Datadeps

Rfam.jl can fetch key files from the server, which can be used with Infernal.jl.
=#

import Rfam

# Path to `Rfam.cm` file containing covariance models of all the families.

Rfam.cm()


text:
module WMesa

using PyCall
using AbstractActuators

export WMesaClient
export move, devposition, setreference, rmove, numaxes, axesnames
export waituntildone, stopmotion, absposition, setabsreference, moveto

struct WMesaClient <: AbstractRobot
    devname::String
    ip::String
    port::Int32
    axis::String
    server::PyObject
end

AbstractActuators.numaxes(dev::WMesaClient) = 1
AbstractActuators.numaxes(::Type{WMesaClient}) = 1
AbstractActuators.axesnames(dev::WMesaClient) = [dev.axis]

function WMesaClient(devname, ip="192.168.0.140", port=9596; axis="θ")
    xmlrpc = pyimport("xmlrpc.client")
    server = xmlrpc.ServerProxy("http://$ip:$port")
    WMesaClient(devname, ip, port, axis, server)
end

AbstractActuators.move(dev::WMesaClient, deg; a=false, r=false, sync=true) =
    dev.server["move"](deg, a, r, sync)

AbstractActuators.moveto(dev::WMesaClient, deg) =
    move(dev, deg[1]; a=false, r=false, sync=true)

AbstractActuators.rmove(dev::WMesaClient, deg; sync=true) =
    dev.server["move"](deg, false, true, sync)

AbstractActuators.devposition(dev::WMesaClient) = dev.server["position"]()
AbstractActuators.absposition(dev::WMesaClient) = dev.server["abs_position"]()

AbstractActuators.setreference(dev::WMesaClient, deg=0) = dev.server["set_reference"](deg)

AbstractActuators.setabsreference(dev::WMesaClient) = dev.server["set_abs_reference"]()

AbstractActuators.waituntildone(dev::WMesaClient) = dev.server["waitUntilDone"]()

AbstractActuators.stopmotion(dev::WMesaClient) = dev.server["stop"]()

end


text:
using SimplexTrees
using Base.Test

include("simplex_array_list.jl")
#include("simplex_tree.jl")


text:
struct Broker
    workers::Dict{WorkerSpec, WorkerPool}
end

Broker() = Broker(Dict())

function request!(broker::Broker, spec::WorkerSpec, method::AbstractString, params)
    withworker!(get!(broker.workers, spec, WorkerPool(spec))) do worker
        return request(worker, method, params)
    end
end

function evalon!(
    broker::Broker,
    spec::WorkerSpec,
    code::AbstractString,
    options::AbstractDict,
)
    @assert !haskey(options, "code")  # what about `:code`?
    params = merge(Dict("code" => code), options)
    return request!(broker, spec, "eval", params)
end

function popspec!(options)
    julia = pop!(options, "julia", nothing)
    julia isa AbstractString ||
        return ArgumentError("Keyword argument `julia` of type `String` is required.")

    project = pop!(options, "project", nothing)
    project isa AbstractString || project === nothing ||
        return ArgumentError("Keyword argument `project` must be a `String`.")

    return WorkerSpec(julia, project)
end

struct EvalAPI
    broker::Broker
end

function (api::EvalAPI)(input)
    if !(input isa AbstractDict)
        return ArgumentError("API `eval` requires keyword argument.")
    end
    # Pass-through unhandled arguments to the worker:
    options = copy(input)

    spec = @return_exception popspec!(options)

    code = pop!(options, "code", nothing)
    code isa AbstractString ||
        return ArgumentError("Keyword argument `code` of type `AbstractString` is required.")

    return evalon!(api.broker, spec, code, options)
end

struct RunAPI
    broker::Broker
end

function (api::RunAPI)(input)
    if !(input isa AbstractDict)
        return ArgumentError("API `run` requires keyword argument.")
    end
    # Pass-through unhandled arguments to the worker:
    options = copy(input)

    script = pop!(options, "script", nothing)
    script isa AbstractString ||
        return ArgumentError("Keyword argument `script` of type `String` is required.")

    @unpack julia, project, code = @return_exception parse_script(script)

    return evalon!(api.broker, WorkerSpec(julia, project), code, options)
end

struct CallMainAPI
    broker::Broker
end

function (api::CallMainAPI)(input)
    if !(input isa AbstractDict)
        return ArgumentError("API `callmain` requires keyword argument.")
    end
    # Pass-through unhandled arguments to the worker:
    options = copy(input)

    spec = @return_exception popspec!(options)
    return request!(api.broker, spec, "callmain", options)
end

struct CallAnyAPI
    broker::Broker
end

function (api::CallAnyAPI)(input)
    if !(input isa AbstractDict)
        return ArgumentError("API `callany` requires keyword argument.")
    end
    # Pass-through unhandled arguments to the worker:
    options = copy(input)

    spec = @return_exception popspec!(options)
    return request!(api.broker, spec, "callany", options)
end

struct HelpAPI
    broker::Broker
end

function (api::HelpAPI)(input)
    if !(input isa AbstractDict)
        return ArgumentError("API `help` requires keyword argument.")
    end
    # Pass-through unhandled arguments to the worker:
    options = copy(input)

    spec = @return_exception popspec!(options)
    return request!(api.broker, spec, "help", options)
end

struct AdHocCLIAPI
    broker::Broker
end

function option_as_kwargname(opt::AbstractString)
    opt = replace(opt, r"^--?" => "")
    opt = replace(opt, "-" => "_")
    return opt
end

function auto_parse_option(val::AbstractString)
    return @something(
        tryparse(Bool, val),
        tryparse(Int, val),
        tryparse(Float64, val),
        val,
    )
end

function parse_adhoc_cli_args(cliargs::AbstractVector)
    args = String[]
    kwargs = Dict{String,Any}()
    for a in cliargs
        if startswith(a, "-")
            if occursin("=", a)
                k, v = split(a, "="; limit=2)
                kwargs[option_as_kwargname(k)] = auto_parse_option(v)
            else
                kwargs[option_as_kwargname(a)] = true
            end
        else
            push!(args, a)
        end
    end
    return args, kwargs
end

function (api::AdHocCLIAPI)(input)
    if !(input isa AbstractDict)
        return ArgumentError("API `adhoccli` requires keyword argument.")
    end
    # Pass-through unhandled arguments to the worker:
    options = copy(input)

    spec = @return_exception popspec!(options)
    cliargs = @return_exception getstrings(options, "args")

    if ["--help"] == cliargs || ["-h"] == cliargs
        return request!(api.broker, spec, "help", options)
    end

    options["args"], options["kwargs"] = parse_adhoc_cli_args(cliargs)
    return request!(api.broker, spec, "callany", options)
end

function serve(server)
    broker = Broker()
    dispatcher = Dict{String, Any}(
        "eval" => EvalAPI(broker),
        "run" => RunAPI(broker),
        "callmain" => CallMainAPI(broker),
        "callany" => CallAnyAPI(broker),
        "help" => HelpAPI(broker),
        "adhoccli" => AdHocCLIAPI(broker),
    )
    JSONRPC.serve(dispatcher, server, JSONRPC.NDJSON)
end

function serve()
    path = joinpath(DEPOT_PATH[1], "jlcli", "socket")
    jlclidir = dirname(path)
    if !isdir(jlclidir)
        mkpath(jlclidir)
    end
    serve(path)
end


text:
using Documenter, PETSc2

makedocs(
  format = :html,
  sitename = "Petsc2.jl",
  pages = [
    "index.md",
    "build.md",
    "init.md",
    "vec.md",
    "vec_interface.md",
    "mat.md",
    "mat_interface.md",
    "constants.md",
    "ksp.md",
    "pc.md",
    "options.md",
    "error.md"
  ]
)


text:
"""
    li2(x::Float64)::Float64

Returns the real dilogarithm of a real number `x` of type `Float64`.

Author: Alexander Voigt

License: MIT

# Example
```julia
li2(1.0)
```
"""
function li2(x::Float64)::Float64
    cp = (
        0.9999999999999999502e+0,
       -2.6883926818565423430e+0,
        2.6477222699473109692e+0,
       -1.1538559607887416355e+0,
        2.0886077795020607837e-1,
       -1.0859777134152463084e-2
    )
    cq = (
        1.0000000000000000000e+0,
       -2.9383926818565635485e+0,
        3.2712093293018635389e+0,
       -1.7076702173954289421e+0,
        4.1596017228400603836e-1,
       -3.9801343754084482956e-2,
        8.2743668974466659035e-4
    )

    # transform to [0, 1/2]
    (y, r, s) = if x < -1.0
        l = log(1.0 - x)
        (1.0/(1.0 - x), -pi*pi/6.0 + l*(0.5*l - log(-x)), 1.0)
    elseif x == -1.0
        return -pi*pi/12.
    elseif x < 0.0
        (x/(x - 1.0), -0.5*log1p(-x)^2, -1.0)
    elseif x == 0.0
        return 0.0
    elseif x < 0.5
        (x, 0.0, 1.0)
    elseif x < 1.0
        (1.0 - x, pi*pi/6.0 - log(x)*log(1.0 - x), -1.0)
    elseif x == 1.0
        return pi*pi/6.0
    elseif x < 2.0
        l = log(x)
        (1.0 - 1.0/x, pi*pi/6.0 - l*(log(1.0 - 1.0/x) + 0.5*l), 1.0)
    else
        (1.0/x, pi*pi/3.0 - 0.5*log(x)^2, -1.0)
    end

    y2 = y*y
    y4 = y2*y2

    p = cp[1] + y * cp[2] + y2 * (cp[3] + y * cp[4]) +
        y4 * (cp[5] + y * cp[6])
    q = cq[1] + y * cq[2] + y2 * (cq[3] + y * cq[4]) +
        y4 * (cq[5] + y * cq[6] + y2 * cq[7])

    r + s*y*p/q
end

"""
    li2(z::ComplexF64)::ComplexF64

Returns the complex dilogarithm of a complex number `z` of type `ComplexF64`.

Author: Alexander Voigt

License: MIT

# Example
```julia
li2(1.0 + 1.0im)
```
"""
function li2(z::ComplexF64)::ComplexF64
    clog(z) = 0.5*log(abs2(z)) + angle(z)*1.0im

    rz = real(z)
    iz = imag(z)

    if iz == 0.0
        if rz <= 1.0
            return li2(rz)
        else # rz > 1.
            return li2(rz) - pi*log(rz)*1.0im
        end
    end

    nz = abs2(z)

    if nz < eps(Float64)
        return z*(1.0 + 0.25*z)
    end

    (u::ComplexF64, rest::ComplexF64, sgn::Float64) = if rz <= 0.5
        if nz > 1.0
            (-clog(1.0 - 1.0 / z), -0.5 * clog(-z)^2 - pi^2 / 6.0, -1.0)
        else # nz <= 1.
            (-clog(1.0 - z), 0.0 + 0.0im, 1.0)
        end
    else # rz > 0.5
        if nz <= 2.0*rz
            l = -clog(z)
            (l, l * clog(1.0 - z) + pi^2 / 6.0, -1.0)
        else # nz > 2.0*rz
            (-clog(1.0 - 1.0 / z), -0.5 * clog(-z)^2 - pi^2 / 6.0, -1.0)
        end
    end

    bf = (
        - 1.0/4.0,
          1.0/36.0,
        - 1.0/3600.0,
          1.0/211680.0,
        - 1.0/10886400.0,
          1.0/526901760.0,
        - 4.0647616451442255e-11,
          8.9216910204564526e-13,
        - 1.9939295860721076e-14,
          4.5189800296199182e-16
    )

    u2::ComplexF64 = u*u
    u4::ComplexF64 = u2*u2

    rest + sgn * (
        u +
        u2 * (bf[1] +
        u  * (bf[2] +
        u2 * (
            bf[3] +
            u2*bf[4] +
            u4*(bf[5] + u2*bf[6]) +
            u4*u4*(bf[7] + u2*bf[8] + u4*(bf[9] + u2*bf[10]))
        )))
    )
end


text:
module FESpaceFactoriesTests

using Test

using Gridap.Geometry
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.FESpaces
using Gridap.CellData
using Gridap.FESpaces: ExtendedFESpace

domain = (0,1,0,1)
partition = (10,10)
model = CartesianDiscreteModel(domain,partition)
order = 3

# From a single reffe

V = FESpace(model,QUAD4)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,QUAD4,conformity=:H1)
@test isa(V,UnconstrainedFESpace)
@test num_free_dofs(V) == 121

V = FESpace(model,QUAD4,conformity=:L2)
@test isa(V,UnconstrainedFESpace)
@test num_free_dofs(V) == 400

V = FESpace(model,QUAD4,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

# From parameter list describing the reffe

reffe = ReferenceFE(lagrangian,Float64,order)

V = FESpace(model,reffe)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,reffe,conformity=:H1)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,reffe,conformity=:L2)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,reffe,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

# From a cell-wise vector of reffes

cell_reffe = ReferenceFE(model,lagrangian,Float64,order)

V = FESpace(model,cell_reffe)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,cell_reffe,conformity=:H1)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,cell_reffe,conformity=:L2)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,cell_reffe,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

# From a CellFE

cell_fe = FiniteElements(PhysicalDomain(),model,lagrangian,Float64,order)

V = FESpace(model,cell_fe)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,cell_fe,conformity=:L2)
@test isa(V,UnconstrainedFESpace)

V = FESpace(model,cell_fe,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

# Now from a restricted model

const R = 0.7

function is_in(coords)
  n = length(coords)
  x = (1/n)*sum(coords)
  d = x[1]^2 + x[2]^2 - R^2
  d < 0
end

Ω = Triangulation(model)
cell_coords = get_cell_coordinates(Ω)
cell_mask = lazy_map(is_in,cell_coords)
model_in = DiscreteModel(model,cell_mask)

# From a single reffe

V = FESpace(model_in,QUAD4)
@test isa(V,ExtendedFESpace)

V = FESpace(model_in,QUAD4,conformity=:L2)
@test isa(V,ExtendedFESpace)

V = FESpace(model_in,QUAD4,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

# From parameter list describing the reffe

reffe = ReferenceFE(lagrangian,Float64,order)

V = FESpace(model_in,reffe)
@test isa(V,ExtendedFESpace)

V = FESpace(model_in,reffe,conformity=:L2)
@test isa(V,ExtendedFESpace)

V = FESpace(model_in,reffe,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

# From a cell-wise vector of reffes

cell_reffe = ReferenceFE(model_in,lagrangian,Float64,order)

V = FESpace(model_in,cell_reffe)
@test isa(V,ExtendedFESpace)

V = FESpace(model_in,cell_reffe,conformity=:L2)
@test isa(V,ExtendedFESpace)

V = FESpace(model_in,cell_reffe,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

# From a CellFE

cell_fe = FiniteElements(PhysicalDomain(),model_in,lagrangian,Float64,order)

V = FESpace(model_in,cell_fe)
@test isa(V,ExtendedFESpace)

V = FESpace(model_in,cell_fe,conformity=:L2)
@test isa(V,ExtendedFESpace)

V = FESpace(model_in,cell_fe,conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

end # module


text:
module TestExamplesMPIP4estMesh3D

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "p4est_3d_dgsem")

@testset "P4estMesh MPI 3D" begin

# Run basic tests
@testset "Examples 3D" begin
  # Linear scalar advection
  @trixi_testset "elixir_advection_basic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
      # Expected errors are exactly the same as with TreeMesh!
      l2   = [0.00016263963870641478],
      linf = [0.0014537194925779984])
  end

  @trixi_testset "elixir_advection_restart.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
      l2   = [0.002590388934758452],
      linf = [0.01840757696885409])
  end

  @trixi_testset "elixir_advection_cubed_sphere.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_cubed_sphere.jl"),
      l2   = [0.002006918015656413],
      linf = [0.027655117058380085])
  end

  # Compressible Euler
  @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic.jl"),
      l2   = [0.0015106060984283647, 0.0014733349038567685, 0.00147333490385685, 0.001473334903856929, 0.0028149479453087093],
      linf = [0.008070806335238156, 0.009007245083113125, 0.009007245083121784, 0.009007245083102688, 0.01562861968368434],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.010380390326164493, 0.006192950051354618, 0.005970674274073704, 0.005965831290564327, 0.02628875593094754],
      linf = [0.3326911600075694, 0.2824952141320467, 0.41401037398065543, 0.45574161423218573, 0.8099577682187109],
      tspan = (0.0, 0.2),
      coverage_override = (polydeg=3,)) # Prevent long compile time in CI
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic_hohqmesh.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic_hohqmesh.jl"),
      l2   = [0.0042023406458005464, 0.004122532789279737, 0.0042448149597303616, 0.0036361316700401765, 0.007389845952982495],
      linf = [0.04530610539892499, 0.02765695110527666, 0.05670295599308606, 0.048396544302230504, 0.1154589758186293])
  end
end

end # P4estMesh MPI

end # module


text:
export problem2

"""
  problem2(N=4e6)

# Problem 2
## Even Fibonacci numbers

By considering the terms in the Fibonacci sequence whose values do not exceed four million, find the sum of the even-valued terms.
"""
function problem2(N=4e6)
  n = 0
  f = 0
  s = 0

  while f < N
    n += 3
    s += f
    f = fib(n)
  end
  return s
end

function fib(n::Int)
  n < 0 && return 0
  n < 2 && return n
  a = 1
  b = 0
  c = 0

  while n > 2
    n -= 1
    c = a
    a += b
    b = c
  end

  return a+b
end

max_stars_repo_name: halleysfifthinc/ProjectEuler

max_issues_repo_name: halleysfifthinc/ProjectEuler

max_forks_repo_name: halleysfifthinc/ProjectEuler


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((1+x^8)^(1/2)*(x+x^17+2*x^9)^-1*(1+2*x^8), x) == :(-1//4*(1+x^8)^-1//2+-1//4*arctanh((1+x^8)^(1/2)))
@test integrate((1+x^2)^-1, x) == :(arctan(x))
@test integrate(x^-1*(1+x^8)^-1//2, x) == :(-1//4*arctanh((1+x^8)^(1/2)))
@test integrate(x*(1+-1*x^3)^-1//2, x) == :(2*(1+-1*x^3)^(1/2)*(1+3^(1/2)+-1x)^-1+-1*3^1//4*((1+3^(1/2)+-1x)^-2*(1+-1x))^-1//2*((1+3^(1/2)+-1x)^-2*(1+x+x^2))^(1/2)*(1+-1*x^3)^-1//2*(2+-1*3^(1/2))^(1/2)*(1+-1x)*Elliptic.E(arcsin((1+3^(1/2)+-1x)^-1*(1+-1x+-1*3^(1/2))),-7+-4*3^(1/2))+2//3*2^(1/2)*3^3//4*((1+3^(1/2)+-1x)^-2*(1+-1x))^-1//2*((1+3^(1/2)+-1x)^-2*(1+x+x^2))^(1/2)*(1+-1*x^3)^-1//2*(1+-1x)*Elliptic.F(arcsin((1+3^(1/2)+-1x)^-1*(1+-1x+-1*3^(1/2))),-7+-4*3^(1/2)))
@test integrate(x^-1*(1+-1*x^3)^-1//2, x) == :(-2//3*arctanh((1+-1*x^3)^(1/2)))
@test integrate(x*(-71+x^4+-96x+10*x^2)^-1//2, x) == :(1//8*log(10001+x^8+-1408*x^3+-128*x^5+20*x^6+54*x^4+3124*x^2+(-71+x^4+-96x+10*x^2)^(1/2)*(781+x^6+-528x+-80*x^3+15*x^4+27*x^2)))
@test integrate(tan(x)^-2*(x+-1*tan(x)), x) == :(-1//2*x^2+-1*x*cot(x))
@test integrate(1+tan(x)^2+x*tan(x), x) == :((1/2)*im*x^2+(1/2)*im*Polylogarithms.polylog(2,-1*exp(2*im*x))+-1*x*log(1+exp(2*im*x))+tan(x))
@test integrate(x^-1*sin(x), x) == :(sinintegral(x))
@test integrate(x^-1*(x+exp(x))^-1//3*(3*(x+exp(x))^1//3+5*x^2+(2*x^2+3x)*exp(x)), x) == :(3*log(x)+3*x*(x+exp(x))^2//3)
@test integrate(x^-1+(x+log(x))^-3//2*(1+x^-1), x) == :(-2*(x+log(x))^-1//2+log(x))
@test integrate((x^3+x*log(x)^2+2*x^2*log(x))^-1*(x^2+log(x)^2+(x+log(x))^(1/2)*(1+x)+2*x*log(x)), x) == :(-2*(x+log(x))^-1//2+log(x))
@test integrate((log(x)^3+-1*x^2*log(x))^-1*(-1*x^2+-1*log(x)+2*log(x)^2), x) == :((1/2)*log(x+log(x))+-1//2*log(x+-1*log(x))+logintegral(x))
@test integrate((4+x^6+-5*x^4+5*x^2)^-1*(6+x^4+-3*x^2), x) == :(-1*arctan(3^(1/2)+-2x)+arctan((1/2)*x*(1+x^4+-3*x^2))+arctan(3^(1/2)+2x))

max_stars_repo_path: test/0 Independent test suites/Bronstein Problems.jl


text:
"""
	MultipleCounts()
	MultipleCounts(counts, proba)

Holds something like the photon counting probabilities with their respective
probability (in order to use them as a single observation). Can be declared
empty as a placeholder.

	Fields:
		- counts::Union{Nothing, Vector{ModeOccupation}, Vector{PartitionOccupancy}},
		- proba::Union{Nothing,Vector{Real}}
"""
mutable struct MultipleCounts

	counts::Union{Nothing, Vector{ModeOccupation}, Vector{PartitionOccupancy}}
	proba::Union{Nothing,Vector{Real}}

	MultipleCounts() = new(nothing,nothing)
	MultipleCounts(counts, proba) = new(counts,proba)

end

"""
	EventProbability(probability::Union{Nothing, Number})
	EventProbability(mc::MultipleCounts)

Holds the probability or probabilities of an [`Event`](@ref).

	Fields:
		- probability::Union{Number,Nothing, MultipleCounts}
		- precision::Union{Number,Nothing}
		- failure_probability::Union{Number,Nothing}
"""
mutable struct EventProbability
    probability::Union{Number,Nothing, MultipleCounts}
    precision::Union{Number,Nothing} # see remarks in conventions
    failure_probability::Union{Number,Nothing}

    function EventProbability(probability::Union{Nothing, Number})

        if probability == nothing
            new(nothing, nothing, nothing)
        else
            try
                probability = clean_proba(probability)
                new(probability,nothing,nothing)
            catch
                error("invalid probability")
            end
        end
    end

	function EventProbability(mc::MultipleCounts)

		try
			mc.proba = clean_pdf(mc.proba)
			new(mc,nothing,nothing)
		catch
			error("invalid probability")
		end

	end

end

"""
	Event{TIn<:InputType, TOut<:OutputMeasurementType}

Event linking an input to an output.

	Fields:
		- input_state::Input{TIn}
		- output_measurement::TOut
		- proba_params::EventProbability
		- interferometer::Interferometer
"""
struct Event{TIn<:InputType, TOut<:OutputMeasurementType}

    input_state::Input{TIn}
    output_measurement::TOut
    proba_params::EventProbability
    interferometer::Interferometer

    function Event{TIn, TOut}(input_state, output_measurement, interferometer::Interferometer, proba_params = nothing) where {TIn<:InputType, TOut<:OutputMeasurementType}

        if proba_params == nothing
            proba_params = EventProbability(nothing)
        end

        new{TIn,TOut}(input_state, output_measurement, proba_params, interferometer)
    end

    Event(i,o,interf,p = nothing) = Event{get_parametric_type(i)[1], get_parametric_type(o)[1]}(i,o,interf,p)

end

function check_probability_empty(ev::Event)
    if ev.proba_params.probability != nothing
		@warn "probability was already set in, rewriting"
	end
end


text:
function jitter!(ec::EChart, series::Int; pctxrange::Real = 0.05, pctyrange::Real = 0)

    #http://blogs.sas.com/content/iml/2011/07/05/jittering-to-prevent-overplotting-in-statistical-graphics.html
    #range method supposed to be from Graphical Methods for Data Analysis
    #Chambers, Cleveland, Kleiner, and Tukey (1983)

    #Calculate range, then use random pct of range
    _y = []
    _x = []

    for datum in ec.series[series].data
        push!(_x, datum[1])
        push!(_y, datum[2])
    end

    xrange = maximum(_x) - minimum(_x)
    yrange = maximum(_y) - minimum(_y)

    #To avoid inexact error, ensure array is float
    typeof(ec.series[series].data) == Array{Array{Int64,1},1} ?
        ec.series[1].data = [convert(Array{Float64,1}, x) for x in ec.series[series].data] :
            nothing

    #Iterate over series to jitter dataset
    for datum in ec.series[series].data

        datum[1] = rand(-pctxrange:0.001:pctxrange) * xrange + datum[1]
        datum[2] = rand(-pctyrange:0.001:pctyrange) * yrange + datum[2]

    end

    return ec
end

function jitter!(ec::EChart; pctxrange::Real = 0.05, pctyrange::Real = 0)

	for n in 1:length(ec.series)
		jitter!(ec, n, pctxrange = pctxrange, pctyrange = pctyrange)
	end

	return ec

end


text:
Base.print(G::FOTF) = show(G)

function Base.show(io::IO, G::FOTF)
    println(typeof(G))
    print("\n")
    numstr = polyshow(G.num, G.nn)
    denstr = polyshow(G.den, G.nd)
    print(numstr)
    print("\n")
    maxlengthstr = max(length(numstr), length(denstr))
    G.ioDelay == 0 ? println("-"^maxlengthstr) : println("-"^maxlengthstr*" exp(-$(G.ioDelay)*s)")
    print(denstr)
end

function polyshow(p, np)
    if typeof(p) <: Number
        return p
    end

    if length(np) == 0
        p=0
        np=0
    end

    ii = sortperm(np, rev=true)
    np = sort(np, rev=true)
    P = ""
    p = p[ii]

    for i = 1:length(p)
        P = P*"+"*"$(p[i])"*"s^"*"{$(np[i])}"
    end

    P = P[2:end]
    P = replace(P, "s^{0}" => "")
    P = replace(P, "s^{0.0}" => "")
    P = replace(P, "^{1}" => "")
    P = replace(P, "+-" => "-")
    P = replace(P, "+1s" => "+s")
    P = replace(P, "-1s" => "-s")

    return P
end


text: include("./WoodburyFactorizations.jl")


text:
# ---
# title: 777. Swap Adjacent in LR String
# id: problem777
# author: Indigo
# date: 2021-01-28
# difficulty: Medium
# categories: Brainteaser
# link: <https://leetcode.com/problems/swap-adjacent-in-lr-string/description/>
# hidden: true
# ---
#
# In a string composed of `'L'`, `'R'`, and `'X'` characters, like
# `"RXXLRXRXL"`, a move consists of either replacing one occurrence of `"XL"`
# with `"LX"`, or replacing one occurrence of `"RX"` with `"XR"`. Given the
# starting string `start` and the ending string `end`, return `True` if and only
# if there exists a sequence of moves to transform one string to the other.
#
#
#
# **Example 1:**
#
#
#
#     Input: start = "RXXLRXRXL", end = "XRLXXRRLX"
#     Output: true
#     Explanation: We can transform start to end following these steps:
#     RXXLRXRXL ->
#     XRXLRXRXL ->
#     XRLXRXRXL ->
#     XRLXXRRXL ->
#     XRLXXRRLX
#
#
# **Example 2:**
#
#
#
#     Input: start = "X", end = "L"
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: start = "LLR", end = "RRL"
#     Output: false
#
#
# **Example 4:**
#
#
#
#     Input: start = "XL", end = "LX"
#     Output: true
#
#
# **Example 5:**
#
#
#
#     Input: start = "XLLR", end = "LXLX"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= start.length <= 104`
#   * `start.length == end.length`
#   * Both `start` and `end` will only consist of characters in `'L'`, `'R'`, and `'X'`.
#
#
## @lc code=start
using LeetCode

function can_transform(s::String, t::String)
    i, j = 1, 1
    len = length(s)
    len != length(t) && return false
    while i <= len || j <= len
        while i <= len && s[i] == 'X'; i += 1; end
        while j <= len && t[j] == 'X' ; j += 1; end
        ((i <= len) ⊻ (j <= len)) && return false
        if i <= len && j <= len
            (s[i] != t[j] || (s[i] == 'L' && i < j) || (s[i] == 'R' && i > j)) && return false
        end
        i += 1; j += 1
    end
    return true
end
## @lc code=end


text:
module ConcaveHull

using NearestNeighbors
using RecipesBase
using LinearAlgebra

include("hull.jl")
include("concave_hull.jl")

export concave_hull, in_hull, area

@recipe function f(h::Hull)
    linewidth --> 2
    label --> ""
    v = copy(h.vertices)
    push!(v,v[1])
    x = [vv[1] for vv in v]
    y = [vv[2] for vv in v]
    x, y
end

end # module


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to the string splitting.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "String splitting" begin
    str = "Test 😅 \e[38;5;231;48;5;243mTest 😅 \e[38;5;201;48;5;243mTest\e[0m"

    result = split_string(str, -10)
    @test result == (
        "",
        "Test 😅 \e[38;5;231;48;5;243mTest 😅 \e[38;5;201;48;5;243mTest\e[0m"
    )

    result = split_string(str, 5)
    @test result == (
        "Test ",
        "😅 \e[38;5;231;48;5;243mTest 😅 \e[38;5;201;48;5;243mTest\e[0m"
    )

    result = split_string(str, 6)
    @test result == (
        "Test  ",
        "  \e[38;5;231;48;5;243mTest 😅 \e[38;5;201;48;5;243mTest\e[0m"
    )

    result = split_string(str, 6)
    @test result == (
        "Test  ",
        "  \e[38;5;231;48;5;243mTest 😅 \e[38;5;201;48;5;243mTest\e[0m"
    )

    result = split_string(str, 7)
    @test result == (
        "Test 😅",
        " \e[38;5;231;48;5;243mTest 😅 \e[38;5;201;48;5;243mTest\e[0m"
    )

    result = split_string(str, 8)
    @test result == (
        "Test 😅 \e[38;5;231;48;5;243m",
        "Test 😅 \e[38;5;201;48;5;243mTest\e[0m"
    )

    result = split_string(str, 13)
    @test result == (
        "Test 😅 \e[38;5;231;48;5;243mTest ",
        "😅 \e[38;5;201;48;5;243mTest\e[0m"
    )

    result = split_string(str, 14)
    @test result == (
        "Test 😅 \e[38;5;231;48;5;243mTest  ",
        "  \e[38;5;201;48;5;243mTest\e[0m"
    )

    result = split_string(str, 15)
    @test result == (
        "Test 😅 \e[38;5;231;48;5;243mTest 😅",
        " \e[38;5;201;48;5;243mTest\e[0m"
    )

    result = split_string(str, 1000)
    @test result == (
        "Test 😅 \e[38;5;231;48;5;243mTest 😅 \e[38;5;201;48;5;243mTest\e[0m",
        ""
    )
end


text:
# This is meant to be run with something like
# mpiexec -np 4 julia mpi.jl

using Test, MonteCarlo
import MPI

let
    function workload(x)
        sleep(x)
        return (MPI.Comm_rank(MPI.COMM_WORLD), x)
    end

    MPI.Init()
    N_workers = MPI.Comm_size(MPI.COMM_WORLD) - 1

    # Test if the scheduler generates this grouping:
    # first group:   [4, 7, ..., 7]
    # second group:  [11, 4, ..., 4]
    # third group:   [_, 4, ..., 4]
    # where each number corresponds to the sleep time of the relevant processes
    # More specifically checks if the "8" gets pairs with the initial "3" by checking
    # the total runtime

    inputs = fill(4, 3N_workers - 1)
    inputs[N_workers + 1] = 11
    inputs[2:N_workers] .= 7

    output = mpi_queue(workload, inputs, verbose = true)

    MPI.Barrier(MPI.COMM_WORLD)

    if MPI.Comm_rank(MPI.COMM_WORLD) == 0
        results = zeros(Int64, N_workers)
        for (id, val) in output
            results[id] += val
        end
        println("Time slept: $results")
        @testset "MPI Scheduler" begin
            @test all(results .== 15)
        end
    end
end


text:
## INSTRUCTIONS:
## RUN IN LINE WITH SIMULATION MAIN FLOW CODE
## REQUIRES "data_file"

### UPDATE FUNCTION FOR DATA FILE ARRAYS

function UpdateArray(data_dictionary,string,alterations)
  # access array
  #println(data_dictionary)
  #println(alterations)
  data_array = data_dictionary[string];
  # set modifications for all conditions here
  for i in 1:Int(size(alterations)[1])
    updates = alterations[Int(i)];
    index = updates[1];
    value = updates[2];
    data_array[Int(index)]=value;
  data_dictionary[string] = data_array;
  end
  return data_array
end

#### set values for multiple variables

# # # TGFB SIGNALING RATES
tgfb12_to_ap1_sp1_gain = 1.5
tgfb3_to_smad_gain_rate = 1.0 # TGFB3_R_phosphorylation_Smad: Smad =(TGFB3_R)=> pSmad
tgfb12_to_smad_rate = 0.2
MAPK_activation_smad_rate = 0.2
# # # BETA-CATEININ RATES
bcat_deg_rate = 0.1 # default 0.1
# bcat_deg_gsk3_rate = 0.05 # default 0.1

bcat_trans_rate = 1.0
pmapk_inact_gsk3_rate = 1.0
# # # BETA-CATEININ COMPLEXING
ecad_complex_gain = 10.0 # 6 or 8
ecad_complex_reverse = 1.0
ecad_bcat_degradation_rate = 0.1
lef1_complex_gain = 2.0
tcf4_complex_gain = 1.0
tgfb3_pik3_rate=1.5
apc_complex_forward_rate=8.0

#ecad_bais_factor = 1.0 # multiply ecadherin-bcat only
ecad_trans_rate = 1.5
ecad_deg_rate=0.5 #0.1

lef1_trans_rate = 1.0

snail_slug_trans_rate = 1.0 # increases lef1 mrna too much

### SELECTED DATA FILE INDEXES

# INDEX VALUES CHANGE WITH UPDATES TO MODEL.NET

# initial conditions
tgfb12_initial_cond_index = 88 # TGFB12
R_initial_cond_index = 87
vegf_initial_cond_index = 93 # VEGF

# TGFB signaling
tgfb3_binding_rate_index = 3 # TGFB3_binding: TGFB3+R =([])=> TGFB3_R
R_deg_rate_index = 55 # degradation_R: R =([])=> []
TGFB12_R_activation_RAS_rate_index = 15 # TGFB12_R_activation_RAS: RAS =(TGFB12_R)=> RAS_GTP
RAS_activation_RAF_rate_index  = 16 # RAS_GTP_phosphorylation_RAF: RAF =(RAS_GTP)=> RAF_P
RAF_activation_MAPK_rate_index = 17 # RAF_P_phosphorylation_MAPK: MAPK =(RAF_P)=> pMAPK
MAPK_activation_AP1_SP1_rate_index = 19 # MAPK_phosphorylation_AP1_SP1: AP1_SP1 =(pMAPK)=> AP1_SP1_P
MAPK_activation_smad_rate_index = 18 # MAPK_phosphorylation_Smad: Smad =(pMAPK)=> pSmad
# Bcatenin complex formation
bcat_ecad_forward_index = 20 # Binding_Ecadherin_Bcatenin: Ecadherin+Bcatenin =([])=> Ecadherin_Bcatenin
bcat_ecad_reverse_index = 21 # Binding_Ecadherin_Bcatenin: Ecadherin_Bcatenin =([])=> Ecadherin+Bcatenin (reverse)
#bcat_apc_reverse_index = 23
bcat_tcf4_forward_index = 24 # TCF4_Bcatenin_complex_formation: TCF4+Bcatenin =([])=> Active_TCF4
bcat_tcf4_reverse_index = 25 # TCF4_Bcatenin_complex_formation: Active_TCF4 =([])=> TCF4+Bcatenin (reverse)
bcat_lef1_forward_index = 26 # LEF1_Bcatenin_complex_formation: LEF1+Bcatenin =([])=> Active_LEF1
bcat_lef1_reverse_index = 27 # LEF1_Bcatenin_complex_formation: Active_LEF1 =([])=> LEF1+Bcatenin (reverse)
bcat_deg_rate_index = 41 # Bcatenin_degradation: Bcatenin =([])=> []
#bcat_deg_gsk3_mediated_index = 67 # Bcatenin_GSK3_MEDIATED_degradation: Bcatenin =(GSK3)=> []
#pmapk_inact_gsk3_index = 68 # 68	GSK3_deactivation_by_pMAPK: GSK3 =(pMAPK)=> []
bcat_trans_rate_index = 163 # translation_mRNA_Bcatenin: mRNA_Bcatenin = Bcatenin

ecad_bcat_degradation_rate_index = 56 # Ecadherin_Bcatenin_degradation: Ecadherin_Bcatenin =([])=> []

# VEGF signaing
PKC_activation_RAF_rate_index = 11 # PKC_activation_RAF: RAF =(Active_PKC)=> RAF_P

snail_slug_trans_rate_index = 179 # induction_gene_SNAIL_SLUG: gene_SNAIL_SLUG = mRNA_SNAIL_SLUG
lef1_trans_rate_index = 167 # translation_mRNA_LEF1: mRNA_LEF1 = LEF1

apc_complex_forward_rate_index = 22 # Binding_GSK3_APC_AXIN: GSK3_APC_AXIN+Bcatenin =([])=> GSK3_APC_AXIN_Bcatenin
tgfb3_pik3_rate_index = 14 # TGFB3_R_phosphorylation_PI3K: PI3K =(TGFB3_R)=> PI3K_P

#ecad bcat balance
ecad_trans_rate_index = 162 # translation_mRNA_Ecadherin: mRNA_Ecadherin = Ecadherin
ecad_deg_rate_index = 42 # Ecadherin_degrdation: Ecadherin =([])=> []

# TGFB dosing
tgfb12_forward_rate_index = 111 # addition_TGFb12: [] =([])=> TGFB12
tgfb12_reverse_rate_index = 112 # addition_TGFb12: TGFB12 =([])=> [] (reverse)
tgfb12_binding_rate_index = 1 # TGFB12_binding: TGFB12+R =([])=> TGFB12_R

# VEGF dosing
vegf_forward_rate_index = 113 # addition_VEGF: [] =([])=> VEGF
vegf_reverse_rate_index = 114 # addition_VEGF: VEGF =([])=> [] (reverse)
vegf_binding_rate_index = 5 # VEGF_binding: VEGF+VEGFR =([])=> VEGF_VEGFR

# TGFB3_Ab dosing
TGFB3_Ab_foward_rate_index = 104 # addition_TGFB3_Ab: [] =([])=> TGFB3_Ab (reverse)
TGFB3_Ab_reverse_rate_index = 103 # addition_TGFB3_Ab: TGFB3_Ab =([])=> []

# U0126 dosing
U0126_foward_rate_index = 102 # addition_U0126: [] =([])=> U0126 (reverse)
U0126_reverse_rate_index = 101 # addition_U0126: U0126 =([])=> []

# TGFB3 dosing
TGFB3_foward_rate_index = 115 # addition_TGFB3: [] =([])=> TGFB3
TGFB3_reverse_rate_index = 116 # addition_TGFB3: TGFB3 =([])=> [] (reverse)
tgfb12_to_smad_rate_index = 12 # TGFB12_R_phosphorylation_Smad: Smad =(TGFB12_R)=> pSmad
tgfb3_to_smad_gain_rate_index = 13 # TGFB3_R_phosphorylation_Smad: Smad =(TGFB3_R)=> pSmad

# Snail slug transfeciton
ss_virus_forward_rate_index = 106 # addition_ss_virus: [] =([])=> ss_virus (reverse)
ss_virus_reverse_rate_index = 107 # addition_ss_virus: ss_virus =([])=> []

# DN smad
# DN SMAD temoprary
#SMAD_induciton_rate_index = 139 # induction_gene_Smad: gene_Smad = mRNA_Smad
# reverse convention based on product term
smad_virus_forward_rate_index = 108 # addition_smad_virus: [] =([])=> smad_virus (reverse)
smad_virus_reverse_rate_index = 107 # addition_smad_virus: smad_virus =([])=> []

# DN LEF1
#lef1_induciton_rate_index = 145
# reverse convention based on product term
LEF1_virus_forward_rate_index = 110 # addition_lef1_virus: [] =([])=> lef1_virus (reverse)
LEF1_virus_reverse_rate_index = 109 # addition_lef1_virus: lef1_virus =([])=> []

yreg_deg_rate_index = 46 # degradation_mRNA_YREG1: mRNA_YREG1 = []

#### THESE SETTINGS ARE FOR STEADY STATE AND ARE DEFAULT ####

# initialize arrays to hold updates for each array in data_file
initial_cond_update_array = [];
rate_constant_update_array = [];
control_parameter_update_array = [];

# initial conditions
push!(initial_cond_update_array,[tgfb12_initial_cond_index,0.0]) # TGFB12 initial conditions
push!(initial_cond_update_array,[R_initial_cond_index,6.0]) # R initial condition
push!(initial_cond_update_array,[vegf_initial_cond_index,0.0]) # VEGF

# rates
# TGFB dosing
push!(rate_constant_update_array,[tgfb12_forward_rate_index,0.0]) # addition_TGFb12: [] =([])=> TGFB12
push!(rate_constant_update_array,[tgfb12_reverse_rate_index,0.0]) # addition_TGFb12: TGFB12 =([])=> [] (reverse)
push!(rate_constant_update_array,[tgfb12_binding_rate_index,3.0]) # TGFB12_binding: TGFB12+R =([])=> TGFB12_R

# VEGF dosing
push!(rate_constant_update_array,[vegf_forward_rate_index,0.0]) # addition_VEGF: [] =([])=> VEGF
push!(rate_constant_update_array,[vegf_reverse_rate_index,0.0]) # addition_VEGF: VEGF =([])=> [] (reverse)
push!(rate_constant_update_array,[vegf_binding_rate_index,3.0]) # VEGF_binding: VEGF+VEGFR =([])=> VEGF_VEGFR

# ecad deg and trans
push!(rate_constant_update_array,[ecad_deg_rate_index,ecad_deg_rate]) #default 0.1
push!(rate_constant_update_array,[ecad_trans_rate_index ,ecad_trans_rate]) #default 1.0

## Inhibitor dosing

# SNAIL_SLUG dosing
# push!(rate_constant_update_array,[snail_slug_forward_rate_index,0.0]) # addition_VEGF: [] =([])=> VEGF
# push!(rate_constant_update_array,[snail_slug_reverse_rate_index,0.0]) # addition_VEGF: VEGF =([])=> [] (reverse)

# TGFB3_Ab dosing
push!(rate_constant_update_array,[TGFB3_Ab_foward_rate_index,0.0])
push!(rate_constant_update_array,[TGFB3_Ab_reverse_rate_index,0.0])

# U0126 dosing
push!(rate_constant_update_array,[U0126_foward_rate_index ,0.0])
push!(rate_constant_update_array,[U0126_reverse_rate_index,0.0])

## Virus dosing

push!(rate_constant_update_array,[ss_virus_forward_rate_index,0.0])
push!(rate_constant_update_array,[ss_virus_reverse_rate_index,0.0])

#push!(rate_constant_update_array,[smad_virus_forward_rate_index,0.0])
#push!(rate_constant_update_array,[smad_virus_reverse_rate_index,0.0])

#push!(rate_constant_update_array,[LEF1_virus_forward_rate_index,0.0])
#push!(rate_constant_update_array,[LEF1_virus_reverse_rate_index,0.0])

# TGFB3 dosing

push!(rate_constant_update_array,[TGFB3_foward_rate_index,0.0])
push!(rate_constant_update_array,[TGFB3_reverse_rate_index,0.0])

##
# TGFB signaling
push!(rate_constant_update_array,[tgfb3_binding_rate_index,3.5]) # TGFB3_binding: TGFB3+R =([])=> TGFB3_R
push!(rate_constant_update_array,[tgfb3_to_smad_gain_rate_index,tgfb3_to_smad_gain_rate]) # TGFB3_binding: TGFB3+R =([])=> TGFB3_R

push!(rate_constant_update_array,[tgfb12_to_smad_rate_index,tgfb12_to_smad_rate]) # TGFB3_binding: TGFB3+R =([])=> TGFB3_R

push!(rate_constant_update_array,[R_deg_rate_index,0.01]) # R_degradation: R =([])=> []
push!(rate_constant_update_array,[TGFB12_R_activation_RAS_rate_index,2.0]) #TGFB12_R_activation_RAS: RAS =(TGFB12_R)=> RAS_GTP
push!(rate_constant_update_array,[RAS_activation_RAF_rate_index,tgfb12_to_ap1_sp1_gain]) # RAS_GTP_phosphorylation_RAF: RAF =(RAS_GTP)=> RAF_P
push!(rate_constant_update_array,[RAF_activation_MAPK_rate_index,tgfb12_to_ap1_sp1_gain]) # RAF_P_phosphorylation_MAPK: MAPK =(RAF_P)=> pMAPK
push!(rate_constant_update_array,[MAPK_activation_AP1_SP1_rate_index,tgfb12_to_ap1_sp1_gain]) # MAPK_phosphorylation_AP1_SP1: AP1_SP1 =(pMAPK)=> AP1_SP1_P
push!(rate_constant_update_array,[MAPK_activation_smad_rate_index,MAPK_activation_smad_rate]) # MAPK_phosphorylation_AP1_SP1: AP1_SP1 =(pMAPK)=> AP1_SP1_P

# VEGF signaing
push!(rate_constant_update_array,[PKC_activation_RAF_rate_index,0.1]) # PKC_activation_RAF: RAF =(Active_PKC)=> RAF_P

# Bcatenin complex formation
push!(rate_constant_update_array,[bcat_ecad_forward_index,ecad_complex_gain]) # 19	Binding_Ecadherin_Bcatenin: Ecadherin+Bcatenin =([])=> Ecadherin_Bcatenin
push!(rate_constant_update_array,[bcat_ecad_reverse_index,ecad_complex_reverse]) # 20	Binding_Ecadherin_Bcatenin: Ecadherin_Bcatenin =([])=> Ecadherin+Bcatenin (reverse)
push!(rate_constant_update_array,[bcat_tcf4_forward_index,tcf4_complex_gain]) # 21	TCF4_Bcatenin_complex_formation: TCF4+Bcatenin =([])=> Active_TCF4
push!(rate_constant_update_array,[bcat_tcf4_reverse_index,1.0]) # 22	TCF4_Bcatenin_complex_formation: Active_TCF4 =([])=> TCF4+Bcatenin (reverse)
push!(rate_constant_update_array,[bcat_lef1_forward_index,lef1_complex_gain]) # 23	LEF1_Bcatenin_complex_formation: LEF1+Bcatenin =([])=> Active_LEF1
push!(rate_constant_update_array,[bcat_lef1_reverse_index,1.0]) # 24	LEF1_Bcatenin_complex_formation: Active_LEF1 =([])=> LEF1+Bcatenin (reverse)

push!(rate_constant_update_array,[bcat_deg_rate_index,bcat_deg_rate]) # Bcatenin_degradation: Bcatenin =([])=> []
#push!(rate_constant_update_array,[bcat_deg_gsk3_mediated_index,bcat_deg_gsk3_rate])
#push!(rate_constant_update_array,[pmapk_inact_gsk3_index,pmapk_inact_gsk3_rate])

push!(rate_constant_update_array,[yreg_deg_rate_index,1.2])
push!(rate_constant_update_array,[bcat_trans_rate_index,bcat_trans_rate])
push!(rate_constant_update_array,[ecad_bcat_degradation_rate_index,ecad_bcat_degradation_rate])

# LEF1
push!(rate_constant_update_array,[lef1_trans_rate_index,lef1_trans_rate])

#
push!(rate_constant_update_array,[apc_complex_forward_rate_index,apc_complex_forward_rate])
push!(rate_constant_update_array,[tgfb3_pik3_rate_index,tgfb3_pik3_rate])

# control parameter array
#num_regs = 15; # number of inducing or repressing intereactions (could have multiple for species)
# do forced update for now
#for i in 1:num_regs
  #push!(control_parameter_update_array[i,1],GAIN_MASTER) # Gain for transcriptional activity
#end
#push!(control_parameter_update_array,[7,GAIN_TCF4_TGFB3])
#push!(control_parameter_update_array,[6,GAIN_BCAT_TGFB3])
#push!(control_parameter_update_array,[5,GAIN_ACTIVE_LEF1_ECAD])

# push modifications to RAM
#println(data_dictionary["INITIAL_CONDITION_ARRAY"])
#println(data_dictionary)
#UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
#UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)
#UpdateArray(data_dictionary,"CONTROL_PARAMETER_ARRAY",control_parameter_update_array)
###  I can't figure out how to get data_dictionary to be passed into this streamed file,
### It would seem that the stream opens a new variable space and then merges it with the main variable space
### Hence, these will need to be called post-stream

#


text:
@with_kw mutable struct LShapedData{T <: Real}
    Q::T = 1e10
    θ::T = -1e10
    iterations::Int = 0
end

@with_kw mutable struct LShapedParameters{T <: Real}
    τ::T = 1e-6
    log::Bool = true
end

struct LShaped{T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver} <: AbstractLShapedSolver{T,A,M,S}
    structuredmodel::JuMP.Model
    solverdata::LShapedData{T}

    # Master
    mastersolver::M
    mastervector::A
    c::A
    x::A
    Q_history::A

    # Subproblems
    nscenarios::Int
    subproblems::Vector{SubProblem{T,A,S}}
    subobjectives::A

    # Cuts
    θs::A
    cuts::Vector{SparseHyperPlane{T}}
    θ_history::A

    # Params
    parameters::LShapedParameters{T}
    progress::ProgressThresh{T}

    function (::Type{LShaped})(model::JuMP.Model,x₀::AbstractVector,mastersolver::AbstractMathProgSolver,subsolver::AbstractMathProgSolver; kw...)
        if nworkers() > 1
            warn("There are worker processes, consider using distributed version of algorithm")
        end
        length(x₀) != model.numCols && error("Incorrect length of starting guess, has ",length(x₀)," should be ",model.numCols)
        !haskey(model.ext,:SP) && error("The provided model is not structured")

        T = promote_type(eltype(x₀),Float32)
        c_ = convert(AbstractVector{T},JuMP.prepAffObjective(model))
        c_ *= model.objSense == :Min ? 1 : -1
        x₀_ = convert(AbstractVector{T},copy(x₀))
        mastervector = convert(AbstractVector{T},copy(x₀))
        A = typeof(x₀_)

        msolver = LQSolver(model,mastersolver)
        M = typeof(msolver)
        S = LQSolver{typeof(LinearQuadraticModel(subsolver)),typeof(subsolver)}
        n = StochasticPrograms.nscenarios(model)

        lshaped = new{T,A,M,S}(model,
                               LShapedData{T}(),
                               msolver,
                               mastervector,
                               c_,
                               x₀_,
                               A(),
                               n,
                               Vector{SubProblem{T,A,S}}(),
                               A(zeros(n)),
                               A(fill(-1e10,n)),
                               Vector{SparseHyperPlane{T}}(),
                               A(),
                               LShapedParameters{T}(;kw...),
                               ProgressThresh(1.0, "L-Shaped Gap "))
        # Initialize solver
        init!(lshaped,subsolver)
        return lshaped
    end
end
LShaped(model::JuMP.Model,mastersolver::AbstractMathProgSolver,subsolver::AbstractMathProgSolver; kw...) = LShaped(model,rand(model.numCols),mastersolver,subsolver; kw...)

function (lshaped::LShaped)()
    # Reset timer
    lshaped.progress.tfirst = lshaped.progress.tlast = time()
    # Start procedure
    while true
        status = iterate!(lshaped)
        if status != :Valid
            return status
        end
    end
end


text:
@testset "Demultiplexer" begin
    function randdna(n)
        return LongDNASeq(rand([DNA_A, DNA_C, DNA_G, DNA_T, DNA_N], n))
    end

    function make_errors(seq, p=0.03)
        seq = copy(seq)
        nucs = DNA['A', 'C', 'G', 'T', 'N']
        i = 1
        while i ≤ lastindex(seq)
            if rand() < p
                r = rand()
                if r < 0.1
                    # insertion
                    insert!(seq, i, rand(nucs))
                elseif r < 0.2
                    # deletion
                    deleteat!(seq, i)
                else
                    # substitution
                    seq[i] = rand(nucs)
                end
            end
            i += 1
        end
        return seq
    end

    @testset "Hamming distance" begin
        barcodes = [dna"ATGG", dna"CAGA", dna"GGAA", dna"TACG"]
        dplxr = Demultiplexer(barcodes, n_max_errors=1, distance=:hamming)

        for i in 1:lastindex(barcodes)
            @test dplxr[i] == barcodes[i]
        end

        @test demultiplex(dplxr, dna"ATGG") === (1, 0)
        @test demultiplex(dplxr, dna"CAGA") === (2, 0)
        @test demultiplex(dplxr, dna"GGAA") === (3, 0)
        @test demultiplex(dplxr, dna"TACG") === (4, 0)

        # every 1bp substitution is recoverable
        for (i, barcode) in enumerate(barcodes)
            # substitution
            for j in 1:lastindex(barcode)
                barcode′ = copy(barcode)
                for nt in [DNA_A, DNA_C, DNA_G, DNA_T, DNA_N]
                    barcode′[j] = nt
                    @test demultiplex(dplxr, barcode′) == (i, count(!=, barcode, barcode′))
                end
            end
        end

        n_ok = 0
        n_ok_with_fallback = 0
        for _ in 1:10_000
            i = rand(1:lastindex(barcodes))
            seq = make_errors(barcodes[i] * randdna(10))
            n_ok += demultiplex(dplxr, seq)[1] == i
            n_ok_with_fallback += demultiplex(dplxr, seq, true)[1] == i
        end
        # empirically, n_ok / 10_000 is ~0.985
        # @show n_ok
        @test n_ok / 10_000 > 0.98
        @test n_ok < n_ok_with_fallback
    end

    @testset "Levenshtein distance" begin
        barcodes = [dna"ATGG", dna"CAGA", dna"GGAA", dna"TACG"]
        dplxr = Demultiplexer(barcodes, n_max_errors=1, distance=:levenshtein)

        for i in 1:lastindex(barcodes)
            @test dplxr[i] == barcodes[i]
        end

        @test demultiplex(dplxr, dna"ATGG") === (1, 0)
        @test demultiplex(dplxr, dna"CAGA") === (2, 0)
        @test demultiplex(dplxr, dna"GGAA") === (3, 0)
        @test demultiplex(dplxr, dna"TACG") === (4, 0)

        # every 1bp substitution/insertion/deletion is recoverable
        for (i, barcode) in enumerate(barcodes)
            # substitution
            for j in 1:lastindex(barcode)
                barcode′ = copy(barcode)
                for nt in [DNA_A, DNA_C, DNA_G, DNA_T, DNA_N]
                    barcode′[j] = nt
                    i′, d = demultiplex(dplxr, barcode′)
                    @test i′ == i && 0 ≤ d ≤ 1
                end
            end
            # insertion
            for j in 1:lastindex(barcode)
                barcode′ = copy(barcode)
                insert!(barcode′, j, rand([DNA_A, DNA_C, DNA_G, DNA_T, DNA_N]))
                i′, d = demultiplex(dplxr, barcode′)
                @test i′ == i && 0 ≤ d ≤ 1
            end
            # deletion
            for j in 1:lastindex(barcode)
                barcode′ = copy(barcode)
                deleteat!(barcode′, j)
                i′, d = demultiplex(dplxr, barcode′)
                @test i′ == i && 0 ≤ d ≤ 1
            end
        end

        n_ok = 0
        n_ok_with_fallback = 0
        for _ in 1:10_000
            i = rand(1:lastindex(barcodes))
            seq = make_errors(barcodes[i] * randdna(10))
            n_ok += demultiplex(dplxr, seq)[1] == i
            n_ok_with_fallback += demultiplex(dplxr, seq, true)[1] == i
        end
        # empirically, n_ok / 10_000 is ~0.995
        # @show n_ok
        @test n_ok / 10_000 > 0.99
        @test n_ok < n_ok_with_fallback
    end
end


text:
# Unit testing for mathfuns.jl

using NumericExtensions
using Base.Test

@test sqr([1, 2, 3]) == [1, 4, 9]
@test_approx_eq rcp([1, 2, 3]) 1.0 / [1, 2, 3]
@test_approx_eq rsqrt([1, 2, 3]) 1.0 / sqrt([1, 2, 3])
@test_approx_eq rcbrt([1, 2, 3]) 1.0 / cbrt([1, 2, 3])

@test_approx_eq logsumexp(2.0, 3.0) log(exp(2.0) + exp(3.0))
@test_approx_eq logsumexp(3.0, 2.0) log(exp(2.0) + exp(3.0))


text:
# Which symbols to exclude from the export
const _EXCLUDE_SYMBOLS = []
function _hide_function(f::Function) # adds the given function to the excluded symbols
    push!(_EXCLUDE_SYMBOLS, Symbol(f))
end

# Create our module level logger (this will get precompiled)
const _LOGGER = Memento.getlogger(@__MODULE__)

"Suppresses information and warning messages output for PowerPlots, for fine grained control use the Memento package"
function silence()
    Memento.info(_LOGGER, "Suppressing information and warning messages for the rest of this session.  Use the Memento package for more fine-grained control of logging.")
    Memento.setlevel!(Memento.getlogger(PowerPlots), "error")
end
_hide_function(silence) # do not export, potential conflict of PowerModels.silence()

"allows the user to set the logging level without the need to add Memento"
function logger_config!(level)
    Memento.config!(Memento.getlogger("PowerPlots"), level)
end
_hide_function(logger_config!) # do not export, potential conflict of PowerModels.logger_config!()

function __init__()
    # Register the module level logger at runtime so that folks can access the logger via `getlogger(PowerPlots)`
    # NOTE: If this line is not included then the precompiled `PowerPlots._LOGGER` won't be registered at runtime.
    Memento.register(_LOGGER)
end


text:
const VarCostArgs = Union{Float64, NTuple{2, Float64}, Vector{NTuple{2, Float64}}}

abstract type OperationalCost <: DeviceParameter end

mutable struct VariableCost{T}
    cost::T
end

VariableCost(; cost) = VariableCost(cost)

function IS.serialize(val::VariableCost)
    return Dict{String, Any}("cost" => serialize(val.cost))
end

# The default implementation can't figure out the variable Union.

function IS.deserialize(::Type{VariableCost}, data::Dict)
    if data["cost"] isa Real
        return VariableCost(Float64(data["cost"]))
    elseif data["cost"][1] isa Array
        variable = Vector{Tuple{Float64, Float64}}()
        for array in data["cost"]
            push!(variable, Tuple{Float64, Float64}(array))
        end
    else
        @assert data["cost"] isa Tuple || data["cost"] isa Array
        variable = Tuple{Float64, Float64}(data["cost"])
    end

    return VariableCost(variable)
end

get_cost(vc::VariableCost) = vc.cost
Base.length(vc::VariableCost) = length(vc.cost)
Base.getindex(vc::VariableCost, ix::Int) = getindex(vc.cost, ix)

function get_breakpoint_upperbounds(vc::VariableCost{T}) where {T}
    throw(ArgumentError("Method not supported for VariableCost using $(T)"))
end

function get_slopes(vc::VariableCost{T}) where {T}
    throw(ArgumentError("Method not supported for VariableCost using $(T)"))
end

"""
Calculates the upper bounds of a variable cost function represented as a collection of piece-wise linear segments.

"""
function get_breakpoint_upperbounds(vc::VariableCost{Vector{NTuple{2, Float64}}})
    return get_breakpoint_upperbounds(get_cost(vc))
end

function get_breakpoint_upperbounds(vc::Vector{NTuple{2, Float64}})
    bp_ubs = Vector{Float64}(undef, length(vc))
    for (ix, component) in enumerate(vc)
        if ix == 1
            bp_ubs[ix] = component[2]
            continue
        end
        bp_ubs[ix] = component[2] - sum(bp_ubs[1:(ix - 1)])
    end
    return bp_ubs
end

"""
Calculates the slopes for the variable cost represented as a piece wise linear cost function. This function returns n - slopes for n - piecewise linear elements in the function. The first element of the return array corresponds to the average cost at the minimum operating point. If your formulation uses n -1 slopes, you can disregard the first component of the array. If the first point in the variable cost has a quantity of 0.0, the first slope returned will be 0.0, otherwise, the first slope represents the trajectory to get from the origin to the first point in the variable cost.

"""
function get_slopes(vc::VariableCost{Vector{NTuple{2, Float64}}})
    return get_slopes(get_cost(vc))
end

function get_slopes(vc::Vector{NTuple{2, Float64}})
    slopes = Vector{Float64}(undef, length(vc))
    previous = (0.0, 0.0)
    for (ix, component) in enumerate(vc)
        if ix == 1
            slopes[ix] = component[2] == 0.0 ? 0.0 : component[1] / component[2] # if the cost component starts at the y-origin make the first slope 0.0
            previous = component
            continue
        end
        slopes[ix] = (component[1] - previous[1]) / (component[2] - previous[2])
        previous = component
    end
    return slopes
end