text:
export  computeSP,
        computeSP2

function computeSP(a_::Vector{Int64},
                   α::Vector{Vector{Float64}},β::Vector{Vector{Float64}},
                   nodes::Vector{Vector{Float64}},weights::Vector{Vector{Float64}},
                   ind::Matrix{Int64};
                   issymmetric::BitArray=falses(length(α)),
                   zerotol::Float64=1e-10)
                  # display("war hier")
    @assert minimum(a_)>=0 "no negative degrees allowed"
    l, p = size(ind); # p-variate basis
    p == 1 && computeSP(a_,α[1],β[1],nodes[1],weights[1];issymmetric=issymmetric[1])
    l -= 1
    # m = length(a_) # m-dimensional scalar product
    @assert maximum(a_)<=l "not enough elements in multi-index (requested: $(maximum(a)), max: $l)"
    @assert length(α)==length(β)==length(nodes)==length(weights)==length(issymmetric)==p "inconsistent number of rec. coefficients and/or nodes/weights"

    a = Vector{Int64}()
    for aa in a_
        !iszero(aa) && push!(a,aa)
    end

    if  iszero(length(a))
        return prod(β[i][1] for i=1:p)
    elseif isone(length(a))
        return 0.
    else
        inds_uni = multi2uni(a,ind)
        val = 1.
        for i=1:p
            # compute univariate scalar product
            @inbounds v = computeSP(inds_uni[i,:],α[i],β[i],nodes[i],weights[i];issymmetric=issymmetric[i])
            isapprox(v,0,atol=zerotol) ? (return 0.) : (val*=v)
        end
    end
    return val
end

function computeSP(a::Vector{Int64},op::Vector{OrthoPoly},q::Vector{Quad},ind::Matrix{Int64})
    @assert length(op)==length(q) "Inconsistent lengh of polynomials and quadrature rules"
    @assert issymmetric.(op)==issymmetric.(q)
    α,β = coeffs(op)
    n,w = nw(q)
    computeSP(a,α,β,n,w,ind;issymmetric=issymmetric.(op))
end

function computeSP(a::Vector{Int64},mOP::MultiOrthoPoly)
    @assert typeof(mOP.uni)==Vector{OrthoPolyQ} "no quadrature rules provided"
    ops = [ mOP.uni[i].op for i=1:length(mOP.uni) ]
    quads = [ mOP.uni[i].quad for i=1:length(mOP.uni) ]
    computeSP(a,ops,quads,mOP.ind)
end

"""
__Univariate__
```
computeSP(a::Vector{Int64},α::Vector{Float64},β::Vector{Float64},nodes::Vector{Float64},weights::Vector{Float64};issymmetric::Bool=false)
computeSP(a::Vector{Int64},op::OrthoPoly,q::Quad;issymmetric=issymmetric(op))
computeSP(a::Vector{Int64},opq::OrthoPolyQ)
```
__Multivariate__
```
computeSP( a::Vector{Int64},
           α::Vector{Vector{Float64}},β::Vector{Vector{Float64}},
           nodes::Vector{Vector{Float64}},weights::Vector{Vector{Float64}},
           ind::Matrix{Int64};
           issymmetric::BitArray=falses(length(α)))
computeSP(a::Vector{Int64},op::Vector{OrthoPoly},q::Vector{Quad},ind::Matrix{Int64})
computeSP(a::Vector{Int64},mOP::MultiOrthoPoly)
```

Computes the scalar product
```math
\\langle \\phi_{a_1},\\phi_{a_2},\\cdots,\\phi_{a_n} \\rangle,
```
where `n = length(a)`.
This requires to provide the recurrence coefficients `(α,β)` and the quadrature rule
`(nodes,weights)`, as well as the multi-index `ind`.
If provided via the keyword `issymmetric`, symmetry of the weight function is exploited.
All computations of the multivariate scalar products resort back to efficient computations
of the univariate scalar products.
Mathematically, this follows from Fubini's theorem.

The function is multiply dispatched to facilitate its use with `OrthoPolyQ` or a
suitable combination of `OrthoPoly` and its quadrature rule `Quad`.

!!! note
    - Zero entries of ``a`` are removed automatically to simplify computations, which follows from
    ```math
    \\langle \\phi_i, \\phi_j, \\phi_0,\\cdots,\\phi_0 \\rangle = \\langle \\phi_i, \\phi_j \\rangle,
    ```
    because `\\phi_0 = 1`.
    - It is checked whether enough quadrature points are supplied to solve the integral exactly.
"""
function computeSP(a_::Vector{Int64},α::Vector{Float64},β::Vector{Float64},nodes::Vector{Float64},weights::Vector{Float64};issymmetric::Bool=false)
    @assert minimum(a_) >= 0 "no negative degrees allowed"
    # remove zero entries
    # this works because ``<Φ_i,Φ_j,Φ_0,...,Φ_0> = <Φ_i,Φ_j>``
    # hence, avoiding unnecessary operations that introduced numerical gibberish
    a = Vector{Int64}()
    for aa in a_
        !iszero(aa) && push!(a,aa)
    end
    # a = filter(!iszero,a_)
    # simplification in case the measure is symmetric w.r.t t=0
    # exploit symmetry of the density, see Theorem 1.17 in W. Gautschi's book
    # and exploit the fact that
    (issymmetric && isodd(sum(a))) && return 0.
    # Gauss quadrature rules have exactness 2N-1
    @assert Int(ceil(0.5*(sum(a)+1)))<=length(nodes) "not enough nodes to integrate exactly ($(length(nodes)) provided, where $(Int(ceil(0.5*(sum(a)+1)))) are needed)."
    if iszero(length(a))
        return first(β)
    elseif isone(length(a))
        return 0.
    elseif length(a) == 2
        @inbounds a[1] == a[2] ? (return computeSP2(first(a),β)[end]) : (return 0.)
    else
        # display("war hier")
        f = ones(Float64,length(nodes))
        @simd for i = 1:length(a)
            @inbounds f = f.*evaluate(a[i],nodes,α,β)
        end
        # f = reduce(*,hcat(map(i->evaluate(a[i],nodes,α,β),1:length(a))...),dims=2)
        # display(f-f_)
        return dot(weights,f)
    end
end
function computeSP(a::Vector{Int64},op::OrthoPoly,q::Quad)
    @assert issymmetric(op)==issymmetric(q) "inconsistent symmetries"
    computeSP(a,op.α,op.β,q.nodes,q.weights;issymmetric=issymmetric(op))
end
computeSP(a::Vector{Int64},opq::OrthoPolyQ) = computeSP(a,opq.op,opq.quad)

"""
    computeSP2(n::Int64,β::Vector{Float64})
    computeSP2(n::Int64,op::OrthoPoly) = computeSP2(n,op.β)
    computeSP2(op::OrthoPoly) = computeSP2(op.deg,op.β)
    computeSP2(opq::OrthoPolyQ) = computeSP2(opq.op)

Computes the `n` *regular* scalar products aka 2-norms of the orthogonal polynomials, namely
```math
\\|ϕ_i\\|^2 = \\langle \\phi_i,\\phi_i\\rangle \\quad \\forall i \\in \\{ 0,\\dots,n \\}.
```
Notice that only the values of `β` of the recurrence coefficients `(α,β)` are required.
The computation is based on equation (1.3.7) from Gautschi, W. "Orthogonal Polynomials: Computation and Approximation".
Whenever there exists an analytic expressions for `β`, this function should be used.

The function is multiply dispatched to facilitate its use with the composite types
`OrthoPoly` and `OrthoPolyQ`.
"""
function computeSP2(n::Int64,β::Vector{Float64})
    @assert n>=0 "can only compute scalar products for non-negative degrees"
    @assert length(β) >= n + 1
    n == 0 && return β[1]
    s = ones(Float64,n+1)
    @inbounds s[1] = β[1]  # order 0computeSP2
    for i=2:n+1
        @inbounds s[i] = s[i-1]*β[i]
    end
    return s
end
computeSP2(n::Int64,op::OrthoPoly) = computeSP2(n,op.β)
computeSP2(op::OrthoPoly) = computeSP2(op.deg,op.β)
computeSP2(opq::OrthoPolyQ) = computeSP2(opq.op)

# function computeSP2(n::Int64,α::Vector{Float64},β::Vector{Float64},nodes::Vector{Float64},weights::Vector{Float64})
#     @assert n>=0 "degree has to be non-negative"
#     @assert n<=length(α) "not enough recurrence coefficients"
#     Nquad, Nrec = length(nodes), length(α)
#     n==0 ? (return β[1]) : (dot(weights,evaluate(n,nodes,α,β).^2))
# end
# computeSP2(n::Int64,opq::OrthoPolyQ) = computeSP2(n,opq.op.α,opq.op.β,opq.quad.nodes,opq.quad.weights)


text:
abstract type DataStorage end

"""
    savedataset(::DataStorage, blocks, data)

Save a data container `data` with sample blocks `blocks` to a storage.
along with metadata on the dataset.
"""
function savedataset(storage::DataStorage, block, data)
    savemetadata(storage, getmetadata(block, data))
    savedataset!(gethandle(storage), storage, block, data)
end

getmetadata(block, data) = Dict("block" => block, "n" => nobs(data))

"""
    savedataset!(handle, ::DataStorage, block, data)

Save a `block` of `data` to a storage's  `handle`.
"""
function savedataset! end

savedataset!(handle, storage, block::FastAI.WrapperBlock, data) =
    savedataset!(handle, storage, FastAI.wrapped(block), data)
"""
    savemetadata(::DataStorage, blocks, data)
"""
function savemetadata end
function loadmetadata end

function gethandle end

"""
    loaddataset(::DataStorage) -> (data, blocks)

Load a dataset
"""
function loaddataset(storage::DataStorage)
    metadata = loadmetadata(storage)
    block = metadata["block"]
    data = loaddataset(gethandle(storage), storage, block)
    return data, block
end

loaddataset(handle, storage, block::FastAI.WrapperBlock) =
    loaddataset(handle, storage, FastAI.wrapped(block))

ext: jl


text:
__precompile__()

module OrbitTomographyAddOns

# Import standard packages
using JLD2
using FileIO
using Interpolations
using ProgressMeter
using Distributions

# Import special packages
using Equilibrium
using GuidingCenterOrbits
using OrbitTomography

greet() = print("Hello World!")

include("add_ons.jl")
export add_noise, get_orbel_volume

end # module


text:
Model_NCS = HA_solver.NCS_growth_params()
@testset "Solve_Smolyak" begin
  #test for NCS_growth_params, check if params make sense
  Model_NCS.β = 2.0
  @test_throws AssertionError HA_solver.Solve_Smolyak(Model_NCS)
  Model_NCS.β = 0.984

  Model_NCS.σ = 0.0
  @test_throws AssertionError HA_solver.Solve_Smolyak(Model_NCS)
  Model_NCS.σ = 2.0

  Model_NCS.rho = -1.0
  @test_throws AssertionError HA_solver.Solve_Smolyak(Model_NCS)
  Model_NCS.rho = 0.979

  Model_NCS.sigma_e = 0.0
  @test_throws AssertionError HA_solver.Solve_Smolyak(Model_NCS)
  Model_NCS.sigma_e = 0.0072

  Model_NCS.α = -0.5
  @test_throws AssertionError HA_solver.Solve_Smolyak(Model_NCS)
  Model_NCS.α = 0.33

  Model_NCS.δ = -0.08
  @test_throws AssertionError HA_solver.Solve_Smolyak(Model_NCS)
  Model_NCS.δ = 0.025

  #check if the model solves, if it does not solve (solution not converge), that will throw error
  Ev, grid = Solve_Smolyak(Model_NCS)
  #if no error if thrown than finish tests
  @test 1>0

end

@testset "simulate_smolyak_NGS" begin
  #test for NCS_growth_params, check if params make sense
  Model_NCS.β = 2.0
  @test_throws AssertionError HA_solver.simulate_smolyak_NGS(Model_NCS, Ev)
  Model_NCS.β = 0.984

  Model_NCS.sigma_e = -10.0
  @test_throws AssertionError HA_solver.simulate_smolyak_NGS(Model_NCS, Ev)
  Model_NCS.sigma_e = 0.0072

  Model_NCS.ρ = -1.0
  @test_throws AssertionError HA_solver.simulate_smolyak_NGS(Model_NCS, Ev)
  Model_NCS.ρ = 0.979

  Model_NCS.α = -0.5
  @test_throws AssertionError HA_solver.simulate_smolyak_NGS(Model_NCS, Ev)
  Model_NCS.α = 0.33

  Model_NCS.δ = -0.08
  @test_throws AssertionError HA_solver.simulate_smolyak_NGS(Model_NCS, Ev)
  Model_NCS.δ = 0.025
  #check if the model solves, if it does not solve (solution not converge), that will throw error
  sim_result = Solve_Smolyak(Model_NCS)
  #check if the simulation results have sense
  @test abs(sim_result[1] - Model_NCS.ρ)<=0.01
  @test abs(sim_result[2] - Model_NCS.sigma_e)<=0.001
  #interest rate for this model should be small
  @test sim_result[3]<=0.05
  #if no error if thrown than finish tests
  @test 1>0

end


text:
"""
    za = agg(z, nas, nat)

Aggregate field `z` to an array `za` of size `(nas,nas,nat)`
"""
function agg(zi, nas, nat)
    nss = size(zi)
    if (length(nss)>=3)
        ns = nss[1]; nt = nss[3]
    else
        ns = nss[1]; nt = 1
    end
    sdim = div(ns, nas)
    tdim = div(nt, nat)
    ir = 1:sdim; jr = ir; kr = 1:tdim
    rs = 0:sdim:(ns-sdim)
    rt = 0:tdim:(nt-tdim)
    z = zeros(Float64, nas, nas, nat)
    n = zeros(Int64, nas, nas, nat)
    for i = 1:sdim
        for j = 1:sdim
            for k = 1:tdim
                n = n+.~isnan.(zi[i.+rs,j.+rs,k.+rt])
		z = z+zi[i.+rs,j.+rs,k.+rt].*.~isnan.(zi[i.+rs,j.+rs,k.+rt])
            end
        end
    end
    z = z./n
end


text:
using Vimes, Pkg, Test

@testset "Vimes.jl" begin
    @testset "Run Vimes on JSON" begin
        tmp_depot = mktempdir()
        tmp_env = mktempdir()
        atexit((tmp_depot) -> rm(tmp_depot; force = true, recursive = true))
        atexit((tmp_env) -> rm(tmp_env; force = true, recursive = true))
        original_depot_path = deepcopy(Base.DEPOT_PATH)
        empty!(Base.DEPOT_PATH)
        pushfirst!(Base.DEPOT_PATH, tmp_depot)
        Pkg.activate(tmp_env; shared = false)
        Pkg.develop("JSON")
        json_directory = joinpath(tmp_depot, "dev", "JSON")
        @test Vimes.go(json_directory; tests = 1) isa Integer
        @test Vimes.go(json_directory; tests = 1, dead = true) isa Integer
        @test Vimes.go(json_directory; tests = 2, procs = 2) isa Integer
        empty!(Base.DEPOT_PATH)
        append!(Base.DEPOT_PATH, original_depot_path)
        rm(tmp_depot; force = true, recursive = true)
        rm(tmp_env; force = true, recursive = true)
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE

"""
    abaqus_download(model_name; dryrun=false)

Download ABAQUS model from Internet. `model_name` is the name of the input
file.

Given some model name from documentation, e.g., `et22sfse`, download that
file to local file system. This function uses environment variables to
determine the download url and place of storage.

In order to use this function, one must set environment variable
`ABAQUS_DOWNLOAD_URL`, which determines a location where to download. For
example, if the path to model is `https://domain.com/v6.14/books/eif/et22sfse.inp`,
`ABAQUS_DOWNLOAD_URL` will be the basename of that path, i.e.,
`https://domain.com/v6.14/books/eif`.

By default, the model will be downloaded to current directory. If that is not
desired, one can set another environment variable `ABAQUS_DOWNLOAD_DIR`, and
in that case the file will be downloaded to that directory.

Function call will return full path to downloaded file or nothing, if download
is failing because of missing environment variable `ABAQUS_DOWNLOAD_DIR`.
"""
function abaqus_download(model_name, env=ENV; dryrun=false)
    path = get(env, "ABAQUS_DOWNLOAD_DIR", "")
    fn = joinpath(path, model_name)
    if isfile(fn)  # already downloaded
        return fn
    end
    if !haskey(env, "ABAQUS_DOWNLOAD_URL")
        error("ABAQUS input file $fn not found and `ABAQUS_DOWNLOAD_URL` not ",
              "set, unable to download file. To enable automatic model ",
              "downloading, set url to models to environment variable
              `ABAQUS_DOWNLOAD_URL`")
    end
    url = joinpath(env["ABAQUS_DOWNLOAD_URL"], model_name)
    @debug("Downloading model $model_name to $fn")
    dryrun || download(url, fn)
    return fn
end


text:
for n = 10 : 1
    print(n)
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Documenter, Atmosphere

makedocs(
    sitename = "Atmosphere",
    modules = [Atmosphere],
    pages   = ["Home" => "index.md",
        "Library" => "library.md"]
)

deploydocs(
    repo   = "github.com/byuflowlab/Atmosphere.jl")


text:
using Images, Colors, IndirectArrays, Test

@testset "ColorizedArray" begin
    intensity = [0.1 0.3; 0.2 0.4]
    labels = IndirectArray([1 2; 2 1], [RGB(1,0,0), RGB(0,0,1)])
    A = ColorizedArray(intensity, labels)
    target = intensity .* labels
    @test eltype(A) == RGB{Float64}
    @test size(A) == (2,2)
    @test axes(A) == (Base.OneTo(2), Base.OneTo(2))
    for i = 1:4
        @test A[i] === target[i]
    end
    for j = 1:2, i = 1:2
        @test A[i,j] === target[i,j]
    end
    for (a,t) in zip(A, target)
        @test a === t
    end
end

nothing


text:
module TestMOIWrapper

using MathOptInterface
using NLopt
using Test

const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_runtests()
    model = MOI.instantiate(NLopt.Optimizer; with_bridge_type = Float64)
    MOI.set(model, MOI.RawOptimizerAttribute("algorithm"), :LD_SLSQP)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            optimal_status = MOI.LOCALLY_SOLVED,
            atol = 1e-2,
            rtol = 1e-2,
            exclude = Any[
                MOI.ConstraintBasisStatus,
                MOI.ConstraintDual,
                MOI.ConstraintName,
                MOI.DualObjectiveValue,
                MOI.ObjectiveBound,
                MOI.NLPBlockDual,
                MOI.SolverVersion,
                MOI.VariableBasisStatus,
                MOI.VariableName,
            ]
        ),
        exclude = String[
            # TODO(odow): investigate failures. A lot of these are probably just
            # suboptimal solutions.
            "test_conic_NormInfinityCone_3",
            "test_conic_NormInfinityCone_INFEASIBLE",
            "test_conic_NormOneCone",
            "test_conic_NormOneCone_INFEASIBLE",
            "test_conic_linear_INFEASIBLE",
            "test_conic_linear_INFEASIBLE_2",
            "test_constraint_qcp_duplicate_diagonal",  # Fails on 32-bit
            "test_infeasible_",
            "test_linear_DUAL_INFEASIBLE",
            "test_linear_INFEASIBLE",
            "test_linear_add_constraints",
            "test_linear_VectorAffineFunction_empty_row",
            "test_model_ScalarAffineFunction_ConstraintName",
            "test_model_duplicate_ScalarAffineFunction_ConstraintName",
            "test_nonlinear_invalid",
            "test_objective_ObjectiveFunction_VariableIndex",
            "test_quadratic_SecondOrderCone_basic",
            "test_quadratic_constraint_integration",
            "test_quadratic_nonconvex_constraint_basic",
            "test_quadratic_nonconvex_constraint_integration",
            "test_unbounded_",
            "test_solve_TerminationStatus_DUAL_INFEASIBLE",
            "test_solve_result_index",
        ]
    )
    return
end

end

TestMOIWrapper.runtests()


text:
# arf_get_d and arf_get_si clash
get_d(x::ArfLike; rnd::Union{arb_rnd,RoundingMode} = RoundNearest) = get_d(x, rnd)
function get_d(x::ArfLike, rnd::Union{arb_rnd,RoundingMode})
    ccall(@libarb(arf_get_d), Float64, (Ref{arf_struct}, arb_rnd), x, rnd)
end
get_si(x::ArfLike; rnd::Union{arb_rnd,RoundingMode} = RoundNearest) = get_si(x, rnd)
function get_si(x::ArfLike, rnd::Union{arb_rnd,RoundingMode})
    ccall(@libarb(arf_get_si), Int, (Ref{arf_struct}, arb_rnd), x, rnd)
end

# arf_mul is given using a #DEFINE which doesn't work in Julia, implement this manually
arbcall"int arf_mul_rnd_any(arf_t z, arf_t x, arf_t y, slong prec, arf_rnd_t rnd)"
mul!(
    res::ArfLike,
    x::ArfLike,
    y::ArfLike,
    prec::Integer,
    rnd::Union{arb_rnd,RoundingMode},
) = mul_rnd_any!(res, x, y, prec, rnd)
mul!(
    res::ArfLike,
    x::ArfLike,
    y::ArfLike;
    prec::Integer = precision(res),
    rnd::Union{arb_rnd,RoundingMode} = RoundNearest,
) = mul!(res, x, y, prec, rnd)


text:
# Definition of all routes of the route tree.

function fade(x, y, speed, δt)
    return x, y + speed * δt
end

function post(x, y, speed, δt; cut = 10.0, angle = 45.0)
    if y < cut
        return x, y + speed *δt
    else
        return x + speed * δt * cosd(angle), y + speed * δt * sind(angle)
    end
end


text: coulomb = Z -> ℓ -> r -> -Z/r + ℓ*(ℓ+1)/2r^2


text:
## ISSUES:
# - Why does CategoricalLikelihood add a zero to the input?
# - No IndependentMOKernel for different kernels on each output
# - LFGP sampler doesn't work with multi-output
# - elbo/expected_loglik needs a more general approach for multi-output GPs

using ApproximateGPs
using GPLikelihoods
using Distributions
using StatsFuns
using ParameterHandling
using LinearAlgebra

# The CategoricalLikelihood in GPLikelihoods currently takes C-1 inputs for C
# classes and defines the last one as 0 - technically this is fine but seems a bit unnatural?
# Also, when using a different kernel per class I don't think this works?
struct CustomCategorical{Tl<:GPLikelihoods.AbstractLink} <: GPLikelihoods.AbstractLikelihood
    invlink::Tl
end

CustomCategorical(l=softmax) = CustomCategorical(Link(l))

(l::CustomCategorical)(f::AbstractVector{<:Real}) = Categorical(l.invlink(f))

function (l::CustomCategorical)(fs::AbstractVector)
    return Product(Categorical.(l.invlink.(fs)))
end

C = 3
N = 100

X = rand(N)
X_mo = KernelFunctions.MOInputIsotopicByOutputs(X, C)

lengthscale = 0.1

k = with_lengthscale(SEKernel(), lengthscale)
mo_kernel = IndependentMOKernel(k)

jitter = 1e-12

lik = CustomCategorical()

gp = GP(mo_kernel)

lgp = LatentGP(gp, lik, jitter)
lgp_fx = LatentGP(gp, lik, jitter)(X_mo)

# # TODO: rand doesn't work - f is a flat vector, needs to be reshaped before
# lik(f) is called This is essentially the same problem as the ELBO below but
# for LatentGPs. Should the dimensionality be defined by the likelihood (would
# need to specify num_classes in CategoricalLikelihood for example)
f_broken, y_broken = rand(lfgp)

f = rand(gp(X_mo, jitter))

rowvecs_y(f, out_dim) = RowVecs(reshape(f, :, out_dim))

rf = rowvecs_y(f, C)
y = rand(lik(rf))

onehot(y) = (sort(unique(y))' .== y)
# y = onehot(y)

M = N  # Number of inducing points

Z = copy(X)[1:M]
Z_mo = KernelFunctions.MOInputIsotopicByOutputs(Z, C)
fz = lgp(Z_mo).fx

q = MvNormal(zeros(M * C), Matrix{Float64}(I, M * C, M * C))

svgp = SVGP(fz, q)

post = posterior(svgp)

#!! A nasty hack to get the ELBO to work.
# The essential problem is that q_f is
# passed as a flat vector, but needs to be reshaped according to the number of
# classes. This approach works if the marginals q_f are independent for each
# class, but in general the marginals might need to be MvNormal.

# Where should this be fixed? Ideally a general solution that lets `marginals`
# work with MOGPs?
function expected_loglik_multiclass(mc::MonteCarlo, y::AbstractVector, q_f::AbstractVector{<:Normal}, lik; num_classes=C)
    q_f_ = reshape(q_f, :, num_classes)
    f_μ = mean.(q_f_)
    fs = f_μ .+ std.(q_f_) .* randn(eltype(f_μ), tuple(size(q_f_)..., mc.n_samples))
    p_ys = mapslices(lik, fs, dims=2)
    lls = loglikelihood.(p_ys, y)
    return sum(lls) / mc.n_samples
end

function ApproximateGPs.expected_loglik(mc::MonteCarlo, y::AbstractVector, q_f::AbstractVector{<:Normal}, lik::CustomCategorical)
    return expected_loglik_multiclass(mc, y, q_f, lik; num_classes=C)
end

elbo(svgp, lgp_fx, py; quadrature=MonteCarlo())


text:
module Orbitals

using Fermi
using Fermi.Geometry
using Fermi.GaussianBasis
using TensorOperations

export AbstractOrbitals, AtomicOrbitals, AbstractRestrictedOrbitals, AbstractUnrestrictedOrbitals
export GeneralRestrictedOrbitals
export RHFOrbitals
export UHFOrbitals

"""
    Fermi.AbstractOrbitals

Abstract type common to all orbitals

_struct tree:_

**AbstractOrbitals**  (Top level)
"""
abstract type AbstractOrbitals end

struct AtomicOrbitals <: AbstractOrbitals
    basisset::BasisSet
end

function AtomicOrbitals()
    return AtomicOrbitals(BasisSet())
end

function AtomicOrbitals(mol:: Molecule, basis::String)
    return AtomicOrbitals(BasisSet(mol, basis))
end

abstract type AbstractRestrictedOrbitals <: AbstractOrbitals end
abstract type AbstractUnrestrictedOrbitals <: AbstractOrbitals end

"""
    Fermi.GeneralRestrictedOrbitals

Dummy orbital struct that can be used to hold any type of Restricted orbitals.
Can be used to read in custom orbitals as

    Orbs = GeneralRestrictedOrbitals(X)

where X is an Array object. The molecule and basis will be deduced from the current options.
Alternatively, one can pass these informations explicitly:

    Orbs = GeneralRestrictedOrbitals(X, molecule=mol, name="myorbitals", basis="cc-pvdz")

`name` and `basis` are Strings, whereas `mol` is a `Fermi.Geometry.Molecule` object.

# Fields

    name      String with a label for the orbital
    basis     String indicating the basis set used to construct the orbitals
    molecule  Molecule object for which the orbitals were constructed
    C         NxN AbstractArray with the AO(lines) → MO(orbitals) coefficients

_struct tree:_

**GeneralRestrictedOrbitals** <: AbstractOrbitals
"""
struct GeneralRestrictedOrbitals{T} <: AbstractRestrictedOrbitals
    molecule::Molecule
    basis::String
    sd_energy::T
    C::AbstractArray{T,2}
end

function GeneralRestrictedOrbitals(C::AbstractArray{T,2}; mol=nothing, basis="undef", sd_energy=zero(T)) where T <: AbstractFloat

    mol === nothing ? mol = Fermi.Geometry.Molecule() : nothing
    basis == "undef" ? basis = Fermi.Options.get("basis") : nothing

    GeneralRestrictedOrbitals{T}(mol, basis, sd_energy, C)
end

"""
    Fermi.HartreeFock.RHFOrbitals

Struct holding information about Restricted Hartree--Fock orbitals

# Fields

    molecule   Molecule object associated with the orbitals
    basis      Basis set used to compute the orbitals
    eps        Orbital energies, i.e. diagonal of the Fock matrix
    C          Coefficients of the AO->MO transformation matrix
"""
struct RHFOrbitals <: AbstractRestrictedOrbitals
    molecule::Molecule
    basis::String
    eps::AbstractArray{Float64,1}
    sd_energy::Float64
    C::AbstractArray{Float64,2}
end

struct UHFOrbitals <: AbstractUnrestrictedOrbitals
    molecule::Molecule
    basis::String
    epsα::AbstractArray{Float64,1}
    epsβ::AbstractArray{Float64,1}
    sd_energy::Float64
    Cα::AbstractArray{Float64,2}
    Cβ::AbstractArray{Float64,2}
end

end # Module


text:
abstract type SDIRKMutableCache <: OrdinaryDiffEqMutableCache end

@cache mutable struct ImplicitEulerCache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
end

function alg_cache(alg::ImplicitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  γ, c = 1, 1
  @iipnlsolve

  atmp = similar(u,uEltypeNoUnits)

  ImplicitEulerCache(u,uprev,uprev2,du1,fsalfirst,k,z,dz,b,tmp,atmp,J,W,uf,jac_config,linsolve,nlsolver)
end

mutable struct ImplicitEulerConstantCache{F,N} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolver::N
end

function alg_cache(alg::ImplicitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  γ, c = 1, 1
  @oopnlsolve
  ImplicitEulerConstantCache(uf,nlsolver)
end

mutable struct ImplicitMidpointConstantCache{F,N} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolver::N
end

function alg_cache(alg::ImplicitMidpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  γ, c = 1//2, 1//2
  @oopnlsolve
  ImplicitMidpointConstantCache(uf,nlsolver)
end

@cache mutable struct ImplicitMidpointCache{uType,rateType,JType,WType,UF,JC,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  b::uType
  tmp::uType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
end

function alg_cache(alg::ImplicitMidpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  γ, c = 1//2, 1//2
  @iipnlsolve
  ImplicitMidpointCache(u,uprev,du1,fsalfirst,k,z,dz,b,tmp,J,W,uf,jac_config,linsolve,nlsolver)
end

mutable struct TrapezoidConstantCache{F,uType,tType,N} <: OrdinaryDiffEqConstantCache
  uf::F
  uprev3::uType
  tprev2::tType
  nlsolver::N
end

function alg_cache(alg::Trapezoid,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  γ, c = 1//2, 1
  @oopnlsolve

  uprev3 = u
  tprev2 = t

  TrapezoidConstantCache(uf,uprev3,tprev2,nlsolver)
end

@cache mutable struct TrapezoidCache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,tType,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  uprev3::uType
  tprev2::tType
  nlsolver::N
end

function alg_cache(alg::Trapezoid,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  γ, c = 1//2, 1
  @iipnlsolve

  uprev3 = zero(u)
  tprev2 = t
  atmp = similar(u,uEltypeNoUnits)

  TrapezoidCache(u,uprev,uprev2,du1,fsalfirst,k,z,dz,b,tmp,atmp,J,W,uf,jac_config,linsolve,uprev3,tprev2,nlsolver)
end

mutable struct TRBDF2ConstantCache{F,Tab,N} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::TRBDF2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tab = TRBDF2Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.d, tab.γ
  @oopnlsolve
  TRBDF2ConstantCache(uf,nlsolver,tab)
end

@cache mutable struct TRBDF2Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,Tab,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  zprev::uType
  zᵧ::uType
  z::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::TRBDF2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = TRBDF2Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.d, tab.γ
  @iipnlsolve

  atmp = similar(u,uEltypeNoUnits); zprev = similar(u); zᵧ = similar(u)

  TRBDF2Cache(u,uprev,du1,fsalfirst,k,zprev,zᵧ,z,dz,b,tmp,atmp,J,
              W,uf,jac_config,linsolve,nlsolver,tab)
end

mutable struct SDIRK2ConstantCache{F,N} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolver::N
end

function alg_cache(alg::SDIRK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  γ, c = 1, 1
  @oopnlsolve
  SDIRK2ConstantCache(uf,nlsolver)
end

@cache mutable struct SDIRK2Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
end

function alg_cache(alg::SDIRK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  γ, c = 1, 1
  @iipnlsolve

  z₁ = similar(u); z₂ = z
  atmp = similar(u,uEltypeNoUnits)

  SDIRK2Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,dz,b,tmp,atmp,J,
              W,uf,jac_config,linsolve,nlsolver)
end

mutable struct SSPSDIRK2ConstantCache{F,N} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolver::N
end

function alg_cache(alg::SSPSDIRK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  γ, c = 1//4, 1//1
  @oopnlsolve
  SSPSDIRK2ConstantCache(uf,nlsolver)
end

@cache mutable struct SSPSDIRK2Cache{uType,rateType,JType,WType,UF,JC,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  dz::uType
  b::uType
  tmp::uType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
end

function alg_cache(alg::SSPSDIRK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  γ, c = 1//4, 1//1
  @iipnlsolve

  z₁ = similar(u); z₂ = z
  atmp = similar(u,uEltypeNoUnits)

  SSPSDIRK2Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,dz,b,tmp,J,
                W,uf,jac_config,linsolve,nlsolver)
end

mutable struct Kvaerno3ConstantCache{UF,Tab,N} <: OrdinaryDiffEqConstantCache
  uf::UF
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Kvaerno3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tab = Kvaerno3Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, 2tab.γ
  @oopnlsolve
  Kvaerno3ConstantCache(uf,nlsolver,tab)
end

@cache mutable struct Kvaerno3Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,Tab,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Kvaerno3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Kvaerno3Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, 2tab.γ
  @iipnlsolve

  z₁ = similar(u); z₂ = similar(u); z₃ = similar(u); z₄ = z
  atmp = similar(u,uEltypeNoUnits)

  Kvaerno3Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,z₃,z₄,dz,b,tmp,atmp,J,
                W,uf,jac_config,linsolve,nlsolver,tab)
end

mutable struct Cash4ConstantCache{F,N,Tab} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Cash4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tab = Cash4Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ,tab.γ
  @oopnlsolve
  Cash4ConstantCache(uf,nlsolver,tab)
end

@cache mutable struct Cash4Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,N,Tab,F} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  z₅::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Cash4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Cash4Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ,tab.γ
  @iipnlsolve

  z₁ = similar(u); z₂ = similar(u); z₃ = similar(u); z₄ = similar(u); z₅ = z
  atmp = similar(u,uEltypeNoUnits)

  Cash4Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,z₃,z₄,z₅,dz,b,tmp,atmp,J,
              W,uf,jac_config,linsolve,nlsolver,tab)
end

mutable struct Hairer4ConstantCache{F,N,Tab} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Union{Hairer4,Hairer42},u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  if alg isa Hairer4
    tab = Hairer4Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  else
    tab = Hairer42Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  end
  γ, c = tab.γ, tab.γ
  @oopnlsolve
  Hairer4ConstantCache(uf,nlsolver,tab)
end

@cache mutable struct Hairer4Cache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,Tab,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType
  z₂::uType
  z₃::uType
  z₄::uType
  z₅::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::Union{Hairer4,Hairer42},u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  if alg isa Hairer4
    tab = Hairer4Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  else # Hairer42
    tab = Hairer42Tableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  end
  γ, c = tab.γ, tab.γ
  @iipnlsolve

  z₁ = similar(u); z₂ = similar(u); z₃ = similar(u); z₄ = similar(u); z₅ = z
  atmp = similar(u,uEltypeNoUnits)

  Hairer4Cache(u,uprev,du1,fsalfirst,k,z₁,z₂,z₃,z₄,z₅,dz,b,tmp,atmp,J,
               W,uf,jac_config,linsolve,nlsolver,tab)
end

@cache mutable struct ESDIRK54I8L2SACache{uType,rateType,uNoUnitsType,JType,WType,UF,JC,Tab,F,N} <: SDIRKMutableCache
  u::uType
  uprev::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z₁::uType; z₂::uType; z₃::uType; z₄::uType; z₅::uType; z₆::uType; z₇::uType; z₈::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::JType
  W::WType
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::ESDIRK54I8L2SA,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = ESDIRK54I8L2SATableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ, tab.γ
  @iipnlsolve

  z₁ = zero(u); z₂ = zero(u); z₃ = zero(u); z₄ = zero(u)
  z₅ = zero(u); z₆ = zero(u); z₇ = zero(u); z₈ = z
  atmp = similar(u,uEltypeNoUnits)

  ESDIRK54I8L2SACache(u,uprev,du1,fsalfirst,k,z₁,z₂,z₃,z₄,z₅,z₆,z₇,z₈,dz,b,tmp,atmp,J,
               W,uf,jac_config,linsolve,nlsolver,tab)
end

mutable struct ESDIRK54I8L2SAConstantCache{F,N,Tab} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolver::N
  tab::Tab
end

function alg_cache(alg::ESDIRK54I8L2SA,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tab = ESDIRK54I8L2SATableau(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  γ, c = tab.γ,tab.γ
  @oopnlsolve
  ESDIRK54I8L2SAConstantCache(uf,nlsolver,tab)
end


text:
function creatematerialdict()
  Dict(
    :steel => (sg=0.283, ea=29.0e6, er=29.0e6, gj=11.5e6),
    :monel => (sg=0.319, ea=26.0e6, er=26.0e6, gj=9.9e6),
    :titanium => (sg=0.164, ea=15.5e6, er=15.5e6, gj=5.9e6),
    :aluminium => (sg=0.1, ea=10.3e6, er=10.3e6, gj=3.9e6),
  )
end


text:
# Code from @carolsatye 20211112
# Minor modifications @PeterMaxwell 20211114 (compiles, but otherwise not tested correctness)

function mm_spec_s1flash_dimspec()::MPCCDimSpec
    dimspec = MPCCDimSpec(21, 2, 2, 18, 20, 1, 0)
    return dimspec
end
# Dimension of spatial variables
# Columns in F (number of complementarities?)
# Number of F functions per complementarity variable (usually 2)
# Number of equality constraints
# Number of inequality constraints
# Number of continuous parameters for PMPCC
# Number of discrete parameters for PMPCC

function mm_spec_s1flash_defn(x, pr, ps)

    # PM20211114: some modifications are a work-around for a bug in Symbolics.jl - https://github.com/JuliaSymbolics/Symbolics.jl/issues/319

    # PM20211114: Constants. PM, modified to include prefix 'c' because F was conflicting with the F constraints
    cF = 1
    cA = [3.97786; 4.00139; 3.93002]
    cB = [1064.84; 1170.874; 1182.774]
    cC = [-41.136; -48.833; -52.532]
    cp = 5
    cz = [0.5; 0.3; 0.2]

    w = Symbolics.scalarize( x[1:3] )
    y = Symbolics.scalarize( x[4:6] )
    L = x[7]
    V = x[8]
    K = Symbolics.scalarize( x[9:11] )
    psat = Symbolics.scalarize( x[12:14] )          #pure component vapor pressure
    log_psat = Symbolics.scalarize( x[15:17] )
    a = x[18]                                       #corrected vapor fraction
    at = x[19]                                      #vapor fration calculated with the Rachford-rice equation
    sv = x[20]                                      #slack for vapor fraction
    sl = x[21]                                      #slack for liquid fraction

    f = 0.5*(a*cF - V)^2

    ce = Vector{Num}(undef, 18)
    ce[1:3] = Symbolics.scalarize( log_psat .- cA .+ cB./(pr[1] .+ cC) )        #Antoine's equation
    ce[4:6] = Symbolics.scalarize( psat .- 10 .^ log_psat )
    ce[7:9] = Symbolics.scalarize( K*cp .- psat )                               #Raoult's law
    ce[10:12] = Symbolics.scalarize( K.*w .- y )
    ce[13:15] = Symbolics.scalarize( L*w .+ V*y .- cF*cz )                      #Balance per component
    ce[16] = Symbolics.scalarize( L + V - cF )                                  #Total balance
    ce[17] = Symbolics.scalarize( sum((cz .* (K .- 1)) ./ (1 .+ at*(K .- 1))) ) #Rachford-Rice equation
    ce[18] = Symbolics.scalarize( a - sv + sl - at )

    ci = Vector{Num}(undef, 20)
    ci[1:3] = Symbolics.scalarize( w )
    ci[4:6] = Symbolics.scalarize( 1 .- w )
    ci[7:9] = Symbolics.scalarize( y )
    ci[10:12] = Symbolics.scalarize( 1 .- y )
    ci[13:15] = Symbolics.scalarize( K )
    ci[16:18] = Symbolics.scalarize( psat )
    ci[19] = Symbolics.scalarize( a )
    ci[20] = Symbolics.scalarize( 1 - a )

    F = Matrix{Num}(undef, 2, 2)
    F[1,1] = Symbolics.scalarize( L )
    F[2,1] = Symbolics.scalarize( sl )
    F[1,2] = Symbolics.scalarize( V )
    F[2,2] = Symbolics.scalarize( sv )

    return MPCCDefinition(f, ce, ci, F)
end

function mm_spec_s1flash_nzmask()::MPCCModelNZMask
    nczmask = MPCCModelNZMask(
        # ce
        [   Set{Int64}([15]),
            Set{Int64}([16]),
            Set{Int64}([17]),
            Set{Int64}([12, 15]),
            Set{Int64}([13, 16]),
            Set{Int64}([14, 17]),
            Set{Int64}([9, 12]),
            Set{Int64}([10, 13]),
            Set{Int64}([11, 14]),
            Set{Int64}([1, 4, 9]),
            Set{Int64}([2, 5, 10]),
            Set{Int64}([3, 6, 11]),
            Set{Int64}([1, 4, 7, 8, 9]),
            Set{Int64}([2, 5, 7, 8, 10]),
            Set{Int64}([3, 6, 7, 8, 111]),
            Set{Int64}([7, 8]),
            Set{Int64}([9, 10, 11, 19]),
            Set{Int64}([18, 19, 20, 21])
        ],
        # ci
        [   Set{Int64}([1]),
            Set{Int64}([2]),
            Set{Int64}([3]),
            Set{Int64}([1]),
            Set{Int64}([2]),
            Set{Int64}([3]),
            Set{Int64}([4]),
            Set{Int64}([5]),
            Set{Int64}([6]),
            Set{Int64}([4]),
            Set{Int64}([5]),
            Set{Int64}([6]),
            Set{Int64}([9]),
            Set{Int64}([10]),
            Set{Int64}([11]),
            Set{Int64}([12]),
            Set{Int64}([13]),
            Set{Int64}([14]),
            Set{Int64}([18]),
            Set{Int64}([18]),
        ],
        # F
        [
            # F1
            [   Set{Int64}([7]), Set{Int64}([8]) ],
            # F2
            [   Set{Int64}([21]), Set{Int64}([20]) ]
        ] )
    return nczmask
end

function mm_spec_s1flash_testvectors()
    return Vector{MPCCModelTestVector}(undef, 0)
end

function mm_spec_s1flash_knownsols()

    return []
end

function mm_spec_s1flash_parameterisations(t)

    defn1 = MPCCParameterisationDefn(
            Vector{Num}([ t ]),
            (380.0, 400.0),               # tspan
            "Standard"
        )

    return  Vector{MPCCParameterisationDefn}( [ defn1 ] )
end


text:
export
    ActionContext,
    ContextFree,
    IntegratedContinuous,

    DriveAction,
    NextState,
    AccelTurnrate,
    AccelDesang,
    LatLonAccel,
    LaneFollowingAccel,
    AccelSteeringAngle,

    BicycleGeom,

    pull_action!,
    propagate

abstract ActionContext

type ContextFree <: ActionContext end
type IntegratedContinuous <: ActionContext
    Δt::Float64 # timestep
    n_integration_steps::Int # number of substeps taken during integration
end

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

abstract DriveAction
Base.length{A<:DriveAction}(a::Type{A}) = error("length not defined for DriveAction $a")
Base.convert{A<:DriveAction}(a::Type{A}, v::Vector{Float64}) = error("convert v → a not implemented for DriveAction $a")
Base.copy!(v::Vector{Float64}, a::DriveAction) = error("copy! not implemented for DriveAction $a")
Base.convert{A<:DriveAction}(::Type{Vector{Float64}}, a::A) = copy!(Array(Float64, length(A)), a)
propagate(veh::Vehicle, action::DriveAction, context::ActionContext, roadway::Roadway) = error("propagate not implemented for DriveAction $action and context $context")

immutable AccelTurnrate <: DriveAction
    a::Float64
    ω::Float64
end
Base.show(io::IO, a::AccelTurnrate) = @printf(io, "AccelTurnrate(%6.3f,%6.3f)", a.a, a.ω)
Base.length(::Type{AccelTurnrate}) = 2
Base.convert(::Type{AccelTurnrate}, v::Vector{Float64}) = AccelTurnrate(v[1], v[2])
function Base.copy!(v::Vector{Float64}, a::AccelTurnrate)
    v[1] = a.a
    v[2] = a.ω
    v
end
function propagate(veh::Vehicle, action::AccelTurnrate, context::IntegratedContinuous, roadway::Roadway)

    a = action.a # accel
    ω = action.ω # turnrate

    x = veh.state.posG.x
    y = veh.state.posG.y
    θ = veh.state.posG.θ
    v = veh.state.v

    δt = context.Δt / context.n_integration_steps

    for i in 1 : context.n_integration_steps
        x += v*cos(θ)*δt
        y += v*sin(θ)*δt
        θ += ω*δt
        v += a*δt
    end

    posG = VecSE2(x, y, θ)
    VehicleState(posG, roadway, v)
end
function Base.get(::Type{AccelTurnrate}, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    accel = get(ACC, rec, roadway, vehicle_index, pastframe)
    turnrate = get(TURNRATEG, rec, roadway, vehicle_index, pastframe)
    AccelTurnrate(accel, turnrate)
end
function pull_action!(::Type{AccelTurnrate}, a::Vector{Float64}, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    a[1] = get(ACC, rec, roadway, vehicle_index, pastframe)
    a[2] = get(TURNRATEG, rec, roadway, vehicle_index, pastframe)
    a
end

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

immutable AccelDesang <: DriveAction
    a::Float64
    ϕdes::Float64
end
Base.show(io::IO, a::AccelDesang) = @printf(io, "AccelDesang(%6.3f,%6.3f)", a.a, a.ϕdes)
Base.length(::Type{AccelDesang}) = 2
Base.convert(::Type{AccelDesang}, v::Vector{Float64}) = AccelDesang(v[1], v[2])
function Base.copy!(v::Vector{Float64}, a::AccelDesang)
    v[1] = a.a
    v[2] = a.ϕdes
    v
end
function propagate(veh::Vehicle, action::AccelDesang, context::IntegratedContinuous, roadway::Roadway)

    a = action.a # accel
    ϕdes = action.ϕdes # desired heading angle

    x = veh.state.posG.x
    y = veh.state.posG.y
    θ = veh.state.posG.θ
    v = veh.state.v

    δt = context.Δt/context.n_integration_steps

    for i in 1 : context.n_integration_steps

        posF = Frenet(VecSE2(x, y, θ), roadway)
        ω = ϕdes - posF.ϕ

        x += v*cos(θ)*δt
        y += v*sin(θ)*δt
        θ += ω*δt
        v += a*δt
    end

    posG = VecSE2(x, y, θ)
    VehicleState(posG, roadway, v)
end

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

"""
    LatLonAccel
Acceleration in the frenet frame
"""
immutable LatLonAccel <: DriveAction
    a_lat::Float64
    a_lon::Float64
end
Base.show(io::IO, a::LatLonAccel) = @printf(io, "LatLonAccel(%6.3f, %6.3f)", a.a_lat, a.a_lon)
Base.length(::Type{LatLonAccel}) = 2
Base.convert(::Type{LatLonAccel}, v::Vector{Float64}) = LatLonAccel(v[1], v[2])
function Base.copy!(v::Vector{Float64}, a::LatLonAccel)
    v[1] = a.a_lat
    v[2] = a.a_lon
    v
end
function propagate(veh::Vehicle, action::LatLonAccel, context::IntegratedContinuous, roadway::Roadway)

    a_lat = action.a_lat
    a_lon = action.a_lon

     v = veh.state.v
     ϕ = veh.state.posF.ϕ
    ds = v*cos(ϕ)
     t = veh.state.posF.t
    dt = v*sin(ϕ)

    ΔT = context.Δt
    ΔT² = ΔT*ΔT
    Δs = ds*ΔT + 0.5*a_lon*ΔT²
    Δt = dt*ΔT + 0.5*a_lat*ΔT²

    ds₂ = ds + a_lon*ΔT
    dt₂ = dt + a_lat*ΔT
    speed₂ = sqrt(dt₂*dt₂ + ds₂*ds₂)
    v₂ = sign(ds₂)*speed₂
    ϕ₂ = atan2(dt₂, ds₂) + (v₂ < 0.0)*π # handle negative speeds

    roadind = move_along(veh.state.posF.roadind, roadway, Δs)
    footpoint = roadway[roadind]
    posG = convert(VecE2, footpoint.pos) + polar(t + Δt, footpoint.pos.θ + π/2)
    posG = VecSE2(posG.x, posG.y, footpoint.pos.θ + ϕ₂)
    VehicleState(posG, roadway, v₂)

    # posF = Frenet(roadind, footpoint.s, t + Δt, ϕ₂)
    # VehicleState(posG, posF, v₂)
end
function Base.get(::Type{LatLonAccel}, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    accel_lat = get(ACCFT, rec, roadway, vehicle_index, pastframe)
    accel_lon = get(ACCFS, rec, roadway, vehicle_index, pastframe)
    LatLonAccel(accel_lat, accel_lon)
end
function pull_action!(::Type{LatLonAccel}, a::Vector{Float64}, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    a[1] = get(ACCFT, rec, roadway, vehicle_index, pastframe)
    a[2] = get(ACCFS, rec, roadway, vehicle_index, pastframe)
    a
end

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

"""
    LaneFollowingAccel
Longitudinal acceleration
"""
immutable LaneFollowingAccel <: DriveAction
    a::Float64
end
Base.show(io::IO, a::LaneFollowingAccel) = @printf(io, "LaneFollowingAccel(%6.3f)", a.a)
Base.length(::Type{LaneFollowingAccel}) = 1
Base.convert(::Type{LaneFollowingAccel}, v::Vector{Float64}) = LaneFollowingAccel(v[1])
function Base.copy!(v::Vector{Float64}, a::LaneFollowingAccel)
    v[1] = a.a
    v
end
function propagate(veh::Vehicle, action::LaneFollowingAccel, context::IntegratedContinuous, roadway::Roadway)

    a_lon = action.a

    ds = veh.state.v

    ΔT = context.Δt
    ΔT² = ΔT*ΔT
    Δs = ds*ΔT + 0.5*a_lon*ΔT²

    v₂ = ds + a_lon*ΔT

    roadind = move_along(veh.state.posF.roadind, roadway, Δs)
    posG = roadway[roadind].pos
    VehicleState(posG, roadway, v₂)
end
###############

immutable NextState <: DriveAction
    s::VehicleState
end
Base.show(io::IO, a::NextState) = print(io, "NextState(", a.s, ")")
Base.length(::Type{NextState}) = 11
function Base.convert(::Type{NextState}, v::Vector{Float64})
    VehicleState(VecSE2(v[1],v[2],v[3]), # x, y, θ
                 Frenet(
                    RoadIndex(
                        CurveIndex(round(Int, v[4]), v[5]),
                        LaneTag(round(Int, v[6]), round(Int, v[7])),
                    ),
                    v[8], v[9], v[10] # s, t, ϕ
                 ),
                 v[11] # speed
    )
end
function Base.copy!(v::Vector{Float64}, a::NextState)
    v[1] = a.s.posG.x
    v[2] = a.s.posG.y
    v[2] = a.s.posG.θ
    v[2] = a.s.posF.roadind.ind.i
    v[2] = a.s.posG.roadind.ind.t
    v[2] = a.s.posG.roadind.tag.segment
    v[2] = a.s.posG.roadind.tag.lane
    v[2] = a.s.posG.s
    v[2] = a.s.posG.t
    v[2] = a.s.posG.ϕ
    v
end
propagate{C<:ActionContext}(veh::Vehicle, action::NextState, context::C, roadway::Roadway) = action.s

###########

@with_kw type BicycleGeom
    wheel_base = 3.0 # distance between front and rear axle [m]
    wheel_base_offset = -1.0 # distance along vehicle from state position where the rear axle is [m]
end

immutable AccelSteeringAngle <: DriveAction
    a::Float64 # accel [m/s²]
    δ::Float64 # steering angle [rad]
end
Base.show(io::IO, a::AccelSteeringAngle) = @printf(io, "AccelSteeringAngle(%6.3f,%6.3f)", a.a, a.δ)
Base.length(::Type{AccelSteeringAngle}) = 2
Base.convert(::Type{AccelSteeringAngle}, v::Vector{Float64}) = AccelSteeringAngle(v[1], v[2])
function Base.copy!(v::Vector{Float64}, a::AccelSteeringAngle)
    v[1] = a.a
    v[2] = a.δ
    v
end
function AutomotiveDrivingModels.propagate(veh::Vehicle, action::AccelSteeringAngle, context::IntegratedContinuous, roadway::Roadway,
    geom::BicycleGeom=BicycleGeom(),
    )

    L = geom.wheel_base
    l = geom.wheel_base_offset

    a = action.a # accel [m/s²]
    δ = action.δ # steering wheel angle [rad]

    x = veh.state.posG.x
    y = veh.state.posG.y
    θ = veh.state.posG.θ
    v = veh.state.v

    Δt = context.Δt

    s = v*Δt + a*Δt*Δt/2 # distance covered
    v′ = v + a*Δt

    if abs(δ) < 0.01 # just drive straight
        posG = veh.state.posG + polar(s, θ)
    else # drive in circle

        R = L/tan(δ) # turn radius

        β = s/R
        xc = x - R*sin(θ) + l*cos(θ)
        yc = y + R*cos(θ) + l*sin(θ)

        θ′ = mod(θ+β, 2π)
        x′ = xc + R*sin(θ+β) - l*cos(θ′)
        y′ = yc - R*cos(θ+β) - l*sin(θ′)

        posG = VecSE2(x′, y′, θ′)
    end

    VehicleState(posG, roadway, v′)
end

# function Base.get(::Type{AccelSteeringAngle}, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0,
#     geom::BicycleGeom=BicycleGeom(),
#     )
#     accel = get(ACC, rec, roadway, vehicle_index, pastframe)
#     turnrate = get(TURNRATEG, rec, roadway, vehicle_index, pastframe)
#     AccelSteeringAngle(accel, turnrate)
# end
# function AutomotiveDrivingModels.pull_action!(::Type{AccelSteeringAngle}, a::Vector{Float64}, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
#     a[1] = get(ACC, rec, roadway, vehicle_index, pastframe)
#     a[2] = get(TURNRATEG, rec, roadway, vehicle_index, pastframe)
#     a
# end


text:
# # Bump On Tail
#
#md # [`notebook`](@__NBVIEWER_ROOT_URL__notebooks/bump_on_tail.ipynb)
#

import Splittings: advection!, UniformMesh
import Splittings: compute_rho, compute_e
import Splittings: CubicSpline, @Strang

using Plots
using LaTeXStrings

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

function push_t!( f, mesh1, v, dt )

    advection!( f, mesh1, v, 0.5dt, CubicSpline(), 1)

end

function push_v!( f, mesh1, mesh2, nrj, dt )

    rho = compute_rho(mesh2, f)
    e   = compute_e(mesh1, rho)
    advection!( f, mesh2, e, dt, CubicSpline(), 2)
    push!(nrj, 0.5*log(sum(e.*e)*mesh1.step))

end

function vlasov_poisson(mesh1  :: UniformMesh,
                        mesh2  :: UniformMesh,
                        f      :: Array{Float64,2},
                        nstep  :: Int64,
                        dt     :: Float64)

    x = mesh1.points
    v = mesh2.points

    nrj = Float64[]
    for istep in 1:nstep

        @Strang(  push_t!( f, mesh1, v, dt ),
                  push_v!( f, mesh1, mesh2, nrj, dt ))

    end
    nrj

end

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

α = 0.03
kx  = 0.3
x1min, x1max = 0.0, 2π / kx
n1, n2 = 32, 64
x2min, x2max = -9., 9.
mesh1 = UniformMesh(x1min, x1max, n1)
mesh2 = UniformMesh(x2min, x2max, n2)
f = zeros(Float64,(mesh1.length,mesh2.length))
for (i,x) in enumerate(mesh1.points), (j,v) in enumerate(mesh2.points)
    f[i,j]  = (1.0+α*cos(kx*x)) / (10*sqrt(2π)) * (9*exp(-0.5*v^2)+2*exp(-2*(v-4.5)^2))
end

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

nstep = 500
t = range(0.0, stop=50.0, length=nstep)
dt = t[2]
@elapsed nrj = vlasov_poisson( mesh1, mesh2, f, nstep, dt)

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

plot(t, nrj, label=L"\frac{1}{2} \log(∫e²dx)")
savefig("bot-plot.png"); nothing # hide

using Test                       #src
@testset "Bump On Tail" begin    #src
@test length(nrj) > 0            #src
end                              #src

#md # ![](bot-plot.png)


text:
using IncrementalInference
using Test

@testset "test the basics" begin

fg = initfg()

addVariable!(fg, :x1, ContinuousScalar)
addVariable!(fg, :x2, ContinuousScalar)
addFactor!(fg, [:x1;:x2], LinearConditional(Normal()), autoinit=false)
addFactor!(fg, [:x2], Prior(Normal()), autoinit=false)

@test exists(fg, :x1)

@test !exists(fg, :l13)

end


text:
"""

    This example show's how to use Term's errors functionality

Term provides functionality to produce styled and informative error
messages and tracestacks. This is not active by default as not everybody
might want this, but it's easier to turn on.

NOTE: this is a preliminary feature, there might be bugs and not all error
types may be correctly represented. While generally useful, it might be best
to hold back on this feature for code that needs to reliably print acurate
error information.
"""

import Term: install_stacktrace
import Term: hLine

install_stacktrace()

print(hLine("Fancy Errors", style = "bold blue"))

"""
Done!

Now error messages will be fancy!

To run this example uncomment one of the lines below to find the
corresponding error message
"""

# ------------  MethodError
# 1 - "a"

# ------------  DomainError
# sqrt(-1)

# ------------  AssertionError
# import Term: Panel
# Panel("mytext", title="this title is waaaay too long!!!")

# ------------  UndefVarError
# println(x)

# ------------  BoundsError
# v = collect(1:10)
# v[20]

# ------------  DivideError
# div(2, 0)

# ------------  StackOverflowError
# a() = b()
# b() = a()
# a()

# ------------  KeyError
# mydict = Dict(:a=>"a", :b=>"b")
# mydict["a"]

# ------------  InexactError
# Int(2.5)

# ------------  UndefKeywordError
# function my_func(;my_arg::Int)
#     return my_arg + 1
# end
# my_func()

# ------------  DimensionMismatch
a = zeros(20, 20)
b = zeros(5, 4)
a .+ b


text:
#=
# TODO: NEED VARIABLE SUBTYPES...

is_convex(::typeof(VARIABLE), x::Interval{T}) where T = true
is_concave(::typeof(VARIABLE), x::Interval{T}) where T = true

is_loglogconvex(::typeof(VARIABLE), x::Interval{T}) where T = true
is_loglogconcave(::typeof(VARIABLE), x::Interval{T}) where T = true

is_increasing(::typeof(VARIABLE), x::Interval{T}) where T = true
is_decreasing(::typeof(VARIABLE), x::Interval{T}) where T = false

is_positive(::typeof(VARIABLE), x::Interval{T}) where T = x >= zero(T)
is_negative(::typeof(VARIABLE), x::Interval{T}) where T = x <= zero(T)
is_locked(::typeof(VARIABLE), x::Interval{T}) where T = one(T) ∉ x
=#

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
using JuliaOSSTemplate
using Test

@test add(3, 4) == 7

max_stars_repo_name: utilForever/julia-oss-template

max_issues_repo_name: utilForever/julia-oss-template

max_forks_repo_name: utilForever/julia-oss-template


text:
_scip = JuMP.optimizer_with_attributes(SCIP.Optimizer, "display/verblevel" => 0)

@testset "run MINLP model" begin
    #       7    9      even arc numbers for
    #   () - d2 - ()    reversed arcs
    #   /1   /3   /5
    #  s1 - () - d1
    #    11   13

    # slightly irregular grid
    nodes = [Node(0,0), Node(45,0), Node(25, 22)]
    demand = [-50, 20, 30]
    bounds = fill(Bounds(60,80), 3)
    diams = [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]]

    topo = Topology([Node(t...) for t in [(0,22), (0,0), (25,22), (25,0), (45,22), (45,0)]],
                    [Arc(t...) for t in [(1,2), (2,1), (3,4), (4,3), (5,6), (6,5),
                                         (1,3), (3,1), (3,5), (5,3),
                                         (2,4), (4,2), (4,6), (6,4)]])

    @testset "low flow: very easy instance" begin
        inst = Instance(nodes, 1*demand, bounds, diams)

        result = optimize(inst, topo, GndStr.MINLP(_scip))
        @test result.status == MOI.OPTIMAL

        zsol = result.solution.zsol
        @test zsol[4,1] == true
        @test zsol[11,1] == true
        @test zsol[13,1] == true
        @test sum(zsol) == 3

        @test result.dualbound ≈ 67.0
    end

    @testset "medium flow: difficult instance" begin
        inst = Instance(nodes, 5*demand, bounds, diams)

        result = optimize(inst, topo, GndStr.MINLP(_scip))
        @test result.status == MOI.OPTIMAL

        zsol = result.solution.zsol
        @test zsol[4,1] == true
        @test zsol[11,2] == true
        @test zsol[13,1] == true
        @test sum(zsol) == 3

        @test result.dualbound ≈ 72.0
    end

    @testset "high flow on triangle: infeasible" begin
        inst3 = Instance([Node(0,0), Node(50,0)],
                         20*[-50, 50],
                         [Bounds(60,80), Bounds(60,80)],
                         [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2)]])
        topo3 = Topology([Node(0,0), Node(50,0), Node(30, 40)],
                         [Arc(1,3), Arc(1,2), Arc(2,3)])

        result = optimize(inst3, topo3, GndStr.MINLP(_scip))
        @test result.status == MOI.INFEASIBLE
        @test result.solution == nothing
        @test result.dualbound == Inf
    end

    @testset "on the example with subproblem/relaxation gap" begin
        #     __s1__  __s2
        #   t3      t4
        nodes = [Node(100,0), Node(300,0), Node(0,0), Node(200,0)]
        arcs = [Arc(1,3), Arc(1,4), Arc(2,4)]
        topo = Topology(nodes, arcs)

        demand = [-600, -400, 400, 600]
        bounds = fill(Bounds(40, 80), 4)
        diams = [Diameter(1.0, 1.0), Diameter(2.0, 2.0)]
        inst = Instance(nodes, demand, bounds, diams, ploss_coeff_nice)

        result = optimize(inst, topo, GndStr.MINLP(_scip))
        @test result.status == MOI.OPTIMAL
        zsol = result.solution.zsol
        @test sum(zsol[:,2]) == 1
        @test sum(zsol[:,1]) == 2 # solution not quite uniqe
        @test result.dualbound ≈ 400.0
    end

    @testset "difficult instance with disconnected candidates" begin
        inst = Instance([Node(0,0), Node(100,0), Node(200,100)],
                        [800, -900, 100],
                        fill(Bounds(40,80), 3),
                        [Diameter(t...) for t in [(1.0, 1.0), (2.0, 3.2)]],
                        ploss_coeff_nice)
        topo = squaregrid(2, 3, 100.0, antiparallel=true)

        result = optimize(inst, topo, GndStr.MINLP(_scip))
        @test result.status == MOI.OPTIMAL

        zsol = result.solution.zsol
        @test sum(zsol) == 3

        @test result.dualbound ≈ 520.0
    end
end


text:
# pseudocode by Baurzhan Muftakhidinov
# Read Eclipse sim software output formatted in Fortran Arrays
# seek(s, 0)

function readlong(file)
    return ntoh(read(file, Int32))
end

function readdouble(file)
    return ntoh(read(file, Float64))
end

function readreal(file)
    return ntoh(read(file, Float32))
end

function readchar(file)
    # 4 chars to form a CHAR
    return strip(String(read(file, 4)))
end

function readstring(file)
    # 8 chars to form a STRING
    return strip(String(read(file, 8)))
end

struct EclSection
    kwd::String
    dim::Int32
    dtype::String
    endi2::Int32
    data::Vector{}
    function EclSection(file::IOStream)
        kwd = readstring(file)
        #println("DBG : Section found: $kwd")
        # Dimension, number of data points
        dim = readlong(file)
        dtype = readchar(file)
        # endianness sign
        endi2 = readlong(file)
        data = []
        while length(data) < dim
            tmp = []
            NBYTES = readlong(file)
            if dtype == "CHAR"
                num = Int(NBYTES/8)
                fun = readstring
            elseif dtype == "INTE"
                num = Int(NBYTES/4)
                fun = readlong
            elseif dtype == "DOUB"
                num = Int(NBYTES/8)
                fun = readdouble
            elseif dtype == "REAL"
                num = Int(NBYTES/4)
                fun = readreal
            elseif dtype == "LOGI"
                num = Int(NBYTES/4)
                fun = readlong
            else
                println("Unknown dtype", dtype)
                exit
            end
            for i in 1:num
                push!(tmp, fun(file))
            end
            NBYTES2 = readlong(file)  # extra endianness sign
            append!(data, tmp)
        end
        new(kwd, dim, dtype, endi2, data)
    end
end

struct EclFile
    data::Vector{EclSection}
    function EclFile(filename::String)
        s = open(filename, "r")
        data = EclSection[]

        while ! eof(s)
            readlong(s)  #endianness sign
            push!(data, EclSection(s))
        end

        #readlong(s)
        close(s)
        new(data)
    end
end

#a = EclFile("NORNE/NORNE_ATW2013.SMSPEC")
#b = EclFile("NORNE/NORNE_ATW2013.UNSMRY")
#c = EclFile("NORNE_ATW2013.UNRST")
#println(c.data[23]) #.data[390:400])


text:
rich_scaled(x::Float64) = 0.95 / (1 + exp(-x))

function run_model(ps::ParameterSetting, scale::Float64, window::Int, numb_updates::Int; init_w=0.5)
    @unpack max_tspan, pr_quick_rec, trans_rate, drop, pr_death, N, λ, hhs = ps
    ppl, grph = create_model(N, λ, hhs)
    s = vcat(dropdims(sum(retrieve_data(ppl), dims=1), dims=1), init_w)
    daily_sus = Vector{Int64}(undef, numb_updates)
    daily_inf = Vector{Int64}(undef, numb_updates)
    daily_rec = Vector{Int64}(undef, numb_updates)
    daily_dead = Vector{Int64}(undef, numb_updates)
    wgts = Vector{Float64}(undef, numb_updates)

    @inbounds for i in 1:window
        model_update!(ppl, grph, init_w, max_tspan, pr_quick_rec, trans_rate, drop, pr_death)
        daily_sus[i] = sum(Int.([ ppl[j].status for j in 1:length(ppl) ] .== :S))
        daily_inf[i] = sum(Int.([ ppl[j].status for j in 1:length(ppl) ] .== :I))
        daily_rec[i] = sum(Int.([ ppl[j].status for j in 1:length(ppl) ] .== :R))
        daily_dead[i] = sum(Int.([ ppl[j].status for j in 1:length(ppl) ] .== :D))
        wgts[i] = init_w
    end

    ps = 0.
    @inbounds for i in window + 1:numb_updates
        x = [ones(window) 1:window]
        y = daily_inf[i - window:i - 1]
        _, cf = x \ y
        ps = ps + (cf * scale)
        w = rich_scaled(ps)
        model_update!(ppl, grph, w, max_tspan, pr_quick_rec, trans_rate, drop, pr_death)
        daily_sus[i] = sum(Int.([ ppl[j].status for j in 1:length(ppl) ] .== :S))
        daily_inf[i] = sum(Int.([ ppl[j].status for j in 1:length(ppl) ] .== :I))
        daily_rec[i] = sum(Int.([ ppl[j].status for j in 1:length(ppl) ] .== :R))
        daily_dead[i] = sum(Int.([ ppl[j].status for j in 1:length(ppl) ] .== :D))
        wgts[i] = w
    end
    return vcat(s', hcat(daily_sus, daily_inf, daily_rec, daily_dead, wgts))
end


text:
using Distributions, Documenter

makedocs(
    format = Documenter.HTML(),
    sitename = "Modelselection",
    modules = [GLM],
    pages = [
        "Home" => "index.md"
    ],
    debug = true,
)

deploydocs(
    repo   = "github.com/ParallelGSReg/ModelSelection.jl.git",
)


text:
mutable struct TranslationParameters

    vmax::Float64
    tau_factor::Float64
    KL::Float64

    function TranslationParameters()
		this = new()
	end
end


text:
# run tests in parallel
include("runtests_setup.jl")

comm = MPI.COMM_WORLD
global const comm_size = MPI.Comm_size(MPI.COMM_WORLD)
global const comm_rank = MPI.Comm_rank(MPI.COMM_WORLD)

# size of the system owned by this process (ie. local sys size)
sys_size = GridapDistributedPETScWrappers.C.PetscInt(3)

# create a 3x3 block for each process
# create these with smallest precision, so they can be promoted
tmp3 = convert(Array{ComplexF64, 1}, [1.0 + 0im; 2.0 + 1.0im; 3.0 + 2.0im])
tmp4 = convert( Array{ComplexF64, 2}, [1.0 + 1im   2 + 2im  3 + 3im; 4 + 4im  5 + 5im 7 + 7im; 7 + 7im 8 + 8im 9 + 9im])

for (i, ST) in enumerate(GridapDistributedPETScWrappers.C.petsc_type)
  if GridapDistributedPETScWrappers.have_petsc[i]
    # @testset "Scalar type $ST" begin # uncomment when nested test results can be printed

    # Have both local and global versions of system for testing
    global A_julia = zeros(ST, sys_size, sys_size)
    global rhs = zeros(ST, sys_size)
    global A_julia_global = zeros(ST, comm_size*sys_size, comm_size*sys_size)
    global rhs_global = zeros(ST, comm_size*sys_size)

    # create right hand side
    for i=1:sys_size
      rhs[i] = convert(ST, RC(tmp3[i]))

      for j=0:(comm_size - 1)  # global
        idx = i + j*sys_size
        rhs_global[idx] = convert(ST, RC(tmp3[i]))
      end

    end

    # create matrix A
    for i=1:sys_size
      for j=1:sys_size
        A_julia[i,j] = convert(ST, RC(tmp4[i,j]))

        for k=0:(comm_size - 1)
          idxm = i + k*sys_size
          idxn = j + k*sys_size
          A_julia_global[idxm, idxn] = convert(ST, RC(tmp4[i,j]) )
        end

      end
    end

    x_julia = A_julia\rhs

    include("vecp.jl")
    #include("c_funcs.jl")
    # end
  end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing

function calc_model(mesh_name; with_volume_load=false, debug_print=false)
    meshfile = Pkg.dir("JuliaFEM")*"/test/testdata/3d_block.med"
    mesh = aster_read_mesh(meshfile, mesh_name)

    block = Problem(Elasticity, "BLOCK", 3)
    block.properties.finite_strain = false
    block.properties.geometric_stiffness = false
    block.elements = create_elements(mesh, "BLOCK")
    update!(block, "youngs modulus", 288.0)
    update!(block, "poissons ratio", 1/3)
    with_volume_load && update!(block, "displacement load 3", 576.0)

    traction = Problem(Elasticity, "traction force", 3)
    traction.properties.finite_strain = false
    traction.properties.geometric_stiffness = false
    traction.elements = create_elements(mesh, "LOAD")
    update!(traction, "displacement traction force 3", 288.0)

    bc = Problem(Dirichlet, "symmetry boundary condition", 3, "displacement")
    symyz = create_elements(mesh, "SYMYZ")
    symxz = create_elements(mesh, "SYMXZ")
    symxy = create_elements(mesh, "SYMXY")
    update!(symyz, "displacement 1", 0.0)
    update!(symxz, "displacement 2", 0.0)
    update!(symxy, "displacement 3", 0.0)
    push!(bc, symyz, symxz, symxy)

    solver = LinearSolver("Solver block problem")
    push!(solver, block, traction, bc)
    solver()

    max_u = maximum(block.assembly.u)
    nu = round(Int, length(block.assembly.u)/3)
    u = reshape(block.assembly.u, 3, nu)
    if debug_print
        f = reshape(full(block.assembly.f), 3, nu)
        dump(round(u', 5))
        dump(round(f', 5))
        info("max |u| = $max_u")
    end
    return block, u
end

@testset "test 3d block HEX8" begin
    block, u = calc_model("BLOCK_HEX8"; with_volume_load=true)
    @test isapprox(maximum(u), 2.0)
end

@testset "test 3d block TET4" begin
#   block, u = calc_model("BLOCK_TET4", :TE4, :TR3; with_volume_load=true)
#   @test isapprox(maximum(u), 2.1329516539440205)
    block, u = calc_model("BLOCK_TET4"; with_volume_load=false)
    @test isapprox(maximum(u), 1.0)
end

@testset "test 3d block TET10" begin
    block, u = calc_model("BLOCK_TET10"; with_volume_load=false)
#    @test isapprox(maximum(u), 2.13656216413056)
    @test isapprox(maximum(u), 1.0)
end


text:
using GeoGreensFunctions

const TESTDIR = @__DIR__
const TESTFILES = filter(x -> startswith(x, "test_") && endswith(x, ".jl"), readdir(TESTDIR))
foreach(x -> include(joinpath(TESTDIR, x)), TESTFILES)


text:
function check_dataframe_has_name(x::DataFrame)
    return size(x)[2] >= 1 && names(x)[1] == "name"
end

function check_dataframe_firstcol(x::DataFrame)
    first = x[!,1]
    return isa(first, AbstractVector{<:AbstractString}) || isa(first, Vector{Nothing})
end

function check_dataframe_in_constructor(x::DataFrame, message::String)
    if !check_dataframe_has_name(x)
        throw(ArgumentError("first column of '" * message * "' should be 'name'"))
    end

    if !check_dataframe_firstcol(x)
        throw(ArgumentError("first column of '" * message * "' should contain strings or nothings"))
    end
end

function check_assay_dimensions(dims::Tuple{Vararg{Int}}, ref::Tuple{Vararg{Int}})
    return length(dims) >= 2 && dims[1] == ref[1] && dims[2] == ref[2]
end

"""
The `SummarizedExperiment` class is a Bioconductor container for matrix-like data with annotations on the rows and columns.
It is the data structure underlying analysis workflows for many genomics data modalities,
ranging from microarrays, bulk and single-cell RNA sequencing, ChIP-seq, epigenomics and beyond.

Any number of arrays (also known as "assays") can be stored in the container,
provided they are assigned to unique names and all have the same extents for the first two dimensions.
This reflects the fact that we often have multiple experimental readouts of the same shape, e.g., raw counts, normalized values, quality metrics.
These assays are held as an `OrderedDict` so the order of their addition is respected.

The row and column annotations are stored as `DataFrame`s, with number of rows equal to the number of assay rows and columns, respectively.
Any number and type of columns may be present in each `DataFrame`,
with the only constraint being that the first column must be a `"name"` column of strings containing the feature/sample names.
If no names are present, the `"name"` column must contain `nothing`s.

Each instance may also contain arbitrary metadata not associated with the rows or columns.
"""
mutable struct SummarizedExperiment
    assays::OrderedDict{String,AbstractArray}
    rowdata::DataFrame
    coldata::DataFrame
    metadata::Dict{String,Any}

    @doc """
        SummarizedExperiment()

    Create an empty `SummarizedExperiment` with no assays and empty row/column annotations.

    # Examples
    ```jldoctest
    julia> using SummarizedExperiments

    julia> SummarizedExperiment()
    0x0 SummarizedExperiment
      assays(0):
      rownames:
      rowdata(1): name
      colnames:
      coldata(1): name
      metadata(0):
    ```
    """
    function SummarizedExperiment()
        dummy = DataFrame(name = String[])
        new(
           OrderedDict{String,AbstractArray}(),
           dummy,
           deepcopy(dummy),
           Dict{String,Any}()
        )
    end

    @doc """
        SummarizedExperiment(assays)

    Create an instance of a `SummarizedExperiment` with the supplied assays.

    All entries of `assays` should have the same extents for the first two dimensions.
    However, they can otherwise have any number of other dimensions.
    Each assay can be of different type.
    `assays` should contain at least one assay matrix.

    For the `coldata` and `rowdata`, an empty `DataFrame` is created with a `"name"` column containing all `nothing`s.

    # Examples
    ```jldoctest
    julia> using SummarizedExperiments

    julia> assays = OrderedDict{String, AbstractArray}(
              "foobar" => [[1,2] [3,4] [5,6]],
              "whee" => [[1.2,2.3] [3.4,4.5] [5.6,7.8]]);

    julia> x = SummarizedExperiment(assays)
    2x3 SummarizedExperiment
      assays(2): foobar whee
      rownames:
      rowdata(1): name
      colnames:
      coldata(1): name
      metadata(0):
    ```
    """
    function SummarizedExperiment(assays::OrderedDict{String, AbstractArray})
        if length(assays) < 1
            throw(ErrorException("expected at least one array in 'assays' if 'rowdata' or 'coldata' are not supplied"))
        end

        refdims = size(first(assays).second)
        if length(refdims) < 2
            throw(DimensionMismatch("'assays' should contain arrays with 2 or more dimensions"))
        end

        first_name = first(assays).first
        for (key, val) in assays
            if !check_assay_dimensions(size(val), refdims)
                throw(DimensionMismatch("'assays[" * repr(key) * "]' and 'assays[" * repr(first_name) *
                    "]' should have the same extents for the first 2 dimensions"))
            end
        end

        rowdata = DataFrame(name = Vector{Nothing}(undef, refdims[1]))
        coldata = DataFrame(name = Vector{Nothing}(undef, refdims[2]))
        new(assays, rowdata, coldata, Dict{String,Any}())
    end

    @doc """
        SummarizedExperiment(assays, rowdata, coldata, metadata = Dict{String,Any}())

    Create an instance of a `SummarizedExperiment` with the supplied assays and the row/column annotations.

    All entries of `assays` should have the same extents for the first two dimensions.
    However, they can otherwise have any number of other dimensions.
    Each assay can be of different type.

    For `rowdata`, the number of rows must be equal to the extent of the first dimension for each entry in `assays`.
    Similarly, for `coldata`, the number of rows must be equal to the extent of the second dimension.
    In both cases, the first column must be called `"name"` and contain a `Vector` of `String`s or `Nothing`s (if no names are available).

    `assays` may also be empty.

    # Examples
    ```jldoctest
    julia> using SummarizedExperiments

    julia> assays = OrderedDict{String, AbstractArray}(
              "foobar" => [[1,2] [3,4] [5,6]],
              "whee" => [[1.2,2.3] [3.4,4.5] [5.6,7.8]]);

    julia> rowdata = DataFrame(
              name = [ "X", "Y" ],
              type = ["protein", "transcript"]);

    julia> coldata = DataFrame(
              name = [ "a", "b", "c" ],
              treatment = ["normal", "drug1", "drug2"]);

    julia> x = SummarizedExperiment(assays, rowdata, coldata)
    2x3 SummarizedExperiment
      assays(2): foobar whee
      rownames: X Y
      rowdata(2): name type
      colnames: a b c
      coldata(2): name treatment
      metadata(0):
    ```
    """
    function SummarizedExperiment(
            assays::OrderedDict{String, AbstractArray},
            rowdata::DataFrame,
            coldata::DataFrame,
            metadata::Dict{String,Any} = Dict{String,Any}()
        )

        refdims = (size(rowdata)[1], size(coldata)[1])

        for (key, val) in assays
            if !check_assay_dimensions(size(val), refdims)
                throw(DimensionMismatch("dimensions of 'assays[" * repr(key) * "]' are not consistent with 'rowdata' or 'coldata'"))
            end
        end

        check_dataframe_in_constructor(rowdata, "rowdata")
        check_dataframe_in_constructor(coldata, "coldata")
        new(assays, rowdata, coldata, metadata)
    end
end


text:
using CountyClusteringNotebook
using Test

@testset "CountyClusteringNotebook.jl" begin
    # Write your tests here.
end


text:
# Basedir is expected at this location relative to home for file transfer
const relBaseDir =
	joinpath("Documents", "projects", "p2019", "college_stratification");

"""
	$(SIGNATURES)

Points to `college_stratification`. All project files are in this directory.

Currently not relocatable.
"""
function base_dir(computer :: Union{Symbol, Computer} = :current)
	return joinpath(home_dir(computer), relBaseDir)
end

"""
	$(SIGNATURES)

Directory for Julia shared code. Perhaps not needed?
"""
function julia_dir(computer :: Union{Symbol, Computer} = :current)
	return joinpath(home_dir(computer),  "Documents",  "julia")
end

"""
	$(SIGNATURES)

Code and other material that get synced to `github` are here.
"""
function project_dir(computer :: Union{Symbol, Computer} = :current)
    return joinpath(base_dir(computer),  "CollegeStrat");
end

# Files for the paper live here.
# On local machine: write directly to Dropbox.
# On remote: these files go into an arbitrary local dir.
function paper_dir(computer :: Union{Symbol, Computer} = :current)
	if is_remote(computer)
		pDir = joinpath(project_dir(computer), "paper");
	else
		pDir = joinpath(dropbox_dir(computer), "lutz", "paper");
	end
	return pDir
end

"""
	$(SIGNATURES)

Base directory for all output files.
"""
function computer_out_dir(computer :: Union{Symbol, Computer} = :current)
    return joinpath(base_dir(computer), "out")
end

"""
	$(SIGNATURES)

Comparisons across cases are stored here.
"""
function global_comparison_dir(computer :: Union{Symbol, Computer} = :current)
    return joinpath(computer_out_dir(computer), "compare")
end

"""
	$(SIGNATURES)

Base directory for all `mat` files (binaries).
"""
function computer_mat_dir(computer :: Union{Symbol, Computer} = :current)
    return joinpath(base_dir(computer), "mat")
end

computer_json_dir(computer :: Union{Symbol, Computer} = :current) =
	joinpath(base_dir(computer), "mat");

"""
	$(SIGNATURES)

Logs need to be outside of the directory that gets `rsync`ed to longleaf.
Otherwise logs of running jobs get overwritten.
"""
function computer_log_dir(computer :: Union{Symbol, Computer} = :current)
	return joinpath(base_dir(computer),  "log")
end

# Put files generated by tests here
function test_dir(computer = :current)
	return joinpath(computer_out_dir(computer), "test_files")
end

## -----------  Files needed for paper
# Live in Dropbox (if local)

dropbox_dir(computer = :current) = "/Users/lutz/Dropbox/Dropout Policies";

"""
	$(SIGNATURES)

Directory where notation files (e.g. notation preamble) are stored.
Local dir, so it works on server.
"""
function notation_dir(computer = :current)
	# if is_remote(computer)
		pDir = joinpath(project_dir(computer), "notation");
	# else
	# 	pDir = joinpath(dropbox_dir(computer), "lutz", "notation");
	# end
	return pDir
end

notation_copy_dir(computer = :current) =
	joinpath(dropbox_dir(computer), "lutz", "notation");

"""
	$(SIGNATURES)

Copy a file, given by a relative path, from `srcDir` to `tgDir`.
Make `tgDir` if needed.
"""
function copy_file(fPath :: AbstractString, srcDir :: AbstractString, tgDir)
    srcPath = joinpath(srcDir, fPath);
    tgPath = joinpath(tgDir, fPath);
    if isfile(srcPath)
        make_dir(tgPath);
        cp(srcPath, tgPath; force = true);
    else
        @warn "Source file not found: $fPath"
    end
end

# -------


text:
using ITensors
using ITensorMPOCompression
using Revise
using Test
using Printf

import ITensorMPOCompression.truncate!
import ITensorMPOCompression.truncate
import ITensorMPOCompression.orthogonalize!

#brute force method to control the default float display format.
# Base.show(io::IO, f::Float64) = @printf(io, "%1.5f", f)

#
#  We need consistent output from randomMPS in order to avoid flakey unit testset
#  So we just pick any seed so that randomMPS (hopefull) gives us the same
#  pseudo-random output for each test run.
#
using Random
Random.Random.seed!(12345);

function test_truncate(makeH,N::Int64,NNN::Int64,ms::matrix_state,epsSVD::Float64,epsrr::Float64,eps::Float64)
    mlr=mirror(ms.lr)
    hx=0.5

    sites = siteinds("SpinHalf", N)
    psi=randomMPS(sites)
    #
    # Make right canonical, then compress to left canonical
    #
    H=makeH(sites,NNN,hx,ms.ul)
    E0l=inner(psi',H,psi)
    @test is_regular_form(H,ms.ul,eps)
    #@show get_Dw(H)
    orthogonalize!(H;orth=ms.lr,epsrr=epsrr)
    #@show get_Dw(H)

    E1l=inner(psi',H,psi)
    RE=abs((E0l-E1l)/E0l)
    #@printf "E0=%1.5f E1=%1.5f rel. error=%.5e  \n" E0l E1l RE
    @test RE<2*eps
    @test is_regular_form(H,ms.ul,eps)
    @test is_canonical(H,ms,eps)

    truncate!(H;orth=mlr,cutoff=epsSVD)
    #@show get_Dw(H)
    ss=truncate!(H;orth=ms.lr,cutoff=epsSVD)
    #@show get_Dw(H)
    @test is_regular_form(H,ms.ul,eps)
    @test is_canonical(H,ms,eps)
    #@show ITensorMPOCompression.min(ss) ITensorMPOCompression.max(ss)
    # make sure the energy in unchanged
    E2l=inner(psi',H,psi)
    RE=abs((E0l-E2l)/E0l)
    @printf "E0=%8.5f Etrunc=%8.5f rel. error=%7.1e RE/espSVD=%6.2f \n" E0l E2l RE RE/epsSVD
    if epsSVD<=1e-12
        @test (RE/epsSVD)<1.0 #typically this will remove any meaning SVs
    else
        @test (RE/epsSVD)<1000.0 #now we may remove real SVs and expect much bigger energy errors
    end
end

@testset "Compress full MPO" begin
    eps=2e-13
    epsSVD=1e-12
    epsrr=1e-12
    ll=matrix_state(lower,left)
    ul=matrix_state(upper,left)
    lr=matrix_state(lower,right)
    ur=matrix_state(upper,right)

    #                                 V=N sites
    #                                   V=Num Nearest Neighbours in H
    test_truncate(make_transIsing_MPO,5,1,ll,epsSVD,epsrr,eps)
    test_truncate(make_transIsing_MPO,5,3,ll,epsSVD,epsrr,eps)
    test_truncate(make_transIsing_MPO,5,3,ul,epsSVD,epsrr,eps)
    test_truncate(make_transIsing_MPO,5,3,lr,epsSVD,epsrr,eps)
    test_truncate(make_transIsing_MPO,5,3,ur,epsSVD,epsrr,eps)
    #test_truncate(make_transIsing_MPO,15,10,ll,epsSVD,epsrr,eps) #know fail, triggers fixing RL_prime
    test_truncate(make_transIsing_MPO,15,10,lr,epsSVD,epsrr,eps)

    epsSVD=.00001
    test_truncate(make_transIsing_MPO,10,7,ll,epsSVD,epsrr,eps)
    test_truncate(make_transIsing_MPO,10,7,lr,epsSVD,epsrr,eps)
    test_truncate(make_transIsing_MPO,10,7,ur,epsSVD,epsrr,eps)
    test_truncate(make_transIsing_MPO,10,7,ul,epsSVD,epsrr,eps)

    #
    # Heisenberg from AutoMPO
    #
    epsSVD=.00001
    test_truncate(make_Heisenberg_AutoMPO,10,7,lr,epsSVD,epsrr,eps)

end

@testset "Test ground states" begin
    eps=3e-13
    epsSVD=1e-12
    epsrr=1e-12
    N=10
    NNN=5
    hx=0.5
    db=ITensors.using_debug_checks()

    sites = siteinds("SpinHalf", N)
    H=make_transIsing_MPO(sites,NNN,hx,lower)

    ITensors.ITensors.disable_debug_checks()
    E0,psi0=fast_GS(H,sites)
    if db  ITensors.ITensors.enable_debug_checks() end
    truncate!(H;orth=left,cutoff=epsSVD,epsrr=epsrr)

    ITensors.ITensors.disable_debug_checks()
    E1,psi1=fast_GS(H,sites)
    if db  ITensors.ITensors.enable_debug_checks() end

    overlap=abs(inner(psi0',psi1))
    RE=abs((E0-E1)/E0)
    @printf "Trans. Ising E0/N=%1.15f E1/N=%1.15f rel. error=%.1e overlap-1.0=%.1e \n" E0/(N-1) E1/(N-1) RE overlap-1.0
    @test abs(E0-E1)<eps
    @test abs(overlap-1.0)<eps

    hx=0.0
    H=make_Heisenberg_AutoMPO(sites,NNN,hx,lower)

    ITensors.ITensors.disable_debug_checks()
    E0,psi0=fast_GS(H,sites)
    if db  ITensors.ITensors.enable_debug_checks() end

    truncate!(H;orth=right,cutoff=epsSVD,epsrr=epsrr)

    ITensors.ITensors.disable_debug_checks()
    E1,psi1=fast_GS(H,sites)
    if db  ITensors.ITensors.enable_debug_checks() end

    overlap=abs(inner(psi0',psi1))
    RE=abs((E0-E1)/E0)
    @printf "Heisenberg   E0/N=%1.15f E1/N=%1.15f rel. error=%.1e overlap-1.0=%.1e \n" E0/(N-1) E1/(N-1) RE overlap-1.0
    @test abs(E0-E1)<eps
    @test abs(overlap-1.0)<eps
end

@testset "Look at bond singular values for large lattices" begin
    NNN=5
    hx=0.5

    @printf "                 max(sv)           min(sv)\n"
    @printf "  N    Dw  left     mid    right   mid\n"
    for N in [6,8,10,12,18,20,40,80]
        sites = siteinds("SpinHalf", N)
        H=make_transIsing_MPO(sites,NNN,hx,lower)
        specs=truncate!(H,cutoff=1e-10)
        imid::Int64=N/2-1
        max_1  =specs[1   ].spectrum[1]
        max_mid=specs[imid].spectrum[1]
        max_N  =specs[N-1 ].spectrum[1]
        min_mid=specs[imid].spectrum[end]
        Dw=maximum(get_Dw(H))
        @printf "%4i %4i %1.5f %1.5f %1.5f %1.5f\n" N Dw  max_1 max_mid max_N min_mid
        @test (max_1-max_N)<1e-10
        @test max_mid<1.0
    end
end
#=
@testset "Test with conserved QNs" begin
    N = 10
    NNN = 1
    hx=0.0 #can't make and sx op with QNs in play
    sites = siteinds("S=1/2",N;conserve_qns=true)
    # to build our own MPO we need to put QNs on all the link indices.
    # from autoMPO we see QN("Sz",0) => 3 as the qn for a link index.
    H=make_transIsing_MPO(sites,NNN,hx,upper)
    #H=make_transIsing_AutoMPO(sites,NNN,hx,lower)
    pprint(H[2],1e-14)
    orthogonalize!(H)
    #pprint(H[2],1e-14)
    #@show H[2]
    #i=Index(QN("Sz",0)=>3;dir=ITensors.In,tags="Link,l=1")
    #@show i
end
 =#


text:
function fibsq04()
#
#  test polynomial suggested by Goedecker
#
    n = 4;
    p = fib(n)
    p = conv(p,p)
    z = [-.7748041132154339;
         -.07637893113374573-.8147036471703865*im; -.07637893113374573+.8147036471703865*im;
         1.927561975482925]
    z = [z  2*ones(n)]
    p, PolyZeros(z)
end


text:
@testset "Horn" begin
    p, q = synthetic_gradient(SynthSphere(), radius = 50)
    Z = convert_gradient(Horn(), p, q)
    @test Z[38,38] ≈ -0.0027751467456787315
    @test Z[1,1] ≈ 0.0
    @test Z[76,76] ≈ 47.05959259743248
end


text:
using SimulatedAnnealing
using Statistics
using Test

@testset "Travelling salesman" begin
    include("../example/travelling_salesman.jl")

    energies = energy.(samples)
    decrement_rules = Dict(
        "constant" => ConstantDecrementRule,
        "AVL" => AVLDecrementRule)
    stop_criteria = Dict(
        "OVG" => () -> OVGCriterion(energy_reference=mean(energies)),
        "SSV" => SSVCriterion)

    for (rulename, rule) in decrement_rules
        for (critname, criterion) in stop_criteria
            @testset "Decrement: $rulename, Stop: $critname" begin
                sa_prob = AnnealingOptimization(
                    energy,
                    propose_candidate,
                    criterion(),
                    rule(),
                    std(energies),
                    samples[1],
                    n*(n - 1))

                best_tour, tour_length = simulated_annealing(sa_prob)
                @test tour_length == 12
            end
        end
    end
end


text:
using PowerSystems
using NLsolve
const PSY = PowerSystems

############### Data Network ########################
include(joinpath(dirname(@__FILE__), "dynamic_test_data.jl"))
include(joinpath(dirname(@__FILE__), "data_utils.jl"))
############### Data Network ########################
threebus_file_dir = joinpath(dirname(@__FILE__), "ThreeBusMulti.raw")
threebus_sys = System(threebus_file_dir, runchecks = false)

function inv_case78(static_device)
    return DynamicInverter(
        name = get_name(static_device),
        ω_ref = 1.0, # ω_ref,
        converter = converter_high_power(), #converter
        outer_control = outer_control(), #outer control
        inner_control = inner_control(), #inner control voltage source
        dc_source = dc_source_lv(), #dc source
        freq_estimator = pll(), #pll
        filter = filt(), #filter
    )
end

function dyn_gen_multi_tg(generator)
    return PSY.DynamicGenerator(
        name = get_name(generator),
        ω_ref = 1.0, # ω_ref,
        machine = machine_classic(), #machine
        shaft = shaft_damping(), #shaft
        avr = avr_none(), #avr
        prime_mover = tg_type2(), #tg
        pss = pss_none(),
    )
end

# Add dynamic generators to the system (each gen is linked through a static one)
for g in get_components(Generator, threebus_sys)
    if get_number(get_bus(g)) == 101
        case_gen = inv_case78(g)
        add_component!(threebus_sys, case_gen, g)
    elseif get_number(get_bus(g)) == 102
        case_gen = dyn_gen_multi_tg(g)
        add_component!(threebus_sys, case_gen, g)
    end
end

#Create Ybus_Fault
fault_branches = deepcopy(collect(get_components(Branch, threebus_sys)))
for br in fault_branches
    if get_name(br) == "1"
        br.r = 4 * br.r
        br.x = 4 * br.x
    end
end
Ybus_fault = PSY.Ybus(fault_branches, get_components(Bus, threebus_sys))[:, :];


text:
module LearningMetaheuristics

end

max_stars_repo_name: cocomoff/LearningMetaheuristics

max_issues_repo_name: cocomoff/LearningMetaheuristics

max_forks_repo_name: cocomoff/LearningMetaheuristics


text:
using LearningHorse.LossFunction

@testset "LossFunction" begin
    sy, st = rand(), rand()
    y, t = rand(20), rand(20)
    my1, my2, my3 = fill!(Array{Float64}(undef, 1, 1), sy), rand(20, 1), rand(20, 2)

    LOSSES = [
        mse, cee, mae, huber, logcosh_loss, poisson, hinge, smooth_hinge
    ]

    for f in LOSSES
        @test typeof(f(sy, st)) <: Number
        @test typeof(f(y, t)) == Float64
        @test typeof(f(y, t, reduction="sum")) == Float64
        @test size(f(y, t, reduction="none")) == (20, )
        @test_throws ArgumentError f(y, t, reduction="NaN")
        @test f(my1, st) == f(sy, st)
        @test f(my2, st) == f(vec(my2), st)
        @test_throws DimensionMismatch f(my3, st)
    end
end


text:
# Unit tests for JUDI linear operators (without PDE solves)
# Philipp Witte (pwitte.slim@gmail.com)
# May 2018
#

using JUDI.TimeModeling, SegyIO, Test, LinearAlgebra

# Example structures

example_info(; n=(120,100), nsrc=2, ntComp=1000) = Info(prod(n), nsrc, ntComp)
example_model(; n=(120,100), d=(10f0, 10f0), o=(0f0, 0f0), m=randn(Float32, n)) = Model(n, d, o, m)

function example_rec_geometry(; nsrc=2, nrec=120)
    xrec = range(50f0, stop=1150f0, length=nrec)
    yrec = 0f0
    zrec = range(50f0, stop=50f0, length=nrec)
    return Geometry(xrec, yrec, zrec; dt=4f0, t=1000f0, nsrc=nsrc)
end

function example_src_geometry(; nsrc=2)
    xrec = range(100f0, stop=1000f0, length=nsrc)
    yrec = 0f0
    zrec = range(50f0, stop=50f0, length=nsrc)
    return Geometry(xrec, yrec, zrec; dt=4f0, t=1000f0, nsrc=nsrc)
end

# Tests

function test_transpose(Op)
    @test isequal(size(Op), size(conj(Op)))
    @test isequal(reverse(size(Op)), size(transpose(Op)))
    @test isequal(reverse(size(Op)), size(adjoint(Op)))
    @test isequal(reverse(size(Op)), size(transpose(Op)))
    return true
end

function test_getindex(Op)
    # requires: Op.info.nsrc == 2
    Op_sub = Op[1]
    @test isequal(Op_sub.info.nsrc, 1)
    @test isequal(Op_sub.model, Op.model)
    @test isequal(size(Op_sub), convert(Tuple{Int64, Int64}, size(Op) ./ 2))

    Op_sub = Op[1:2]
    @test isequal(Op_sub.info.nsrc, 2)
    @test isequal(Op_sub.model, Op.model)
    @test isequal(size(Op_sub), size(Op))
    return true
end

########################################################### judiModeling ###############################################

@testset "judiModeling Unit Test" begin

    info = example_info()
    model = example_model()
    F_forward = judiModeling(info, model; options=Options())
    F_adjoint = judiModelingAdjoint(info, model; options=Options())

    @test isequal(typeof(F_forward), judiModeling{Float32, Float32})
    @test isequal(typeof(F_adjoint), judiModelingAdjoint{Float32, Float32})

    # conj, transpose, adjoint
    @test test_transpose(F_forward)
    @test test_transpose(F_adjoint)

    # get index
    @test test_getindex(F_forward)
    @test test_getindex(F_adjoint)

end

############################################################# judiPDE ##################################################

@testset "judiPDE Unit Test" begin

    info = example_info()
    model = example_model()
    rec_geometry = example_rec_geometry()

    PDE_forward = judiPDE("PDE", info, model, rec_geometry; options=Options())
    PDE_adjoint = judiPDEadjoint("PDEadjoint", info, model, rec_geometry; options=Options())

    @test isequal(typeof(PDE_forward), judiPDE{Float32, Float32})
    @test isequal(typeof(PDE_adjoint), judiPDEadjoint{Float32, Float32})

    # conj, transpose, adjoint
    @test test_transpose(PDE_forward)
    @test test_transpose(PDE_adjoint)

    # get index
    @test test_getindex(PDE_forward)
    @test test_getindex(PDE_adjoint)

    # Multiplication w/ judiProjection
    src_geometry = example_src_geometry()
    Ps = judiProjection(info, src_geometry)

    PDE = PDE_forward*transpose(Ps)
    @test isequal(typeof(PDE), judiPDEfull{Float32, Float32})
    @test isequal(PDE.recGeometry, rec_geometry)
    @test isequal(PDE.srcGeometry, src_geometry)

    PDEad = PDE_adjoint*transpose(Ps)
    @test isequal(typeof(PDEad), judiPDEfull{Float32, Float32})
    @test isequal(PDEad.srcGeometry, rec_geometry)
    @test isequal(PDEad.recGeometry, src_geometry)

end

######################################################### judiPDEfull ##################################################

@testset "judiPDEfull Unit Test" begin

    info = example_info()
    model = example_model()
    rec_geometry = example_rec_geometry()
    src_geometry = example_src_geometry()
    PDE = judiModeling(info, model, src_geometry, rec_geometry; options=Options())

    @test isequal(typeof(PDE), judiPDEfull{Float32, Float32})
    @test isequal(PDE.recGeometry, rec_geometry)
    @test isequal(PDE.srcGeometry, src_geometry)

    @test test_transpose(PDE)
    @test test_getindex(PDE)

end

######################################################## judiJacobian ##################################################

@testset "judiJacobian Unit Test" begin

    info = example_info()
    model = example_model()
    rec_geometry = example_rec_geometry()
    src_geometry = example_src_geometry()
    wavelet = randn(Float32, src_geometry.nt[1])
    PDE = judiModeling(info, model, src_geometry, rec_geometry; options=Options())
    q = judiVector(src_geometry, wavelet)

    J = judiJacobian(PDE, q)

    @test isequal(typeof(J), judiJacobian{Float32, Float32})
    @test isequal(J.recGeometry, rec_geometry)
    @test isequal(J.srcGeometry, src_geometry)
    @test isequal(size(J)[2], prod(model.n))
    @test test_transpose(J)

    # get index
    J_sub = J[1]
    @test isequal(J_sub.info.nsrc, 1)
    @test isequal(J_sub.model, J.model)
    @test isequal(size(J_sub)[1], Int(size(J)[1]/2))

    J_sub = J[1:2]
    @test isequal(J_sub.info.nsrc, 2)
    @test isequal(J_sub.model, J.model)
    @test isequal(size(J_sub), size(J))

end

####################################################### judiProjection #################################################

@testset "judiProjection Unit Test" begin

    info = example_info()
    rec_geometry = example_rec_geometry()
    src_geometry = example_src_geometry()
    wavelet = randn(Float32, src_geometry.nt[1])
    q = judiVector(src_geometry, wavelet)

    Pr = judiProjection(info, rec_geometry)
    Ps = judiProjection(info, src_geometry)

    @test isequal(typeof(Pr), judiProjection{Float32, Float32})
    @test isequal(typeof(Ps), judiProjection{Float32, Float32})
    @test isequal(Pr.geometry, rec_geometry)
    @test isequal(Ps.geometry, src_geometry)

    @test test_transpose(Pr)
    @test test_transpose(Ps)

    Pr_sub = Pr[1]
    @test isequal(Pr_sub.info.nsrc, 1)
    @test isequal(size(Pr_sub), convert(Tuple{Int64, Int64}, size(Pr) ./ 2))

    Pr_sub = Pr[1:2]
    @test isequal(Pr_sub.info.nsrc, 2)
    @test isequal(size(Pr_sub), size(Pr))

    RHS = transpose(Ps)*q
    @test isequal(typeof(RHS), judiRHS{Float32})
    @test isequal(RHS.geometry, q.geometry)

end