text:
using StaticArrays, ForwardDiff, Requires

export DiscreteDS, DiscreteDS1D, evolve, evolve!, timeseries, dimension, jacobian

abstract type DiscreteDynamicalSystem <: DynamicalSystem end
#######################################################################################
#                                     Constructors                                    #
#######################################################################################
function test_functions(u0, eom, jac)
  length(size(u0)) == 1 || throw(ArgumentError("Initial condition must an AbstractVector"))
  D = length(u0)
  su0 = SVector{D}(u0); sun = eom(u0);
  length(sun) == length(u0) ||
  throw(DimensionMismatch("E.o.m. does not give same sized vector as initial condition"))
  if !issubtype((typeof(sun)), SVector)
    throw(ArgumentError("E.o.m. should create an SVector (from StaticArrays)"))
  end
  J1 = jac(u0); J2 = jac(SVector{length(u0)}(u0))
  if !issubtype((typeof(J1)), SMatrix) || !issubtype((typeof(J2)), SMatrix)
    throw(ArgumentError("Jacobian function should create an SMatrix (from StaticArrays)!"))
  end
  return true
end
function test_functions(u0, eom)
  jac = (x) -> ForwardDiff.jacobian(eom, x)
  test_discrete(u0, eom, fd_jac)
end

"""
    DiscreteDS(state, eom [, jacob]) <: DynamicalSystem
`D`-dimensional discrete dynamical system (used for `D ≤ 10`).
# Fields:
* `state::SVector{D}` : Current state-vector of the system, stored in the data format
  of `StaticArray`'s `SVector`.
* `eom::F` (function) : The function that represents the system's equations of motion
  (also called vector field). The function is of the format: `eom(u) -> SVector`
  which means that given a state-vector `u` it returns an `SVector` containing the
  next state.
* `jacob::J` (function) : A function that calculates the system's jacobian matrix,
  based on the format: `jacob(u) -> SMatrix` which means that given a state-vector
  `u` it returns an `SMatrix` containing the Jacobian at that state.
  If the `jacob` is not provided by the user, it is created with *tremendous* efficiency
  using the module `ForwardDiff`. Most of the time, for low dimensional systems, this
  Jacobian is within a few % of speed of a user-defined one.
"""
mutable struct DiscreteDS{D, T<:Real, F, J} <: DiscreteDynamicalSystem
  state::SVector{D,T}
  eom::F
  jacob::J
end
# constructor without jacobian (uses ForwardDiff)
function DiscreteDS(u0::AbstractVector, eom)
  su0 = SVector{length(u0)}(u0)
  @inline ForwardDiff_jac(x) = ForwardDiff.jacobian(eom, x)
  # test_functions(su0, eom, ForwardDiff_jac)
  return DiscreteDS(su0, eom, ForwardDiff_jac)
end
function DiscreteDS(u0::AbstractVector, eom, jac)
  su0 = SVector{length(u0)}(u0)
  # test_functions(su0, eom, jac)
  return DiscreteDS(su0, eom, jac)
end

"""
    DiscreteDS1D(state, eom [, deriv]) <: DynamicalSystem
One-dimensional discrete dynamical system.
# Fields:
* `state::Real` : Current state of the system.
* `eom::F` (function) : The function that represents the system's equation of motion:
  `eom(x) -> Real`.
* `deriv::D` (function) : A function that calculates the system's derivative given
  a state: `deriv(x) -> Real`. If it is not provided by the user
  it is created automatically using the module `ForwardDiff`.
"""
mutable struct DiscreteDS1D{S<:Real, F, D} <: DiscreteDynamicalSystem
  state::S
  eom::F
  deriv::D
end
function DiscreteDS1D(x0, eom)
  fd_deriv(x) = ForwardDiff.derivative(eom, x)
  DiscreteDS1D(x0, eom, fd_deriv)
end

dimension(::DiscreteDS{D, T, F, J})  where {D<:ANY, T<:ANY, F<:ANY, J<:ANY} = D
dimension(::DiscreteDS1D) = 1
jacobian(ds::DynamicalSystem) = ds.jacob(ds.state)
#######################################################################################
#                                 System Evolution                                    #
#######################################################################################
"""
```julia
evolve([state, ] ds::DynamicalSystem, T=1; diff_eq_kwargs = Dict()) -> new_state
```
Evolve a `state` (or the system's state) under the dynamics
of `ds` for total "time" `T`. For discrete systems `T` corresponds to steps and
thus it must be integer. Returns the final state after evolution.

The **keyword** argument `diff_eq_kwargs` (applicable only in `ContinuousDS`)
is a dictionary `Dict{Symbol, ANY}`
of keyword arguments
passed into the `solve` of the `DifferentialEquations.jl` package,
for example `Dict(:abstol => 1e-9)`.
If you want to specify a solver,
do so by using the symbol `:solver`, e.g.:
`Dict(:solver => DP5(), :maxiters => 1e9)`. This requires you to have been first
`using OrdinaryDiffEq` or `using DifferentialEquations` to access the solvers.

This function *does not store* any information about intermediate steps.
Use `timeseries` if you want to produce timeseries of the system.
"""
function evolve(ds::DiscreteDynamicalSystem, N::Int = 1)
  st = ds.state
  st = evolve(st, ds, N)
end
function evolve(state, ds::DiscreteDynamicalSystem, N::Int = 1)
  f = ds.eom
  for i in 1:N
    state = f(state)
  end
  return state
end

"""
```julia
evolve!(ds::DynamicalSystem, T; diff_eq_kwargs = Dict()) -> ds
```
Evolve (in-place) a dynamical system for total "time" `T`, setting the final
state as the system's state.
"""
function evolve!(ds::DiscreteDynamicalSystem, N::Int = 1)
  st = ds.state
  ds.state = evolve(st, ds, N)
  return ds
end

"""
```julia
timeseries(ds::DynamicalSystem, T; kwargs...)
```
Create a matrix that will contain the timeseries of the sytem, after evolving it
for time `T` (`D` is the system dimensionality). *Each column corresponds to
one dynamic variable.*

For the discrete case, `T` is an integer and a `T×D` matrix is returned. For the
continuous case, a `K×D` matrix is returned, with `K = length(0:dt:T)` with
`0:dt:T` representing the time vector.
# Keywords:
* `mutate = true` : whether to update the dynamical system's state with the
  final state of the timeseries.
* `dt = 0.05` : (only for continuous) Time step of value output during the solving
  of the continuous system.
* `diff_eq_kwargs = Dict()` : (only for continuous) A dictionary `Dict{Symbol, ANY}`
  of keyword arguments
  passed into the `solve` of the `DifferentialEquations.jl` package,
  for example `Dict(:abstol => 1e-9)`. If you want to specify a solver,
  do so by using the symbol `:solver`, e.g.:
  `Dict(:solver => DP5(), :maxiters => 1e9)`. This requires you to have been first
  `using OrdinaryDiffEq` to access the solvers.
"""
function timeseries(ds::DiscreteDS, N::Real; mutate = true)
  st = ds.state
  T = eltype(st)
  D = length(st)
  ts = Array{T}(N, D)
  f = ds.eom
  ts[1,:] .= ds.state
  for i in 2:N
    st = f(st)
    ts[i, :] .= st
  end
  if mutate
    ds.state = ts[end, :]
  end
  return ts
end

function timeseries(ds::DiscreteDS1D, N::Int; mutate = true)
  x = deepcopy(ds.state)
  f = ds.eom
  ts = Vector{eltype(x)}(N)
  ts[1] = x
  for i in 2:N
    x = f(x)
    ts[i] = x
  end
  if mutate
    ds.state = x
  end
  return ts
end

#######################################################################################
#                                 Pretty-Printing                                     #
#######################################################################################
import Base.show
function Base.show(io::IO, s::DiscreteDS{N, S, F, J}) where {N<:ANY, S<:ANY, F<:ANY, J<:ANY}
  print(io, "$N-dimensional discrete dynamical system:\n",
  "state: $(s.state)\n", "e.o.m.: $F\n", "jacobian: $J")
end

@require Juno begin
  function Juno.render(i::Juno.Inline, s::DiscreteDS{N, S, F, J}) where {N<:ANY, S<:ANY, F<:ANY, J<:ANY}
    t = Juno.render(i, Juno.defaultrepr(s))
    t[:head] = Juno.render(i, Text("$N-dimensional discrete dynamical system"))
    t
  end
end

# 1-D
function Base.show(io::IO, s::DiscreteDS1D{S, F, J}) where {S<:ANY, F<:ANY, J<:ANY}
  print(io, "1-dimensional discrete dynamical system:\n",
  "state: $(s.state)\n", "e.o.m.: $F\n", "jacobian: $J")
end
@require Juno begin
  function Juno.render(i::Juno.Inline, s::DiscreteDS1D{S, F, J}) where {S<:ANY, F<:ANY, J<:ANY}
    t = Juno.render(i, Juno.defaultrepr(s))
    t[:head] = Juno.render(i, Text("1-dimensional discrete dynamical system"))
    t
  end
end


text:
export heisenberg, transverse_ising

"""
    heisenberg(nbit::Int; periodic::Bool=true)

1D Heisenberg hamiltonian defined as ``\\sum_{i=1}^{n} X_{i}X_{i+1} + Y_{i}Y_{i+1} + Z_{i}Z_{i+1}``, where ``n`` is specified by `nbit`.
`periodic` means the boundary condition is periodic.

References
----------------------
* de Oliveira, Mário J. "Ground-state properties of the spin-1/2 antiferromagnetic Heisenberg chain obtained by use of a Monte Carlo method." Physical Review B 48.9 (1993): 6141-6143.
"""
function heisenberg(nbit::Int; periodic::Bool=true)
    map(1:(periodic ? nbit : nbit-1)) do i
        j=i%nbit+1
        repeat(nbit,X,(i,j)) + repeat(nbit, Y, (i,j)) + repeat(nbit, Z, (i,j))
    end |> sum
end

"""
    transverse_ising(nbit::Int, h::Number; periodic::Bool=true)

1D transverse Ising hamiltonian defined as ``\\sum_{i=1}^{n} hX_{i} + Z_{i}Z_{i+1}``, where ``n`` is specified by `nbit`.
`periodic` means the boundary condition is periodic.
"""
function transverse_ising(nbit::Int, h::Number; periodic::Bool=true)
    ising_term = map(1:(periodic ? nbit : nbit-1)) do i
        repeat(nbit,Z,(i,i%nbit+1))
    end |> sum
    ising_term + h*sum(map(i->put(nbit,i=>X), 1:nbit))
end

ext: jl


text:
"""
    box_fill_outline(pt_topleft, colfill;
                     height = row_height(), width = EM, luminfac = 2,
                     opacity = missing)

Draw a filled box.
Default height adapts to the current 'Top API' font, which is used for text tables.
Outline color luminosity is 200% the luminosity of fill color, if possible.
For dark backgrounds, consider luminfac = 0.5.
By default, opacity is not changed from the current state. Range 0-1, where 1 is opaque.
"""
function box_fill_outline(pt_topleft, colfill;
                          height = row_height(), width = EM, luminfac = 2, boxopacity = missing)
    gsave()
    !ismissing(boxopacity) && setopacity(boxopacity)
    sethue(colfill)
    box(pt_topleft, pt_topleft + (width, height), :fill)
    # outline
    lu = get_current_lumin()
    sethue(color_with_lumin(colfill, min(100, luminfac * lu)))
    !ismissing(boxopacity) && setopacity(boxopacity)
    box(pt_topleft, pt_topleft + (width, height), :stroke)
    grestore()
end
"""
label_boxed(pt, str; colfill = color_with_lumin(get_current_RGBA(), 0.5 * get_current_lumin()),
               luminfac = 1.0, boxopacity = 0.75, valign = :bottom)

Using the 'Toy API. For keyword arguments, see 'box_fill_outline'. Could be extended with 'text' arguments.
"""
function label_boxed(pt, str; colfill = color_with_lumin(get_current_RGBA(), 0.5 * get_current_lumin()),
               luminfac = 1.0, boxopacity = 0.75, valign = :bottom)
    @layer begin
        box_fill_outline(pt, colfill; luminfac, height = -FS / 2, width = pixelwidth(str), boxopacity)
        #sethue(PALETTE[10])
        text(str, pt; valign)
    end
end

"""
    setlabel(endpoint, vx, vy, α)

Add a label at the end point, rounded absolute length with unit.
"""
function setlabel(endpoint, vx, vy, α)
    Δvx = (vx * cos(α), vx *sin(α))
    Δvy = (-vy * sin(α), vy *cos(α))
    Δv = Δvx + Δvy
    l = hypot(Δv)
    rounded = round(unit(l), l, digits = 1)
    direction = atan(Δv[2], Δv[1]) |> °
    labdir = direction + 90°
    labpos = if -90° <= labdir < -70°
        :S
    elseif labdir < -20°
        :SE
    elseif labdir < 20°
        :E
    elseif labdir < 70°
        :NE
    elseif labdir < 110°
        :N
    elseif labdir < 160°
        :NW
    elseif labdir < 200°
        :W
    elseif labdir < 250°
        :SW
    else
        :S
    end
    label(string(rounded) , labpos, endpoint, offset = div(EM, 4))
end


text:
# Autogenerated wrapper script for libjulia_jll for armv7l-linux-musleabihf-cxx11-julia_version+1.8.0
export libjulia

using LibUnwind_jll
using LibUV_jll
JLLWrappers.@generate_wrapper_header("libjulia")
JLLWrappers.@declare_library_product(libjulia, "libjulia.so.1")
function __init__()
    JLLWrappers.@generate_init_header(LibUnwind_jll, LibUV_jll)
    JLLWrappers.@init_library_product(
        libjulia,
        "lib/libjulia.so",
        nothing,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module TestInterpreter
    using Base.Test
    using Brainfuck

    type MyIO <: IO
        bytes::Vector{Uint8}
    end
    Base.write(io::MyIO, val::Uint8) = push!(io.bytes, val)
    Base.read(io::MyIO, ::Type{Uint8}) = shift!(io.bytes)

    text = readall(Pkg.dir("Brainfuck", "test", "programs", "print.bfk"))
    ops = Brainfuck.lexandparse(text)

    io1 = MyIO(Uint8[])
    write(io1, 'a')
    write(io1, 'b')
    io2 = IOBuffer()
    Brainfuck.interpret(ops, io_in = io1, io_out = io2)
    @test takebuf_array(io2) == [0x62, 0x61, 0x63, 0x62]

    io1 = MyIO(Uint8[])
    write(io1, 'a')
    write(io1, 'b')
    io2 = IOBuffer()
    Brainfuck.brainfuck(text, io_in = io1, io_out = io2)
    @test takebuf_array(io2) == [0x62, 0x61, 0x63, 0x62]
end


text:
using Test, LinearAlgebra, SparseArrays, Random
using QDLDL
rng = Random.MersenneTwister(1312)

@testset "non-quasidefinite" begin

  m = 20

  A = random_psd(m)
  A[:,10] .= 0
  A[10,:] .= 0

  @test_throws ErrorException qdldl(A)

end

nothing


text:
"""
   function create_interpolation_grid(xbnd, xres, thegrid)
Given boundaries of a domain, and a resolution in each direction, create
an interpolation grid. This is where the interpolation functions for
a set of variables will be evaluated.
"""
function create_interpolation_grid(xbnd, xres, thegrid)
    x1g = collect(range(xbnd[1, 1], xbnd[2, 1], step = xres[1]))
    x2g = collect(range(xbnd[1, 2], xbnd[2, 2], step = xres[2]))
    x3g = collect(range(xbnd[1, 3], xbnd[2, 3], step = xres[3]))
    intrp_brck = ClimateMachine.InterpolationBrick(thegrid, xbnd, x1g, x2g, x3g)
    return intrp_brck
end

"""
   function interpolate_variables(objects, brick)
Create an interpolation function from data and evaluate that
function on the interpolation brick passed.
"""
function interpolate_variables(objects, brick)
    i_objects = []
    for object in objects
        nvars = size(object.data, 2)
        i_object = Array{FT}(undef, brick.Npl, nvars)
        ClimateMachine.interpolate_local!(brick, object.data, i_object)
        push!(i_objects, i_object)
    end

    return i_objects
end


text:
using Test

@time begin
    @time include("test0.jl")
    @time include("test2.jl")
    @time include("test_limiter.jl")
end


text:
"
    MPSComoving(leftstate,window,rightstate)

    muteable window of tensors on top of an infinite chain
"
mutable struct MPSComoving{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbstractMPS
    left_gs::InfiniteMPS{A,B}
    window::FiniteMPS{A,B}
    right_gs::InfiniteMPS{A,B}

    function MPSComoving(left::InfiniteMPS{A,B},window::FiniteMPS{A,B},right::InfiniteMPS{A,B}) where {A<:GenericMPSTensor,B<:MPSBondTensor}
        #todo:insert checks
        return new{A,B}(left,window,right);
    end
end

#allow construction with an array of tensors
MPSComoving(left::InfiniteMPS,site_tensors::Array,right::InfiniteMPS) = MPSComoving(left,FiniteMPS(site_tensors),right)

# allow construction with one large tensorkit space
MPSComoving(f, elt,P::ProductSpace, args...; kwargs...) = MPSComoving(f, elt, collect(P), args...; kwargs...)

# allow construction given only a physical space and length
MPSComoving(f,elt, N::Int, V::VectorSpace, args...; kwargs...) = MPSComoving(f, elt,fill(V, N), args...; kwargs...)

function MPSComoving(f, elt, physspaces::Vector{<:Union{S,CompositeSpace{S}}}, maxvirtspace::S,
                    leftgs::M,
                    rightgs::M) where {S<:ElementarySpace,M<:InfiniteMPS}

    left = virtualspace(leftgs,0);
    right = virtualspace(rightgs,length(physspaces));
    window = FiniteMPS(f,elt,physspaces,maxvirtspace,left=left,right=right);
    MPSComoving(leftgs,window,rightgs)
end
function MPSComoving(f,elt,
                    physspaces::Vector{<:Union{S,CompositeSpace{S}}},
                    virtspaces::Vector{S},
                    leftgs::M,
                    rightgs::M) where {S<:ElementarySpace,M<:InfiniteMPS}
    N = length(physspaces)
    length(virtspaces) == N+1 || throw(DimensionMismatch())

    tensors = [TensorMap(f, elt,virtspaces[n] ⊗ physspaces[n], virtspaces[n+1]) for n=1:N]

    return MPSComoving(leftgs,FiniteMPS(tensors,overwrite=true),rightgs)
end

#take a window from an infinite mps, and use that as state
function MPSComoving(state::InfiniteMPS{A,B},len::Int) where {A,B}
    CLs = Vector{Union{Missing,B}}(missing,len+1)
    ALs = Vector{Union{Missing,A}}(missing,len)
    ARs = Vector{Union{Missing,A}}(missing,len)
    ACs = Vector{Union{Missing,A}}(missing,len)

    ALs.= state.AL[1:len];
    ARs.= state.AR[1:len];
    ACs.= state.AC[1:len];
    CLs.= state.CR[0:len];

    MPSComoving(state,FiniteMPS(ALs,ARs,ACs,CLs),state)
end

Base.copy(state::MPSComoving{A,B}) where {A,B} = MPSComoving(state.left_gs,copy(state.window),state.right_gs);

Base.length(state::MPSComoving) = length(state.window)
Base.size(psi::MPSComoving, i...) = size(psi.window, i...)

Base.eltype(::Type{MPSComoving{Mtype,Vtype}}) where {Mtype<:GenericMPSTensor,Vtype<:MPSBondTensor} = Mtype
site_type(::Type{MPSComoving{Mtype,Vtype}}) where {Mtype<:GenericMPSTensor,Vtype<:MPSBondTensor} = Mtype
bond_type(::Type{MPSComoving{Mtype,Vtype}}) where {Mtype<:GenericMPSTensor,Vtype<:MPSBondTensor} = Vtype

TensorKit.space(psi::MPSComoving{<:MPSTensor}, n::Integer) = space(psi.AC[n], 2)
virtualspace(psi::MPSComoving, n::Integer) =
    n < length(psi) ? _firstspace(psi.AC[n+1]) : dual(_lastspace(psi.AC[n]))

r_RR(state::MPSComoving)=r_RR(state.right_gs,length(state))
l_LL(state::MPSComoving)=l_LL(state.left_gs,1)

function Base.getproperty(psi::MPSComoving,prop::Symbol)
    if prop == :AL
        return ALView(psi)
    elseif prop == :AR
        return ARView(psi)
    elseif prop == :AC
        return ACView(psi)
    elseif prop == :CR
        return CRView(psi)
    else
        return getfield(psi,prop)
    end
end

max_Ds(f::MPSComoving) = max_Ds(f.window)

Base.:*(psi::MPSComoving, a::Number) = rmul!(copy(psi),a)
Base.:*(a::Number, psi::MPSComoving) = lmul!(a,copy(psi))

function TensorKit.lmul!(a::Number, psi::MPSComoving)
    psi.window.ACs .*=a;
    psi.window.CLs .*=a;
    return psi
end

function TensorKit.rmul!(psi::MPSComoving,a::Number)
    psi.window.ACs .*=a;
    psi.window.CLs .*=a;
    return psi
end

function TensorKit.dot(psi1::MPSComoving, psi2::MPSComoving)
    length(psi1) == length(psi2) || throw(ArgumentError("MPS with different length"))
    psi1.left_gs == psi2.left_gs || throw(ArgumentError("left InfiniteMPS is different"))
    psi1.right_gs == psi2.right_gs || throw(ArgumentError("right InfiniteMPS is different"))

    ρr = transfer_right(r_RR(psi2),psi2.AR[2:end],psi1.AR[2:end]);
    return tr(_transpose_front(psi1.AC[1])' * _transpose_front(psi2.AC[1]) * ρr)
end

TensorKit.norm(psi::MPSComoving) = norm(psi.AC[1])

TensorKit.normalize!(psi::MPSComoving) = rmul!(psi, 1/norm(psi))
TensorKit.normalize(psi::MPSComoving) = normalize!(copy(psi))


text:
module uJSON
    using LegacyStrings

	const ujsonlib = Libdl.find_library(["ujsonlib"],[Pkg.dir("uJSON", "deps")])
	const TYPES = Any # Union(Dict, Array, String, Number, Bool, Nothing) # Types it may encounter
	const KEY_TYPES = Union{UTF32String} # Types it may encounter as object keys

	type UltraObject
	    # array to put JSON data in
	    array::Array{TYPES, 1}
	    # array to hold keys in
	    route::Array{TYPES, 1}
	    # whether we are in a dict (as opposed to an array)
	    in_dict::Bool
	    # working object used as a reference to where we're putting data now
	    working_obj::TYPES
	    UltraObject() = new(TYPES[], TYPES[], false, nothing)
	end

	function set_last!{T}(uo::UltraObject, value::T, key::AbstractString)
	    if uo.in_dict
	        setindex!(uo.working_obj, value, key)
	    else
	        push!(uo.working_obj, value)
	    end
	end

	function get_string(key_::Ptr{Int32}, key_length_::Ptr{Int32})
        utf32(key_, unsafe_load(key_length_))
	end

	function get_key(uobj_::Ptr{Void}, key_::Ptr{Int32}, key_length_::Ptr{Int32})
	    uo = unsafe_pointer_to_objref(uobj_)::UltraObject
	    # seems to be slightly quicker if we don't call get_string
        key = uo.in_dict ?  utf32(key_, unsafe_load(key_length_)): ""
		return uo, key
	end

	function startnew(uobj_::Ptr{Void},
	                  key_::Ptr{Int32},
	                  key_length_::Ptr{Int32},
	                  is_dict_::Ptr{Int32})
	    uo, key = get_key(uobj_, key_, key_length_)

	    is_dict = Bool(unsafe_load(is_dict_))
	    new_item = is_dict ? Dict{KEY_TYPES, TYPES}() : TYPES[]

	    if uo.working_obj == nothing
	        uo.working_obj = uo.array
	    end
	    set_last!(uo, new_item, key)
	    push!(uo.route, new_item)
	    uo.in_dict = is_dict
	    uo.working_obj = new_item
	    return nothing
	end
	const startnew_c = cfunction(startnew, Void, (Ptr{Void},
	                                              Ptr{Int32},
	                                              Ptr{Int32},
	                                              Ptr{Int32}))

	function exitob(uobj_::Ptr{Void})
	    uo = unsafe_pointer_to_objref(uobj_)::UltraObject
	    pop!(uo.route)
	    uo.working_obj = length(uo.route) > 0 ? last(uo.route) : nothing
	    uo.in_dict = isa(uo.working_obj, Dict)
	    return nothing
	end
	const exitob_c = cfunction(exitob, Void, (Ptr{Void},))

	# null, Bool, int
	function addnbi(uobj_::Ptr{Void},
	                key_::Ptr{Int32},
	                key_length_::Ptr{Int32},
	                value_::Ptr{Int64},
	                value_type_::Ptr{Int32})
	    uo, key = get_key(uobj_, key_, key_length_)

	    value_type = unsafe_load(value_type_)
	    if value_type == -1
	        value = nothing
	    elseif value_type == 0
	        value = Bool(unsafe_load(value_))
	    elseif value_type == 1
	        value = unsafe_load(value_)
	    else
	        error("value type unknown: ", value_type)
	    end
	    set_last!(uo, value, key)
	    return nothing
	end
	const addnbi_c = cfunction(addnbi,
	                           Void,
	                           (Ptr{Void},
	                            Ptr{Int32},
	                            Ptr{Int32},
	                            Ptr{Int64},
	                            Ptr{Int32}))

	function adddouble(uobj_::Ptr{Void},
	                   key_::Ptr{Int32},
	                   key_length_::Ptr{Int32},
	                   value_::Ptr{Float64})
	    uo, key = get_key(uobj_, key_, key_length_)

	    value = unsafe_load(value_)::Float64
	    set_last!(uo, value, key)
	    return nothing
	end
	const adddouble_c = cfunction(adddouble, Void, (Ptr{Void}, Ptr{Int32}, Ptr{Int32}, Ptr{Float64}))

	function addstring(uobj_::Ptr{Void},
	                   key_::Ptr{Int32},
	                   key_length_::Ptr{Int32},
	                   value_::Ptr{Int32},
	                   value_length_::Ptr{Int32})
	    uo, key = get_key(uobj_, key_, key_length_)

	    value = get_string(value_, value_length_)
	    set_last!(uo, value, key)
	    return nothing
	end
	const addstring_c = cfunction(addstring,
	                              Void,
	                              (Ptr{Void},
	                               Ptr{Int32},
	                               Ptr{Int32},
	                               Ptr{Int32},
	                               Ptr{Int32}))

# not used as reading the string in is
	# function parsefile(filename::String)
	#     uo = UltraObject()

	#     result = ccall( (:process_file, ujsonlib),
	#                     Int32,
	#                     (Ptr{UInt8}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Any),
	#                     pointer(filename), startnew_c, exit_ob_c, addnbi_c, add_double_c, add_string_c, uo)
	#     # the first (and only) item in the base array is the actual data structure
	#     if result != 1
	#         error("error processing JSON")
	#     end
	#     if length(uo.array) == 0
	#         error("no JSON found")
	#     end
	#     return uo.array[1]
	# end

    function parsefile(filename::String)
        parse(readstring(filename))
    end

	function parse(io::IO)
		str = read(io, String)
		return parse(str)
	end

	function parse(str::String)
	    uo = UltraObject()
	    result = ccall( (:process_string, ujsonlib),
	                    Int32,
	                    (Ptr{UInt8}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Any),
	                    str, startnew_c, exitob_c, addnbi_c, adddouble_c, addstring_c, uo)
	    # the first (and only) item in the base array is the actual data structure
	    if result != 1
	        error("error processing JSON")
	    end
	    if length(uo.array) == 0
	        error("no JSON found")
	    end
	    return uo.array[1]
	end
end


text:
using AWSBatch: JobState, SUBMITTED, PENDING, RUNNABLE, STARTING, RUNNING, SUCCEEDED, FAILED

@testset "JobState" begin
    @testset "parse" begin
        @test length(instances(JobState)) == 7
        @test parse(JobState, "SUBMITTED") == SUBMITTED
        @test parse(JobState, "PENDING") == PENDING
        @test parse(JobState, "RUNNABLE") == RUNNABLE
        @test parse(JobState, "STARTING") == STARTING
        @test parse(JobState, "RUNNING") == RUNNING
        @test parse(JobState, "SUCCEEDED") == SUCCEEDED
        @test parse(JobState, "FAILED") == FAILED
    end

    @testset "order" begin
        @test SUBMITTED < PENDING < RUNNABLE < STARTING < RUNNING < SUCCEEDED < FAILED
    end
end


text:
module CYCLOPS_v6_2a_MultiCoreModule_Smooth

using StatsBase
using MultivariateStats
using CYCLOPS_v6_2a_AutoEncoderModule_multi
using CYCLOPS_v6_2a_PreNPostprocessModule

export backgroundmetrics_global_eigen
export multicore_backgroundmetrics_global_eigen
export multicore_backgroundstatistics_global_eigen
export smoothness_measures, circ_diff, circ_diff_phases

###########################################################################
#  Note this is written to take adventage of a 6 core device
###########################################################################

function circ_diff(data::Array{Float64,2})
       	circd=hcat(diff(data,2),data[:,1]-data[:,end])
       	circd
end

###########################################################################
function circ_diff_phases(data::Array{Float64,1})
       	circd=[diff(data),data[1]+ (2*pi-data[end])]
       	circd
end

###########################################################################
###########################################################################
###########################################################################
function circ_diff_phases2(data)
	ans=zeros(length(data))

	ans[1]=((2*pi+data[2])-data[end])/2

		for i in 2:length(data)-1
			ans[i]=(data[i+1]-data[i-1])/2
		end

	ans[end]=((2*pi+data[1])-data[end-1])/2

	ans
end

###########################################################################
function smoothness_measures(l_seeddata::Array{Float64,2},l_eigendata::Array{Float64,2},estimated_phaselist)
	estimated_phaselist1s=mod.(estimated_phaselist,2*pi)
	use_order_c=sortperm(estimated_phaselist1s)
	use_order_l=sortperm(vec(l_eigendata[1,:]))

	new_phaselist=vec(estimated_phaselist1s[use_order_c])
	circ_seeddata=l_seeddata[:,use_order_c]
	circ_eigendata=l_eigendata[:,use_order_c]

	lin_seeddata=l_seeddata[:,use_order_l]
	lin_eigendata=l_eigendata[:,use_order_l]

	nsamp=size(l_seeddata)[2]

	gdiffs				=circ_diff(circ_seeddata)
	gdiffs2				=gdiffs .* gdiffs
	gdiffsm				=sqrt.(sum(gdiffs2,1))
	num					=sum(gdiffsm)/nsamp

	gdiffs				=diff(lin_seeddata,2)
	gdiffs2				=gdiffs .* gdiffs
	gdiffsm				=sqrt.(sum(gdiffs2,1))
	denom				=sum(gdiffsm)/(nsamp-1)
	measure_raw			=num/denom

	######################################################
	gdiffs				=circ_diff(circ_seeddata)
	gdiffs2				=gdiffs .* gdiffs
	gdiffsm				=(sum(gdiffs2,1))
	num					=sum(gdiffsm)/nsamp

	gdiffs				=diff(lin_seeddata,2)
	gdiffs2				=gdiffs .* gdiffs
	gdiffsm				=(sum(gdiffs2,1))
	denom				=sum(gdiffsm)/(nsamp-1)
	measure_raw_sq		=num/denom
	######################################################
	gdiffs				=circ_diff(circ_eigendata)
	gdiffs2				=gdiffs .* gdiffs
	gdiffsm				=sqrt.(sum(gdiffs2,1))
	num					=sum(gdiffsm)/nsamp

	gdiffs				=diff(lin_eigendata,2)
	gdiffs2				=gdiffs .* gdiffs
	gdiffsm				=sqrt.(sum(gdiffs2,1))
	denom				=sum(gdiffsm)/(nsamp-1)
	measure_eigen		=num/denom

	######################################################
	gdiffs				=circ_diff(circ_eigendata)
	gdiffs2				=gdiffs .* gdiffs
	gdiffsm				=(sum(gdiffs2,1))
	num					=sum(gdiffsm)/nsamp

	gdiffs				=diff(lin_eigendata,2)
	gdiffs2				=gdiffs .* gdiffs
	gdiffsm				=(sum(gdiffs2,1))
	denom				=sum(gdiffsm)/(nsamp-1)
	measure_eigen_sq	=num/denom
    ###########################################################################

	[measure_eigen,measure_eigen_sq,measure_raw,measure_raw_sq]

end

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

function backgroundmetrics_global_eigen(seed_ldata::Array{Float64,2},ESize::Integer,N_best::Integer,N_runs::Integer)

	background_global_metrics=zeros(N_runs,3)
	netsize=ESize

	for count=1:N_runs
		best_error=10E40;

		shuffled_seed_data=Row_Shuffle(seed_ldata)
		EigenShuffledData=GetNEigenGenes(shuffled_seed_data,ESize)
		estimated_phaselist,bestnet,variance_global_metrics=CYCLOPS_Order(ESize,EigenShuffledData,N_best);
		background_global_metrics[count,:]=variance_global_metrics
	end
	background_global_metrics
end

###########################################################################
function multicore_backgroundmetrics_global_eigen(seed_ldata::Array{Float64,2},ESize::Integer,N_best::Integer,N_trials::Integer)

	a1=@spawn backgroundmetrics_global_eigen(seed_ldata,ESize,N_best,Int(N_trials/5));
	a2=@spawn backgroundmetrics_global_eigen(seed_ldata,ESize,N_best,Int(N_trials/5));
	a3=@spawn backgroundmetrics_global_eigen(seed_ldata,ESize,N_best,Int(N_trials/5));
	a4=@spawn backgroundmetrics_global_eigen(seed_ldata,ESize,N_best,Int(N_trials/5));
	a5=@spawn backgroundmetrics_global_eigen(seed_ldata,ESize,N_best,Int(N_trials/5));

	global1=fetch(a1);
	global2=fetch(a2);
	global3=fetch(a3);
	global4=fetch(a4);
	global5=fetch(a5);

	global_back=vcat(global1,global2,global3,global4,global5);
	global_back
end

###########################################################################
function multicore_backgroundstatistics_global_eigen(seed_ldata::Array{Float64,2},ESize::Integer,N_best::Integer,N_trials::Integer,testmetrics)

	backgroundmetrics=multicore_backgroundmetrics_global_eigen(seed_ldata,ESize,N_best,N_trials)

	gm1=backgroundmetrics[:,1]
	gm2=backgroundmetrics[:,2]
	gm3=backgroundmetrics[:,3]

	nless1=length(gm1[gm1 .<= testmetrics[1]]);
	nless2=length(gm2[gm2 .<= testmetrics[2]]);
	nless3=length(gm3[gm3 .<= testmetrics[3]]);

	p1=nless1/N_trials
	p2=nless2/N_trials
	p3=nless3/N_trials

	[p1,p2,p3]
	###Note p3 is the statistical measure to be used########
end

###########################################################################
end

max_stars_repo_name: AdithyaVenkateshMohan/CYCLOPS

max_issues_repo_name: AdithyaVenkateshMohan/CYCLOPS

max_forks_repo_name: AdithyaVenkateshMohan/CYCLOPS


text:
function test_simovici_jaroszewicz_distance()
  v = 6

  u = [1; 1; 1; 1; 1; 1]
  nu = [6; 0; 0; 0; 0; 0]
  ku = 1

  a = [1; 1; 1; 1; 2; 2]
  na = [4; 2; 0; 0; 0; 0]
  ka = 2

  b = [1; 2; 3; 1; 2; 3]
  nb = [2; 2; 2; 0; 0; 0]
  kb = 3

  c = [1; 1; 2; 3; 4; 3]
  nc = [2; 1; 2; 1; 0; 0]
  kc = 4

  d = [1; 2; 3; 4; 5; 1]
  nd = [2; 1; 1; 1; 1; 0]
  kd = 5

  t = [1; 2; 3; 4; 5; 6]
  nt = [1; 1; 1; 1; 1; 1]
  kt = 6

  # β = 2.0
  @test Kpax3.distsj(u, nu, ku, u, nu, ku, v, β=2.0) == 0.0
  @test isapprox(Kpax3.distsj(u, nu, ku, a, na, ka, v, β=2.0), 0.5333333333333333)
  @test isapprox(Kpax3.distsj(u, nu, ku, b, nb, kb, v, β=2.0), 0.7999999999999999)
  @test isapprox(Kpax3.distsj(u, nu, ku, c, nc, kc, v, β=2.0), 0.8666666666666666)
  @test isapprox(Kpax3.distsj(u, nu, ku, d, nd, kd, v, β=2.0), 0.9333333333333333)
  @test isapprox(Kpax3.distsj(u, nu, ku, t, nt, kt, v, β=2.0), 1.0)

  @test Kpax3.distsjlower(nu, ku, v, β=2.0) == 0.0
  @test isapprox(Kpax3.distsjlower(na, ka, v, β=2.0), 0.5333333333333333)
  @test isapprox(Kpax3.distsjlower(nb, kb, v, β=2.0), 0.7999999999999999)
  @test isapprox(Kpax3.distsjlower(nc, kc, v, β=2.0), 0.8666666666666666)
  @test isapprox(Kpax3.distsjlower(nd, kd, v, β=2.0), 0.9333333333333333)
  @test isapprox(Kpax3.distsjlower(nt, kt, v, β=2.0), 1.0)

  @test Kpax3.distsj(a, na, ka, a, na, ka, v, β=2.0) == 0.0
  @test isapprox(Kpax3.distsj(a, na, ka, u, nu, ku, v, β=2.0), 0.5333333333333333)
  @test isapprox(Kpax3.distsj(a, na, ka, b, nb, kb, v, β=2.0), 0.5333333333333333)
  @test isapprox(Kpax3.distsj(a, na, ka, c, nc, kc, v, β=2.0), 0.4666666666666667)
  @test isapprox(Kpax3.distsj(a, na, ka, d, nd, kd, v, β=2.0), 0.5333333333333333)
  @test isapprox(Kpax3.distsj(a, na, ka, t, nt, kt, v, β=2.0), 0.4666666666666667)

  @test Kpax3.distsj(b, nb, kb, b, nb, kb, v, β=2.0) == 0.0
  @test isapprox(Kpax3.distsj(b, nb, kb, u, nu, ku, v, β=2.0), 0.7999999999999999)
  @test isapprox(Kpax3.distsj(b, nb, kb, a, na, ka, v, β=2.0), 0.5333333333333333)
  @test isapprox(Kpax3.distsj(b, nb, kb, c, nc, kc, v, β=2.0), 0.3333333333333333)
  @test isapprox(Kpax3.distsj(b, nb, kb, d, nd, kd, v, β=2.0), 0.2666666666666667)
  @test isapprox(Kpax3.distsj(b, nb, kb, t, nt, kt, v, β=2.0), 0.2)

  @test Kpax3.distsj(c, nc, kc, c, nc, kc, v, β=2.0) == 0.0
  @test isapprox(Kpax3.distsj(c, nc, kc, u, nu, ku, v, β=2.0), 0.8666666666666666)
  @test isapprox(Kpax3.distsj(c, nc, kc, a, na, ka, v, β=2.0), 0.4666666666666667)
  @test isapprox(Kpax3.distsj(c, nc, kc, b, nb, kb, v, β=2.0), 0.3333333333333333)
  @test isapprox(Kpax3.distsj(c, nc, kc, d, nd, kd, v, β=2.0), 0.2)
  @test isapprox(Kpax3.distsj(c, nc, kc, t, nt, kt, v, β=2.0), 0.1333333333333333)

  @test Kpax3.distsj(d, nd, kd, d, nd, kd, v, β=2.0) == 0.0
  @test isapprox(Kpax3.distsj(d, nd, kd, u, nu, ku, v, β=2.0), 0.9333333333333333)
  @test isapprox(Kpax3.distsj(d, nd, kd, a, na, ka, v, β=2.0), 0.5333333333333333)
  @test isapprox(Kpax3.distsj(d, nd, kd, b, nb, kb, v, β=2.0), 0.2666666666666667)
  @test isapprox(Kpax3.distsj(d, nd, kd, c, nc, kc, v, β=2.0), 0.2)
  @test isapprox(Kpax3.distsj(d, nd, kd, t, nt, kt, v, β=2.0), 0.0666666666666667)

  @test Kpax3.distsj(t, nt, kt, t, nt, kt, v, β=2.0) == 0.0
  @test isapprox(Kpax3.distsj(t, nt, kt, u, nu, ku, v, β=2.0), 1.0)
  @test isapprox(Kpax3.distsj(t, nt, kt, a, na, ka, v, β=2.0), 0.4666666666666667)
  @test isapprox(Kpax3.distsj(t, nt, kt, b, nb, kb, v, β=2.0), 0.2)
  @test isapprox(Kpax3.distsj(t, nt, kt, c, nc, kc, v, β=2.0), 0.1333333333333333)
  @test isapprox(Kpax3.distsj(t, nt, kt, d, nd, kd, v, β=2.0), 0.0666666666666667)

  @test Kpax3.distsjupper(nt, kt, v, β=2.0) == 0.0
  @test isapprox(Kpax3.distsjupper(nu, ku, v, β=2.0), 1.0)
  @test isapprox(Kpax3.distsjupper(na, ka, v, β=2.0), 0.4666666666666667)
  @test isapprox(Kpax3.distsjupper(nb, kb, v, β=2.0), 0.2)
  @test isapprox(Kpax3.distsjupper(nc, kc, v, β=2.0), 0.1333333333333333)
  @test isapprox(Kpax3.distsjupper(nd, kd, v, β=2.0), 0.0666666666666667)

  @test Kpax3.distsj(u, nu, ku, t, nt, kt, v, β=2.0) <= (Kpax3.distsj(u, nu, ku, b, nb, kb, v, β=2.0) + Kpax3.distsj(b, nb, kb, t, nt, kt, v, β=2.0))
  @test Kpax3.distsj(a, na, ka, b, nb, kb, v, β=2.0) <= (Kpax3.distsj(a, na, ka, c, nc, kc, v, β=2.0) + Kpax3.distsj(c, nc, kc, b, nb, kb, v, β=2.0))
  @test Kpax3.distsj(b, nb, kb, c, nc, kc, v, β=2.0) <= (Kpax3.distsj(b, nb, kb, d, nd, kd, v, β=2.0) + Kpax3.distsj(d, nd, kd, c, nc, kc, v, β=2.0))
  @test Kpax3.distsj(c, nc, kc, d, nd, kd, v, β=2.0) <= (Kpax3.distsj(c, nc, kc, a, na, ka, v, β=2.0) + Kpax3.distsj(a, na, ka, d, nd, kd, v, β=2.0))
  @test Kpax3.distsj(a, na, ka, b, nb, kb, v, β=2.0) <= (Kpax3.distsj(a, na, ka, u, nu, ku, v, β=2.0) + Kpax3.distsj(u, nu, ku, b, nb, kb, v, β=2.0))
  @test Kpax3.distsj(a, na, ka, b, nb, kb, v, β=2.0) <= (Kpax3.distsj(a, na, ka, t, nt, kt, v, β=2.0) + Kpax3.distsj(t, nt, kt, b, nb, kb, v, β=2.0))
  @test Kpax3.distsj(a, na, ka, d, nd, kd, v, β=2.0) <= (Kpax3.distsj(a, na, ka, u, nu, ku, v, β=2.0) + Kpax3.distsj(u, nu, ku, d, nd, kd, v, β=2.0))
  @test Kpax3.distsj(a, na, ka, d, nd, kd, v, β=2.0) <= (Kpax3.distsj(a, na, ka, t, nt, kt, v, β=2.0) + Kpax3.distsj(t, nt, kt, d, nd, kd, v, β=2.0))

  # β = 3.0
  @test Kpax3.distsj(u, nu, ku, u, nu, ku, v, β=3.0) == 0.0
  @test isapprox(Kpax3.distsj(u, nu, ku, a, na, ka, v, β=3.0), 0.6857142857142857)
  @test isapprox(Kpax3.distsj(u, nu, ku, b, nb, kb, v, β=3.0), 0.9142857142857143)
  @test isapprox(Kpax3.distsj(u, nu, ku, c, nc, kc, v, β=3.0), 0.9428571428571429)
  @test isapprox(Kpax3.distsj(u, nu, ku, d, nd, kd, v, β=3.0), 0.9714285714285714)
  @test isapprox(Kpax3.distsj(u, nu, ku, t, nt, kt, v, β=3.0), 1.0)

  @test Kpax3.distsjlower(nu, ku, v, β=3.0) == 0.0
  @test isapprox(Kpax3.distsjlower(na, ka, v, β=3.0), 0.68571428571428572063)
  @test isapprox(Kpax3.distsjlower(nb, kb, v, β=3.0), 0.91428571428571425717)
  @test isapprox(Kpax3.distsjlower(nc, kc, v, β=3.0), 0.94285714285714294913)
  @test isapprox(Kpax3.distsjlower(nd, kd, v, β=3.0), 0.97142857142857141906)
  @test isapprox(Kpax3.distsjlower(nt, kt, v, β=3.0), 1.0)

  @test Kpax3.distsj(a, na, ka, a, na, ka, v, β=3.0) == 0.0
  @test isapprox(Kpax3.distsj(a, na, ka, u, nu, ku, v, β=3.0), 0.6857142857142857)
  @test isapprox(Kpax3.distsj(a, na, ka, b, nb, kb, v, β=3.0), 0.3428571428571429)
  @test isapprox(Kpax3.distsj(a, na, ka, c, nc, kc, v, β=3.0), 0.3142857142857143)
  @test isapprox(Kpax3.distsj(a, na, ka, d, nd, kd, v, β=3.0), 0.3428571428571429)
  @test isapprox(Kpax3.distsj(a, na, ka, t, nt, kt, v, β=3.0), 0.3142857142857143)

  @test Kpax3.distsj(b, nb, kb, b, nb, kb, v, β=3.0) == 0.0
  @test isapprox(Kpax3.distsj(b, nb, kb, u, nu, ku, v, β=3.0), 0.9142857142857143)
  @test isapprox(Kpax3.distsj(b, nb, kb, a, na, ka, v, β=3.0), 0.3428571428571429)
  @test isapprox(Kpax3.distsj(b, nb, kb, c, nc, kc, v, β=3.0), 0.1428571428571429)
  @test isapprox(Kpax3.distsj(b, nb, kb, d, nd, kd, v, β=3.0), 0.1142857142857143)
  @test isapprox(Kpax3.distsj(b, nb, kb, t, nt, kt, v, β=3.0), 0.0857142857142857)

  @test Kpax3.distsj(c, nc, kc, c, nc, kc, v, β=3.0) == 0.0
  @test isapprox(Kpax3.distsj(c, nc, kc, u, nu, ku, v, β=3.0), 0.9428571428571429)
  @test isapprox(Kpax3.distsj(c, nc, kc, a, na, ka, v, β=3.0), 0.3142857142857143)
  @test isapprox(Kpax3.distsj(c, nc, kc, b, nb, kb, v, β=3.0), 0.1428571428571429)
  @test isapprox(Kpax3.distsj(c, nc, kc, d, nd, kd, v, β=3.0), 0.0857142857142857)
  @test isapprox(Kpax3.distsj(c, nc, kc, t, nt, kt, v, β=3.0), 0.0571428571428571)

  @test Kpax3.distsj(d, nd, kd, d, nd, kd, v, β=3.0) == 0.0
  @test isapprox(Kpax3.distsj(d, nd, kd, u, nu, ku, v, β=3.0), 0.9714285714285714)
  @test isapprox(Kpax3.distsj(d, nd, kd, a, na, ka, v, β=3.0), 0.3428571428571429)
  @test isapprox(Kpax3.distsj(d, nd, kd, b, nb, kb, v, β=3.0), 0.1142857142857143)
  @test isapprox(Kpax3.distsj(d, nd, kd, c, nc, kc, v, β=3.0), 0.0857142857142857)
  @test isapprox(Kpax3.distsj(d, nd, kd, t, nt, kt, v, β=3.0), 0.0285714285714286)

  @test Kpax3.distsj(t, nt, kt, t, nt, kt, v, β=3.0) == 0.0
  @test isapprox(Kpax3.distsj(t, nt, kt, u, nu, ku, v, β=3.0), 1.0)
  @test isapprox(Kpax3.distsj(t, nt, kt, a, na, ka, v, β=3.0), 0.3142857142857143)
  @test isapprox(Kpax3.distsj(t, nt, kt, b, nb, kb, v, β=3.0), 0.0857142857142857)
  @test isapprox(Kpax3.distsj(t, nt, kt, c, nc, kc, v, β=3.0), 0.0571428571428571)
  @test isapprox(Kpax3.distsj(t, nt, kt, d, nd, kd, v, β=3.0), 0.0285714285714286)

  @test Kpax3.distsjupper(nt, kt, v, β=3.0) == 0.0
  @test isapprox(Kpax3.distsjupper(nu, ku, v, β=3.0), 1.0)
  @test isapprox(Kpax3.distsjupper(na, ka, v, β=3.0), 0.3142857142857143)
  @test isapprox(Kpax3.distsjupper(nb, kb, v, β=3.0), 0.0857142857142857)
  @test isapprox(Kpax3.distsjupper(nc, kc, v, β=3.0), 0.0571428571428571)
  @test isapprox(Kpax3.distsjupper(nd, kd, v, β=3.0), 0.0285714285714286)

  @test Kpax3.distsj(u, nu, ku, t, nt, kt, v, β=3.0) <= (Kpax3.distsj(u, nu, ku, b, nb, kb, v, β=3.0) + Kpax3.distsj(b, nb, kb, t, nt, kt, v, β=3.0))
  @test Kpax3.distsj(a, na, ka, b, nb, kb, v, β=3.0) <= (Kpax3.distsj(a, na, ka, c, nc, kc, v, β=3.0) + Kpax3.distsj(c, nc, kc, b, nb, kb, v, β=3.0))
  @test Kpax3.distsj(b, nb, kb, c, nc, kc, v, β=3.0) <= (Kpax3.distsj(b, nb, kb, d, nd, kd, v, β=3.0) + Kpax3.distsj(d, nd, kd, c, nc, kc, v, β=3.0))
  @test Kpax3.distsj(c, nc, kc, d, nd, kd, v, β=3.0) <= (Kpax3.distsj(c, nc, kc, a, na, ka, v, β=3.0) + Kpax3.distsj(a, na, ka, d, nd, kd, v, β=3.0))
  @test Kpax3.distsj(a, na, ka, b, nb, kb, v, β=3.0) <= (Kpax3.distsj(a, na, ka, u, nu, ku, v, β=3.0) + Kpax3.distsj(u, nu, ku, b, nb, kb, v, β=3.0))
  @test Kpax3.distsj(a, na, ka, b, nb, kb, v, β=3.0) <= (Kpax3.distsj(a, na, ka, t, nt, kt, v, β=3.0) + Kpax3.distsj(t, nt, kt, b, nb, kb, v, β=3.0))
  @test Kpax3.distsj(a, na, ka, d, nd, kd, v, β=3.0) <= (Kpax3.distsj(a, na, ka, u, nu, ku, v, β=3.0) + Kpax3.distsj(u, nu, ku, d, nd, kd, v, β=3.0))
  @test Kpax3.distsj(a, na, ka, d, nd, kd, v, β=3.0) <= (Kpax3.distsj(a, na, ka, t, nt, kt, v, β=3.0) + Kpax3.distsj(t, nt, kt, d, nd, kd, v, β=3.0))

  nothing
end

test_simovici_jaroszewicz_distance()


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule eccodes_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("eccodes")
JLLWrappers.@generate_main_file("eccodes", UUID("b04048ba-5ccd-5610-b3f6-85129a548705"))
end  # module eccodes_jll


text:
function FieldOfViewDetector(arg0::PVCoordinatesProvider, arg1::FieldOfView)
    return FieldOfViewDetector((PVCoordinatesProvider, FieldOfView), arg0, arg1)
end

function FieldOfViewDetector(arg0::PVCoordinatesProvider, arg1::jdouble, arg2::VisibilityTrigger, arg3::FieldOfView)
    return FieldOfViewDetector((PVCoordinatesProvider, jdouble, VisibilityTrigger, FieldOfView), arg0, arg1, arg2, arg3)
end

function g(obj::FieldOfViewDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_fov(obj::FieldOfViewDetector)
    return jcall(obj, "getFOV", FieldOfView, ())
end

function get_pv_target(obj::FieldOfViewDetector)
    return jcall(obj, "getPVTarget", PVCoordinatesProvider, ())
end


text:
"""
	function beCount(nbar::Real)

Returns Bose-Einstein distributed counts for average count rate nbar
"""
function beCount(nbar::Real)
	p = 1/(nbar+1)
	f = p*nbar
	s = p
	r = rand()
	count = 0
	while r>s
		count += 1
		p *= f
		s += p
	end
	return count
end

export poissonCount


text:
using ParametricCurves
using StaticArrays
using Test

@testset "ParametricCurves.jl" begin
    p = ParametricCurves.ParametricCurve((t -> cos(t), t -> sin(t)))
    @test ParametricCurves.tangent(p, 0.) ≈ SVector{2,Float64}(0, 1)
    @test ParametricCurves.normal(p, 0.) ≈ SVector{2,Float64}(-1, 0)
    @test ParametricCurves.radius(p, 0) ≈ 1.0
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

struct Footer <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Footer) = (:version, :schema, :dictionaries, :recordBatches, :custom_metadata)

function Base.getproperty(x::Footer, field::Symbol)
    if field === :version
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), MetadataVersion)
        return MetadataVersions.V1
    elseif field === :schema
        o = FlatBuffers.offset(x, 6)
        if o != 0
            y = FlatBuffers.indirect(x, o + FlatBuffers.pos(x))
            return FlatBuffers.init(Schema, FlatBuffers.bytes(x), y)
        end
    elseif field === :dictionaries
        o = FlatBuffers.offset(x, 8)
        if o != 0
            return FlatBuffers.Array{Block}(x, o)
        end
    elseif field === :recordBatches
        o = FlatBuffers.offset(x, 10)
        if o != 0
            return FlatBuffers.Array{Block}(x, o)
        end
    elseif field === :custom_metadata
        o = FlatBuffers.offset(x, 12)
        if o != 0
            return FlatBuffers.Array{KeyValue}(x, o)
        end
    end
    return nothing
end

footerStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 4)
footerAddVersion(b::FlatBuffers.Builder, version::MetadataVersion) = FlatBuffers.prependslot!(b, 0, version, 0)
footerAddSchema(b::FlatBuffers.Builder, schema::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 1, schema, 0)
footerAddDictionaries(b::FlatBuffers.Builder, dictionaries::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 2, dictionaries, 0)
footerStartDictionariesVector(b::FlatBuffers.Builder, numelems) = FlatBuffers.startvector!(b, 24, numelems, 8)
footerAddRecordBatches(b::FlatBuffers.Builder, recordbatches::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 3, recordbatches, 0)
footerStartRecordBatchesVector(b::FlatBuffers.Builder, numelems) = FlatBuffers.startvector!(b, 24, numelems, 8)
footerEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct Block <: FlatBuffers.Struct
    bytes::Vector{UInt8}
    pos::Base.Int
end

FlatBuffers.structsizeof(::Base.Type{Block}) = 24

Base.propertynames(x::Block) = (:offset, :metaDataLength, :bodyLength)

function Base.getproperty(x::Block, field::Symbol)
    if field === :offset
        return FlatBuffers.get(x, FlatBuffers.pos(x), Int64)
    elseif field === :metaDataLength
        return FlatBuffers.get(x, FlatBuffers.pos(x) + 8, Int32)
    elseif field === :bodyLength
        return FlatBuffers.get(x, FlatBuffers.pos(x) + 16, Int64)
    end
    return nothing
end

function createBlock(b::FlatBuffers.Builder, offset::Int64, metadatalength::Int32, bodylength::Int64)
    FlatBuffers.prep!(b, 8, 24)
    prepend!(b, bodylength)
    FlatBuffers.pad!(b, 4)
    prepend!(b, metadatalength)
    prepend!(b, offset)
    return FlatBuffers.offset(b)
end

max_stars_repo_name: jonalm/arrow-julia

max_issues_repo_name: jonalm/arrow-julia

max_forks_repo_name: jonalm/arrow-julia


text:
# ------------------------------------------------------------
# Parallel -> Algorithm is run in parallel
# ------------------------------------------------------------
@define_trait Parallel

@define_traitfn Parallel init_subproblems!(lshaped::AbstractLShapedSolver{F,T,A,M,S}, subsolver::SubSolver) where {F, T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver} = begin
    function init_subproblems!(lshaped::AbstractLShapedSolver{F,T,A,M,S}, subsolver::SubSolver, !Parallel) where {F, T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver}
        # Prepare the subproblems
        load_subproblems!(lshaped, scenarioproblems(lshaped.stochasticprogram), subsolver)
        append!(lshaped.subobjectives, zeros(nbundles(lshaped)))
        return lshaped
    end

    function init_subproblems!(lshaped::AbstractLShapedSolver{F,T,A,M,S}, subsolver::SubSolver, Parallel) where {F, T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver}
        @unpack κ = lshaped.parameters
        # Load initial decision
        put!(lshaped.decisions, 1, lshaped.x)
        # Create subproblems on worker processes
        @sync begin
            for w in workers()
                lshaped.subworkers[w-1] = RemoteChannel(() -> Channel{Vector{SubProblem{F,T,A,S}}}(1), w)
                @async load_worker!(scenarioproblems(lshaped.stochasticprogram), w, lshaped.subworkers[w-1], lshaped.x, subsolver, lshaped.parameters.bundle)
            end
            # Prepare memory
            push!(lshaped.subobjectives, zeros(nbundles(lshaped)))
            push!(lshaped.finished, 0)
            log_val = lshaped.parameters.log
            lshaped.parameters.log = false
            log!(lshaped)
            lshaped.parameters.log = log_val
        end
        # Prepare work channels
        for w in workers()
            lshaped.work[w-1] = RemoteChannel(() -> Channel{Int}(10), w)
            put!(lshaped.work[w-1], 1)
        end
        return lshaped
    end
end

@define_traitfn Parallel iterate!(lshaped::AbstractLShapedSolver) = begin
    function iterate!(lshaped::AbstractLShapedSolver, !Parallel)
        iterate_nominal!(lshaped)
    end

    function iterate!(lshaped::AbstractLShapedSolver, Parallel)
        iterate_parallel!(lshaped)
    end
end

@define_traitfn Parallel nbundles(lshaped::AbstractLShapedSolver) = begin
    function nbundles(lshaped::AbstractLShapedSolver, !Parallel)
        return ceil(Int, lshaped.nscenarios/lshaped.parameters.bundle)
    end

    function nbundles(lshaped::AbstractLShapedSolver, Parallel)
        return nbundles(lshaped, scenarioproblems(lshaped.stochasticprogram))
    end
end

function nbundles(lshaped::AbstractLShapedSolver, sp::StochasticPrograms.ScenarioProblems)
    jobsize = ceil(Int, lshaped.nscenarios/nworkers())
    n = ceil(Int, jobsize/lshaped.parameters.bundle)
    remainder = lshaped.nscenarios-(nworkers()-1)*jobsize
    bundlerem = ceil(Int, remainder/lshaped.parameters.bundle)
    return n*(nworkers()-1)+bundlerem
end

function nbundles(lshaped::AbstractLShapedSolver, sp::StochasticPrograms.DScenarioProblems)
    return sum([ceil(Int, nscen/lshaped.parameters.bundle) for nscen in sp.scenario_distribution])
end

@define_traitfn Parallel init_workers!(lshaped::AbstractLShapedSolver) = begin
    function init_workers!(lshaped::AbstractLShapedSolver, Parallel)
        for w in workers()
            lshaped.active_workers[w-1] = remotecall(work_on_subproblems!,
                                                     w,
                                                     lshaped.subworkers[w-1],
                                                     lshaped.work[w-1],
                                                     lshaped.cutqueue,
                                                     lshaped.decisions,
                                                     lshaped.parameters.bundle)
        end
        return nothing
    end
end

@define_traitfn Parallel close_workers!(lshaped::AbstractLShapedSolver) = begin
    function close_workers!(lshaped::AbstractLShapedSolver, Parallel)
        map((w)->close(w), lshaped.work)
        map(wait, lshaped.active_workers)
        return nothing
    end
end

@define_traitfn Parallel calculate_objective_value(lshaped::AbstractLShapedSolver, x::AbstractVector) = begin
    function calculate_objective_value(lshaped::AbstractLShapedSolver, x::AbstractVector, !Parallel)
        return lshaped.c⋅x + sum([subproblem.π*subproblem(x) for subproblem in lshaped.subproblems])
    end

    function calculate_objective_value(lshaped::AbstractLShapedSolver, x::AbstractVector, Parallel)
        Qs = Vector{Float64}(undef, nworkers())
        @sync begin
            for (w,worker) in enumerate(lshaped.subworkers)
                @async Qs[w] = remotecall_fetch(calculate_subobjective, w+1, worker, x)
            end
        end
        return lshaped.c⋅x + sum(Qs)
    end
end

@define_traitfn Parallel fill_submodels!(lshaped::AbstractLShapedSolver, scenarioproblems::StochasticPrograms.ScenarioProblems) = begin
    function fill_submodels!(lshaped::AbstractLShapedSolver, scenarioproblems::StochasticPrograms.ScenarioProblems, !Parallel)
        for (i, submodel) in enumerate(scenarioproblems.problems)
            lshaped.subproblems[i](decision(lshaped))
            fill_submodel!(submodel, lshaped.subproblems[i])
        end
    end

    function fill_submodels!(lshaped::AbstractLShapedSolver, scenarioproblems::StochasticPrograms.ScenarioProblems, Parallel)
        j = 0
        @sync begin
            for w in workers()
                n = remotecall_fetch((sw)->length(fetch(sw)), w, lshaped.subworkers[w-1])
                for i = 1:n
                    k = i+j
                    @async fill_submodel!(scenarioproblems.problems[k],remotecall_fetch((sw,i,x)->begin
                        sp = fetch(sw)[i]
                        sp(x)
                        get_solution(sp)
                    end,
                    w,
                    lshaped.subworkers[w-1],
                    i,
                    decision(lshaped))...)
                end
                j += n
            end
        end
    end
end

@define_traitfn Parallel fill_submodels!(lshaped::AbstractLShapedSolver, scenarioproblems::StochasticPrograms.DScenarioProblems) = begin
    function fill_submodels!(lshaped::AbstractLShapedSolver, scenarioproblems::StochasticPrograms.DScenarioProblems, !Parallel)
        j = 0
        @sync begin
            for w in workers()
                n = remotecall_fetch((sp)->length(fetch(sp).problems), w, scenarioproblems[w-1])
                for i in 1:n
                    k = i+j
                    lshaped.subproblems[k](decision(lshaped))
                    @async remotecall_fetch((sp,i,x,μ,λ,C) -> fill_submodel!(fetch(sp).problems[i],x,μ,λ,C),
                                            w,
                                            scenarioproblems[w-1],
                                            i,
                                            get_solution(lshaped.subproblems[k])...)
                end
                j += n
            end
        end
    end

    function fill_submodels!(lshaped::AbstractLShapedSolver, scenarioproblems::StochasticPrograms.DScenarioProblems, Parallel)
        @sync begin
            for w in workers()
                @async remotecall_fetch(fill_submodels!,
                                        w,
                                        lshaped.subworkers[w-1],
                                        decision(lshaped),
                                        scenarioproblems[w-1])
            end
        end
    end
end

# Parallel routines #
# ======================================================================== #
mutable struct DecisionChannel{A <: AbstractArray} <: AbstractChannel{A}
    decisions::Dict{Int,A}
    cond_take::Condition
    DecisionChannel(decisions::Dict{Int,A}) where A <: AbstractArray = new{A}(decisions, Condition())
end

function put!(channel::DecisionChannel, t, x)
    channel.decisions[t] = copy(x)
    notify(channel.cond_take)
    return channel
end

function take!(channel::DecisionChannel, t)
    x = fetch(channel, t)
    delete!(channel.decisions, t)
    return x
end

isready(channel::DecisionChannel) = length(channel.decisions) > 1
isready(channel::DecisionChannel, t) = haskey(channel.decisions, t)

function fetch(channel::DecisionChannel, t)
    wait(channel, t)
    return channel.decisions[t]
end

function wait(channel::DecisionChannel, t)
    while !isready(channel, t)
        wait(channel.cond_take)
    end
end

SubWorker{F,T,A,S} = RemoteChannel{Channel{Vector{SubProblem{F,T,A,S}}}}
ScenarioProblems{D,S} = RemoteChannel{Channel{StochasticPrograms.ScenarioProblems{D,S}}}
Work = RemoteChannel{Channel{Int}}
Decisions{A} = RemoteChannel{DecisionChannel{A}}
QCut{T} = Tuple{Int,T,SparseHyperPlane{T}}
CutQueue{T} = RemoteChannel{Channel{QCut{T}}}

function load_subproblems!(lshaped::AbstractLShapedSolver{F,T,A}, scenarioproblems::StochasticPrograms.ScenarioProblems, subsolver::MPB.AbstractMathProgSolver) where {F, T <: Real, A <: AbstractVector}
    id = 1
    for i = 1:lshaped.nscenarios
        m = subproblem(scenarioproblems, i)
        y₀ = convert(A, rand(m.numCols))
        push!(lshaped.subproblems,SubProblem(m,
                                             parentmodel(scenarioproblems),
                                             id,
                                             probability(scenario(scenarioproblems,i)),
                                             copy(lshaped.x),
                                             y₀,
                                             subsolver,
                                             F))
        if i % lshaped.parameters.bundle == 0
            id += 1
        end
    end
    return lshaped
end

function load_subproblems!(lshaped::AbstractLShapedSolver{F,T,A}, scenarioproblems::StochasticPrograms.DScenarioProblems, subsolver::MPB.AbstractMathProgSolver) where {F, T <: Real, A <: AbstractVector}
    id = 1
    for i = 1:lshaped.nscenarios
        m = subproblem(scenarioproblems, i)
        y₀ = convert(A, rand(m.numCols))
        push!(lshaped.subproblems,SubProblem(m,
                                             id,
                                             probability(scenario(scenarioproblems,i)),
                                             copy(lshaped.x),
                                             y₀,
                                             masterterms(scenarioproblems,i),
                                             subsolver,
                                             F))
        if i % lshaped.parameters.bundle == 0
            id += 1
        end
    end
    return lshaped
end

function load_worker!(sp::StochasticPrograms.ScenarioProblems,
                      w::Integer,
                      worker::SubWorker,
                      x::AbstractVector,
                      subsolver::SubSolver,
                      bundlesize::Integer)
    n = StochasticPrograms.nscenarios(sp)
    (nscen, extra) = divrem(n, nworkers())
    prev = [nscen + (extra + 2 - p > 0) for p in 2:(w-1)]
    start = isempty(prev) ? 1 : sum(prev) + 1
    stop = min(start + nscen + (extra + 2 - w > 0) - 1, n)
    prev = [begin
            jobsize = nscen + (extra + 2 - p > 0)
            ceil(Int, jobsize/bundlesize)
            end for p in 2:(w-1)]
    start_id = isempty(prev) ? 1 : sum(prev) + 1
    πs = [probability(sp.scenarios[i]) for i = start:stop]
    return remotecall_fetch(init_subworker!,
                            w,
                            worker,
                            sp.parent,
                            sp.problems[start:stop],
                            πs,
                            x,
                            subsolver,
                            bundlesize,
                            start_id)
end

function load_worker!(sp::StochasticPrograms.DScenarioProblems,
                      w::Integer,
                      worker::SubWorker,
                      x::AbstractVector,
                      subsolver::SubSolver,
                      bundlesize::Integer)
    prev = [ceil(Int,sp.scenario_distribution[p-1]/bundlesize) for p in 2:(w-1)]
    start_id = isempty(prev) ? 1 : sum(prev)+1
    return remotecall_fetch(init_subworker!,
                            w,
                            worker,
                            sp[w-1],
                            x,
                            subsolver,
                            bundlesize,
                            start_id)
end

function init_subworker!(subworker::SubWorker{F,T,A,S},
                         parent::JuMP.Model,
                         submodels::Vector{JuMP.Model},
                         πs::A,
                         x::A,
                         subsolver::SubSolver,
                         bundlesize::Integer,
                         start_id::Integer) where {F, T <: Real, A <: AbstractArray, S <: LQSolver}
    subproblems = Vector{SubProblem{F,T,A,S}}(undef, length(submodels))
    id = start_id
    for (i,submodel) = enumerate(submodels)
        y₀ = convert(A, rand(submodel.numCols))
        subproblems[i] = SubProblem(submodel, parent, id, πs[i], x, y₀, get_solver(subsolver), F)
        if i % bundlesize == 0
            id += 1
        end
    end
    put!(subworker, subproblems)
    return nothing
end

function init_subworker!(subworker::SubWorker{F,T,A,S},
                         scenarioproblems::ScenarioProblems,
                         x::A,
                         subsolver::SubSolver,
                         bundlesize::Integer,
                         start_id::Integer) where {F, T <: Real, A <: AbstractArray, S <: LQSolver}
    sp = fetch(scenarioproblems)
    subproblems = Vector{SubProblem{F,T,A,S}}(undef, StochasticPrograms.nsubproblems(sp))
    id = start_id
    for (i,submodel) = enumerate(sp.problems)
        y₀ = convert(A, rand(sp.problems[i].numCols))
        subproblems[i] = SubProblem(submodel, sp.parent, id, probability(sp.scenarios[i]), x, y₀, get_solver(subsolver), F)
        if i % bundlesize == 0
            id += 1
        end
    end
    put!(subworker, subproblems)
    return nothing
end

function work_on_subproblems!(subworker::SubWorker{F,T,A,S},
                              work::Work,
                              cuts::CutQueue{T},
                              decisions::Decisions{A},
                              bundlesize::Int) where {F, T <: Real, A <: AbstractArray, S <: LQSolver}
    subproblems::Vector{SubProblem{F,T,A,S}} = fetch(subworker)
    if isempty(subproblems)
       # Workers has nothing do to, return.
       return
    end
    @sync while true
        t::Int = try
            wait(work)
            take!(work)
        catch err
            if err isa InvalidStateException
                # Master closed the work channel. Worker finished
                return
            end
        end
        if t == -1
            # Worker finished
            return
        end
        x::A = fetch(decisions,t)
        if bundlesize == 1
            for subproblem in subproblems
                @async begin
                    update_subproblem!(subproblem, x)
                    cut = subproblem()
                    Q::T = cut(x)
                    put!(cuts, (t,Q,cut))
                end
            end
        else
            njobs = ceil(Int, length(subproblems)/bundlesize)
            for i = 1:njobs
                @async begin
                    cut_bundle = CutBundle(T)
                    for subproblem in subproblems[(i-1)*bundlesize+1:min(i*bundlesize,length(subproblems))]
                        update_subproblem!(subproblem, x)
                        cut::SparseHyperPlane{T} = subproblem()
                        _add_cut(cuts, t, cut,x)
                        _add_to_bundle!(cut_bundle, cut, x)
                    end
                    if cut_bundle.q < Inf
                        put!(cuts, (t,cut_bundle.q,aggregate!(cut_bundle)))
                    end
                end
            end
        end
    end
end

function _add_cut(cuts::CutQueue, t::Integer, cut::HyperPlane, x::AbstractArray)
    put!(cuts, (t,cut(x),cut))
end
function _add_cut(cuts::CutQueue, t::Integer, cut::HyperPlane{OptimalityCut}, x::AbstractArray)
    nothing
end
function _add_to_bundle!(bundle::CutBundle, cut::HyperPlane, x::AbstractArray)
    bundle.q += cut(x)
end
function _add_to_bundle!(bundle::CutBundle, cut::HyperPlane{OptimalityCut}, x::AbstractArray)
    push!(bundle.cuts, cut)
    bundle.q += cut(x)
end

function calculate_subobjective(subworker::SubWorker{F,T,A,S}, x::A) where {F, T <: Real, A <: AbstractArray, S <: LQSolver}
    subproblems::Vector{SubProblem{F,T,A,S}} = fetch(subworker)
    if length(subproblems) > 0
        return sum([subproblem.π*subproblem(x) for subproblem in subproblems])
    else
        return zero(T)
    end
end

function fill_submodels!(subworker::SubWorker{F,T,A,S},
                         x::A,
                         scenarioproblems::ScenarioProblems) where {F, T <: Real, A <: AbstractArray, S <: LQSolver}
    sp = fetch(scenarioproblems)
    subproblems::Vector{SubProblem{F,T,A,S}} = fetch(subworker)
    for (i, submodel) in enumerate(sp.problems)
        subproblems[i](x)
        fill_submodel!(submodel, subproblems[i])
    end
end

function fill_submodel!(submodel::JuMP.Model, x::AbstractVector, μ::AbstractVector, λ::AbstractVector, C::Real)
    submodel.colVal = x
    submodel.redCosts = μ
    submodel.linconstrDuals = λ
    submodel.objVal = C
    submodel.objVal *= submodel.objSense == :Min ? 1 : -1
end

function fill_submodel!(submodel::JuMP.Model, subproblem::SubProblem)
    fill_submodel!(submodel, get_solution(subproblem)...)
end

function iterate_parallel!(lshaped::AbstractLShapedSolver{F,T}) where {F, T <: Real}
    wait(lshaped.cutqueue)
    while isready(lshaped.cutqueue)
        # Add new cuts from subworkers
        t::Int,Q::T,cut::SparseHyperPlane{T} = take!(lshaped.cutqueue)
        if Q == Inf && !F
            @warn "Subproblem $(cut.id) is infeasible, aborting procedure."
            return :Infeasible
        end
        if !bounded(cut)
            map((w,aw)->!isready(aw) && put!(w,-1), lshaped.work, lshaped.active_workers)
            @warn "Subproblem $(cut.id) is unbounded, aborting procedure."
            return :Unbounded
        end
        add_cut!(lshaped, cut, lshaped.subobjectives[t], Q)
        update_objective!(lshaped, cut)
        lshaped.finished[t] += 1
        if lshaped.finished[t] == nbundles(lshaped)
            lshaped.solverdata.timestamp = t
            lshaped.x[:] = fetch(lshaped.decisions, t)
            lshaped.Q_history[t] = current_objective_value(lshaped, lshaped.subobjectives[t])
            lshaped.solverdata.Q = lshaped.Q_history[t]
            lshaped.solverdata.θ = t > 1 ? lshaped.θ_history[t-1] : -1e10
            take_step!(lshaped)
            lshaped.solverdata.θ = lshaped.θ_history[t]
            # Check if optimal
            if check_optimality(lshaped)
                # Optimal, tell workers to stop
                map((w,aw)->!isready(aw) && put!(w,t), lshaped.work, lshaped.active_workers)
                map((w,aw)->!isready(aw) && put!(w,-1), lshaped.work, lshaped.active_workers)
                # Final log
                log!(lshaped, lshaped.solverdata.iterations)
                return :Optimal
            end
        end
    end
    # Resolve master
    t = lshaped.solverdata.iterations
    if lshaped.finished[t] >= lshaped.parameters.κ*nbundles(lshaped) && length(lshaped.cuts) >= nbundles(lshaped)
        try
            solve_problem!(lshaped, lshaped.mastersolver)
        catch
            # Master problem could not be solved for some reason.
            @unpack Q,θ = lshaped.solverdata
            gap = abs(θ-Q)/(abs(Q)+1e-10)
            @warn "Master problem could not be solved, solver returned status $(status(lshaped.mastersolver)). The following relative tolerance was reached: $(@sprintf("%.1e",gap)). Aborting procedure."
            map((w,aw)->!isready(aw) && put!(w,-1), lshaped.work, lshaped.active_workers)
            return :StoppedPrematurely
        end
        if status(lshaped.mastersolver) == :Infeasible
            @warn "Master is infeasible. Aborting procedure."
            map((w,aw)->!isready(aw) && put!(w,-1), lshaped.work, lshaped.active_workers)
            return :Infeasible
        end
        # Update master solution
        update_solution!(lshaped)
        θ = calculate_estimate(lshaped)
        if t > 1 && abs(θ-lshaped.θ_history[t-1]) <= 10*lshaped.parameters.τ*abs(1e-10+θ) && lshaped.finished[t] != nbundles(lshaped)
            # Not enough new information in master. Repeat iterate
            return :Valid
        end
        lshaped.solverdata.θ = θ
        lshaped.θ_history[t] = θ
        # Project (if applicable)
        project!(lshaped)
        # If all work is finished at this timestamp, check optimality
        if lshaped.finished[t] == nbundles(lshaped)
            # Check if optimal
            if check_optimality(lshaped)
                # Optimal, tell workers to stop
                map((w,aw)->!isready(aw) && put!(w,t), lshaped.work, lshaped.active_workers)
                map((w,aw)->!isready(aw) && put!(w,-1), lshaped.work, lshaped.active_workers)
                # Final log
                log!(lshaped, t)
                return :Optimal
            end
        end
        # Log progress at current timestamp
        log_regularization!(lshaped, t)
        # Send new decision vector to workers
        put!(lshaped.decisions, t+1, lshaped.x)
        map((w,aw)->!isready(aw) && put!(w,t+1), lshaped.work, lshaped.active_workers)
        # Prepare memory for next iteration
        push!(lshaped.subobjectives, zeros(nbundles(lshaped)))
        push!(lshaped.finished, 0)
        # Log progress
        log!(lshaped)
        lshaped.θ_history[t+1] = -Inf
    end
    # Just return a valid status for this iteration
    return :Valid
end


text:
#==============================================================================================
Some conventions:
    All instructions are based on methods where arguments are passed from an external source
    (i.e. they do not come from the CPU's memory).

    Note that many functions are explicitly inlined, but probably none of them have to be.

    **These calls do not affect the program counter.**

    TODO warn about the stack doing weird things
==============================================================================================#

# these can be passed to instructions to specify addressing modes when appropriate
abstract type AddressingMode end
abstract type DirectMode <: AddressingMode end
abstract type IndirectMode <: AddressingMode end
struct Direct <: DirectMode end
struct DirectX <: DirectMode end
struct DirectY <: DirectMode end
struct IndirectX <: IndirectMode end
struct IndirectY <: IndirectMode end

# some instructions only have these modes
const DirectXModes = Union{Direct,DirectX}
const DirectYModes = Union{Direct,DirectY}

# utility functions used frequently in instructions
@inline checkNflag!(c::CPU, val::UInt8) = (c.flags.N = val ≥ 0x80)
@inline checkZflag!(c::CPU, val::UInt8) = (c.flags.Z = val == 0x00)
@inline checkCflag!(c::CPU, val1::UInt8, val2::UInt8) = (c.flags.C = overflow(val1, val2))

# pointers that occur in different addressing modes
pointer(::Type{Direct}, ptr::Π, c::CPU, m::Memory) = ptr
pointer(::Type{DirectX}, ptr::Π, c::CPU, m::Memory) = ptr + c.X
pointer(::Type{DirectY}, ptr::Π, c::CPU, m::Memory) = ptr + c.Y
pointer(::Type{IndirectX}, ptr::Π, c::CPU, m::Memory) = Π(deref(UInt16, ptr + c.X, m))
pointer(::Type{IndirectY}, ptr::Π, c::CPU, m::Memory) = Π(deref(UInt16, ptr, m) + c.Y)

# these dereference in a way appropriate for Indirect
function deref(::Type{T}, ptr::Π, c::CPU, m::Memory) where {T<:AddressingMode}
    deref(pointer(T, ptr, c, m), m)
end

function store!(::Type{T}, ptr::Π, c::CPU, m::Memory, val::UInt8) where {T<:AddressingMode}
    store!(m, pointer(T, ptr, c, m), val)
end

stackpush!(c::CPU, m::Memory, val::UInt8) = (store!(m, stackpointer(c), val); c.SP -= 1)
stackpush!(c::CPU, m::Memory, val::UInt16) = (storeback!(m, stackpointer(c), val); c.SP -= 2)

stackpull!(c::CPU, m::Memory) = (ξ = deref(stackpointer(c), val); c.SP += 1; ξ)
stackpull!(::Type{UInt16}, c::CPU, m::Memory) = (stackpull!(c, m) + 0x0100 * stackpull!(c, m))

# this is for the weird arithmetic of the branching instructions (\boxplus)
(⊞)(x::Unsigned, y::Unsigned) = unsigned(signed(x) + signed(y))

#==============================================================================================
    <load, store instructions>
==============================================================================================#
#-----------------------LDA, LDX, LDY----------------------------------------------
# template for ld instructions in immediate mode
function ld!(c::CPU, reg::Symbol, val::UInt8)
    v = store!(c, reg, val)
    checkNflag!(c, v)
    checkZflag!(c, v)
    v
end

function ld!(c::CPU, m::Memory, reg::Symbol, ::Type{T}, ptr::Π) where {T<:AddressingMode}
    ld!(c, reg, deref(T, ptr, c, m))
end
ld!(c::CPU, m::Memory, reg::Symbol, ptr::Π) = ld!(c, m, reg, Direct, ptr)

# TODO  TESTING!!!! REVERT!!!!
ldaa!(c::CPU, v::AbstractVector{UInt8}) = (c.A = v[1])
lda!(c::CPU, val::UInt8) = ld!(c, :A, val)
ldx!(c::CPU, val::UInt8) = ld!(c, :X, val)
ldy!(c::CPU, val::UInt8) = ld!(c, :Y, val)

lda!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:AddressingMode} = ld!(c, m, :A, T, ptr)
ldx!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:AddressingMode} = ld!(c, m, :X, T, ptr)
ldy!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:AddressingMode} = ld!(c, m, :Y, T, ptr)

lda!(c::CPU, m::Memory, ptr::Π) = ld!(c, m, :A, Direct, ptr)
ldx!(c::CPU, m::Memory, ptr::Π) = ld!(c, m, :X, Direct, ptr)
ldy!(c::CPU, m::Memory, ptr::Π) = ld!(c, m, :Y, Direct, ptr)

#---------------------STA, STX, STY------------------------------------------------
# template for st in zero page or absolute mode
function st!(c::CPU, m::Memory, reg::Symbol, ::Type{T}, ptr::Π) where T<:AddressingMode
    store!(m, pointer(T, ptr, c, m), fetch(c, reg))
end
st!(c::CPU, m::Memory, reg::Symbol, ptr::Π) = st!(c, m, reg, Direct, ptr)

sta!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:AddressingMode} = st!(c, m, :A, T, ptr)
stx!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:AddressingMode} = st!(c, m, :X, T, ptr)
sty!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:AddressingMode} = st!(c, m, :Y, T, ptr)

sta!(c::CPU, m::Memory, ptr::Π) = st!(c, m, :A, Direct, ptr)
stx!(c::CPU, m::Memory, ptr::Π) = st!(c, m, :X, Direct, ptr)
sty!(c::CPU, m::Memory, ptr::Π) = st!(c, m, :Y, Direct, ptr)

#==============================================================================================
    </load, store instructions>
==============================================================================================#

#==============================================================================================
    <register transfers>
==============================================================================================#
function t!(c::CPU, reg1::Symbol, reg2::Symbol)
    v = store!(c, reg2, fetch(c, reg1))
    checkNflag!(c, v)
    checkZflag!(c, v)
end

tax!(c::CPU) = t!(c, :A, :X)
tay!(c::CPU) = t!(c, :A, :Y)
txa!(c::CPU) = t!(c, :X, :A)
tya!(c::CPU) = t!(c, :Y, :A)

#==============================================================================================
    </register transfers>
==============================================================================================#

# TODO check effect of other instructions on stack!!!
#==============================================================================================
    <stack operations>
==============================================================================================#
function tsx!(c::CPU)
    c.X = c.SP
    checkNflag!(c, c.X)
    checkZflag!(c, c.X)
    c.X
end

txs!(c::CPU) = (c.SP = c.X)

pha!(c::CPU, m::Memory) = stackpush!(c, m, c.A)

php!(c::CPU, m::Memory) = stackpush!(c, m, c.P)

function pla!(c::CPU, m::Memory)  # pull stack onto A
    c.A = stackpull!(c, m)
    checkNflag!(c, c.A)
    checkZflag!(c, c.A)
    c.A
end

# in example this doesn't set the blank or B flags for some reason... wtf?
plp!(c::CPU, m::Memory) = (c.P = stackpull!(c, m))

#==============================================================================================
    </stack operations>
==============================================================================================#

#==============================================================================================
    <logical operations>
==============================================================================================#
# template for logical operations
function logical!(logical_op::Function, c::CPU, val::UInt8)
    c.A = logical_op(c.A, val)
    checkNflag!(c, c.A)
    checkZflag!(c, c.A)
    c.A
end

and!(c::CPU, val::UInt8) = logical!(&, c, val)
function and!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:AddressingMode}
    and!(c, deref(T, ptr, c, m))
end
and!(c::CPU, m::Memory, ptr::Π) = and!(c, m, Direct, ptr)

eor!(c::CPU, val::UInt8) = logical!(⊻, c, val)
function eor!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:AddressingMode}
    eor!(c, deref(T, ptr, c, m))
end
eor!(c::CPU, m::Memory, ptr::Π) = eor!(c, m, Direct, ptr)

ora!(c::CPU, val::UInt8) = logical!(|, c, val)
function ora!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:AddressingMode}
    ora!(c, deref(T, ptr, c, m))
end
ora!(c::CPU, m::Memory, ptr::Π) = ora!(c, m, Direct, ptr)

function bit!(c::CPU, m::Memory, ptr::Π)
    v = deref(ptr, m)
    c.A & v == 0 ? status!(c, :Z, true) : status!(c, :Z, false)
    status!(c, :V, Bool((0x40 & v) >> 6))
    status!(c, :N, Bool((0x80 & v) >> 7))
    v
end
bit!(c::CPU, m::Memory, ::Type{Direct}, ptr::Π) = bit!(c, m, ptr)

#==============================================================================================
    </logical operations>
==============================================================================================#

#==============================================================================================
    <arithmetic>
==============================================================================================#
function adc!(c::CPU, val::UInt8)
    val += UInt8(status(c, :C))
    checkCflag!(c, c.A, val)
    c.A += val
    checkNflag!(c, c.A)
    checkZflag!(c, c.A)
    c.A
end

function adc!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:AddressingMode}
    adc!(c, deref(pointer(T, ptr, c, m), m))
end
adc!(c::CPU, m::Memory, ptr::Π) = adc!(c, m, Direct, ptr)

function sbc!(c::CPU, val::UInt8)
    val -= ~UInt8(status(c, :C))
    checkCflag!(c, c.A, val)
    c.A -= val
    checkNflag!(c, c.A)
    checkZflag!(c, c.A)
    c.A
end

function sbc!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:AddressingMode}
    sbc!(c, deref(pointer(T, ptr, c, m), m))
end
sbc!(c::CPU, m::Memory, ptr::Π) = sbc!(c, m, Direct, ptr)

function compare!(c::CPU, reg::Symbol, val::UInt8)
    r = fetch(c, reg)
    r ≥ val ? status!(c, :C, true) : status!(c, :C, false)
    ξ = r - val
    checkNflag!(c, ξ)
    checkZflag!(c, ξ)
    r
end
function compare!(c::CPU, m::Memory, reg::Symbol, ::Type{T}, ptr::Π) where {T<:AddressingMode}
    compare!(c, reg, deref(pointer(T, ptr, c, m), m))
end
compare!(c::CPU, m::Memory, reg::Symbol, ptr::Π) = compare!(c, m, reg, Direct, ptr)

cmp!(c::CPU, val::UInt8) = compare!(c, :A, val)
function cmp!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:AddressingMode}
    compare!(c, m, :A, T, ptr)
end
cmp!(c::CPU, m::Memory, ptr::Π) = compare!(c, m, :A, ptr)

# this instruction only has direct mode
cpx!(c::CPU, val::UInt8) = compare!(c, :X, val)
cpx!(c::CPU, m::Memory, ::Type{Direct}, ptr::Π) = compare!(c, m, :X, T, ptr)
cpx!(c::CPU, m::Memory, ptr::Π) = compare!(c, m, :X, ptr)

# this instruction only has direct mode
cpy!(c::CPU, val::UInt8) = compare!(c, :Y, val)
cpy!(c::CPU, m::Memory, ::Type{Direct}, ptr::Π) = compare!(c, m, :Y, T, ptr)
cpy!(c::CPU, m::Memory, ptr::Π) = compare!(c, m, :Y, ptr)
#==============================================================================================
    </arithmetic>
==============================================================================================#

#==============================================================================================
    <increments and decrements>
==============================================================================================#
# note that the result of this needs to be stored
function increment!(c::CPU, val::UInt8)
    ξ = val + 1
    checkZflag!(c, ξ)
    checkNflag!(c, ξ)
    ξ
end

function inc!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:DirectXModes}
    C.x = increment!(c, deref(pointer(T, ptr, c, m), m))
end
inc!(c::CPU, m::Memory, ptr::Π) = inc!(c, m, Direct, ptr)

inx!(c::CPU) = (c.X = increment!(c, c.X))

iny!(c::CPU) = (c.Y = increment!(c, c.Y))

function decrement!(c::CPU, val::UInt8)
    ξ = val - 1
    checkZflag!(c, ξ)
    checkNflag!(c, ξ)
    ξ
end

function dec!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:DirectXModes}
    C.x = decrement!(c, deref(pointer(T, ptr, c, m), m))
end
dec!(c::CPU, m::Memory, ptr::Π) = dec!(c, m, Direct, ptr)

dex!(c::CPU, m::Memory, ptr::Π) = (c.X = decrement!(c, c.X))

dey!(c::CPU, m::Memory, ptr::Π) = (c.Y = decrement!(c, c.Y))

#==============================================================================================
    </increments and decrements>
==============================================================================================#

#==============================================================================================
    <shifts>
==============================================================================================#
# checkZflag! must be called from outside!
# WTF? why is this inconsistent with the other ops???
function arithmetic_shiftleft!(c::CPU, val::UInt8)
    val & 0x80 > 0x00 ? status!(c, :C, true) : status!(c, :C, false)
    ξ = val << 1
    checkNflag!(c, ξ)
    ξ
end

asl!(c::CPU) = (c.A = arithmetic_shiftleft!(c, c.A); checkZflag!(c, c.A); c.A)
function asl!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:DirectXModes}
    ξ = arithmetic_shiftleft!(c, deref(pointer(T, ptr, c, m), m))
    store!(m, ptr, ξ)
    # we neglect checking A == 0 since it didn't change
    ξ
end
asl!(c::CPU, m::Memory, ptr::Π) = asl!(c, m, Direct, ptr)

function logical_shiftright!(c::CPU, val::UInt8)
    val & 0x01 > 0x00 ? status!(c, :C, true) : status!(c, :C, false)
    ξ = val >> 1
    checkZflag!(c, ξ)
    checkNflag!(c, ξ)
    ξ
end

lsr!(c::CPU) = (c.A = logical_shiftright!(c, c.A))
function lsr!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:DirectXModes}
    ξ = logical_shiftright!(c, deref(pointer(T, ptr, c, m), m))
    store!(m, ptr, ξ)
    ξ
end
lsr!(c::CPU, m::Memory, ptr::Π) = lsr!(c, m, Direct, ptr)

# again, checkZflag! must be called from outside!
function rotateleft!(c::CPU, val::UInt8)
    val & 0x80 > 0x00 ? status!(c, :C, true) : status!(c, :C, false)
    ξ = (val << 1) + (0x01 * convert(UInt8, status(c, :C)))
    checkNflag!(c, ξ)
    ξ
end

rol!(c::CPU) = (c.A = rotateleft!(c, c.A); checkZflag!(c, c.A); c.A)
function rol!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:DirectXModes}
    ξ = rotateleft!(c, deref(pointer(T, ptr, c, m), m))
    store!(m, ptr, ξ)
    # we neglect checking A == 0 since it didn't change
    ξ
end
rol!(c::CPU, m::Memory, ptr::Π) = rol!(c, m, Direct, ptr)

# again, checkZflag! must be called from outside!
function rotateright!(c::CPU, val::UInt8)
    val & 0x01 > 0x00 ? status!(c, :C, true) : status!(c, :C, false)
    ξ = (val >> 1) + (0x80 * convert(UInt8, status(c, :C)))
    checkNflag!(c, ξ)
    ξ
end

ror!(c::CPU) = (c.A = rotateright!(c, c.A); checkZflag!(c.A); c.A)
function ror!(c::CPU, m::Memory, ::Type{T}, ptr::Π) where {T<:DirectXModes}
    ξ = rotateright!(c, deref(pointer(T, ptr, c, m), m))
    store!(m, ptr, ξ)
    # we neglect checking A == 0 since it didn't change
    ξ
end
ror!(c::CPU, m::Memory, ptr::Π) = ror!(c, m, Direct, ptr)
#==============================================================================================
    </shifts>
==============================================================================================#

#==============================================================================================
    <jumps and calls>

    TODO be __very__ careful about how you implement the byte counting of these!!!
==============================================================================================#
# note that this takes UInt16 arguments
jmp!(c::CPU, val::UInt16) = (c.PC = val)

# note, the assembly syntax for this function is really fucking weird
# confusingly, the references call this indirect mode
jmp!(c::CPU, m::Memory, ::Type{Direct}, ptr::Π{UInt16}) = jmp!(c, deref(UInt16, ptr, m))
jmp!(c::CPU, m::Memory, ptr::Π{UInt16}) = jmp!(c, m, Direct, ptr)

function jsr!(c::CPU, m::Memory, val::UInt16)
    stackpush!(c, m, c.PC + 0x0002)
    c.PC = val
end

rts!(c::CPU, m::Memory) = (c.PC = stackpull!(UInt16, c, m))
#==============================================================================================
    </jumps and calls>
==============================================================================================#

#==============================================================================================
    <branches>

    # TODO again be __very__ careful about how you implement byte counting for these!!!
==============================================================================================#
bcc!(c::CPU, m::Memory, val::UInt8) = (!status(c, :C) && (c.PC = c.PC ⊞ val); val)

bcs!(c::CPU, m::Memory, val::UInt8) = (status(c, :C) && (c.PC = c.PC ⊞ val); val)

beq!(c::CPU, m::Memory, val::UInt8) = (status(c, :Z) && (c.PC = c.PC ⊞ val); val)

bmi!(c::CPU, m::Memory, val::UInt8) = (status(c, :N) && (c.PC = c.PC ⊞ val); val)

bne!(c::CPU, m::Memory, val::UInt8) = (!status(c, :Z) && (c.PC = c.PC ⊞ val); val)

bpl!(c::CPU, m::Memory, val::UInt8) = (!status(c, :N) && (c.PC = c.PC ⊞ val); val)

bvc!(c::CPU, m::Memory, val::UInt8) = (!status(c, :V) && (c.PC = c.PC ⊞ val); val)

bvs!(c::CPU, m::Memory, val::UInt8) = (status(c, :V) && (c.PC = c.PC ⊞ val); val)
#==============================================================================================
    </branches>
==============================================================================================#

#==============================================================================================
    <status flag changes>
==============================================================================================#
clc!(c::CPU) = status!(c, :C, false)

cld!(c::CPU) = status!(c, :D, false)

cli!(c::CPU) = status!(c, :I, false)

clv!(c::CPU) = status!(c, :V, false)

sec!(c::CPU) = status!(c, :C, true)

sed!(c::CPU) = status!(c, :D, true)

sei!(c::CPU) = status!(c, :I, true)

#==============================================================================================
    <status flag changes>
==============================================================================================#

#==============================================================================================
    <system functions>
==============================================================================================#
brk!(c::CPU) = status!(c, :B, true)

nop!(c::CPU) = ()

# TODO this is likely fucked up
# __NOTE!!__ that PC is also incremented the normal way afterwards from reading the instruction
function rti!(c::CPU, m::Memory)
    ptr = stackpointer(c) + 0x01  # are you sure about this?
    c.P = deref(ptr, m)
    ptr = ptr + 0x01
    pc_small = deref(ptr, m)
    ptr = ptr + 0x01
    pc_big = deref(ptr, m)
    c.PC = bincat(pc_small, pc_big)
    c.SP += 0x03  # not completely confident in this either
end

#==============================================================================================
    </system functions>
==============================================================================================#


text:
# Uses the same optimiser as the base instance. Only works if that solver
# understands lazy constraints and SOCP.
struct OSSBExactNaive <: OSSBOptimisationAlgorithm
  solver
end

function optimise_graves_lai(instance::CombinatorialInstance{T},
                       algo::OSSBExactNaive,
                       θ::Dict{T, Float64}; with_trace::Bool=true) where T
  # No need for an LP formulation! Just need an optimisation solver, given in the `OSSBExactNaive` object.
  t0 = time_ns()

  # TODO: make it an external function that can be overridden by combinatorial instances.
  # TODO: probably wrong when order has some importance, like in paths.
  # Generate a list of solutions. Probably highly inefficient, but it should work for any kind of combinatorial
  # instance, as long as is_feasible and is_partially_acceptable are correctly implemented.
  # Technique: for each item, consider that it may be a part of a solution or not. Hence, build all combinations,
  # stop one branch of computations as soon as a set of arms can no longer lead to a solution (something is already
  # violating a constraint).
  # Initially: either take the first item or not. Maybe the first item cannot lead to a feasible solution.
  potential_solutions = Vector{T}[T[], T[first(all_arm_indices(instance))]]
  filter!(p -> is_partially_acceptable(instance, p), potential_solutions)

  solutions = Vector{T}[]
  for i in all_arm_indices(instance)[2:end]
    # Rebuild the set of potential solutions in this loop: for each existing element, try adding i or keeping the same thing.
    partial_solutions = copy(potential_solutions)
    deleteat!(potential_solutions, 1:length(potential_solutions)) # Empty the vector, so that it can be filled again.

    for p in partial_solutions
      # If this is already a feasible solution, put it aside.
      if length(p) > 0 && is_feasible(instance, p) && !(p in solutions)
        push!(solutions, p)
      end

      # The current potential solution is already known to be... a potential solution.
      push!(potential_solutions, p)

      # The other part of the recursion: try adding the new item i.
      p2 = [p..., i]
      if is_feasible(instance, p2)
        push!(solutions, p2)
      end
      if is_partially_acceptable(instance, p2)
        push!(potential_solutions, p2)
      end
    end
  end

  potential_solutions = nothing # Free some memory.

  # Precompute useful results.
  d = dimension(instance)
  x_max = solve_linear(copy(instance), θ)
  r_max = sum(θ[arm] for arm in x_max)
  Δx(x) = r_max - sum(θ[i] for i in x)
  I = _ossb_confusing_arms(instance, θ)
  # I = all_arm_indices(instance)

  # Complete formulation:
  #   min   ∑ η_x Δ_x
  #   s.t.  η_x ≥ 0                 ∀ x ∈ X
  #         t_i = ∑ η_x x_i         ∀ i ∈ {1, 2… d}
  #         ∑ x_i / t_i ≤ Δ_x^2     ∀ x ∈ X

  # Create the optimisation program.
  m = Model(algo.solver)

  @variable(m, η[x in solutions] >= 0)
  @variable(m, t[i in I] >= 0)
  @variable(m, inv_t[i in I] >= 0)

  @constraint(m, c_t[i in I],       t[i] == sum(η[x] for x in solutions if i in x)) # This line
  # supposes that all arms are used in at least one solution, otherwise the line will fail
  # (reduction on an empty collection). This is a safe assumption.
  @constraint(m, c_inv_t[i in I],   [t[i], inv_t[i], sqrt(2)] in RotatedSecondOrderCone())

  # Add the gap constraints and build the objective. Ignore the optimum solutions.
  obj = AffExpr()
  for x in solutions
    Δx² = (Δx(x))^2
    if Δx² >= 1.0e-5 # TODO: \Delta_min?
      # Always add this solution to the objective function.
      add_to_expression!(obj, Δx(x), η[x]) # Performance-wise, better than repeated +=.

      # Generate a constraint only if it is not vacuous (i.e. not 0 <= Δx²).
      if length([i for i in x if i in I]) > 0
        @constraint(m, sum(inv_t[i] for i in x if i in I) <= Δx²)
      end
    end
  end
  @objective(m, Min, obj)

  # Imagine that the previous operations can be performed without exhausting memory...
  set_silent(m)
  optimize!(m)

  if termination_status(m) != MOI.OPTIMAL
    @show termination_status(m)
    println(m)
    error(termination_status(m))
  end

  η_val = [value(η[x]) for x in solutions]
  res = GravesLaiResults(objective_value(m), solutions, η_val)
  if with_trace
    return res, GravesLaiDetails((time_ns() - t0) / 1_000_000, objective_value(m))
  else
    return res
  end
end