text:
"""
Differential Evolution MCMC object
* `n_groups`: number of groups of particles. Default = 4
* `Np`: number of particles per chain. Default = number of parameters * 3 if
    priors are passed
* `burnin`: number of burnin iterations. Default = 1000
* `α`: migration probability. Default = .1
* `β`: mutation probability. Default = .1
* `ϵ`: noise in crossover step. Default = .001
* `σ`: standard deviation of noise added to parameters for mutation. Default = .05
* `κ`: recombination with probability (1-κ) during crossover. Default = 1.0
* `visualize`: generate scatter plot of particle locations on each iteration. Default=false
* `iter`: current iteration
* `progress`: show progress meter if true. Default = true

Constructor signature:

```@example
DE(;n_groups=4, priors=nothing, Np=num_parms(priors)*3,burnin=1000,
    α=.1, β=.1, ϵ=.001, σ=.05, bounds, visualize=false, iter=1, progress=false)
```
References:

* Ter Braak, C. J. A Markov Chain Monte Carlo version of the genetic algorithm Differential Evolution: easy Bayesian computing for real parameter spaces.

* Turner, B. M., Sederberg, P. B., Brown, S. D., & Steyvers, M. (2013). A method for efficiently sampling from distributions with correlated dimensions. Psychological methods, 18(3), 368.

* Turner, B. M., & Sederberg, P. B. (2012). Approximate Bayesian computation with differential evolution. Journal of Mathematical Psychology, 56(5), 375-385.
"""
mutable struct DE{T1} <: AbstractSampler
    n_groups::Int64
    Np::Int64
    burnin::Int64
    α::Float64
    β::Float64
    ϵ::Float64
    σ::Float64
    κ::Float64
    bounds::T1
    iter::Int64
    progress::Bool
end

function DE(;n_groups=4, priors=nothing, Np=num_parms(priors)*3, burnin=1000, α=.1, β=.1, ϵ=.001, σ=.05,
    κ=1.0, bounds, progress=true)
    return DE(n_groups, Np, burnin, α, β, ϵ, σ, κ, bounds, 1, progress)
end

"""
A model object containing the log likelihood function and prior distributions
* `priors`: prior distributions
* `model`: log likelihood function
* `names`: parameter names
"""
struct DEModel{F,L,T} <: AbstractModel where {F<:Function,L,T}
    priors::L
    model::F
    names::T
end

function DEModel(;priors, model, names=String.(keys(priors)))
    priors′ = values(priors)
    return DEModel(priors′, model, names)
 end

"""
Computes the number of parameters based on scalars or vectors
"""
function num_parms(priors)
    if isnothing(priors)
        error("Np undefined. Define Np=x in constructor or pass priors for default Np")
    end
    n = 0
    for p in priors
        length(p) == 1 ? n += 1 : (n += p[2])
    end
    return n
end

"""
* `θ`: a vector of parameters
* `samples`: a 2-dimensional array containing all acccepted proposals
* `accept`: proposal acceptance. 1: accept, 0: reject
* `weight`: particle weight based on model fit (currently posterior log likelihood)
* `lp`: a vector of log posterior probabilities associated with each accepted proposal
"""
mutable struct Particle{T}
    Θ::Vector{T}
    samples::Array{T,2}
    accept::Vector{Bool}
    weight::Float64
    lp::Vector{Float64}
end

Base.broadcastable(x::Particle) = Ref(x)

function Particle(;Θ=[.0], samples=Array{eltype(Θ),2}(undef,1,1),
    accept=Bool[],weight=0.0)
    Particle(Θ, samples, accept, weight, Float64[])
end


text:
# a DisjointSet is a structure that stores counts of the number of observations that
# fall into each element of the power set of the given sets. Note that a PowerSet
# would be similar, but (for example) the slot for A would contain A&B, whereas for
# a DisjointSet, the slot for A is for observations of (A and not-B and not-C), etc.
mutable struct DisjointSet
	labels
	counts

	function DisjointSet(dat::AbstractMatrix, labs::AbstractVector)
		# dat must be 2-D and the num of cols must be equal to the length of the labels
		@assert size(dat,2) == size(labs,1)
		nrow = size(dat,1); ncol = size(dat,2)

		pssize = 2^ncol
		cts = zeros(Int64, pssize)

		for row in 1:nrow
			# which cell to increment? start with 0 and add powers of 2 from right to left
			# then add 1 because Julia is 1-based
			i = 0
			for col in 1:ncol
				if (dat[row, col])
					i += 2^(ncol-col)
				end
			end
			cts[i+1] += 1
		end
		new(labs, cts)
	end
end

# an EulerState is just a Vector{Float64}
const EulerState = Vector{Float64}

mutable struct EulerObject
	nparams
	labels
	lb
	ub
	sizes
	target
	specs
	evalfn
end

mutable struct EulerSpec
	shape
	clamp
	statepos

	function EulerSpec(shape, clamp, statepos)
		# this constructor enforces invariants
		if in(shape, [:circle, :square, :triangle])
			@assert length(clamp) == 2
			@assert length(statepos) == 2
		elseif in(shape, [:rectangle])
			@assert length(clamp) == 3
			@assert length(statepos) == 3
		else
			error("Unknown EulerSpec shape: ", string(shape))
		end
		new(shape, clamp, statepos)
	end
end

EulerSpec() = EulerSpec(:circle)
function EulerSpec(shape)
	if in(shape, [:circle, :square, :triangle])
		EulerSpec(shape, [NaN, NaN], [0, 0])
	elseif in(shape, [:rectangle])
		EulerSpec(shape, [NaN, NaN, NaN], [0, 0, 0])
	else
		error("Unknown EulerSpec shape: ", string(shape))
	end
end

ext: jl


text:
using NeMo
using Base.Test

@test delconf(newconf()) == 0
@test synapseswo() == 0
## TODO: add all tests


text:
using SphericalFourierBesselDecompositions
using Documenter

DocMeta.setdocmeta!(SphericalFourierBesselDecompositions, :DocTestSetup, :(using SphericalFourierBesselDecompositions); recursive=true)

magickopts = `-density 300 -background none -define icon:auto-resize=256,128,96,64,48,32,16`
cmd = `convert $(@__DIR__)/src/assets/favicon.svg $magickopts $(@__DIR__)/src/assets/favicon.ico`
@show cmd
run(cmd)

makedocs(;
    modules=[SphericalFourierBesselDecompositions],
    authors="Henry Gebhardt <henry.s.gebhardt@jpl.nasa.gov> and contributors",
    repo="https://github.com/hsgg/SphericalFourierBesselDecompositions.jl/blob/{commit}{path}#{line}",
    sitename="SphericalFourierBesselDecompositions.jl",
    #sitename="[SphericalFourierBesselDecompositions.jl](https://github.com/hsgg/SphericalFourierBesselDecompositions.jl)",
    #sitename="SphericalFourierBesselDecompositions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://hsgg.github.io/SphericalFourierBesselDecompositions.jl",
	assets = ["assets/favicon.ico"],
    ),
    pages=[
        "Home" => "index.md",
        "Tutorial: Catalogue and Window" => "tutorial_catalog.md",
        "Reference" => "reference.md",
        "Index" => "myindex.md",
    ],
)

deploydocs(;
    repo="github.com/hsgg/SphericalFourierBesselDecompositions.jl",
)


text:
mutable struct VAR
    endogenous::Matrix{Float64}
    exogenous::Matrix{Float64}
    endocoefs::Array{Float64, 3}
    exocoefs::Matrix{Float64}
end

eigen(  obj::VAR) = eigen(  companion(obj))
eigvals(obj::VAR) = eigvals(companion(obj))

function companion(obj::VAR)
    p = size(obj.endogenous, 2)
    k = size(obj.endocoefs, 3)
    return [reshape(obj.endocoefs, p, p*k); Matrix{Float64}(I, (k-1)*p, k*p)]
end

function simulate(obj::VAR, innovations::Matrix{Float64}, init::Matrix{Float64})
    iT, p = size(innovations)
    k = size(obj.endocoefs, 3)
    if size(obj.endocoefs, 1) != p
         throw(ArgumentError("Wrong dimensions"))
    end
    if size(init, 2) != p
         throw(ArgumentError("Wrong dimensions"))
    end
    if size(init, 1) != k
        throw(ArgumentError("Wrong dimensions"))
    end
    Y = zeros(iT + k, p)
    Y[1:k, :] = init
    for t = 1:iT
        for i = 1:k
            Y[t + k, :] += (obj.endocoefs[:, :, i]*Y[t + k - i, :] + innovations[t, :])
        end
    end
    return Y
end
simulate(obj::VAR, innovations::Matrix{Float64}) = simulate(obj, innovations, obj.endogenous[1:size(obj.endocoefs, 3),:])
simulate(obj::VAR, iT::Integer) = simulate(obj, randn(iT, size(obj.endogenous, 2)))

function convert(::Type{VAR}, obj::CivecmI1)
    p = size(obj.endogenous, 2)
    endocoefs = Array{Float64}(undef, p, p, obj.lags)
    endocoefs[:, :, 1] = obj.α*obj.β[1:p, :]' + I
    if obj.lags > 1
        endocoefs[:, :, 1] += obj.Γ[:,1:p]
        endocoefs[:, :, obj.lags] = -obj.Γ[:,(obj.lags-2)*p+1:(obj.lags-1)*p]
        for i = 1:obj.lags - 2
            endocoefs[:,:,i+1] = obj.Γ[:,p*(obj.lags-i-1)+1:p*(obj.lags-i)] - obj.Γ[:,p*(obj.lags-i-2)+1:p*(obj.lags-i-1)]
        end
    end
    return VAR(obj.endogenous, obj.exogenous, endocoefs, zeros(0,0))
end

function convert(::Type{VAR}, obj::CivecmI2)
    p = size(obj.endogenous, 2)
    endocoefs = Array{Float64}(p, p, obj.lags)
    endocoefs[:,:,1] = 2I + (obj.α*ρ(obj)'*τ(obj)' + obj.α*δ(obj)'*obj.τ⊥' + obj.ζt'*τ(obj)')[1:p,1:p]
    endocoefs[:,:,2] = -I - (obj.α*δ(obj)'*obj.τ⊥' + obj.ζt'*τ(obj)')[1:p,1:p]
    if obj.lags > 2
        Ψ = (obj.Z3\(obj.Z0 - obj.Z2*obj.τ*ρ(obj)*obj.α' - obj.Z1*(obj.τ⊥*δ(obj)*obj.α' + obj.τ*obj.ζt)))'
        endocoefs[:,:,1] += Ψ[:,1:p]
        endocoefs[:,:,obj.lags] = -Ψ[:,(obj.lags-3)*p+1:(obj.lags-2)*p]
        if obj.lags > 3
            endocoefs[:,:,2] += Ψ[:,p+1:2p] - 2*Ψ[:,1:p]
            endocoefs[:,:,obj.lags-1] = Ψ[:,(obj.lags-4)*p+1:(obj.lags-3)*p] - 2*Ψ[:,(obj.lags-3)*p+1:(obj.lags-2)*p]
        end
        for i = 1:obj.lags - 4
            endocoefs[:,:,i+2] = Ψ[:,p*(i-1)+1:p*i] - 2*Ψ[:,p*i+1:p*(i+1)] + Ψ[:,p*(i+1)+1:p*(i+2)]
        end
    end
    return VAR(obj.endogenous, obj.exogenous, endocoefs, zeros(0,0))
end


text:
using MAT

@testset "simulateDSM" begin
    # MATLAB output
    test = matread("resources/simulateDSM.mat")

    # ntf input test
    u1 = 0.5 * sin.(2*pi*85/8192 * (0:8191))
    u2 = 0.5 * sin.(2*pi*1045/8192 * (0:8191))
    u3 = 8 * sin.(2*pi*146/8192 * (0:8191))
    u4 = 16 * undbv(-3) * sin.(2*pi*164/16384 * (0:16383))
    u5 = 16 * undbv(-30) * sin.(2*pi*164/16384 * (0:16383))
    ntf1 = synthesizeNTF(5, 32, 1)
    ntf2 = synthesizeNTF(8, 64, 1, 1.5, 0.125)
    ntf3 = synthesizeNTF(7, 8, 1, 2)
    ntf4 = synthesizeNTF(7, 8, 1, 8)
    ntf5 = synthesizeNTF(6, 25, 1, 4)
    u    = [  u1,   u2,   u3,   u3,   u4,   u5]
    ntf  = [ntf1, ntf2, ntf3, ntf4, ntf5, ntf5]
    nlev = [   2,    2,   17,   17,   17,   17]
    for i = 1:length(u)
        v, _, _, y = simulateDSM(u[i], ntf[i], nlev[i])
        @test vec(test["v_$i"]) == v
        # Allow higher error tolerance because of differences in MATLAB and
        # Julia state space models. This also why xn and xmax are ignored;
        # the x vectors aren't the same. Default tolerance is sqrt(eps).
        @test isapprox(vec(test["y_$i"]), y, rtol=10*sqrt(eps(Float64)))
    end

    # ABCD input test
    u = range(0, 0.6; length=30)
    ntf = synthesizeNTF(5, 42, 1)
    a, g, b, c = realizeNTF(ntf)
    ABCD = stuffABCD(a, g, b, c)
    for i = 1:length(u)
        v, xn, xmax, y = simulateDSM(u[i]*ones(10000), ABCD)
        @test vec(test["v"][i:i,:]) == v
        @test isapprox(test["xn"][i,:,:], xn)
        @test isapprox(test["xmax"][i,:,:], xmax)
        @test isapprox(test["y"][i,:,:], y)
    end
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

# # @typegraph example
# In this notebook, we will be going through a tutorial on how to use the `@typegraph` macro to extract the mapping of how type are transformed through functions.

using SemanticModels.ModelTools

include("../../src/modeltools/typegraph.jl");

# After loading `@typegraph` in to our workspace, we simply begin the extraction by calling the macro immediately followed by an `Expr` which will return an edge list which is easily passed through to `MetaGraphs.jl` to visualize the transformations that are taking place throughout the code.
#
# *_To learn more about `Expr` & metaprogramming, we recommend looking at the offcial [Julia docs](@https://docs.julialang.org/en/v1.0/manual/metaprogramming/) on the topic._

edgelist = @typegraph begin
    import Base: count

    """    AgentModel

    the root type for agent based models.

    See also: StateModel
    """
    abstract type AgentModel end

    """     StateModel

    holds the components of an agent based simulation using finite state machines.

    - states: a collection of distinct states an agent can occupy
    - agents: a collection of states `aᵢ = agents[i] ∈ states` indicating that agent `i` is in state `aᵢ`
    - transitions: the functions `f: states -> states`
    """
    mutable struct StateModel <: AgentModel
        states
        agents
        transitions
        loads
    end

    # +
    # Determine number of agents in a given state
    function count(sm::StateModel, state)
        return length(filter(a->(a==state), sm.agents))
    end

    # counts the number of agents
    function count(sm::StateModel)
        return length(sm.agents)
    end
    # -

    # stateload computes the fraction of agents in each state
    # it is used by tick! to update the statemodel for computing the probability of infection.
    function stateload(sm::StateModel, state::Symbol)
        return (count(sm, state)+1)/(count(sm)+1)
    end

    #     tick!(sm::StateModel)
    #
    # performs the operations that need to happen once per time step. Things like caching values that need to be computed at the beginning of the timestep like the distribution of the agents' states.

    function tick!(sm::StateModel)
        sm.loads = map(s->stateload(sm, s), sm.states)
        return sm.loads
    end

    #     step!(sm::StateModel, n=1)
    #
    # advance the simulation by `n` ticks of time.
    # This is an in-place operation that modifies the current state of the simulation.
    #

    function step!(sm::StateModel, n)
        for s in 1:n
          tick!(sm)
          for (i, a) in enumerate(sm.agents)
              sm.agents[i] = sm.transitions[a](sm, i,a)
          end
            describe(sm)
        end
        return sm
    end

    """    describe(sm::StateModel)

    summarize the state of the simulation for presentation or analysis.
    """
    function describe(sm::StateModel)
        counts = zeros(Int, size(sm.states))
        d = Dict{eltype(sm.states), Int}()
        for (i,s) in enumerate(sm.states)
            d[s] = i
        end
        for a in sm.agents
            i = d[a]
            counts[i] += 1
        end
        return collect(map(x->Pair(x...), zip(sm.states, counts)))
    end;

    # ## Run the model
    #
    # This script defines a basic agent based model of disease spread called SIRS. Each agent is in one of 3 states
    #
    # 1. $S$ Susceptible
    # 2. $I$ Infected
    # 3. $R$ Recovered
    #
    # <img src="https://docs.google.com/drawings/d/e/2PACX-1vSeA7mAQ-795lLVxCWXzbkFQaFOHMpwtB121psFV_2cSUyXPyKMtvDjssia82JvQRXS08p6FAMr1hj1/pub?w=1031&h=309">
    #
    # The agents go from `S->I`, `I-R`, and `R->S` based on random numbers. The probability of S-> is dependent on the fraction of agents in state :I. The probability of recovering is a constant ρ, and the disease confers some temporary immunity with probability μ.
    #
    # This script has an entrypoint to call it so that you can include this file and run as many simulations as you want. The intended use case is to repeatedly call `main` and accumulate the return values into an array for later analysis.

    function main(nsteps)
        n = 20
        a = fill(:S, n)
        ρ = 0.5 + randn(Float64)/4 # chance of recovery
        μ = 0.5 # chance of immunity
        T = Dict(
            :S=>(x...)->rand(Float64) < stateload(x[1], :I) ? :I : :S,
            :I=>(x...)->rand(Float64) < ρ ? :I : :R,
            :R=>(x...)->rand(Float64) < μ ? :R : :S,
        )

        sam = StateModel([:S, :I, :R], a, T, zeros(Float64,3))
        newsam = step!(deepcopy(sam), nsteps)
        @show newsam.agents
        counts = describe(newsam)
        return newsam, counts
    end;

    # +
    # # An example of how to run this thing.
    main(10);
    # -

    end;

# In the above example, we have a simple agent based simulation where we defined new stucts to collect the singltonian type information for our simulation. To reduce the nosie from collecting the iterations throughout the runtime of our model we need to collect the unique calls through our our example which contain the relevant information.

E = unique((f.func, f.args, f.ret) for f in edgelist)

# ## visualizing the edges
#
# Now that we have extracted the relevant type information, we want to visualize these transformations in a knowledge graph.

using MetaGraphs;
using LightGraphs;

g = MetaDiGraph();
set_indexing_prop!(g,:label);

for e in E
    try
        g[e[2],:label]
    catch
        add_vertex!(g,:label,e[2]) # add ags
    end

    try
        g[e[3],:label]
    catch
        add_vertex!(g,:label,e[3]) # add rets
    end

    try
        add_edge!(g,g[e[2],:label],g[e[3],:label],:label,e[1]) # add func edges
    catch
        nothing
    end
end

savegraph("exampletypegraph.dot",g,DOTFormat());

run(`dot -Tsvg -O exampletypegraph.dot`);

s = read("exampletypegraph.dot.svg",String);
display("image/svg+xml",s)

# +
# TODO use a different script and repeat, extraction, drawing, and rendering

# +
# TODO discuss difference in type graphs


text:
using LaTeXTabulars, Test, LaTeXStrings

# for testing
using LaTeXTabulars: latex_cell

"Normalize whitespace, for more convenient testing."
squash_whitespace(string) = strip(replace(string, r"[ \n\t]+" => " "))

@test squash_whitespace(" something  \n with line  breaks \n  and   stuff \n") ==
    "something with line breaks and stuff"

"Comparison using normalized whitespace. For testing."
≅(a, b) = squash_whitespace(a) == squash_whitespace(b)

@testset "tabular" begin
    tb = Tabular("lcl")
    tlines = [Rule(:top),
              [L"\alpha", L"\beta", "sum"],
              Rule(:mid),
              [1, 2, 3],
              Rule(),           # a nice \hline to make it ugly
              [4.0 "5" "six";   # a matrix
               7 8 9],
              (CMidRule(1, 2), CMidRule("lr", 1, 1)), # just to test tuples
              [MultiColumn(2, :c, "centered")],       # ragged!
              Rule(:bottom)]
    tlatex = raw"\begin{tabular}{lcl}
                 \toprule
                 $\alpha$ & $\beta$ & sum \\
                 \midrule
                 1 & 2 & 3 \\
                 \hline
                 4.0 & 5 & six \\
                 7 & 8 & 9 \\ \cmidrule{1-2} \cmidrule(lr){1-1}
                 \multicolumn{2}{c}{centered} \\
                 \bottomrule
                 \end{tabular}"

    tlatex = replace(tlatex, "\r\n"=>"\n")
    @test latex_tabular(String, tb, tlines) ≅ tlatex
    tmp = tempname()
    latex_tabular(tmp, tb, tlines)
    @test isfile(tmp) && read(tmp, String) ≅ tlatex
    @test read(tmp, String) ≅ tlatex
end

@test_throws ArgumentError latex_cell(stdout, MultiColumn(2, :BAD, ""))
@test_throws ArgumentError CMidRule(3, 1)     # not ≤
@test_throws MethodError latex_cell(stdout, ("un", "supported"))
@test_throws MethodError CMidRule(1, 1, 1, 2) # invalid types

@testset "longtable" begin
    lt = LongTable("rrr", ["alpha", "beta", "gamma"])
    tlines = [[1 2 3 ;
               4.0 "5" "six"],
              Rule(:h)]
    tlatex = raw"\begin{longtable}[c]{rrr}
                 \hline
                 alpha & beta & gamma \\
                 \hline
                 \endfirsthead
                 \multicolumn{3}{l}
                 {{\bfseries \tablename\ \thetable{} --- continued from previous page}} \\
                 \hline
                 alpha & beta & gamma \\
                 \hline
                 \endhead
                 \hline
                 \multicolumn{3}{r}{{\bfseries Continued on next page}} \\
                 \hline
                 \endfoot
                 \hline
                 \endlastfoot
                 1 & 2 & 3 \\
                 4.0 & 5 & six \\
                 \hline
                 \end{longtable}"

    tlatex = replace(tlatex, "\r\n"=>"\n")
    @test latex_tabular(String, lt, tlines) ≅ tlatex
    tmp = tempname()
    latex_tabular(tmp, lt, tlines)
    @test isfile(tmp) && read(tmp, String) ≅ tlatex
    @test read(tmp, String) ≅ tlatex
end


text:
__precompile__(false)

module Buildkite

import HTTP
import JSON

export BuildkiteAPI

abstract type AbstractBuildkiteAPI end

Base.@kwdef struct BuildkiteAPI <: AbstractBuildkiteAPI
    base_url::HTTP.URI = HTTP.URI("https://api.buildkite.com/v2/")
    access_token::String = ""
end

function buildkite_api_uri(api::BuildkiteAPI, path)
    HTTP.URIs.merge(api.base_url, path = api.base_url.path * path)
end

function buildkite_request(api::BuildkiteAPI, request_method, endpoint;
                           handle_error = true,
                           headers = Dict(),
                           params = Dict(),
                           allowredirects::Bool = true,
                           idle_timeout = 20,
                           status_exception = false)
    api_endpoint = buildkite_api_uri(api, endpoint)
    request_headers = convert(Dict{String, String}, headers)
    if !haskey(request_headers, "Authorization")
        request_headers["Authorization"] = "Bearer $(api.access_token)"
    end
    if !haskey(request_headers, "User-Agent")
        request_headers["User-Agent"] = "Buildkite-jl"
    end
    if request_method == HTTP.get
        api_uri = HTTP.URIs.merge(api_endpoint, query = params)
        println("DEBUG: ", api_uri)
        r = request_method(api_uri, request_headers,
                           redirect = allowredirects, status_exception = false, idle_timeout=idle_timeout)
    else
        api_uri = string(api_uri)
        r = request_method(api_uri, request_headers, JSON.json(params),
                           redirect=allowredirects,
                           status_exception=status_exception,
                           idle_timeout=idle_timeout)
    end
    if handle_error
        #handle_response_error(r)
    end
    return r
end

# REST primitives

function buildkite_get(api::BuildkiteAPI, endpoint = ""; options...)
    buildkite_request(api, HTTP.get, endpoint; options...)
end

function buildkite_post(api::BuildkiteAPI, endpoint = ""; options...)
    buildkite_request(api, HTTP.post, endpoint; options...)
end

function buildkite_put(api::BuildkiteAPI, endpoint = ""; options...)
    buildkite_request(api, HTTP.put, endpoint; options...)
end

function buildkite_delete(api::BuildkiteAPI, endpoint = ""; options...)
    buildkite_request(api, HTTP.delete, endpoint; options...)
end

function buildkite_patch(api::BuildkiteAPI, endpoint = ""; options...)
    buildkite_request(api, HTTP.patch, endpoint; options...)
end

function buildkite_get_json(api::BuildkiteAPI, endpoint = ""; options...)
    JSON.parse(HTTP.payload(buildkite_get(api, endpoint; options...), String))
end

function buildkite_post_json(api::BuildkiteAPI, endpoint = ""; options...)
    JSON.parse(HTTP.payload(buildkite_post(api, endpoint; options...), String))
end

function buildkite_put_json(api::BuildkiteAPI, endpoint = ""; options...)
    JSON.parse(HTTP.payload(buildkite_put(api, endpoint; options...), String))
end

function buildkite_delete_json(api::BuildkiteAPI, endpoint = ""; options...)
    JSON.parse(HTTP.payload(buildkite_delete(api, endpoint; options...), String))
end

function buildkite_patch_json(api::BuildkiteAPI, endpoint = ""; options...)
    JSON.parse(HTTP.payload(buildkite_patch(api, endpoint; options...), String))
end

function hello_world()
    base_url = HTTP.URI("https://api.buildkite.com")
    r = buildkite_get_json(BuildkiteAPI(base_url=base_url), "")
    return r["response"]
end

# organization api

struct Organization
    api::BuildkiteAPI
    data::Dict
end

function organization(api::BuildkiteAPI, name)
    return Organization(api, buildkite_get_json(api, "organizations/$(lowercase(name))"))
end

# pipelines api

struct Pipeline
    api::BuildkiteAPI
    data::Dict
end

function pipelines(api::BuildkiteAPI, organization; page=0, pagination=false)
    query_params = Dict("page" => page)
    endpoint = "organizations/$(lowercase(organization))/pipelines"
    return [Pipeline(api, p) for p in buildkite_get_json(api, endpoint; params = query_params)]
end

function pipelines(org::Buildkite.Organization)
    return pipelines(org.api, org.data["name"])
end

# build api

struct Build
    api::BuildkiteAPI
    data::Dict
end

function builds(api::BuildkiteAPI; state=nothing)
    query_params = Dict("state" => state)
    endpoint = "builds"
    return [Build(api, b) for b in buildkite_get_json(api, endpoint; params = query_params)]
end

build_state(b::Build) = b.data["state"]

end


text:
"""
    AbsAttrib{T} <: AbstractArray{T,1}

Abstract type of attribute. It is a subtype of AbstractArray.
"""
abstract type AbsAttrib{T} <: AbstractArray{T,1} end

Base.IndexStyle(::Type{<:AbsAttrib}) = IndexLinear()

"""
    VecAttrib{T} <: AbsAttrib{T}

Attribute stores a quantity of type `T`, one value per shape.
The value is indexed to the serial number of the object.

# Example
Attribute to allow access to the locations of the vertices.
```
using Test
using StaticArrays
using MeshCore: VecAttrib, ShapeColl, P1, attribute
xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0]
N, T = size(xyz, 2), eltype(xyz)
locs =  VecAttrib([SVector{N,T}(xyz[i, :]) for i in 1:size(xyz, 1)])
vertices = ShapeColl(P1, size(xyz, 1))
vertices.attributes["geom"] = locs
a = attribute(vertices, "geom")
@test a[2] == [633.3333333333334, 0.0]
```
Attribute to label the vertices.
```
using Test
using StaticArrays
using MeshCore: VecAttrib, ShapeColl, P1, attribute
a = VecAttrib([1 for i in 1:size(xyz, 1)])
vertices = ShapeColl(P1, size(xyz, 1))
vertices.attributes["label"] = a
@test a[3] == 1
```

The vector attribute supports the methods `size`, `getindex`, `setindex!`.
"""
struct VecAttrib{T} <: AbsAttrib{T}
    v::Vector{T}
end

Base.size(a::A) where {A<:VecAttrib} = (length(a.v),)
Base.getindex(a::A, i::Int) where {A<:VecAttrib} = a.v[i]
Base.getindex(a::A, I::Vararg{Int,N}) where {A<:VecAttrib,N} = a.v[[I...]]
Base.setindex!(a::A, v, i::Int) where {A<:VecAttrib} = (a.v[i] = v)
Base.setindex!(a::A, v, I::Vararg{Int,N}) where {A<:VecAttrib,N} = (a.v[[I...]] .= v)

"""
    FunAttrib{T,M,F} <: AbsAttrib{T}

Attribute provides access to a quantity of type `T`, one value per shape.
The value is indexed to the serial number of the object.

# Example
Attribute to allow access to the locations of the vertices.
```
using Test
using StaticArrays
using MeshCore: VecAttrib, ShapeColl, P1, attribute
xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0]
N, T = size(xyz, 2), eltype(xyz)
flocs = FunAttrib(SVector{N,T}([0.0 0.0]), size(xyz, 1), i -> SVector{N,T}(xyz[i, :]))
@test flocs[10] == [633.3333333333334, 800.0]
```
Attribute to label the vertices.
```
using Test
using StaticArrays
using MeshCore: VecAttrib, ShapeColl, P1, attribute
xyz = [0.0 0.0; 633.3333333333334 0.0; 1266.6666666666667 0.0]
a = FunAttrib(0, size(xyz, 1), i -> i)
vertices = ShapeColl(P1, size(xyz, 1))
vertices.attributes["label"] = a
@test a[3] == 1
```

The function attribute supports the methods `size`, `getindex`.
It does not support `setindex!`: there is no way to change the data.
"""
struct FunAttrib{T,M,F} <: AbsAttrib{T}
    z::T
    m::M
    f::F
end

Base.size(a::A) where {A<:FunAttrib} = (a.m,)
Base.getindex(a::A, i::Int) where {A<:FunAttrib} = a.f(i)
Base.getindex(a::A, I::Vararg{Int,N}) where {A<:FunAttrib,N} = a.f.([I...])
# The `show` method may need some more work. It doesn't seem to play nice with @show .
Base.show(io::IOContext{Base.TTY}, unused::MIME{Symbol("text/plain")}, a::X) where {X<:FunAttrib} = let
    print(io, "$(a.m)-element FunAttrib")
end

"""
    datavaluetype(AbsAttrib{T}) where {T}

Retrieve the type of the attribute  value.
"""
datavaluetype(a::AbsAttrib{T}) where {T} = T


text:
################
# Script for getting MNIST data, selecting ones and sevens, shuffling
# and storing that subset in the data folder.
#
# Run this script standalone (instead of using your dev kernel and without
# using the .juliarc.jl of this project), like
# cd data
# julia preparemnist.jl
################

# Pkg.add("MNIST")
import MNIST

function preparemnist(Xy, func1, func2, p_X, p_y, maxclasssize, seed)
    X, y = Xy

    srand(seed)

    # split into classes
    X_a = X[:,func1(y)]    # X[:,(y.==7.0)]
    X_b = X[:,func2(y)]    # X[:,(y.!=7.0)]

    # cleanup
    X = nothing
    y = nothing
    gc()

    classsize = min(size(X_a)[2], size(X_b)[2], maxclasssize)
    println("Stratifying data to have 2 * $classsize examples (dimension $(size(X_a)[1]))...")

    # subsample
    X_a = X_a[:, randperm(size(X_a)[2])[1:classsize]]
    X_b = X_b[:, randperm(size(X_b)[2])[1:classsize]]

    # normalize and transpose for convenience
    X_a = X_a' ./ 255.0
    X_b = X_b' ./ 255.0

    # bring together
    y = vcat(repeat([1.0], inner=[classsize]), repeat([-1.0], inner=[classsize]))
    X = vcat(X_a, X_b)
    X_a = nothing
    X_b = nothing
    gc()

    println("Shuffling...")
    perm = randperm(2*classsize)
    X = X[perm,:]
    y = y[perm]

    println("Storing data into $p_X and $p_y...")
    writedlm(p_X, X)
    writedlm(p_y, y)
end

################
# folders
################
for folder in ["./seven_vs_all", "../results", "../results/seven_vs_all",
               "./seven_vs_six", "../results/seven_vs_six"]
    try
        mkdir(folder)
    end
end

################
# splits
################
for seed in 1337:(1337+9)
    println("\n", repeat("-", 10), "Seed $seed", repeat("-", 10))
    try
        mkdir("./seven_vs_all/$seed")
    end
    try
        mkdir("./seven_vs_six/$seed")
    end
    try
        mkdir("../results/seven_vs_all/$seed")
    end
    try
        mkdir("../results/seven_vs_six/$seed")
    end

    println("\nPreparing 7-vs-all TEST data...")
    preparemnist(MNIST.traindata(), (y -> y .== 7.0), (y -> y .!= 7.0),
                 "./seven_vs_all/$seed/X_test.dlm", "./seven_vs_all/$seed/y_test.dlm", 2500, seed)
    gc()

    println("\nPreparing 7-vs-all TRAINING data...")
    preparemnist(MNIST.testdata(), (y -> y .== 7.0), (y -> y .!= 7.0),
                 "./seven_vs_all/$seed/X_train.dlm", "./seven_vs_all/$seed/y_train.dlm", 1000, seed)
    gc()

    println("\nPreparing 7-vs-six TEST data...")
    preparemnist(MNIST.traindata(), (y -> y .== 7.0), (y -> y .== 6.0),
                 "./seven_vs_six/$seed/X_test.dlm", "./seven_vs_six/$seed/y_test.dlm", 2500, seed)
    gc()

    println("\nPreparing 7-vs-six TRAINING data...")
    preparemnist(MNIST.testdata(), (y -> y .== 7.0), (y -> y .== 6.0),
                 "./seven_vs_six/$seed/X_train.dlm", "./seven_vs_six/$seed/y_train.dlm", 1000, seed)
    gc()
end

println("\nDone.\n")

return nothing

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
using OCToolbox

using LinearAlgebra
using QuantumInformation

i2 = Matrix{Complex{Float64}}(I, 2, 2)

# we define our initial and final states
Ψ = [1+0.0im, 0.0]
ρ = [0+0.0im, 1.0]

N_ensemble = 10
ΔRange = 10  * 2π
H_drift = collect(range(-ΔRange, ΔRange, length = N_ensemble)) .* [sz]
H_ctrl = [π * sx, π * sy]

# we set up a functional for a robust pulse
function fn(controls)
    controls = complex.(real.(controls))
    err = 0
    for i = 1:N_ensemble
        U = pw_full_evolution(H_drift[i], H_ctrl, controls, Δt, i2)
        err += C2(ρ, U * Ψ)
    end
    err
end

K = 2
N = 20
T = 10
Δt = T/N

control_guess = rand(K, N).*0.001
fn(control_guess)

using Zygote
Zygote.gradient(fn, control_guess)

o = GRAPE(fn, control_guess, K, N)

fn(o.minimizer)


text:
#Based upon alg 2 in https://ieeexplore.ieee.org/document/8573778

"""
$(TYPEDEF)
`SR3` is an optimizer framework introduced [by Zheng et. al., 2018](https://ieeexplore.ieee.org/document/8573778) and used within
[Champion et. al., 2019](https://arxiv.org/abs/1906.10612). `SR3` contains a sparsification parameter `λ`, a relaxation `ν`.
It solves the following problem
```math
\\argmin_{x, w} \\frac{1}{2} \\| Ax-b\\|_2 + \\lambda R(w) + \\frac{\\nu}{2}\\|x-w\\|_2
```
Where `R` is a proximal operator and the result is given by `w`.

# Fields
$(FIELDS)

# Example
```julia
opt = SR3()
opt = SR3(1e-2)
opt = SR3(1e-3, 1.0)
opt = SR3(1e-3, 1.0, SoftThreshold())
```
## Note
Opposed to the original formulation, we use `ν` as a relaxation parameter,
as given in [Champion et. al., 2019](https://arxiv.org/abs/1906.10612). In the standard case of
hard thresholding the sparsity is interpreted as `λ = threshold^2 / 2`, otherwise `λ = threshold`.
"""
mutable struct SR3{T, V, P <: AbstractProximalOperator} <: AbstractOptimizer{T}
   """Sparsity threshold"""
   λ::T
   """Relaxation parameter"""
   ν::V
   """Proximal operator"""
   R::P

   function SR3(threshold::T = 1e-1, ν::V = 1.0, R::P = HardThreshold()) where {T,V <: Number,P <: AbstractProximalOperator}
      @assert all(threshold .> zero(eltype(threshold))) "Threshold must be positive definite"
      @assert ν > zero(V) "Relaxation must be positive definite"

       λ = isa(R, HardThreshold) ? threshold.^2 /2 : threshold
       return new{typeof(λ), V, P}(λ, ν, R)
   end

   function SR3(threshold::T , R::P) where {T,P <: AbstractProximalOperator}
      @assert all(threshold .> zero(eltype(threshold))) "Threshold must be positive definite"
       λ = isa(R, HardThreshold) ? threshold.^2 /2 : threshold
       ν = one(eltype(λ))
       return new{typeof(λ), eltype(λ), P}(λ, ν, R)
   end
end

Base.summary(::SR3) = "SR3"

function (opt::SR3{T,V,R})(X, A, Y, λ::V = first(opt.λ);
    maxiter::Int64 = maximum(size(A)), abstol::V = eps(eltype(T)), progress = nothing, kwargs...)  where {T, V, R}

   n, m = size(A)
   ν = opt.ν
   W = copy(X)

   # Init matrices
   H = A'*A+I(m)*ν
   H = cholesky!(H)
   X̂ = A'*Y

   w_i = similar(W)
   w_i .= W
   iters = 0

   iters = 0
   converged = false

   xzero = zero(eltype(X))
   obj = xzero
   sparsity = xzero
   conv_measure = xzero

   _progress = isa(progress, Progress)
   initial_prog = _progress ? progress.counter : 0

   @views while (iters < maxiter) && !converged
       iters += 1

       # Solve ridge regression
       X .= H \ (X̂ .+ W*ν)
       #ldiv!(X, H, X̂ .+ W*ν)
       # Proximal
       opt.R(W, X, λ)

       conv_measure = norm(w_i .- W, 2)

       if _progress
           obj = norm(Y - A*X, 2)
           sparsity = norm(X, 0, λ)

           ProgressMeter.next!(
           progress;
           showvalues = [
               (:Threshold, λ), (:Objective, obj), (:Sparsity, sparsity),
               (:Convergence, conv_measure)
           ]
           )
       end

       if conv_measure < abstol
           converged = true

           #if _progress

            #   ProgressMeter.update!(
            #   progress,
            #   initial_prog + maxiter
            #   )
           #end

       else
           w_i .= W
       end
   end
   # We really search for W here
   X .= W
   @views clip_by_threshold!(X, λ)
   return
end


text:
using TensorOperations

function make4(a1,a2,a3,a4,a12,a34,a1234)
   t1 = HTensorNode(a1)
   t2 = HTensorNode(a2)
   t3 = HTensorNode(a3)
   t4 = HTensorNode(a4)
   t12 = HTensorNode(a12, [t1,t2])
   t34 = HTensorNode(a34, [t3,t4])
   t1234 = HTensorNode(a1234, [t12,t34])
   return t1234
end

# B[i,j,k] = A[a,b,k] U[i,a] V[j,b]
function AxUV(A, U, V)
   T = eltype(A)
   na, nb, nk = size(A)
   ni = size(U,1)
   nj = size(V,1)
   @assert size(U,2) == na
   @assert size(V,2) == nb
   tmp = Array{T}(ni, nb)
   B = zeros(T, ni, nj, nk)
   for k=1:nk
      fill!(tmp, zero(T))
      for b=1:nb, a=1:na, i=1:ni
         tmp[i,b] += A[a,b,k] * U[i,a]
      end
      for b=1:nb, j=1:nj, i=1:ni
         B[i,j,k] += tmp[i,b] * V[j,b]
      end
   end
   return B
end

function cost(u1,u2,u3,u4,u12,u34,u1234,target)
   a12 = AxUV(u12, u1, u2)
   a34 = AxUV(u34, u3, u4)
   n1, n2, n12 = size(a12)
   n3, n4, n34 = size(a34)
   a1234 = AxUV( u1234, reshape(a12, n1*n2, n12), reshape(a34, n3*n4, n34))
   guess = reshape(a1234, n1, n2, n3, n4, 1)
   return sum(abs2, guess-target)
end


text:
# Sequential assignment. Multiple locations

## ------------- Switches

function StructLH.describe(switches :: EntryDecisionSwitches{F1}) where F1
    nLoc = n_locations(switches);
    capStr = limited_capacity(switches)  ?  "limited"  :  "unlimited";
    return [
        "College Entry"  "Sequential. $nLoc locations";
        "Capacities"  capStr]
end

"""
    $(SIGNATURES)

Constructor for `EntryDecisionSwitches`: One location. Equal type mass. No capacity constraints by default.
"""
function make_entry_switches_oneloc(J, nc; entryPrefScale = 1.0,
    typeMass = 1.0, capacity_cV = nothing)

    F1 = typeof(entryPrefScale);
    if isnothing(capacity_cV)
        capacity_clM = fill(F1(CapacityInf), nc, 1);
    else
        capacity_clM = matrix_from_vector(capacity_cV);
    end
    switches = EntryDecisionSwitches{F1}(
        nTypes = J, nColleges = nc, nLocations = 1,
        typeMass_jlM = fill(F1(typeMass), J, 1),
        capacity_clM = capacity_clM,
        entryPrefScale = entryPrefScale,
        valueLocal = zero(F1), calValueLocal = false);
    @assert validate_es(switches)
    return switches
end

function validate_es(switches :: EntryDecisionSwitches{F1}) where F1
    isValid = true;
    J = n_types(switches);
    nc = n_colleges(switches);
    nl = n_locations(switches);
    if nl == 1
        if switches.calValueLocal  ||  !(switches.valueLocal == zero(F1))
            isValid = false;
            @warn  "Cannot have local value with one location"
        end
    end
    return isValid
end

## ---------------  Entry Decision constructor

function init_entry_decision(objId :: ObjectId,
    switches :: EntryDecisionSwitches{F1},
    st :: SymbolTable) where F1

    pEntryPref = init_entry_prefscale(switches, st);
    pValueLocal = init_value_local(switches, st);
    pvec = ParamVector(objId, [pEntryPref, pValueLocal]);
    return EntryDecision(objId, pvec, ModelParams.value(pEntryPref),
        ModelParams.value(pValueLocal), switches)
end

function init_value_local(switches :: EntryDecisionSwitches{F1},
    st :: SymbolTable) where F1
    if n_locations(switches) == 1
        # Cannot have local value
        valueLocal = zero(F1);
        calValueLocal = false;
    else
        valueLocal = switches.valueLocal;
        calValueLocal = switches.calValueLocal;
    end

    ub = max(F1(5.0), valueLocal + 2.0);
    pEntryPref = Param(:valueLocal, LatexLH.description(st, :uLocal), latex(st, :uLocal), valueLocal, valueLocal, F1(0.0), ub,
        calValueLocal);
    return pEntryPref
end

## ------------  Entry decisions

"""
    $(SIGNATURES)

Compute entry probabilities and expected values at entry from admission rule and entry decision objects.

It is possible to solve for a subset of types. Only those listed in `rank_jV` will be solved. But keep in mind that capacity constraints may then not bind.

The type dimension of the `EntryResults` object matches that of the `EntryDecision`.

Optional: solve each student's entry decision without preference shocks, but continue to update enrollments with preference shocks. This is useful to diagnose how important preference shocks are.

Note that multiple locations only matter if they are not identical. If all colleges are available in all locations (and not full), the fraction going local (conditional on entry) is only a function of the number of locations and the value of going local. `vCollege_jcM` and `vWork_jV` do not matter. This is easy to check analytically.
"""
function entry_decisions(entryS :: EntryDecision{F1},
    admissionS :: AbstractAdmissionsRule{I1, F1},
    admProbFct :: AF1,
    vWork_jV :: Vector{F1}, vCollege_jcM :: Matrix{F1},
    endowPctV :: Vector{F1},
    rank_jV :: Vector{I2};
    prefShocks :: Bool = true) where
    {AF1 <: AbstractAdmProbFct{<: Real}, I1, I2 <: Integer, F1}

    # nTypes = n_types(entryS);
    nc = n_colleges(entryS);
    nl = n_locations(entryS);
    er = EntryResults(entryS.switches);

    # Loop over students in order of ranking
    for j in rank_jV
        # Loop over student locations in ascending order
        # We are now processing student (j, l)
        for l = 1 : nl
            full_clM = (enrollment_cl(er) .>= capacities(entryS));

            # This is the standard one-step entry decision, but with colleges from all locations stacked.
            entryProb_clM, er.eVal_jlM[j, l], entryProbBest_clM =
                entry_decisions_one_student(
                    entryS, admissionS, admProbFct,
                    vWork_jV[j], vCollege_jcM[j,:],
                    endowPctV[j], full_clM, l; prefShocks = true);

            # Record enrollment
            typeMass = type_mass_jl(entryS, j, l);
            er.enroll_clM .+= typeMass .* entryProb_clM;
            er.enrollLocal_clM[:,l] .+= typeMass .* entryProb_clM[:,l];
            er.enrollBest_clM .+= typeMass .* entryProbBest_clM;

            if !prefShocks
                # Solve again without pref shocks
                # But enrollment is determined by original problem
                entryProb_clM, er.eVal_jlM[j, l], entryProbBest_clM =
                    entry_decisions_one_student(
                        entryS, admissionS, admProbFct,
                        vWork_jV[j], vCollege_jcM[j,:],
                        endowPctV[j], full_clM, l; prefShocks = false);
            end

            # Record entry probs
            entryProb_cV = vec(sum(entryProb_clM, dims = 2));
            make_valid_probs!(entryProb_cV);
            # Fraction entering any `c` college; local or not.
            er.fracEnter_jlcM[j,l,:] .= entryProb_cV;

            # Probability of entering college `c` as the best college
            entryProbBest_cV = vec(sum(entryProbBest_clM, dims = 2));
            make_valid_probs!(entryProbBest_cV);
            er.fracEnterBest_jlcM[j,l,:] .= entryProbBest_cV;
            # Fraction entering a local `c` college.
            er.fracLocal_jlcM[j,l,:] = entryProb_clM[:,l];
        end
    end

    @assert validate_er(er; validateFracLocal = prefShocks)
    return er
end

# # College enrollments of student `j` in `l`
# function enrollment_cl(entryS :: EntryDecision{F1},
#     entryProb_clM :: Matrix{F1}, j :: Integer, l :: Integer) where F1

#     return entryProb_clM .* type_mass_jl(entryS, j, l)
# end

# --------------


text:
#
#	gaussian
#
# process:
gaus_rf(ν::T where T<:Real) = exp(-ν^2/2)
# field:
gaus_rf(V::Vector{T} where T<:Real) = gaus_rf(norm(V))

#
#	Rect / sinc
#
# process:
rect_rf(ν::T where T<:Real) = abs(ν) < 1 ? 1 : 0
sinc_rf(t::T where T<:Real) = (2/sqrt(2pi)) * (t!=0 ? sin(t)/t : 1)
# field:
rect_rf(V::Vector{T} where T<:Real) = rect_rf(norm(V))
sinc_rf(V::Vector{T} where T<:Real) = sinc_rf(norm(V))


text:
using Test
using Random, SparseArrays, LinearAlgebra, StableRNGs
using OpenQuantumSystems
import DelayDiffEq
import QuantumOpticsBase

@testset "master ansatz" begin

    Random.seed!(StableRNG(0), 1)

    D(op1::Array, op2::Array) = abs(norm(op1 - op2))
    D(x1::StateVector, x2::StateVector) = norm(x2 - x1)
    D(op1::AbstractOperator, op2::AbstractOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))
    D(op1::AbstractSuperOperator, op2::AbstractSuperOperator) =
        abs(tracedistance_nh(dense(op1), dense(op2)))

    # TODO: change to macro
    HR = 0.01
    shift = (2.0 * HR)^0.5
    mode1 = Mode(300., shift)
    mol1 = Molecule([mode1], 3, [12500., 12750.])
    mol2 = Molecule([mode1], 3, [12500., 12800.])
    agg = Aggregate([mol1, mol2])
    aggInds = getIndices(agg)
    vibindices = getVibIndices(agg, aggInds)
    aggIndsLen = length(aggInds)
    basis = GenericBasis([aggIndsLen])
    FCFact = getFranckCondonFactors(agg, aggInds)
    FCProd = getFCProd(agg, FCFact, aggInds, vibindices)

    Ham_sys = getAggHamSystemSmall(agg)
    Ham_I = getAggHamInteraction(agg, aggInds, FCFact)
    Ham_0 = getAggHamSystemBath(agg, aggInds, FCFact)
    Ham_B = getAggHamBathBig(agg)
    Ham = getAggHamiltonian(agg, aggInds, FCFact)

    Ham_0_lambda, Ham_0_S = eigen(Ham_0.data)
    Ham_0_Sinv = inv(Ham_0_S)
    Ham_0_lambda = diagm(Ham_0_lambda)

    data = Matrix(Hermitian(rand(ComplexF64, aggIndsLen, aggIndsLen)))
    rho0 = DenseOperator(basis, basis, data)
    normalize!(rho0)
    # tests have to be quick enough
    t_max = 0.001
    t_count = 10
    t0 = 0.
    t_step = (t_max - t0) / (t_count)
    tspan = [t0:t_step:t_max;]

    T = 300
    mu_array = [[2, 1]]
    W01 = thermal_state(T, mu_array, Ham, vibindices, aggInds; diagonalize=false, diagonal=true)
    rho_traced1 = trace_bath(W01.data, agg, FCProd, aggInds, vibindices)
    mu_array = [[1, 2]]
    W02 = thermal_state(T, mu_array, Ham, vibindices, aggInds; diagonalize=false, diagonal=true)
    rho_traced2 = trace_bath(W02.data, agg, FCProd, aggInds, vibindices)

    W0 = 1.0*W01 + 0.0*W02
    W0_bath = get_rho_bath(W0, agg, FCProd, aggInds, vibindices)
    rho0 = trace_bath(W0, agg, FCProd, aggInds, vibindices)

    W0 = DenseOperator(W0.basis_l, W0.basis_r, complex(W0.data))
    W0_bath = DenseOperator(W0_bath.basis_l, W0_bath.basis_r, complex(W0_bath.data))
    rho0 = DenseOperator(rho0.basis_l, rho0.basis_r, complex(rho0.data))

    p = (Ham_0, Ham_I, Ham_0_lambda, Ham_0_S, Ham_0_Sinv, Ham_B, W0, W0_bath, agg, FCProd, aggInds, vibindices, ComplexF64)
    Tspan, rho_t = master_ansatz(
        rho0,
        tspan,
        p;
        reltol = 1.0e-4,
        abstol = 1.0e-4,
        alg = DelayDiffEq.MethodOfSteps(DelayDiffEq.Tsit5())
    )
    rho_prev = deepcopy(rho0)
    for t_i = 2:length(rho_t)
        t = Tspan[t_i]
        rho_I = rho_t[t_i]
        U_op_S = evolutionOperator(Ham_sys, t)
        rho = U_op_S * rho_I * U_op_S'
    end

end # testset


text:
module CineFiles

using ColorTypes: Color, N0f8, N6f10, N4f12, N0f16, gray, Gray, AbstractGray
using ColorVectorSpace
using LRUCache
using StaticArrays: SizedMatrix, SizedVector
using FixedPointNumbers: Normed, reinterpret, nbitsfrac

import Base: eltype, length, size, getindex, firstindex, lastindex, iterate, read, read!
export CineFile

include("LUT.jl")

abstract type RawFrame end
abstract type BinaryData end

rawcounts(v::AbstractGray) = reinterpret(gray(v))

struct Levels{T}
    white::T
    black::T
end

function Levels{T}(W::S, B::S) where {T<:AbstractGray, S<:Int32}
    scale = rawcounts(T(1))
    Levels{T}(T(W / scale), T(B / scale))
end

struct Unpacked{T} <: RawFrame
    width::Int
    height::Int
    tmp::Array{T,2}
    unpacked::Array{T,2}
    levels::Levels
    Unpacked{T}(width, height, whiteL, blackL) where {T<:Color} = new{T}(
        width,
        height,
        Array{T,2}(undef, width, height),
        Array{T,2}(undef, width, height),
        Levels{T}(whiteL, blackL),
    )
end

struct Packed{T,S} <: RawFrame
    width::Int
    height::Int
    tmp::Array{S,2}
    unpacked::Array{T,2}
    pack::Vector{UInt8}
    levels::Levels
end

function Packed{T,S}(width, height, whiteL, blackL) where {T<:Gray,S<:Gray}
    true_bits = nbitsfrac(eltype(T)) * width * height
    true_bytes = ceil(Int, true_bits / 8)
    Packed{T,S}(
        width,
        height,
        Array{S,2}(undef, width, height),
        Array{T,2}(undef, width, height),
        Vector{UInt8}(undef, true_bytes),
        Levels{T}(whiteL, blackL),
    )
end

function Packed{T}(width, height, whiteL, blackL) where {T<:Gray}
    true_bits = nbitsfrac(eltype(T)) * width * height
    true_bytes = ceil(Int, true_bits / 8)
    Packed{T,T}(
        width,
        height,
        Array{T,2}(undef, width, height),
        Array{T,2}(undef, width, height),
        Vector{UInt8}(undef, true_bytes),
        Levels{T}(whiteL, blackL),
    )
end

struct Time64 <: BinaryData
    fractions::UInt32
    seconds::UInt32
end

struct WBgain <: BinaryData
    R::Float32
    B::Float32
end

struct IMfilter <: BinaryData
    Dim::Int32
    Shifts::Int32
    Bias::Int32
    Coef::SizedMatrix{5,5,Int32,2,Matrix{Int32}}
end

struct Rect <: BinaryData
    Left::Int32
    Top::Int32
    Right::Int32
    Bottom::Int32
end

struct TimeCode <: BinaryData
    framesU::UInt8
    framesT::UInt8
    dropFrameFlag::UInt8
    colorFrameFlag::UInt8
    secondsU::UInt8
    secondsT::UInt8
    flag1::UInt8
    minutesU::UInt8
    minutesT::UInt8
    flag2::UInt8
    hoursU::UInt8
    hoursT::UInt8
    flag3::UInt8
    flag4::UInt8
    userBitData::UInt32
end

struct SetupHeader <: BinaryData
    FrameRate16::UInt16
    Shutter16::UInt16
    PostTrigger16::UInt16
    FrameDelay16::UInt16
    AspectRatio::UInt16
    Res7::UInt16
    Res8::UInt16
    Res9::UInt8
    Res10::UInt8
    Res11::UInt8
    TrigFrame::UInt8
    Res12::UInt8
    DescriptionOld::SizedVector{121,Int8}
    Mark::UInt16
    Length::UInt16
    Res13::UInt16
    SigOption::UInt16
    BinChannels::Int16
    SamplesPerImage::UInt8
    BinName::SizedMatrix{8,11,Int8}
    AnaOption::UInt16
    AnaChannels::Int16
    Res6::UInt8
    AnaBoard::UInt8
    ChOption::SizedVector{8,Int16}
    AnaGain::SizedVector{8,Float32}
    AnaUnit::SizedMatrix{8,6,Int8}
    AnaName::SizedMatrix{8,11,Int8}
    lFirstImage::Int32
    dwImageCount::UInt32
    nQFactor::Int16
    wCineFileType::UInt16
    szCinePath::SizedMatrix{4,65,Int8}
    Res14::UInt16
    Res15::UInt8
    Res16::UInt8
    Res17::UInt16
    Res18::Float64
    Res19::Float64
    Res20::UInt16
    Res1::Int32
    Res2::Int32
    Res3::Int32
    ImWidth::UInt16
    ImHeight::UInt16
    EDRShutter16::UInt16
    Serial::UInt32
    Saturation::Int32
    Res5::UInt8
    AutoExposure::UInt32
    bFlipH::Int32
    bFlipV::Int32
    Grid::UInt32
    FrameRate::UInt32
    Shutter::UInt32
    EDRShutter::UInt32
    PostTrigger::UInt32
    FrameDelay::UInt32
    bEnableColor::Int32
    CameraVersion::UInt32
    FirmwareVersion::UInt32
    SoftwareVersion::UInt32
    RecordingTimeZone::Int32
    CFA::UInt32
    Bright::Int32
    Contrast::Int32
    Gamma::Int32
    Res21::UInt32
    AutoExpLevel::UInt32
    AutoExpSpeed::UInt32
    AutoExpRect::Rect
    WBGain::SizedVector{4,WBgain}
    Rotate::Int32
    WBView::WBgain
    RealBPP::UInt32
    Conv8Min::UInt32
    Conv8Max::UInt32
    FilterCode::Int32
    FilterParam::Int32
    UF::IMfilter
    BlackCalSVer::UInt32
    WhiteCalSVer::UInt32
    GrayCalSVer::UInt32
    bStampTime::Int32
    SoundDest::UInt32
    FRPSteps::UInt32
    FRPImgNr::SizedVector{16,Int32}
    FRPRate::SizedVector{16,UInt32}
    FRPExp::SizedVector{16,UInt32}
    MCCnt::Int32
    MCPercent::SizedVector{64,Float32}
    CICalib::UInt32
    CalibWidth::UInt32
    CalibHeight::UInt32
    CalibRate::UInt32
    CalibExp::UInt32
    CalibEDR::UInt32
    CalibTemp::UInt32
    HeadSerial::SizedVector{4,UInt32}
    RangeCode::UInt32
    RangeSize::UInt32
    Decimation::UInt32
    MasterSerial::UInt32
    Sensor::UInt32
    ShutterNs::UInt32
    EDRShutterNs::UInt32
    FrameDelayNs::UInt32
    ImPosXAcq::UInt32
    ImPosYAcq::UInt32
    ImWidthAcq::UInt32
    ImHeightAcq::UInt32
    Description::SizedVector{4096,Int8}
    RisingEdge::Int32
    FilterTime::UInt32
    LongReady::Int32
    ShutterOff::Int32
    Res4::SizedVector{16,UInt8}
    bMetaWB::Int32
    Hue::Int32
    BlackLevel::Int32
    WhiteLevel::Int32
    LensDescription::SizedVector{256,Int8}
    LensAperture::Float32
    LensFocusDistance::Float32
    LensFocalLength::Float32
    fOffset::Float32
    fGain::Float32
    fSaturation::Float32
    fHue::Float32
    fGamma::Float32
    fGammaR::Float32
    fGammaB::Float32
    fFlare::Float32
    fPedestalR::Float32
    fPedestalG::Float32
    fPedestalB::Float32
    fChroma::Float32
    ToneLabel::SizedVector{256,Int8}
    TonePoints::Int32
    fTone::SizedVector{64,Float32}
    UserMatrixLabel::SizedVector{256,Int8}
    EnableMatrices::Int32
    fUserMatrix::SizedVector{9,Float32}
    EnableCrop::Int32
    CropRect::Rect
    EnableResample::Int32
    ResampleWidth::UInt32
    ResampleHeight::UInt32
    fGain16_8::Float32
    FRPShape::SizedVector{16,UInt32}
    TrigTC::TimeCode
    fPbRate::Float32
    fTcRate::Float32
    CineName::SizedVector{256,UInt8}
end

struct CineFileHeader <: BinaryData
    Type::UInt16
    HeaderSize::UInt16
    Compression::UInt16
    Version::UInt16
    FirstMovieImage::Int32
    TotalImageCount::UInt32
    FirstImageNum::Int32
    ImageCount::UInt32
    OffImageHeader::UInt32
    OffSetup::UInt32
    OffImageOffsets::UInt32
    TriggerTime::Time64
end

struct BitmapInfoHeader <: BinaryData
    Size::UInt32
    Width::Int32
    Height::Int32
    Planes::UInt16
    BitCount::UInt16
    Compression::UInt32
    SizeImage::UInt32
    XPelsPerMeter::Int32
    YPelsPerMeter::Int32
    ClrUsed::UInt32
    ClrImportant::UInt32
    FPS::UInt16
end

struct CineHeader{T}
    cine::CineFileHeader
    bitmap::BitmapInfoHeader
    setup::SetupHeader
    imglocs::Vector{Int64}
    imgoffset::Int
    dt::Vector{Float64}
    raw::RawFrame
end

struct CineFile{T}
    path::String
    header::CineHeader{T}
    data::LRU{Int,Array{T,2}}
end

function unpack!(pack::Vector{UInt8}, unpacked::Array{Gray{N6f10}})
    pack16 = UInt16.(pack)
    indx = CartesianIndices(unpacked)[:, end:-1:1]
    unpacked[indx[1:4:end]] =
        reinterpret.(N6f10, (pack16[1:5:end] .<< 2) .| (pack16[2:5:end] .>> 6))
    unpacked[indx[2:4:end]] =
        reinterpret.(
            N6f10,
            ((pack16[2:5:end] .& 0b00111111) .<< 4) .| (pack16[3:5:end] .>> 4),
        )
    unpacked[indx[3:4:end]] =
        reinterpret.(
            N6f10,
            ((pack16[3:5:end] .& 0b00001111) .<< 6) .| (pack16[4:5:end] .>> 2),
        )
    unpacked[indx[4:4:end]] =
        reinterpret.(N6f10, ((pack16[4:5:end] .& 0b00000011) .<< 8) .| pack16[5:5:end])
end

function unpack!(pack::Vector{UInt8}, unpacked::Array{Gray{N4f12}})
    @warn "unpacking 12 bit images untested"
    pack16 = UInt16.(pack)
    indx = CartesianIndices(unpacked)[:, end:-1:1]
    unpacked[indx[1:2:end]] =
        reinterpret.(N6f10, (pack16[1:3:end] .<< 4) .| (pack16[2:3:end] .>> 4))
    unpacked[indx[2:2:end]] =
        reinterpret.(N6f10, ((pack16[2:3:end] .& 0b00001111) .<< 8) .| (pack16[3:3:end]))
end

function linearize!(raw_data::RawFrame)
    Wlevel = raw_data.levels.white
    Blevel = raw_data.levels.black
    raw_data.tmp .= raw_data.unpacked
    raw_data.tmp[raw_data.tmp.>Wlevel] .= Wlevel
    raw_data.tmp[raw_data.tmp.<Blevel] .= Blevel
    raw_data.tmp .= (raw_data.tmp .- Blevel) ./ (Wlevel - Blevel)
end

function linearize!(raw_data::Packed{Gray{N6f10}})
    Blevel, Wlevel = Gray{N4f12}.((64, 4064) ./ (2^12))
    raw_data.tmp .= lookup.(raw_data.unpacked, Ref(CINE_LUT))
    raw_data.tmp[raw_data.tmp.>Wlevel] .= Wlevel
    raw_data.tmp[raw_data.tmp.<Blevel] .= Blevel
    raw_data.tmp .= (raw_data.tmp .- Blevel) ./ (Wlevel - Blevel)
end

function CineHeader(fname)
    open(fname) do f
        read(f, UInt16) == UInt(18755) || error(basename(fname), " is not a .cine file")
        seek(f, 0)
        cine = read(f, CineFileHeader)
        seek(f, cine.OffImageHeader)
        bitmap = read(f, BitmapInfoHeader)

        seek(f, cine.OffSetup)
        setup = read(f, SetupHeader)

        bitsized = (bitmap.BitCount == 8) || (bitmap.BitCount == 24) ? UInt8 : UInt16

        if (bitmap.BitCount == 8) || (bitmap.BitCount == 16)
            bittype = Gray{Normed{bitsized,Int(setup.RealBPP)}}
        elseif (bitmap.BitCount == 24) || (bitmap.BitCount == 36)
            bittype = BGR{Normed{bitsized,Int(setup.RealBPP)}}
        end

        cine.ImageCount > 0 || error("no images exist in file")

        seek(f, cine.OffImageOffsets)
        imglocs = read!(f, Array{Int64}(undef, cine.ImageCount))

        seekstart(f)
        imgoffset = 1
        while read(f, UInt16) != 1002
            imgoffset += 1
        end

        dt = zeros(cine.ImageCount)
        skip(f, 2)
        for i in eachindex(dt)
            fracstart = read(f, UInt32)
            secstart = read(f, UInt32)
            dt[i] =
                (secstart - cine.TriggerTime.seconds) +
                ((fracstart - cine.TriggerTime.fractions) / 2^32)
        end

        # Packing options
        if bitmap.Compression == 0
            raw = Unpacked{bittype}(
                bitmap.Width,
                bitmap.Height,
                setup.WhiteLevel,
                setup.BlackLevel,
            )
            return CineHeader{bittype}(cine, bitmap, setup, imglocs, imgoffset, dt, raw)
        elseif bitmap.Compression == 256
            raw = Packed{bittype,Gray{N6f10}}(
                bitmap.Width,
                bitmap.Height,
                setup.WhiteLevel,
                setup.BlackLevel,
            )
            return CineHeader{Gray{N4f12}}(cine, bitmap, setup, imglocs, imgoffset, dt, raw)
        elseif bitmap.Compression == 1024
            raw = Packed{bittype,Gray{N4f12}}(
                bitmap.Width,
                bitmap.Height,
                setup.WhiteLevel,
                setup.BlackLevel,
            )
            return CineHeader{bittype}(cine, bitmap, setup, imglocs, imgoffset, dt, raw)
        end
    end
end

function readframe!(f::IO, frame, h, frameidx)
    1 <= frameidx <= h.cine.ImageCount || BoundsError(h.dt, frameidx)
    seek(f, h.imglocs[frameidx])
    skip(f, read(f, UInt32) - 4)
    read!(f, h.raw)
    linearize!(h.raw)
    # non-allocating version of rotl90
    ind1, ind2 = axes(h.raw.tmp)
    n = first(ind2) + last(ind2)
    for i in axes(h.raw.tmp, 1), j in ind2
        frame[n-j, i] = h.raw.tmp[i, j]
    end
    return frame
end

function read!(f::IOStream, frame::Packed)
    read!(f, frame.pack)
    unpack!(frame.pack, frame.unpacked)
end

read!(f::IOStream, frame::Unpacked) = read!(f, frame.unpacked)
# TODO: implement with lower compilation overhead
read(f::IOStream, S::Type{T}) where {T<:BinaryData} = S(read.(Ref(f), S.types)...)

readframe!(filename, frame, h, frameidx) =
    open(f -> readframe!(f, frame, h, frameidx), filename)
readframe(f, h, frameidx) =
    readframe!(f, similar(h.raw.tmp, size(h.raw.tmp, 2), size(h.raw.tmp, 1)), h, frameidx)

"""
    CineFile(filepath, cachelimit=0.25)

Load header information and create a frame cache for a Phantom .cine file.
    `cachelimit` sets the maximum size of the frame cache as a fraction of
    your system's free RAM. Indexing and iteration is supported for CineFiles.

"""
function CineFile(filepath, cachelimit = 0.25)
    header = CineHeader(filepath)
    framesize = Base.summarysize(header.raw.tmp)
    maxcachedframes = ceil(Int, cachelimit * Sys.free_memory() / framesize)
    data = LRU{Int,Array{eltype(header),2}}(maxsize = maxcachedframes)
    return CineFile(filepath, header, data)
end

Base.eltype(h::CineHeader{T}) where {T} = T
Base.length(cf::CineFile) = length(cf.header.dt)
Base.eltype(cf::CineFile) = typeof(cf.header.raw.tmp)
Base.size(cf::CineFile) = (length(cf), size(cf.header.raw.tmp, 2), size(cf.header.raw.tmp, 1))

function Base.getindex(cf::CineFile, idx::Int)
    get!(cf.data, idx) do
        readframe(cf.path, cf.header, idx)
    end
end

Base.firstindex(cf::CineFile) = 1
Base.lastindex(cf::CineFile) = length(cf)
Base.getindex(cf::CineFile, I) = [cf[i] for i in I]
Base.iterate(cf::CineFile, state = 1) = state > length(cf) ? nothing : (cf[state], state + 1)

# precompile as the final step of the module definition:
if ccall(:jl_generating_output, Cint, ()) == 1   # if we're precompiling the package
    let
        precompile(read, (IOStream, SetupHeader))

        # CineFile(joinpath(dirname(@__DIR__), "test/data/8bpp.cine"))[1]
    end
end

end


text:
using IBAN
using Base.Test

## valid IBANs
test_iban1 = "GB82 WEST 1234 5698 7654 32"

@test is_valid_iban(test_iban1)
@test is_valid_iban(replace(test_iban1, " ", ""))
@test is_valid_iban(replace(test_iban1, " ", ","))

### invalid IBANs

## wrong bban
test_wban1 = "GB82 WEST 2134 5698 7654 32"
@test !is_valid_iban(test_wban1)

## unknown country
test_wban2 = "GX82 WEST 1234 5698 7654 32"
@test !is_valid_iban(test_wban2)

## wrong bban length
test_wban3 = "GB82 WEST 1234 5698 7654 3200"
@test !is_valid_iban(test_wban3)

## wrong global re
test_wban4 = "GEXX WEST 1234 5698 7654 32"
@test !is_valid_iban(test_wban4)

## example fake, valid IBANs
valid_ibans = [
              "AL90 2081 1008 0000 0010 3953 1801",
              "AT13 1490 0220 1001 0999",
              "BE68 8440 1037 0034",
              "CZ42 0100 0000 1955 0503 0267",
              "EE34 2200 2210 3412 6658",
              "FI98 1428 3500 1711 41",
              "FR76 3006 6100 4100 0105 7380 116",
              "DE12 5001 0517 0648 4898 90",
              "HU29 1170 8001 2054 7794 0000 0000",
              "IE92 BOFI 9000 1710 0279 52",
              "IT68 D030 0203 2800 0040 0162 854",
              "LI10 0880 0000 0201 7630 6",
              "LU76 1111 0008 7296 0000",
              "MT98 MMEB 4409 3000 0000 0902 7293 051",
              "NL18 ABNA 0484 8698 68",
              "PT50 0035 0683 0000 0007 8431 1",
              "SE64 1200 0000 0121 7014 5230",
              "CH39 0870 4016 0754 7300 7",
              "SK93 1111 0000 0010 5736 1004",
              "SI56 0310 0100 1300 933",
              "ES10 2090 3200 5000 4104 5040",
              "GB96 MIDL 4027 1522 8598 82",
              "AD12 0001 2030 2003 5910 0100",
              "AT61 1904 3002 3457 3201",
              "BE68 5390 0754 7034",
              "CY17 0020 0128 0000 0012 0052 7600",
              "CZ65 0800 0000 1920 0014 5399",
              "DK50 0040 0440 1162 43",
              "EE38 2200 2210 2014 5685",
              "FI21 1234 5600 0007 85",
              "FR14 2004 1010 0505 0001 3M02 606",
              "DE89 3704 0044 0532 0130 00",
              "GI75 NWBK 0000 0000 7099 453",
              "GR16 0110 1250 0000 0001 2300 695",
              "HU42 1177 3016 1111 1018 0000 0000",
              "IS14 0159 2600 7654 5510 7303 39",
              "IE29 AIBK 9311 5212 3456 78",
              "IT60 X054 2811 1010 0000 0123 456",
              "LV80 BANK 0000 4351 9500 1",
              "LT12 1000 0111 0100 1000",
              "LU28 0019 4006 4475 0000",
              "NL91 ABNA 0417 1643 00",
              "NO93 8601 1117 947",
              "PL27 1140 2004 0000 3002 0135 5387",
              "PT50 0002 0123 1234 5678 9015 4",
              "SK31 1200 0000 1987 4263 7541",
              "SI56 1910 0000 0123 438",
              "ES80 2310 0001 1800 0001 2345",
              "SE35 5000 0000 0549 1000 0003",
              "CH39 0070 0115 2018 4917 3",
              "GB29 NWBK 6016 1331 9268 19"]

for iban in valid_ibans
    @test is_valid_iban(iban)
end

## REXPs
bban_GB_re = r"[A-Z]{4}[0-9]{14}"
@test ismatch(bban_GB_re, "WEST12345698765432")

bban_gb_rexp = IBAN.construct_bban_re(IBAN.bban_format["GB"])
@test bban_GB_re.pattern == bban_gb_rexp.pattern

@test IBAN.bban_length(IBAN.bban_format["GB"]) == 18
@test IBAN.iban_length(IBAN.bban_format["GB"]) == 22

## split IBAN
country, checksum, bban = split_iban(test_iban1)
@test country == "GB"
@test checksum == "82"

## wrong global re
test_wban4 = "GEXX WEST 1234 5698 7654 32"
@test_throws ErrorException split_iban(test_wban4)

## pos2number
@test IBAN.pos2number('A') == 10
@test IBAN.pos2number('B') == 11
@test IBAN.pos2number('Z') == 35
@test IBAN.pos2number('1') == 1
@test IBAN.pos2number('9') == 9
@test_throws ErrorException IBAN.pos2number('!')


text:
function run_scopf_contigency_cuts(ini_file::String, model_type::Type, optimizer; scenario_id::String="", kwargs...)
    goc_data = parse_goc_files(ini_file, scenario_id=scenario_id)
    network = build_pm_model(goc_data)
    return run_scopf_contigency_cuts(network, model_type, optimizer; kwargs...)
end

"""
Solves a SCOPF problem by iteratively checking for violated contingencies and
resolving until a fixed-point is reached
"""
function run_scopf_contigency_cuts(network::Dict{String,<:Any}, model_type::Type, optimizer; max_iter::Int=100, time_limit::Float64=Inf)
    if _IM.ismultinetwork(network)
        error(_LOGGER, "run_scopf_contigency_cuts can only be used on single networks")
    end

    time_start = time()

    network_base = deepcopy(network)
    network_active = deepcopy(network)

    gen_contingencies = network_base["gen_contingencies"]
    branch_contingencies = network_base["branch_contingencies"]

    network_active["gen_contingencies"] = []
    network_active["branch_contingencies"] = []

    multinetwork = build_scopf_multinetwork(network_active)

    result = run_scopf(multinetwork, model_type, optimizer)
    if !(result["termination_status"] == _PM.OPTIMAL || result["termination_status"] == _PM.LOCALLY_SOLVED || result["termination_status"] == _PM.ALMOST_LOCALLY_SOLVED)
        error(_LOGGER, "base-case SCOPF solve failed in run_scopf_contigency_cuts, status $(result["termination_status"])")
    end
    #_PM.print_summary(result["solution"])
    solution = result["solution"]["nw"]["0"]
    solution["per_unit"] = result["solution"]["per_unit"]

    _PM.update_data!(network_base, solution)
    _PM.update_data!(network_active, solution)

    result["iterations"] = 0

    iteration = 1
    contingencies_found = 1
    while contingencies_found > 0
        time_start_iteration = time()

        contingencies = check_contingency_violations(network_base, contingency_limit=iteration)
        #println(contingencies)

        contingencies_found = 0
        #append!(network_active["gen_contingencies"], contingencies.gen_contingencies)
        for cont in contingencies.gen_contingencies
            if cont in network_active["gen_contingencies"]
                warn(_LOGGER, "generator contingency $(cont.label) is active but not secure")
            else
                push!(network_active["gen_contingencies"], cont)
                contingencies_found += 1
            end
        end

        #append!(network_active["branch_contingencies"], contingencies.branch_contingencies)
        for cont in contingencies.branch_contingencies
            if cont in network_active["branch_contingencies"]
                warn(_LOGGER, "branch contingency $(cont.label) is active but not secure")
            else
                push!(network_active["branch_contingencies"], cont)
                contingencies_found += 1
            end
        end

        if contingencies_found <= 0
            info(_LOGGER, "no new violated contingencies found, scopf fixed-point reached")
            break
        else
            info(_LOGGER, "found $(contingencies_found) new contingencies with violations")
        end

        info(_LOGGER, "active contingencies: gen $(length(network_active["gen_contingencies"])), branch $(length(network_active["branch_contingencies"]))")

        time_solve_start = time()
        multinetwork = build_scopf_multinetwork(network_active)
        result = run_scopf(multinetwork, model_type, optimizer)
        if !(result["termination_status"] == _PM.OPTIMAL || result["termination_status"] == _PM.LOCALLY_SOLVED || result["termination_status"] == _PM.ALMOST_LOCALLY_SOLVED)
            warn(_LOGGER, "scopf solve failed with status $(result["termination_status"]), terminating fixed-point early")
            break
        end
        # for (nw,nw_sol) in result["solution"]["nw"]
        #     if nw != "0"
        #         println(nw, " ", nw_sol["delta"])
        #     end
        # end
        info(_LOGGER, "objective: $(result["objective"])")
        solution = result["solution"]["nw"]["0"]
        solution["per_unit"] = result["solution"]["per_unit"]

        _PM.update_data!(network_base, solution)
        _PM.update_data!(network_active, solution)

        time_iteration = time() - time_start_iteration
        time_remaining = time_limit - (time() - time_start)
        if time_remaining < time_iteration
            warn(_LOGGER, "insufficent time for next iteration, time remaining $(time_remaining), estimated iteration time $(time_iteration)")
            break
        end
        iteration += 1
    end

    result["solution"] = solution
    result["iterations"] = iteration

    return result
end

function run_scopf_ptdf_cuts(ini_file::String, model_type::Type, optimizer; scenario_id::String="", kwargs...)
    goc_data = parse_goc_files(ini_file, scenario_id=scenario_id)
    network = build_pm_model(goc_data)
    return run_scopf_ptdf_cuts!(network, model_type, optimizer; kwargs...)
end

"""
Solves a SCOPF problem by iteratively checking for violated branch flow
constraints in contingencies and resolving until a fixed-point is reached.

The base-case model is formulation agnostic.  The flow cuts are based on PTDF
and utilize the DC Power Flow assumption.
"""
function run_scopf_ptdf_cuts!(network::Dict{String,<:Any}, model_type::Type, optimizer; max_iter::Int=100, time_limit::Float64=Inf)
    if _IM.ismultinetwork(network)
        error(_LOGGER, "run_scopf_ptdf_cuts can only be used on single networks")
    end

    time_start = time()

    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = _PM.run_opf(network, model_type, optimizer)
    if !(result["termination_status"] == _PM.OPTIMAL || result["termination_status"] == _PM.LOCALLY_SOLVED || result["termination_status"] == _PM.ALMOST_LOCALLY_SOLVED)
        error(_LOGGER, "base-case OPF solve failed in run_scopf_ptdf_cuts, status $(result["termination_status"])")
    end
    info(_LOGGER, "objective: $(result["objective"])")
    _PM.update_data!(network, result["solution"])

    result["iterations"] = 0

    iteration = 1
    cuts_found = 1
    while cuts_found > 0
        time_start_iteration = time()

        cuts = check_contingencies_branch_power(network, total_cut_limit=iteration, gen_flow_cuts=[], branch_flow_cuts=[])

        cuts_found = length(cuts.gen_cuts)+length(cuts.branch_cuts)
        if cuts_found <= 0
            info(_LOGGER, "no violated cuts found scopf fixed-point reached")
            break
        else
            info(_LOGGER, "found $(cuts_found) branch flow violations")
        end

        append!(network["gen_flow_cuts"], cuts.gen_cuts)
        append!(network["branch_flow_cuts"], cuts.branch_cuts)

        info(_LOGGER, "active cuts: gen $(length(network["gen_flow_cuts"])), branch $(length(network["branch_flow_cuts"]))")

        time_solve_start = time()
        result = run_scopf_cuts(network, model_type, optimizer)
        if !(result["termination_status"] == _PM.OPTIMAL || result["termination_status"] == _PM.LOCALLY_SOLVED || result["termination_status"] == _PM.ALMOST_LOCALLY_SOLVED)
            warn(_LOGGER, "scopf solve failed with status $(result["termination_status"]), terminating fixed-point early")
            break
        end
        info(_LOGGER, "objective: $(result["objective"])")
        _PM.update_data!(network, result["solution"])

        time_iteration = time() - time_start_iteration
        time_remaining = time_limit - (time() - time_start)
        if time_remaining < time_iteration
            warn(_LOGGER, "insufficent time for next iteration, time remaining $(time_remaining), estimated iteration time $(time_iteration)")
            break
        end
        iteration += 1
    end

    result["iterations"] = iteration
    return result
end


text:
module RandomMatrices

using Combinatorics
using GSL
using SpecialFunctions, FastGaussQuadrature
using LinearAlgebra

import Base: isinf, rand, convert
import Distributions: ContinuousUnivariateDistribution,
    ContinuousMatrixDistribution,
    Beta, Chi,
    cdf, pdf, entropy, insupport, mean, median, modes, kurtosis, skewness, std, var, moment

export bidrand,    #Generate random bidiagonal matrix
       tridrand,   #Generate random tridiagonal matrix
       sprand,     #Generate random sparse matrix
       eigvalrand, #Generate random set of eigenvalues
       eigvaljpdf, #Eigenvalue joint probability density
       cumulant, freecumulant, cf, mgf,
       cdf, pdf, entropy, insupport, mean, median, modes, kurtosis, skewness, std, var, moment

const Dim2 = Tuple{Int, Int} #Dimensions of a rectangular matrix

#Classical Gaussian matrix ensembles
include("GaussianEnsembles.jl")

# Classical univariate distributions
include("densities/Semicircle.jl")
include("densities/TracyWidom.jl")

# Ginibre
include("Ginibre.jl")

# determinantal point processes
include("dpp.jl")

#Generating matrices of Haar measure
include("Haar.jl")
include("HaarMeasure.jl")

#Fast histogrammer for matrix eigenvalues - hist_eig
include("FastHistogram.jl")

#Formal power series
include("FormalPowerSeries.jl")

#Statistical tests based on random matrix theory
include("StatisticalTests.jl")

#Stochastic processes
include("StochasticProcess.jl")

end


text:
module TestDispatch

using Test, JETTest, MacroTools
import JETTest: RuntimeDispatchReport, module_filter
import JET: get_reports

# DispatchAnalyzer
# ================

# simple cases
# ------------

# `f(::Int)` a concrete call and just type stable and anything shouldn't be reported
f(a) = a
f(a::Number) = a
let
    analyzer, = @analyze_dispatch f(10)  # should be ok
    @test isempty(get_reports(analyzer))
end

# if the argument type isn't well typed, compiler can't determine which method to call,
# and it will lead to runtime dispatch
let
    analyzer, = analyze_dispatch((Vector{Any},)) do ary
        f(ary[1]) # runtime dispatch !
    end
    @test length(get_reports(analyzer)) == 1
    r = first(get_reports(analyzer))
    @test isa(r, RuntimeDispatchReport)
end

# if we annotate `@noinline` to a function, then its call won't be inlined and will be
# dispatched runtime
@inline   g1(a) = return a
@noinline g2(a) = return a
let
    analyzer, = analyze_dispatch((Vector{Any},)) do ary
        a = ary[1]
        g1(a) # this call should be statically resolved and inlined
        g2(a) # this call should be statically resolved but not inlined, and will be dispatched
    end

    # NOTE the following test is line-sensitive !
    @test length(get_reports(analyzer)) == 1
    r = first(get_reports(analyzer))
    @test isa(r, RuntimeDispatchReport)
    @test any(r.vst) do vf
        vf.file === Symbol(@__FILE__) &&
        vf.line == (@__LINE__) - 9
    end
end

# real-world targets
# ------------------

# the `unoptimize_throw_blocks` configuration disables optimizations on "throw blocks" by default,
# but `DispatchAnalyzer` ignores problems from them, so we don't get error reports here
let
    analyzer, = @analyze_dispatch sin(10)
    @test isempty(get_reports(analyzer))
end

# filter
# ------

# problem: when ∑1/n exceeds 30 ?
function compute(x)
    r = 1
    s = 0.0
    n = 1
    @time while r < x
        s += 1/n
        if s ≥ r
            println("round $r/$x has been finished") # we're not interested type-instabilities within this call
            r += 1
        end
        n += 1
    end
    return n, s
end

# we will get bunch of reports from the `println` call
let
    analyzer, = @analyze_dispatch compute(30)
    @test !isempty(get_reports(analyzer))
end
# if we use different `frame_filter`, a fresh analysis should run
let
    analyzer, = @analyze_dispatch frame_filter=module_filter(@__MODULE__) compute(30)
    @test isempty(get_reports(analyzer))
end

# skip_nonconcrete_calls
# ----------------------

const F1_DEFINITION_LINE = @__LINE__
f1(a) = sin(a)
f1(a::Number) = cos(a)
let
    # by default, we'd ignore error reports from `f1` calls
    analyzer, = analyze_dispatch((Vector{Any},)) do ary
        f1(ary[1]) # runtime dispatch !
    end
    @test length(get_reports(analyzer)) == 1
    r = first(get_reports(analyzer))
    @test isa(r, RuntimeDispatchReport)
end
let
    # when the `skip_nonconcrete_calls` configuration is turned off, we will get error reports
    # from those non-concrete calls of `f1`
    analyzer, = analyze_dispatch((Vector{Any},); skip_nonconcrete_calls=false) do ary
        f1(ary[1]) # runtime dispatch !
    end
    @test length(get_reports(analyzer)) ≥ 3
    @test any(get_reports(analyzer)) do r
        isa(r, RuntimeDispatchReport) &&
        last(r.vst).file === Symbol(@__FILE__) && last(r.vst).line == (@__LINE__) - 5 # report for `f1(ary[1])`
    end
    @test any(get_reports(analyzer)) do r
        isa(r, RuntimeDispatchReport) &&
        last(r.vst).file === Symbol(@__FILE__) && last(r.vst).line == F1_DEFINITION_LINE+1 # report for `sin(a)` within `f1(a)`
    end
    @test any(get_reports(analyzer)) do r
        isa(r, RuntimeDispatchReport) &&
        last(r.vst).file === Symbol(@__FILE__) && last(r.vst).line == F1_DEFINITION_LINE+2 # report for `cos(a)` within `f1(a::Number)`
    end
end

# Test integration
# ================

# runs `f()` in an isolated testset, so that test functions don't influence the current test suite
function with_isolated_testset(f)
    ts = Test.DefaultTestSet("isolated")
    Test.push_testset(ts)
    try
        mktemp() do path, io
            redirect_stdout(io) do
                f()
            end
        end
    catch err
        rethrow(err)
    finally
        Test.pop_testset()
    end
    return ts
end

# do a test that is equivalent to `macroex` but with `test_nodispatch`
macro addfntest(macroex)
    funcex = MacroTools.postwalk(macroex) do x
        pat = :(@test_nodispatch(args__))
        if @capture(x, $pat)
            if isexpr(x, :macrocall)
                kwargs = map(x->Expr(:kw, x.args...), x.args[findall(x->isexpr(x,:(=)), x.args)])
                testex = x.args[findfirst(x->isexpr(x,:call), x.args)::Int]
                return :(test_nodispatch(Base.typesof($(testex.args...)); $(kwargs...)))
            end
        end
        return x
    end

    return quote
        $macroex
        $funcex
    end
end

# positive case
@addfntest let
    ts = with_isolated_testset() do
        @test_nodispatch f(10) # ok
    end
    @test ts.n_passed == 1
end

# negative case
# NOTE the following test is line-sensitive !
g = 10
ff(a) = f(g) # untyped, thus runtime dispatch
let
    ts = with_isolated_testset() do
        @test_nodispatch ff(10)
    end
    @test ts.n_passed == 0
    @test length(ts.results) == 1
    r = ts.results[1]
    @test isa(r, Test.Fail)
    @test r.source === LineNumberNode((@__LINE__)-6, @__FILE__)
end

# actual error within `@test_nodispatch`
# NOTE the following test is line-sensitive !
let
    ts = with_isolated_testset() do
        @test_nodispatch undefvar(10)
    end
    @test ts.n_passed == 0
    @test length(ts.results) == 1
    r = ts.results[1]
    @test isa(r, Test.Error)
    @test r.source === LineNumberNode((@__LINE__)-6, @__FILE__)
    @test occursin("UndefVarError", r.value)
end

# supports `skip` keyword
@addfntest let
    ts = with_isolated_testset() do
        @test_nodispatch skip=true ff(10)
    end
    @test ts.n_passed == 0
    @test length(ts.results) == 1
    r = ts.results[1]
    @test isa(r, Test.Broken)
    @test r.test_type === :skipped

    ts = with_isolated_testset() do
        skip = true
        @test_nodispatch skip=skip ff(10)
    end
    @test ts.n_passed == 0
    @test length(ts.results) == 1
    r = ts.results[1]
    @test isa(r, Test.Broken)
    @test r.test_type === :skipped
end

# supports `broken` keyword
@addfntest let
    ts = with_isolated_testset() do
        @test_nodispatch broken=true ff(10)
    end
    @test ts.n_passed == 0
    @test length(ts.results) == 1
    r = ts.results[1]
    @test isa(r, Test.Broken)
    @test r.test_type === :test_nodispatch

    ts = with_isolated_testset() do
        broken = true
        @test_nodispatch broken=broken ff(10)
    end
    @test ts.n_passed == 0
    @test length(ts.results) == 1
    r = ts.results[1]
    @test isa(r, Test.Broken)
    @test r.test_type === :test_nodispatch

    ts = with_isolated_testset() do
        @test_nodispatch broken=true f(10)
    end
    @test ts.n_passed == 0
    @test length(ts.results) == 1
    r = ts.results[1]
    @test isa(r, Test.Error)
    @test r.test_type === :test_unbroken
end

# supports JET's configurations
@addfntest let
    ts = with_isolated_testset() do
        @test_nodispatch sin(10) # ok
    end
    @test ts.n_passed == 1

    ts = with_isolated_testset() do
        @test_nodispatch sin(10) # still be ok
    end
    @test ts.n_passed == 1

    ts = with_isolated_testset() do
        @test_nodispatch skip_unoptimized_throw_blocks=false sin(10)
    end
    @test ts.n_passed == 0
    @test length(ts.results) == 1
    r = ts.results[1]
    @test isa(r, Test.Fail)

    ts = with_isolated_testset() do
        @test_nodispatch skip_unoptimized_throw_blocks=false broken=true sin(10)
    end
    @test ts.n_passed == 0
    @test length(ts.results) == 1
    r = ts.results[1]
    @test isa(r, Test.Broken)
end

end # module TestDispatch


text:
"""
    OrdinalMarginLoss <: SupervisedLoss

Modifies a margin loss `loss` to be used on an ordinal domain with
number of levels `N`. It treats each level as an integer between
1 and `N`, inclusive, and penalizes output according to the sum of
level thresholds crossed relative to target

Assumes target is encoded in an Index encoding scheme where levels are
numbered between 1 and `N`
"""
struct OrdinalMarginLoss{L<:MarginLoss} <: SupervisedLoss
    loss::L
    N::Int
    OrdinalMarginLoss{L}(loss::L, N::Int) where {L<:MarginLoss} = new{L}(loss, N)
end

function OrdinalMarginLoss(loss::L, N::Int) where {L<:MarginLoss}
    OrdinalMarginLoss{L}(loss, N)
end

@generated function (::Type{T})(N::Int, args...) where {T<:OrdinalMarginLoss}
    L = typeof(T) == UnionAll ? T.var.ub : T.parameters[1]
    :(OrdinalMarginLoss($L(args...), N))
end

for fun in (:value, :deriv, :deriv2)
    @eval @fastmath function ($fun)(
            l::OrdinalMarginLoss,
            target::Number,
            output::Number)
        not_target = Int(1 != target)
        sgn = sign(target - 1)
        retval = not_target * ($fun)(l.loss, sgn, output - 1)
        for t = 2:l.N
            not_target = Int(t != target)
            sgn = sign(target - t)
            retval += not_target * ($fun)(l.loss, sgn, output - t)
        end
        retval
    end
end

for prop in [:isminimizable, :isdifferentiable,
             :istwicedifferentiable,
             :isconvex, :isstrictlyconvex,
             :isstronglyconvex, :isnemitski,
             :isunivfishercons, :isfishercons,
             :islipschitzcont, :islocallylipschitzcont]
    @eval ($prop)(l::OrdinalMarginLoss) = ($prop)(l.loss)
end

for fun in (:isdifferentiable, :istwicedifferentiable)
    @eval function ($fun)(
            l::OrdinalMarginLoss,
            target::Number,
            output::Number)
        for t = 1:target - 1
            ($fun)(l.loss, output - t) || return false
        end
        for t = target + 1:l.N
            ($fun)(l.loss, t - output) || return false
        end
        return true
    end
end


text:
t₀ = 0.
q₀ = [1.]
p₀ = [1.]
x₀ = [1., 1.]
λ₀ = [0.]

function f_ode(x, fx)
    fx[1] = x[1]
end

ode  = ODE{Float64, typeof(f_ode)}(1, 1, f_ode, t₀, q₀)
ode1 = ODE(f_ode, t₀, q₀)
ode2 = ODE(f_ode, q₀)

@test ode == ode1
@test ode == ode2

function f_pode(q, p, fp)
    fq[1] = q[1]
end

function g_pode(q, p, fp)
    fp[1] = 2p[1]
end

pode  = PODE{Float64, typeof(f_pode), typeof(g_pode)}(1, 1, f_pode, g_pode, t₀, q₀, p₀)
pode1 = PODE(f_pode, g_pode, t₀, q₀, p₀)
pode2 = PODE(f_pode, g_pode, q₀, p₀)

@test pode == pode1
@test pode == pode2

sode  = SODE{Float64}(1, 1, q₀, p₀, t₀, do_nothing, do_nothing, f_pode, g_pode)
sode1 = SODE(q₀, p₀, t₀; v=f_pode, f=g_pode)
sode2 = SODE(q₀, p₀; v=f_pode, f=g_pode)

@test sode == sode1
@test sode == sode2
# @test pode1 == pode2

function f_dae(x, fx)
    fx[1] = x[1]
    fx[2] = x[2]
end

function u_dae(x, λ, fu)
    fu[1] = +λ[1]
    fu[1] = -λ[1]
end

function ϕ_dae(x, λ, fϕ)
    fϕ[1] = x[2] - x[1]
end

dae  = DAE{Float64}(2, 1, f_dae, u_dae, ϕ_dae, t₀, x₀, λ₀)
dae1 = DAE(f_dae, u_dae, ϕ_dae, t₀, x₀, λ₀)
dae2 = DAE(f_dae, u_dae, ϕ_dae, x₀, λ₀)

@test dae == dae1
@test dae == dae2

function f_pdae(q, p, fp)
    fq[1] = q[1]
end

function g_pdae(q, p, fp)
    fp[1] = p[1]
end

function u_pdae(q, p, λ, fu)
    fu[1] = +λ[1]
end

function v_pdae(q, p, λ, fv)
    fv[1] = -λ[1]
end

function ϕ_pdae(q, p, λ, fϕ)
    fϕ[1] = p[1] - q[1]
end

pdae  = PDAE{Float64}(1, 1, f_pdae, g_pdae, u_pdae, v_pdae, ϕ_pdae, t₀, q₀, p₀, λ₀)
pdae1 = PDAE(f_pdae, g_pdae, u_pdae, v_pdae, ϕ_pdae, t₀, q₀, p₀, λ₀)
pdae2 = PDAE(f_pdae, g_pdae, u_pdae, v_pdae, ϕ_pdae, q₀, p₀, λ₀)

@test pdae == pdae1
@test pdae == pdae2


text:
### Document
abstract type AbstractDocument end

struct LdaDocument <: AbstractDocument
  # this is a fully observed data from the LDA model
  theta::Array{Float64,1} # the topic probs for the doc
  z::Array{Int64,1} # the topic for each word
  terms::Array{Int64,1} # the word tokens

  LdaDocument(alpha::Array{Float64,1},
                  P::Array{Float64,2},
                  N::Int64) = begin # length of the doc
    d = new(
      Array{Float64,1}(undef,size(P,2)),
      Array{Int64,1}(undef,N),
      Array{Int64,1}(undef,N)
    )
    GenerateDoc(d,alpha,P)
    return d
  end

  LdaDocument(theta::Array{Float64,1},
              z::Array{Int64,1},
              terms::Array{Int64,1}) = begin
    d = new(theta,z,N)
    return d
  end
end

function GenerateDoc(doc::LdaDocument,
                     alpha::Array{Float64,1},
                     Phi::Array{Float64,2})
    dd = Dirichlet(alpha)
    doc.theta .= vec(rand(dd,1))
    cat = Categorical(vec(doc.theta))
    doc.z .= rand(cat,length(doc))
    for i in 1:length(doc)
        @inbounds dc = Categorical(Phi[:,doc.z[i]])
        @inbounds doc.terms[i] = rand(dc,1)[1]
    end
    return
end

struct Document <: AbstractDocument
  #this is actual data, where only the terms are observed
  terms::Array{Int64,1} # the word tokens
  Document(terms::Array{Int64,1}) = new(terms)
end

function length(doc::AbstractDocument)
  return size(doc.terms,1)
end

### Corpus
abstract type AbstractCorpus end

struct LdaCorpus <: AbstractCorpus
  # this is a fully observed data from the LDA model
  docs::Array{LdaDocument,1}
  alpha::Array{Float64,1}
  beta::Array{Float64,1}
  Phi::Array{Float64,2}
  weights::Array{Array{Float64,1}} # only unweighted terms supported
  lexicon::Array{String,1}

  LdaCorpus(k::Int64,
            lexLength::Int64,
            corpLambda::Int64,
            corpLength::Int64,
            scaleK::Float64,
            scaleL::Float64) = begin # length of the doc
    w = Array{Array{Float64,1},1}(undef,corpLength)
    lex = string.([1:1:lexLength;]) # there is no
    a = fill(scaleK,k) # scale parameter for the Dirichlet topic prior
    b = fill(scaleL,lexLength) # scale parameter for the Dirichlet token prior
    dl = Poisson(corpLambda)
    docLengths = rand(dl,corpLength) # the lengths of the docs in the corpus
    db = Dirichlet(b)
    P = rand(db,k) # the Dirichlet token prior, containing one lexLength vector for each k
    d = Array{LdaDocument,1}(undef,corpLength)
    for i in 1:corpLength
      w[i] = ones(docLengths[i])
      @inbounds d[i] = LdaDocument(a,P,docLengths[i])
    end
    return new(d, a, b, P, w, lex)
  end

  LdaCorpus(docs::Array{LdaDocument,1}, # the documents
            alpha::Array{Float64,1},
            beta::Array{Float64,1},
            Phi::Array{Float64,2},
            weights::Array{Float64,1},
            lexicon::Array{String,1}) = begin
    c = new(docs,alpha,beta,Phi,weights)
    return c
  end
end

function CorpusTopics(corpus::LdaCorpus)
  cat(dims=2,map(i->vec(i.theta), corpus.docs)...) # get a 2d array of (document wise) mixed membership for the corpus
end

function AllTopics(corpus::LdaCorpus)
  alltopics = []
  for i in 1:length(corpus)
    append!(alltopics,corpus.docs[i].z)
  end
  return convert(Array{Int,1},alltopics)
end

struct Corpus <: AbstractCorpus
  docs::Array{Document,1}
  weights::Array{Array{Float64,1},1}
  lexicon::Array{String,1}

  Corpus(docs::Array{Document,1},
         weights::Array{Array{Float64,1},1},
         lexicon::Array{String,1}) = begin
    return new(
      docs,
      weights,
      lexicon
    )
  end

  Corpus(docs::Array{Document,1},
         lexicon::Array{String,1}) = begin
    return new(
      docs,
      map(x -> ones(Float64,length(x)), docs), # no weights
      lexicon
    )
  end
end

function length(corpus::AbstractCorpus)
  return length(corpus.docs)
end

# Expand  a term:count pair into a <count>-length sequence [term, term, ....]
function termToWordSequence(term::AbstractString)
  parts = split(term, ":")
  fill(parse(Int64, parts[1]) + 1, parse(Int64, parts[2]))
end

function readDocs(stream)
    corpus = readlines(stream)
    docs = Array{Document,1}(undef,length(corpus))
    for i in 1:length(corpus)
      @inbounds terms = split(corpus[i], " ")[2:end]
      @inbounds docs[i] = Document(termToWordSequence(terms[1]))
      for ii in 2:length(terms)
        @inbounds append!(docs[i].terms, termToWordSequence(terms[ii]))
      end
    end
    return docs
end

function readLexicon(stream)
  lines = readlines(stream)
  chomped = map(chomp, convert(Array{AbstractString,1}, lines))
  convert(Array{String,1},chomped) # convert from substrings
end


text: Dict{String, Any}("bestdayever" => Dict{String, Any}("value" => "1987-07-05", "type" => "date-local"))


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie
CairoMakie.activate!() #HIDE
let
    x = -1:0.02:1
    y = -1.5:0.02:2
    egg(x,y) = x^2 + y^2/(1.4 + y/5)^2
    segg = [egg(x,y) for x in x, y in y]
    fig = Figure(resolution = (470, 550))
    ax = Axis(fig, xlabel = "x", ylabel = "y", backgroundcolor = :black,
    xgridstyle=:dash, ygridstyle=:dash, xgridcolor = :grey, ygridcolor = :grey)
    cl =contour!(x, y, segg, linewidth = 0.85,colormap = :viridis,
                levels = 0:0.02:1)
    cbar = Colorbar(fig, cl, label ="egg-l", labelpadding = 0, width = 15,
                ticksize=15, tickalign = 1, height = Relative(1))
    fig[1, 1] = ax
    fig[1, 2] = cbar
    colgap!(fig.layout, 7)
    fig
    save(joinpath(@__DIR__, "output", "contour1.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
#Fichier qui permet de créer les différents PLNE pour le problème du VRP

using JuMP
using CPLEX
#using GLPK

function createVRP_MTZ(params, nodes, demands, costs, t)

    model = Model(CPLEX.Optimizer)
    #model = Model(GLPK.Optimizer)

    #=

    m = nb de véhicule
    Q = capacité d'un véhicule

    =#

    m = params["k"]
    n = params["n"]
    Q = params["Q"]

    ####################VARIABLES#######################
    @variable(model, x[i=0:n, j=0:n], Bin) #Variable binaire x_(i,j) pour chaque arête
    for i in 0:n
        delete(model, x[i, i]) #on enlève les variables qui correspondent aux arêtes en trop (les (i, i))
    end

    @variable(model, 0 <= w[i=1:n] <= Q) #Variables w_i

    ####################CONTRAINTES####################
    @constraint(model, sum(x[0, j] for j in 1:n) <= m)#contrainte 6
    @constraint(model, sum(x[i, 0] for i in 1:n) <= m)#contrainte 7

    for i in 1:n

        nodesIndexWithoutI = filter(e -> e != i, 0:n)

        @constraint(model, sum(x[i,j] for j in nodesIndexWithoutI) == 1)#contraintes 8
        @constraint(model, sum(x[j,i] for j in nodesIndexWithoutI) == 1)#contraintes 9

    end

    for i in 1:n

        nodesIndexWithoutI = filter(e -> e != i, 1:n)
        for j in nodesIndexWithoutI
            @constraint(model, w[i] - w[j] >= demands[i, t] - (Q+demands[i, t]) * (1 - x[i, j]))#contrainte 10
        end
    end

    ####################FONCTION OBJECTIF####################
    @objective(model, Min, sum(edgeCost * x[edge[1], edge[2]] for (edge, edgeCost) in costs))

    return model

end


text:
"""
    AbstractLayout

Supertype for specifying the layout of new instances of arrays.
"""
abstract type AbstractLayout{T} end

""" MemoryLayout """
abstract type MemoryLayout{T} <: AbstractLayout{T} end

struct DynamicMemory{T} <: MemoryLayout{T}
    length::Int
end

struct FixedMemory{T} <: MemoryLayout{T}
    length::Int
end

struct StaticMemory{T,N} <: MemoryLayout{T} end

MemoryLayout{T}(length::Int) where {T} = FixedMemory{T}(length)
MemoryLayout{T}(::StaticInt{N}) where {T,N} = StaticMemory{T,N}()

""" AbstractArrayLayout """
abstract type AbstractArrayLayout{T,N,L} <: AbstractLayout{T} end

struct StrideLayout{T,N,L<:MemoryLayout{T},B,C,S,I<:Tuple{Vararg{Any,N}}} <: AbstractArrayLayout{T,N,L}
    layout::L
    contiguous_batch_size::B
    contiguous_axis::C
    strides::S
    indices::I

    function StrideLayout(x::AbstractArray{T,N}) where {T,N}
        inds = map(ArrayInterface.indices, axes(x))
        m = MemoryLayout{T}(prod(map(static_length, inds)))
        b = ArrayInterface.contiguous_batch_size(x)
        c = ArrayInterface.contiguous_axis(x)
        s = ArrayInterface.strides(x)
        return new{T,N,typeof(m),typeof(b),typeof(c),typeof(s),typeof(inds)}(
            m, b, c, s, inds
        )
    end
end

struct ArrayLayout{T,N,L<:AbstractLayout{T},I<:Tuple{Vararg{Any,N}}} <: AbstractArrayLayout{T,N,L}
    layout::L
    indices::I
end

abstract type SparseLayout{T,N,L} <: AbstractArrayLayout{T,N,L} end

struct DiagonalLayout{T,L<:ArrayLayout{T,1}} <: SparseLayout{T,2,L}
    layout::L
end

""" layout """
function layout(x::AbstractArray{T,N}) where {T,N}
    if ArrayInterface.defines_strides(x)
        return StrideLayout(x)
    else
        return ArrayLayout(x)
    end
    inds = map(ArrayInterface.indices, axes(x))
    return ArrayLayout(MemoryLayout{T}(prod(map(static_length, inds))), inds)
end
function layout(x::AbstractVector{T}) where {T}
    if ArrayInterface.can_change_size(x)
        return DynamicMemory(length(x))
    else
        return MemoryLayout{T}(length(x))
    end
end
layout(x::Diagonal) = DiagonalLayout(layout(x.diag))

""" layout """
layout(x::DynamicMemory, ::Type{T}) where {T} = DynamicMemory{T}(x.length)
layout(x::FixedMemory, ::Type{T}) where {T} = FixedMemory{T}(x.length)
layout(x::StaticMemory{T1,N}, ::Type{T2}) where {T1,T2,N} = StaticMemory{T2,N}()
layout(x::ArrayLayout, ::Type{T}) where {T} = ArrayLayout(layout(x.layout, T), x.axes)
layout(x::DiagonalLayout, ::Type{T}) where {T} = DiagonalLayout(layout(x.layout, T))
layout(x::AbstractArray, ::Type{T}) where {T} = layout(layout(x), T)

""" combine_layouts """
@inline function combine_layouts(bc::AbstractBroadcasted)
    typed_layout = Fix2(layout, Base.Broadcast.combine_eltypes(bc.f, bc.args))
    return broadcast_layout(map(typed_layout, bc.args)...)
end

# This is the where new subtypes of `AbstractLayout` need to specify how they combine with
# eachother
""" broadcast_layout(x::AbstractLayout, y::AbstractLayout) """
broadcast_layout(x) = x
broadcast_layout(x, y, zs...) = broadcast_layout(broadcast_layout(x, y), zs...)
function broadcast_layout(x::DiagonalLayout, y::DiagonalLayout)
    return DiagonalLayout(broadcast_layout(x.layout, y.layout))
end
function broadcast_layout(x::AbstractArrayLayout, y::AbstractArrayLayout)
    return ArrayLayout(
        combine_layout(memory_layout(x), memory_layout(y)),
        _bcs(axes(x), axes(y))
    )
end

_bcs(::Tuple{}, ::Tuple{}) = ()
_bcs(::Tuple{}, y::Tuple) = (first(y), _bcs((), tail(y))...)
_bcs(shape::Tuple, ::Tuple{}) = (shape[1], _bcs(tail(shape), ())...)
function _bcs(x::Tuple, y::Tuple)
    return (ArrayInterface.broadcast_axis(first(x), first(y)), _bcs(tail(x), tail(y))...)
end

# TODO
""" vcat_layouts """

""" hcat_layouts """

""" cat_layouts """


text:
include("dimDiff.jl")
include("PV.jl")
include("vertical_flux.jl")

include("layer.jl")
include("smoother.jl")