text:
#=

  Breaking News puzzle in ConstraintSolver.jl

  From http://brownbuffalo.sourceforge.net/BreakingNewsClues.html
  """
  Title: Breaking News
  Author: Faith Johnson
  Publication: Dell Logic Puzzles
  Issue: April, 1998
  Page: 9
  Stars: 1

  The Daily Galaxy sent its four best reporters
      (Corey, Jimmy, Lois, and Perry)
  to different locations
      (Bayonne, New Hope, Port Charles, and South Amboy)
  to cover four breaking news events
      (30-pound baby, blimp launching, skyscraper dedication, and
       beached whale).
  Their editor is trying to remember where each of the reporters is.
  Can you match the name of each reporter with the place he or she
  was sent, and the event that each covered?

  1. The 30-pound baby wasn't born in South Amboy or New Hope.
  2. Jimmy didn't go to Port Charles.
  3. The blimp launching and the skyscraper dedication were covered,
     in some order, by Lois and the reporter who was sent to Port Charles.
  4. South Amboy was not the site of either the beached whale or the
     skyscraper dedication.
  5. Bayonne is either the place that Corey went or the place where
     the whale was beached, or both.

  Determine: Reporter -- Location -- Story
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function breaking_news(print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    n = 4
    reporters = 1:n
    Corey, Jimmy, Lois, Perry = reporters

    @variable(model, 1 <= locations[1:n] <= n, Int)
    Bayonne, NewHope, PortCharles, SouthAmboy = locations

    @variable(model, 1 <= news[1:n] <= n, Int)
    Baby, Blimp, Skyscraper, Whale = news

    @constraint(model, locations in CS.AllDifferent())
    @constraint(model, news in CS.AllDifferent())

    # 1. The 30-pound baby wasn"t born in South Amboy or New Hope.
    @constraint(model, Baby != SouthAmboy)
    @constraint(model, Baby != NewHope)

    # 2. Jimmy didn"t go to Port Charles.
    @constraint(model, Jimmy != PortCharles)

    # 3. The blimp launching and the skyscraper dedication were covered,
    #    in some order, by Lois and the reporter who was sent to
    #    Port Charles.
    #=
    (
        (Blimp == Lois && Skyscraper == PortCharles)
        ||
        (Skyscraper == Lois && Blimp == PortCharles)
    )
    =#
    @constraint(model, Lois != PortCharles)
    @constraint(model, (Blimp == Lois && Skyscraper == PortCharles) || (Skyscraper == Lois && Blimp == PortCharles))

    # 4. South Amboy was not the site of either the beached whale or the
    #    skyscraper dedication.
    @constraint(model, SouthAmboy != Whale)
    @constraint(model, SouthAmboy != Skyscraper)

    # 5. Bayonne is either the place that Corey went or the place where
    #    the whale was beached, or both.
    c5 = @variable(model, [1:2],Bin)
    #=
    (
        Bayonne == Corey || Bayonne == Whale
    )
    =#
    @constraint(model, c5[1] := {Bayonne == Corey})
    @constraint(model, c5[2] := {Bayonne == Whale})
    @constraint(model, sum(c5) == 1)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                locations_val = convert.(Integer,JuMP.value.(locations; result=sol))
                news_val = convert.(Integer,JuMP.value.(news; result=sol))
                println("locations:$locations_val")
                println("news:$news_val")
                println()

            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

@time breaking_news()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module Issue980

import MathOptInterface
const MOI = MathOptInterface

abstract type StaticArray end
(::Type{SA})(x...) where {SA<:StaticArray} = SA(x)

function crash()
    model = MOI.Utilities.Model{Float64}()
    bridged = MOI.Bridges.LazyBridgeOptimizer(MOI.Utilities.Model{Float64}())
    return MOI.copy_to(bridged, model)
end

crash()

end

ext: jl

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
# Description:
# Dimensions: 2
#
# The Cross-in-Tray function has multiple global minima. It is shown here with a smaller domain in the second plot, so that its characteristic "cross" will be visible.
#
# Input Domain:
# The function is usually evaluated on the square xi ∈ [-10, 10], for all i = 1, 2.
#
# Global Minima:
# https://www.sfu.ca/~ssurjano/crossit3a.png

function crossit(xx)

    x1 = xx[1];
    x2 = xx[2];

    fact1 = sin(x1) * sin(x2);
    fact2 = exp(abs(100 - sqrt(x1^2 + x2^2) / pi));

    y = -0.0001 * (abs(fact1 * fact2) + 1)^0.1;

    return [y]
end


text:
function update_omega!(s::StateFS, c::ConstantsFS, d::DataFS, t::TunersFS)
  Q = length(s.omega)
  for q in 1:Q
    update_omega!(q, s, c, d, t)
  end
end

function update_omega!(q::Integer, s::StateFS, c::ConstantsFS, d::DataFS,
                       t::TunersFS)

  function loglike(om_q::Float64)::Float64
    om = copy(s.omega)  # make a copy
    om[q] = om_q
    p = compute_p(d.X, om)

    ll = 0.
    for i in 1:d.data.I
      lpdf(rik::Bool)::Float64 = logpdf(Bernoulli(p[i]), rik)
      for k in 1:c.constants.K
        ll += lpdf(s.r[i, k])
      end
    end

    return ll
  end

  logprior(om_q::Float64)::Float64 = logpdf(c.omega_prior, om_q)
  log_prob(om_q::Float64)::Float64 = loglike(om_q) + logprior(om_q)

  s.omega[q] = MCMC.metropolisAdaptive(s.omega[q], log_prob, t.omega[q])
end

# Omega = p when using Cell means model for X.
function update_omega_cellmeans!(s::StateFS, c::ConstantsFS, d::DataFS,
                                 t::TunersFS)
  Q = length(s.omega)
  for q in 1:Q
    update_omega_cellmeans!(q, s, c, d, t)
  end
end

function update_omega_cellmeans!(q::Integer, s::StateFS, c::ConstantsFS,
                                 d::DataFS, t::TunersFS)
  a_new, b_new = params(c.p_prior)

  xi = findfirst(z -> z == 1, d.X[i, :])
  K = c.constants.K
  for i in 1:d.data.I
    if xi == q
      ri_sum = sum(s.r[i, :])
      a_new += ri_sum
      b_new += (K - ri_sum)
    end
  end

  s.omega[q] = rand(Beta(a_new, b_new))
end


text:
import Base: AbstractSet, isempty, length, in, issubset, start, next, done,
             push!, ==, copy, intersect!, intersect, union!, union

struct CustomSet{T} <: AbstractSet{T}
    elements::Array{T, 1}
end

isempty(s::CustomSet) = isempty(s.elements)
length(s::CustomSet) = length(s.elements)
in(element, s::CustomSet) = in(element, s.elements)  # this also defines issubset(::CustomSet, ::CustomSet)
==(s1::CustomSet, s2::CustomSet) = issubset(s1, s2) && issubset(s2, s1)
copy(s::CustomSet) = CustomSet(copy(s.elements))
push!(s::CustomSet, element) = push!(s.elements, element)

# Iterator protocol
start(::CustomSet) = 1
next(s::CustomSet, state) = s.elements[state], state + 1
done(s::CustomSet, state) = state > length(s)

function disjoint(s1::CustomSet, s2::CustomSet)
    for element in s1
        element in s2 && return false
    end
    return true
end

function complement!(s1::CustomSet, s2::CustomSet)
    for (i, element) in enumerate(s1)
        element in s2 && deleteat!(s1.elements, i)
    end
    return s1
end
complement(s1::CustomSet, s2::CustomSet) = complement!(copy(s1), s2)

function intersect!(s1::CustomSet, s2::CustomSet)
    tbd = Int[]
    for (i, element) in enumerate(s1)
        element in s2 || push!(tbd, i)
    end
    deleteat!(s1.elements, tbd)
    return s1
end
intersect(s1::CustomSet, s2::CustomSet) = intersect!(copy(s1), s2)

function union!(s1::CustomSet, s2::CustomSet)
    for element in s2
        element in s1 || push!(s1, element)
    end
    return s1
end
union(s1::CustomSet, s2::CustomSet) = union!(copy(s1), s2)


text:
# The number of components in the color prior.
const D = 2

# The number of types of celestial objects (here, stars and galaxies).
const Ia = 2

type CatalogEntry
    pos::Vector{Float64}
    is_star::Bool
    star_fluxes::Vector{Float64}
    gal_fluxes::Vector{Float64}
    gal_frac_dev::Float64
    gal_ab::Float64
    gal_angle::Float64
    gal_scale::Float64
    objid::String
    thing_id::Int
end

"""
Parameters of a single normal component of a galaxy.

Attributes:
  etaBar: The weight of the galaxy component
  nuBar: The scale of the galaxy component
"""
immutable GalaxyComponent
    etaBar::Float64
    nuBar::Float64
end

typealias GalaxyPrototype Vector{GalaxyComponent}

"""
Pre-defined shapes for galaxies.

Returns:
  dev_prototype: An array of GalaxyComponent for de Vaucouleurs galaxy types
  exp_prototype: An array of GalaxyComponent for exponenttial galaxy types
"""
function get_galaxy_prototypes()
    dev_amp = [
        4.26347652e-2, 2.40127183e-1, 6.85907632e-1, 1.51937350,
        2.83627243, 4.46467501, 5.72440830, 5.60989349]
    dev_amp /= sum(dev_amp)
    dev_var = [
        2.23759216e-4, 1.00220099e-3, 4.18731126e-3, 1.69432589e-2,
        6.84850479e-2, 2.87207080e-1, 1.33320254, 8.40215071]

    exp_amp = [
        2.34853813e-3, 3.07995260e-2, 2.23364214e-1,
        1.17949102, 4.33873750, 5.99820770]
    exp_amp /= sum(exp_amp)
    exp_var = [
        1.20078965e-3, 8.84526493e-3, 3.91463084e-2,
        1.39976817e-1, 4.60962500e-1, 1.50159566]

    # Adjustments to the effective radius hard-coded above.
    # (The effective radius is the distance from the center containing half
    # the light.)
    effective_radii = [1.078031, 0.928896]
    dev_var /= effective_radii[1]^2
    exp_var /= effective_radii[2]^2

    exp_prototype = [GalaxyComponent(exp_amp[j], exp_var[j]) for j in 1:6]
    dev_prototype = [GalaxyComponent(dev_amp[j], dev_var[j]) for j in 1:8]
    (dev_prototype, exp_prototype)
end

const galaxy_prototypes = get_galaxy_prototypes()

immutable PriorParams
    a::Vector{Float64}  # formerly Phi
    r_mean::Vector{Float64}
    r_var::Vector{Float64}
    k::Matrix{Float64}  # formerly Xi
    c_mean::Array{Float64, 3} # formerly Omega
    c_cov::Array{Float64, 4} # formerly Lambda
end

function load_prior()
    # set a = [.99, .01] if stars are underrepresented
    # due to the greater flexibility of the galaxy model
    #a = [0.28, 0.72]
    a = [0.99, 0.01]
    r_mean = Array(Float64, Ia)
    r_var = Array(Float64, Ia)
    k = Array(Float64, D, Ia)
    c_mean = Array(Float64, B - 1, D, Ia)
    c_cov = Array(Float64, B - 1, B - 1, D, Ia)

    v05 = VERSION >= v"0.5.0-dev" ? "-v05" : ""
    stars_file = open(joinpath(cfgdir, "stars$D$v05.dat"))
    r_fit1, k[:, 1], c_mean[:,:,1], c_cov[:,:,:,1] = deserialize(stars_file)
    close(stars_file)

    gals_file = open(joinpath(cfgdir, "gals$D$v05.dat"))
    r_fit2, k[:, 2], c_mean[:,:,2], c_cov[:,:,:,2] = deserialize(gals_file)
    close(gals_file)

    # These "magic numbers" have been in use for a while.
    # They were initially gamma parameters, and now they are log normal
    # parameters.  TODO: Get rid of these and use an empirical prior.
    # r = [0.47 1.28; 1/0.012 1/0.11] # These were gamma (shape, scale)

    mean_brightness = [0.47 / 0.012, 1.28 / 0.11 ]
    var_brightness = [0.47 / (0.012 ^ 2), 1.28 / (0.11 ^ 2) ]

    # The prior contains parameters of a lognormal distribution with
    # the desired means.
    r_var = log.(var_brightness ./ (mean_brightness .^ 2) + 1)
    r_mean = log.(mean_brightness) - 0.5 * r_var
    PriorParams(a, r_mean, r_var, k, c_mean, c_cov)
end

const prior = load_prior()

"""
Return a default-initialized VariationalParams object.
"""
function init_source(init_pos::Vector{Float64})
    ret = Array(Float64, length(CanonicalParams))
    ret[ids.a[2, 1]] = 0.5
    ret[ids.a[1, 1]] = 1.0 - ret[ids.a[2, 1]]
    ret[ids.u[1]] = init_pos[1]
    ret[ids.u[2]] = init_pos[2]
    ret[ids.r1] = log(2.0)
    ret[ids.r2] = 1e-3
    ret[ids.e_dev] = 0.5
    ret[ids.e_axis] = 0.5
    ret[ids.e_angle] = 0.
    ret[ids.e_scale] = 1.
    ret[ids.k] = 1. / size(ids.k, 1)
    ret[ids.c1] = 0.
    ret[ids.c2] =  1e-2
    ret
end

"""
Return a VariationalParams object initialized form a catalog entry.
"""
function init_source(ce::CatalogEntry)
    # TODO: sync this up with the transform bounds
    ret = init_source(ce.pos)

    ret[ids.a[1, 1]] = ce.is_star ? 0.8: 0.2
    ret[ids.a[2, 1]] = ce.is_star ? 0.2: 0.8

    ret[ids.r1[1]] = log(max(0.1, ce.star_fluxes[3]))
    ret[ids.r1[2]] = log(max(0.1, ce.gal_fluxes[3]))

    function get_color(c2, c1)
        c2 > 0 && c1 > 0 ? min(max(log(c2 / c1), -9.), 9.) :
            c2 > 0 && c1 <= 0 ? 3.0 :
                c2 <= 0 && c1 > 0 ? -3.0 : 0.0
    end

    function get_colors(raw_fluxes)
        [get_color(raw_fluxes[c+1], raw_fluxes[c]) for c in 1:4]
    end

    ret[ids.c1[:, 1]] = get_colors(ce.star_fluxes)
    ret[ids.c1[:, 2]] = get_colors(ce.gal_fluxes)

    ret[ids.e_dev] = min(max(ce.gal_frac_dev, 0.015), 0.985)

    ret[ids.e_axis] = ce.is_star ? .8 : min(max(ce.gal_ab, 0.015), 0.985)
    ret[ids.e_angle] = ce.gal_angle
    ret[ids.e_scale] = ce.is_star ? 0.2 : max(ce.gal_scale, 0.2)

    ret
end


text:
using  MCMCChains, Test

# https://github.com/TuringLang/AdvancedMH.jl/pull/63
@testset "order of parameters" begin
    chains = Chains(rand(2, 10, 4), vcat(["μ[$i]" for i in 1:9], :lp), (internals=[:lp],))
    @test names(chains, :parameters) == [Symbol("μ[$i]") for i in 1:9]
end

@testset "describe sections" begin

    a3d = rand(500, 9, 4)

    cnames = ["lp__"  , "accept_stat__", "stepsize__" , "treedepth__" ,
        "n_leapfrog__" , "divergent__"  , "energy__", "sigma", "mu" ]

    pi = filter(p -> length(p) > 2 && p[end-1:end] == "__", cnames)
    p = filter(p -> !(p in  pi), cnames)

    c = Chains(a3d,
        cnames,
            Dict(
                :parameters => p,
                :internals => pi
                )
        )

    c2 = Chains(a3d,
        cnames,
            Dict(
                :parameters => p,
                :internals => pi
            )
        )

    tmpdir = tempdir()
    open(joinpath(tmpdir, "sections_test1.txt"), "w") do io
        describe(io, c, sections=:parameters);
        describe(io, c, sections=:internals);
        describe(io, c);
    end

    open(joinpath(tmpdir, "sections_test2.txt"), "w") do io
        describe(io, c2, sections=:parameters);
        describe(io, c2, sections=:internals);
        describe(io, c2);
    end

    @test open(f->read(f, String), joinpath(tmpdir, "sections_test1.txt")) ==
        open(f->read(f, String), joinpath(tmpdir, "sections_test2.txt"))
end


text:
# Probabilistic Principal Component Analysis

"""Probabilistic PCA type"""
struct PPCA{T<:Real}
    mean::Vector{T}       # sample mean: of length d (mean can be empty, which indicates zero mean)
    W::Matrix{T}          # weight matrix: of size d x p
    σ²::T                 # residual variance
end

## properties

indim(M::PPCA) = size(M.W, 1)
outdim(M::PPCA) = size(M.W, 2)

mean(M::PPCA) = fullmean(indim(M), M.mean)
projection(M::PPCA) = svd(M.W).U # recover principle components from the weight matrix
var(M::PPCA) = M.σ²
loadings(M::PPCA) = M.W

## use

function transform(m::PPCA{T}, x::AbstractVecOrMat{T}) where {T<:Real}
    xn = centralize(x, m.mean)
    W  = m.W
    n = outdim(m)
    M = W'W .+ m.σ² * Matrix{T}(I, n, n)
    return inv(M)*m.W'*xn
end

function reconstruct(m::PPCA{T}, z::AbstractVecOrMat{T}) where {T<:Real}
    W  = m.W
    WTW = W'W
    n = outdim(m)
    M  = WTW .+ var(m) * Matrix{T}(I, n, n)
    return W*inv(WTW)*M*z .+ mean(m)
end

## show

function Base.show(io::IO, M::PPCA)
    print(io, "Probabilistic PCA(indim = $(indim(M)), outdim = $(outdim(M)), σ² = $(var(M)))")
end

## core algorithms

function ppcaml(Z::AbstractMatrix{T}, mean::Vector{T};
                tol::Real=1.0e-6, # convergence tolerance
                maxoutdim::Int=size(Z,1)-1) where {T<:Real}

    check_pcaparams(size(Z,1), mean, maxoutdim, 1.)

    d = size(Z,1)

    # SVD decomposition
    Svd = svd(Z)
    λ = Svd.S
    ord = sortperm(λ; rev=true)
    V = λ[ord]

    # filter 0 eigenvalues and adjust number of latent dimensions
    idxs = findall(V .< tol)
    l = length(idxs)
    l = l == 0 ? maxoutdim : l

    # variance "loss" in the projection
    σ² = sum(V[l+1:end])/(d-l)

    @inbounds for i in 1:l
        V[i] = sqrt(V[i] - σ²)
    end

    W = Svd.U[:,ord[1:l]]*diagm(0 => V[1:l])

    return PPCA(mean, W, σ²)
end

function ppcaem(S::AbstractMatrix{T}, mean::Vector{T}, n::Int;
                maxoutdim::Int=size(S,1)-1,
                tol::Real=1.0e-6,   # convergence tolerance
                maxiter::Integer=1000) where {T<:Real}

    check_pcaparams(size(S,1), mean, maxoutdim, 1.)

    d = size(S,1)
    q = maxoutdim
    Iq = Matrix{T}(I, q, q)
    Id = Matrix{T}(I, d, d)
    W = Matrix{T}(I, d, q)
    σ² = zero(T)
    M⁻¹ = inv(W'W .+ σ² * Iq)

    i = 1
    L_old = 0.
    chg = NaN
    converged = false
    while i < maxiter
        # EM-steps
        SW = S*W
        W⁺  = SW*inv(σ²*Iq + M⁻¹*W'*SW)
        σ²⁺ = tr(S - SW*M⁻¹*(W⁺)')/d
        # new parameters
        W = W⁺
        σ² = σ²⁺
        # log likelihood
        C = W*W'.+ σ²*Id
        M⁻¹ = inv(W'*W .+ σ²*Iq)
        C⁻¹ = (Id - W*M⁻¹*W')/σ²
        L = (-n/2)*(log(det(C)) + tr(C⁻¹*S))  # (-n/2)*d*log(2π) omitted

        @debug "Likelihood" iter=i L=L ΔL=abs(L_old - L)
        chg = abs(L_old - L)
        if chg < tol
            converged = true
            break
        end
        L_old = L
        i += 1
    end
    converged || throw(ConvergenceException(maxiter, chg, oftype(chg, tol)))

    return PPCA(mean, W, σ²)
end

function bayespca(S::AbstractMatrix{T}, mean::Vector{T}, n::Int;
                 maxoutdim::Int=size(S,1)-1,
                 tol::Real=1.0e-6,   # convergence tolerance
                 maxiter::Integer=1000) where {T<:Real}

    check_pcaparams(size(S,1), mean, maxoutdim, 1.)

    d = size(S,1)
    q = maxoutdim
    Iq = Matrix{T}(I, q, q)
    Id = Matrix{T}(I, d, d)
    W = Matrix{T}(I, d, q)
    wnorm = zeros(T, q)
    σ² = zero(T)
    M = W'*W .+ σ²*Iq
    M⁻¹ = inv(M)
    α = zeros(T, q)

    i = 1
    chg = NaN
    L_old = 0.
    converged = false
    while i < maxiter
        # EM-steps
        SW = S*W
        W⁺  = SW*inv(σ²*(Iq+diagm(0=>α)*M/n) + M⁻¹*W'*SW)
        σ²⁺ = tr(S - SW*M⁻¹*(W⁺)')/d
        # new parameters
        W = W⁺
        σ² = σ²⁺
        @inbounds for j in 1:q
            wnorm[j] = norm(W[:,j])^2
            α[j] = wnorm[j] < eps() ? maxintfloat(Float64) : d/wnorm[j]
        end
        # log likelihood
        C = W*W'.+ σ²*Id
        M = W'*W .+ σ²*Iq
        M⁻¹ = inv(M)
        C⁻¹ = (Id - W*M⁻¹*W')/σ²
        L = (-n/2)*(log(det(C)) + tr(C⁻¹*S))  # (-n/2)*d*log(2π) omitted

        @debug "Likelihood" iter=i L=L ΔL=abs(L_old - L)
        chg = abs(L_old - L)
        if chg < tol
            converged = true
            break
        end
        L_old = L
        i += 1
    end
    converged || throw(ConvergenceException(maxiter, chg, oftype(chg, tol)))

    return PPCA(mean, W[:,wnorm .> 0.], σ²)
end

## interface functions

function fit(::Type{PPCA}, X::AbstractMatrix{T};
             method::Symbol=:ml,
             maxoutdim::Int=size(X,1)-1,
             mean=nothing,
             tol::Real=1.0e-6,   # convergence tolerance
             maxiter::Integer=1000) where {T<:Real}

    @assert !SparseArrays.issparse(X) "Use Kernel PCA for sparce arrays"

    d, n = size(X)

    # process mean
    mv = preprocess_mean(X, mean)
    if !(isempty(mv) || length(mv) == d)
        throw(DimensionMismatch("Dimensions of weight matrix and mean are inconsistent."))
    end

    if method == :ml
        Z = centralize(X, mv)
        M = ppcaml(Z, mv, maxoutdim=maxoutdim, tol=tol)
    elseif method == :em || method == :bayes
        S = covm(X, isempty(mv) ? 0 : mv, 2)
        if method == :em
            M = ppcaem(S, mv, n, maxoutdim=maxoutdim, tol=tol, maxiter=maxiter)
        elseif method == :bayes
            M = bayespca(S, mv, n, maxoutdim=maxoutdim, tol=tol, maxiter=maxiter)
        end
    else
        throw(ArgumentError("Invalid method name $(method)"))
    end

    return M::PPCA
end


text:
mutable struct GebOrganism
    key::UInt64

    genotype::String
    rules::Vector{Rule}

    coordinates::Tuple{Float64, Float64}
    direction::UInt16

    network::Network

    parent_genotypes::Tuple{String, String}

    time_birth::UInt64
    age::UInt64

    daughters::Vector{GebOrganism}

    function GebOrganism(genotype, coordinates)
        GebOrganism(0, genotype, coordinates)
    end
    function GebOrganism(key, genotype, coordinates)
        GebOrganism(key, genotype, coordinates, "", "", 0)
    end
    function GebOrganism(genotype, coordinates, parent_1_genotype, parent_2_genotype, time)
        GebOrganism(0, genotype, coordinates, parent_1_genotype, parent_2_genotype, time)
    end
    function GebOrganism(key, genotype, coordinates, parent_1_genotype, parent_2_genotype, time)
        rules = get_filtered_rules(genotype)
        direction = rand(0:359)
        network = get_axiom_network()

        new(key, genotype, rules, coordinates, direction, network,
            (parent_1_genotype, parent_2_genotype), time, 0, [])
    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
#=
Dash port of Shiny telephones by region example:

https://shiny.rstudio.com/gallery/telephones-by-region.html
=#

using Dash, DashBootstrapComponents
using CSV, DataFrames, HTTP

data_path = "https://cdn.opensource.faculty.ai/world-phones/data.csv";
data = DataFrame(CSV.File(HTTP.get(data_path).body));

app = dash(external_stylesheets = [dbc_themes.BOOTSTRAP]);

controls = dbc_card(
    [
        dbc_formgroup([
            dbc_label("Region"),
            dcc_dropdown(
                id = "region-selector",
                options = [
                    Dict("label" => col, "value" => col) for
                    col in names(data) if col != "Year"
                ],
                value = "S.Amer",
            ),
        ]),
        html_hr(),
        html_p("Data from AT&T (1961) The World's Telephones.", className = "text-muted"),
    ],
    body = true,
)

app.layout = dbc_container(
    [
        html_h1("Telephones by region"),
        html_hr(),
        dbc_row(
            [dbc_col(controls, md = 4), dbc_col(dcc_graph(id = "phones-graph"), md = 8)],
            align = "center",
        ),
    ],
    fluid = true,
);

callback!(
    app,
    Output("phones-graph", "figure"),
    Input("region-selector", "value"),
) do region
    fig_data = [(y = data[:, region], type = "bar")]
    fig_layout = (
        xaxis = (tickvals = [i for i = 0:6], ticktext = data[:, "Year"], title = "Year"),
        yaxis = (title = "Number of Telephones",),
    )

    return (data = fig_data, layout = fig_layout)

end;

run_server(app, "0.0.0.0", 8888, debug = true);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
# Make c library
run(`make -C v4lcapture/`)

# set linker path to library
libv4lpath = joinpath(dirname(@__FILE__),"v4lcapture")

#create a file to push to the DL_LOAD_PATH when module loads,
f = open("loadpath.jl", "w")

write(f,"# This is an automatically generated file, do not edit.\n")
write(f,"using Libdl\n")
write(f,"push!(Libdl.DL_LOAD_PATH, \"$(libv4lpath)\")")

close(f)


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
config = MOIT.TestConfig()

bridged_mock = MOIB.Variable.Vectorize{Float64}(mock)

@testset "get scalar constraint" begin
    x, cx = MOI.add_constrained_variable(bridged_mock, MOI.GreaterThan(1.0))
    fx = MOI.SingleVariable(x)
    func = 2.0 * fx
    set = MOI.GreaterThan(5.0)
    err = MOI.ScalarFunctionConstantNotZero{
        Float64, typeof(func), typeof(set)}(1.0)
    @test_throws err MOI.add_constraint(bridged_mock, func + 1.0, set)

    c = MOI.add_constraint(bridged_mock, func, set)
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), c) ≈ func
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), c) == set
    MOI.set(bridged_mock, MOI.ConstraintName(), c, "c")

    @testset "Mock model" begin
        MOI.set(mock, MOI.VariableName(),
                MOI.get(mock, MOI.ListOfVariableIndices()), ["y"])
        MOI.set(mock, MOI.ConstraintName(),
                MOI.get(mock, MOI.ListOfConstraintIndices{
                    MOI.VectorOfVariables, MOI.Nonnegatives}()),
                ["cy"])
        s = """
        variables: y
        cy: [y] in MathOptInterface.Nonnegatives(1)
        c: 2.0y >= 3.0
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(mock, model, ["y"], ["cy", "c"])
    end
    @testset "Bridged model" begin
        MOI.set(bridged_mock, MOI.VariableName(), x, "x")
        MOI.set(bridged_mock, MOI.ConstraintName(), cx, "cx")
        s = """
        variables: x
        cx: x >= 1.0
        c: 2.0x >= 5.0
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(bridged_mock, model, ["x"], ["cx", "c"])
    end
end

@testset "exp3 with add_constrained_variable for `y`" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [log(5), 0.0],
                          (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [0.0],
                          (MOI.VectorAffineFunction{Float64}, MOI.ExponentialCone)   => [[-1.0, log(5)-1, 1/5]])

    MOI.empty!(bridged_mock)
    x = MOI.add_variable(bridged_mock)
    @test MOI.get(bridged_mock, MOI.NumberOfVariables()) == 1
    fx = MOI.SingleVariable(x)
    xc = MOI.add_constraint(bridged_mock, 2.0fx, MOI.LessThan(4.0))
    y, yc = MOI.add_constrained_variable(bridged_mock, MOI.LessThan(5.0))
    @test yc.value == y.value == -1
    @test MOI.get(bridged_mock, MOI.NumberOfVariables()) == 2
    @test length(MOI.get(bridged_mock, MOI.ListOfVariableIndices())) == 2
    @test Set(MOI.get(bridged_mock, MOI.ListOfVariableIndices())) == Set([x, y])
    fy = MOI.SingleVariable(y)
    ec = MOI.add_constraint(bridged_mock,
                            MOIU.operate(vcat, Float64, fx, 1.0, fy),
                            MOI.ExponentialCone())

    MOI.optimize!(bridged_mock)
    @test MOI.get(bridged_mock, MOI.VariablePrimal(), x) ≈ log(5)
    @test MOI.get(bridged_mock, MOI.VariablePrimal(), y) ≈ 5.0
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), xc) ≈ 2log(5)
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), yc) ≈ 5
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), ec) ≈ [log(5), 1., 5.0]
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), xc) ≈ 0.0
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), yc) ≈ -1/5
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), ec) ≈ [-1., log(5)-1, 1/5]

    err = ErrorException(
        "Cannot add two `SingleVariable`-in-`MathOptInterface.LessThan{Float64}`" *
        " on the same variable MathOptInterface.VariableIndex(-1)."
    )
    @test_throws err MOI.add_constraint(bridged_mock, MOI.SingleVariable(y), MOI.LessThan(4.0))

    cis = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{Float64}, MOI.ExponentialCone}())
    @test length(cis) == 1

    @testset "get `UnknownVariableAttribute``" begin
        err = ArgumentError(
            "Variable bridge of type `MathOptInterface.Bridges.Variable.VectorizeBridge{Float64,MathOptInterface.Nonpositives}`" *
            " does not support accessing the attribute `MathOptInterface.Test.UnknownVariableAttribute()`."
        )
        @test_throws err MOI.get(bridged_mock, MOIT.UnknownVariableAttribute(), y)
    end

    @testset "set `ConstraintSet`" begin
        ci = MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}(y.value)
        attr = MOI.ConstraintSet()
        err = MOI.SetAttributeNotAllowed(attr,
            "The variable `MathOptInterface.VariableIndex(12345676)` is bridged by the `VectorizeBridge`.")
        @test_throws err MOI.set(bridged_mock, attr, ci, MOI.LessThan(4.0))
    end

    @testset "MultirowChange" begin
        change = MOI.MultirowChange(y, [(3, 0.0)])
        message = "The change $change" *
            " contains variables bridged into a function with nonzero constant."
        err = MOI.ModifyConstraintNotAllowed(cis[1], change, message)
        @test_throws err MOI.modify(bridged_mock, cis[1], change)
    end

    @testset "ScalarCoefficientChange" begin
        change = MOI.ScalarCoefficientChange(y, 0.0)
        attr = MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()
        message = "The change MathOptInterface.ScalarCoefficientChange{Float64}(MathOptInterface.VariableIndex(-1), 0.0)" *
            " contains variables bridged into a function with nonzero constant."
        err = MOI.ModifyObjectiveNotAllowed(change, message)
        @test_throws err MOI.modify(bridged_mock, attr, change)
    end

    MOI.set(bridged_mock, MOI.VariableName(), x, "x")
    MOI.set(bridged_mock, MOI.ConstraintName(), xc, "xc")
    MOI.set(bridged_mock, MOI.ConstraintName(), ec, "ec")
    z = MOI.get(mock, MOI.ListOfVariableIndices())[2]
    @testset "Mock model" begin
        MOI.set(mock, MOI.VariableName(), z, "z")
        MOI.set(mock, MOI.ConstraintName(),
                MOI.get(mock, MOI.ListOfConstraintIndices{
                    MOI.VectorOfVariables, MOI.Nonpositives}()),
                ["zc"])
        s = """
        variables: x, z
        zc: [z] in MathOptInterface.Nonpositives(1)
        xc: 2.0x <= 4.0
        ec: [x, 1.0, z + 5.0] in MathOptInterface.ExponentialCone()
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(mock, model, ["x", "z"], ["zc", "xc", "ec"])
    end
    @testset "Bridged model" begin
        MOI.set(bridged_mock, MOI.VariableName(), y, "y")
        MOI.set(bridged_mock, MOI.ConstraintName(), yc, "yc")
        s = """
        variables: x, y
        yc: y <= 5.0
        xc: 2.0x <= 4.0
        ec: [x, 1.0, y] in MathOptInterface.ExponentialCone()
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(bridged_mock, model, ["x", "y"], ["yc", "xc", "ec"])
    end

    @test MOI.supports(bridged_mock, MOI.VariablePrimalStart(), MOI.VariableIndex)
    @test MOI.supports(bridged_mock, MOI.VariablePrimalStart(), typeof(MOIB.bridge(bridged_mock, y)))
    MOI.set(bridged_mock, MOI.VariablePrimalStart(), y, 1.0)
    @test MOI.get(mock, MOI.VariablePrimalStart(), z) == -4
    @test MOI.get(bridged_mock, MOI.VariablePrimalStart(), y) == 1

    test_delete_bridged_variable(bridged_mock, y, MOI.LessThan{Float64}, 2, (
        (MOI.VectorOfVariables, MOI.Nonpositives, 0),
    ))
end


text:
# Test stacking of SAC traces
using Test
using SAC

normal_dist(x, μ, σ) = exp(-(x - μ)/2σ^2)/sqrt(2σ^2*π)

@testset "Stacking" begin
    # Create five traces which are just a normal distribution and stack the traces
    # with shifts of the means, meaning the stack should be just a normal distribution
    # itself.
    let N = 5, μ = 1:5, σ = 0.1, delta = 0.01, b = 0.0, e = 10,
            npts = round(Int, (e - b)/delta) + 1

        local S = Array{SACtr}(undef, N)
        for i in 1:N
            S[i] = SACtr(delta, npts, b)
            S[i].t[:] = normal_dist.(time(S[i]), μ[i], σ)
            S[i].a = μ[i]
            S[i].ka = "mu"
            SAC.update_headers!(S[i])
        end

        # Create the stack
        local t1 = -0.5
        local t2 = 0.5
        local st = stack(S, (-0.5,0.5), collect(μ), ones(N))
        @test st.t ≈ normal_dist.(t1:delta:t2, 0, σ)
    end
end


text:
export SplitMergeMove

struct SplitMergeMove <: InvMcmcMove
    ν::Int
    ind_split::Vector{Int}
    ind_merge::Vector{Int}
    counts::Vector{Int}
    function SplitMergeMove(
        ;ν::Int=1
        )
        ind_split = zeros(Int, ν)
        ind_merge = zeros(Int, 2ν)
        new(ν, ind_split, ind_merge, [0,0])
    end
end

Base.show(io::IO, x::SplitMergeMove) = print(io, "SplitMergeMove(ν=$(x.ν))")

function multi_split_prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    move::SplitMergeMove,
    pointers::InteractionSequence{Int}
    )

    ν = move.ν
    N = length(S_curr)
    ε = rand(1:min(N, ν))
    ind_split = view(move.ind_split, 1:ε)
    ind_loc = view(move.ind_merge, 1:2ε)

    store = Vector{Int}[]

    log_ratio = 0.0

    StatsBase.seqsample_a!(1:N, ind_split)
    StatsBase.seqsample_a!(1:(N+ε), ind_loc)

    for i in Iterators.reverse(ind_split)
        @assert length(S_prop[i])>1 "Cannot split path of length one."
        tmp1 = popat!(S_prop, i)
        n = length(tmp1)
        k = rand(1:(n-1)) # Split location
        log_ratio += log(n-1)
        tmp2 = popfirst!(pointers)
        copy!(tmp2, tmp1[(k+1):end])
        resize!(tmp1, k)
        pushfirst!(store, tmp2)
        pushfirst!(store, tmp1)
    end

    for i in ind_loc
        tmp = popfirst!(store)
        insert!(S_prop, i, tmp)
    end

    log_ratio += log(min(N, ν)) - log(min(floor(length(S_prop/2)), ν))

    return log_ratio

end

function multi_merge_prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    move::SplitMergeMove,
    pointers::InteractionSequence{Int}
    )

    ν = move.ν
    N = length(S_curr)
    ε = rand(1:min(floor(Int, N/2), ν))

    ind_merge = view(move.ind_merge, 1:2ε)
    ind_loc = view(move.ind_split, 1:ε)

    store = Vector{Int}[]

    log_ratio = 0.0

    StatsBase.seqsample_a!(1:N, ind_merge)
    StatsBase.seqsample_a!(1:(N-ε), ind_loc)

    i = length(ind_merge)
    while i > 0
        # Merge i-1 and i
        tmp2 = popat!(S_prop, ind_merge[i])
        tmp1 = popat!(S_prop, ind_merge[i-1])
        append!(tmp1, tmp2)
        pushfirst!(pointers, tmp2)
        pushfirst!(store, tmp1)
        log_ratio -= log(length(tmp1)-1)
        i -= 2
    end

    for i in ind_loc
        tmp = popfirst!(store)
        insert!(S_prop, i, tmp)
    end

    log_ratio += log(min(floor(Int, N/2), ν)) - log(min(length(S_prop), ν))

    return log_ratio
end

# function log_binomial_ratio_split_merge(N, τ, ε)
#     # \binom{N-τ}{ε}/\binom{N}{ε}
#     out = 0.0
#     if τ < ε
#         for j in 0:(τ-1)
#             out += log(N-ε-j) - log(N-j)
#         end
#     else
#         for j in 0:(ε-1)
#             out += log(N-τ-j) - log(N-j)
#         end
#     end
#     return out
# end

# Special case (when we have paths of length one present)
function special_multi_split_prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    move::SplitMergeMove,
    pointers::InteractionSequence{Int}
    )

    ν = move.ν
    N = length(S_curr)
    ind_ok_len = findall(x->length(x)>1, S_curr) # Good paths (length > 1)
    τ = N-length(ind_ok_len)                     # Number of bad paths
    ε = rand(1:min(N-τ, ν))
    ind_split = view(move.ind_split, 1:ε)
    ind_loc = view(move.ind_merge, 1:2ε)

    store = Vector{Int}[]

    log_ratio = 0.0

    StatsBase.seqsample_a!(ind_ok_len, ind_split)  # Sample subsequence of good paths
    StatsBase.seqsample_a!(1:(N+ε), ind_loc)

    for i in Iterators.reverse(ind_split)
        @assert length(S_prop[i])>1 "Cannot split path of length one."
        tmp1 = popat!(S_prop, i)
        n = length(tmp1)
        k = rand(1:(n-1)) # Split location
        log_ratio += log(n-1)
        tmp2 = popfirst!(pointers)
        copy!(tmp2, tmp1[(k+1):end])
        resize!(tmp1, k)
        pushfirst!(store, tmp2)
        pushfirst!(store, tmp1)
    end

    for i in ind_loc
        tmp = popfirst!(store)
        insert!(S_prop, i, tmp)
    end

    log_ratio += log(min(N-τ, ν)) - log(min(floor(length(S_prop/2)), ν))

    log_ratio += log_binomial_ratio_split_merge(N, τ, ε)

    return log_ratio

end

function special_multi_merge_prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    move::SplitMergeMove,
    pointers::InteractionSequence{Int}
    )

    ν = move.ν
    N = length(S_curr)
    ε = rand(1:min(floor(Int, N/2), ν))

    ind_merge = view(move.ind_merge, 1:2ε)
    ind_loc = view(move.ind_split, 1:ε)

    store = Vector{Int}[]

    log_ratio = 0.0

    StatsBase.seqsample_a!(1:N, ind_merge)
    StatsBase.seqsample_a!(1:(N-ε), ind_loc)

    i = length(ind_merge)
    while i > 0
        # Merge i-1 and i
        tmp2 = popat!(S_prop, ind_merge[i])
        tmp1 = popat!(S_prop, ind_merge[i-1])
        append!(tmp1, tmp2)
        pushfirst!(pointers, tmp2)
        pushfirst!(store, tmp1)
        log_ratio -= log(length(tmp1)-1)
        i -= 2
    end

    for i in ind_loc
        tmp = popfirst!(store)
        insert!(S_prop, i, tmp)
    end

    M = length(S_prop)
    τ = sum(length(x)==1 for x in S_prop) # Find number of length one paths
    log_ratio -= log_binomial_ratio_split_merge(M, τ, ε)

    log_ratio += log(min(floor(Int, N/2), ν)) - log(min(M-τ, ν))

    return log_ratio
end

function prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    move::SplitMergeMove,
    pointers::InteractionSequence{Int},
    V::UnitRange
    )

    log_ratio = 0.0
    τ = sum(length(x)==1 for x in S_curr)
    N = length(S_curr)

    # If there is only one path of lenth one, we cannot do anything.
    # So just reject instantly
    if (N==τ==1)
        return 0, suff_stat_curr
    # If one path of length > 1 we split with probability one
    elseif N==1
        do_split = true
        log_ratio += multi_split_prop_sample!(
            S_curr, S_prop,
            move,
            pointers
        )
        # If S_prop has all length one paths, we would come back with probability one
        # hence the following...
        log_ratio += all(length(x)==1 for x in S_prop) ? 0.0 : log(0.5)
    # If no length one paths, use std proposal
    elseif τ==0
        do_split = rand(Bool)
        if do_split
        log_ratio += multi_split_prop_sample!(
            S_curr, S_prop,
            move,
            pointers
        )
        else
            log_ratio += multi_merge_prop_sample!(
                S_curr, S_prop,
                move,
                pointers
            )
        end
    # If all paths are length one, do merge with probability one
    elseif τ==N
        do_split = false
        log_ratio += special_multi_merge_prop_sample!(
            S_curr, S_prop,
            move,
            pointers
        )
        log_ratio += log(0.5) # Accounts for uneven choice between this and split move
    # Else to use special proposal function
    else
        do_split = rand(Bool)
        if do_split
            log_ratio += special_multi_split_prop_sample!(
                S_curr, S_prop,
                move,
                pointers
            )
        else
            log_ratio += special_multi_merge_prop_sample!(
                S_curr, S_prop,
                move,
                pointers
            )
        end
    end
    return log_ratio
end


text:
MLJModelInterface.save(::MLJTuning.EitherTunedModel, fitresult::Machine) =
    serializable(fitresult)

function MLJModelInterface.restore(::MLJTuning.EitherTunedModel, fitresult)
    fitresult.fitresult = restore(fitresult.model, fitresult.fitresult)
    return fitresult
end

max_stars_repo_name: alan-turing-institute/MLJTuning

max_issues_repo_name: alan-turing-institute/MLJTuning

max_forks_repo_name: alan-turing-institute/MLJTuning


text:
DATA_FILE = "02/data"

# This solution is going to be a lot longer than needed. I want to explore the
# use of structs and generic functions, and the relationship between "command"
# and "function" seemed perfect.

# I use two structs: Location1 and Location2. These are used in part 1 and
# part 2, respectively. Thanks to Julia's multiple dispatch, I don't need to
# redefine functions between the two parts.

## Shared code for part 1 and 2

VALID_MANOEUVRES = ["forward", "up", "down"]

commands = readlines(DATA_FILE)

# Location is the abstract supertype of Location1 and Location2
abstract type Location end
aoc_answer(loc::Location) = loc.horizontal * loc.depth

"""
    manoeuvre!(loc, command)

Interpret the command and use it to manoeuvre the given location. This
function dispatches based on `loc` through the given command. A command
such as "forward 1" is interpreted as `forward!(loc, 1)`.

# Examples
```julia-repl
julia> l = Location1() # initial position (0, 0)
julia> manoeuvre(l, "forward 1")
# new position: (1, 0)
```
"""
function manoeuvre!(loc::Location, command::String)
    command_split = split(command, " ")

    manoeuvre_direction = command_split[1]
    if manoeuvre_direction ∉ VALID_MANOEUVRES
        throw(ArgumentError(
            manoeuvre_direction *
            "is not a valid direction. Must be one of: " *
            join(VALID_COMMANDS, ", ")
        ))
    end
    manoeuvre_function = getfield(Main, Symbol(manoeuvre_direction * "!"))

    manoeuvre_amount = parse(Int64, command_split[2])

    manoeuvre_function(loc, manoeuvre_amount)
end

# Part 1

mutable struct Location1 <: Location
    horizontal::Int64
    depth::Int64
end
Location1() = Location1(0, 0)

function forward!(loc::Location1, x::Int64)
    loc.horizontal += x
end

function up!(loc::Location1, x::Int64)
    loc.depth -= x
    loc.depth = loc.depth < 0 ? 0 : loc.depth
end

function down!(loc::Location1, x::Int64)
    loc.depth += x
end

l1 = Location1()
for command in commands
    manoeuvre!(l1, command)
end
aoc_answer(l1)

## Part 2

mutable struct Location2 <: Location
    horizontal::Int64
    depth::Int64
    aim::Int64
end
Location2() = Location2(0, 0, 0)

function up!(loc::Location2, x::Int64)
    loc.aim -= x
end

function down!(loc::Location2, x::Int64)
    loc.aim += x
end

function forward!(loc::Location2, x::Int64)
    loc.horizontal += x
    loc.depth += loc.aim * x
    loc.depth = loc.depth < 0 ? 0 : loc.depth
end

l2 = Location2()
for command in commands
    manoeuvre!(l2, command)
end
aoc_answer(l2)


text:
"A circular or spherical specie of particle with homogenious material properties."
mutable struct Specie{T<:Real}
  ρ::T # density
  r::T # radius
  c::Complex{T} # sound speed
  num_density::T # number density
end

"A homogenious background material."
mutable struct Medium{T}
  ρ::T # density
  c::Complex{T} # sound speed
end

"Returns the volume fraction of the specie."
volume_fraction(sp::Specie; dim = 2) = (dim == 2 ? 1 : 4/3) * sp.r^dim * pi * sp.num_density

"Returns pressure wave speed, when β is the adiabatic bulk modulus/"
p_speed(ρ,β,μ) = sqrt(β+4*μ/3)/sqrt(ρ)

function Specie(ρ::T, r::T, c=one(Complex{T}); volfrac::T=0.1*one(T)) where T <: AbstractFloat
  Specie{T}(ρ,r,c,volfrac/(T(pi)*r^2.0))
end

function Specie(T=Float64;ρ = 1.0, r = 1.0, c= 1.0+0.0im, volfrac=0.1)
  if T <: Int T = Float64 end
  Specie{T}(T(ρ),T(r),Complex{T}(c),T(volfrac/(pi*r^2.0)))
end
# function Specie(T=Float64; ρ = zero(T), r = one(T), c=one(Complex{T}), volfrac=0.1*one(T))
#   Specie{T}(T(ρ),T(r),c,T(volfrac/(T(pi)*r^2.0)))
# end

# function Specie(T; ρ = zero(Float64), r = one(Float64), c=one(Complex{Float64}), volfrac=0.1*one(Float64))
#   Specie{Float64}(Float64(ρ),Float64(r),c,Float64(volfrac/(Float64(pi)*r^2.0)))
# end

Medium(;ρ=1.0, c=1.0+0.0im) = Medium{typeof(ρ)}(ρ,Complex{typeof(ρ)}(c))
# Medium(;ρ::T = 1.0, c::Union{R,Complex{T}} = 1.0+0.0im) = Medium(ρ,Complex{T}(c))


text:
function define_compiler(machine)
    define(machine, "create",   op_create)
    define(machine, ":",        op_colon)
    define(machine, ";",        op_semicolon,   true)
    define(machine, ",",        op_comma,       true)
    define(machine, "'",        op_apostrophe,  true)
    define(machine, "[",        op_leftbracket, true)
    define(machine, "]",        op_rightbracket)
    define(machine, "immediate", op_immediate,  true)
    define(machine, "see",      op_see)
    define(machine, "execute",  op_execute)
    define(machine, "postpone", op_postpone,    true)
    define(machine, "postponed", op_postponed,  true)
    define(machine, "mark",     op_mark,        true) # non-std ( -- i ) Push the current index in the just compiled definition
    define(machine, "slot",     op_slot,        true) # non-std ( -- ) Add a slot to the just compiled definition to fill later
    define(machine, "store",    op_store,       true) # non-std ( i v -- ) Store a value in the just compiled definition at i
end

function _compile(engine, wordstr)
    word = codeof(engine, wordstr)
    if isnothing(word)
        print(engine.out, "$(wordstr)?\n")
        return false
    elseif word isa ExecutionToken && word.immediate
        execute_codeword(engine, word) # ???? branching?
        return true
    else
        push!(engine.latest.code, word)
        return true
    end
end

function op_word(machine, parent, myidx)
    wordstr = word(machine.input)[1]
    machine.latest = machine.dictionary[wordstr] = ExecutionToken(wordstr)
    return 1
end

function op_create(machine, parent, myidx)
    wordstr = word(machine.input)[1]
    machine.latest = machine.dictionary[wordstr] = ExecutionToken(wordstr)
    return 1
end

function op_colon(machine, parent, myidx)
    if machine.mode != MODE_INTERPRET
        return print(machine.out, "Colon is only allowed while interpreting.")
    end
    machine.mode = MODE_COMPILE
    return op_create(machine, parent, myidx)
end

function op_semicolon(machine, parent, myidx)
    if machine.mode != MODE_COMPILE
        return print(machine.out, "Semicolon is only allowed while compiling.")
    end
    machine.mode = MODE_INTERPRET
    return 1
end

function op_comma(machine, parent, myidx)
    push!(machine.latest.code, pop!(machine.stack))
    return 1
end

function op_apostrophe(machine, parent, myidx)
    push!(machine.stack, codeof(machine, word(machine.input)[1]))
    return 1
end

function op_leftbracket(machine, parent, myidx)
    machine.mode = MODE_INTERPRET
    return 1
end

function op_rightbracket(machine, parent, myidx)
    machine.mode = MODE_COMPILE
    return 1
end

function op_immediate(machine, parent, myidx)
    if isnothing(machine.latest)
        return print(machine.out, "Invalid IMMEDIATE.")
    end
    machine.latest.immediate = true
    return 1
end

function op_see(machine, parent, myidx)
    what = codeof(machine, word(machine.input)[1])
    println(machine.out, ": $(what.name) " * join(what.code, ' ') * (what.immediate ? " ; immediate" : " ;"))
    return 1
end

function op_execute(machine, _...)
    execute_codeword(machine, pop!(machine.stack))
    return 1
end

function op_postpone(machine, parent, myidx)
    push!(machine.latest.code, codeof(machine, "postponed"))
    push!(machine.latest.code, codeof(machine, word(machine.input)[1]))
    return 1
end

function op_postponed(machine, parent, myidx)
    push!(machine.latest.code, parent.code[myidx + 1])
    return 2
end

function op_mark(machine, parent, myidx)
    push!(machine.stack, length(machine.latest.code))
    return 1
end

function op_slot(machine, parent, myidx)
    push!(machine.latest.code, 0)
    return 1
end

function op_store(machine, parent, myidx)
    val = pop!(machine.stack)
    pos = pop!(machine.stack)::Number
    machine.latest.code[pos] = val
    return 1
end


text:
using Soss

m = @model Prior,x,σ begin
    β ~ Prior
    yhat = β .* x
    y ~ For(eachindex(x)) do j
        Normal(yhat[j], σ)
    end
end

x = randn(3);
Prior = Normal()
truth = rand(m(Prior=Prior, x=x, σ=2.0));

pairs(truth)

post = dynamicHMC(m(Prior=Prior, x=x, σ=2.0), (y=truth.y,));
particles(post)

pred = predictive(m, :β)

predictive(m,:yhat)


text:
const _spa_so_filename =
                string(
                    @__DIR__, "/spa",
                    isapple ? "_macos" : islinux ? "_linux" : "_win",
                    ".so",
                )

_check_spa() = _check_so(_spa_so_filename)
_check_spa()

# export _get_za
# that azimuth is measured westward from south.
_get_za(zenithout::Ref{Float64},azimuthout::Ref{Float64},
        year::Int32,month::Int32,day::Int32,
        hour::Int32,minute::Int32,second::Float64,
        timezone::Float64, delta_ut1::Float64,delta_t::Float64,
        longitude::Float64,latitude::Float64,elevation::Float64,
        pressure::Float64, temperature::Float64,slope::Float64,
        azm_rotation::Float64,atmos_refract::Float64,
        function_int::Int32,) =
            ccall((:get_za, _spa_so_filename),
                Nothing,
                (Int32, Int32, Int32, Int32, Int32, Float64,
                Float64, Float64, Float64,
                Float64, Float64, Float64,
                Float64, Float64, Float64,
                Float64, Float64, Int32,
                Ref{Cdouble}, Ref{Cdouble},),
                year, month,  day,  hour,  minute,  second,
                timezone,  delta_ut1,  delta_t,
                longitude,  latitude,  elevation,
                pressure,  temperature,  slope,
                azm_rotation,  atmos_refract,
                function_int,
                zenithout, azimuthout, )

# export SpaInput

@with_kw mutable struct SpaInput
    year::Int32
    month::Int32
    day::Int32
    hour::Int32
    minute::Int32
    second::Float64
    delta_ut1::Float64 = 0.0
    delta_t::Float64 = 66.6
    timezone::Float64 = 0.0
    longitude::Float64 = 0.0
    latitude::Float64 = 0.0
    elevation::Float64 = 350000.0
    pressure::Float64 = 66.6
    temperature::Float64 = 166.6
    slope::Float64 = 0.0
    azm_rotation::Float64 = 0.0
    atmos_refract::Float64 = 0.5667
    function_int::Int32 = 0

end
function SpaInput(dt::DateTime, pos::LLA)
    SpaInput(year=Dates.year(dt),
            month=Dates.month(dt),
            day=Dates.day(dt),
            hour=Dates.hour(dt),
            minute=Dates.minute(dt),
            second=Dates.second(dt) + Dates.millisecond(dt) / 1000000.0,
            latitude=pos.lat, longitude=pos.lon, elevation=pos.alt,
            )
end
# export get_za!
function get_za!(zenithout::Ref{Float64} ,
                azimuthout::Ref{Float64},
                spa::SpaInput, )::Nothing
    _get_za(zenithout, azimuthout,
            spa.year, spa.month, spa.day, spa.hour, spa.minute, spa.second,
            spa.timezone, spa.delta_ut1, spa.delta_t,
            spa.longitude, spa.latitude, spa.elevation,
            spa.pressure, spa.temperature, spa.slope,
            spa.azm_rotation, spa.atmos_refract,
            spa.function_int,)
    return
end

function test_spa()

    spa = SpaInput(year=2003, month=10, day=17,
            hour=12, minute=0, second=30,
            timezone=0.0, delta_ut1=0, delta_t=67,
            longitude=0.0, latitude=0.0, elevation=1860.14,
            pressure=66.6, temperature=166.6, slope=0.0,
            azm_rotation=0.0, atmos_refract=0.5667,
            function_int=0,)

    zenithout = Ref{Cdouble}(0.0)
    azimuthout = Ref{Cdouble}(0.0)
    get_za!(zenithout, azimuthout, spa, )
    # println("  ", zenithout[], "\t", azimuthout[])
    # 9.936110030111635     202.08736805795803
    zenithout[], azimuthout[]
    # @test isapprox(zenithout[], 9.936110030111635)
    # @test isapprox(azimuthout[],  202.08736805795803)
end


text:
# Functions defined for the distributed primal-dual interior-point algorithm presented in
# the paper "Distributed primal–dual interior-point methods for solving tree-structured coupled convex problems
# using message-passing" by Sina Khoshfetrat Pakazad, Anders Hansson, Martin S. Andersen & Isak Nielsen

# (c): Shervin Parvini Ahmadi

@everywhere function sendto(p::Int; args...)   # This function is taken from the ParallelDataTransfer library on Github and is slightly modified
      for (nm, val) in args
          wait(@spawnat(p, eval(Main, Expr(:(=), nm, val))))
      end
end

@everywhere type Node
  Nodename::Int64                        # Name Nodes with Integers (useful for debugging)
  ProcsId::Int64                         # The processor in which the Node is supposed to live
  var::Vector{Int64}                     # Index of variables
  parent::Union{Node,Void}               # Parent Node
  children::Union{Vector{Node},Void}     # Child Nodes
  ChProcsId::Vector{Int64}               # The processors in which the children Nodes are supposed to live
  Sep                                    # Separator of variables w.r.t variables defined on child Nodes
  Res::Vector{Int64}                     # Residual of variables w.r.t variables defined on parent Node
  H::Array                               # Quadratic term (square matrix) in the objective function
  g::Array                               # Linear term (column vector) in the objective function
  A::Array                               # x's coefficient in the equality constraint
  b::Array                               # Constant term in the equality constraint
  A_I::Array                             # x's coefficient in the inequality constraint
  b_I::Array                             # Constant term in the inequality constraint
  x::Array                               # x
  Delta_x::Array                         # Delata_x
  lambda::Array                          # lambda (the multiplier in Lagrangian which corresponds to the inequality constraint)
  Delta_lambda::Array                    # Delta_lambda
  v::Array                               # v (the multiplier in Lagrangian which corresponds to the equality constraint)
  Delta_v::Array                         # Delta_v
  res_cent::Array                        # Needs to be stored in the upward pass stage in order to be used to calculate Delta_lambda later
  H1::Array                              # Partitioning in the upward pass stage which needs to be stored to be used in the downward-pass stage
  H2::Array                              # Partitioning in the upward pass stage which needs to be stored to be used in the downward-pass stage
  h1::Array                              # Partitioning in the upward pass stage which needs to be stored to be used in the downward-pass stage
  h2::Array                              # Partitioning in the upward pass stage which needs to be stored to be used in the downward-pass stage
end

function setparent(N::Node)
  if ~(N.children == nothing)
    for i=1:length(N.children)
      N.children[i].parent = N;
      setparent(N.children[i]);
    end
  end
end

function setchildrenprocsid(N::Node)
  if ~(N.children == nothing)
    for i=1:length(N.children)
      push!(N.children[i].parent.ChProcsId,N.children[i].ProcsId);
      setchildrenprocsid(N.children[i]);
    end
  end
end

function SepRes(N::Node)

  if N.parent == Void()
    N.Res = N.var
  else
    for i in N.var
      if isempty(find(N.parent.var.==i))
        push!(N.Res,i)
      end
    end
  end
  if ~(N.children==Void())
    for i in 1:length(N.children)
      push!(N.Sep,[]);
    end
    for i in N.var
      for j in 1:length(N.children)
        if ~isempty(find(N.children[j].var.==i))
          push!(N.Sep[j],i)
        end
      end
    end
  end

  if ~(N.children == nothing)
    for k=1:length(N.children)
      SepRes(N.children[k])
    end
  end

end

function Preprocessing_EquCons(N::Node)
  E_RC=zeros(length(N.Res),length(N.var)); for iii=1:length(N.var); if any(N.Res.==N.var[iii]);E_RC[find(N.Res.==N.var[iii]),iii]=1;end;end
  if ~(N.children == nothing)
    for k=1:length(N.children)
      received_message =  Preprocessing_EquCons(N.children[k])
      if ~(received_message == nothing)
        E_SC=zeros(length(N.Sep[k]),length(N.var)); for iii=1:length(N.var); if any(N.Sep[k].==N.var[iii]);E_SC[find(N.Sep[k].==N.var[iii]),iii]=1;end;end
        N.A = [N.A; received_message[1] * E_SC];
        N.b = [N.b; received_message[2]];
      end
    end
  end
  if rank(N.A * E_RC') < size(N.A,1)
    if ~(N.parent == nothing)
      E_SC=zeros(length(setdiff(N.var,N.Res)),length(N.var)); for iii=1:length(N.var); if any(setdiff(N.var,N.Res).==N.var[iii]);E_SC[find(setdiff(N.var,N.Res).==N.var[iii]),iii]=1;end;end
      QR = qrfact(N.A * E_RC');
      trans_A = QR[:Q]'*[N.A*E_RC' N.A*E_SC'];
      trans_b = QR[:Q]' * N.b;
      msg_A = trans_A[rank(N.A * E_RC')+1:end,length(N.Res)+1:end];
      msg_b = trans_b[rank(N.A * E_RC')+1:end,:];
      rank_A_RC = rank(N.A * E_RC')
      N.A = trans_A[1:rank_A_RC,:] * [E_RC;E_SC];
      N.b = trans_b[1:rank_A_RC,:];
      return(msg_A,msg_b)
    else
      QR = qrfact(N.A);
      trans_A = QR[:Q]' * N.A;
      trans_b = QR[:Q]' * N.b;
      rank_A = rank(N.A)
      N.A = trans_A[1:rank_A,:];
      N.b = trans_b[1:rank_A,:];
      return(nothing)
    end
  else
    return(nothing)
  end
end

function find_m_p(N::Node)
  m=size(N.A_I,1);
  p=size(N.A,1);
  if ~(N.children == nothing)
    for i=1:length(N.children)
      m_p = find_m_p(N.children[i])
      m += m_p[1]; p += m_p[2];
    end
  end
  return(m,p)
end

@everywhere function init_lambda_v_x(N::Node,x)
  N.lambda = ones(size(N.A_I,1),1)
  N.v = ones(size(N.A,1),1)
  N.x = x[N.var]
  N.Delta_x = Matrix(length(N.var),1)
  if ~(N.children == nothing)
    for i=1:length(N.children)
      init_lambda_v_x(N.children[i],x)
    end
  end
end

@everywhere function ResidualPDMessagePassing(N::Node,x_loc,lambda,v,t1)
  res_primal = N.A * x_loc - N.b
  if ~isempty(N.A_I)
    res_cent = -diagm(vec(lambda)) * (N.A_I * x_loc - N.b_I) - 1/t1*ones(length(lambda),1)
    res_dual = N.H*x_loc + N.g + N.A_I'*lambda + N.A'*v
  else
    res_dual =  N.H*x_loc + N.g  + N.A'*v
    res_cent = [];
  end

  return(res_primal,res_dual,res_cent)
end

@everywhere function upwardpass(N::Node,t1)
  E_RC=zeros(length(N.Res),length(N.var)); for iii=1:length(N.var); if any(N.Res.==N.var[iii]);E_RC[find(N.Res.==N.var[iii]),iii]=1;end;end
  (res_primal,res_dual,N.res_cent) = ResidualPDMessagePassing(N,N.x,N.lambda,N.v,t1)
  H_msg = N.H
  if ~isempty(N.A_I)
    for j=1:size(N.A_I,1)
      H_msg -= (N.lambda[j]./(N.A_I[j,:]'*N.x-N.b_I[j])).*N.A_I[j,:]*N.A_I[j,:]'
    end
    r_msg = res_dual + N.A_I'*diagm(vec(1./(N.A_I*N.x-N.b_I)))*N.res_cent
  else
    r_msg = res_dual
  end

  if ~(N.children == nothing)

      if any(N.ChProcsId.!=myid())
        abroad_children_indices = find(N.ChProcsId.!=myid())
        msg_Q2_q2 = Vector{Future}(length(abroad_children_indices))
        @sync for (idx, pid) in enumerate(N.ChProcsId[abroad_children_indices])
          @async msg_Q2_q2[idx] = @spawnat pid upwardpass_wraptree(t1)
        end

        if any(N.ChProcsId.==myid())
          for i in find(N.ChProcsId.==myid())
            E_SC=zeros(length(N.Sep[i]),length(N.var)); for iii=1:length(N.var); if any(N.Sep[i].==N.var[iii]);E_SC[find(N.Sep[i].==N.var[iii]),iii]=1;end;end ###
            (msg_Q,msg_q) = upwardpass(N.children[i],t1)
            H_msg += E_SC' * msg_Q * E_SC;
            r_msg += E_SC' * msg_q;
          end
        end

        for (idx,i) in enumerate(abroad_children_indices)
          (msg_Q2,msg_q2) = fetch(msg_Q2_q2[idx])
          E_SC=zeros(length(N.Sep[i]),length(N.var)); for iii=1:length(N.var); if any(N.Sep[i].==N.var[iii]);E_SC[find(N.Sep[i].==N.var[iii]),iii]=1;end;end
          H_msg += E_SC' * msg_Q2 * E_SC;
          r_msg += E_SC' * msg_q2;
        end

    else
      for i=1:length(N.children)
        E_SC=zeros(length(N.Sep[i]),length(N.var)); for iii=1:length(N.var); if any(N.Sep[i].==N.var[iii]);E_SC[find(N.Sep[i].==N.var[iii]),iii]=1;end;end
        (msg_Q,msg_q) = upwardpass(N.children[i],t1)
        H_msg += E_SC' * msg_Q * E_SC;
        r_msg += E_SC' * msg_q;
      end
    end
  end

  E_SC=zeros(length(setdiff(N.var,N.Res)),length(N.var)); for iii=1:length(N.var); if any(setdiff(N.var,N.Res).==N.var[iii]);E_SC[find(setdiff(N.var,N.Res).==N.var[iii]),iii]=1;end;end

  Q_yy = E_SC * H_msg * E_SC';
  q_y = E_SC * r_msg;
  A_y = N.A * E_SC';
  Q_zy= E_RC * H_msg * E_SC';
  Q_zz = E_RC * H_msg * E_RC';
  q_z = E_RC * r_msg;
  A_z = N.A * E_RC';

  O = [Q_zz A_z';A_z zeros(size(A_z,1),size(A_z,1))];

  H_opt = -O \ [Q_zy;A_y];
  N.H1 = H_opt[1:length(N.Res),:];
  N.H2 = H_opt[length(N.Res)+1:end,:];

  h_opt = -O \ [q_z;res_primal];
  N.h1 = h_opt[1:length(N.Res),:];
  N.h2 = h_opt[length(N.Res)+1:end,:];

  msg_Q = N.H1' * Q_zz * N.H1 + N.H1' * Q_zy + Q_zy' * N.H1 + Q_yy;
  msg_q = N.H1' * Q_zz * N.h1 + N.H1' * q_z + Q_zy' * N.h1 + q_y;
  return(msg_Q,msg_q)

end

@everywhere function upwardpass_wraptree(t1)
  get_true_tree = subtree
  return(upwardpass(get_true_tree,t1))
end

@everywhere function downwardpass(N::Node,Delta_x_send)

  if setdiff(N.var,N.Res)==[]
    N.Delta_x =  N.h1
    N.Delta_v =  N.h2
  else
    N.Delta_x[findin(N.var,N.Res)] = N.H1 * Delta_x_send + N.h1
    N.Delta_x[findin(N.var,setdiff(N.var,N.Res))] = Delta_x_send
    N.Delta_v = N.H2 * Delta_x_send + N.h2
  end

  if ~(N.children == nothing)
    if any(N.ChProcsId.!=myid())
      abroad_children_indices = find(N.ChProcsId.!=myid())
      confirm = Vector{Future}(length(abroad_children_indices))
      for (idx, pid) in enumerate(N.ChProcsId[abroad_children_indices])
        to_send1 = N.Delta_x[findin(N.var,N.Sep[abroad_children_indices[idx]])]
		confirm[idx] = @spawnat pid downwardpass_wraptree(to_send1)
      end

      if any(N.ChProcsId.==myid())
        for i in find(N.ChProcsId.==myid())
          to_send2 = N.Delta_x[findin(N.var,N.Sep[i])]
		  downwardpass(N.children[i],to_send2)
        end
      end
      @sync for idx = 1:length(abroad_children_indices)
        @async wait(confirm[idx])
      end

    else
      for i=1:length(N.children)
        to_send3 = N.Delta_x[findin(N.var,N.Sep[i])]
		downwardpass(N.children[i],to_send3)
      end
    end
  end
end

@everywhere function downwardpass_wraptree(Delta_x_send)
  get_true_tree = subtree
  return(downwardpass(get_true_tree,Delta_x_send))
end

function MessagePassingDirection(N::Node,t1,n)
  upwardpass(N,t1)
  Delta_x_tot = downwardpass(N,nothing)
end

@everywhere function  find_Delta_lambda(N::Node)
  N.Delta_lambda = -diagm(vec(1./(N.A_I*N.x-N.b_I)))*(diagm(vec(N.lambda))*N.A_I*N.Delta_x - N.res_cent)
  if ~(N.children == nothing)
    if any(N.ChProcsId.!=myid())
      abroad_children_indices = find(N.ChProcsId.!=myid())
      confirm = Vector{Future}(length(abroad_children_indices))
      @sync for (idx, pid) in enumerate(N.ChProcsId[abroad_children_indices])
          @async confirm[idx] = @spawnat pid find_Delta_lambda_wraptree()
      end
      if any(N.ChProcsId.==myid())
        for i in find(N.ChProcsId.==myid())
          find_Delta_lambda(N.children[i])
        end
      end
      @sync for idx = 1:length(abroad_children_indices)
        @async wait(confirm[idx])
      end

    else
      for i=1:length(N.children)
        find_Delta_lambda(N.children[i])
      end
    end
  end
  return(nothing)
end

@everywhere function find_Delta_lambda_wraptree()
  get_true_tree = subtree
  return(find_Delta_lambda(get_true_tree))
end

@everywhere function stepsize_firststage(N::Node,Beta)
  Temp = -N.lambda./N.Delta_lambda
  t_d = 1
  if ~isempty(Temp[N.Delta_lambda.<0])
    t_d = min(1,minimum(Temp[N.Delta_lambda.<0]))
  end
  t = 0.99 * t_d
  while sum(N.A_I*(N.x+t*N.Delta_x)-N.b_I.>=0)>0
    t = Beta * t;
  end

  if ~(N.children == nothing)
    if any(N.ChProcsId.!=myid())
      abroad_children_indices = find(N.ChProcsId.!=myid())
      t_msg = Vector{Future}(length(abroad_children_indices))
      @sync for (idx, pid) in enumerate(N.ChProcsId[abroad_children_indices])
        @async t_msg[idx] = @spawnat pid stepsize_firststage_wraptree(Beta)
      end
      if any(N.ChProcsId.==myid())
        for i in find(N.ChProcsId.==myid())
          t = min(t,stepsize_firststage(N.children[i],Beta))
        end
      end
      for idx=1:length(abroad_children_indices)
        t = min(t,fetch(t_msg[idx]))
      end

    else
      for i=1:length(N.children)
        t = min(t,stepsize_firststage(N.children[i],Beta))
      end
    end

  end
  return(t)
end

@everywhere function stepsize_firststage_wraptree(Beta)
  get_true_tree = subtree
  return(stepsize_firststage(get_true_tree,Beta))
end

@everywhere function TotalResiduals(N::Node,t1,n)
  (res_primal_tot,res_dual1,res_cent_tot) = ResidualPDMessagePassing(N,N.x,N.lambda,N.v,t1)
  res_dual_tot=zeros(n,1);
  kk=1;for ii in N.var; res_dual_tot[ii,1]+=res_dual1[kk,1];kk+=1;end
  if ~(N.children == nothing)
    if any(N.ChProcsId.!=myid())
      abroad_children_indices = find(N.ChProcsId.!=myid())
      res_msg_reciept = Vector{Future}(length(abroad_children_indices))
      @sync for (idx, pid) in enumerate(N.ChProcsId[abroad_children_indices])
        @async res_msg_reciept[idx] = @spawnat pid TotalResiduals_wraptree(t1,n)
      end
      if any(N.ChProcsId.==myid())
        for i in find(N.ChProcsId.==myid())
          (res_primal,res_dual2,res_cent) = TotalResiduals(N.children[i],t1,n)
          res_primal_tot = [res_primal_tot;res_primal]
          res_dual_tot += res_dual2
          res_cent_tot = [res_cent_tot;res_cent]
        end
      end
      for idx=1:length(abroad_children_indices)
        res_msg = fetch(res_msg_reciept[idx])
        res_primal_tot = [res_primal_tot;res_msg[1]]
        res_dual_tot += res_msg[2]
        res_cent_tot = [res_cent_tot;res_msg[3]]
      end

    else
      for i=1:length(N.children)
        (res_primal,res_dual2,res_cent) = TotalResiduals(N.children[i],t1,n)
        res_primal_tot = [res_primal_tot;res_primal]
        res_dual_tot += res_dual2
        res_cent_tot = [res_cent_tot;res_cent]
      end
    end
  end
  return(res_primal_tot,res_dual_tot,res_cent_tot)
end

@everywhere function TotalResiduals_wraptree(t1,n)
  get_true_tree = subtree
  return(TotalResiduals(get_true_tree,t1,n))
end

@everywhere function TotalResiduals_Ahead(N::Node,t1,n,t_step)
  (res_primal_tot,res_dual1,res_cent_tot) = ResidualPDMessagePassing(N,N.x+t_step*N.Delta_x,N.lambda+t_step*N.Delta_lambda,N.v+t_step*N.Delta_v,t1)
  res_dual_tot=zeros(n,1);
  kk=1;for ii in N.var; res_dual_tot[ii,1]+=res_dual1[kk,1];kk+=1;end
  if ~(N.children == nothing)
    if any(N.ChProcsId.!=myid())
      abroad_children_indices = find(N.ChProcsId.!=myid())
      res_msg_reciept = Vector{Future}(length(abroad_children_indices))
      @sync for (idx, pid) in enumerate(N.ChProcsId[abroad_children_indices])
        @async res_msg_reciept[idx] = @spawnat pid TotalResiduals_Ahead_wraptree(t1,n,t_step)
      end
      if any(N.ChProcsId.==myid())
        for i in find(N.ChProcsId.==myid())
          (res_primal,res_dual2,res_cent) = TotalResiduals_Ahead(N.children[i],t1,n,t_step)
          res_primal_tot = [res_primal_tot;res_primal]
          res_dual_tot += res_dual2
          res_cent_tot = [res_cent_tot;res_cent]
        end
      end
      for idx=1:length(abroad_children_indices)
        res_msg = fetch(res_msg_reciept[idx])
        res_primal_tot = [res_primal_tot;res_msg[1]]
        res_dual_tot += res_msg[2]
        res_cent_tot = [res_cent_tot;res_msg[3]]
      end
    else
      for i=1:length(N.children)
        (res_primal,res_dual2,res_cent) = TotalResiduals_Ahead(N.children[i],t1,n,t_step)
        res_primal_tot = [res_primal_tot;res_primal]
        res_dual_tot += res_dual2
        res_cent_tot = [res_cent_tot;res_cent]
      end
    end

  end
  return(res_primal_tot,res_dual_tot,res_cent_tot)
end

@everywhere function TotalResiduals_Ahead_wraptree(t1,n,t_step)
  get_true_tree = subtree
  return(TotalResiduals_Ahead(get_true_tree,t1,n,t_step))
end

function MessagePassingStepSize(N::Node,t1,alpha,Beta,n)
  t_final = stepsize_firststage(N,Beta)

  residual1 = TotalResiduals(N,t1,n)
  res = norm([residual1[1];residual1[2];residual1[3]])^2
  residual2 = TotalResiduals_Ahead(N,t1,n,t_final)
  res_t = norm([residual2[1];residual2[2];residual2[3]])^2

  back_no = 0;
  if (res_t > (1-alpha*t_final)^2*res)
    back_no = 1;
  end
  return(t_final,back_no)
end

@everywhere function find_surduality(N::Node)
  sur_duality = -( (N.A_I*N.x - N.b_I)'*N.lambda );
  if ~(N.children == nothing)
    if any(N.ChProcsId.!=myid())
      abroad_children_indices = find(N.ChProcsId.!=myid())
      sur_duality_msg_reciept = Vector{Future}(length(abroad_children_indices))
      @sync for (idx, pid) in enumerate(N.ChProcsId[abroad_children_indices])
        @async sur_duality_msg_reciept[idx] = @spawnat pid find_surduality_wraptree()
      end
      if any(N.ChProcsId.==myid())
        for i in find(N.ChProcsId.==myid())
          sur_duality += find_surduality(N.children[i])
        end
      end
      for idx=1:length(abroad_children_indices)
        sur_duality += fetch(sur_duality_msg_reciept[idx])
      end

    else
      for i=1:length(N.children)
        sur_duality += find_surduality(N.children[i])
      end
    end

  end
  sur_duality = reshape(sur_duality,1)[1]
  return(sur_duality)
end

@everywhere function find_surduality_wraptree()
  get_true_tree = subtree
  return(find_surduality(get_true_tree))
end

@everywhere function Update_lambda_v_x(N::Node,stepsize)
  N.lambda += stepsize * N.Delta_lambda
  N.v += stepsize * N.Delta_v
  N.x += stepsize * N.Delta_x
  if ~(N.children == nothing)
    if any(N.ChProcsId.!=myid())
      abroad_children_indices = find(N.ChProcsId.!=myid())
      confirm = Vector{Future}(length(abroad_children_indices))
      @sync for (idx, pid) in enumerate(N.ChProcsId[abroad_children_indices])
        @async confirm[idx] = @spawnat pid Update_lambda_v_x_wraptree(stepsize)
      end
      if any(N.ChProcsId.==myid())
        for i in find(N.ChProcsId.==myid())
          Update_lambda_v_x(N.children[i],stepsize)
        end
      end
      @sync for idx = 1:length(abroad_children_indices)
        @async wait(confirm[idx])
      end

    else
      for i=1:length(N.children)
        Update_lambda_v_x(N.children[i],stepsize)
      end
    end

  end
  return(nothing)
end

@everywhere function Update_lambda_v_x_wraptree(stepsize)
  get_true_tree = subtree
  return(Update_lambda_v_x(get_true_tree,stepsize))
end

max_stars_repo_name: ShervinParvini/DistPrimDualIntPoint


text:
# an AbstractIndex is a thing that can be used to look up ordered things by name, but that
# will also accept a position or set of positions or range or other things and pass them
# through cleanly.
# an Index is the usual implementation.
# a SimpleIndex only works if the things are integer indexes, which is weird.
typealias Indices Union(Real, AbstractVector{Real})

abstract AbstractIndex

type Index <: AbstractIndex   # an OrderedDict would be nice here...
    lookup::Dict{Symbol, Indices}      # name => names array position
    names::Vector{Symbol}
end
function Index(x::Vector{Symbol})
    for n in x
        if !is_valid_identifier(n)
            error("Names must be valid identifiers.")
        end
    end
    x = make_unique(x)
    Index(Dict{Symbol, Indices}(tuple(x...), tuple([1:length(x)]...)), x)
end
Index() = Index(Dict{Symbol, Indices}(), Symbol[])
Base.length(x::Index) = length(x.names)
Base.names(x::Index) = copy(x.names)
Base.copy(x::Index) = Index(copy(x.lookup), copy(x.names))
Base.deepcopy(x::Index) = Index(deepcopy(x.lookup), deepcopy(x.names))
Base.isequal(x::Index, y::Index) = isequal(x.lookup, y.lookup) && isequal(x.names, y.names)
Base.(:(==))(x::Index, y::Index) = isequal(x, y)

function names!(x::Index, nm::Vector{Symbol})
    if length(nm) != length(x)
        error("Lengths don't match.")
    end
    for n in nm
        if !is_valid_identifier(n)
            error("Names must be valid identifiers.")
        end
    end
    for i in 1:length(nm)
        delete!(x.lookup, x.names[i])
        x.lookup[nm[i]] = i
    end
    x.names = nm
    return x
end

function rename!(x::Index, nms)
    for (from, to) in nms
        if haskey(x, from)
            if haskey(x, to)
                error("Tried renaming $from to $to, when $to already exists in the Index.")
            end
            if !is_valid_identifier(to)
                error("Names must be valid identifiers.")
            end
            x.lookup[to] = col = pop!(x.lookup, from)
            if !isa(col, Array)
                x.names[col] = to
            end
        end
    end
    return x
end

rename!(x::Index, from, to) = rename!(x, zip(from, to))
rename!(x::Index, from::Symbol, to::Symbol) = rename!(x, ((from, to),))
rename!(x::Index, f::Function) = rename!(x, [(x,f(x)) for x in x.names])

rename(x::Index, args...) = rename!(copy(x), args...)

Base.haskey(x::Index, key::Symbol) = haskey(x.lookup, key)
Base.haskey(x::Index, key::Real) = 1 <= key <= length(x.names)
Base.keys(x::Index) = names(x)

function Base.push!(x::Index, nm::Symbol)
    x.lookup[nm] = length(x) + 1
    push!(x.names, nm)
    return x
end

function Base.delete!(x::Index, idx::Integer)
    # reset the lookup's beyond the deleted item
    for i in (idx + 1):length(x.names)
        x.lookup[x.names[i]] = i - 1
    end
    delete!(x.lookup, x.names[idx])
    splice!(x.names, idx)
    return x
end

function Base.delete!(x::Index, nm::Symbol)
    if !haskey(x.lookup, nm)
        return x
    end
    idx = x.lookup[nm]
    return delete!(x, idx)
end

Base.getindex(x::Index, idx::Symbol) = x.lookup[idx]
Base.getindex(x::AbstractIndex, idx::Real) = int(idx)
Base.getindex(x::AbstractIndex, idx::AbstractDataVector{Bool}) = getindex(x, array(idx, false))
Base.getindex{T}(x::AbstractIndex, idx::AbstractDataVector{T}) = getindex(x, dropna(idx))
Base.getindex(x::AbstractIndex, idx::AbstractVector{Bool}) = find(idx)
Base.getindex(x::AbstractIndex, idx::Ranges) = [idx]
Base.getindex{T <: Real}(x::AbstractIndex, idx::AbstractVector{T}) = convert(Vector{Int}, idx)
Base.getindex(x::AbstractIndex, idx::AbstractVector{Symbol}) = [[x.lookup[i] for i in idx]...]

type SimpleIndex <: AbstractIndex
    length::Integer
end
SimpleIndex() = SimpleIndex(0)
Base.length(x::SimpleIndex) = x.length
Base.names(x::SimpleIndex) = nothing


text:
export dual
export truncate, ambo, bevel, expand, snub
export kiss, join, meta, ortho, gyro

function dual end

function truncate end
function ambo end
function bevel end
function expand end
function snub end

function kiss end
function join end
function meta end
function ortho end
function gyro end


text:
strrep(s::String, r::UnitRange) = [repeat(s, n) for n in r]
naValues = vcat(strrep("*", 1:23), strrep("#", 1:23), "NAME WITHHELD BY AGENCY", "NAME WITHHELD BY OPM", "NAME UNKNOWN", "UNSP", "<NA>", "000000", "999999", "")

#    ~FieldName,~Length, ~Type,
CPDFStandardFormat = DataFrame(["SCT_TABLE_ID"  2  String;
    "DATA_CODE" 9  String;
    "CODE_USE_FROM 1"  6  DateFormat("yyyymm");
    "CODE_USE_UNTIL 1"  6 DateFormat("yyyymm");
    "CODE_USE_FROM 2"  6  DateFormat("yyyymm");
    "CODE_USE_UNTIL 2"  6  DateFormat("yyyymm");
    "TRANSLATION_1"  40  String;
    "TRANSLATION_IN_USE_FROM 1"  6  DateFormat("yyyymm");
    "TRANSLATION_IN_USE_UNTIL 1"  6  DateFormat("yyyymm");
    "TRANSLATION_2"  40  String;
    "TRANSLATION_IN_USE_FROM 2"  6  DateFormat("yyyymm");
    "TRANSLATION_IN_USE_UNTIL 2"  6  DateFormat("yyyymm");
    "TRANSLATION_3"  40  String;
    "TRANSLATION_IN_USE_FROM 3"  6  DateFormat("yyyymm");
    "TRANSLATION_IN_USE_UNTIL 3"  6  DateFormat("yyyymm");
    "TRANSLATION_4"  40  String;
    "TRANSLATION_IN_USE_FROM 4"  6  DateFormat("yyyymm");
    "TRANSLATION_IN_USE_UNTIL 4"  6  DateFormat("yyyymm");
    "TRANSLATION_5"  40  String;
    "TRANSLATION_IN_USE_FROM 5"  6  DateFormat("yyyymm");
    "TRANSLATION_IN_USE_UNTIL 5"  6  DateFormat("yyyymm");
    "TRANSLATION_6"  40  String;
    "TRANSLATION_IN_USE_FROM 6"  6  DateFormat("yyyymm");
    "TRANSLATION_IN_USE_UNTIL 6"  6 DateFormat("yyyymm")])

# Used to use a tibble, now use julia vectors...
#    ~FieldName,~Length,~Type,
file_format = DataFrame(["PSEUDO_ID" 9 Int;
    "EMPLOYEE_NAME" 23 String;
    "FILE_DATE" 8  DateFormat("yyyymmdd");
    "AGENCY" 2 String;
    "SUB_AGENCY" 2 String;
    "DUTY_STATION" 9 String;
    "AGE" 6 String;
    "EDUCATION_LEVEL" 2 String;
    "PAY_PLAN" 2 String;
    "GRADE" 2 String;
    "LOS_LEVEL" 6 String;
    "OCCUPATION" 4 String;
    "OCCUPATIONAL_CATEGORY" 1 String;
    "ADJUSTED_BASIC_PAY" 6 Int;
    "SUPERVISORY_STATUS" 1 String;
    "TYPE_OF_APPOINTMENT" 2 String;
    "WORK_SCHEDULE" 1 String;
    "NSFTP_IND"  1 String])

#file = raw"C:\Users\matth\OneDrive\Documents\julia\FWF\test\testfiles\Status_Non_DoD_1981_03.txt"
#FWF.read(file, convert(Array{Int},OPMFedFormat[:x2]),
#        header=convert(Array{String}, OPMFedFormat[:x1]), types=convert(Array{Union{Type, DateFormat}},OPMFedFormat[:x3]),
#        missings=naValues)

#FWF.read(IOBuffer("abc12310102017\ndef45610112017\n"), [3,3,8], types=[String,Int,DateFormat("mmddyyyy")])

# This function will parse and return a tibble of a file with the assumption it is a FedScope fixed width
# file.  Results will be unknown if the file is not a FedScope fixed width file.
# This function also has the ability to expand collapsed fields using helper functions.
# Parameters:
# path: The path to the file that will be read.
# post2014: boolean flag that signals the file is after 2014 and thus contains "Years Since Degree"
# incYearsSinceDegree: Will include the years since degree column as NA if not present.
function opm_parse_fwf(path, incYearsSinceDegree = true)
    # removed post 2014, files we deal with were always pre-2014
    data = file_format

    #colsObj = mapply(function(x, y) { y }, pull(data, FieldName), pull(data, Type), SIMPLIFY = FALSE, USE.NAMES = TRUE)
    result = FWF.read(path, convert(Array{Int},data[:x2]),
        header=convert(Array{String}, data[:x1]), types=convert(Array{Union{Type, DateFormat}},data[:x3]), missings=naValues)
    #if (incYearsSinceDegree)
    #    result = result %>% add_column(`YEARS SINCE DEGREE` = NA, .after = "AGE") %>% mutate(`YEARS SINCE DEGREE` = as.integer(`YEARS SINCE DEGREE`))
    #    end
    return result
end

# This function will read in the cpdf format file named sctdata.  This is a static data table of
# normalized values from the main dataset.
function opm_parse_sctfile_fwf(path)
    data = CPDFStandardFormat
    # Ugly line of code to convert a tbl to parameter list
    #colsObj <- mapply(function(x, y) { y }, pull(data, FieldName), pull(data, Type), SIMPLIFY = FALSE, USE.NAMES = TRUE)
    result = FWF.read(path, convert(Array{Int},data[:x2]),
        header=convert(Array{String}, data[:x1]), types=convert(Array{Union{Type, DateFormat}},data[:x3]), missings=naValues)
    return result
end
# sct = opm_parse_sctfile_fwf(raw"C:\Users\matth\Desktop\data\1973-09-to-2014-06\SCTFILE.TXT")


text:
"""
    mutable struct Scenarios <: Forecast
        name::String
        resolution::Dates.Period
        scenario_count::Int64
        data::Union{
            SortedDict{Dates.DateTime, Matrix{CONSTANT}},
            SortedDict{Dates.DateTime, Matrix{POLYNOMIAL}},
            SortedDict{Dates.DateTime, Matrix{PWL}},
        }
        scaling_factor_multiplier::Union{Nothing, Function}
        internal::InfrastructureSystemsInternal
    end

A Discrete Scenario Based time series for a particular data field in a Component.

# Arguments
- `name::String`: user-defined name
- `resolution::Dates.Period`: forecast resolution
- `scenario_count::Int64`: Number of scenarios
- `data::Union{SortedDict{Dates.DateTime, Matrix{CONSTANT}}, SortedDict{Dates.DateTime, Matrix{POLYNOMIAL}}, SortedDict{Dates.DateTime, Matrix{PWL}}}`: timestamp - scalingfactor
- `scaling_factor_multiplier::Union{Nothing, Function}`: Applicable when the time series
  data are scaling factors. Called on the associated component to convert the values.
- `internal::InfrastructureSystemsInternal`
"""
mutable struct Scenarios <: Forecast
    "user-defined name"
    name::String
    "timestamp - scalingfactor"
    data::Union{
        SortedDict{Dates.DateTime, Matrix{CONSTANT}},
        SortedDict{Dates.DateTime, Matrix{POLYNOMIAL}},
        SortedDict{Dates.DateTime, Matrix{PWL}},
    }
    "Number of scenarios"
    scenario_count::Int64
    "forecast resolution"
    resolution::Dates.Period
    "Applicable when the time series data are scaling factors. Called on the associated component to convert the values."
    scaling_factor_multiplier::Union{Nothing, Function}
    internal::InfrastructureSystemsInternal
end

function Scenarios(;
    name,
    data,
    scenario_count,
    resolution,
    scaling_factor_multiplier = nothing,
    normalization_factor = 1.0,
    internal = InfrastructureSystemsInternal(),
)
    data = handle_normalization_factor(convert_data(data), normalization_factor)
    Scenarios(name, data, scenario_count, resolution, scaling_factor_multiplier, internal)
end

"""
Construct Scenarios from a SortedDict of Arrays.

# Arguments
- `name::AbstractString`: user-defined name
- `input_data::AbstractDict{Dates.DateTime, Matrix{Float64}}`: time series data.
- `resolution::Dates.Period`: The resolution of the forecast in Dates.Period`
- `normalization_factor::NormalizationFactor = 1.0`: optional normalization factor to apply
  to each data entry
- `scaling_factor_multiplier::Union{Nothing, Function} = nothing`: If the data are scaling
  factors then this function will be called on the component and applied to the data when
  [`get_time_series_array`](@ref) is called.
"""
function Scenarios(
    name::AbstractString,
    input_data::AbstractDict,
    resolution::Dates.Period;
    normalization_factor::NormalizationFactor = 1.0,
    scaling_factor_multiplier::Union{Nothing, Function} = nothing,
)
    scenario_count = size(first(values(input_data)))[2]

    return Scenarios(
        name = name,
        data = input_data,
        scenario_count = scenario_count,
        resolution = resolution,
        scaling_factor_multiplier = scaling_factor_multiplier,
        normalization_factor = normalization_factor,
        internal = InfrastructureSystemsInternal(),
    )
end

"""
Construct Scenarios from a Dict of TimeArrays.

# Arguments
- `name::AbstractString`: user-defined name
- `input_data::AbstractDict{Dates.DateTime, TimeSeries.TimeArray}`: time series data.
- `normalization_factor::NormalizationFactor = 1.0`: optional normalization factor to apply
  to each data entry
- `scaling_factor_multiplier::Union{Nothing, Function} = nothing`: If the data are scaling
  factors then this function will be called on the component and applied to the data when
  [`get_time_series_array`](@ref) is called.
- `timestamp = :timestamp`: If the values are DataFrames is passed then this must be the column name that
  contains timestamps.
"""
function Scenarios(
    name::AbstractString,
    input_data::AbstractDict{Dates.DateTime, <:TimeSeries.TimeArray};
    normalization_factor::NormalizationFactor = 1.0,
    scaling_factor_multiplier::Union{Nothing, Function} = nothing,
)
    data = SortedDict{Dates.DateTime, Matrix{Float64}}()
    resolution =
        TimeSeries.timestamp(first(values(input_data)))[2] -
        TimeSeries.timestamp(first(values(input_data)))[1]
    for (k, v) in input_data
        data[k] = TimeSeries.values(v)
    end

    return Scenarios(
        name,
        data,
        resolution;
        normalization_factor = normalization_factor,
        scaling_factor_multiplier = scaling_factor_multiplier,
    )
end

function Scenarios(ts_metadata::ScenariosMetadata, data::SortedDict)
    return Scenarios(
        name = get_name(ts_metadata),
        scenario_count = get_scenario_count(ts_metadata),
        resolution = get_resolution(ts_metadata),
        data = data,
        scaling_factor_multiplier = get_scaling_factor_multiplier(ts_metadata),
        internal = InfrastructureSystemsInternal(get_time_series_uuid(ts_metadata)),
    )
end

function Scenarios(info::TimeSeriesParsedInfo)
    return Scenarios(
        info.name,
        info.data,
        info.resolution;
        normalization_factor = info.normalization_factor,
        scaling_factor_multiplier = info.scaling_factor_multiplier,
    )
end

function ScenariosMetadata(time_series::Scenarios)
    return ScenariosMetadata(
        get_name(time_series),
        get_resolution(time_series),
        get_initial_timestamp(time_series),
        get_interval(time_series),
        get_scenario_count(time_series),
        get_count(time_series),
        get_uuid(time_series),
        get_horizon(time_series),
        get_scaling_factor_multiplier(time_series),
    )
end

function get_array_for_hdf(forecast::Scenarios)
    interval_count = get_count(forecast)
    scenario_count = get_scenario_count(forecast)
    horizon = get_horizon(forecast)
    data = get_data(forecast)

    data_for_hdf = Array{Float64, 3}(undef, scenario_count, horizon, interval_count)
    for (ix, f) in enumerate(values(data))
        data_for_hdf[:, :, ix] = transpose(f)
    end
    return data_for_hdf
end

function get_horizon(forecast::Scenarios)
    return size(first(values(get_data(forecast))))[1]
end

"""Get [`Scenarios`](@ref) `name`."""
get_name(value::Scenarios) = value.name
"""Get [`Scenarios`](@ref) `resolution`."""
get_resolution(value::Scenarios) = value.resolution
"""Get [`Scenarios`](@ref) `scenario_count`."""
get_scenario_count(value::Scenarios) = value.scenario_count
"""Get [`Scenarios`](@ref) `data`."""
get_data(value::Scenarios) = value.data
"""Get [`Scenarios`](@ref) `scaling_factor_multiplier`."""
get_scaling_factor_multiplier(value::Scenarios) = value.scaling_factor_multiplier
"""Get [`Scenarios`](@ref) `internal`."""
get_internal(value::Scenarios) = value.internal

"""Set [`Scenarios`](@ref) `name`."""
set_name!(value::Scenarios, val) = value.name = val
"""Set [`Scenarios`](@ref) `resolution`."""
set_resolution!(value::Scenarios, val) = value.resolution = val
"""Set [`Scenarios`](@ref) `scenario_count`."""
set_scenario_count!(value::Scenarios, val) = value.scenario_count = val
"""Set [`Scenarios`](@ref) `data`."""
set_data!(value::Scenarios, val) = value.data = val
"""Set [`Scenarios`](@ref) `scaling_factor_multiplier`."""
set_scaling_factor_multiplier!(value::Scenarios, val) =
    value.scaling_factor_multiplier = val
"""Set [`Scenarios`](@ref) `internal`."""
set_internal!(value::Scenarios, val) = value.internal = val

eltype_data(forecast::Scenarios) = eltype_data_common(forecast)
get_count(forecast::Scenarios) = get_count_common(forecast)
get_initial_times(forecast::Scenarios) = get_initial_times_common(forecast)
get_initial_timestamp(forecast::Scenarios) = get_initial_timestamp_common(forecast)
get_interval(forecast::Scenarios) = get_interval_common(forecast)
get_window(f::Scenarios, initial_time::Dates.DateTime; len = nothing) =
    get_window_common(f, initial_time; len = len)
iterate_windows(forecast::Scenarios) = iterate_windows_common(forecast)


text:
module TestRecording

using StopWatches
using StopWatchesBase: Record, StopWatchesBase
using Test

using ..Utils: @spawn

function test_allocation_hack()
    swb = StopWatchesBase
    @test swb._sample_record1() === swb._sample_record1()::Record
    @test swb._sample_record2() === swb._sample_record2()::Record
    @test swb._sample_record3() === swb._sample_record3()::Record
    @test swb._sample_record1().location.fullname[1] === fullname(swb)[1]
    @test swb._sample_record1().location.tag === nothing
    @test swb._sample_record2().location.tag === :sample_record2
    @test swb._sample_record1().location.line + 1 == swb._sample_record2().location.line
end

function record_once()
    ans = @stopwatch :record_once begin
        y = 1
    end
    return (ans, y)
end

function test_record_once()
    @test record_once() == (1, 1)
end

function no_tag()
    ans = @stopwatch y = 1
    return (ans, y)
end

function test_no_tag()
    @test no_tag() == (1, 1)
end

function test_record_in_spawn()
    tasks = map(1:10) do i
        @spawn begin
            ans = @stopwatch :record_in_spawn begin
                y = i * 10
            end
            return (ans, y)
        end
    end
    results = fetch.(tasks)
    @test first.(results) == (1:10) .* 10
    @test last.(results) == (1:10) .* 10
end

end  # module


text:
module Cambrian

using Random
using Logging
using Statistics
using YAML
using JSON
import Formatting
import Dates

include("config.jl")
include("logger.jl")
include("individual.jl")
include("evolution.jl")
include("selection.jl")
include("evaluation.jl")
include("step.jl")
include("oneplus.jl")
include("GA.jl")

end


text:
__precompile__(true)

module Format

import Base.show
using Compat
using Compat.Printf
const V6_COMPAT = VERSION < v"0.7.0-DEV"

_stdout() = @static V6_COMPAT ? STDOUT : stdout
_codeunits(s) = Vector{UInt8}(@static V6_COMPAT ? s : codeunits(s))
@static if V6_COMPAT
    const m_eval = eval
else
    m_eval(expr) = Core.eval(@__MODULE__, expr)
end

export FormatSpec, FormatExpr, printfmt, printfmtln, format, generate_formatter
export pyfmt, cfmt, fmt
export fmt_default, fmt_default!, reset!, default_spec, default_spec!

# Deal with mess from #16058
# Later, use Strs package!
isdefined(Main, :ASCIIStr) || (const ASCIIStr = String)
isdefined(Main, :UTF8Str)  || (const UTF8Str = String)
isdefined(Main, :AbstractChar) || (const AbstractChar = Char)

include("cformat.jl" )
include("fmtspec.jl")
include("fmtcore.jl")
include("formatexpr.jl")
include("fmt.jl")

end # module Format


text:
module Catalyst

using Reexport, ModelingToolkit
@reexport using ModelingToolkit
import MacroTools
import Base: (==), merge!, merge
using Latexify

const ExprValues = Union{Expr,Symbol,Float64,Int}

include("expression_utils.jl")
include("reaction_network.jl")

# reaction network macro
export @reaction_network, @add_reactions
export @reaction_func

# functions to query network properties
include("networkapi.jl")
export species, params, reactions, speciesmap, paramsmap, numspecies, numreactions, numparams
export make_empty_network, addspecies!, addparam!, addreaction!
export dependants, dependents

# for Latex printing of ReactionSystems
include("latexify_recipes.jl")

# for making and saving graphs
import Base.Iterators: flatten
using Catlab.Graphics.Graphviz
import Catlab.Graphics.Graphviz: Graph, Edge
include("graphs.jl")
export Graph, savegraph

end # module


text:
module FeedbackChains

using Flux
import Flux: functor
import Base: getindex, show
using MacroTools: @forward

using ..Splitters
using ..AbstractMergers
using ..AbstractFeedbackNets
export FeedbackChain

"""
    FeedbackChain{T<:Tuple}

Tuple-like structure similar to a Flux.Chain with support for `Splitter`s and `Merger`s.
"""
struct FeedbackChain{T<:Tuple} <: AbstractFeedbackNet
    layers::T
    FeedbackChain(xs...) = new{typeof(xs)}(xs)
end

# empty feedback chain just returns the input
function (c::FeedbackChain{Tuple{}})(h, x)
    return h, x
end # function (c::FeedbackChain{Tuple{}})

"""
    (c::FeedbackChain)(h, x)

Apply a `FeedbackChain` to input `x` with hidden state `h`. `h` should take the
form of a dictionary mapping `Splitter` names to states.
"""
function (c::FeedbackChain)(h, x)
    newh = Dict{String, Any}()
    for layer ∈ c.layers
        if layer isa Splitter
            newh[splitname(layer)] = x
        elseif layer isa AbstractMerger
            x = layer(x, h)
        else
            x = layer(x)
        end
    end
    return newh, x
end # function (c::FeedbackChain)

# These overloads ensure that a FeedbackChain behaves as Flux expects, e.g.,
# when moving to gpu or collecting parameters.
functor(f::FeedbackChain) = f.layers, ls -> FeedbackChain(ls...)
#children(c::FeedbackChain) = c.layers
#mapchildren(f, c::FeedbackChain) = FeedbackChain(f.(c.layers)...)

# These overloads ensure that indexing / slicing etc. work with FeedbackChains
@forward FeedbackChain.layers Base.getindex, Base.length, Base.first, Base.last,
         Base.iterate, Base.lastindex
getindex(c::FeedbackChain, i::AbstractArray) = FeedbackChain(c.layers[i]...)

function show(io::IO, c::FeedbackChain)
    print(io, "FeedbackChain(")
    join(io, c.layers, ", ")
    print(io, ")")
end # function show
end # module FeedbackChains


text:
@testset "Without Pluto" begin
    using PlutoHooks

    ref = @use_ref(1)
    @test ref[] == 1

    x = 2
    @use_effect([x]) do
        ref[] = x
        () -> (ref[] = 9999)
    end
    @test ref[] == 2
    # cleanup never called without pluto ✓
    @test ref[] != 9999

    state, setstate = @use_state(5)
    @test state == 5
    @test_nowarn setstate(99)
    # setstate does nothing without pluto ✓
    @test state == 5

    y = 7
    result = @use_deps([y]) do
        ref2 = @use_ref(1)
        ref2[] = y
    end
    @test result == 7

    result = @use_memo([]) do
      result * y
    end

    @test result == 49
end