text:
#=

Handling of discrete states and parameters for ESSE

Ignacio Quintero Mächler

t(-_-t)

October 30 2017

=#

"""
    struct Sgh
      g::Set{Int64}
      h::Int64
    end

Composite type for Geographical State with `g` areas & `h` hidden states
"""
struct Sgh
  g::Set{Int64}
  h::Int64
end

"""
    isequal(x::Sgh, y::Sgh)

Compares equality between two `Sgh` types.
"""
isSghequal(x::Sgh, y::Sgh) = x.g == y.g && x.h == y.h

"""
    create_states(k::Int64, h::Int64)

Create EGeoHiSSE states
"""
function create_states(k::Int64, h::Int64)

  # create individual areas subsets
  gs = Array{Set{Int64}, 1}(undef, k)
  for i = Base.OneTo(k)
    gs[i] = Set([i])
  end
  gS = sets(gs)
  popfirst!(gS)            # remove empty
  sort!(gS, by = length)   # arrange by range size

  # add hidden states and create Sgh objects
  S = Array{Sgh, 1}()
  for i = 0:(h-1), j = gS
    push!(S, Sgh(j, i))
  end

  return S
end

"""
  sets(x::Array{String,1})

Create all sets from each string (letters) in x.
"""
function sets(x::Array{String,1})
  ss = [""]
  for elem = x, j = eachindex(ss)
    push!(ss, ss[j]*elem)
  end

  return ss
end

"""
  sets(x::Array{SubString{String},1})

Create all subsets from each string (letters) in x.
"""
function sets(x::Array{SubString{String},1})
  ss = [""]
  for elem = x, j = eachindex(ss)
    push!(ss, ss[j]*elem)
  end

  return ss
end

"""
  sets(x::Vector{Set{Int64}})

Create all sets from a vector of sets of length `1`.
"""
function sets(x::Vector{Set{Int64}})

  ss = Vector{Set{Int64}}()
  push!(ss, Set{Int64}())
  for i = x, j = eachindex(ss)
    push!(ss, union(ss[j], i))
  end

  return ss
end

"""
   vicsubsets(x::String)

Create all vicariance subsets from a string of areas.
"""
function vicsubsets(x::String)
  lx  = lastindex(x)
  sx  = split(x, "")
  ss  = sets(sx)
  sort!(ss, by = length)
  pop!(ss)
  popfirst!(ss)
  ls = lastindex(ss)

  vss = Array{NTuple{2, String},1}()
  for i = Base.OneTo(ls)
    push!(vss, (ss[i], ss[ls-i+1]))
  end

  return vss
end

"""
    vicsubsets(x::Set{Int64})

Create all vicariance subsets from a Set of areas.
"""
function vicsubsets(x::Set{Int64})
  lx  = length(x)

  ss = Vector{Set{Int64}}()
  push!(ss, Set{Int64}())
  for i = x, j = eachindex(ss)
    push!(ss, union(ss[j], i))
  end
  sort!(ss, by = length)
  pop!(ss)
  popfirst!(ss)

  ls = lastindex(ss)
  vss = Array{NTuple{2, Set{Int64}},1}()
  for i = Base.OneTo(ls)
    push!(vss, (ss[i], ss[ls-i+1]))
  end

  return vss
end

"""
    build_par_names(k::Int64, h::Int64, ny::Int64, model::Array{Int64,1})

Build dictionary for parameter names and indexes for EGeoHiSSE for
`k` areas, `h` hidden states and `ny` covariates.
"""
function build_par_names(k    ::Int64,
                         h    ::Int64,
                         ny   ::Int64,
                         model::NTuple{3,Bool})

  # generate individual area names
  ia = Array{Char,1}(undef,0)
  for i = 0:(k-1)
    push!(ia, ('A' + i)::Char)
  end

  ## build parameters name
  par_nams = Array{String,1}(undef,0)

  # add λ names but only one between-regions speciation rate
  for j = 0:(h-1)
    for a::Char = ia
      push!(par_nams, "lambda_"*string(a)*"_"*string(j))
    end
    lastindex(ia) > 1 && push!(par_nams, "lambda_W_"*string(j))
  end

  # add μ names for endemics
  for j = 0:(h-1), i = ia
    push!(par_nams, "mu_"*string(i)*"_"*string(j))
  end

  # add area gains `g`
  # transitions can only through **one area** transition
  for i = 0:(h-1), a = ia, b = ia
    a == b && continue
    push!(par_nams, "gain_"*string(a)*string(b)*"_"*string(i))
  end

  # add area looses `l`
  # transitions can only through **one area** transition
  for i = 0:(h-1), a = ia
    push!(par_nams, "loss_"*string(a)*"_"*string(i))
  end

  # add q between hidden states
  for j = 0:(h-1), i = 0:(h-1)
    j == i && continue
    push!(par_nams, "q_"*string(j)*string(i))
  end

  yppar = ny == 1 ? 1 : div(ny,
    model[1]*k +
    model[2]*k +
    model[3]*k*(k-1))

  # add betas
  if model[1]
    for i = 0:(h-1), a = ia, l = Base.OneTo(yppar)
      push!(par_nams,
        "beta_lambda_"*string(l)*"_"*string(a)*"_"*string(i))
    end
  end
  if model[2]
    for i = 0:(h-1), a = ia, l = Base.OneTo(yppar)
      push!(par_nams,
        "beta_mu_"*string(l)*"_"*string(a)*"_"*string(i))
    end
  end
  if model[3]
    for i = 0:(h-1), a = ia, b = ia, l = Base.OneTo(yppar)
      a == b && continue
      push!(par_nams,
        "beta_q_"*string(l)*"_"*string(a)*string(b)*"_"*string(i))
    end
  end

  pardic::Dict{String, Int64} = Dict(par_nams[i]::String => i::Int64
    for i = Base.OneTo(lastindex(par_nams)))

  return pardic::Dict{String, Int64}
end

"""
    build_par_names(k::Int64, T::Bool)

Build dictionary for parameter names and indexes for EHiSSE.
"""
function build_par_names(k::Int64, T::Bool)

  par_nams = String[]
  # build parameters name
  for i in 0:(k-1)
    push!(par_nams, "lambda$i")
  end

  # add μ names
  for i in 0:(k-1)
    push!(par_nams, "mu$i")
  end

  # add q names
  for j in 0:(k-1), i in 0:(k-1)
    if i == j
      continue
    end
    push!(par_nams, "q$j$i")
  end

  # add betas
  if T
    for i in 0:(k-1)
      push!(par_nams, "beta$i")
    end
  else
    # add betas
    for j in 0:(k-1), i in 0:(k-1)
      if i == j
        continue
      end
      push!(par_nams, "beta$j$i")
    end
  end

  pardic = Dict(par_nams[i] => i for i in Base.OneTo(length(par_nams)))

  return pardic
end

"""
    set_constraints(constraints::NTuple{N,String},
                    pardic     ::Dict{String,Int64},
                    h          ::Int64)

Make a Dictionary linking parameter that are to be the same.
"""
function set_constraints(constraints::NTuple{N,String},
                         pardic     ::Dict{String,Int64},
                         k          ::Int64,
                         h          ::Int64,
                         ny         ::Int64,
                         model      ::NTuple{3,Bool}) where {N}

  # Dict of hidden state correspondence
  hsc = dict_hscor(k, h, ny, model)

  dcp0  = Dict{Int64,Int64}()  # constrains dictionary for base hidden state 0
  dcfp0 = Dict{Int64,Int64}()  # constrains dictionary for other hidden states
  zp    = Set(Int64[])         # zeros for base hidden state 0
  zfp   = Set(Int64[])         # zeros for other hidden states

  for c in constraints

    spl = map(x -> strip(x), split(c, '='))

    if occursin("q", spl[1])
      dcp0[pardic[spl[2]]] = pardic[spl[1]]
      continue
    end

    # if no equality
    if length(spl) < 2
      @warn "No equality found in $c"
      continue
    end

    # for parameters fixed to `0`
    if isequal(spl[end], "0")
      for i in Base.OneTo(length(spl)-1)
        hsi  = split(spl[i], "_")[end]
        if hsi == "0"
          push!(zp, pardic[spl[i]])
        else
          push!(zfp, pardic[spl[i]])
        end
      end

    # for constraints
    else
      for i in Base.OneTo(length(spl)-1)

        sp1 = spl[i]
        sp2 = spl[i+1]

        # divide
        spl1 = split(sp1, "_")
        spl2 = split(sp2, "_")

        # concatenate without hidden state
        spc1 = spc2 = ""
        for s in Base.OneTo(lastindex(spl1)-1)
          spc1 *= spl1[s]*"_"
        end
        for s in Base.OneTo(lastindex(spl2)-1)
          spc2 *= spl2[s]*"_"
        end

        # convert to hidden states
        hs1 = parse(Int64, spl1[end])
        hs2 = parse(Int64, spl2[end])

        # minmax for hidden states
        mnh, mxh = minmax(hs1, hs2)

        # which the maximum
        spmn, spmx = hs1 < hs2 ? (sp1, sp2) : (sp2, sp1)

        # set non-shared hidden states to 0
        wp = pardic[spmx]
        for j in (mnh+1):mxh
          push!(zfp, wp)
          wp = hsc[wp]
        end

        if spc1 != spc2
          # set shared hidden states equal
          for j in 0:mnh
            if iszero(j)
              dcp0[pardic[spc2*"0"]] = pardic[spc1*"0"]
            else
              dcfp0[pardic[spc2*"$j"]] = pardic[spc1*"$j"]
            end
          end
        end
      end
    end
  end

  # check double feedback loops
  dcp = Dict{Int64,Int64}()
  for (key, value) in dcp0
    if in(value, keys(dcp)) && key == dcp[value]
      continue
    end
    dcp[key] = value
  end

  dcfp = Dict{Int64,Int64}()
  for (key, value) in dcfp0
    if in(value, keys(dcfp)) && key == dcfp[value]
      continue
    end
    dcfp[key] = value
  end

  return dcp, dcfp, zp, zfp
end

"""
    dict_hscor(k::Int64, h::Int64, ny::Int64)

Create dictionary of hidden states correspondence.
"""
function dict_hscor(k::Int64, h::Int64, ny::Int64, model::NTuple{3, Bool})

 # number of covariates
  yppar = ny == 1 ? 1 : div(ny,
    model[1]*k +
    model[2]*k +
    model[3]*k*(k-1))

  # starting indices for models 2 and 3
  m2s = model[1]*k*h*yppar
  m3s = m2s + model[2]*k*h*yppar

  hsc = Dict{Int64,Int64}()

  for j in 1:(h-1)
    for i in 1:k

      # speciation
      hsc[(k+1)*j + i] = (k+1)*(j-1) + i

      # extinction
      s = (k+1)*h
      hsc[s + k*j + i] = s + k*(j-1) + i

      # gain
      s = (2k+1)*h
      for a in 1:(k-1)
        hsc[s + k*(k-1)*j + a + (k-1)*(i-1)] = s + k*(k-1)*(j-1) + a + (k-1)*(i-1)
      end

      # loss
      s = (2k+1 + k*(k-1))*h
      hsc[s + k*j + i] = s + k*(j-1) + i

      # betas
      s = (3k+1+k*(k-1))*h + h*(h-1)

      if model[1]
        for l = Base.OneTo(yppar)
          hsc[s + k*j + i + (l-1)*(i-1)] = s + k*(j-1) + i + (l-1)*(i-1)
        end
      end
      if model[2]
        for l = Base.OneTo(yppar)
          hsc[s + k*j + m2s + i + (l-1)*(i-1)] = s + k*(j-1) + m2s + i + (l-1)*(i-1)
        end
      end

      if model[3]
        for a = 1:(k-1), l = Base.OneTo(yppar)
          hsc[s + k*(k-1)*j + m3s + a + (a-1)*yypar + (i-1)*(k-1)] =
            s + k*(k-1)*(j-1) + m3s + a + (a-1)*yypar + (i-1)*(k-1)
        end
      end
    end

    # between-region speciation
    hsc[(k+1)*j + (k+1)] = (k+1)*(j-1) + (k+1)
  end

  return hsc
end

"""
  set_multivariate(mvpars::NTuple{N,String},
                   pardic::Dict{String,Int64}) where {N}

Returns all vectors for multivariate updates given mvpars.
"""
function set_multivariate(mvpars::NTuple{N,String},
                          pardic::Dict{String,Int64}) where {N}

  mvps = Array{Int64,1}[]
  for c in mvpars
    spl = map(x -> strip(x), split(c, '='))
    si = spl[1]

    if si == ""
      break
    end

    for (k,v) in pardic
      # if matches a parameter
      if occursin(Regex("^$si"), k)
        mvp = [v]
        # get area of parameter
        ar = match(r"[A-Z]", k).match
        occursin('W', ar) && continue
        # get hidden state
        hd = match(r"[0-9]$", k).match
        # get other parameters in c that has this area and hidden state
        for sj in spl
          si == sj && continue
          for (kj, vj) in pardic
            if occursin(Regex("^$sj"), kj) && occursin(ar*"_"*hd, kj)
              push!(mvp, vj)
            end
          end
        end
        push!(mvps, sort!(mvp))
      end
    end
  end

  return mvps
end

"""
    states_to_values(tipst::Dict{Int64,Int64}, S::Array{Sgh,1}, k::Int64, h::Int64)

Transform numbered tip_values to array with 1s and 0s
"""
function states_to_values(tipst::Dict{Int64,Int64}, S::Array{Sgh,1}, k::Int64)

  tip_val = Dict{Int64,Array{Float64,1}}()
  for (key, val) in tipst
    push!(tip_val, key => zeros(k))
    for s in S[val].g
      setindex!(tip_val[key], 1.0, s)
    end
  end

  return tip_val
end

"""
    states_to_values(tipst::Dict{Int64,Int64}, k::Int64)

Transform numbered tip_values to array with 1s and 0s
"""
function states_to_values(tipst::Dict{Int64,Int64}, k::Int64)
  tip_val = Dict{Int64,Array{Float64,1}}()
  for (key, val) in tipst
    push!(tip_val, key => zeros(k))
    setindex!(tip_val[key], 1.0, val)
  end

  return tip_val
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
## String utilities for vectors.

"""
Non-standard string literal `sz`.

Convert a space-separated string to an occupation vector.
"""
macro sz_str(x)
    Int[parse(Int, n) for n in split(x)]
end

"""
    to_str(v::Vector{Int})

Convert an occupation vector `v` to a space-separated string.
"""
to_str(v::AbstractVector{Int}) = join(v, " ")

ext: jl


text:
"""
  Methods for the manipulation of geometric data.

  # AUTHORSHIP
    * Author    : Eduardo J Alvarez
    * Email     : Edo.AlvarezR@gmail.com
    * Created   : Aug 2017
    * License   : MIT License
"""
module GeometricTools

using Printf
using LinearAlgebra
using Requires
using Statistics
import PyCall
import Dierckx
import Roots
import QuadGK
import HDF5

const patch = PyCall.PyNULL()
const module_path = splitdir(@__FILE__)[1]      # Path to this module
                                                # Type of multidiscretize input
const multidisctype = Vector{Tuple{Float64,Int64,Float64,Bool}}

function __init__()
    copy!(patch, PyCall.pyimport_conda("matplotlib.patches", "patches"))
    @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" include("GeometricTools_plotting.jl")
end

for header_name in ["vtk", "geometry", "misc", "gridabstract", "airfoil",
                    "surfacing", "plot3d", "vtk_shapes", "conics",
                    "statistics", "linearalgebra", "xdmf",
                    "DEPRECATED"]
    include("GeometricTools_"*header_name*".jl")
end

end # END OF MODULE


text:
using Test
using Petri
import Petri: stripnullterms, fluxes, odefunc, quotesplat
using ModelingToolkit
using MacroTools
import MacroTools: postwalk

@testset "Equality" begin
    @variables S, E, I, R
    sir = Model([S,I,R], [(S+I, 2I), (I,R)])
    @test sir == sir
    seir = Model([S,E, I,R], [(S+I, E+I), (I,R)])
    x = Model([S, E, I,R], [(S+I, 2I), (I,R)])
    y = Model([S, E, I,R], [(S+I, E+I), (I,R), (R,S)])
    @test sir != seir
    @test seir != x
    @test seir != y
end

include("stochastic.jl")

mutable struct SIRState{T,F}
    S::T
    I::T
    R::T
    β::F
    γ::F
    μ::F
end

function test1()
    @variables S, I, R, β, γ, μ
    N = +(S,I,R)
    ϕ = [(S > 0) * (I > 0),
         I > 0,
         R > 0]

    Δ = [(S~S-1, I~I+1),
         (I~I-1, R~R+1),
         (R~R-1, S~S+1)]

    Λ = [β*S*I/N,
         γ*I,
         μ*R]

    m = Model([S,I,R], Δ, Λ, ϕ)
    p = Problem(m, SIRState(100, 1, 0, 0.5, 0.15, 0.05), 150)
end

p1 = test1()

m = p1.m
map(m.Λ) do λ
    body, args = @show Petri.funcbody(λ)
end

S, I, R = m.S

Δin  = [S+I, I, R]
Δout = [2I, R, S]

Δinm = [1 1 0;
        0 1 0;
        0 0 1]

Δoutm = [0 2 0;
         0 0 1;
         1 0 0]

du = (Δoutm - Δinm)'m.Λ

answer = map(enumerate(du)) do (i, ex)
    body, args = Petri.funcbody(ex)
    body′ = stripnullterms(body)
    state = m.S[i].op.name
    :(du.$(state) = $body′)
end
f = Petri.funckit(:f, (:du, :state, :p, :t), quote $(answer...)end )

include("ode.jl")


text:
function generaCalles(ps_predio, ps_publico, anchoEspacioPublico)

    conjuntoAnchoEspPublico = unique(anchoEspacioPublico)

    ps_predio_ = polyShape.polyExpand(ps_predio, .01)
    ps_calles_ = polyShape.polyDifference_v3(ps_publico, ps_predio_)

    ps_calles = PolyShape([], 0)
    for i = 1:ps_calles_.NumRegions
        V_i = ps_calles_.Vertices[i]
        ps_calle_i = PolyShape([V_i], 1)
        for k = 1:size(V_i,1)
            q = V_i[k,:]
            if k <= size(V_i,1)-2
                p1 = V_i[k+1,:]
                p2 = V_i[k+2,:]
            elseif k == size(V_i,1)-1
                p1 = V_i[k+1,:]
                p2 = V_i[1,:]
            elseif k == size(V_i,1)
                p1 = V_i[1,:]
                p2 = V_i[2,:]
            end
            dist = poly2D.distPointLine(q, p1, p2)

            for j = conjuntoAnchoEspPublico
                if abs(dist - j) <= .111
                    ps_calle_i = polyShape.polyExpandSides(ps_calle_i, 100, k)
                end
            end
        end
        ps_calles.Vertices = push!(ps_calles.Vertices, ps_calle_i.Vertices[1])
    end

    ps_calles.NumRegions = length(ps_calles.Vertices)

    return ps_calles

end


text:
module ScalarSolvers

using JuMP
using NLPModels
using LinearAlgebra
using Printf
using State, Stopping
using Polynomials

include("includes.jl")

include("solvers.jl")

end # module


text:
using LightGraphs, MatrixDepot, SparseArrays, Arpack, MatrixNetworks

matrix_network_datasets()
A = load_matrix_network("dfs_example")
Acc,p = largest_component(A)

biggy = sprand(Int,
  10^6, 10^6,
  0.000001)
biggy_sim = (biggy + biggy')/2

using LinearAlgebra

biggy_eigs = eigs(biggy)

biggys


text:
using CSV, DataFrames, Dates, LinearAlgebra, Distributions, DelimitedFiles, SharedArrays, Plots
gr()

import Base.Threads.@spawn
# Custom package
using Jevo

# Get date to append to output file
date = Dates.format(Dates.today(), "yyyy_mm_dd")

# Define Parameters
l0 = 10
ϵ0 = 0.7
n = 4
N = 100
fl = 1/2N
f0 = 400/2N
κ_arr = [0, 0.2, 0.5, 1, 2, 5, 10, 20]#, 20, 45]
F = Jevo.fermi_fitness(epsilon=ϵ0, f0=f0, fl=fl, l0=l0)
emat = ϵ0* (ones(n, n) - Matrix{Float64}(I, n, n))

function make_histogram(data)
    x = unique(data)
    y = [sum([d == t for d in data]) for t in x]
    return x, y/sum(y)
end

function Run(N, l, emat, F, κ, l0, steps)
    E_arr = zeros(Float64, steps)
    L_arr = zeros(Float64, steps)
    pop = Jevo.mono_pop(N=N, l=l)
    Jevo.initiate!(pop, emat)
    for r in 1:steps
        Jevo.bp_substitution!(pop, emat, F)
        if rand() < κ/N
            Jevo.driver_mutation!(pop)
        end
        if rand() < 1/10
            Jevo.l_substitution!(pop, emat, F)
        end
        E_arr[r] = Jevo.get_energy(pop, emat)
        L_arr[r] = length(pop.seqs)
        if E_arr[r]/(L_arr[r]*ϵ0) > (Jevo.γ_0(n))
            Jevo.initiate!(pop, emat)
        end
    end
    return E_arr, L_arr
end

# Predefine output arrays
k_list = Array{Array{Float64, 1}, 1}()
l_list = Array{Array{Float64, 1}, 1}()
Q_list = Array{Array{Float64, 1}, 1}()
κ_list = Array{Array{Float64, 1}, 1}()

for i in 1:length(κ_arr)
    push!(l_list, Float64[])
    push!(k_list, Float64[])
    push!(Q_list, Float64[])
    push!(κ_list, Float64[])
end

# uses all availabe cores
Threads.@threads for i in 1:length(κ_arr)
    E, L = Run(N, 40, emat, F, κ_arr[i], l0, 10^8)
    pairs, hist = make_histogram([(x, y) for (x, y) in zip(E[10^6:10^8], L[10^6:10^8])])
    k_list[i] = [pair[1] for pair in pairs]
    l_list[i] = [pair[2] for pair in pairs]
    Q_list[i] = hist
    κ_list[i] = ones(length(pairs)) .* κ_arr[i]
    println("Run κ = $(κ_arr[i]) done.")
end

df = DataFrame(Q=[(Q_list...)...], k=[(k_list...)...], l=[(l_list...)...], kappa=[(κ_list...)...])

CSV.write(date * "_results.csv", df)


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct BinaryOpReq <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function BinaryOpReq(; kwargs...)
        obj = new(meta(BinaryOpReq), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end #type BinaryOpReq
function meta(::Type{BinaryOpReq})
    allflds = Pair{Symbol,Union{Type,String}}[:i1 => Int64, :i2 => Int64]
    meta(ProtoMeta(BinaryOpReq), BinaryOpReq, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
end
function Base.getproperty(obj::BinaryOpReq, name::Symbol)
    if name === :i1
        return (obj.__protobuf_jl_internal_values[name])::Int64
    elseif name === :i2
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

mutable struct BinaryOpResp <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function BinaryOpResp(; kwargs...)
        obj = new(meta(BinaryOpResp), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end #type BinaryOpResp
function meta(::Type{BinaryOpResp})
    allflds = Pair{Symbol,Union{Type,String}}[:result => Int64]
    meta(ProtoMeta(BinaryOpResp), BinaryOpResp, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
end
function Base.getproperty(obj::BinaryOpResp, name::Symbol)
    if name === :result
        return (obj.__protobuf_jl_internal_values[name])::Int64
    else
        getfield(obj, name)
    end
end

# service methods for TestMath
const _TestMath_methods = MethodDescriptor[
        MethodDescriptor("Mul", 1, BinaryOpReq, BinaryOpResp),
        MethodDescriptor("Add", 2, BinaryOpReq, BinaryOpResp)
    ] # const _TestMath_methods
const _TestMath_desc = ServiceDescriptor("TestMath", 1, _TestMath_methods)

TestMath(impl::Module) = ProtoService(_TestMath_desc, impl)

mutable struct TestMathStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    TestMathStub(channel::ProtoRpcChannel) = new(ProtoServiceStub(_TestMath_desc, channel))
end # type TestMathStub

mutable struct TestMathBlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    TestMathBlockingStub(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_TestMath_desc, channel))
end # type TestMathBlockingStub

Mul(stub::TestMathStub, controller::ProtoRpcController, inp::BinaryOpReq, done::Function) = call_method(stub.impl, _TestMath_methods[1], controller, inp, done)
Mul(stub::TestMathBlockingStub, controller::ProtoRpcController, inp::BinaryOpReq) = call_method(stub.impl, _TestMath_methods[1], controller, inp)

Add(stub::TestMathStub, controller::ProtoRpcController, inp::BinaryOpReq, done::Function) = call_method(stub.impl, _TestMath_methods[2], controller, inp, done)
Add(stub::TestMathBlockingStub, controller::ProtoRpcController, inp::BinaryOpReq) = call_method(stub.impl, _TestMath_methods[2], controller, inp)

export BinaryOpReq, BinaryOpResp, TestMath, TestMathStub, TestMathBlockingStub, Mul, Add


text:
function position_mod(x::Float64, L::Float64)
	if x < 0.0
		x += L
	elseif x > L
		x -= L
	end

	return x
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Author: Paula B. Staudt and Rafael P. Soares
#* $Id$
#*------------------------------------------------------------------
type batch_isothermic
	batch_isothermic()=begin
		new(
			batch_basic(),
			[
				:(T = _base_1.Inlet.T),
				:(P = _base_1.Inlet.P),
			],
			[
				"Isotermic","Isobaric",
			],
		)
	end
	_base_1::batch_basic
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	attributes::Dict{Symbol,Any}
end
export batch_isothermic
function setEquationFlow(in::batch_isothermic)
	addEquation(1)
	addEquation(2)
end
function atributes(in::batch_isothermic,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/batch"
	fields[:Brief]="Model of a isothermal batch reactor"
	fields[:Info]="
Assumptions
* isothermic
"
	drive!(fields,_)
	return fields
end
batch_isothermic(_::Dict{Symbol,Any})=begin
	newModel=batch_isothermic()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
module TestFilter

using DSP

include("dspfilter.jl")

# Write your package code here.

function add_one(x)
    return x + 1
end

function add_two(x)
    return x + 2
end

export add_one,
    DSPFilter,
    init,
    run_filter

end # module


text:
using GAIO

# the Henon map
a, b = 1.4, 0.3
f((x,y)) = (1 - a*x^2 + y, b*x)

center, radius = (0, 0), (3, 3)
P = BoxPartition(Box(center, radius))
F = BoxMap(f, P)
A = relative_attractor(F, P[:], steps = 16)

plot(A)


text:
export gradient!

function gradient!(top::Var; debug=false)
    sorted = topsort(top)
    isgraddefined(top) || initgrad!(top)
    for v in sorted
        f = v.creator
        f isa SymbolFunc && continue
        debug && println(typeof(f))
        gys = [output.value.grad for output in f.outputs]
        gxs = backward(v.creator, gys...)
        gxs isa Tuple || (gxs = (gxs,))
        for (x, gx) in zip(f.args, gxs)
            debug && println("\t", typeof(gx), ":", size(gx))
            isgraddefined(x) && !isnothing(gx) ? x.grad = x.grad + gx : x.grad = gx
        end
    end
    return sorted[end]
end


text:
##############################################################################
##
## FixedEffectLinearMap
##
## Model matrix of categorical variables
## mutiplied by diag(1/sqrt(∑w * interaction^2, ..., ∑w * interaction^2) (Jacobi preconditoner)
##
## We define these methods used in lsmr! (duck typing):
## eltyp
## size
## mul!
##
##############################################################################
mutable struct FixedEffectLinearMapCPU{T}
	fes::Vector{<:FixedEffect}
	scales::Vector{<:AbstractVector}
	caches::Vector{<:AbstractVector}
	tmp::Vector{Union{Nothing, <:AbstractVector}}
	nthreads::Int
end

function FixedEffectLinearMapCPU{T}(fes::Vector{<:FixedEffect}, ::Type{Val{:cpu}}, nthreads) where {T}
	scales = [zeros(T, fe.n) for fe in fes]
	caches = [zeros(T, length(fes[1].interaction)) for fe in fes]
	fecoefs = [[zeros(T, fe.n) for _ in 1:nthreads] for fe in fes]
	return FixedEffectLinearMapCPU{T}(fes, scales, caches, fecoefs, nthreads)
end

LinearAlgebra.adjoint(fem::FixedEffectLinearMapCPU) = Adjoint(fem)

function Base.size(fem::FixedEffectLinearMapCPU, dim::Integer)
	(dim == 1) ? length(fem.fes[1].refs) : (dim == 2) ? sum(fe.n for fe in fem.fes) : 1
end

Base.eltype(x::FixedEffectLinearMapCPU{T}) where {T} = T

function LinearAlgebra.mul!(fecoefs::FixedEffectCoefficients,
	Cfem::Adjoint{T, FixedEffectLinearMapCPU{T}},
	y::AbstractVector, α::Number, β::Number) where {T}
	fem = adjoint(Cfem)
	rmul!(fecoefs, β)
	for (fecoef, fe, cache, tmp) in zip(fecoefs.x, fem.fes, fem.caches, fem.tmp)
		gather!(fecoef, fe.refs, α, y, cache, tmp, fem.nthreads)
	end
	return fecoefs
end

function gather!(fecoef::AbstractVector, refs::AbstractVector, α::Number,
	y::AbstractVector, cache::AbstractVector, tmp::AbstractVector, nthreads::Integer)
	n_each = div(length(y), nthreads)
	Threads.@threads for t in 1:nthreads
		fill!(tmp[t], 0.0)
		gather!(tmp[t], refs, α, y, cache, ((t - 1) * n_each + 1):(t * n_each))
	end
	for x in tmp
		fecoef .+= x
	end
	gather!(fecoef, refs, α, y, cache, (nthreads * n_each + 1):length(y))
end

function gather!(fecoef::AbstractVector, refs::AbstractVector, α::Number,
	y::AbstractVector, cache::AbstractVector, irange::AbstractRange)
	@inbounds @simd for i in irange
		fecoef[refs[i]] += α * y[i] * cache[i]
	end
end

function LinearAlgebra.mul!(y::AbstractVector, fem::FixedEffectLinearMapCPU,
			  fecoefs::FixedEffectCoefficients, α::Number, β::Number)
	rmul!(y, β)
	for (fecoef, fe, cache) in zip(fecoefs.x, fem.fes, fem.caches)
		scatter!(y, α, fecoef, fe.refs, cache, fem.nthreads)
	end
	return y
end

function scatter!(y::AbstractVector, α::Number, fecoef::AbstractVector,
	refs::AbstractVector, cache::AbstractVector, nthreads::Integer)
	n_each = div(length(y), nthreads)
	Threads.@threads for t in 1:nthreads
		scatter!(y, α, fecoef, refs, cache, ((t - 1) * n_each + 1):(t * n_each))
	end
	scatter!(y, α, fecoef, refs, cache, (nthreads * n_each + 1):length(y))
end

function scatter!(y::AbstractVector, α::Number, fecoef::AbstractVector,
	refs::AbstractVector, cache::AbstractVector, irange::AbstractRange)
	@inbounds @simd for i in irange
		y[i] += α * fecoef[refs[i]] * cache[i]
	end
end

##############################################################################
##
## Implement AbstractFixedEffectSolver interface
##
##############################################################################

mutable struct FixedEffectSolverCPU{T} <: AbstractFixedEffectSolver{T}
	m::FixedEffectLinearMapCPU{T}
	weights::AbstractVector
	b::AbstractVector{T}
	r::AbstractVector{T}
	x::FixedEffectCoefficients{<: AbstractVector{T}}
	v::FixedEffectCoefficients{<: AbstractVector{T}}
	h::FixedEffectCoefficients{<: AbstractVector{T}}
	hbar::FixedEffectCoefficients{<: AbstractVector{T}}
end

function AbstractFixedEffectSolver{T}(fes::Vector{<:FixedEffect}, weights::AbstractWeights, ::Type{Val{:cpu}}, nthreads = Threads.nthreads()) where {T}
	m = FixedEffectLinearMapCPU{T}(fes, Val{:cpu}, nthreads)
	b = zeros(T, length(weights))
	r = zeros(T, length(weights))
	x = FixedEffectCoefficients([zeros(T, fe.n) for fe in fes])
	v = FixedEffectCoefficients([zeros(T, fe.n) for fe in fes])
	h = FixedEffectCoefficients([zeros(T, fe.n) for fe in fes])
	hbar = FixedEffectCoefficients([zeros(T, fe.n) for fe in fes])
	return update_weights!(FixedEffectSolverCPU(m, weights, b, r, x, v, h, hbar), weights)
end

function update_weights!(feM::FixedEffectSolverCPU, weights::AbstractWeights)
	for (scale, fe) in zip(feM.m.scales, feM.m.fes)
		scale!(scale, fe.refs, fe.interaction, weights)
	end
	for (cache, scale, fe) in zip(feM.m.caches, feM.m.scales, feM.m.fes)
		cache!(cache, fe.refs, fe.interaction, weights, scale)
	end
	feM.weights = weights
	return feM
end

function scale!(scale::AbstractVector, refs::AbstractVector, interaction::AbstractVector, weights::AbstractVector)
        fill!(scale, 0)
	@inbounds @simd for i in eachindex(refs)
		scale[refs[i]] += abs2(interaction[i]) * weights[i]
	end
	# Case of interaction variatble equal to zero in the category (issue #97)
	for i in 1:length(scale)
	    scale[i] = scale[i] > 0 ? (1 / sqrt(scale[i])) : 0.0
	end
end

function cache!(cache::AbstractVector, refs::AbstractVector, interaction::AbstractVector, weights::AbstractVector, scale::AbstractVector)
	@inbounds @simd for i in eachindex(cache)
		cache[i] = interaction[i] * sqrt(weights[i]) * scale[refs[i]]
	end
end

function solve_residuals!(r::AbstractVector, feM::FixedEffectSolverCPU{T}; tol::Real = sqrt(eps(T)), maxiter::Integer = 100_000) where {T}
	copyto!(feM.r, r)
	if !(feM.weights isa UnitWeights)
		feM.r .*=  sqrt.(feM.weights)
	end
	copyto!(feM.b, feM.r)
	mul!(feM.x, feM.m', feM.b, 1, 0)
        if length(feM.x.x) > 1
            x, ch = lsmr!(feM.x, feM.m, feM.b, feM.v, feM.h, feM.hbar; atol = tol, btol = tol, maxiter = maxiter)
	    iter, converged = ch.mvps + 1, ch.isconverged
        else
            iter, converged = 1, true
        end
	mul!(feM.r, feM.m, feM.x, -1, 1)
	if !(feM.weights isa UnitWeights)
		feM.r ./=  sqrt.(feM.weights)
	end
	copyto!(r, feM.r)
	return r, iter, converged
end

function solve_residuals!(X::AbstractMatrix, feM::FixedEffects.FixedEffectSolverCPU; progress_bar = true, kwargs...)
    iterations = Int[]
    convergeds = Bool[]
    bar = MiniProgressBar(header = "Demean Variables:", color = Base.info_color(), percentage = false, max = size(X, 2))
    for j in 1:size(X, 2)
    	v0 = time()
        _, iteration, converged = solve_residuals!(view(X, :, j), feM; kwargs...)
        v1 = time()
        # remove progress_bar if estimated time lower than 2sec
	    if progress_bar && (j == 1) && ((v1 - v0) * size(X, 2) <= 2)
	    	progress_bar = false
	    end
    	if progress_bar
    		bar.current = j
    	    showprogress(stdout, bar)
    	end
        push!(iterations, iteration)
        push!(convergeds, converged)
    end
    if progress_bar
    	end_progress(stdout, bar)
    end
    return X, iterations, convergeds
end

function solve_coefficients!(r::AbstractVector, feM::FixedEffectSolverCPU{T}; tol::Real = sqrt(eps(T)), maxiter::Integer = 100_000) where {T}
	copyto!(feM.b, r)
	if !(feM.weights isa UnitWeights)
		feM.b .*=  sqrt.(feM.weights)
	end
	fill!(feM.x, 0)
	x, ch = lsmr!(feM.x, feM.m, feM.b, feM.v, feM.h, feM.hbar; atol = tol, btol = tol, maxiter = maxiter)
	for (x, scale) in zip(feM.x.x, feM.m.scales)
		x .*=  scale
	end
	x = Vector{eltype(r)}[x for x in feM.x.x]
	full(normalize!(x, feM.m.fes; tol = tol, maxiter = maxiter), feM.m.fes), div(ch.mvps, 2), ch.isconverged
end


text:
using ProgressMeter

# Given an exact solution, test the convergence of the propagator, as
# a function of the amount of time steps.
function convergence(exact, y0, tmax, Na, Nb, propagator::MagnusPropagator)
    Ns = round(Int, logspace(Na, Nb, 100))
    errors = Vector{Float64}(length(Ns))
    times = Vector{Float64}(length(Ns))
    @showprogress for j = eachindex(Ns)
        tic()
        r = integrate(y0, tmax, Ns[j], propagator)
        times[j] = r[:milliseconds]
        errors[j] = sqrt(sumabs2(exact-r[:V]))
    end
    Ns,errors,times
end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()
using PackageCompiler
create_app(
  ".",
  joinpath(@__DIR__, "..", "build"),
  #precompile_execution_file=joinpath(@__DIR__, "run_example.jl"),
  filter_stdlibs=true,
  audit=true,
  force=true,
  )


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

using Statistics, Random
using StatsBase, Distributions, ArraysOfArrays, IntervalSets

@testset "NamedTupleDist" begin
    dist = @inferred NamedTupleDist(
        a = 5, b = Weibull(2, 1),
        c = -4..5,
        d = MvNormal([1.2 0.5; 0.5 2.1]),
        x = [3 4; 2 5],
        e = [Normal(1.1, 0.2)]
    )

    @test typeof(@inferred varshape(dist)) <: NamedTupleShape

    shape = varshape(dist)
    shapes = map(varshape, dist)

    for k in keys(dist)
        @test getproperty(shapes, k) == varshape(getproperty(dist, k))
    end

    X_unshaped = [0.2, -0.4, 0.3, -0.5, 0.9]
    X_shaped = shape(X_unshaped)
    @test (@inferred logpdf(unshaped(dist), X_unshaped)) == logpdf(Weibull(2, 1), 0.2) + logpdf(Uniform(-4, 5), -0.4) + logpdf(MvNormal([1.2 0.5; 0.5 2.1]), [0.3, -0.5]) + logpdf(Normal(1.1, 0.2), 0.9)
    @test (@inferred logpdf(dist, X_shaped)) == logpdf(unshaped(dist), X_unshaped)
    @test (@inferred logpdf(dist, X_shaped[])) == logpdf(unshaped(dist), X_unshaped)

    @test (@inferred mode(unshaped(dist))) == [mode(dist.b), 0.5, 0.0, 0.0, 1.1]
    @test (@inferred mode(dist)) == shape(mode(unshaped(dist)))[]

    @test (@inferred mean(unshaped(dist))) == [mean(dist.b), 0.5, 0.0, 0.0, 1.1]
    @test (@inferred mean(dist)) == shape(mean(unshaped(dist)))[]

    @test @inferred(var(unshaped(dist))) ≈ [var(dist.b), 6.75, 1.2, 2.1, 0.04]
    @test @inferred(var(dist)) == (a = 0, b = var(dist.b), c = var(dist.c), d = var(dist.d), x = var(dist.x), e = var(dist.e))

    @test begin
        ref_cov =
            [var(dist.b)  0.0   0.0  0.0 0.0;
             0.0  6.75  0.0  0.0 0.0;
             0.0  0.0   1.2  0.5 0.0;
             0.0  0.0   0.5  2.1 0.0;
             0.0  0.0   0.0  0.0 0.04 ]

        (@inferred cov(unshaped(dist))) ≈ ref_cov
    end

    @test @inferred(rand(dist)) isa NamedTuple
    @test pdf(dist, rand(dist)) > 0
    @test @inferred(rand(dist, ())) isa ShapedAsNT
    @test pdf(dist, rand(dist, ())) > 0
    @test @inferred(rand(dist, 100)) isa ShapedAsNTArray
    @test all(x -> x > 0, pdf.(Ref(dist), rand(dist, 10^3)))

    let X = X = varshape(dist).(nestedview(Array{eltype(unshaped(dist))}(undef, length(unshaped(dist)), 14)))
        @test @inferred(rand!(dist, view(X, 1))) == view(X, 1)
        @test @inferred(rand!(dist, X)) === X
    end

    let X = varshape(dist).([Array{Float64}(undef, totalndof(varshape(dist))) for i in 1:11])
        @test @inferred(rand!(dist, view(X, 1))) == view(X, 1)
        @test @inferred(rand!(dist, X)) === X
    end

    testrng() = MersenneTwister(0xaef035069e01e678)

    let X = rand(unshaped(dist), 10), Xn = nestedview(X)
        @test @inferred(Distributions._pdf(unshaped(dist), Xn[1])) == @inferred(pdf(unshaped(dist), Xn[1]))
        @test @inferred(pdf(unshaped(dist), X)) == @inferred(broadcast(pdf, Ref(unshaped(dist)), Xn))

        @test @inferred(Distributions._logpdf(unshaped(dist), Xn[1])) == @inferred(logpdf(unshaped(dist), Xn[1]))
        @test @inferred(logpdf(unshaped(dist), X)) == @inferred(broadcast(logpdf, Ref(unshaped(dist)), Xn))

        @test @inferred(insupport(unshaped(dist), Xn[1])) == true
        @test @inferred(insupport(unshaped(dist), fill(-Inf, length(Xn)))) == false
        @test @inferred(insupport(unshaped(dist), X)) == fill(true, length(Xn))
    end

    @test @inferred(rand(unshaped(dist))) isa Vector{Float64}
    @test shape(@inferred(rand(testrng(), unshaped(dist))))[] == @inferred(rand(testrng(), dist, ()))[] == @inferred(rand(testrng(), dist))
    @test @inferred(rand(unshaped(dist), 10^3)) isa Matrix{Float64}
    @test shape.(nestedview(@inferred(rand(testrng(), unshaped(dist), 10^3)))) == @inferred(rand(testrng(), dist, 10^3))

    propnames = propertynames(dist, true)
    @test propnames == (:a, :b, :c, :d, :x, :e, :_internal_distributions, :_internal_shape)
    @test @inferred(keys(dist)) == propertynames(dist, false)
    internaldists = getproperty(dist, :_internal_distributions)
    internalshape = getproperty(dist, :_internal_shape)
    @test all(i -> typeof(getindex(internaldists, i)) == typeof(getproperty(dist, keys(dist)[i])), 1:length(internaldists))
    @test all(key -> getproperty(getproperty(internalshape, key), :shape) == valshape(getproperty(internalshape, key)), keys(internalshape))

    @test @inferred(convert(NamedTupleDist, (x = 5, z = Normal()))) == NamedTupleDist(x = 5, z = Normal())
    @test @inferred(merge((a = 42,), NamedTupleDist(x = 5, z = Normal()))) == (a = 42, x = ConstValueDist(5), z = Normal())
    @test @inferred(NamedTupleDist((;dist...))) == dist
    @test @inferred(merge(dist)) === dist
    @test @inferred(merge(
        NamedTupleDist(x = Normal(), y = 42),
        NamedTupleDist(y = Normal(4, 5)),
        NamedTupleDist(z = Exponential(3), a = 4.2),
    )) == NamedTupleDist(x = Normal(), y = Normal(4, 5), z = Exponential(3), a = 4.2)
    @test @inferred(merge(NamedTupleDist(a = 42,), (x = 5, z = Normal()))) == NamedTupleDist(a = 42, x = ConstValueDist(5), z = Normal())
    @test @inferred(merge(
        NamedTupleDist(x = Normal(), y = 42),
        (y = Normal(4, 5),),
        NamedTupleDist(z = Exponential(3), a = 4.2),
    )) == NamedTupleDist(x = Normal(), y = Normal(4, 5), z = Exponential(3), a = 4.2)
end


text:
#=
simple integer knapsack problem with convex objective

minₓ f(c ⊙ x) :
b' x ≤ B
x ≥ 1
x ∈ ℤⁿ

where c > 0, b > 0, B > 0, and f : ℝⁿ → ℝ is a convex spectral function
=#

struct Knapsack <: ExampleInstance
    n::Int
    f::VecSpecExt
    relax_int::Bool # relax integrality
end

function build(inst::Knapsack)
    n = inst.n

    # generate parameters
    c = 2 * rand(n)
    b = rand(n)
    B = 3 * sum(b)

    # build model
    model = JuMP.Model()
    JuMP.@variable(model, x[1:n] >= 1)
    JuMP.@constraint(model, dot(b, x) <= B)

    if !inst.relax_int
        JuMP.set_integer.(x)
    end

    JuMP.@variable(model, y)
    JuMP.@objective(model, Min, y)
    add_homog_spectral(inst.f, n, vcat(y, c .* x), model)

    # save for use in tests
    model.ext[:x] = x
    model.ext[:c] = c

    return model
end

function test_extra(inst::Knapsack, model::JuMP.Model)
    stat = JuMP.termination_status(model)
    @test stat in OPT_OR_LIMIT
    (stat == MOI.OPTIMAL) || return

    tol = eps()^0.2
    x = JuMP.value.(model.ext[:x])
    if !inst.relax_int
        @test x ≈ round.(Int, x) atol = tol rtol = tol
    end
    @test all(>(1 - tol), x)
    # check conic feasibility
    c = JuMP.value.(model.ext[:c])
    obj_result = get_val(pos_only(c .* x), inst.f)
    @test JuMP.objective_value(model) ≈ obj_result atol = tol rtol = tol
    return
end


text:
struct ParticleSwarm{T} <: ZerothOrderOptimizer
    lower::Vector{T}
    upper::Vector{T}
    n_particles::Int
end

"""
# Particle Swarm
## Constructor
```julia
ParticleSwarm(; lower = [],
                upper = [],
                n_particles = 0)
```

The constructor takes 3 keywords:
* `lower = []`, a vector of lower bounds, unbounded below if empty or `Inf`'s
* `upper = []`, a vector of upper bounds, unbounded above if empty or `Inf`'s
* `n_particles = 0`, the number of particles in the swarm, defaults to least three

## Description
The Particle Swarm implementation in Optim.jl is the so-called Adaptive Particle
Swarm algorithm in [1]. It attempts to improve global coverage and convergence by
switching between four evolutionary states: exploration, exploitation, convergence,
and jumping out. In the jumping out state it intentially tries to take the best
particle and move it away from its (potentially and probably) local optimum, to
improve the ability to find a global optimum. Of course, this comes a the cost
of slower convergence, but hopefully converges to the global optimum as a result.

Note, that convergence is never assessed for ParticleSwarm. It will run until it
reaches the maximum number of iterations set in Optim.Options(iterations=x)`.

## References
- [1] Zhan, Zhang, and Chung. Adaptive particle swarm optimization, IEEE Transactions on Systems, Man, and Cybernetics, Part B: CyberneticsVolume 39, Issue 6 (2009): 1362-1381
"""
ParticleSwarm(; lower = [], upper = [], n_particles = 0) = ParticleSwarm(lower, upper, n_particles)

Base.summary(::ParticleSwarm) = "Particle Swarm"

mutable struct ParticleSwarmState{Tx,T} <: ZerothOrderState
    x::Tx
    iteration::Int
    lower::Tx
    upper::Tx
    c1::T # Weight variable; currently not exposed to users
    c2::T # Weight variable; currently not exposed to users
    w::T  # Weight variable; currently not exposed to users
    limit_search_space::Bool
    n_particles::Int
    X
    V
    X_best
    score::Vector{T}
    best_score::Vector{T}
    x_learn
    current_state
    iterations::Int
end

function initial_state(method::ParticleSwarm, options, d, initial_x::AbstractArray{T}) where T
    #=
    Variable X represents the whole swarm of solutions with
    the columns being the individual particles (= solutions to
    the optimization problem.)
    In each iteration the cost function is evaluated for all
    particles. For the next iteration all particles "move"
    towards their own historically best and the global historically
    best solution. The weighing coefficients c1 and c2 define how much
    towards the global or individual best solution they are pulled.

    In each iteration there is a check for an additional special
    solution which consists of the historically global best solution
    where one randomly chosen parameter is modified. This helps
    the swarm jumping out of local minima.
    =#
    n = length(initial_x)
    # TODO do we even need a lower((n) upper    ) that is different from method.lower(upper)
    # do some checks on input parameters
    @assert length(method.lower) == length(method.upper) "lower and upper must be of same length."
    if length(method.lower) > 0
        lower = copyto!(similar(initial_x), copy(method.lower))
        upper = copyto!(similar(initial_x), copy(method.upper))
        limit_search_space = true
        @assert length(lower) == length(initial_x) "limits must be of same length as x_initial."
        @assert all(upper .> lower) "upper must be greater than lower"
    else
        lower = copy(initial_x)
        upper = copy(initial_x)
        limit_search_space = false
    end

    if method.n_particles > 0
        if method.n_particles < 3
          @warn("Number of particles is set to 3 (minimum required)")
          n_particles = 3
        else
          n_particles = method.n_particles
        end
    else
      # user did not define number of particles
       n_particles = maximum([3, length(initial_x)])
    end
    c1 = T(2)
    c2 = T(2)
    w = T(1)

    X = Array{T,2}(undef, n, n_particles)
    V = Array{T,2}(undef, n, n_particles)
    X_best = Array{T,2}(undef, n, n_particles)
    dx = zeros(T, n)
    score = zeros(T, n_particles)
    x = copy(initial_x)
    best_score = zeros(T, n_particles)
    x_learn = copy(initial_x)

    current_state = 0

    value!!(d, initial_x)
    score[1] = value(d)

    # if search space is limited, spread the initial population
    # uniformly over the whole search space
    if limit_search_space
        for i in 1:n_particles
            for j in 1:n
                ww = upper[j] - lower[j]
                X[j, i] = lower[j] + ww * rand(T)
                X_best[j, i] = X[j, i]
                V[j, i] = ww * (rand(T) * T(2) - T(1)) / 10
            end
        end
    else
        for i in 1:n_particles
            for j in 1:n
                if i == 1
                    if abs(initial_x[i]) > T(0)
                        dx[j] = abs(initial_x[i])
                    else
                        dx[j] = T(1)
                    end
                end
                X[j, i] = initial_x[j] + dx[j] * rand(T)
                X_best[j, i] = X[j, i]
                V[j, i] = abs(X[j, i]) * (rand(T) * T(2) - T(1))
            end
        end
    end

    for j in 1:n
        X[j, 1] = initial_x[j]
        X_best[j, 1] = initial_x[j]
    end

    for i in 2:n_particles
        score[i] = value(d, X[:, i])
    end

    ParticleSwarmState(
        x,
        0,
        lower,
        upper,
        c1,
        c2,
        w,
        limit_search_space,
        n_particles,
        X,
        V,
        X_best,
        score,
        best_score,
        x_learn,
        0,
        options.iterations)
end

function update_state!(f, state::ParticleSwarmState{T}, method::ParticleSwarm) where T
    n = length(state.x)

    if state.iteration == 0
        copyto!(state.best_score, state.score)
        f.F = Base.minimum(state.score)
    end
    f.F = housekeeping!(state.score,
                              state.best_score,
                              state.X,
                              state.X_best,
                              state.x,
                              value(f),
                              state.n_particles)
    # Elitist Learning:
    # find a new solution named 'x_learn' which is the current best
    # solution with one randomly picked variable being modified.
    # Replace the current worst solution in X with x_learn
    # if x_learn presents the new best solution.
    # In all other cases discard x_learn.
    # This helps jumping out of local minima.
    worst_score, i_worst = findmax(state.score)
    for k in 1:n
        state.x_learn[k] = state.x[k]
    end
    random_index = rand(1:n)
    random_value = randn()
    sigma_learn = 1 - (1 - 0.1) * state.iteration / state.iterations

    r3 = randn() * sigma_learn

    if state.limit_search_space
        state.x_learn[random_index] = state.x_learn[random_index] + (state.upper[random_index] - state.lower[random_index]) / 3.0 * r3
    else
        state.x_learn[random_index] = state.x_learn[random_index] + state.x_learn[random_index] * r3
    end

    if state.limit_search_space
        if state.x_learn[random_index] < state.lower[random_index]
            state.x_learn[random_index] = state.lower[random_index]
        elseif state.x_learn[random_index] > state.upper[random_index]
            state.x_learn[random_index] = state.upper[random_index]
        end
    end

    score_learn = value(f, state.x_learn)
    if score_learn < f.F
        f.F = score_learn * 1.0
        for j in 1:n
            state.X_best[j, i_worst] = state.x_learn[j]
            state.X[j, i_worst] = state.x_learn[j]
            state.x[j] = state.x_learn[j]
        end
        state.score[i_worst] = score_learn
        state.best_score[i_worst] = score_learn
    end

    # TODO find a better name for _f (look inthe paper, it might be called f there)
    state.current_state, _f = get_swarm_state(state.X, state.score, state.x, state.current_state)
    state.w, state.c1, state.c2 = update_swarm_params!(state.c1, state.c2, state.w, state.current_state, _f)
    update_swarm!(state.X, state.X_best, state.x, n, state.n_particles, state.V, state.w, state.c1, state.c2)

    if state.limit_search_space
        limit_X!(state.X, state.lower, state.upper, state.n_particles, n)
    end
    compute_cost!(f, state.n_particles, state.X, state.score)

    state.iteration += 1
    false
end

function update_swarm!(X::AbstractArray{Tx}, X_best, best_point, n, n_particles, V,
                       w, c1, c2) where Tx
  # compute new positions for the swarm particles
  for i in 1:n_particles
      for j in 1:n
          r1 = rand(Tx)
          r2 = rand(Tx)
          vx = X_best[j, i] - X[j, i]
          vg = best_point[j] - X[j, i]
          V[j, i] = V[j, i]*w + c1*r1*vx + c2*r2*vg
          X[j, i] = X[j, i] + V[j, i]
      end
    end
end

function get_mu_1(f::Tx) where Tx
    if Tx(0) <= f <= Tx(4)/10
        return Tx(0)
    elseif Tx(4)/10 < f <= Tx(6)/10
        return Tx(5) * f - Tx(2)
    elseif Tx(6)/10 < f <= Tx(7)/10
        return Tx(1)
    elseif Tx(7)/10 < f <= Tx(8)/10
        return -Tx(10) * f + Tx(8)
    else
        return Tx(0)
    end
end

function get_mu_2(f::Tx) where Tx
    if Tx(0) <= f <= Tx(2)/10
        return Tx(0)
    elseif Tx(2)/10 < f <= Tx(3)/10
        return Tx(10) * f - Tx(2)
    elseif Tx(3)/10 < f <= Tx(4)/10
        return Tx(1)
    elseif Tx(4)/10 < f <= Tx(6)/10
        return -Tx(5) * f + Tx(3)
    else
        return Tx(0)
    end
end

function get_mu_3(f::Tx) where Tx
    if Tx(0) <= f <= Tx(1)/10
        return Tx(1)
    elseif Tx(1)/10 < f <= Tx(3)/10
        return -Tx(5) * f + Tx(3)/2
    else
        return Tx(0)
    end
end

function get_mu_4(f::Tx) where Tx
    if Tx(0) <= f <= Tx(7)/10
        return Tx(0)
    elseif Tx(7)/10 < f <= Tx(9)/10
        return Tx(5) * f - Tx(7)/2
    else
        return Tx(1)
    end
end

function get_swarm_state(X::AbstractArray{Tx}, score, best_point, previous_state) where Tx
    # swarm can be in 4 different states, depending on which
    # the weighing factors c1 and c2 are adapted.
    # New state is not only depending on the current swarm state,
    # but also from the previous state.
    n, n_particles = size(X)
    f_best, i_best = findmin(score)
    d = zeros(Tx, n_particles)
    for i in 1:n_particles
        dd = Tx(0.0)
        for k in 1:n_particles
            for dim in 1:n
                @inbounds ddd = (X[dim, i] - X[dim, k])
                dd += ddd * ddd
            end
        end
        d[i] = sqrt(dd)
    end
    dg = d[i_best]
    dmin = Base.minimum(d)
    dmax = Base.maximum(d)

    f = (dg - dmin) / max(dmax - dmin, sqrt(eps(Tx)))

    mu = zeros(Tx, 4)
    mu[1] = get_mu_1(f)
    mu[2] = get_mu_2(f)
    mu[3] = get_mu_3(f)
    mu[4] = get_mu_4(f)
    best_mu, i_best_mu = findmax(mu)
    current_state = 0

    if previous_state == 0
        current_state = i_best_mu
    elseif previous_state == 1
        if mu[1] > 0
            current_state = 1
        else
          if mu[2] > 0
              current_state = 2
          elseif mu[4] > 0
              current_state = 4
          else
              current_state = 3
          end
        end
    elseif previous_state == 2
        if mu[2] > 0
            current_state = 2
        else
          if mu[3] > 0
              current_state = 3
          elseif mu[1] > 0
              current_state = 1
          else
              current_state = 4
          end
        end
    elseif previous_state == 3
        if mu[3] > 0
            current_state = 3
        else
          if mu[4] > 0
              current_state = 4
          elseif mu[2] > 0
              current_state = 2
          else
              current_state = 1
          end
        end
    elseif previous_state == 4
        if mu[4] > 0
            current_state = 4
        else
            if mu[1] > 0
                current_state = 1
            elseif mu[2] > 0
                current_state = 2
            else
                current_state = 3
            end
        end
    end
    return current_state, f
end

function update_swarm_params!(c1, c2, w, current_state, f::T) where T

    delta_c1 = T(5)/100 + rand(T) / T(20)
    delta_c2 = T(5)/100 + rand(T) / T(20)

    if current_state == 1
        c1 += delta_c1
        c2 -= delta_c2
    elseif current_state == 2
        c1 += delta_c1 / 2
        c2 -= delta_c2 / 2
    elseif current_state == 3
        c1 += delta_c1 / 2
        c2 += delta_c2 / 2
    elseif current_state == 4
        c1 -= delta_c1
        c2 -= delta_c2
    end

    if c1 < T(3)/2
        c1 = T(3)/2
    elseif c1 > T(5)/2
        c1 = T(5)/2
    end

    if c2 < T(3)/2
        c2 = T(5)/2
    elseif c2 > T(5)/2
        c2 = T(5)/2
    end

    if c1 + c2 > T(4)
        c_total = c1 + c2
        c1 = c1 / c_total * 4
        c2 = c2 / c_total * 4
    end

    w = 1 / (1 + T(3)/2 * exp(-T(26)/10 * f))
    return w, c1, c2
end

function housekeeping!(score, best_score, X, X_best, best_point,
                       F, n_particles)
    n = size(X, 1)
    for i in 1:n_particles
        if score[i] <= best_score[i]
            best_score[i] = score[i]
            for k in 1:n
                X_best[k, i] = X[k, i]
            end
            if score[i] <= F
                for k in 1:n
                  	best_point[k] = X[k, i]
                end
              	F = score[i]
            end
        end
    end
    return F
end

function limit_X!(X, lower, upper, n_particles, n)
    # limit X values to boundaries
    for i in 1:n_particles
        for j in 1:n
            if X[j, i] < lower[j]
              	X[j, i] = lower[j]
            elseif X[j, i] > upper[j]
              	X[j, i] = upper[j]
            end
        end
    end
    nothing
end

function compute_cost!(f,
                       n_particles::Int,
                       X::Matrix,
                       score::Vector)

    for i in 1:n_particles
        score[i] = value(f, X[:, i])
    end
    nothing
end


text:
export prox_l2s!

function prox_l2s!(
  x   ::Vector{TF},
  rho ::TF,
  m   ::Vector{TF}
  ) where {TF<:Real}

  x .= (x .* rho .+ m) ./ (rho .+ 1.0);
  return x
end


text:
using Plots, QuadGK, ColorSchemes

# Vector e.g. [1, 4, 5] or collected Range e.g. collect(1:5)
inputstates = [9,15,30]
sort!(inputstates)

function constant(n)
    1/sqrt(size(inputstates, 1))
    #2/(n*pi)*(1-cos(n*pi/2))
end

function calculate(states)
    dx = 0.1
    dt = 1/30
    elapse = 10
    a = 10
    m = 0.2

    h = 1.054571817
    xran = 0:dx:a
    yran = 0:dx:elapse
    lbls = ["real" "imaginary"]
    clrs = ColorSchemes.GnBu_9

    Cₙ = [constant(i) for i in states]

    Ψₙ(x) = [sin(pi*x*i/a)*sqrt(2/a) for i in states]
    dΨₙ(x) = [cos(pi*x*i/a)*sqrt(2/a)*pi*i/a for i in states]

    φₙ(t) = [exp(-im*(i*pi)^2*h*t/(2*m*a^2)) for i in states]

    radii = [[0,abs(i)] for i in Cₙ]

    Eₙ = [(i*pi*h)^2/2*m*a^2 for i in states]

    anim = @animate for ti in 0:dt:elapse

            Ψ(x) = sum(Ψₙ(x).*φₙ(ti).*Cₙ)
            dΨ(x) = sum(dΨₙ(x).*φₙ(ti).*Cₙ)

            ref(x) = real(Ψ(x))
            imf(x) = imag(Ψ(x))
            ρ(x) = abs2(Ψ(x))
            J(x) = real(im*h/(2*m)*(Ψ(x)*conj(dΨ(x))-conj(Ψ(x))*dΨ(x)))
            xΨ(x) = x*ρ(x)
            x²Ψ(x) = x^2*ρ(x)
            expcX = first(quadgk(xΨ, 0, a))
            expcX² = first(quadgk(x²Ψ, 0, a))
            σₓ = sqrt(expcX²-expcX^2)

            theta = [[0,i] for i in angle.(φₙ(ti))]

            p1 = plot(xran, [ref, imf], ylims = [-1,1], color = [clrs[5] clrs[3]], labels = lbls, title = "Ψ(x, t)")
            p2 = plot(xran, ρ, ylims = [0,1], color = clrs[6], labels = false, title = "∣Ψ(x, t)∣²", fill = true)
            p3 = plot(theta, radii, proj =:polar, lims = [0,1], showaxis = false, labels = false, palette = [clrs[4], clrs[5], clrs[6], clrs[7], clrs[8]], title = "φ(t)")
            scatter!(p2, [expcX], [0.5], color = clrs[5], xerror = σₓ, markerstrokecolor=clrs[8], labels = false)

            p5 = plot(xran, J, ylims = [-1.5, 1.5], labels = false, title = "J(x,t)", color = clrs[6])

            l1 = @layout[a b; c{0.4h} d]

            plot(p1, p2, p3, p5, fontfamily = :serif, layout = l1)

    end

    Eplot = vcat(((states[begin]-1)*pi*h)^2/2*m*a^2, Eₙ, ((states[end]+1)*pi*h)^2/2*m*a^2)
    Cplot = vcat(0, abs2.(Cₙ), 0)
    Psi(x, t) = abs2(sum(Ψₙ(x).*φₙ(t).*Cₙ))
    l2 = @layout[a{0.3w} [b; c]]

    p4 = plot(Cplot, Eplot, xlims = [0, 1], yticks = :log10, markershapes = :circle, ylabel = "E", xlabel = "∣C∣²", labels = false, line = clrs[4], markercolors = clrs[4])
    hline!(p4, [sum(abs2.(Cₙ).*Eₙ)], color = clrs[9], label = false)
    p8 = surface(xran, yran, Psi, camera = [15, 75], seriescolor = cgrad(clrs), ylabel = "t", zlabel = "∣Ψ(x, t)∣²")
    p6 = contourf(xran, yran, Psi, seriescolor = cgrad(clrs), ylabel = "t")

    stats = plot(p4, p8, p6, layout = l2, fontfamily = :serif)

    savefig(stats, "isw1_$inputstates _stats.png")

    gif(anim, "isw1_$inputstates.gif", fps = 1/dt)
end

@time calculate(inputstates)


text:
#### PH Types Used in Interface ####

"""
Struct containing all information needed by PH about a variable (along with the variable id) to run the PH algorithm.

**Fields**

* `name::String` : User specified name of the variable
* `is_integer::Bool` : True if the variable is an integer or binary
"""
struct VariableInfo
    name::String
    is_integer::Bool # true if variable is integer or binary
end

#### PH Subproblem Interface ####

"""
Abstract type for ProgressiveHedging subproblem types.

A concrete subtype handles all details of creating, solving and updating subproblems for ProgressiveHedging (PH).  In particular, it handles all interaction with any modeling language that the subproblem is written in.  See `JuMPSubproblem` for an implementation that uses JuMP as the optimization modeling language.

Variables and their values are identified and exchanged between PH and a Subproblem type using the `VariableID` type.  A unique `VariableID` is associated with each variable in the subproblem by the Subproblem implementation by using the `ScenarioID` of the subproblem as well as the `StageID` to which this variable belongs.  This combination uniquely identifies a node in the scenario tree to which the variable can be associated.  Variables associated with the same node in a scenario tree and sharing the same name are assumed to be consensus variables whose optimal value is determined by PH.  The final component of a `VariableID` is an `Index` which is just a counter assigned to a variable to differentiate it from other variables at the same node.  See `VariableID` type for more details.

Any concrete subtype **must** implement the following functions:
* [`add_ph_objective_terms`](@ref)
* [`objective_value`](@ref)
* [`report_values`](@ref)
* [`report_variable_info`](@ref)
* [`solve_subproblem`](@ref)
* [`update_ph_terms`](@ref)
See help strings on each function for details on arguments, returned objects and expected performance of each function. See JuMPSubproblem for an example using JuMP.

To use `warm_start=true` the concrete subtype must also implement
* [`warm_start`](@ref)
See help on `warm_start` for more information.  See JuMPSubproblem for an example using JuMP.

To use the extensive form functionality, the concrete subtype must implement
* [`ef_copy_model`](@ref)
* [`ef_node_dict_constructor`](@ref)
See the help on the functions for more details. See JuMPSubproblem for an example using JuMP. Note that the extensive form model is always constructed as a `JuMP.Model` object.

To use penalty parameter types other than `ScalarPenaltyParameter`, concrete subproblem types also need to implement
* [`report_penalty_info`](@ref)
See the help on individual penalty parameter types and `report_penalty_info` for more details.

To use lower bound computation functionality, the concrete subtype must also implement
* [`add_lagrange_terms`](@ref)
* [`update_lagrange_terms`](@ref)
See the help on these functions for more details.
"""
abstract type AbstractSubproblem end

## Required Interface Functions ##

# The below functions must be implemented by any new subtype of AbstractSubproblem.
# Note they all default to throwing an `UnimplementedError`.

"""
    add_ph_objective_terms(as::AbstractSubproblem,
                           vids::Vector{VariableID},
                           r::Union{Float64,Dict{VariableID,Float64}}
                           )::Dict{VariableID,Float64}

Create model variables for Lagrange multipliers and hat variables and add Lagrange and proximal terms to the objective function.

Returns mapping of the variable to any information needed from the subproblem to compute the penalty parameter. Only used if `is_subproblem_dependent(typeof(r))` returns `true`.

**Arguments**

* `as::AbstractSubproblem` : subproblem object (replace with appropriate type)
* `vids::Vector{VariableID}` : list of `VariableIDs` which need ph terms created
* `r::AbstractPenaltyParameter` : penalty parameter on quadratic term
"""
function add_ph_objective_terms(as::AbstractSubproblem,
                                vids::Vector{VariableID},
                                r::Union{Float64,Dict{VariableID,Float64}},
                                )::Nothing
    throw(UnimplementedError("add_ph_objective_terms is unimplemented"))
end

"""
    objective_value(as::AbstractSubproblem)::Float64

Return the objective value of the solved subproblem.

**Arguments**

* `as::AbstractSubproblem` : subproblem object (replace with appropriate type)
"""
function objective_value(as::AbstractSubproblem)::Float64
    throw(UnimplementedError("objective_value is unimplemented"))
end

"""
    report_penalty_info(as::AbstractSubproblem,
                        pp<:AbstractPenaltyParameter,
                        )::Dict{VariableID,Float64}

Returns mapping of variable to a value used to compute the penalty parameter.

This function must be implemented for any concrete penalty parameter type for which `is_subproblem_dependent` returns true and any concrete subproblem type wishing to use that penalty parameter. The mapping must contain all subproblem values needed for processing by `process_penalty_subproblem`. The returned mapping is handed unaltered to `process_penalty_subproblem` as the `subproblem_dict` argument.

**Arguments**

*`as::AbstractSubproblem` : subproblem object (replace with appropriate type)
*`pp<:AbstractPenaltyParameter` : penalty parameter type
"""
function report_penalty_info(as::AbstractSubproblem,
                             pp::Type{P},
                             )::Dict{VariableID,Float64} where P <: AbstractPenaltyParameter
    throw(UnimplementedError("report_penalty_info is unimplmented for subproblem of type $(typeof(as)) and penalty parameter of type $(pp)."))
end

"""
    report_variable_info(as::AbstractSubproblem,
                         st::ScenarioTree
                         )::Dict{VariableID, VariableInfo}

Assign `VariableID`s to all model variables and build a map from those ids to the required variable information (e.g., variable name).  See `VariableInfo` help for details on required variable information.

**Arguments**

* `as::AbstractSubproblem` : subproblem object (replace with appropriate type)
* `st::ScenarioTree` : scenario tree for the entire PH problem
"""
function report_variable_info(as::AbstractSubproblem,
                              st::ScenarioTree
                              )::Dict{VariableID, VariableInfo}
    throw(UnimplementedError("report_variable_info is unimplemented"))
end

"""
    report_values(as::AbstractSubproblem,
                  vars::Vector{VariableID}
                  )::Dict{VariableID, Float64}

Return the variable values specified by `vars`.

**Arguments**

* `as::AbstractSubproblem` : subproblem object (replace with appropriate type)
* `vars::Vector{VariableID}` : collection of `VariableID`s to gather values for
"""
function report_values(as::AbstractSubproblem,
                       vars::Vector{VariableID}
                       )::Dict{VariableID, Float64}
    throw(UnimplementedError("report_values is unimplemented"))
end

"""
    solve_subproblem(as::AbstractSubproblem)::MOI.TerminationStatusCode

Solve the subproblem specified by `as` and return the status code.

**Arguments**

* `as::AbstractSubproblem` : subproblem object (replace with appropriate type)
"""
function solve_subproblem(as::AbstractSubproblem)::MOI.TerminationStatusCode
    throw(UnimplementedError("solve_subproblem is unimplemented"))
end

"""
    update_ph_terms(as::AbstractSubproblem,
                    w_vals::Dict{VariableID,Float64},
                    xhat_vals::Dict{VariableID,Float64}
                    )::Nothing

Update the values of the PH variables in this subproblem with those given by `w_vals` and `xhat_vals`.

**Arguments**

* `as::AbstractSubproblem` : subproblem object (replace with appropriate type)
* `w_vals::Dict{VariableID,Float64}` : Values to update Lagrangian variables with
* `xhat_vals::Dict{VariableID,Float64}` : Values to update proximal variables with
"""
function update_ph_terms(as::AbstractSubproblem,
                         w_vals::Dict{VariableID,Float64},
                         xhat_vals::Dict{VariableID,Float64}
                         )::Nothing
    throw(UnimplementedError("update_ph_terms is unimplemented"))
end

"""
    warm_start(as::AbstractSubproblem)::Nothing

Use the values of previous solves non-PH variables as starting points of the next solve.

**Arguments**

* `as::AbstractSubproblem` : subproblem object (replace with appropriate type)
"""
function warm_start(as::AbstractSubproblem)::Nothing
    throw(UnimplementedError("warm_start is unimplemented"))
end

## Optional Interface Functions ##

"""
    add_lagrange_terms(as::AbstractSubproblem,
                       vids::Vector{VariableID}
                       )::Nothing

Create model variables for Lagrange multipliers and adds the corresponding terms to the objective function.

**Arguments**

* `as::AbstractSubproblem` : subproblem object (replace with appropriate type)
* `vids::Vector{VariableID}` : list of `VariableIDs` which need Lagrange terms

See also [`add_ph_objective_terms`](@ref).
"""
function add_lagrange_terms(as::AbstractSubproblem,
                            vids::Vector{VariableID}
                            )::Nothing
    throw(UnimplementedError("add_lagrange_terms is unimplemented"))
end

"""
    update_lagrange_terms(as::AbstractSubproblem,
                          w_vals::Dict{VariableID,Float64}
                          )::Nothing

Update the values of the dual variables in this subproblem with those given by `w_vals`.

**Arguments**

* `as::AbstractSubproblem` : subproblem object (replace with appropriate type)
* `w_vals::Dict{VariableID,Float64}` : values to update Lagrange dual variables with

See also [`update_ph_terms`](@ref).
"""
function update_lagrange_terms(as::AbstractSubproblem,
                               w_vals::Dict{VariableID,Float64}
                               )::Nothing
    throw(UnimplementedError("update_lagrange_terms is unimplemented"))
end

# These functions are required only if an extensive form construction is desired

"""
    ef_copy_model(destination::JuMP.Model,
                  original::AbstractSubproblem,
                  scid::ScenarioID,
                  scen_tree::ScenarioTree,
                  node_dict::Dict{NodeID, Any}
                  )

Copy the subproblem described by `original` to the extensive form model `destination`.

**Arguments**

* `destination::JuMP.Model` : extensive form of the model that is being built
* `original::AbstractSubproblem` : subproblem object (replace with appropriate type)
* `scid::ScenarioID` : `ScenarioID` corresponding to this subproblem
* `scen_tree::ScenarioTree` : scenario tree for the entire PH problem
* `node_dict::Dict{NodeID, Any}` : dictionary for transferring nodal information from one submodel to another
"""
function ef_copy_model(destination::JuMP.Model,
                       original::AbstractSubproblem,
                       scid::ScenarioID,
                       scen_tree::ScenarioTree,
                       node_dict::Dict{NodeID, Any}
                       )
    throw(UnimplementedError("ef_copy_model is unimplemented"))
end

"""
    ef_node_dict_constructor(::Type{S}) where S <: AbstractSubproblem

Construct dictionary that is used to carry information between subproblems for `ef_copy_model`.

**Arguments**

* `::Type{S}` : Subproblem type
"""
function ef_node_dict_constructor(::Type{S}) where S <: AbstractSubproblem
    throw(UnimplementedError("ef_node_dict_constructor is unimplemented"))
end


text:
function subbac!(u::Matrix{Float64}, b::Vector{Float64})
# Back-substitution on an Upper Triangle
  n = size(u, 1)
  for i in n:-1:1
    total = b[i]
    if i < n
      for j in i+1:n
        total -= u[i, j] * b[j]
      end
    end
    b[i] = total / u[i, i]
  end
end


text:
using COSMO, Test, LinearAlgebra, Statistics, Random

# This test is precision agnostic and will be run with TestFloat precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end
mutable struct TestProblem
    P
    q
    constraints
end

@testset "Simple Tests" begin

for TestFloat in UnitTestFloats
    tol = TestFloat(1e-3)
    rng = Random.MersenneTwister(41)

    function simpleQP(Type::Type{T}) where {T <: AbstractFloat}
        P = T[4. 1;1 2]
        q = T[1; 1.]
        A = T[1. 1;1 0; 0 1]
        l = T[1.;0;0]
        u = T[1.;0.7;0.7]

        constraint1 = COSMO.Constraint(-A, u, COSMO.Nonnegatives)
        constraint2 = COSMO.Constraint(A, -l, COSMO.Nonnegatives)
        constraints = [constraint1; constraint2]
        return TestProblem(P, q, constraints)
    end

    @testset "Simple Tests (T = $(TestFloat))" begin

        @testset "Simple QP" begin
            p = simpleQP(TestFloat)
            model = COSMO.Model{TestFloat}()
            assemble!(model, p.P, p.q, p.constraints, settings = COSMO.Settings{TestFloat}(verbose=true))

            res = COSMO.optimize!(model);

            @test res.status == :Solved
            @test isapprox(norm(res.x - TestFloat[0.3; 0.7]), zero(TestFloat), atol=tol)
            @test isapprox(res.obj_val, TestFloat(1.8800000298331538), atol=tol)

            @test typeof(model) == COSMO.Model{TestFloat}
            @test typeof(model.p) == COSMO.ProblemData{TestFloat}
            @test typeof(model.settings) == COSMO.Settings{TestFloat}
            @test typeof(res) == COSMO.Result{TestFloat}
            @test typeof(res.x) <: AbstractVector{TestFloat}
            @test typeof(res.y) <: AbstractVector{TestFloat}
            @test typeof(res.s) <: AbstractVector{TestFloat}
        end

        @testset "update_max_iter" begin
            p = simpleQP(TestFloat)
            settings = COSMO.Settings{TestFloat}(max_iter=20)
            model = COSMO.Model{TestFloat}()
            assemble!(model, p.P, p.q, p.constraints, settings = settings)

            res = COSMO.optimize!(model);
            @test res.status == :Max_iter_reached
        end

        @testset "update_check_termination" begin
             p = simpleQP(TestFloat)
            settings = COSMO.Settings{TestFloat}(check_termination=100000)
            model = COSMO.Model{TestFloat}()
            assemble!(model, p.P, p.q, p.constraints, settings = settings)

            res = COSMO.optimize!(model);

            @test res.status == :Max_iter_reached

        end

        @testset "timelimit" begin
            p = simpleQP(TestFloat)
            settings = COSMO.Settings{TestFloat}(time_limit = 1, check_termination = 100000000, max_iter = 10000000)
            model = COSMO.Model{TestFloat}()
            assemble!(model, p.P, p.q, p.constraints, settings = settings)

            res = COSMO.optimize!(model);
            @test res.status == :Time_limit_reached
        end

         @testset "warm_starting" begin
            p = simpleQP(TestFloat)
            settings = COSMO.Settings{TestFloat}(check_termination = 1)
            model = COSMO.Model{TestFloat}()
            assemble!(model, p.P, p.q, p.constraints, settings = settings)

            res1 = COSMO.optimize!(model);
            n = 2
            m = 6
            x0 = res1.x
            y0 = res1.y
            s0 = res1.s
            model = COSMO.Model{TestFloat}()
            assemble!(model, p.P, p.q, p.constraints, settings = settings)
            COSMO.warm_start_primal!(model, x0)
            COSMO.warm_start_dual!(model, y0)

            # check that warm starting was applied correctly
            @test model.vars.x == x0
            @test model.vars.μ == -y0
            @test norm(model.vars.s.data - s0, 2) < 1e-4

            # actually, warm start with a bit of noise
            COSMO.warm_start_primal!(model, x0 .+ TestFloat(0.01) * rand(rng, TestFloat, n))
            COSMO.warm_start_dual!(model, y0 .+ TestFloat(0.01) * rand(rng, TestFloat, m))

            res2 = COSMO.optimize!(model);

            @test res1.status == :Solved && res2.status == :Solved && res2.iter < res1.iter

            @test_throws DimensionMismatch COSMO.warm_start_primal!(model, rand(TestFloat,4))
            @test_throws DimensionMismatch COSMO.warm_start_dual!(model, rand(TestFloat, 2))
        end

    end
end # UnitTestFloats
end
nothing


text:
using SafeTestsets

@safetestset "Recursive Search" begin
    using TwoDotsModels, TwoDots, Test
    game = Game()
    map(x->x.color = "green", game.dots)
    game.dots[1,1].color = "red"
    game.dots[2,1].color = "red"
    game.dots[2,2].color = "red"
    game.dots[2,3].color = "red"

    game.dots[3,1].color = "red"
    game.dots[4,1].color = "red"
    game.dots[4,2].color = "red"
    game.dots[4,3].color = "red"
    options = search(game.dots, game.dots[1,1])

    check(dot, r, c) = dot.row == r && dot.col == c
    indices = [
        (1,1),
        (2,1),
        (3,1),
        (4,1),
        (4,2),
        (4,3)
    ]
    @test map((v,t)->check(v, t...), options[1], indices) |> all

    indices = [
        (1,1),
        (2,1),
        (3,1),
        (4,1),
        (4,2),
    ]
    @test map((v,t)->check(v, t...), options[2], indices) |> all

    indices = [
        (1,1),
        (2,1),
        (3,1),
        (4,1),
    ]
    @test map((v,t)->check(v, t...), options[3], indices) |> all

    indices = [
        (1,1),
        (2,1),
        (3,1),
    ]
    @test map((v,t)->check(v, t...), options[4], indices) |> all

    indices = [
        (1,1),
        (2,1),
        (2,2),
        (2,3),
    ]
    @test map((v,t)->check(v, t...), options[5], indices) |> all

    indices = [
        (1,1),
        (2,1),
        (2,2),
    ]
    @test map((v,t)->check(v, t...), options[6], indices) |> all

    indices = [
        (1,1),
        (2,1),
    ]
    @test map((v,t)->check(v, t...), options[7], indices) |> all

    @test length(options) == 7
end


text:
module TimeMarching

# To do
# - Make use of in-place rhs operators, using same form as OrdinaryDiffEq

import ViscousFlow:@get

#using Compat
#using Compat: hasmethod

using LinearAlgebra

import LinearAlgebra: ldiv!

export System, Constrained, Unconstrained, RK, IFRK, IFHERK, r₁, r₂, B₂, B₁ᵀ,
          plan_constraints

"Abstract type for a system of ODEs"
abstract type System{C} end

const Constrained = true
const Unconstrained = false

# Functions that get extended by individual systems
function r₁ end
function r₂ end
function B₂ end
function B₁ᵀ end
function plan_constraints end

struct RKParams{N}
  c::Vector{Float64}
  a::Matrix{Float64}
end

using ..SaddlePointSystems

const RK31 = RKParams{3}([0.5, 1.0, 1.0],
                      [1/2        0        0
                       √3/3 (3-√3)/3        0
                       (3+√3)/6    -√3/3 (3+√3)/6])

const Euler = RKParams{1}([1.0],ones(1,1))

include("timemarching/rk.jl")
include("timemarching/ifrk.jl")
include("timemarching/ifherk.jl")

end