text:
# # Nonnegative over a variety

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Getting started/circle.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Getting started/circle.ipynb)

# The polynomial ``1 - y^2`` is nonnegative for all ``y`` in the unit circle.
# This can be verified using Sum-of-Squares.

using Test #src
using DynamicPolynomials
using SumOfSquares
@polyvar x y
S = @set x^2 + y^2 == 1

# We need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v0.23.0/installation/#Supported-solvers) for a list of the available choices.
# The domain over which the nonnegativity of ``1 - y^2`` should be certified
# is specified through the `domain` keyword argument.

import CSDP
model = SOSModel(CSDP.Optimizer)
set_silent(model)
con_ref = @constraint(model, 1 - y^2 >= 0, domain = S)
optimize!(model)

# We can see that the model was feasible:

@test JuMP.termination_status(model) == MOI.OPTIMAL #src
@test JuMP.primal_status(model) == MOI.FEASIBLE_POINT #src
solution_summary(model)

# The certificate can be obtained as follows:

dec = sos_decomposition(con_ref, 1e-6) #src
@test length(dec) == 1 #src
@test first(dec) ≈ x rtol = 1e-6 #src
sos_decomposition(con_ref, 1e-6)

# It returns ``x^2`` which is a valid certificate as:
# $$ 1 - y^2 \equiv x^2 \pmod{x^2 + y^2 - 1} $$

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
function init_observable_mappings!(m::SmetsWouters)

    observables = OrderedDict{Symbol,Observable}()
    population_mnemonic = get(get_setting(m, :population_mnemonic))

    ############################################################################
    ## 1. GDP
    ############################################################################
    gdp_fwd_transform = function (levels)
        # FROM: Level of nominal GDP (FRED :GDP series)
        # TO:   Quarter-to-quarter percent change of real, per-capita GDP, adjusted for population smoothing

        levels[!,:temp] = percapita(m, :GDP, levels)
        gdp = 1000 * nominal_to_real(:temp, levels)
        oneqtrpctchange(gdp)
    end

    gdp_rev_transform = loggrowthtopct_annualized_percapita

    observables[:obs_gdp] = Observable(:obs_gdp, [:GDP__FRED, population_mnemonic, :GDPDEF__FRED],
                                       gdp_fwd_transform, gdp_rev_transform,
                                       "Real GDP Growth", "Real GDP Growth Per Capita")

    ############################################################################
    ## 2. Hours per-capita
    ############################################################################

    hrs_fwd_transform = function (levels)
        # FROM: Average weekly hours (AWHNONAG) & civilian employment (CE16OV)
        # TO:   log (3 * per-capita weekly hours / 100)
        # Note: Not sure why the 3 is there.

        levels[!,:temp] = levels[!,:AWHNONAG] .* levels[!,:CE16OV]
        weeklyhours = percapita(m, :temp, levels)
        100*log.(3 * weeklyhours / 100)
    end

    hrs_rev_transform = logleveltopct_annualized_percapita

    observables[:obs_hours] = Observable(:obs_hours, [:AWHNONAG__FRED, :CE16OV__FRED],
                                         hrs_fwd_transform, hrs_rev_transform,
                                         "Hours Per Capita", "Log Hours Per Capita")

    ############################################################################
    ## 3. Wages
    ############################################################################

    wages_fwd_transform = function (levels)
        # FROM: Nominal compensation per hour (:COMPNFB from FRED)
        # TO: quarter to quarter percent change of real compensation (using GDP deflator)

        oneqtrpctchange(nominal_to_real(:COMPNFB, levels))
    end

    wages_rev_transform = loggrowthtopct_annualized

    observables[:obs_wages] = Observable(:obs_wages, [:COMPNFB__FRED, :GDPDEF__FRED],
                                         wages_fwd_transform, wages_rev_transform,
                                         "Percent Change in Wages",
                                         "Q-to-Q Percent Change of Real Compensation (using GDP deflator)")

    ############################################################################
    ## 4. GDP Deflator
    ############################################################################

    gdpdeflator_fwd_transform = function (levels)
        # FROM: GDP deflator (index)
        # TO:   Approximate quarter-to-quarter percent change of gdp deflator,
        #       i.e.  quarterly gdp deflator inflation

        oneqtrpctchange(levels[!,:GDPDEF])
    end

    gdpdeflator_rev_transform = loggrowthtopct_annualized

    observables[:obs_gdpdeflator] = Observable(:obs_gdpdeflator, [:GDPDEF__FRED],
                                               gdpdeflator_fwd_transform, gdpdeflator_rev_transform,
                                               "GDP Deflator",
                                               "Q-to-Q Percent Change of GDP Deflator")

    ############################################################################
    ## 5. Nominal short-term interest rate (3 months)
    ############################################################################

    nominalrate_fwd_transform = function (levels)
        # FROM: Nominal effective federal funds rate (aggregate daily data at a
        #       quarterly frequency at an annual rate)
        # TO:   Nominal effective fed funds rate, at a quarterly rate

        annualtoquarter(levels[!,:DFF])
    end

    nominalrate_rev_transform = quartertoannual

    observables[:obs_nominalrate] = Observable(:obs_nominalrate, [:DFF__FRED],
                                               nominalrate_fwd_transform, nominalrate_rev_transform,
                                               "Nominal FFR",
                                               "Nominal Effective Fed Funds Rate")

    ############################################################################
    ## 6. Consumption
    ############################################################################

    consumption_fwd_transform = function (levels)
        # FROM: Nominal consumption
        # TO:   Real consumption, approximate quarter-to-quarter percent change,
        #       per capita, adjusted for population filtering

        levels[!,:temp] = percapita(m, :PCE, levels)
        cons = 1000 * nominal_to_real(:temp, levels)
        oneqtrpctchange(cons)
    end

    consumption_rev_transform = loggrowthtopct_annualized_percapita

    observables[:obs_consumption] = Observable(:obs_consumption, [:PCE__FRED, population_mnemonic],
                                               consumption_fwd_transform, consumption_rev_transform,
                                               "Consumption growth per capita",
                                               "Consumption growth adjusted for population filtering")

    ############################################################################
    ## 7. Investment growth per capita
    ############################################################################

    investment_fwd_transform = function (levels)
        # FROM: Nominal investment
        # INTO: Real investment, approximate quarter-to-quarter percent change,
        #       per capita, adjusted for population filtering

        levels[!,:temp] = percapita(m, :FPI, levels)
        inv = 10000 * nominal_to_real(:temp, levels)
        oneqtrpctchange(inv)
    end

    investment_rev_transform  = loggrowthtopct_annualized_percapita

    observables[:obs_investment] = Observable(:obs_investment, [:FPI__FRED, population_mnemonic],
                                              investment_fwd_transform, investment_rev_transform,
                                              "Real Investment per capita",
                                              "Real investment per capita, adjusted for population filtering")

    ############################################################################
    # 8 - (8 + n_anticipated_shocks)
    ############################################################################

    for i = 1:n_anticipated_shocks(m)
        # FROM: OIS expectations of $i-period-ahead interest rates at a quarterly rate
        # TO:   Same

        ant_fwd_transform = function (levels)
            levels[!, Symbol("ant$i")]
        end

        ant_rev_transform = quartertoannual

        observables[Symbol("obs_nominalrate$i")] = Observable(Symbol("obs_ant$i"), [Symbol("ant$(i)__OIS")],
                                                      ant_fwd_transform, ant_rev_transform,
                                                      "Anticipated Shock $i",
                                                      "$i-period ahead anticipated monetary policy shock")
    end

    m.observable_mappings = observables
end

ext: jl


text:
using nbody
using Test

@testset "nbody.jl" begin
end

@testset "utils.jl" begin
    @testset "same_size" begin
        @test  nbody.same_size([1,2,3], [2,3,4], [5,6,7])
        @test !nbody.same_size([1,2,3], [2,3,4], [5,6])
    end
end

@testset "bodies.jl" begin
    @test_throws DimensionMismatch Bodies(
        [1.0, 2.0],
        ([1.0, 2.0], [1.0, 2.0], [1.0, 2.0]),
        ([1.0, 2.0], [1.0, 2.0], [1.0, 2.0]),
        ([1.0, 2.0], [1.0, 2.0], [1.0]))
end


text:
"""
The semagrams logo, rendered using our homegrown SVG library.
"""
module Logo

using ..SVG
export LOGO

julia_blue = "#4063d8"
julia_green = "#389826"
julia_red = "#cb3c33"
julia_purple = "#9558b2"

circles = [[-1,0], [0,1], [1,0], [0,-1]]
wires = [(1,1,julia_green),(1,-1,julia_red),(-1,-1,julia_purple),(-1,1,julia_blue)]
stroke_width = 0.11
r = 0.5
portr = 0.16

scale=50

offset = 1 + r + portr

function make_box(cx,cy)
  @svg(circle(r=r,cx=cx,cy=cy,style="stroke-width:$stroke_width;fill:none;stroke:black"))
end

function make_wire(x,y,color)
  src = [x, y * r]
  towards = src .+ [0, y * 0.5r]
  tgt = [x * r, y]
  from = tgt .+ [x * 0.5r, 0]
  @svg g() do
    circle(r=portr,cx=src[1],cy=src[2],style="fill:$color;stroke:none")
    circle(r=portr,cx=tgt[1],cy=tgt[2],style="fill:$color;stroke:none")
    path(d="M $(svgpointstring([src';])) C $(svgpointstring([towards';from';tgt']))",
         style="stroke-width:$stroke_width;stroke:$color;fill:none")
  end
end

LOGO = @svg svg(version="1.1", baseProfile="full",
                width=2 * offset * scale, height=2 * offset * scale,
                xmlns="http://www.w3.org/2000/svg") do
  g(transform="scale($scale $scale) translate($offset $offset)") do
    $([make_box(c...) for c in circles]...)
    $([make_wire(w...) for w in wires]...)
  end
end

end


text:
using Random: rand

mutable struct RandomModel <: ALifeModel
    num_organisms
    num_genotypes
    current_organisms::Vector{Organism}
    current_time

    RandomModel(num_organisms, num_genotypes) = new(num_organisms, num_genotypes, [], 0)
end

struct IntGenotype <: Genotype
    genotype::Int
end

function init!(model::RandomModel)
    model.current_organisms = _get_random_organism(model)
end

function iterate!(model::RandomModel)
    current_organisms = model.current_organisms
    model.current_organisms = _get_random_organism(model)
    model.current_time += 1
    return current_organisms
end

function _get_random_organism(model::RandomModel)
    new_genotypes = rand(1:model.num_genotypes, model.num_organisms)

    if model.current_organisms == []
        return [Organism(IntGenotype(new), Ancestor(), model.current_time, model.current_time+1)
                    for new in new_genotypes]
    else
        return [Organism(IntGenotype(new), old, model.current_time, model.current_time+1)
                    for (new, old) in zip(new_genotypes, model.current_organisms)]
    end
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
"""
Creates an equation matrix in which rows correspond to atoms
and columns correspond to compounds.

# Examples
```jldoctest
julia> equationmatrix(ce"H2 + Cl2 → HCl")
2×3 Array{Int64,2}:
 2  0  1
 0  2  1
```
"""
function equationmatrix(equation::ChemEquation{T}) where T
    vect = elements(equation)
    charged = hascharge(equation)
    if charged && "e" ∉ vect
        push!(vect, "")
    end

    mat = zeros(T, length(vect), length(equation.tuples))
    for (j, compoundtuple) ∈ enumerate(equation.tuples)
        compound = compoundtuple[1]
        for (element, k) ∈ compound.tuples
            i = findfirst(isequal(element), vect)
            mat[i, j] = k
        end
        if charged && hascharge(compound)
            mat[end, j] = compound.charge
        end
    end

    return mat
end

"Wrapper function for [`balancematrix`](@ref)."
function _balancematrix(equation)
    -nullspacex(equationmatrix(equation))
end

"""
Balances an equation matrix using the *nullspace method*.
Returns an array in which each column represents a solution.

# References
- [Thorne (2009)](https://arxiv.org/ftp/arxiv/papers/1110/1110.4321.pdf)
"""
balancematrix(equation::ChemEquation) = _balancematrix(equation)

"""
Same as [`balancematrix(::ChemEquation)`](@ref),
but for a chemical equation with integer coefficients.

By default, the solutions of integer matrices are displayed as integers.
If `fractions` is true, they will be displayed as rational fractions instead.
"""
function balancematrix(equation::ChemEquation{T}; fractions=false) where T<:IntegerX
    mat = _balancematrix(equation)
    if !fractions
        mat ./= gcd(mat)
        T.(mat)
    end
    return mat
end

"""
Balances the coefficients of a chemical equation.
If the equation cannot be balanced, an error is thrown.

!!! info
    The original equation is not modified.

For equations with integer coefficients, the minimal integer solution is default.
If `fractions` is true, it will be displayed as rational fractions instead.

# Examples
```jldoctest
julia> balance(ce"Cr2O7{-2} + H{+} + {-} = Cr{+3} + H2O")
Cr2O7{-2} + 14 H{+} + 6 e = 2 Cr{+3} + 7 H2O

julia> balance(ce"Fe + Cl2 = FeCl3", fractions=true)
Fe + 3//2 Cl2 = FeCl3

julia> balance(ChemEquation{Rational}("H2 + Cl2 = HCl"))
1//2 H2 + 1//2 Cl2 = HCl
```
"""
function balance(equation::ChemEquation; fractions=false)
    if fractions
        eq = ChemEquation{Rational}(equation.tuples)
        mat = balancematrix(equation, fractions=true)
    else
        eq = ChemEquation(equation.tuples)
        mat = balancematrix(equation)
    end
    num = size(mat, 2)
    if num == 1
        for (i, k) ∈ enumerate(mat)
            eq.tuples[i] = (eq.tuples[i][1], k)
        end
    elseif num > 1
        error("Chemical equation $equation can be balanced in infinitely many ways")
    else
        error("Chemical equation $equation cannot be balanced")
    end
    return eq
end


text:
using EcologicalNetworks
using StatsPlots
using Plots
using Statistics

# U = UnipartiteNetwork(rand(Bool, (10,10)))
#U = nz_stream_foodweb()[5]

null_sample = rand(null2(B_A), 10_000)
 simplify!.(null_sample)

is_a_valid_sample(x, N) = (richness(x) == richness(N)) && (links(x) == links(N))
is_a_valid_sample_of_U(x) = is_a_valid_sample(x, B_A)
# filter!(is_a_valid_sample_of_U, null_sample)

expected_spectral_radius = ρ.(null_sample)
# expected_mod = Q(brim(lp(null_sample...))...)

null_graph = density!(
    (ρ(B_A) .- expected_spectral_radius)./std(expected_spectral_radius),
    frame=:zerolines)
plot!(title = "Modèle Null Type 2 d'Emboîtement de l'Archipel", xlabel = "x", ylabel = "f(x)", leg = false)

png("null.png")

# density!(
#     ( Q_A .- expected_mod)./std(expected_mod),
#     frame=:zerolines)

max_stars_repo_name: kiristern/reseaux

max_issues_repo_name: kiristern/reseaux

max_forks_repo_name: kiristern/reseaux


text:
#
# Check if structure is in mmCIF format and returns the list of
# desired fields if that is the case
#
function check_mmCIF(file)

  local ifield, mmCIF, line

  mmCIF_fields = Indexes_mmCIF_fields()

  pdbfile = open(file,"r")

  mmCIF = false

  while ! eof(pdbfile)
    try
      line = strip(readline(pdbfile))
      if length(line) < 5
        continue
      end
    catch
      break
    end
    if line[1:5] == "loop_"
      mmCIF = true
      ifield = 0
      while ! eof(pdbfile)
        try
          line = strip(readline(pdbfile))
        catch
          break
        end
        if line[1:1] == "#"
          continue
        end
        if ifield > 0
          ifield = ifield + 1
          if occursin("_atom_site.label_atom_id",line) ; mmCIF_fields.name = ifield ; end
          if occursin("_atom_site.id",line) ; mmCIF_fields.index = ifield ; end
          if occursin("_atom_site.label_comp_id",line) ; mmCIF_fields.resname = ifield ; end
          if occursin("_atom_site.label_asym_id",line) ; mmCIF_fields.chain = ifield ; end
          if occursin("_atom_site.label_seq_id",line) ; mmCIF_fields.resnum = ifield ; end
          if occursin("_atom_site.Cartn_x",line) ; mmCIF_fields.x = ifield ; end
          if occursin("_atom_site.Cartn_y",line) ; mmCIF_fields.y = ifield ; end
          if occursin("_atom_site.Cartn_z",line) ; mmCIF_fields.z = ifield ; end
          if occursin("_atom_site.B_iso_or_equiv",line) ; mmCIF_fields.beta = ifield ; end
          if occursin("_atom_site.occupancy",line) ; mmCIF_fields.occup = ifield ; end
          if occursin("ATOM",line) && occursin("HETATOM",line)
            break
          end
        end
        if occursin("_atom_site.group_PDB",line)
          ifield = 1
        end
      end
    end
  end
  close(pdbfile)

  if mmCIF
    if mmCIF_fields.name == 0 ||
       mmCIF_fields.resname == 0 ||
       mmCIF_fields.chain == 0 ||
       mmCIF_fields.resnum == 0 ||
       mmCIF_fields.x == 0 ||
       mmCIF_fields.y == 0 ||
       mmCIF_fields.z == 0
      error(" ERROR: Could not find all necessary fields in mmCIF file. ")
    end
  end

  return mmCIF, mmCIF_fields

end


text:
module TransferEntropy
include("core.jl")
include("mutualinfo/interface.jl")
include("transferentropy/interface.jl")
end # module


text:
commands = readlines("day14/input.txt")

function decode_addresses(mask, value)
    value_bits = bitstring(value)
    out_bits = []
    for i in 0:35
        if mask[end-i] == '0'
            pushfirst!(out_bits, value_bits[end-i])
        else
            pushfirst!(out_bits, mask[end-i])
        end
    end
    out_bitstring = join(out_bits)
    addresses_out = [out_bitstring]
    while contains(addresses_out[1], 'X')
        address = popfirst!(addresses_out)
        idx = findfirst('X', address)
        push!(addresses_out, address[1:idx-1]*"0"*address[idx+1:end], address[1:idx-1]*"1"*address[idx+1:end])
    end
    return addresses_out
end

current_mask = ""
memory = Dict{Int, Int}()
for instruction in commands
    if contains(instruction, "[")
        address, value = match(r"mem\[(\d+)\] = (\d+)", instruction).captures
        for address in decode_addresses(current_mask, parse(Int, address))
            global memory[parse(Int, address; base=2)] = parse(Int, value)
        end
    else
        global current_mask = match(r"mask = (\w+)", instruction).captures[1]
    end
end

println(sum(v for v in values(memory)))


text:
#Read in the file and run latentgmm
#Lanfeng Pan, Oct 21, 2015

#change the path to the new folder if there is new version of data
datapath = Pkg.dir("LatentGaussianMixtureModel", "examples")

#adding all available cpu cores, utilizing the parallel computing
addprocs(2)

#These lines loads all the functions
Pkg.update()
import LatentGaussianMixtureModel
println("The loaded package version is", LatentGaussianMixtureModel.PKGVERSION)

#LatentGaussianMixtureModel is our package
#import LatentGaussianMixtureModel
import Distributions, StatsBase, GaussianMixtureTest
@everywhere using Distributions, StatsBase, GaussianMixtureTest

#Read in the patients covariates X

X = readcsv(joinpath(datapath, "X.csv"));
groupindex_raw = readcsv(joinpath(datapath, "groupindex.csv"));
Y_raw = readcsv(joinpath(datapath, "Y.csv"));

N, J = size(X)

#read in the groupindex for transplant center and convert to integer vector
levelsdictionary = levelsmap(groupindex_raw);
groupindex  = ones(Int64, N);
for i in 1:N
    groupindex[i] = levelsdictionary[groupindex_raw[i]]
end
#after the convert, maximum value in groupindex should be the number of transplant centers
n = maximum(groupindex);

# read in the binary response Y and convert to Bool Vector
#Not sure what Y looks like, so I list several possibilities here.

Y = Array{Bool}(N);
for i in 1:N
    if Y_raw[i] in ["Yes","YES","y", "Y", 1, "1", 1.0, "1.0", "true", "TRUE", "True"]
        Y[i] = true
    else
        Y[i] = false
    end
end

#Test if X, Y, groupindex are in the same length
@assert length(Y) == length(groupindex) == N

X = X .- mean(X, 1);

####-----------------------------
## ** Part One: Hypothesis Test

#The recommended ntrials=25. Setting it to some smaller number can save much time if we just want to test if the code is working.
#If there is error, set debuginfo=true to see more information

#set the number of components We want to test
# Null hypothesis: C = 1
lr1=LatentGaussianMixtureModel.EMtest(X, Y, groupindex, 1, ntrials=5, vtau=[.5, .3, .1;], debuginfo=true)

#If the reject C=1, further test C=2
lr2=LatentGaussianMixtureModel.EMtest(X, Y, groupindex, 2, ntrials=5, debuginfo=true)

####-------------------------------------------------
## **Part Two: Model Fitting**

#After decide the number of components, try model fitting
C=2

## Fitting
wi, mu, sigmas, betas, ml_C = LatentGaussianMixtureModel.latentgmm(X, Y, groupindex, C; maxiteration=1000, an=1/n, debuginfo=false, tol=.001)

###------------
##Only if we need to try multiple initial values
wi, mu, sigmas, betas, ml_C = LatentGaussianMixtureModel.latentgmmrepeat(X, Y, groupindex, C, ntrials=5)
###-----------------

## If you want to specify your own initial values
#[0.7, .3] is the weight pi
#[-1.0, -0.1] is the means mu
#[0.4, 0.3] is the standard deviation sigma
wi, mu, sigmas, betas, ml_C = LatentGaussianMixtureModel.latentgmm(X, Y, groupindex, 2, ones(J), [0.7, 0.3], [-1.0, -0.1], [0.4, 0.3]; maxiteration=1000, an=1/n, debuginfo=false, tol=.001)

# Print the predicted gamma
gammaprediction = LatentGaussianMixtureModel.predictgamma(X, Y, groupindex, wi, mu, sigmas, betas);
#writecsv("gammaprediction.csv", gammaprediction)
p=LatentGaussianMixtureModel.coefpvalue(X, Y, groupindex, wi, mu, sigmas, betas)

####-------------------------------------------------
## Part Three: False Dicovery Rate
#
CNull = findmax(wi)[2]
clFDR, rejectid = LatentGaussianMixtureModel.FDR(X, Y, groupindex, wi, mu, sigmas, betas, [CNull;], alphalevel=0.05)
println("The rejected transplant centers are:", rejectid)
println("Their probabiity of belong to majority is:", round(clFDR[rejectid], 4))

####--------------------------------------
## How to save the current work

using JLD
save(joinpath(datapath, "saveall_$(LatentGaussianMixtureModel.PKGVERSION).jld"), "wi", wi, "mu", mu, "sigmas", sigmas, "betas", betas, "X", X, "Y", Y ,"groupindex", groupindex, "lr1", lr1, "lr2", lr2, "gammaprediction", gammaprediction, "clFDR", clFDR, "rejectid", rejectid, "PKGVERSION", LatentGaussianMixtureModel.PKGVERSION)

##Warning! Please load all the packages first before laod the jld file.
#Or some of the data may not be able to recover.
using JLD
import Distributions, StatsBase
@everywhere using Distributions, StatsBase
@load joinpath(datapath, "saveall_$(LatentGaussianMixtureModel.PKGVERSION).jld")

############----------------------
## Plot a graph
#Pkg.add("KernelEstimator")
using RCall, KernelEstimator
mhat = MixtureModel(map((u, v) -> Normal(u, v), mu, sigmas), wi)
xs = linspace(minimum(gammaprediction)-0.5, maximum(gammaprediction)+1., 400);
denhat = pdf.(mhat, xs);
denpredict = kerneldensity(gammaprediction, xeval=xs)
den1 = probs(mhat)[1].*pdf.(mhat.components[1], xs);
den2 =  probs(mhat)[2] .* pdf.(mhat.components[2], xs);

##Send the data to R
@rput gammaprediction xs denhat den1 den2 denpredict

# regular R code inside the triple quotes.
rprint("""
#pdf("denhat_mpm.pdf", width=10, height=6.18)
plot(xs, denhat, lwd=4, type="l", xlab=expression(gamma), ylab="")
lines(xs, den1, lwd=3, lty=2, col="blue")
lines(xs, den2, lwd=3, lty=2, col="red")
#rug(gammaprediction)
#lines(xs, denpredict, lty=4, lwd=3, col="blue")
legend("topright",c("Estimated Density", "Component 1", "Component 2"), lty=c(1,2,2), lwd=c(4, 3, 3), col=c("black", "blue", "red"))
#dev.off()
NULL
""")


text:
module FermiDiracIntegrals

using Polylogarithms

export F

"""
Complete Fermi-Dirac-integral

Arguments:
* j
* x

Formula:
``{F_j(x) = \\frac{1}{\\Gamma(j+1)} \\int_0^{\\infty}{\\frac{t^j}{\\exp(t-x)+1}dt}}``

Implementation:
Using the polylogarithm
"""
function F(j,x)
    -polylog(j+1,-exp(x))
end

"""
Approximation of the complete Fermi-Dirac-integral for j = 1/2

Checked for a relative tolerance of 3% in the range x = -100:0.1:100

Speed: 100 times faster than the polylog version

Source:
J. S. Blakemore: Approximations for Fermi-Dirac Integrals. Solid-State Electronics, 25(11):1067-1076, 1982.
"""
function F(::Val{1/2},x)
    if x < 1.3
        1/(exp(-x)+0.27)
    else
        4/3/sqrt(pi)*(x^2+pi^2/6)^(3/4)
    end
end

end


text:
module TestRastrigin
    using Evolutionary
    using Base.Test

    function test_result(result::Vector, fitness::Float64, N::Int, tol::Float64)
        if round(fitness) == 0
            @test_approx_eq_eps result zeros(N) tol
            @test_approx_eq_eps fitness 0.0 tol
        else
            warn("Found local minimum!!!")
            @test sum(round(abs(result))) < N
        end
    end

    # Objective function
    function rastrigin(x::Vector{Float64})
        n = length(x)
        return 10n + sum([ x[i]^2 - 10cos(2π*x[i]) for i in 1:n ])
    end

    # Parameters
    N = 3
    P = 100

    # Testing: (μ/μ_I,λ)-σ-Self-Adaptation-ES
    # with non-isotropic mutation operator y' := y + (σ_1 N_1(0, 1), ..., σ_N N_N(0, 1))
    result, fitness, cnt = es( rastrigin, N;
            initStrategy = strategy(σ = .5ones(N), τ = 1/sqrt(2*N), τ0 = 1/sqrt(N)),
            recombination = average, srecombination = averageSigmaN,
            mutation = anisotropic, smutation = anisotropicSigma,
            selection=:comma,
            μ = 15, λ = P,
            iterations = 1000)
    println("(15/15,$(P))-σ-SA-ES-AS => F: $(fitness), C: $(cnt), OBJ: $(result)")
    test_result(result, fitness, N, 1e-1)

    # Testing: CMA-ES
    result, fitness, cnt = cmaes( rastrigin, N; μ = 15, λ = P, tol = 1e-8)
    println("(15/15,$(P))-CMA-ES => F: $(fitness), C: $(cnt), OBJ: $(result)")
    test_result(result, fitness, N, 1e-1)

    # Testing: GA
    result, fitness, cnt = ga( rastrigin, N;
            populationSize = P,
            ɛ = 0.1,
            selection = sus,
            crossover = intermediate(0.25),
            mutation = domainrange(fill(1.0,N)),
            tol = 1e-5)
    println("GA(p=$(P),x=.8,μ=.1,ɛ=0.1) => F: $(fitness), C: $(cnt), OBJ: $(result)")
    test_result(result, fitness, N, 1e-1)

end


text:
#
# common.jl -
#
# Common functions for testing.
#

"""

```julia
floating_point_type(A, B, ...)
```

yields the floating-point type for operations between arrays `A`, `B`, ...

"""
floating_point_type(args::AbstractArray...) =
    float(real(promote_type(map(eltype, args)...)))

"""

```julia
relative_precision(A, B, ...)
```

yields the worst of the relative precisions of the element types of arrays `A`,
`B`, ...

"""
relative_precision(args::AbstractArray...) = max(map(relative_precision, args)...)
relative_precision(A::AbstractArray{T}) where {T} = eps(float(real(T)))

"""

```julia
test_api(P, A, x, y; atol=0, rtol=sqrt(eps(relative_precision(x,y))))
```

test LazyAlgebra API for mapping `P(A)` using variables `x` and `y`.

"""
function test_api(::Type{P}, A::Mapping, x0::AbstractArray, y0::AbstractArray;
                  rtol::Real=sqrt(relative_precision(x0,y0)),
                  atol::Real=0) where {P<:Union{Direct,InverseAdjoint}}
    x = vcopy(x0)
    y = vcopy(y0)
    z = apply(P, A, x)
    @test x == x0
    T = floating_point_type(x, y, z)
    for α in (0, 1, -1,  2.71, π),
        β in (0, 1, -1, -1.33, Base.MathConstants.φ),
        scratch in (false, true)
        @test apply!(α, P, A, x, scratch, β, vcopy(y)) ≈
            T(α)*z + T(β)*y  atol=atol rtol=rtol
        if scratch
            vcopy!(x, x0)
        else
            @test x == x0
        end
    end
end

function test_api(::Type{P}, A::Mapping, x0::AbstractArray, y0::AbstractArray;
                  rtol::Real=sqrt(relative_precision(x0,y0)),
                  atol::Real=0) where {P<:Union{Adjoint,Inverse}}
    x = vcopy(x0)
    y = vcopy(y0)
    z = apply(P, A, y)
    @test y == y0
    T = floating_point_type(x, y, z)
    for α in (0, 1, -1,  2.71, π),
        β in (0, 1, -1, -1.33, Base.MathConstants.φ),
        scratch in (false, true)
        @test apply!(α, P, A, y, scratch, β, vcopy(x)) ≈
            T(α)*z + T(β)*x  atol=atol rtol=rtol
        if scratch
            vcopy!(y, y0)
        else
            @test y == y0
        end
    end
end


text:
using ODE, Test, LinearAlgebra

tols = [5e-2, 1e-2, 1e-2]

solvers = [
           ## Non-stiff
           # fixed step
           ODE.ode1,
           ODE.ode2_midpoint,
           ODE.ode2_heun,
           ODE.ode4,
           ODE.ode4ms,
           ODE.ode5ms,
           # adaptive
#           ODE.ode21, # this fails on Travis with 0.4?! TODO revert once fixed.
           ODE.ode23,
           ODE.ode45_dp,
           ODE.ode45_fe,
           ODE.ode78,

           ## Stiff
           # fixed-step
           ODE.ode4s_s,
           ODE.ode4s_kr,
           # adaptive
           ODE.ode23s]

dtypes = [Float32, Float64, BigFloat]
for solver in solvers
    println("using $solver")
    for (tol,T) in zip(tols, dtypes)
        if solver==ODE.ode45_fe && T==Float32
            # For some reason ode45_fe hangs with Float32!
            continue
        end

        # dy
        # -- = 6 ==> y = 6t
        # dt
        t,y=solver((t,y)->T(6), T(0), T[0:.1:1;])
        @test maximum(abs.(y-6t)) < tol
        @test eltype(t)==T

        # dy
        # -- = 2t ==> y = t.^2
        # dt
        t,y=solver((t,y)->T(2)*t, T(0), T[0:.001:1;])
        @test maximum(abs.(y-t.^2)) < tol
        @test eltype(t)==T

        # dy
        # -- = y ==> y = y0*e.^t
        # dt
        t,y=solver((t,y)->y, T(1), T[0:.001:1;])
        @test maximum(abs.(y-ℯ.^t)) < tol
        @test eltype(t)==T

        t,y=solver((t,y)->y, T(1), T[1:-.001:0;])
        @test maximum(abs.(y .- ℯ .^ (t.-1))) < tol
        @test eltype(t)==T

        # dv       dw
        # -- = -w, -- = v ==> v = v0*cos(t) - w0*sin(t), w = w0*cos(t) + v0*sin(t)
        # dt       dt
        #
        # y = [v, w]
        t,y=solver((t,y)->[-y[2]; y[1]], T[1, 2], T[0:.001:2*pi;])
        ys = hcat(y...)'   # convert Vector{Vector{T}} to Matrix{T}
        @test maximum(abs.(ys-[cos.(t)-2*sin.(t) 2*cos.(t)+sin.(t)])) < tol
        @test eltype(t)==T
    end
end

# Test negative starting times ODE.ode23s
@test length(ODE.ode23s((t,y)->[-y[2]; y[1]], [1., 2.], [-5., 0])[1]) > 1

# rober testcase from http://www.unige.ch/~hairer/testset/testset.html
let
    println("ROBER test case")
    function f(t, y)
        ydot = similar(y)
        ydot[1] = -0.04*y[1] + 1.0e4*y[2]*y[3]
        ydot[3] = 3.0e7*y[2]*y[2]
        ydot[2] = -ydot[1] - ydot[3]
        ydot
    end
    t = [0., 1e11]
    t,y = ode23s(f, [1.0, 0.0, 0.0], t; abstol=1e-8, reltol=1e-8,
                                        maxstep=1e11/10, minstep=1e11/1e18)

    refsol = [0.2083340149701255e-07,
              0.8333360770334713e-13,
              0.9999999791665050] # reference solution at tspan[2]
    @test LinearAlgebra.norm(refsol-y[end], Inf) < 2e-10
end
include("interface-tests.jl")
include("common.jl")
println("All looks OK")


text:
module SeqSpace

using GZip
using BSON: @save
using LinearAlgebra: norm, svd, Diagonal
using Statistics: quantile, std
using Flux, Zygote
using ProgressMeter

import BSON

include("io.jl")
include("rank.jl")
include("model.jl")
include("voronoi.jl")
include("pointcloud.jl")
include("generate.jl")
include("normalize.jl")
include("manifold.jl")
include("infer.jl")
include("scrna.jl")

using .PointCloud, .DataIO, .SoftRank, .ML

export Result, HyperParams
export linearprojection, fitmodel, extendfit
export marshal, unmarshal

# ------------------------------------------------------------------------
# globals

# ------------------------------------------------------------------------
# types

"""
    mutable struct HyperParams
        dₒ :: Int
        Ws :: Array{Int,1}
        BN :: Array{Int,1}
        DO :: Array{Int,1}
        N  :: Int
        δ  :: Int
        η  :: Float64
        B  :: Int
        V  :: Int
        k  :: Int
        γₓ :: Float32
        γᵤ :: Float32
        g  :: Function
    end

HyperParams is a collection of parameters that specify the network architecture and training hyperparameters of the autoencoder.
`dₒ` is the desired output dimensionality of the encoding layer
`Ws` is a collection of the network layer widths. The number of entries controls the depth. The decoder is mirror-symmetric.
`BN` is the collection of layers that will be _followed_ by batch normalization.
`DO` is the collection of layers that will be _followed_ by dropout.
`N`  is the number of epochs to train against
`δ`  is the number of epochs that will be sampled for logging
`η`  is the learning rate
`B`  is the batch size
`V`  is the number of points to partition for validation purposes, i.e. won't be training against
`k`  is the number of neighbors to use to estimate geodesics
`γₓ` is the prefactor of distance soft rank loss
`γᵤ` is the prefactor of uniform density loss
`g ` is the metric given to latent space
"""
mutable struct HyperParams
    dₒ :: Int          # output dimensionality
    Ws :: Array{Int,1} # (latent) layer widths
    BN :: Array{Int,1} # (latent) layers followed by batch normalization
    DO :: Array{Int,1} # (latent) layers followed by drop outs
    N  :: Int          # number of epochs to run
    δ  :: Int          # epoch subsample factor for logging
    η  :: Float64      # learning rate
    B  :: Int          # batch size
    V  :: Int          # number of points to partition for validation
    k  :: Int          # number of neighbors to use to estimate geodesics
    γₓ :: Float32      # prefactor of distance soft rank loss
    γᵤ :: Float32      # prefactor of uniform density loss
    g  :: Function     # metric given to latent space
end

"""
    euclidean²(x)

Generate the matrix of pairwise distances between vectors `x`, assuming the Euclidean metric.
`x` is assumed to be ``d \times N`` where ``d`` denotes the dimensionality of the vector and ``N`` denotes the number.
"""
const euclidean²(x) = sum( (x[d,:]' .- x[d,:]).^2 for d in 1:size(x,1) )

"""
    cylinders²(x)

Generate the matrix of pairwise distances between vectors `x`, assuming the points are distributed on a cylinder.
`x` is assumed to be ``d \times N`` where ``d`` denotes the dimensionality of the vector and ``N`` denotes the number.
The first coordinate of `x` is assumed to be the polar coordinate.
"""
function cylinder²(x)
    c = cos.(π.*(x[1,:]))
    s = sin.(π.*(x[1,:]))

    return (c' .- c).^2 .+ (s' .- s).^2 .+ euclidean²(x[2:end,:])
end

HyperParams(; dₒ=2, Ws=Int[], BN=Int[], DO=Int[], N=200, δ=10, η=1e-3, B=64, V=1, k=12, γₓ=1, γᵤ=1e-1, g=euclidean²) = HyperParams(dₒ, Ws, BN, DO, N, δ, η, B, V, k, γₓ, γᵤ, g)

"""
    struct Result
        param :: HyperParams
        loss  :: NamedTuple{(:train, :valid), Tuple{Array{Float64,1},Array{Float64,1}} }
        model
    end

Store the output of a trained autoencoder.
`param` stores the input hyperparameters used to design and train the neural network.
`loss` traces the dynamics of the optimization found during training.
`model` represents the learned pullback and pushforward functions.
"""
struct Result
    param :: HyperParams
    loss  :: NamedTuple{(:train, :valid), Tuple{Array{Float64,1},Array{Float64,1}} }
    model
end

"""
    marshal(r::Result)

Serialize a trained autoencoder to binary format suitable for disk storage.
Store parameters of model as contiguous array
"""
function marshal(r::Result)
    # trainable parameters
    ω₁ = r.model.pullback    |> cpu |> Flux.params |> collect
    ω₂ = r.model.pushforward |> cpu |> Flux.params |> collect

    # untrainables
    β₁ = [ (μ=layer.μ,σ²=layer.σ²) for layer in r.model.pullback.layers    if isa(layer,Flux.BatchNorm) ]
    β₂ = [ (μ=layer.μ,σ²=layer.σ²) for layer in r.model.pushforward.layers if isa(layer,Flux.BatchNorm) ]

    return Result(r.param,r.loss,
            (
                pullback=(
                    params=ω₁,
                    batchs=β₁,
                ),
                pushforward=(
                    params=ω₂,
                    batchs=β₂,
                ),
                size=size(r.model.pullback.layers[1].weight,2)
            )
    )
end

"""
    unmarshal(r::Result)

Deserialize a trained autoencoder from binary format to semantic format.
Represents model as a collection of functors.
"""
function unmarshal(r)
    autoencoder = model(r.model.size, r.param.dₒ;
          Ws         = r.param.Ws,
          normalizes = r.param.BN,
          dropouts   = r.param.DO
    )

    Flux.loadparams!(autoencoder.pullback,    r.model.pullback.params)
    Flux.loadparams!(autoencoder.pushforward, r.model.pushforward.params)
    Flux.trainmode!(autoencoder.identity, false)

    i = 1
    for layer in autoencoder.pullback.layers
        if isa(layer, Flux.BatchNorm)
            layer.μ  = r.model.pullback.batchs[i].μ
            layer.σ² = r.model.pullback.batchs[i].σ²
            i += 1
        end
    end

    i = 1
    for layer in autoencoder.pushforward.layers
        if isa(layer, Flux.BatchNorm)
            layer.μ  = r.model.pushforward.batchs[i].μ
            layer.σ² = r.model.pushforward.batchs[i].σ²
            i += 1
        end
    end

    param = HyperParams(;
        dₒ = r.param.dₒ,
        Ws = r.param.Ws,
        BN = r.param.BN,
        DO = r.param.DO,
        N  = r.param.N,
        δ  = r.param.δ,
        η  = r.param.η,
        B  = r.param.B,
        V  = r.param.V,
        k  = r.param.k,
        γₓ = r.param.γₓ,
        γᵤ = r.param.γᵤ,
    )

    return Result(param, r.loss, autoencoder)
end

# ------------------------------------------------------------------------
# utility functions

# α := rescale data by
# δ := subsample data by
function pointcloud(;α=1, δ=1)
    verts, _ = open("$root/gut/mesh_apical_stab_000153.ply") do io
        read_ply(io)
    end

    return α*vcat(
        map(v->v.x, verts)',
        map(v->v.y, verts)',
        map(v->v.z, verts)'
    )[:,1:δ:end]
end

function expression(;raw=false)
    scrna, genes, _  = if raw
        GZip.open("$root/dvex/dge_raw.txt.gz") do io
            read_matrix(io; named_cols=false, named_rows=true)
        end
    else
        GZip.open("$root/dvex/dge_normalized.txt.gz") do io
            read_matrix(io; named_cols=true, named_rows=true)
        end
    end

    return scrna, genes
end

# ------------------------------------------------------------------------
# i/o

# assumes a BSON i/o
function load(io::IO)
    database = BSON.parse(io)
    result, input = database[:result], database[:in]
end

mean(x) = length(x) > 0 ? sum(x) / length(x) : 0
mean(x::Matrix;dims=1) = sum(x;dims=dims) / size(x,dims)

function cor(x, y)
    μ = (
        x=mean(x),
        y=mean(y)
    )
    var = (
       x=mean(x.^2) .- μ.x^2,
       y=mean(y.^2) .- μ.y^2
    )

    return (mean(x.*y) .- μ.x.*μ.y) / sqrt(var.x*var.y)
end

"""
    buildloss(model, D², param)

Return a loss function used to train a neural network `model` according to input hyperparameters `param`.
`model` is a object with three fields, `pullback`, `pushforward`, and `identity`.
`pullback` and `pushforward` refers to the encoder and decoder layers respectively, while the identity is the composition.
`D²` is a matrix of pairwise distances that will be used as a quenched hyperparameter in the distance soft rank loss.
"""
function buildloss(model, D², param)
    return function(x, i::T, output::Bool) where T <: AbstractArray{Int,1}
        z = model.pullback(x)
        y = model.pushforward(z)

        # reconstruction loss
        ϵᵣ = sum(sum((x.-y).^2, dims=2)) / sum(sum(x.^2,dims=2))

        # distance softranks
        Dz² = param.g(z)
        Dx² = D²[i,i]

        dx, dz = PointCloud.upper_tri(Dx²), PointCloud.upper_tri(Dz²)
        rx, rz = softrank(dx ./ mean(dx)), softrank(dz ./ mean(dz))
        ϵₓ = 1 - cor(1 .- rx, 1 .- rz)
        #=
        ϵₓ = mean(
            let
                dx, dz = Dx²[:,j], Dz²[:,j]
                rx, rz = softrank(dx ./ mean(dx)), softrank(dz ./ mean(dz))
                1 - cor((1 .- rx).^2, (1 .- rz).^2)
            end for j ∈ 1:size(Dx²,2)
        )
        =#

        ϵᵤ = mean(
            let
                zₛ = sort(z[d,:])
                mean( ( (2*i/length(zₛ)-1) - s)^2 for (i,s) in enumerate(zₛ) )
            end for d ∈ 1:size(z,1)
        )

        #=
        ϵᵤ = let
            a = Voronoi.volumes(z)
            std(a) / mean(a)
        end
        =#

        if output
            println(stderr, "ϵᵣ=$(ϵᵣ), ϵₓ=$(ϵₓ), ϵᵤ=$(ϵᵤ)")
        end

        return ϵᵣ + param.γₓ*ϵₓ + param.γᵤ*ϵᵤ
    end
end

# ------------------------------------------------------------------------
# main functions

"""
    linearprojection(x, d; Δ=1, Λ=nothing)

Project an empirical distance matrix `x` onto `d` top principal components.
Centers the result to have zero mean.
Returns the projection, as well as a function to transform a projected vector back to the embedding space.
Ignores top `Δ` principal components.
If `Λ` is not nothing, assumes it is a precomputed SVD decomposition.
"""
function linearprojection(x, d; Δ=1, Λ=nothing)
    Λ = isnothing(Λ) ? svd(x) : Λ

    ι = (1:d) .+ Δ
    ψ = Diagonal(Λ.S[ι])*Λ.Vt[ι,:]
    μ = mean(ψ, dims=2)

    x₀ = (Δ > 0) ? Λ.U[:,1:Δ]*Diagonal(Λ.S[1:Δ])*Λ.Vt[1:Δ,:] : 0

    embed(x)   = (x₀ .+ (Λ.U[:,ι]*(x.+μ)))
    embed(x,i) = (x₀ .+ (Λ.U[:,ι]*(x.+μ[i])))

    return (
        embed      = embed,
        projection = (ψ .- μ),
    )
end

"""
    fitmodel(data, param; D²=nothing, chatty=true)

Train an autoencoder model, specified with `param` hyperparams, to fit `data`.
`data` is assumed to be sized ``d \times N`` where ``d`` and ``N`` are dimensionality and cardinality respectively.
If not nothing, `D²` is assumed to be a precomputed distance matrix of point cloud `data`.
If `chatty` is true, function will print to `stdout`.
Returns a `Result` type.
"""
function fitmodel(data, param; D²=nothing, chatty=true)
    D² = isnothing(D²) ? geodesics(data, param.k).^2 : D²

    M = model(size(data,1), param.dₒ;
          Ws         = param.Ws,
          normalizes = param.BN,
          dropouts   = param.DO
    )

    nvalid = size(data,2) - ((size(data,2)÷param.B)-param.V)*param.B
    batch, index = validate(data, nvalid)

    loss = buildloss(M, D², param)
    E    = (
        train = Float64[],
        valid = Float64[]
    )

    progress = Progress(param.N; desc=">training model", output=stderr)
    log = (n) -> begin
        if (n-1) % param.δ == 0
            push!(E.train,loss(batch.train, index.train, chatty))
            push!(E.valid,loss(batch.valid, index.valid, chatty))
        end

        next!(progress)

        nothing
    end

    Flux.trainmode!(M.identity, true)
    train!(M, batch.train, index.train, loss;
        η   = param.η,
        B   = param.B,
        N   = param.N,
        log = log
    )
    Flux.trainmode!(M.identity, false)

    return Result(param, E, M), (batch=batch, index=index, D²=D², log=log)
end

"""
    extendfit(result::Result, input, epochs)

Retrain model within `result` on `input` data for `epochs` more iterations.
Returns a new `Result`.
"""
function extendfit(result::Result, input, epochs)
    loss = buildloss(result.model, input.D², result.param)
    train!(result.model, input.y.train, input.index.train, loss;
        η   = result.param.η,
        B   = result.param.B,
        N   = epochs,
        log = input.log
    )

    return result, input
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using Test

# -----------------------------------------------------------------------------------------
# Test the tess_set_debug_param(inst, name, Int) method.
@testset "tess_set_debug_param(inst, name, $type)" for type in [ UInt8, Int16, Int32, Int64, Int]
    local inst = TessInst("eng", datadir)

    @test tess_set_debug_param(inst, "matcher_debug_level", type(3)) == true
    @test tess_set_debug_param(inst, "classify_class_pruner_threshold", type(100)) == true
    @test tess_get_param(inst, "classify_class_pruner_threshold", type) == type(229)

    tess_delete!(inst)

    local err  = "Instance has been freed."
    @test (@test_logs (:error, err) tess_set_debug_param(inst, "matcher_debug_level", type(100))) == false
end

# -----------------------------------------------------------------------------------------
# Test the tess_set_debug_param(inst, name, Bool) method.
@testset "tess_set_debug_param(inst, name, Bool)" begin
    local inst = TessInst("eng", datadir)

    @test tess_set_debug_param(inst, "tessedit_display_outwords", true) == true
    @test tess_get_param(inst, "tessedit_display_outwords", Bool) == true
    @test tess_set_debug_param(inst, "tessedit_unrej_any_wd", true) == true
    @test tess_get_param(inst, "tessedit_unrej_any_wd", Bool) == false

    tess_delete!(inst)

    local err  = "Instance has been freed."
    @test (@test_logs (:error, err) tess_set_debug_param(inst, "matcher_debug_flags", true)) == false
end

# -----------------------------------------------------------------------------------------
# Test the tess_set_debug_param(inst, name, String) method.
@testset "tess_set_debug_param(inst, name, String)" begin
    local inst = TessInst("eng", datadir)

    @test tess_set_debug_param(inst, "debug_file", "test.log") == true
    @test tess_get_param(inst, "debug_file", String) == "test.log"
    @test tess_set_debug_param(inst, "page_separator", "\n* * *\n") == true
    @test tess_get_param(inst, "page_separator", String) == "\f"

    tess_delete!(inst)

    local err  = "Instance has been freed."
    @test (@test_logs (:error, err) tess_set_debug_param(inst, "debug_file", "abc")) == false
end


text:
function pointcoords(geom::AbstractGeometry)
    @assert geotype(geom) == :Point
    [tuple(coordinates(geom)...)]
end

function multipointcoords(geom::AbstractGeometry)
    @assert geotype(geom) == :MultiPoint
    coords = coordinates(geom)
    first.(coords), last.(coords)
end

function linestringcoords(geom::AbstractGeometry)
    @assert geotype(geom) == :LineString
    coords = coordinates(geom)
    first.(coords), last.(coords)
end

function multilinestringcoords(geom::AbstractGeometry)
    @assert geotype(geom) == :MultiLineString
    x, y = Float64[], Float64[]
    for line in coordinates(geom)
        append!(x, first.(line)); push!(x, NaN)
        append!(y, last.(line)); push!(y, NaN)
    end
    x, y
end

function polygoncoords(geom::AbstractGeometry)
    @assert geotype(geom) == :Polygon
    ring = first(coordinates(geom)) # currently doesn't plot holes
    first.(ring), last.(ring)
end

function multipolygoncoords(geom::AbstractGeometry)
    @assert geotype(geom) == :MultiPolygon
    x, y = Float64[], Float64[]
    for poly in coordinates(geom)
        ring = first(coordinates(poly)) # currently doesn't plot holes
        append!(x, first.(ring)); push!(x, NaN)
        append!(y, last.(ring)); push!(y, NaN)
    end
    x, y
end

shapecoords(geom::AbstractPoint) = pointcoords(geom)
shapecoords(geom::AbstractMultiPoint) = multipointcoords(geom)
shapecoords(geom::AbstractLineString) = linestringcoords(geom)
shapecoords(geom::AbstractMultiLineString) = multilinestringcoords(geom)
shapecoords(geom::AbstractPolygon) = polygoncoords(geom)
shapecoords(geom::AbstractMultiPolygon) = multipolygoncoords(geom)

function shapecoords(geom::AbstractGeometry)
    gtype = geotype(geom)
    if gtype == :Point
        return pointcoords(geom)
    elseif gtype == :MultiPoint
        return multipointcoords(geom)
    elseif gtype == :LineString
        return linestringcoords(geom)
    elseif gtype == :MultiLineString
        return multilinestringcoords(geom)
    elseif gtype == :Polygon
        return polygoncoords(geom)
    elseif gtype == :MultiPolygon
        return multipolygoncoords(geom)
    else
        warn("unknown geometry type: $gtype")
    end
end

RecipesBase.@recipe f(geom::AbstractPoint) = (
    aspect_ratio --> 1;
    seriestype --> :scatter;
    legend --> :false;
    shapecoords(geom)
)

RecipesBase.@recipe f(geom::AbstractMultiPoint) = (
    aspect_ratio --> 1;
    seriestype --> :scatter;
    legend --> :false;
    shapecoords(geom)
)

RecipesBase.@recipe f(geom::AbstractLineString) = (
    aspect_ratio --> 1;
    seriestype --> :path;
    legend --> :false;
    shapecoords(geom)
)

RecipesBase.@recipe f(geom::AbstractMultiLineString) = (
    aspect_ratio --> 1;
    seriestype --> :path;
    legend --> :false;
    shapecoords(geom)
)

RecipesBase.@recipe f(geom::AbstractPolygon) = (
    aspect_ratio --> 1;
    seriestype --> :shape;
    legend --> :false;
    shapecoords(geom)
)

RecipesBase.@recipe f(geom::AbstractMultiPolygon) = (
    aspect_ratio --> 1;
    seriestype --> :shape;
    legend --> :false;
    shapecoords(geom)
)

RecipesBase.@recipe function f(geom::AbstractGeometry)
    aspect_ratio --> 1
    legend --> :false
    gtype = geotype(geom)
    if gtype == :Point || gtype == :MultiPoint
        seriestype := :scatter
    elseif gtype == :LineString || gtype == :MultiLineString
        seriestype := :path
    elseif gtype == :Polygon || gtype == :MultiPolygon
        seriestype := :shape
    else
        @warn("unknown geometry type: $gtype")
    end
    shapecoords(geom)
end

RecipesBase.@recipe function f(geom::Vector{<:Union{Missing, AbstractGeometry}})
    aspect_ratio --> 1
    legend --> :false
    for g in skipmissing(geom)
        @series begin
            gtype = geotype(g)
            if gtype == :Point || gtype == :MultiPoint
                seriestype := :scatter
            elseif gtype == :LineString || gtype == :MultiLineString
                seriestype := :path
            elseif gtype == :Polygon || gtype == :MultiPolygon
                seriestype := :shape
            else
                @warn("unknown geometry type: $gtype")
            end
            shapecoords(g)
        end
    end
end

RecipesBase.@recipe f(feature::AbstractFeature) = geometry(feature)
RecipesBase.@recipe f(features::Vector{<:AbstractFeature}) = geometry.(features)
RecipesBase.@recipe f(collection::AbstractFeatureCollection) = features(collection)
RecipesBase.@recipe f(collection::AbstractGeometryCollection) = geometries(collection)


text:
# redis client in julia
# wraps the stable hiredis C library

module HiRedis

using Logging

logger = SimpleLogger(stdout, Logging.Warn)
global_logger(logger)

const REDIS_ERR = -1
const REDIS_OK = 0

const REDIS_REPLY_STRING = 1
const REDIS_REPLY_ARRAY = 2
const REDIS_REPLY_INTEGER = 3
const REDIS_REPLY_NIL = 4
const REDIS_REPLY_STATUS = 5
const REDIS_REPLY_ERROR = 6

global redisContext = 0
global pipelinedCommandCount = 0

struct RedisReadTask
    rtype::Int32
    elements::Int32
    idx::Int32
    obj::Ptr{Nothing}
    parent::Ptr{RedisReadTask}
    privdata::Ptr{Nothing}
end

function create_string(task::Ptr{RedisReadTask}, str::Ptr{UInt8}, len::UInt)
    # not implemented
    ret::Ptr{Nothing} = 0
    ret
end

function create_array(task::Ptr{RedisReadTask}, len::Int32)
    # not implemented
    ret::Ptr{Nothing} = 0
    ret
end

function create_integer(task::Ptr{RedisReadTask}, int::Int)
    # not implemented
    ret::Ptr{Nothing} = 0
    ret
end

function create_nil(task::Ptr{RedisReadTask})
    # not implemented
    ret::Ptr{Nothing} = 0
    ret
end

function free_object(obj::Ptr{Nothing})
    # not implemented
    ret::Nothing = 0
    ret
end

const create_string_c = @cfunction(create_string, Ptr{Nothing}, (Ptr{RedisReadTask}, Ptr{UInt8}, UInt))

const create_array_c = @cfunction(create_array, Ptr{Nothing}, (Ptr{RedisReadTask}, Int32))

const create_integer_c = @cfunction(create_integer, Ptr{Nothing}, (Ptr{RedisReadTask}, Int))

const create_nil_c = @cfunction(create_nil, Ptr{Nothing}, (Ptr{RedisReadTask},))

const free_object_c = @cfunction(free_object, Nothing, (Ptr{Nothing},))

struct RedisReplyObjectFunctions
    create_string_c
    create_array_c
    create_integer_c
    create_nil_c
    free_object_c
end

struct RedisReader
    err::Int32
    errstr::Ptr{UInt8}
    buf::Ptr{UInt8}
    pos::UInt
    len::UInt
    maxbuf::UInt
    rstack::Array{RedisReadTask,1}
    ridx::Int32
    reply::Ptr{Nothing}
    fn::Ptr{RedisReplyObjectFunctions}
    privdata::Ptr{Nothing}
end

struct RedisContext
    err::Int32
    errstr::Ptr{UInt8}
    fd::Int32
    flags::Int32
    obuf::Ptr{UInt8}
    reader::Ptr{RedisReader}
end

struct RedisReply
    rtype::Int32                  # REDIS_REPLY_*
    integer::UInt64               # The integer when type is REDIS_REPLY_INTEGER
    len::Int32                    # Length of string
    str::Ptr{UInt8}               # Used for both REDIS_REPLY_ERROR and REDIS_REPLY_STRING
    elements::UInt                # number of elements, for REDIS_REPLY_ARRAY
    element::Ptr{Ptr{RedisReply}} # elements vector for REDIS_REPLY_ARRAY
end

function start_session(host::String = "127.0.0.1", port::Int = 6379)
    global redisContext = ccall((:redisConnect, "libhiredis"), Ptr{RedisContext}, (Ptr{UInt8}, Int32), host, port)
end

function end_session()
    if redisContext != 0 # isdefined(:redisContext)
        ccall((:redisFree, "libhiredis"), Nothing, (Ptr{RedisContext},), redisContext::Ptr{RedisContext})
    end
end

"""Free memory allocated to objects returned from hiredis"""
function free_reply_object(redisReply)
    ccall((:freeReplyObject, "libhiredis"), Nothing, (Ptr{RedisReply},), redisReply)
end

"""
Appends commands to an output buffer. Pipelining is sending a batch of commands
to redis to be processed in bulk. It cuts down the number of network requests.
"""
function pipeline_command(command::String)
    if redisContext == 0 # !isdefined(:redisContext)
        start_session()
    end
    @debug(string("RedisClient.pipeline_command: ", command))
    global pipelinedCommandCount += 1
    ccall((:redisAppendCommand, "libhiredis"), Int32, (Ptr{RedisContext}, Ptr{UInt8}), redisContext::Ptr{RedisContext}, command)
end

"""
In a blocking context, this function first checks if there are unconsumed
replies to return and returns one if so. Otherwise, it flushes the output
buffer to the socket and reads until it has a reply.
"""
function call_get_reply(redisReply::Array{Ptr{RedisReply},1})
    ccall((:redisGetReply, "libhiredis"), Int32, (Ptr{RedisContext}, Ptr{Ptr{RedisReply}}), redisContext::Ptr{RedisContext}, redisReply)
end

"""
Calls call_get_reply until the pipelinedCommandCount is 0 or an error
is returned. Adds the results from each reply to an Array, then returns the
Array.
"""
function get_reply()
    redisReply = Array(Ptr{RedisReply},1)
    results = Any[]
    while pipelinedCommandCount::Int > 0 && call_get_reply(redisReply) == REDIS_OK
        push!(results, get_result(redisReply[1]))
        global pipelinedCommandCount -= 1
    end
    results
end

"""
Converts the reply object from hiredis into a String, int, or Array
as appropriate the the reply type.
"""
function get_result(redisReply::Ptr{RedisReply})
    r = unsafe_load(redisReply)
    if r.rtype == REDIS_REPLY_ERROR
        error(unsafe_string(r.str))
    end
    ret::Any = nothing
    if r.rtype == REDIS_REPLY_STRING
        ret = unsafe_string(r.str)
    elseif r.rtype == REDIS_REPLY_INTEGER
        ret = Int(r.integer)
    elseif r.rtype == REDIS_REPLY_ARRAY
        n = Int(r.elements)
        results = String[]
        replies = unsafe_wrap(Array, r.element, n)
        for i in 1:n
            ri = unsafe_load(replies[i])
            push!(results, unsafe_string(ri.str))
        end
        ret = results
    end
    free_reply_object(redisReply)
    ret
end

"""Pipelines a block of ordinary blocking calls."""
macro pipeline(expr::Expr)
    Expr(:block, map(x ->
        begin
            if x.args[1] in names(HiRedis)
                args = copy(x.args)
                push!(args, Expr(:kw, :pipeline, true))
                Expr(x.head, args...)
            else
                x
            end
        end, filter(x -> typeof(x) == Expr, expr.args))...)
end

"""Issues a blocking command to hiredis."""
function do_command(command::String)
    if redisContext == 0 # !isdefined(:redisContext)
#         error("redisContext not defined. Please call RedisClient.start_session.")
        start_session()
    end
    @debug(string("RedisClient.do_command: ", command))
    redisReply = ccall((:redisvCommand, "libhiredis"), Ptr{RedisReply}, (Ptr{RedisContext}, Ptr{UInt8}), redisContext::Ptr{RedisContext}, command)
    get_result(redisReply)
end

"""Issues a blocking command to hiredis, accepting command arguments as an Array."""
function do_command(argv::AbstractVector{<: Any})
    if redisContext == 0 # !isdefined(:redisContext)
        start_session()
    end
    redisReply = ccall((:redisCommandArgv, "libhiredis"), Ptr{RedisReply}, (Ptr{RedisContext}, Int32, Ptr{Ptr{UInt8}}, Ptr{UInt}), redisContext::Ptr{RedisContext}, length(argv), argv, C_NULL)
    get_result(redisReply)
end

"""Switches between blocking and pipelined command execution according to flag."""
function docommand(cmd::String, pipeline::Bool)
    (pipeline || (pipelinedCommandCount::Int > 0)) ? pipeline_command(cmd) : do_command(cmd)
end

# @doc "Set the string value of a key." Dict{Any,Any}(
#     :params => Dict(
#         :key           => "The key to set",
#         :value         => "The value to set",
#         :ex            => "Set the specified expire time, in seconds",
#         :px            => "Set the specified expire time, in milliseconds",
#         :pipeline      => "A flag to indicate that the command should be pipelined"
#     )) ->
"""
Set the string value of a key.

Params:

* key - The key to set
* value - The value to set
* ex - Set the specified expire time, in seconds
* px - Set the specified expire time, in milliseconds
* pipeline - A flag to indicate that the command should be pipelined
"""
function kvset(key::String, value::Any; ex::Int=0, px::Int=0, pipeline::Bool=false)
    cmd = string("SET ", key, " ", value)
    if ex > 0
        cmd = string(cmd, " EX ", ex)
    elseif px > 0
        cmd = string(cmd, " PX ", ex)
    end
    docommand(cmd, pipeline)
end

"""Get the value of a key."""
function kvget(key::String; pipeline::Bool=false)
    docommand(string("GET ", key), pipeline)
end

"""Increment the integer value of a key by one."""
function incr(key::String; pipeline::Bool=false)
    docommand(string("INCR ", key), pipeline)
end

"""Increment the integer value of a key by the given amount."""
function incrby(key::String, by::Int; pipeline::Bool=false)
    docommand(string("INCRBY ", key, " ", by), pipeline)
end

"""Delete a key."""
function del(key::String)
    docommand(string("DEL ", key))
end

"""Determine if a key exists."""
function exists(key::String)
    docommand(string("EXISTS ", key))
end

"""Find all keys matching the given pattern."""
function getkeys(pattern::String)
    docommand(string("KEYS ", pattern))
end

"""Return a serialized version of the value stored at the specified key."""
function rdump(key::String)
    docommand(string("DUMP ", key))
end

"""Determine the type stored at key."""
function rtype(key::String)
    docommand(string("TYPE ", key))
end

"""Set the string value of a hash field."""
function hset(key::String, attr_name::String, attr_value::Any; pipeline::Bool=false)
    #TODO do_command(["HSET %s %s %s", key, attr_name, string(attr_value)])
    docommand(string("HSET ", key, " ", attr_name, " ", attr_value), pipeline)
end

"""Get the value of a hash field."""
function hget(key::String, attr_name::String; pipeline::Bool=false)
    docommand(string("HGET ", key, " ", attr_name), pipeline)
end

"""
Set multiple hash fields to multiple values. A variable number of arguments
follow key in `field` `value` format, e.g.:
    `hmset("myhash", "field1", "Hello", "field2", "World")`
"""
function hmset(key::String, argv::Any...; pipeline::Bool=false)
    cmd = string("HMSET ", key)
    for arg in argv
        cmd = string(cmd, " ", string(arg))
    end
    docommand(cmd, pipeline)
end

"""
Set multiple hash fields to multiple values. Fields and values are provided
as an Array.
"""
function hmset(key::String, attrs::AbstractVector{<: Any}; pipeline::Bool=false)
    cmd = string("HMSET ", key)
    for attr in attrs
        cmd = string(cmd, " ", string(attr))
    end
    docommand(cmd, pipeline)
end

"""
Set multiple hash fields to multiple values. Fields and values are provided
as a Dict.
"""
function hmset(key::String, attrs::Dict{String,Any}; pipeline::Bool=false)
    cmd = string("HMSET ", key)
    for (field, val) in attrs
        cmd = string(cmd, " ", field, " ", string(val))
    end
    docommand(cmd, pipeline)
end

"""
Get the values of all the given hash fields. Multiple fields are provided
as additional arguments.
"""
function hmget(key::String, argv::String...; pipeline::Bool=false)
    cmd = string("HMGET ", key)
    for arg in argv
        cmd = string(cmd, " ", arg)
    end
    docommand(cmd, pipeline)
end

"""
Get the values of all the given hash fields. Multiple fields are provided
as an Array.
"""
function hmget(key::String, fields::Array{String,1}; pipeline::Bool=false)
    cmd = string("HMGET ", key)
    for field in fields
        cmd = string(cmd, " ", field)
    end
    docommand(cmd, pipeline)
end

"""Get all the fields and values in a hash."""
function hgetall(key::String)
    reply::Array{Any,1} = do_command(string("HGETALL ", key))
    n = length(reply)
    dict = Dict{String,Any}()
    if n > 1 && mod(n, 2) == 0
        for i = 1:2:n
            dict[reply[i]] = reply[i + 1]
        end
    end
    dict
end

"""
Delete one or more hash fields. Multiple fields are provided
as additional arguments.
"""
function hdel(key::String, argv::String...; pipeline::Bool=false)
    cmd = string("HDEL ", key)
    for arg in argv
        cmd = string(cmd, " ", arg)
    end
    docommand(cmd, pipeline)
end

"""
Delete one or more hash fields. Multiple fields are provided
as an Array.
"""
function hdel(key::String, fields::Array{String,1}; pipeline::Bool=false)
    cmd = string("HDEL ", key)
    for field in fields
        cmd = string(cmd, " ", field)
    end
    docommand(cmd, pipeline)
end

"""Determine if a hash field exists."""
function hexists(key::String, field::String; pipeline::Bool=false)
    docommand(string("HEXISTS ", key, " ", field), pipeline)
end

"""Get all the fields in a hash."""
function hkeys(key::String; pipeline::Bool=false)
    docommand(string("HKEYS ", key), pipeline)
end

"""Get all the values in a hash."""
function hvals(key::String; pipeline::Bool=false)
    docommand(string("HVALS ", key), pipeline)
end

"""Get the number of fields in a hash."""
function hlen(key::String; pipeline::Bool=false)
    docommand(string("HLEN ", key), pipeline)
end

"""Increment the integer value of a hash field by the given number."""
function hincrby(key::String, field::String, increment::Int; pipeline::Bool=false)
    docommand(string("HINCRBY ", key, " ", field, " ", increment), pipeline)
end

"""Increment the float value of a hash field by the given number."""
function hincrby(key::String, field::String, increment::Float64; pipeline::Bool=false)
    docommand(string("HINCRBYFLOAT ", key, " ", field, " ", increment), pipeline)
end

"""Add one or more members to a set."""
function sadd(key::String, argv::Any...; pipeline::Bool=false)
    cmd = string("SADD ", key)
    for arg in argv
        # if isa(arg, Array) || isa(arg, Range)
        if isa(arg, Array)
            for member in arg
                cmd = string(cmd, " ", string(member))
            end
        else
            cmd = string(cmd, " ", string(arg))
        end
    end
    docommand(cmd, pipeline)
end

"""Get all the members in a set."""
function smembers(key::String; pipeline::Bool=false)
    docommand(string("SMEMBERS ", key), pipeline)
end

"""Determine if a given value is a member of a set."""
function sismember(key::String, member::Any; pipeline::Bool=false)
    docommand(string("SISMEMBER ", key, " ", string(member)), pipeline)
end

"""Get the number of members in a set."""
function scard(key::String; pipeline::Bool=false)
    docommand(string("SCARD ", key), pipeline)
end

"""
Remove one or more members from a set. Members are provided
as additional arguments.
"""
function srem(key::String, argv::Any...; pipeline::Bool=false)
    cmd = string("SREM ", key)
    for arg in argv
        cmd = string(cmd, " ", string(arg))
    end
    docommand(cmd, pipeline)
end

"""
Remove one or more members from a set. Members are provided
as an Array.
"""
function srem(key::String, members::AbstractVector{<: Any}; pipeline::Bool=false)
    cmd = string("SREM ", key)
    for member in members
        cmd = string(cmd, " ", string(member))
    end
    docommand(cmd, pipeline)
end

"""Subtract multiple sets. Multiple sets are provided as additional arguments."""
function sdiff(key::String, argv::String...; pipeline::Bool=false)
    cmd = string("SDIFF ", key)
    for arg in argv
        cmd = string(cmd, " ", arg)
    end
    docommand(cmd, pipeline)
end

"""Subtract multiple sets. Multiple sets are provided as an Array."""
function sdiff(keys::Array{String,1}; pipeline::Bool=false)
    cmd = "SDIFF"
    for key in keys
        cmd = string(cmd, " ", key)
    end
    docommand(cmd, pipeline)
end

"""Intersect multiple sets. Multiple sets are provided as additional arguments."""
function sinter(key::String, argv::String...; pipeline::Bool=false)
    cmd = string("SINTER ", key)
    for arg in argv
        cmd = string(cmd, " ", arg)
    end
    docommand(cmd, pipeline)
end

"""Intersect multiple sets. Multiple sets are provided as an Array."""
function sinter(keys::Array{String,1}; pipeline::Bool=false)
    cmd = "SINTER"
    for key in keys
        cmd = string(cmd, " ", key)
    end
    docommand(cmd, pipeline)
end

"""Add multiple sets. Multiple sets are provided as additional arguments."""
function sunion(key::String, argv::String...; pipeline::Bool=false)
    cmd = string("SUNION ", key)
    for arg in argv
        cmd = string(cmd, " ", arg)
    end
    docommand(cmd, pipeline)
end

"""Add multiple sets. Multiple sets are provided as an Array."""
function sunion(keys::Array{String,1}; pipeline::Bool=false)
    cmd = "SUNION"
    for key in keys
        cmd = string(cmd, " ", key)
    end
    docommand(cmd, pipeline)
end

"""Change the selected database for the current connection."""
function selectdb(db::Int)
    do_command(string("SELECT ", db))
end

"""Remove all keys from the current database."""
function flushdb()
    do_command("FLUSHDB")
end

"""Remove all keys from all databases."""
function flushall()
    do_command("FLUSHALL")
end

export start_session, end_session,                                                  # session
    kvset, kvget, incr, incrby, del, exists, getkeys, rdump, rtype,                 # key-value
    hset, hget, hmset, hmget, hgetall, hdel, hexists, hkeys, hvals, hlen, hincrby,  # hash sets
    sadd, smembers, sismember, scard, srem, sdiff, sinter, sunion,                  # sets
    selectdb, flushdb, flushall,                                                    # management
    @pipeline, get_reply,                                                           # pipelining
    do_command, pipeline_command                                                    # generic

end # module


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

"""
    FirstInfeasibleBranching()

Branching rule that always selects the first fractional variable.
"""
struct FirstInfeasibleBranching <: VariableBranchingRule end

function find_branching_var(
    rule::FirstInfeasibleBranching,
    node::Node,
    pool::NodePool,
)::Variable
    return node.fractional_variables[1]
end

"""
    LeastInfeasibleBranching()

Branching strategy that select the fractional variable whose value is the closest
to an integral value.
"""
struct LeastInfeasibleBranching <: VariableBranchingRule end

function find_branching_var(
    rule::LeastInfeasibleBranching,
    node::Node,
    pool::NodePool,
)::Variable
    scores = [max(v - floor(v), ceil(v) - v) for v in node.fractional_values]
    _, max_offset = findmax(scores)
    return node.fractional_variables[max_offset]
end

"""
    MostInfeasibleBranching()

Branching strategy that selects the fractional variable whose value is closest
to 1/2.
"""
struct MostInfeasibleBranching <: VariableBranchingRule end

function find_branching_var(
    rule::MostInfeasibleBranching,
    node::Node,
    pool::NodePool,
)::Variable
    scores = [min(v - floor(v), ceil(v) - v) for v in node.fractional_values]
    _, max_offset = findmax(scores)
    return node.fractional_variables[max_offset]
end


text:
# These reproduce the behaviour of DS9 according to http://ds9.si.edu/doc/ref/how.html

"""
    logstretch(num,a=1000)

A log-stretch as defined by the SAO DS9 application: http://ds9.si.edu/doc/ref/how.html
"""
logstretch(x,a=1000) = log(a*x+1)/log(a)
"""
    powstretch(num, a=1000)

A power-stretch as defined by the SAO DS9 application: http://ds9.si.edu/doc/ref/how.html
"""
powstretch(x,a=1000) = (a^x - 1)/a
"""
    sqrtstretch(num)

A square root stretch (simply defined as Base.sqrt)
"""
sqrtstretch(x) = sqrt(x)
"""
    squarestretch(num)

A squarestretch-stretch as defined by the SAO DS9 application: http://ds9.si.edu/doc/ref/how.html
"""
squarestretch(x) = x^2
"""
    asinhstretch(num)

A hyperbolic arcsin stretch as defined by the SAO DS9 application: http://ds9.si.edu/doc/ref/how.html.
"""
asinhstretch(x) = asinh(10x)/3
"""
    sinhstretch(num)

A hyperbolic sin stretch as defined by the SAO DS9 application: http://ds9.si.edu/doc/ref/how.html
"""
sinhstretch(x) = sinh(3x)/10
# These additional stretches reproduce behaviour from astropy
"""
    logstretch(num,a=1000)

A power distance stretch as defined by astropy.
"""
powerdiststretch(x, a=1000) = (a^x - 1) / (a - 1)

"""
    Percent(99.5)

Returns a callable that calculates display limits that include the given
percent of the image data.
Reproduces the behaviour of the SAO DS9 scale menu.

Example:
```julia
julia> imview(img, clims=Percent(90))
```
This will set the limits to be the 5th percentile to the 95th percentile.
"""
struct Percent
    perc::Float64
    trim::Float64
    Percent(percentage::Number) = new(Float64(percentage), (1 - percentage/100)/2)
end
(p::Percent)(data::AbstractArray) = quantile(vec(data), (p.trim, 1-p.trim))
(p::Percent)(data) = p(collect(data))
Base.broadcastable(p::Percent) = Ref(p)
Base.show(io::IO, p::Percent; kwargs...) = print(io, "Percent($(p.perc))", kwargs...)

"""
    Zscale(options)(data)

Wraps PlotUtils.zscale in a callable with default parameters.
This is a common algorithm for agressively stretching astronomical data
to see faint structure that originated in IRAF: https://iraf.net/forum/viewtopic.php?showtopic=134139
but is now seen in many other applications/libraries (DS9, Astropy, etc.)

Usage:
```
imview(img, clims=Zscale())
implot(img, clims=Zscale(contrast=0.1))
```

Default parameters:
```
nsamples::Int=1000
contrast::Float64=0.25
max_reject::Float64=0.5
min_npixels::Float64=5
k_rej::Float64=2.5
max_iterations::Int=5
```
"""
Base.@kwdef struct Zscale
    nsamples::Int=1000
    contrast::Float64=0.25
    max_reject::Float64=0.5
    min_npixels::Float64=5
    k_rej::Float64=2.5
    max_iterations::Int=5
end
(z::Zscale)(data::AbstractArray) = PlotUtils.zscale(vec(data), z.nsamples; z.contrast, z.max_reject, z.min_npixels, z.k_rej, z.max_iterations)
(z::Zscale)(data) = z(collect(data))
Base.show(io::IO, z::Zscale; kwargs...) = print(io, "Zscale()", kwargs...)
Base.broadcastable(z::Zscale) = Ref(z)

const _default_cmap  = Base.RefValue{Union{Symbol,Nothing}}(:magma)
const _default_clims = Base.RefValue{Any}(Percent(99.5))
const _default_stretch  = Base.RefValue{Any}(identity)

"""
    set_cmap!(cmap::Symbol)
    set_cmap!(cmap::Nothing)

Alter the default color map used to display images when using
`imview` or displaying an AstroImageMat.
"""
function set_cmap!(cmap)
    # Ensure it's valid
    _lookup_cmap(cmap)
    _default_cmap[] = cmap
end
"""
    set_clims!(clims::Tuple)
    set_clims!(clims::Callable)

Alter the default limits used to display images when using
`imview` or displaying an AstroImageMat.
"""
function set_clims!(clims)
    _default_clims[] = clims
end
"""
    set_stretch!(stretch::Function)

Alter the default value stretch functio used to display images when using
`imview` or displaying an AstroImageMat.
"""
function set_stretch!(stretch)
    _default_stretch[] = stretch
end

# Helper to iterate over data skipping missing and non-finite values.
skipmissingnan(itr) = Iterators.filter(el->!ismissing(el) && isfinite(el), itr)

# Convert argument into a colorscheme or AbstractColorList which allow converting
# from numerical data into colors.
function _lookup_cmap(cmap::Symbol)
    if cmap ∉ keys(ColorSchemes.colorschemes)
        error("$cmap not found in ColorSchemes.colorschemes. See: https://juliagraphics.github.io/ColorSchemes.jl/stable/catalogue/")
    end
    return ColorSchemes.colorschemes[cmap]
end
_lookup_cmap(::Nothing) = nothing
_lookup_cmap(acl::AbstractColorList) = acl
_lookup_cmap(colorant::Colorant) = PlotUtils.cgrad([:black, colorant])
_lookup_cmap(colorant::String) = PlotUtils.cgrad([:black, colorant])

function _resolve_clims(img::AbstractArray, clims)
    # Tuple or abstract array
    if typeof(clims) <: AbstractArray || typeof(clims) <: Tuple
        if length(clims) != 2
            error("clims must have exactly two values if provided.")
        end
        imgmin = first(clims)
        imgmax = last(clims)
    # Or as a callable that computes them given an iterator
    else
        imgmin, imgmax = clims(skipmissingnan(img))
    end

    return imgmin, imgmax
end

"""
    imview(img; clims=Percent(99.5), stretch=identity, cmap=:magma, contrast=1.0, bias=0.5)

Create a read only view of an array or AstroImageMat mapping its data values
to Colors according to `clims`, `stretch`, and `cmap`.

The data is first clamped to `clims`, which can either be a tuple of (min, max)
values or a function accepting an iterator of pixel values that returns (min, max).
By default, `clims=extrema` i.e. the minimum and maximum of `img`.
Convenient functions to use for `clims` are:
`extrema`, `zscale`, and `percent(p)`

Next, the data is rescaled to [0,1] and remapped according to the function `stretch`.
Stretch can be any monotonic fuction mapping values in the range [0,1] to some range [a,b].
Note that `log(0)` is not defined so is not directly supported.
For a list of convenient stretch functions, see:
`logstretch`, `powstretch`, `squarestretch`, `asinhstretch`, `sinhstretch`, `powerdiststretch`

Finally the data is mapped to RGB values according to `cmap`. If cmap is `nothing`,
grayscale is used. ColorSchemes.jl defines hundreds of colormaps. A few nice ones for
images include: `:viridis`, `:magma`, `:plasma`, `:thermal`, and `:turbo`.

Crucially, this function returns a view over the underlying data. If `img` is updated
then those changes will be reflected by this view with the exception of `clims` which
is not recalculated.

Note: if clims or stretch is a function, the pixel values passed in are first filtered
to remove non-finite or missing values.

### Defaults
The default values of `clims`, `stretch`, and `cmap` are `extrema`, `identity`, and `nothing`
respectively.
You may alter these defaults using `AstroImages.set_clims!`,  `AstroImages.set_stretch!`, and
`AstroImages.set_cmap!`.

### Automatic Display
Arrays wrapped by `AstroImageMat()` get displayed as images automatically by calling
`imview` on them with the default settings when using displays that support showing PNG images.

### Missing data
Pixels that are `NaN` or `missing` will be displayed as transparent when `cmap` is set
or black if.
+/- Inf will be displayed as black or white respectively.

### Exporting Images
The view returned by `imview` can be saved using general `FileIO.save` methods.
Example:
```julia
v = imview(data, cmap=:magma, stretch=asinhstretch, clims=percent(95))
save("output.png", v)
```
"""
function imview(
    img::AbstractArray{T};
    clims=_default_clims[],
    stretch=_default_stretch[],
    cmap=_default_cmap[],
    contrast=1.0,
    bias=0.5
) where {T}

    # Create flipped view of to match conventions of other programs.
    # Origin is centre of pixel (1,1) at bottom left.
    if ndims(img) == 2
        imgT = view(
            permutedims(img,(2,1)),
            reverse(axes(img,2)),
            :,
        )
    elseif ndims(img) >= 3
        newdims = (2,1, 3:ndims(img)...)
        ds = Tuple(((:) for _ in 2:ndims(img)))
        imgT = view(
            permutedims(img,newdims),
            reverse(axes(img,2)),
            ds...,
        )
    else
        imgT = img
    end
    isempt = isempty(imgT)
    if isempt
        @warn "imview called with empty argument"
        return fill(RGBA{N0f8}(0,0,0,0), 1,1)
    end
    # Users will occaisionally pass in data that is 0D, filled with NaN, or filled with missing.
    # We still need to do something reasonable in those caes.
    nonempty = any(x-> !ismissing(x) && isfinite(x), imgT)
    if !nonempty
        @warn "imview called with all missing or non-finite values"
        return map(px->RGBA{N0f8}(0,0,0,0), imgT)
    end

    imgmin, imgmax = _resolve_clims(imgT, clims)
    normed = clampednormedview(imgT, (imgmin, imgmax))
    return _imview(imgT, normed, stretch, _lookup_cmap(cmap), contrast, bias)
end

# Unwrap AstroImages before view, then rebuild.
# We have to permute the dimensions of the image to get the origin at the bottom left.
# But we don't want this to affect the dimensions of the array.
# Also, this reduces the number of methods we need to compile for imview by standardizing types
# earlier on. The compiled code for showing an array is the same as an array wrapped by an
# AstroImage, except for one unwrapping step.
function imview(
    img::AstroImage;
    kwargs...
)
    return shareheader(img, imview(parent(img); kwargs...))
end

# Special handling for complex images
"""
    imview(img::AbstractArray{<:Complex}; ...)

When applied to an image with complex values, display the magnitude
of the pixels using `imview` and display the phase angle as a panel below
using a cyclical color map.
For more customatization, you can create a view like this yourself:
```julia
vcat(
    imview(abs.(img)),
    imview(angle.(img)),
)
```
"""
function imview(img::AbstractArray{T}; kwargs...) where {T<:Complex}
    mag_view = imview(abs.(img); kwargs...)
    angle_view = imview(angle.(img), clims=(-pi, pi), stretch=identity, cmap=:cyclic_mygbm_30_95_c78_n256_s25)
    vcat(mag_view,angle_view)
end

function _imview(img, normed::AbstractArray{T}, stretch, cmap, contrast, bias) where T

    function colormap(pixr, pixn)::RGBA{N0f8}
        if ismissing(pixr) || !isfinite(pixr) || ismissing(pixn) || !isfinite(pixn)
            # We check pixr in addition to pixn because we want to preserve if the pixels
            # are +-Inf
            stretched = pixr
        else
            stretched = (stretch(pixn) - bias)*contrast+0.5
        end

        # We treat NaN/missing values as transparent
        pix= if ismissing(stretched) || isnan(stretched)
            RGBA{N0f8}(0,0,0,0)
        # We treat Inf values as white / -Inf as black
        elseif isinf(stretched)
            if stretched > 0
                RGBA{N0f8}(1,1,1,1)
            else
                RGBA{N0f8}(0,0,0,1)
            end
        else
            if isnothing(cmap)
                # true/false used as numerical values to prevent unucessary promotion
                s = clamp(stretched, false, true)
                RGBA{N0f8}(s,s,s,1)
            else
                # Look up colormap
                RGBA{N0f8}(get(cmap, stretched, (false, true)))
            end
        end
        return pix
    end
    mapper = mappedarray(colormap, img, normed)

    return maybe_copyheader(img, mapper)
end

"""
    imview_colorbar(img; orientation=:vertical)
Create a colorbar for a given image matching how it is displayed by
`imview`. Returns an image.

`orientation` can be `:vertical` or `:horizontal`.
"""
function imview_colorbar(
    img::AbstractArray;
    orientation=:vertical,
    clims=_default_clims[],
    stretch=_default_stretch[],
    cmap=_default_cmap[],
    contrast=1,
    bias=0.5
)
    imgmin, imgmax = _resolve_clims(img, clims)
    cbpixlen = 100
    data = repeat(range(imgmin, imgmax, length=cbpixlen), 1,10)
    if orientation == :vertical
        data = data'
    elseif orientation == :horizontal
        data = data
    else
        error("Unsupported orientation for colorbar \"$orientation\"")
    end

    # # Stretch the colors:
    # # Construct the image to use as a colorbar
    # cbimg =  imview(data; clims=(imgmin,imgmax), stretch, cmap, k_min=3)
    # # And the colorbar tick locations & labels
    # ticks, cbmin, cbmax = optimize_ticks(imgmin, imgmax)
    # # Now map these to pixel locations through streching and colorlimits:
    # stretchmin = stretch(zero(eltype(data)))
    # stretchmax = stretch(one(eltype(data)))
    # normedticks = clampednormedview(ticks, (imgmin, imgmax))
    # ticksloc = map(ticks,normedticks) do tick, tickn
    #     return cbpixlen * tickn
    # end

    # Strech the ticks
    # Construct the image to use as a colorbar
    cbimg = imview(data; clims=(imgmin,imgmax), stretch=identity, cmap, contrast, bias)
    # And the colorbar tick locations & labels
    ticks, _, _ = optimize_ticks(Float64(imgmin), Float64(imgmax), k_min=3)
    # Now map these to pixel locations through streching and colorlimits:
    stretchmin = stretch(zero(eltype(data)))
    stretchmax = stretch(oneunit(eltype(data)))
    normedticks = clampednormedview(ticks, (imgmin, imgmax))
    ticksloc = map(normedticks) do tickn
        stretched = stretch(tickn)
        stretchednormed = (stretched - stretchmin) * (stretchmax - stretchmin)
        return cbpixlen * stretchednormed
    end
    ticklabels = map(ticks) do t
        @sprintf("%4g", t)
    end
    return cbimg, (ticksloc, ticklabels)
end


text:
export ℝ, ℕ, 𝕋, Each, round_ceil, round_floor
export SysDynamics, SysDynamicsLinear, SysDynamicsLTI, discretize
export ObsDynamics, StaticObsDynamics
export sys_forward, limit_control, sys_A, sys_B, sys_w, obs_forward
export DelayModel, short_delay_names, total_delay, msg_queue_length
export WorldDynamics, RegretLossModel
export CentralControl, CentralControlStateless, init_state, control_one, control_all
export Controller, control!, write_logs
export ControllerFramework, make_controllers, MsgQueue
export Timed, attime, TimedQueue, pushpop!

const ℝ = Float64  # use 64-bit precision for real numbers
const ℕ = Int64
"The set of discrete times\n"
const 𝕋 = ℕ

round_ceil(x::ℝ)::ℕ = ℕ(ceil(x))
round_floor(x::ℝ)::ℕ = ℕ(floor(x))

"""
`Each{T}` is a vector containing elements of type `T`,
one for each agent in the fleet.
"""
Each{T} = Vector{T}

"""
A state dynamics model.

 - Should implement: `sys_forward`
 - Can optionally implement: `limit_control`
"""
abstract type SysDynamics end

"Simulate one time step forward using the system dynamics.\n"
function sys_forward(dy::SysDynamics, x::X, u, t::𝕋)::X where {X} @require_impl end

"""
    Restrict controls output by π into some valid range. The default implementation
    simply returns the original control unmodified.
"""
function limit_control(dy::SysDynamics, u::U, x, t)::U where {U}
    u
end

"""
A time-variant linear dynamics of the form ``x(t+1) = A(t) x(t) + B(t) u(t) + w(t)``.
Should implement 3 methods: `sys_A`, `sys_B`, and `sys_w`.
"""
abstract type SysDynamicsLinear <: SysDynamics end

function sys_A(sys::SysDynamicsLinear, t) @require_impl end
function sys_B(sys::SysDynamicsLinear, t) @require_impl end
function sys_w(sys::SysDynamicsLinear, t) @require_impl end

function sys_forward(dy::SysDynamicsLinear, x::X, u, t::𝕋)::X where {X}
    sys_A(dy, t) * x + sys_B(dy, t) * u + (sys_w(dy, t)::X)
end

"""
Get the corresponding system matrices from a linear system.
"""
sys_A, sys_B, sys_w

"""
An LTI Dynamics of the form ``x(t+1) = A x(t) + B u(t) + w(t)``.
"""
struct SysDynamicsLTI{MA,MB,NF <: Function} <: SysDynamicsLinear
    A::MA
    B::MB
    w::NF
end

sys_A(s::SysDynamicsLTI, t) = s.A
sys_B(s::SysDynamicsLTI, t) = s.B
sys_w(s::SysDynamicsLTI, t) = s.w(t)

function sys_forward(dy::SysDynamicsLTI, x::X, u, t::𝕋)::X where {X}
    dy.A * x + dy.B * u + convert(X, dy.w(t))
end

"""
Converts a continous LTI system, given in the form of (A, B), into a
discrete-time LTI system, given in the form (A′, B′).
"""
function discretize(A::TA, B::TB, delta_t::ℝ)::Tuple{TA,TB} where {TA,TB}
    A_int, err = quadgk(t -> exp(A .* t), 0.0, delta_t)
    # println("discretization error: $err")
    A′ = exp(A .* delta_t)
    B′ = A_int * B
    A′, B′
end

abstract type ObsDynamics end

"""
Assume the observation remains the same.
"""
struct StaticObsDynamics <: ObsDynamics end

"""
    obs_forward(dy::ObsDynamics, x, z::Z, t::𝕋)::Z where Z

Forward predict new observation from old observation and current state. Note that
the OneVision will only work well if `z` does not (or weakly) depend on `x`.
"""
function obs_forward(dy::ObsDynamics, x, z::Z, t::𝕋)::Z where {Z} @require_impl end

function obs_forward(
    dy::StaticObsDynamics, x, z, t::𝕋
) z end

struct DelayModel
    "observation delay: Tx"
    obs::𝕋
    "actuation delay: Tu"
    act::𝕋
    "communication delay: Tc"
    com::𝕋
    "control interaval, i.e., the time between adjacent control steps"
    ΔT::𝕋
end

total_delay(dm::DelayModel) = dm.obs + dm.act + dm.com

@inline short_delay_names(dm::DelayModel) =
    (Tx = dm.obs, Tu = dm.act, Tc = dm.com, Ta = total_delay(dm), ΔT = dm.ΔT)

DelayModel(;obs, act, com, ΔT = 1) = begin
    @assert com ≥ 1 "Communication delay should be at least 1, but got: $com."
    @assert obs ≥ 0 && act ≥ 0 && ΔT ≥ 1
    DelayModel(obs, act, com, ΔT)
end

msg_queue_length(dm::DelayModel) = dm.com

"""
A centralized controller with actuation type `U` and an internal state of type `S`.
`S` could be mutable and the framework will make a deep copy of the state when needed.

Should implement `control_one` and `init_state`.
"""
abstract type CentralControl{U,S} end

function control_one(
    π::CentralControl{U,S}, s::S, xs, zs, t::𝕋, id::ℕ
)::U  where {U,S}
    @require_impl
end

function control_all(
    π::CentralControl{U,S}, s::S, xs,zs, t::𝕋, ids::AbstractVector{ℕ}
)::AbstractVector{U} where {U,S}
    map(i -> control_one(π, s, xs, zs, t, i), ids)
end

function init_state(f::CentralControl{U,S}, t0)::S where {U,S}
    @require_impl
end

"""
A stateless central control. Only needs to implement the following function:
```
    control_one(f, xs, zs, t, id)
```
"""
abstract type CentralControlStateless{U} <: CentralControl{U,Nothing} end

function control_one(f::CentralControlStateless, xs, zs, t::𝕋, id::ℕ) @require_impl end

function control_one(
    f::CentralControlStateless{U}, ::Nothing, xs, zs, t::𝕋, id::ℕ
)::U  where {U}
    control_one(f, xs, zs, t, id)
end

function init_state(::CentralControlStateless, t0) nothing end

function init_state(::CentralControl{U,SymbolMaps.SymbolMap}, t0) where U
    SymbolMaps.SymbolMap()
end

"A distributed controller, typically generated from a controller framework.\n"
abstract type Controller{X,Z,U,Msg,Log} end

function control!(
    ctrl::Controller{X,Z,U,Msg},
    x::X,
    obs::Z,
    msgs::Each{Msg},
)::Tuple{U,Each{Msg}} where {X,Z,U,Msg}
    @require_impl
end

"""
This function will be called at the end of the simulation.
"""
function write_logs(
    ctrl::Controller{X,Z,U,Msg,Log} where {X,Z,U,Msg}
)::Dict{𝕋,Log} where Log
    @require_impl
end

struct WorldDynamics{N,XDynamics <: Tuple,ZDynamics <: Tuple}
    dynamics::XDynamics
    obs_dynamics::ZDynamics

    WorldDynamics(info::Each{<:Tuple{SysDynamics,ObsDynamics}}) = begin
        ds = Tuple(first.(info))
        os = Tuple(last.(info))
        new{length(info),typeof(ds),typeof(os)}(ds, os)
    end
end

struct RegretLossModel{N,X,U,S,XDynamics,UDynamics}
    central::CentralControl{U,S}
    world_model::WorldDynamics{N, XDynamics, UDynamics}
    x_weights::SVector{N, X}
    u_weights::SVector{N, U}
end

MsgQueue{Msg} = FixedQueue{Each{Msg}}

abstract type ControllerFramework{X,Z,U,Msg,Log} end

"""
Returns a tuple of `(controllers, message queues)`
"""
function make_controllers(
    framework::ControllerFramework{X,Z,U,Msg},
    init_status::Each{Tuple{X,Z,U}},
    init_t::𝕋,
) where {X,Z,U,Msg}
    @require_impl
end

struct Timed{X}
    time::𝕋
    value::X
end

Base.getindex(timed::Timed, time::𝕋) = begin
    @assert timed.time == time "expect time = $time, got: $timed."
    timed.value
end

Base.convert(::Type{Timed{Y}}, x::Timed{X}) where {X, Y} =
    if X == Y
        x
    else
        Timed{Y}(x.time, convert(Y, x.value))
    end

@inline attime(t::𝕋) = v -> Timed(t, v)

mutable struct TimedQueue{X}
    in_time::Int
    queue::FixedQueue{Timed{X}}
end

"""
Creates a TimedQueue with all initial elements having the same value.
"""
function TimedQueue(value::X, t1::𝕋, t2::𝕋; eltype::Type{E} = X) where {X, E}
    in_time = t2+1
    queue = FixedQueue([Timed{E}(t, value) for t in t1:t2])
    TimedQueue(in_time, queue)
end

Base.length(q::TimedQueue) = length(q.queue)
Base.isempty(q::TimedQueue) = isempty(q.queue)
Base.getindex(q::TimedQueue, t::𝕋) = q.queue[t]

function pushpop!(q::TimedQueue{T}, x::Timed) where T
    @assert x.time == q.in_time "Queue in_time: $(q.in_time), x = $x"
    q.in_time += 1
    pushpop!(q.queue, x)
end


text:
struct CapacityExpansionScenario{T₁,T₂} <: AbstractScenario
    probability::Probability
    ρ::PriceCurve{T₁,Float64}
    Q̃::InflowSequence{T₂,typeof(Skellefteälven),Float64}

    function CapacityExpansionScenario(ρ::PriceCurve{T₁,Float64}, Q̃::InflowSequence{T₂,typeof(Skellefteälven),Float64}; probability::AbstractFloat = 1.0) where {T₁,T₂}
        T₁ == 24*T₂ || error("Incompatible price period $T₁ and flow period $T₂.")
        return new{T₁,T₂}(Probability(probability), ρ, Q̃)
    end
end

function Base.zero(::Type{CapacityExpansionScenario{T₁,T₂}}) where {T₁,T₂}
    return CapacityExpansionScenario(zero(PriceCurve{T₁,Float64}), zero(InflowSequence{T₂,typeof(Skellefteälven),Float64}); probability = 1.0)
end

function StochasticPrograms.expected(ξ₁::CapacityExpansionScenario, ξ₂::CapacityExpansionScenario)
    ρ = mean([ξ₁.ρ, ξ₂.ρ])
    Q̃ = mean([ξ₁.Q̃, ξ₂.Q̃])
    expected = CapacityExpansionScenario(ρ, Q̃)
    return ExpectedScenario(expected)
end

struct RecurrentCapacityExpansionSampler{T₁,T₂} <: AbstractSampler{CapacityExpansionScenario{T₁,T₂}}
    date::Date
    plants::PlantCollection
    price_forecaster::Forecaster
    flow_forecaster::Forecaster
    electricity_price_rate::Float64
    horizon::Horizon

    function RecurrentCapacityExpansionSampler(date::Date,
                                               plants::PlantCollection,
                                               price_forecaster::Forecaster,
                                               flow_forecaster::Forecaster,
                                               electricity_price_rate::Float64,
                                               horizon::Horizon)
        T₁ = num_hours(horizon)
        T₂ = num_days(horizon)
        return new{T₁,T₂}(date, plants, price_forecaster, flow_forecaster, electricity_price_rate, horizon)
    end
end

function (sampler::RecurrentCapacityExpansionSampler)()
    price_rate = sampler.electricity_price_rate * rand()
    price_curve = forecast(sampler.price_forecaster, month(sampler.date))
    while !all(price_curve .>= 0)
        price_curve = forecast(sampler.price_forecaster, month(sampler.date))
    end
    for d = 1:num_days(sampler.horizon)-1
        date = sampler.date + Dates.Day(d)
        years_passed = Dates.value(Dates.Year(date) - Dates.Year(sampler.date))
        daily_curve = (1 + price_rate)^years_passed * forecast(sampler.price_forecaster, price_curve[end], month(date))
        while !all(daily_curve .>= 0)
            daily_curve = (1 + price_rate)^years_passed * forecast(sampler.price_forecaster, price_curve[end], month(date))
        end
        append!(price_curve, daily_curve[:])
    end
    flows = forecast(sampler.flow_forecaster, week(sampler.date))
    while !all(flows[:, 2] .>= 0)
        flows = forecast(sampler.flow_forecaster, week(sampler.date))
    end
    for w = 1:num_weeks(sampler.horizon)-1
        date = sampler.date + Dates.Week(w)
        weekly_flows = forecast(sampler.flow_forecaster, flows[:,end], week(date))
        while !all(weekly_flows[:, 2] .>= 0)
            weekly_flows = forecast(sampler.flow_forecaster, flows[:,end], week(date))
        end
        flows = hcat(flows, weekly_flows)
    end
    if size(flows, 2) < num_days(sampler.horizon)
        date = sampler.date + Dates.Week(num_weeks(sampler.horizon)-1)
        weekly_flows = forecast(sampler.flow_forecaster, flows[:,end], week(date))
        while !all(weekly_flows[:, 2] .>= 0)
            weekly_flows = forecast(sampler.flow_forecaster, flows[:,end], week(date))
        end
        flows = hcat(flows, weekly_flows[:,1:num_days(sampler.horizon)-size(flows,2)])
    end
    return CapacityExpansionScenario(PriceCurve(price_curve), InflowSequence(sampler.plants, flows))
end


text:
__precompile__()
module TaaSSim

#optimisation
using JuMP

max_stars_repo_name: jprendo/TaaSSimTest

max_issues_repo_name: jprendo/TaaSSimTest

max_forks_repo_name: jprendo/TaaSSimTest