text:
using EzXML
using HTTP: request

"""
    readosm(filename)

`readosm` has only one argument `filename`, taking a string of the pbf-file path and name.
It returns an object containing the OSM data.
"""
function readosm(filename::String)::Map
    return readxmldoc(EzXML.readxml(open(filename, "r")))
end

"""
    queryoverpass(bbox)

`queryoverpass` has only one argument `bbox`.
It returns an object containing the OSM data.
"""
function queryoverpass(bbox::BBox; kwargs...)::Map
    osmdata = queryoverpass("$(bbox.bottom_lat),$(bbox.left_lon),$(bbox.top_lat),$(bbox.right_lon)", kwargs...)
    return osmdata
end

"""
    queryoverpass(lonlat, radius)

`queryoverpass` has two arguments, `lonlat` and `radius`.
It returns an object containing the OSM data.
"""
function queryoverpass(lonlat::LatLon, radius::Real; kwargs...)
    osmdata = queryoverpass("around:$radius,$(lonlat.lat),$(lonlatl.lon)", kwargs...)
    return osmdata
end

"""
    queryoverpass(bounds)

`queryoverpass` has only one argument `bounds`.
It returns an object containing the OSM data.
"""
function queryoverpass(bounds::String; timeout::Int64=25)::Map
    query = """
        [out:xml][timeout:$timeout];
        (
            node($bounds);
            way($bounds);
            relation($bounds);
        );
        out body;
        >;
        out skel qt;
    """
    result = request(
        "GET",
        "https://overpass-api.de/api/interpreter",
        query=Dict("data" => query)
    )
    return readxmldoc(EzXML.readxml(IOBuffer(result.body)))
end

"""
The argument is a xml document, it returns an osm object.
"""
function readxmldoc(xmldoc::EzXML.Document)::Map
    osmdata = Map()
    # Iterate over child elements.
    for xmlnode in EzXML.eachelement(root(xmldoc))
        elname = EzXML.nodename(xmlnode)
        if elname == "bounds"
            osmdata.meta[:bbox] = osmbound(xmlnode)
        elseif elname == "node"
            k, v = osmnode(xmlnode)
            osmdata.nodes[k] = v
        elseif elname == "way"
            k, v = osmway(xmlnode)
            osmdata.ways[k] = v
        elseif elname == "relation"
            k, v = osmrelation(xmlnode)
            osmdata.relations[k] = v
        else
            merge!(osmdata.meta, osmunknown(xmlnode))
        end
    end
    return osmdata
end

function osmbound(xmlnode::EzXML.Node)::BBox
    bbox = BBox(
        parse(Float64, xmlnode["minlat"]),
        parse(Float64, xmlnode["minlon"]),
        parse(Float64, xmlnode["maxlat"]),
        parse(Float64, xmlnode["maxlon"])
    )
    return bbox
end

function osmnode(xmlnode::EzXML.Node)::Tuple{Int64,Node}
    id = parse(Int64, xmlnode["id"])
    latlon = LatLon(
        parse(Float64, xmlnode["lat"]),
        parse(Float64, xmlnode["lon"])
    )
    tags = nothing
    if EzXML.haselement(xmlnode)
        # Iterate over child elements.
        for subxmlnode in EzXML.eachelement(xmlnode)
            elname = EzXML.nodename(subxmlnode)
            if elname == "tag"
                if tags === nothing
                    tags = Dict{Symbol,String}()
                end
                tags[Symbol(subxmlnode["k"])] = subxmlnode["v"]
            end
        end
    end
    return id, Node(latlon, tags)
end

function osmway(xmlnode::EzXML.Node)::Tuple{Int64,Way}
    id = parse(Int64, xmlnode["id"])
    refs = []
    tags = nothing
    if EzXML.haselement(xmlnode)
        # Iterate over child elements.
        for subxmlnode in EzXML.eachelement(xmlnode)
            elname = EzXML.nodename(subxmlnode)
            if elname == "nd"
                push!(refs, parse(Int64, subxmlnode["ref"]))
            elseif elname == "tag"
                if tags === nothing
                    tags = Dict{Symbol,String}()
                end
                tags[Symbol(subxmlnode["k"])] = subxmlnode["v"]
            end
        end
    end
    return id, Way(refs, tags)
end

function osmrelation(xmlnode::EzXML.Node)::Tuple{Int64,Relation}
    id = parse(Int64, xmlnode["id"])
    refs = Int64[]
    types = Symbol[]
    roles = String[]
    tags = nothing
    if EzXML.haselement(xmlnode)
        # Iterate over child elements.
        for subxmlnode in EzXML.eachelement(xmlnode)
            elname = EzXML.nodename(subxmlnode)
            if elname == "member"
                push!(refs, parse(Int64, subxmlnode["ref"]))
                push!(types, Symbol(subxmlnode["type"]))
                push!(roles, subxmlnode["role"])
            elseif elname == "tag"
                if tags === nothing
                    tags = Dict{Symbol,String}()
                end
                tags[Symbol(subxmlnode["k"])] = subxmlnode["v"]
            end
        end
    end
    return id, Relation(refs, types, roles, tags)
end

function osmunknown(xmlnode::EzXML.Node)::Dict{Symbol,Any}
    out = Dict{Symbol,Any}()
    for kv in EzXML.attributes(xmlnode)
        out[Symbol(EzXML.nodename(kv))] = EzXML.nodecontent(kv)
    end
    if EzXML.haselement(xmlnode)
        for subxmlnode in EzXML.eachelement(xmlnode)
            elname = EzXML.nodename(subxmlnode)
            out[Symbol(elname)] = osmunknown(subxmlnode)
        end
    end
    return out
end


text:
module CompilerUtils

using Core.Compiler:
    CodeInfo,
    NativeInterpreter,
    OptimizationParams,
    OptimizationState,
    compact!,
    convert_to_ircode,
    copy,
    copy_exprargs,
    slot2reg

fallback_insert_new_nodes(ir) = compact!(copy(ir), true)

const insert_new_nodes = try
    Core.Compiler.insert_new_nodes
catch
    fallback_insert_new_nodes
end

function is_intermediate_value_type(t)
    @nospecialize t
    return t isa Union{Core.Compiler.DelayedTyp,Core.Compiler.NotFound}
end

function ircode_from_codeinfo(ci::CodeInfo)
    # Making a copy here, as, e.g., `convert_to_ircode` mutates `ci`:
    ci = copy(ci)

    for (i, t) in pairs(ci.ssavaluetypes)
        if is_intermediate_value_type(t)
            ci.ssavaluetypes[i] = Any
        end
    end

    linfo = ci.parent  # MethodInstance
    interp = NativeInterpreter()
    params = OptimizationParams(interp)
    opt = OptimizationState(linfo, copy(ci), params, interp)
    nargs = Int(opt.nargs) - 1
    preserve_coverage = false
    code = copy_exprargs(ci.code)
    ir = convert_to_ircode(ci, code, preserve_coverage, nargs, opt)
    ir = slot2reg(ir, ci, nargs, opt)
    ir = insert_new_nodes(ir)
    return ir
end

end  # module

ext: jl


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

p1 = plot(x=[1,2,3], y=[4,5,6])
p2 = plot(x=[1,2,3], y=[6,7,8])
p3 = plot(x=[5,7,8], y=[8,9,10])
p4 = plot(x=[5,7,8], y=[10,11,12])
gridstack([p1 p2; p3 p4])


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using BenchmarkTools
using CAGD

function generateCubeGrids(n=1, m=1, p=1)
    V,(VV,EV,FV,CV) = Lar.cuboidGrid([n,m,p],true)
    mybox = (V,CV,FV,EV)

    twocubs = Lar.Struct([mybox, Lar.t(.3,.4,.5), Lar.r(pi/5,0,0), Lar.r(0,0,pi/12), mybox])

    V,CV,FV,EV = Lar.struct2lar(twocubs)

    cop_EV = convert(Lar.ChainOp, Lar.coboundary_0(EV::Lar.Cells))
    cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
    model = CAGD.Model(V)
    CAGD.addModelCells!(model, 1, cop_EV)
    CAGD.addModelCells!(model, 2, cop_FE)

    return model
end

m = generateCubeGrids(1, 1, 1)

@btime CAGD.spatial_arrangement(m)


text:
# This script should be called from runtests.jl as
#   ./simplesum.jl 12 13 25

@testset "Basic" begin
    ans = parse(Int64, ARGS[1]) + parse(Int64, ARGS[2])
    @test ans == parse(Int64, ARGS[3])
end


text:
export xor, log_loss, log_fitness, fitness_xor, xor_dataset

function xor(a::Float64, b::Float64)
    if a + b == 1
        1.
    else
        0.
    end
end

"Array to array XOR"
function xor(a::Array{Float64})
    v = a[1]
    for i in 2:length(a)
        v = xor(v, a[i])
    end
    [v]
end

function rand_bin(len::Int64)
    1.0 * rand(0:1, len)
end

function xor_dataset(len::Int64, n_records::Int64)
    X::Array{Array{Float64}}=[]
    y::Array{Array{Float64}}=[]
    for i in 1:n_records
        l = rand_bin(len)
        push!(X, l)
        push!(y, xor(l))
    end
    X, y
end

function log_loss(y_true::Float64, y_pred::Float64)
    y_pred = maximum([minimum([y_pred, 1-10^-15]), 10^-15])
    if y_true == 1.0
        -log(y_pred)
    else
        -log(1-y_pred)
    end
end

function log_loss(y_true::Int64, y_pred::Int64)
    log_loss(float(y_true), float(y_pred))
end

function log_loss(y_true::Array, y_pred::Array)
    log_loss(y_true[1], y_pred[1])
end

function log_fitness(y_true, y_pred)
    max_loss = log_loss(y_true, 1 .- y_true)
    max_loss - log_loss(y_true, y_pred)
end

function fitness_xor(indiv::NEATIndiv, len::Int64=2)
    X, y = xor_dataset(len, 100)
    y_pred = process(indiv, X)
    [sum(log_fitness.(y, y_pred)) / length(X)]
end


text:
module PerformanceRNN

using Flux
using MIDI
using NoteSequences
using NoteSequences.PerformanceRepr

include("layers.jl")
include("generate.jl")
include("pretrain.jl")

end


text:
"""
    energy_balance(object::T,meteo::Atmosphere,constants = Constants())
    energy_balance(object::Array{AbstractComponentModel},meteo::Atmosphere,constants = Constants())
    energy_balance!(object::AbstractComponentModel,meteo::Atmosphere,constants = Constants())

Computes the energy balance of a component based on the type of the model it was parameterized
with in `object.energy`.

At the moment, two models are implemented in the package:

- [`Monteith`](@ref): the model found in Monteith and Unsworth (2013)
- `Missing`: if no computation of the energy balance is needed

# Arguments

- `object::Union{AbstractComponentModel,Array{AbstractComponentModel},
Dict{String,AbstractComponentModel}}`: a [`Component`](@ref) struct, or a Dict/Array of.
- `meteo::Union{Atmosphere,Weather}`: meteorology structure, see [`Atmosphere`](@ref) or
[`Weather`](@ref)
- `constants = Constants()`: physical constants. See [`Constants`](@ref) for more details

# Note

Some models need input values for some variables. For example [`Monteith`](@ref) requires a
value for `Rₛ`, `d` and `skyFraction`. If you read the models from a file, you can
use [`init_status!`](@ref) (see examples).

# Examples

```julia
meteo = Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65)

# Using the model of Monteith and Unsworth (2013) for energy, Farquhar et al. (1980) for
# photosynthesis, and Medlyn et al. (2011) for stomatal conductance:
leaf = LeafModels(energy = Monteith(),
            photosynthesis = Fvcb(),
            stomatal_conductance = Medlyn(0.03, 12.0),
            Rₛ = 13.747, skyFraction = 1.0, PPFD = 1500.0, d = 0.03)

energy_balance(leaf,meteo)

# Using a model file:
model = read_model("a-model-file.yml")

# An example model file is available here:
# "https://raw.githubusercontent.com/VEZY/PlantBiophysics/main/test/inputs/models/plant_coffee.yml"

# Initialising the mandatory variables:
init_status!(model, Rₛ = 13.747, skyFraction = 1.0, PPFD = 1500.0, Tₗ = 25.0, d = 0.03)

# NB: To know which variables has to be initialised according to the models used, you can use
# `to_initialise(leaf)`

# Running a simulation for all component types in the same scene:
energy_balance!(model, meteo)

model["Leaf"].status.Rn
model["Leaf"].status.A
model["Leaf"].status.Cᵢ
```

# References

Duursma, R. A., et B. E. Medlyn. 2012. « MAESPA: a model to study interactions between water
limitation, environmental drivers and vegetation function at tree and stand levels, with an
example application to [CO2] × drought interactions ». Geoscientific Model Development 5 (4):
919‑40. https://doi.org/10.5194/gmd-5-919-2012.

Monteith, John L., et Mike H. Unsworth. 2013. « Chapter 13 - Steady-State Heat Balance: (i)
Water Surfaces, Soil, and Vegetation ». In Principles of Environmental Physics (Fourth Edition),
edited by John L. Monteith et Mike H. Unsworth, 217‑47. Boston: Academic Press.

Schymanski, Stanislaus J., et Dani Or. 2017. « Leaf-Scale Experiments Reveal an Important
Omission in the Penman–Monteith Equation ». Hydrology and Earth System Sciences 21 (2): 685‑706.
https://doi.org/10.5194/hess-21-685-2017.

Vezy, Rémi, Mathias Christina, Olivier Roupsard, Yann Nouvellon, Remko Duursma, Belinda Medlyn,
Maxime Soma, et al. 2018. « Measuring and modelling energy partitioning in canopies of varying
complexity using MAESPA model ». Agricultural and Forest Meteorology 253‑254 (printemps): 203‑17.
https://doi.org/10.1016/j.agrformet.2018.02.005.
"""
function energy_balance!(object::AbstractComponentModel, meteo::Atmosphere, constants = Constants())
    is_init = is_initialised(object)
    !is_init && error("Some variables must be initialized before simulation")

    net_radiation!(object, meteo, constants)
    return nothing
end

# Same as above but non-mutating
function energy_balance(object::AbstractComponentModel, meteo::Atmosphere, constants = Constants())
    object_tmp = copy(object)
    energy_balance!(object_tmp, meteo, constants)
    return object_tmp.status
end

# energy_balance over several objects (e.g. all leaves of a plant) in an Array
function energy_balance!(object::O, meteo::Atmosphere, constants = Constants()) where O <: AbstractArray{<:AbstractComponentModel}

    for i in values(object)
        energy_balance!(i, meteo, constants)
    end

    return nothing
end

# same as the above but non-mutating
function energy_balance(
    object::O,
    meteo::M,
    constants = Constants()
    ) where {O <: Union{AbstractArray{<:AbstractComponentModel},AbstractDict{N,<:AbstractComponentModel} where N},M <: Union{Atmosphere,Weather}}

    # Copy the objects only once before the computation for performance reasons:
    object_tmp = copy(object)

    # Computation:
    energy_balance!(object_tmp, meteo, constants)

    # --- Extracting the outputs: ---

    # Pre-allocating the outputs:
    output = [i.status for i in object_tmp]

    for (i, obj) in enumerate(object_tmp)
        output[i] = obj.status
    end

    output = DataFrame([NamedTuple(i) for i in output])

    return output
end

# energy_balance over several objects (e.g. all leaves of a plant) in a kind of Dict.
function energy_balance!(object::O, meteo::Atmosphere, constants = Constants()) where O <: AbstractDict{N,<:AbstractComponentModel} where N

    for (k, v) in object
        energy_balance!(v, meteo, constants)
    end

    return nothing
end

# same as the above but non-mutating. In this case we add a column with the component name 😃
function energy_balance(
    object::O,
    meteo::M,
    constants = Constants()
    ) where {O <: Union{AbstractArray{<:AbstractComponentModel},AbstractDict{N,<:AbstractComponentModel} where N},M <: Union{Atmosphere,Weather}}

    # Copy the objects only once before the computation for performance reasons:
    object_tmp = copy(object)

    # Computation:
    energy_balance!(object_tmp, meteo, constants)

    # --- Extracting the outputs: ---

    # Pre-allocating the outputs:
    output = Dict([k => v.status for (k, v) in object_tmp])

    for (k, v) in object_tmp
        output[k] = v.status
    end

    output = DataFrame([(NamedTuple(v)..., component = k) for (k, v) in output])

    return output
end

# energy_balance over several meteo time steps (called Weather) and possibly several components.
# Only allowed for components given as a subtype of AbstractDict to track components names in the
# outputs
function energy_balance!(
    object::T,
    meteo::Weather,
    constants = Constants()
    ) where T <: Union{AbstractDict{N,<:AbstractComponentModel} where N}

    # Pre-allocating the time-step outputs:
    timestep_tmp = Dict([k => v.status for (k, v) in object])

    # Pre-allocating the general DataFrame with the first time-step results:
    energy_balance!(object, meteo.data[1], constants)

    for (k, v) in object
        timestep_tmp[k] = v.status
    end

    output_timestep = DataFrame([(NamedTuple(v)..., component = k) for (k, v) in timestep_tmp])

    # Actually pre-allocating the DF:
    output = repeat(output_timestep, length(meteo.data))
    output.time_step = repeat(1:length(meteo.data), inner = size(output_timestep, 1))

    # Computing for all following time-steps:
    for (i, meteo_i) in enumerate(meteo.data[2:end])
        energy_balance!(object, meteo_i, constants)

        for (k, v) in object
            timestep_tmp[k] = v.status
        end

        output_timestep = DataFrame([(NamedTuple(v)..., component = k) for (k, v) in timestep_tmp])

        # Update the values of the global output:
        output[output.time_step .== i,Not(:time_step)] = output_timestep
    end

    return output
end

# energy_balance over several meteo time steps (same as above) but non-mutating
function energy_balance(
    object::T,
    meteo::Weather,
    constants = Constants()
    ) where T <: Union{AbstractDict{N,<:AbstractComponentModel} where N}

    object_tmp = copy(object)

    return energy_balance!(object_tmp, meteo, constants)
end


text:
module Phione

export s2tw

include("words.jl")

const EOFChar = Char(0xA1) # End Of File Char

function destruct!(chn::Channel{Char}, str::AbstractString)
    temp = iterate(str)
    temp ≡ nothing && error("destruct!: string collection must be nonempty")

    char, idx = temp
    put!(chn, char)
    while true
        temp = iterate(str, idx)
        temp ≡ nothing && (put!(chn, EOFChar), break)
        char, idx = temp
        put!(chn, char)
    end
    return nothing
end

char2phr!(chn::Channel{Char}, char::Char) = haskey(WORDS, char) ? put!(chn, WORDS[char]) : put!(chn, char)

function char2phr!(cnChar::Channel{Char}, twChar::Channel{Char})
    while true
        letter = take!(cnChar)
        letter ≡ EOFChar && (put!(twChar, EOFChar), break)
        isdone = false
        if haskey(PHRASELENGTH, letter)
            maxlen = PHRASELENGTH[letter]
            phrase = string(letter)
            counts = 1
            while counts < maxlen && fetch(cnChar) ≠ EOFChar
                phrase *= take!(cnChar)
                counts += 1
                if haskey(PHRASES, phrase)
                    phrase = PHRASES[phrase]
                    isdone = true
                    break
                end
            end
            fun! = isdone ? put! : char2phr!
            temp = iterate(phrase)
            char, idx = temp
            fun!(twChar, char)
            while true
                temp = iterate(phrase, idx)
                temp ≡ nothing && break
                char, idx = temp
                fun!(twChar, char)
            end
        else
            char2phr!(twChar, letter)
        end
    end
    close(cnChar)
    return nothing
end

function assemble!(chn::Channel{Char})
    ret = ""
    while true
        temp = take!(chn)
        temp ≡ EOFChar && (close(chn), break)
        ret *= temp
    end
    return ret
end

function s2tw(str::AbstractString)
    cnChar = Channel{Char}(16)
    twChar = Channel{Char}(16)
    @async destruct!(cnChar, str)
    @async char2phr!(cnChar, twChar)
    return fetch(@async assemble!(twChar))::String
end

end # module


text:
using Test
using OMEinsum
using CUDA
using DoubleFloats

CUDA.allowscalar(false)

@testset "cuda einsum" begin
    a = [randn(fill(3, i)...) for i=1:4]
    ca = a .|> CuArray
    ixs = ((1,2), (2,3))
    xs = (ca[2], ca[2])
    @test loop_einsum!(EinCode(ixs, (1,3)), xs, zeros(3,3) |> CuArray, OMEinsum.get_size_dict(ixs, xs)) ≈ ca[2]*ca[2]
    for f in [ein"ij,jk->ik", ein"ii->", ein"ijj ->i", ein"ij,ik,il->jkl", ein"ii->i", ein"ijl->i", ein"i->ii", ein"ij,jk,kl->il", ein"ij,ij,ij -> ij"]
        cins = map(ix->ca[length(ix)], OMEinsum.getixs(f))
        ins = map(ix->a[length(ix)], OMEinsum.getixs(f))
        @test f(cins...) isa DenseCuArray
        @test loop_einsum(f, cins, OMEinsum.get_size_dict(OMEinsum.getixs(f), cins)) |> Array ≈ f(ins...)
        @test f(ins...) ≈ Array(f(cins...))
    end
end

@testset "fallback - getindex IR error" begin
    a = rand(ComplexF64,2,2,2)
    ca = CuArray(a);
    @test Array(ein"npu,por,dom,lmn -> urdl"(ca,ca,ca,ca)) ≈ ein"npu,por,dom,lmn -> urdl"(a,a,a,a)
end

@testset "isbatchmul" begin
    for (ixs, iy) in [(((1,2), (2,3)), (1,3)), (((1,2,3), (2,3)), (1,3)),
                        (((7,1,2,3), (2,4,3,7)), (1,4,3)),
                        (((3,), (3,)), (3,)), (((3,1), (3,)), (3,1))
                        ]
        xs = ([randn(ComplexF64, fill(4,length(ix))...) |> CuArray for ix in ixs]...,)
        @test EinCode(ixs, iy)(xs...) |> Array ≈ loop_einsum(EinCode(ixs, iy), xs, OMEinsum.get_size_dict(ixs, xs)) |> Array
    end
end

@testset "doublefloats" begin
    D = 2
    T = CuArray(rand(Double64, D, D, D, D, D, D))
    U = CuArray(rand(Double64, D, D, D))

    code = ein"abewcd,bfixgh,ajeycd,jfizgh->wxyz"
    xs = (T,T,T,T)
    M = code(xs...)
    @test M |> Array ≈ loop_einsum(code, xs, OMEinsum.get_size_dict(OMEinsum.getixs(code), xs)) |> Array

    code = ein"(ubcdef,fjz),dhx,(bvghij,eiy),cgw->uvwxyz"
    _code = ein"ubcdef,fjz,dhx,bvghij,eiy,cgw->uvwxyz"
    xs = (T,U,U,T,U,U)
    M = code(xs...)
    # mapreducedim! calls to dynamic tuple splatting.
    @test_broken M |> Array ≈ loop_einsum(_code, xs, OMEinsum.get_size_dict(OMEinsum.getixs(_code), xs)) |> Array
end

@testset "binary einsum" begin
    for code in [
            ein"ij->",  # sum
            ein"ij->j", # sum
            ein"ii->",  # tr
            ein"ii->i",  # diag
            ein"ijk->kij",   # permutedims
            ein"a->aaaa",   # ~diag
            ein"ac->acc",   # ~diag
            ein"->ii",   # ~tr
            ein"i->ik",   # ~sum
            ein"->ik",   # ~sum
            ein"illljkk->kijjcc",   # general
            # binary
            ein",->",
            ein"ijl,jl->il",
            ein"ab,bc->ac",
            ein"abb,bc->ac",  # with diag in
            ein"ab,bc->acc",  # with diag out
            ein"ab,bce->ac",  # with sum in
            ein"ab,bc->ace",  # with sum out
            ein"bal,bcl->lcae",  # with perm in
            ein"bal,bcl->ca",  # with multi-edge in
            ein"bal,bc->lca",  # with multi-edge out
            ein"ddebal,bcf->lcac",  # with all
        ]
        @info code
        D = 2
        xs = [length(ix)==0 ? CUDA.fill(1.2) : CUDA.rand(Float64, fill(D, length(ix))...) for ix in OMEinsum.getixs(code)]
        size_dict = Dict(zip(('a', 'b', 'c', 'd', 'e', 'f','i','j','k','l'), ntuple(x->D, 10)))
        res = einsum(code, (xs...,), size_dict)
        @test Array(res) ≈ loop_einsum(code, (Array.(xs)...,), size_dict)
        @test Array(res) ≈ Array(loop_einsum(code, (xs...,), size_dict))
    end
end

@testset "binary rules" begin
    for (code, a, b) in [
        (ein"j,j->", randn(10), randn(10)),
        (ein"i,->i", randn(10), fill(2.0, ())),
        (ein",->", fill(2.0,()), fill(2.0, ())),
        (ein"il,kl->ikl", randn(10, 10), randn(10, 10)),
        ]
        res0 = code(a, b)
        res1 = code(CuArray(a), CuArray(b))
        @test res1 isa CuArray
        @test res0 ≈ Array(res1)
    end
end


text:
using Transliterate
using Test

@testset "transliterate" begin
    @test transliterate("абвгдђежзијклљмнњопрстћуфхцчџш") == "abvgddjezhzijklljmnnjoprstcufhczchdzsh"
    @test transliterate("Привет") == "Privet"
    @test transliterate("Foo ÿ") == "Foo y"
	@test transliterate("Hællæ, hva skjera?") == "Haellae, hva skjera?"
    @test transliterate("Déjà Vu!") == "Deja Vu!"

    @test transliterate("ث س و"; languages=["ar"]) == "th s w"
    @test transliterate("ћ ћ ћ ჰ"; languages=["sr", "ka"]) == "c c c h"

    custom_replacements = Dict("≠" => "not equal", "∉" => "not in")
    custom_replacements_copy = deepcopy(custom_replacements)
    @test transliterate("≠ ∉"; custom_replacements) == "not equal not in"
    @test custom_replacements == custom_replacements_copy
end


text:
# 'cstrt_hansen_ngon' provides the vertices coordinates (a,b) of a convex
# equilateral small n-gon for n = 2^s and s >= 3
function cstrt_hansen_ngon(n)
if mod(log2(n),1)==0 && n>=8
    # initialization
    F(t) = (cos(t)^2+sin((n/2-2)*t)^2-2*cos(t)*sin((n/2-2)*t)^2)/(2*cos(t)^2*(1+sin((n/2-2)*t)))-4*sin(t/2)^2;
    z = find_zero(F,(2*asin(sin(pi/n)/2),pi/n));
    # construction
    a = zeros(n-1); b = zeros(n-1);
    for i = 1:Int(n/4-1)
        a[Int(mod(i*(n/2-1),n))] = (-1)^(i-1)*sin(2*i*z)/2/cos(z);
        a[Int(n-mod(i*(n/2-1),n))] = -a[Int(mod(i*(n/2-1),n))];
        b[Int(mod(i*(n/2-1),n))] = (1+(-1)^(i-1)*cos(2*i*z))/2/cos(z);
        b[Int(n-mod(i*(n/2-1),n))] = b[Int(mod(i*(n/2-1),n))];
    end
    a[Int(mod((n/4)*(n/2-1),n)+1)] = a[Int(mod((n/4-1)*(n/2-1),n))]+(-1)^(n/4-1)*sin((n/2-2)*z);
    a[Int(n-mod((n/4)*(n/2-1),n)-1)] = -a[Int(mod((n/4)*(n/2-1),n)+1)];
    b[Int(mod((n/4)*(n/2-1),n)+1)] = b[Int(mod((n/4-1)*(n/2-1),n))]+(-1)^(n/4-1)*cos((n/2-2)*z);
    b[Int(n-mod((n/4)*(n/2-1),n)-1)] = b[Int(mod((n/4)*(n/2-1),n)+1)];
    a[Int(n/4)] = 1/2; a[Int(3*n/4)] = -a[Int(n/4)];
    b[Int(n/4)] = (2*cos((n/4-1)*z)-cos(z)*(cos((n/4-1)*z)-sin((n/4-1)*z)))/
    (2*cos(z)*(cos((n/4-1)*z)+sin((n/4-1)*z)));
    b[Int(3*n/4)] = b[Int(n/4)];
    if n>8
        for i = 1:Int(n/4-2)
            a[Int(mod(i*(n/2-1)+n/2,n))] = a[Int(mod(i*(n/2-1),n))] + (-1)^i*sin(2*i*z);
            a[Int(n-mod(i*(n/2-1)+n/2,n))] = -a[Int(mod(i*(n/2-1)+n/2,n))];
            b[Int(mod(i*(n/2-1)+n/2,n))] = b[Int(mod(i*(n/2-1),n))] + (-1)^i*cos(2*i*z);
            b[Int(n-mod(i*(n/2-1)+n/2,n))] = b[Int(mod(i*(n/2-1)+n/2,n))];
        end
    end
    a[Int(n/2)] = 0; b[Int(n/2)] = 1;
    return(a,b)
end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
module LavoisierCore

using ForwardDiff
using DiffResults
#using Sobol
using MappedArrays
using LinearAlgebra
using Roots
#using BenchmarkTools
using Unitful
#using Distributions
using PositiveFactorizations

include("utils.jl")
include("core.jl")
include("solver_core.jl")
include("IAPWS95.jl")
include("gerg2008.jl")
include("solver5.jl")

export IAPWS95,GERG2008
export AbstractHelmholtzModel
export _helmholtzn,    core_helmholtz, core_gibbs,    core_grad_vt,    core_grad_x,
    core_fg_x,    core_fg_vt,    core_dfdv,    core_dfdt,    core_hessian_vt,
    core_fgh_vt,    core_pressure,    compressibility_factor,    core_entropy,
    core_enthalpy,    core_internal_energy,    core_isochoric_heat_capacity,
    core_isobaric_heat_capacity,    core_sound_speed,
    pressure,     entropy,    enthalpy,    internal_energy,    isochoric_heat_capacity,
    isobaric_heat_capacity,    sound_speed
export HelmholtzPhase, helmholtz_phase, volume, moles, mass, temperature
export mol_fraction, mol_number, mass_fraction, mass_number
export pt_flash, mass_density, mol_density, mass_volume, mol_volume
export core_mass_density, core_mol_density, core_mass_volume, core_mol_volume

export Gernert
end

#end

#degreesoffredom(materialAbstractMaterial)
#material.state

#function eqsolve!(materialAbstractMaterial;kwargs)
#degreesofFreedom(material) == 0 && error("Non-zero degrees of freedom.")
#end

##function
#function equilibrium(thermoAbstractThermoModel;stateAbstractState;kwargs)
#end

##moistAir

#end # module

#idea

# data = addcompounds(["CH4","O2","H2","CO2"])
#model = PengRobinson(data,alphaTSoave(),VdWMixtureRule())
#model = PengRobinson(data), if you dont want parameters, it shall have a sane default

#what about transport variables?, those are dependent of the temperature, pressure, etc,
#maybe using models
#but there are not dependent of the model
#X1 = [0.7,0.2,0.05,0.05]
#state1 = pvtxstate(T=400u"K",P=30u"bar",X1)

#material1 = newmaterial(state1,model)

#eqsolve!(material1)
#material2 = eqsolve(material1)
#entalphy(material1) #save cache of results in the form of the state vector and the value
#entropy(material2)  #probably using CAPE open definitions here

#another example, moist air, forking for psycro.jl to use his functions

# data = addcompounds(["H2O","N2","O2"])
#MoistAirModel = ASHRAEMoistAir(data), #the system should check if the compounds are correct, it not, then throw an error
#this can be done in the database, for specific models
#state2 = moiststate(T=273u"K",wetBulb = 291.15u"K")
#moistair1 = newmaterial(model)


text:
abstract type AbstractCFDateTime <: Dates.TimeType
end

const RegTime = Union{Dates.Millisecond,Dates.Second,Dates.Minute,Dates.Hour,Dates.Day}

for (CFDateTime,calendar) in [(:DateTimeStandard,"standard"),
                              (:DateTimeJulian,"julian"),
                              (:DateTimeProlepticGregorian,"prolepticgregorian"),
                              (:DateTimeAllLeap,"allleap"),
                              (:DateTimeNoLeap,"noleap"),
                              (:DateTime360Day,"360day")]
    @eval begin
        struct $CFDateTime <: AbstractCFDateTime
            instant::UTInstant{Millisecond}
            $CFDateTime(instant::UTInstant{Millisecond}) = new(instant)
        end
    end
end


text: # intentionally empty as this is just the client part


text:
using JuMP
using ECOS
using LinearAlgebra
using Random

Random.seed!(1234);

@testset "Test SecondOrderCone formulations" begin
    u0 = rand(10)
    p = rand(10)
    q = rand();

    optimal_value_soc = 0.7413476010760025
    @testset "Primal formulation" begin
        model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
        @variable(model, u[1:10])
        @variable(model, t)
        @objective(model, Min, t)
        @constraint(model, [t, (u - u0)...] in SecondOrderCone())
        @constraint(model, u' * p == q)
        optimize!(model)
        @test JuMP.termination_status(model) == MOI.OPTIMAL
        @test objective_value(model) ≈ optimal_value_soc atol=1e-6
    end
    @testset "Dual formulation" begin
        e1 = [1, zeros(10)...]
        dual_model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
        @variable(dual_model, y1 <= 0)
        @variable(dual_model, y2[1:11])
        @objective(dual_model, Max, q * y1 + dot(vcat(0, u0), y2))
        @constraint(dual_model, e1 - [0, p...] .* y1 - y2 .== 0)
        @constraint(dual_model, y2 in SecondOrderCone())
        optimize!(dual_model)
        @test JuMP.termination_status(dual_model) == MOI.OPTIMAL
        @test objective_value(dual_model) ≈ optimal_value_soc atol=1e-6
    end
end

@testset "Test RotatedSecondOrderCone formulation" begin
    u0 = rand(10)
    p = rand(10)
    q = rand();
    model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
    @variable(model, u[1:10])
    @variable(model, t)
    @objective(model, Min, t)
    @constraint(model, [t, 0.5, (u - u0)...] in RotatedSecondOrderCone())
    @constraint(model, u' * p == q)
    optimize!(model)
    @test JuMP.termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) ≈ 0.5495962595148015 atol=1e-6
end

@testset "Test ExponentialCone formulations" begin
    n = 15;
    m = 10;
    A = randn(m, n);
    b = rand(m, 1);

    model = Model(optimizer_with_attributes(ECOS.Optimizer, "printlevel" => 0))
    @variable(model, t[1:n])
    @variable(model, x[1:n])
    @objective(model, Max, sum(t))
    @constraint(model, sum(x) == 1)
    @constraint(model, A * x .<= b )
    # Cannot use the exponential cone directly in JuMP, hence we use MOI to specify the set.
    @constraint(model, con[i = 1:n], [t[i], x[i], 1.0] in MOI.ExponentialCone())

    optimize!(model);
    @test objective_value(model) ≈ 2.708050207925179
    @test JuMP.termination_status(model) == MOI.OPTIMAL
end


text:
using TensorNetworkCodes
using Test

@testset "five_qubit_code" begin
    code = five_qubit_code()
    @test verify(code)
    @test num_qubits(code) == 5
    @test Int(length(code.logicals) / 2) == 1
    @test find_distance_logicals(code)[1] == 3
end

@testset "five_qubit_surface_code" begin
    code = five_qubit_surface_code()
    @test verify(code)
    @test num_qubits(code) == 5
    @test Int(length(code.logicals) / 2) == 1
    @test find_distance_logicals(code)[1] == 2
end

@testset "steane_code" begin
    code = steane_code()
    @test verify(code)
    @test num_qubits(code) == 7
    @test Int(length(code.logicals) / 2) == 1
    @test find_distance_logicals(code)[1] == 3
end

@testset "random_code" begin
    code = random_code(9, 2)
    @test verify(code)
    @test num_qubits(code) == 9
    @test Int(length(code.logicals) / 2) == 2
    # invalid k > n
    @test_throws ErrorException random_code(5, 6)
end

@testset "random_stabilizer_state" begin
    code = random_stabilizer_state(6)
    @test verify(code)
    @test num_qubits(code) == 6
    @test Int(length(code.logicals) / 2) == 0
end


text:
module UtilsTest
using AlphaVantage
using HTTP
using Test

@testset "Utils" begin
    resp = HTTP.Messages.Response(
        200,
        [
            Pair("Content-Type", "application/json")
        ];
        body="""{"Note": "API limit exceeded"}"""
    )
    @test_throws Exception AlphaVantage._check_api_limit(resp)
end

end


text:
using Base.Test
using AutoEnvs

function runtests()
    include("test_debug_envs.jl")
    include("test_ngsim_utils.jl")
    include("test_ngsim_env.jl")
    include("test_vectorized_ngsim_env.jl")
    include("test_multiagent_ngsim_env.jl")
end

@time runtests()


text:
### Common types

# type parameters:
#   T: element type
#   N: the number of dimensions
#   M: contiguous rank
#
# Note: M is the maximum dimension such that
# slices up to this dimension are contiguous.
#
# For example, given a 3D contiguous array a,
# the contiguous rank of a is 3, that of
# view(a, :, :, 1:3) is 2, and that of
# view(a, :, 1:2, 1:3) is 1, etc.
#
# This rank plays a key role in efficient
# linear indexing and type-stable subview
# calculation.
#
abstract type StridedArrayView{T,N,M} <: DenseArray{T,N} end
abstract type ArrayView{T,N,M} <: StridedArrayView{T,N,M} end
abstract type UnsafeArrayView{T,N,M} <: StridedArrayView{T,N,M} end

# a type for indicating contiguous rank (statically)
mutable struct ContRank{M} end

## auxiliary union types to simplify method definition
## (for internal use only)

const Subs = Union{Real,Colon,AbstractRange}
const SubsNC = Union{Real,AbstractRange}
const SubsRange = Union{Colon,AbstractRange}

### Common methods

iscontiguous(a::DenseArray) = false
iscontiguous(a::Array) = true

offset(a::Array) = 0

contiguousrank(a::StridedArrayView{T,N,M}) where {T,N,M} = M

contrank(a::Array{T,N}) where {T,N} = ContRank{N}
contrank(a::StridedArrayView{T,N,M}) where {T,N,M} = ContRank{M}

length(a::StridedArrayView) = a.len
size(a::StridedArrayView) = a.shp

getdim(s::NTuple{N,Int}, d::Integer) where {N} = (1 <= d <= N ? s[d] : 1)
size(a::StridedArrayView{T,N}, d::Integer) where {T,N} = getdim(size(a), d)

## pointer conversion

unsafe_convert(::Type{Ptr{T}}, a::StridedArrayView{T}) where {T} = pointer(a)

## Create similar array

similar(a::StridedArrayView, ::Type{T}, dims::Dims) where {T} = Array{T}(undef, dims)


text:
# ============================================================================ #
function make_board(str::String)
    # new line feeds cause errors in stat on windows when testing if str is a
    # file...
    if !occursin(r"[\r\n]+", str) && isfile(str)
        raw = read(str, String)
    else
        raw = str
    end

    RAW_STR[1] = raw

    lines = map(x->split(x, r"\s+"), decomment(split(strip(raw), r"\r|\r\n|\n")))

    nrow = length(lines)
    ncol = mapreduce(length, max, lines)
    nel = sum(length, lines)

    brd = init_board(ncol, nrow)

    player_location = [0, 0]
    facing = East

    for r = 1:nrow
        for c = 1:length(lines[r])
            if !isempty(lines[r][c]) && (strip(lines[r][c]) != "x")
                m = match(r"^(\w+)((?:\*?|\^\w?))$", lines[r][c])
                if m == nothing
                    warn("Invalid square definition at [$(c), $(r)]: $(lines[r][c])")
                else
                    set_square(brd, c, nrow - r + 1, m[1])
                    if m[2] == "*"
                        set_coin(brd, c, nrow - r + 1)
                    elseif !isempty(m[2]) && m[2][1] == '^'
                        player_location .= [c, nrow - r + 1]
                        if c == ncol
                            # change default if we're in the last column
                            facing = West
                        end
                        if length(m[2]) > 1
                            if m[2][2] == 'n'
                                facing = North
                            elseif m[2][2] == 'e'
                                facing = East
                            elseif m[2][2] == 's'
                                facing = South
                            elseif m[2][2] == 'w'
                                facing = West
                            end
                        end
                    end
                end
            end
        end
    end

    # player not specified, put them in the first valid square
    if sum(player_location) < 1
        player_location = [first(keys(brd.squares))...]
    end

    set_player(brd, player_location..., facing)

    return brd
end
# ============================================================================ #


text:
## Generic types and methods

"""
    PSF

Abstract type.  Specific psf types will inherit from PSF
"""
abstract type PSF end

"""
    BAMFState

Abstract type.  Specific state types will inherit from BAMFState.

States hold the current state of the RJMCMC Chain (i.e. θ) and can also be used to
define a true underlying state that generated the data.
"""
abstract type BAMFState end

"""
    BAMFData

BAMFData is an abstract type.  Specific data types will inherit from BAMFData.

The most common example is 'directdetection'
"""
abstract type BAMFData end

"""
    StateFlatBg<: BAMFState

Abstract type with children that implement CPU and GPU variants
"""
abstract type StateFlatBg <: BAMFState end

mutable struct StateFlatBg_CUDA <: StateFlatBg # this gets saved in chain
    n::Int
    x::CuArray{Float32}
    y::CuArray{Float32}
    photons::CuArray{Float32}
    bg::Float32
end

"""
    StateFlatBg_CPU(n::Int, x::Vector{Float32}, y::Vector{Float32}, photons::Vector{Float32}, bg::Float32)

State of a model that has x,y positions, total integrated intensity and a flat background offset.
"""
mutable struct StateFlatBg_CPU <: StateFlatBg # this gets saved in chain
    n::Int
    x::Vector{Float32}
    y::Vector{Float32}
    photons::Vector{Float32}
    bg::Float32
end

mutable struct StateFlatBg_Results <: StateFlatBg # this gets saved in chain
    n::Int
    x::Vector{Float32}
    y::Vector{Float32}
    photons::Vector{Float32}
    σ_x::Vector{Float32}
    σ_y::Vector{Float32}
    σ_photons::Vector{Float32}
    bg::Float32
end

StateFlatBg() = StateFlatBg_CPU(0, [0], [0], [0], 0)
StateFlatBg(n::Int) = StateFlatBg_CPU(n, Vector{Float32}(undef,n),Vector{Float32}(undef,n), Vector{Float32}(undef,n), 0)
StateFlatBg(n::Int, x::Vector{Float32},y::Vector{Float32}, photons::Vector{Float32},bg::Float32)=StateFlatBg_CPU(n, x,y,photons, bg)
StateFlatBg_CUDA(n::Int) = StateFlatBg_CPU(n, CuArray{Float32}(undef,n), CuArray{Float32}(undef,n), CuArray{Float32}(undef,n), 0)
StateFlatBg_Results(n::Int) = StateFlatBg_Results(n, Vector{Float32}(undef,n),
    Vector{Float32}(undef,n), Vector{Float32}(undef,n),Vector{Float32}(undef,n),Vector{Float32}(undef,n), Vector{Float32}(undef,n), 0)

function StateFlatBg_CUDACopy!(myempty,myfull)
    idx=threadIdx().x
    myempty[idx]=myfull[idx]
    return nothing
end
function StateFlatBg(sf::StateFlatBg_CUDA) # make a deep copy
    s = StateFlatBg_CUDA(sf.n)
    @cuda threads=s.n StateFlatBg_CUDACopy!(s.x,sf.x)
    @cuda threads=s.n StateFlatBg_CUDACopy!(s.y,sf.y)
    @cuda threads=s.n StateFlatBg_CUDACopy!(s.photons,sf.photons)
    return s
end
function StateFlatBg(sf::StateFlatBg_CPU) # make a deep copy
    s = StateFlatBg(sf.n)
    for nn=1:s.n
        s.x[nn]=sf.x[nn]
        s.y[nn]=sf.y[nn]
        s.photons[nn]=sf.photons[nn]
    end
    s.bg=sf.bg
    return s
end

function addemitter!(s::StateFlatBg,x::Float32,y::Float32,photons::Float32)
    s.n=s.n+1
    push!(s.x,x)
    push!(s.y,y)
    push!(s.photons,photons)
end

function removeemitter!(s::StateFlatBg,idx::Int)
    s.n=s.n-1
    deleteat!(s.x,idx)
    deleteat!(s.y,idx)
    deleteat!(s.photons,idx)
end

function findclosest(s::StateFlatBg,idx::Int)
    if (s.n<2)||(idx<1) return 0 end

    mindis=1f5 #big
    nn=1
    for nn=1:s.n
        if nn!=idx
            dis=(s.x[idx]-s.x[nn])^2+(s.y[idx]-s.y[nn])
            if dis<mindis
                mindis=dis
                idx=nn
            end
        end
    end
    return nn
end

function findother(s::StateFlatBg,idx::Int)
    if (s.n<2)||(idx<1) return 0 end

    nn=Int(ceil(s.n*rand()))
    while nn==idx
        nn=Int(ceil(s.n*rand()))
    end

    return nn
end

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

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

## RJPrior -----------------
mutable struct RJPrior
    sz::Int
    θ_start::Float32
    θ_step::Float32
    pdf::Vector{Float32}
end
RJPrior()=RJPrior(1,1,1,[1])

function priorrnd(rjp::RJPrior)
    r=rand()
    mysum=rjp.pdf[1]
    nn=1
    while (mysum<r)&&(nn<rjp.sz)
        nn+=1
        mysum+=rjp.pdf[nn]
    end
    return Float32((nn-1)*rjp.θ_step+rjp.θ_start)
end

""""
    priorpdf(rjp::RJPrior,θ)

calculate PDF(θ)
"""
function priorpdf(rjp::RJPrior,θ)
nn=round((θ-rjp.θ_start)/rjp.θ_step)
nn=Int(max(1,min(nn,rjp.sz)))
return rjp.pdf[nn]
end

"""
    RJStruct

Holds the data to be analyzed, the priors, the PSF model, and parameters used in the RJMCMC steps.
"""
mutable struct RJStruct # contains data and all static info for Direct Detection passed to BAMF functions
    data::BAMFData
    psf::MicroscopePSFs.PSF
    prior_photons::Distributions.UnivariateDistribution
    prior_background::Distributions.UnivariateDistribution
    bndpixels::Real
    σ_xy::Real
    σ_split::Real
    σ_photons::Real
    σ_background::Real
    modeldata::BAMFData
    testdata::BAMFData
end

# RJStruct(sz,psf,xy_std,I_std,split_std) =
#     RJStruct(sz, psf, xy_std, I_std, split_std,ArrayDD(sz),2,RJPrior(),ArrayDD(sz),ArrayDD(sz))
# RJStruct(sz,psf,xy_std,I_std,split_std,data::BAMFData) =
#     RJStruct(sz, psf, xy_std, I_std, split_std,data,2,RJPrior(),deepcopy(data),deepcopy(data))
# RJStruct(sz,psf,xy_std,I_std,split_std,data::BAMFData,bndpixels,prior_photons) =
#     RJStruct(sz, psf, xy_std, I_std, split_std,data,bndpixels,prior_photons,deepcopy(data),deepcopy(data))

"""
    genBAMFData(rjs::RJStruct)

generate an empty BAMFData stucture for data type in RJStruct
"""
function genBAMFData(rjs::RJStruct)
    return genBAMFData(rjs.data)
end


text:
function suc_model(uc_system, optimizer, initial_conditions::Dict=Dict{String, Any}(), step::Int64=1)

    m = JuMP.Model(optimizer)
    #Time Information
    time_periods = PowerSystems.get_forecasts_horizon(suc_system)
    time_periods_set = 1:time_periods
    data_first_step = PowerSystems.get_forecast_initial_times(suc_system)[step]
    renewable_forecast_sample = first(PowerSystems.get_component_forecasts(PowerSystems.RenewableDispatch, suc_system, data_first_step))
    scenario_count = renewable_forecast_sample.scenario_count
    scenario_set = 1:scenario_count
    minutes_per_period = Dates.Minute(PowerSystems.get_forecasts_resolution(suc_system))/Dates.Minute(1)

    # Thermal Generation
    thermal_generators = PowerSystems.get_components(PowerSystems.ThermalStandard, suc_system)
    gen_pwl_points = Dict(PowerSystems.get_name(g) => 1:length(g.op_cost.variable) for g in thermal_generators)
    thermal_gen_names = [PowerSystems.get_name(g) for g in thermal_generators]

    JuMP.@variable(m, cg[thermal_gen_names,scenario_set,time_periods_set])
    JuMP.@variable(m, pg[thermal_gen_names,scenario_set,time_periods_set] >= 0)
    JuMP.@variable(m, rg[thermal_gen_names,scenario_set,time_periods_set] >= 0)
    JuMP.@variable(m, ug[thermal_gen_names,time_periods_set], binary=true)
    JuMP.@variable(m, vg[thermal_gen_names,time_periods_set], binary=true)
    JuMP.@variable(m, wg[thermal_gen_names,time_periods_set], binary=true)
    JuMP.@variable(m, 0 <= lambda_lg[ g in thermal_gen_names, scenario_set, gen_pwl_points[g], time_periods_set] <= 1)

    # Renewable Generation
    renewable_forecasts = PowerSystems.get_component_forecasts(PowerSystems.RenewableDispatch, suc_system, data_first_step)
    renewable_gen_names = PowerSystems.get_forecast_component_name.(renewable_forecasts)
    JuMP.@variable(m, pw[renewable_gen_names,scenario_set, time_periods_set] >= 0)

    # Loads
    fix_load_forecasts =  PowerSystems.get_component_forecasts(PowerSystems.PowerLoad, suc_system, data_first_step)

    interruptible_load_forecasts =  PowerSystems.get_component_forecasts(PowerSystems.InterruptibleLoad, suc_system, data_first_step)
    interruptible_load_names = PowerSystems.get_forecast_component_name.(interruptible_load_forecasts)
    JuMP.@variable(m, pl[ interruptible_load_names, scenario_set, time_periods_set] >= 0)

    #Objective
    JuMP.@objective(m, Min,
    sum(
        sum(
        (PowerSystems.get_op_cost(g) |> PowerSystems.get_fixed )*ug[PowerSystems.get_name(g),t] +
        (PowerSystems.get_op_cost(g) |> PowerSystems.get_startup )*vg[PowerSystems.get_name(g),t] +
        1/scenario_count*
            sum(
                cg[PowerSystems.get_name(g),s,t] -
                5000*sum(PowerSystems.get_component(il).op_cost.variable.cost[2]*pl[PowerSystems.get_forecast_component_name(il), s, t] for il in interruptible_load_forecasts) -
                5000*sum(pw[PowerSystems.get_forecast_component_name(ren), s, t] for ren in renewable_forecasts)
            for s in scenario_set)
        for g in thermal_generators)
    for t in time_periods_set)
    )

       # Constraints for first time period that require initial conditions
       for g in thermal_generators
        name = PowerSystems.get_name(g)
        gen_ini_cond = get(initial_conditions, name, Dict{Symbol, Any}())
        power_output_t0 = get(gen_ini_cond, :power_output_t0, PowerSystems.get_activepower(g))
        unit_on_t0 = get(initial_conditions, :unit_on_t0, 1.0*(power_output_t0 > 0))
        time_down_t0 = get(initial_conditions, :time_down_t0, 999.0*(1.0 - (power_output_t0 > 0)))
        time_up_t0 = get(initial_conditions, :time_up_t0, 999.0*(power_output_t0 > 0))
        activepowerlimits = PowerSystems.get_tech(g) |> PowerSystems.get_activepowerlimits
        time_minimum = PowerSystems.get_tech(g) |> PowerSystems.get_timelimits
        ramplimits = PowerSystems.get_tech(g) |> PowerSystems.get_ramplimits

        #Commitment Constraints
        if unit_on_t0 > 0
            JuMP.@constraint(m, sum( (ug[name,t]-1) for t in 1:min(time_periods, time_minimum.up - time_up_t0) ) == 0)
        else
           JuMP.@constraint(m, sum( ug[name,t] for t in 1:min(time_periods, time_minimum.down - time_down_t0) ) == 0)
        end

        JuMP.@constraint(m, ug[name,1] - unit_on_t0 == vg[name,1] - wg[name,1])

        for s in scenario_set
            # Ramp Constraints
            JuMP.@constraint(m, pg[name, s, 1] + rg[name, s, 1] - unit_on_t0*(power_output_t0 - activepowerlimits.min) <= ramplimits.up*minutes_per_period)

            JuMP.@constraint(m, unit_on_t0*(power_output_t0 - activepowerlimits.min) - pg[name, s, 1] <= ramplimits.down*minutes_per_period)
        end

        # Shut Down Ramp constraint.
        JuMP.@constraint(m, unit_on_t0*(power_output_t0 - activepowerlimits.min) <= unit_on_t0*(activepowerlimits.max - activepowerlimits.min) - max(0, activepowerlimits.max - ramplimits.down*minutes_per_period)*wg[name,1])

    end

    for t in time_periods_set
        for s in scenario_set
            # Energy Balance Constraint
            JuMP.@constraint(m,
                sum( pg[PowerSystems.get_name(g),s,t] + g.tech.activepowerlimits.min*ug[PowerSystems.get_name(g),t] for g in thermal_generators) +
                sum( pw[PowerSystems.get_forecast_component_name(g),s,t] for g in renewable_forecasts)
                == sum(PowerSystems.get_component(load).maxactivepower*PowerSystems.get_forecast_value(load, t) for load in fix_load_forecasts) +
                sum(pl[PowerSystems.get_forecast_component_name(l),s,t] for l in interruptible_load_forecasts)
            )
        end

        # InterruptibleLoad Upper Bound
        for il in interruptible_load_forecasts, s in scenario_set
            load_value = PowerSystems.get_component(il).maxactivepower*PowerSystems.get_forecast_value(il, t)
            JuMP.set_upper_bound(pl[ PowerSystems.get_forecast_component_name(il), s, t], load_value)
        end

        for reserve in PowerSystems.get_component_forecasts(PowerSystems.StaticReserve, suc_system, data_first_step), s in scenario_set
            JuMP.@constraint(m, sum(rg[PowerSystems.get_name(g), s, t] for g in thermal_generators) >= PowerSystems.get_component(reserve).requirement*PowerSystems.get_forecast_value(reserve, t)) # (3)
        end

         for g in thermal_generators
            name = PowerSystems.get_name(g)
            power_output_t0 = PowerSystems.get_activepower(g)
            unit_on_t0 = (power_output_t0 > 0) ? 1.0 : 0.0
            activepowerlimits = PowerSystems.get_tech(g) |> PowerSystems.get_activepowerlimits
            time_minimum = PowerSystems.get_tech(g) |> PowerSystems.get_timelimits
            ramplimits = PowerSystems.get_tech(g) |> PowerSystems.get_ramplimits
            piecewise_production = PowerSystems.get_op_cost(g) |> PowerSystems.get_variable

            if t > 1
                JuMP.@constraint(m, ug[name,t] - ug[name,t-1] == vg[name,t] - wg[name,t]) # (12)
                for s in scenario_set
                    JuMP.@constraint(m, pg[name, s, t] + rg[name, s, t] - pg[name, s, t-1] <= ramplimits.up*minutes_per_period) # (19)
                    JuMP.@constraint(m, pg[name, s, t-1] - pg[name, s, t] <= ramplimits.down*minutes_per_period) # (20)
                end
            end

           if t >= time_minimum.up || t == time_periods
                JuMP.@constraint(m, sum( vg[name,t2] for t2 in (t-min(time_minimum.up,time_periods)+1):t) <= ug[name,t])  # (13)
            end

            if t >= time_minimum.down || t == time_periods
                JuMP.@constraint(m, sum( wg[name,t2] for t2 in (t-min(time_minimum.down,time_periods)+1):t) <= 1 - ug[name,t])  # (14)
            end

            #Shut down and Start up ramps are 3x faster than regular ramps.
            for s in scenario_set
                JuMP.@constraint(m, pg[name, s, t] + rg[name, s, t] <= (activepowerlimits.max - activepowerlimits.min)*ug[name,t] - max(0, (activepowerlimits.max - 3*ramplimits.up*minutes_per_period))*vg[name,t]) # (17)

                if t < time_periods
                    JuMP.@constraint(m, pg[name, s, t] + rg[name, s, t] <= (activepowerlimits.max - activepowerlimits.min)*ug[name,t]  - max(0, (activepowerlimits.max - 3*ramplimits.down*minutes_per_period))*wg[name,t+1]) # (18)
                end

                JuMP.@constraint(m, pg[name, s, t] == sum((piecewise_production[l][1] - piecewise_production[1][1])*lambda_lg[name, s, l,t] for l in gen_pwl_points[name])) # (21)
                JuMP.@constraint(m, cg[name, s, t] == sum((piecewise_production[l][2] - piecewise_production[1][2])*lambda_lg[name, s, l,t] for l in gen_pwl_points[name])) # (22)
                JuMP.@constraint(m, ug[name,t] == sum(lambda_lg[name, s, l,t] for l in gen_pwl_points[name])) # (23)

            end
        end

        for rgen in renewable_forecasts
            name = PowerSystems.get_forecast_component_name(rgen)
            ub = rgen.component.tech.rating*PowerSystems.get_forecast_value(rgen,t)
            for s in scenario_set
                JuMP.set_upper_bound(pw[name, s,  t], ub[s])
            end
        end

    end

    return m

end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Compute the lighting condition of a satellite:
#       1) Sunlight,
#       2) penumbra, or
#       3) umbra.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Longo, C. R. O., Rickman, S. L (1995). Method for the Calculation of
#       Spacecraft Umbra and Penumbra Shadow Terminator Points. NASA Technical
#       Paper 3547.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export SAT_LIGHTING_SUNLIGHT, SAT_LIGHTING_UMBRA, SAT_LIGHTING_PENUMBRA
export satellite_lighting_condition

# Constants.
const SAT_LIGHTING_SUNLIGHT = 0
const SAT_LIGHTING_PENUMBRA = 1
const SAT_LIGHTING_UMBRA = 2

"""
    satellite_lighting_condition(r_i::AbstractVector, s_i::AbstractVector)

Compute the satellite lighting condition given the Sun unitary vector `s_i` [m]
and the satellite position vector `r_i` [m].

# Returns

* `SAT_LIGHTING_SUNLIGHT`: Satellite is under sunlight.
* `SAT_LIGHTING_PENUMBRA`: Satellite is at penumbra region.
* `SAT_LIGHTING_UMBRA`: Satellite is at umbra region.

"""
function satellite_lighting_condition(r_i::AbstractVector, s_i::AbstractVector)
    # Norme of the sun position vector.
    norm_s_i = norm(s_i)

    # Umbra section [1].
    Xu = R0*norm_s_i/(Rs-R0)
    Au = asin(R0/Xu)

    # Penumbra section [1].
    Xp = R0*norm_s_i/(Rs+R0)
    Ap = asin(R0/Xp)

    # Projection of the satellite position vector on the sun direction [1].
    rs_i = dot(r_i, s_i)*s_i/(norm_s_i^2)

    # Distance of the umbral cone and the spacecraft [1].
    delta_i = r_i - rs_i

    # Location of the umbral cone terminator at the projected spacecraft
    # location [1].
    ep_i = (Xu - norm(rs_i))*tan(Au)

    # Location of the penumbral cone terminator at the projected spacecraft
    # location [1].
    Kp_i = (Xp + norm(rs_i))*tan(Ap)

    # Check if the satellite is under sunlight or umbra / penumbra.
    if dot(rs_i, s_i/norm_s_i) < 0
        if norm(delta_i) > Kp_i
            return SAT_LIGHTING_SUNLIGHT
        else
            if norm(delta_i) < ep_i
                return SAT_LIGHTING_UMBRA
            else
                return SAT_LIGHTING_PENUMBRA
            end
        end
    else
        return SAT_LIGHTING_SUNLIGHT
    end
end


text:
const Position = Tuple{Int,Int}
const FindClosestResult = Union{Position,Nothing}

"Calculate the manhattan distance between two `Position`s"
manhattan_distance(a::Position, b::Position) = abs.(a .- b) |> sum

"""
    closest_coin(goblin::Position, coins::Set{Position})::FindClosestResult

Given a `Position` representing the goblin and a list of `Position`s
representing the coins that have not yet been collected, return the coin
nearest to the goblin's current position. If more than one coin is nearest
the goblin, return `nothing`.
"""
function closest_coin(goblin::Position, coins::Set{Position})::FindClosestResult
    found_coin   = nothing      # The nearest coin found so far
    found_times  = 0            # How many times the current closest coin has been found
    min_distance = typemax(Int) # Distance to the nearest coin found so far

    for coin in coins
        distance = manhattan_distance(goblin, coin)

        # We've seen this distance before...
        if distance == min_distance
            found_times += 1
        end

        # The new closest coin
        if distance < min_distance
            found_coin = coin
            found_times = 1
            min_distance = distance
        end
    end

    return found_times > 1 ? nothing : found_coin
end

"""
    path(goblin::Position, coins::Vector{Position})

Given a `Position` representing the goblin and a list of `Position`s
representing coins, return a list of `Position`s, in order, that represent
the path taken by the goblin from nearest coin to nearest coin. Returns an
empty list if the goblin ever finds itself stuck between two equidistant coins.
"""
function path(goblin::Position, coins::Vector{Position})
    found_path    = []
    coin_set      = Set(coins)
    curr_position = goblin

    while !isempty(coin_set)
        curr_position = closest_coin(curr_position, coin_set)
        isnothing(curr_position) && return []
        push!(found_path, curr_position)
        delete!(coin_set, curr_position)
    end

    return found_path
end

#=------------------------------------------------------------------------------
| Test Cases
------------------------------------------------------------------------------=#

using Test

@testset "Should find a path when positions are negative" begin
    goblin = (1, 1)
    coins  = [(1, 2), (-5, -5), (1, 3)]
    @test path(goblin, coins) == [(1, 2), (1, 3), (-5, -5)]
end

@testset "Should find a path when positions are large numbers" begin
    goblin = (34, 51)
    coins  = [(1050, -34), (-8809080, 2), (-46788, -67683)]
    @test path(goblin, coins) == [(1050, -34), (-46788, -67683), (-8809080, 2)]
end

@testset "Should find a path when positions are all positive" begin
    goblin = (0, 2)
    coins  = [(0, 4), (1, 0), (2, 0), (3, 2), (2, 4)]
    @test path(goblin, coins) == [(0, 4), (2, 4), (3, 2), (2, 0), (1, 0)]
end

@testset "Should return an empty list when the goblin freaks out." begin
    goblin = (100, 100)
    coins  = [(99, 100), (101, 100), (100, 99), (100, 101)]
    @test path(goblin, coins) == []
end


text:
using ITensors,
      Test

include("util.jl")

@testset "MPS Basics" begin

  N = 10
  sites = [Index(2,"Site") for n=1:N]
  psi = MPS(sites)
  @test length(psi) == N
  @test length(MPS()) == 0

  str = split(sprint(show, psi), '\n')
  @test str[1] == "MPS"
  @test length(str) == length(psi) + 2

  @test siteind(psi,2) == sites[2]
  @test hasind(psi[3],linkind(psi,2))
  @test hasind(psi[3],linkind(psi,3))

  psi[1] = ITensor(sites[1])
  @test hasind(psi[1],sites[1])

  @testset "productMPS" begin
    @testset "vector of string input" begin
      sites = siteinds("S=1/2",N)
      state = fill("",N)
      for j=1:N
        state[j] = isodd(j) ? "Up" : "Dn"
      end
      psi = productMPS(sites,state)
      for j=1:N
        sign = isodd(j) ? +1.0 : -1.0
        @test (psi[j]*op(sites,"Sz",j)*dag(prime(psi[j],"Site")))[] ≈ sign/2
      end
      @test_throws DimensionMismatch productMPS(sites, fill("", N - 1))
    end

    @testset "vector of int input" begin
      sites = siteinds("S=1/2",N)
      state = fill(0,N)
      for j=1:N
        state[j] = isodd(j) ? 1 : 2
      end
      psi = productMPS(sites,state)
      for j=1:N
        sign = isodd(j) ? +1.0 : -1.0
        @test (psi[j]*op(sites,"Sz",j)*dag(prime(psi[j],"Site")))[] ≈ sign/2
      end
    end
    @testset "vector of ivals input" begin
      sites  = siteinds("S=1/2",N)
      states = fill(0,N)
      for j=1:N
        states[j] = isodd(j) ? 1 : 2
      end
      ivals  = [state(sites[n],states[n]) for n=1:length(sites)]
      psi = productMPS(ivals)
      for j=1:N
        sign = isodd(j) ? +1.0 : -1.0
        @test (psi[j]*op(sites,"Sz",j)*dag(prime(psi[j],"Site")))[] ≈ sign/2
      end
    end
  end

  @testset "randomMPS with chi==1" begin
    phi = randomMPS(sites)

    @test maxlinkdim(phi) == 1

    @test hasind(phi[1],sites[1])
    @test norm(phi[1])≈1.0

    @test hasind(phi[4],sites[4])
    @test norm(phi[4])≈1.0
  end

  @testset "inner different MPS" begin
    phi = randomMPS(sites)
    psi = randomMPS(sites)
    phipsi = dag(phi[1])*psi[1]
    for j = 2:N
      phipsi *= dag(phi[j])*psi[j]
    end
    @test phipsi[] ≈ inner(phi,psi)

    badsites = [Index(2) for n=1:N+1]
    badpsi = randomMPS(badsites)
    @test_throws DimensionMismatch inner(phi,badpsi)
  end

  @testset "inner same MPS" begin
    psi = randomMPS(sites)
    psidag = dag(psi)
    ITensors.primelinkinds!(psidag)
    psipsi = psidag[1]*psi[1]
    for j = 2:N
      psipsi *= psidag[j]*psi[j]
    end
    @test psipsi[] ≈ inner(psi,psi)
  end

  @testset "scaling MPS" begin
    psi = randomMPS(sites)
    twopsidag = 2.0*dag(psi)
    ITensors.primelinkinds!(twopsidag)
    @test inner(twopsidag, psi) ≈ 2.0*inner(psi,psi)
  end

  @testset "flip sign of MPS" begin
    psi = randomMPS(sites)
    minuspsidag = -dag(psi)
    ITensors.primelinkinds!(minuspsidag)
    @test inner(minuspsidag, psi) ≈ -inner(psi,psi)
  end

  @testset "add MPS" begin
    psi = randomMPS(sites)
    phi = deepcopy(psi)
    xi = add(psi, phi)
    @test inner(xi, xi) ≈ 4.0 * inner(psi, psi)
    # sum of many MPSs
    Ks = [randomMPS(sites) for i in 1:3]
    K12  = add(Ks[1], Ks[2])
    K123 = add(K12, Ks[3])
    @test inner(sum(Ks), K123) ≈ inner(K123,K123)
  end

  @testset "+ MPS" begin
    psi = randomMPS(sites)
    phi = deepcopy(psi)
    xi = psi + phi
    @test inner(xi, xi) ≈ 4.0 * inner(psi, psi)
    # sum of many MPSs
    Ks = [randomMPS(sites) for i in 1:3]
    K12  = Ks[1] + Ks[2]
    K123 = K12 + Ks[3]
    @test inner(sum(Ks), K123) ≈ inner(K123,K123)
  end

  sites = siteinds(2,N)
  psi = MPS(sites)
  @test length(psi) == N # just make sure this works
  @test length(siteinds(psi)) == N

  psi = randomMPS(sites)
  orthogonalize!(psi, N-1)
  @test ITensors.leftlim(psi) == N-2
  @test ITensors.rightlim(psi) == N
  orthogonalize!(psi, 2)
  @test ITensors.leftlim(psi) == 1
  @test ITensors.rightlim(psi) == 3
  psi = randomMPS(sites)
  ITensors.setrightlim!(psi, N+1) # do this to test qr
                                  # from rightmost tensor
  orthogonalize!(psi, div(N, 2))
  @test ITensors.leftlim(psi) == div(N, 2) - 1
  @test ITensors.rightlim(psi) == div(N, 2) + 1

  @test_throws ErrorException linkind(MPS(N, fill(ITensor(), N), 0, N + 1), 1)

  @testset "replacebond!" begin
  # make sure factorization preserves the bond index tags
    psi = randomMPS(sites)
    phi = psi[1]*psi[2]
    bondindtags = tags(linkind(psi,1))
    replacebond!(psi,1,phi)
    @test tags(linkind(psi,1)) == bondindtags

    # check that replacebond! updates llim and rlim properly
    orthogonalize!(psi,5)
    phi = psi[5]*psi[6]
    replacebond!(psi, 5, phi; ortho = "left")
    @test ITensors.leftlim(psi) == 5
    @test ITensors.rightlim(psi) == 7

    phi = psi[5]*psi[6]
    replacebond!(psi, 5, phi; ortho = "right")
    @test ITensors.leftlim(psi) == 4
    @test ITensors.rightlim(psi)==6

    ITensors.setleftlim!(psi, 3)
    ITensors.setrightlim!(psi, 7)
    phi = psi[5]*psi[6]
    replacebond!(psi, 5, phi; ortho = "left")
    @test ITensors.leftlim(psi) == 3
    @test ITensors.rightlim(psi) == 7
  end

end

@testset "orthogonalize! with QNs" begin
  N = 8
  sites = siteinds("S=1/2",N, conserve_qns=true)
  init_state = [isodd(n) ? "Up" : "Dn" for n=1:N]
  psi0 = productMPS(sites,init_state)
  orthogonalize!(psi0,4)
  @test ITensors.leftlim(psi0) == 3
  @test ITensors.rightlim(psi0) == 5
end

# Helper function for making MPS
function basicRandomMPS(N::Int;dim=4)
  sites = [Index(2,"Site") for n=1:N]
  M = MPS(sites)
  links = [Index(dim,"n=$(n-1),Link") for n=1:N+1]
  for n=1:N
    M[n] = randomITensor(links[n],sites[n],links[n+1])
  end
  M[1] *= delta(links[1])
  M[N] *= delta(links[N+1])
  M[1] /= sqrt(inner(M,M))
  return M
end

@testset "MPS gauging and truncation" begin

  N = 30

  @testset "orthogonalize! method" begin
    c = 12
    M = basicRandomMPS(N)
    orthogonalize!(M,c)

    @test ITensors.leftlim(M) == c-1
    @test ITensors.rightlim(M) == c+1

    # Test for left-orthogonality
    L = M[1]*prime(M[1],"Link")
    l = linkind(M,1)
    @test norm(L-delta(l,l')) < 1E-12
    for j=2:c-1
      L = L*M[j]*prime(M[j],"Link")
      l = linkind(M,j)
      @test norm(L-delta(l,l')) < 1E-12
    end

    # Test for right-orthogonality
    R = M[N]*prime(M[N],"Link")
    r = linkind(M,N-1)
    @test norm(R-delta(r,r')) < 1E-12
    for j in reverse(c+1:N-1)
      R = R*M[j]*prime(M[j],"Link")
      r = linkind(M,j-1)
      @test norm(R-delta(r,r')) < 1E-12
    end

    @test norm(M[c]) ≈ 1.0
  end

  @testset "truncate! method" begin
    M = basicRandomMPS(N;dim=10)
    M0 = copy(M)
    truncate!(M;maxdim=5)

    @test ITensors.rightlim(M) == 2

    # Test for right-orthogonality
    R = M[N]*prime(M[N],"Link")
    r = linkind(M,N-1)
    @test norm(R-delta(r,r')) < 1E-12
    for j in reverse(2:N-1)
      R = R*M[j]*prime(M[j],"Link")
      r = linkind(M,j-1)
      @test norm(R-delta(r,r')) < 1E-12
    end

    @test inner(M,M0) > 0.1
  end

end

@testset "Other MPS methods" begin

  @testset "sample! method" begin
    N = 10
    sites = [Index(3,"Site,n=$n") for n=1:N]
    psi = makeRandomMPS(sites,chi=3)
    nrm2 = inner(psi,psi)
    psi[1] *= (1.0/sqrt(nrm2))

    s = sample!(psi)

    @test length(s) == N
    for n=1:N
      @test 1 <= s[n] <= 3
    end

    # Throws becase not orthogonalized to site 1:
    orthogonalize!(psi,3)
    @test_throws ErrorException sample(psi)

    # Throws becase not normalized
    orthogonalize!(psi,1)
    psi[1] *= (5.0/norm(psi[1]))
    @test_throws ErrorException sample(psi)

    # Works when ortho & normalized:
    orthogonalize!(psi,1)
    psi[1] *= (1.0/norm(psi[1]))
    s = sample(psi)
    @test length(s) == N
  end

  @testset "randomMPS with chi > 1" begin
    N = 20
    chi = 8
    sites = siteinds(2,N)
    M = randomMPS(sites,chi)

    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2

    @test norm(M[1]) ≈ 1.0

    @test maxlinkdim(M) == chi

    # Test for right-orthogonality
    R = M[N]*prime(M[N],"Link")
    r = linkind(M,N-1)
    @test norm(R-delta(r,r')) < 1E-10
    for j in reverse(2:N-1)
      R = R*M[j]*prime(M[j],"Link")
      r = linkind(M,j-1)
      @test norm(R-delta(r,r')) < 1E-10
    end
  end

  @testset "randomMPS from initial state (QN case)" begin
    N = 20
    chi = 8
    sites = siteinds("S=1/2",N;conserve_qns=true)

    # Make flux-zero random MPS
    state = [isodd(n) ? 1 : 2 for n=1:N]
    M = randomMPS(sites,state,chi)
    @test flux(M) == QN("Sz",0)

    @test ITensors.leftlim(M) == 0
    @test ITensors.rightlim(M) == 2

    @test norm(M[1]) ≈ 1.0
    @test inner(M,M) ≈ 1.0

    @test maxlinkdim(M) == chi

    # Test making random MPS with different flux
    state[1] = 2
    M = randomMPS(sites,state,chi)
    @test flux(M) == QN("Sz",-2)
    state[3] = 2
    M = randomMPS(sites,state,chi)
    @test flux(M) == QN("Sz",-4)
  end

end

nothing


text:
@timedtestset "cross" begin
    @timedtestset "Correctness tests" begin
        rng, P, Q, gpc = MersenneTwister(123456), 2, 3, GPC()

        f1 = wrap(GP(sin, SEKernel()), gpc)
        f2 = wrap(GP(cos, SEKernel()), gpc)
        f3 = cross([f1, f2])
        f4 = cross([f1])
        f5 = cross([f2])

        x1 = collect(range(-5.0, 5.0; length=P))
        x2 = collect(range(-5.0, 5.0; length=Q))
        x3 = BlockData([x1, x2])
        x4 = BlockData([x1])
        x5 = BlockData([x2])

        # mean
        @test mean(f3(x3)) == vcat(mean(f1(x1)), mean(f2(x2)))
        @test mean(f4(x4)) == mean(f1(x1))
        @test mean(f5(x5)) == mean(f2(x2))

        # cov
        @test cov(f3(x3)) ≈ vcat(
            hcat(cov(f1(x1)), cov(f1(x1), f2(x2))),
            hcat(cov(f2(x2), f1(x1)), cov(f2(x2))),
        )
        @test cov(f3(x3), f3(x3)) == cov(f3(x3))
        @test cov(f4(x4), f4(x4)) == cov(f4(x4))
        @test cov(f5(x5), f5(x5)) == cov(f5(x5))
        @test cov(f4(x4)) == cov(f1(x1))
        @test cov(f5(x5)) == cov(f2(x2))

        # cross-cov (with block)
        @test cov(f1(x1), f2(x2)) == cov(f4(x4), f5(x5))
        @test cov(f2(x2), f1(x1)) == cov(f5(x5), f4(x4))
        @test cov(f3(x3), f4(x4)) == vcat(cov(f1(x1)), cov(f2(x2), f1(x1)))

        @test cov(f5(x5), f3(x3)) == hcat(cov(f2(x2), f1(x1)), cov(f2(x2)))

        # cross-cov (with non-block)
        @test cov(f4(x4)) == cov(f1(x1))
        @test cov(f5(x5)) == cov(f2(x2))
        @test cov(f3(x3), f4(x4)) == cov(f3(x3), f1(x1))
        @test cov(f5(x5), f3(x3)) == cov(f2(x2), f3(x3))
    end
    @timedtestset "Standardised Tests" begin
        rng, P, Q = MersenneTwister(123456), 3, 5
        x0_1, x0_2 = collect(range(-1.0, 1.0; length=P)), collect(range(2.0, 4.0; length=P))
        x1_1, x1_2 = randn(rng, Q), randn(rng, Q)
        x0, x1 = BlockData([x0_1, x0_2]), BlockData([x1_1, x1_2])
        x2, x3 = randn(rng, P), randn(2P)

        gpc = GPC()
        f1 = wrap(GP(sin, SEKernel()), gpc)
        f2 = wrap(GP(cos, SEKernel()), gpc)
        f3 = cross([f1, f2])
        abstractgp_interface_tests(f3, f1, x0, x1, x2, x3)

        # x1, x2 = collect(range(-2.0, 2.0; length=5)), collect(range(1.2, 1.5; length=4))
        # z1, z2 = collect(range(-1.5, 0.75; length=3)), collect(range(0.89, 2.0; length=4))
        # standard_1D_tests(
        #     MersenneTwister(123456),
        #     Dict(:l1=>0.5, :l2=>2.3),
        #     θ->begin
        #         gpc = GPC()
        #         f1 = θ[:l1] * GP(sin, SqExponentialKernel(), gpc)
        #         f2 = θ[:l2] * GP(cos, SqExponentialKernel(), gpc)
        #         f3 = cross([f1, f2])
        #         return f3, f3
        #     end,
        #     BlockData([x1, x2]),
        #     BlockData([z1, z2]),
        # )
    end
    @timedtestset "finites_to_block" begin
        rng, P, Q = MersenneTwister(123456), 3, 5
        xp, xq = randn(rng, P), randn(rng, Q)
        gpc = GPC()
        f_ = wrap(GP(sin, SEKernel()), gpc)
        g_ = wrap(GP(cos, SEKernel()), gpc)
        h_x = Stheno.finites_to_block([f_(xp, 1.0), g_(xq, 1.0)])
    end
end


text:
"""
    basemap(cmd0::String=""; kwargs...)

Plot base maps and frames.

Full option list at [`psbasemap`]($(GMTdoc)basemap.html)

Parameters
----------

- $(GMT.opt_J)
- $(GMT.opt_R)
- **A** | **polygon** :: [Type => Str | []]

    No plotting is performed. Instead, we determine the geographical coordinates of the polygon
    outline for the (possibly oblique) rectangular map domain.
    ($(GMTdoc)basemap.html#a)
- $(GMT.opt_B)
- **D** | **inset** | **inset_box** :: [Type => Str]

    Draw a simple map insert box on the map. Requires -F.
    ($(GMTdoc)basemap.html#d)
- **F** | **box** :: [Type => Str]

    Without further options, draws a rectangular border around any map insert (D), map scale (L)
    or map rose (T)
    ($(GMTdoc)basemap.html#f)
- $(GMT.opt_Jz)
- **L** | **map_scale** :: [Type => Str]

    Draw a map scale.
    ($(GMTdoc)basemap.html#l)
- $(GMT.opt_P)
- **Td** | **rose** :: [Type => Str]

    Draws a map directional rose on the map at the location defined by the reference and anchor points.
    ($(GMTdoc)basemap.html#t)
- **Tm** | **compass** :: [Type => Str]

    Draws a map magnetic rose on the map at the location defined by the reference and anchor points.
    ($(GMTdoc)basemap.html#t)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bo)
- $(GMT.opt_f)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function basemap(cmd0::String="", arg1=nothing; first=true, kwargs...)

	gmt_proggy = (IamModern[1]) ? "basemap "  : "psbasemap "
    length(kwargs) == 0 && return monolitic("psbasemap", cmd0, arg1)

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX" * split(def_fig_size, '/')[1] * "/0")
	cmd, = parse_common_opts(d, cmd, [:F :UVXY :JZ :bo :c :f :p :t :params], first)
    cmd  = parse_these_opts(cmd, d, [[:A :polygon]])
	cmd  = parse_Td(d, cmd)
	cmd  = parse_Tm(d, cmd)
	cmd  = parse_L(d, cmd)
	opt_D = parse_type_anchor(d, "", [:D :inset :inset_box],
	                        (map=("-g", arg2str, 1), outside=("J", nothing, 1), inside=("j", nothing, 1), norm=("-n", arg2str, 1), paper=("-x", arg2str, 1), anchor=("", arg2str, 2), width="+w", size="+w", justify="+j", offset=("+o", arg2str), save="+s", translate="_+t", units="_+u"), 'j')
    (!IamModern[1] && opt_D != "") && (cmd *= opt_D)
    (IamModern[1] && opt_D != "") && @warn("The `inset` option is not available in modern mode. Please use the `inset()` function.")
	_cmd, K = finish_PS_nested(d, [gmt_proggy * cmd], K)
	finish_PS_module(d, _cmd, "", K, O, true, arg1)
end

# ---------------------------------------------------------------------------------------------------
basemap!(cmd0::String="", arg1=nothing; kw...) = basemap(cmd0, arg1; first=false, kw...)

const psbasemap  = basemap 		# Alias
const psbasemap! = basemap!		# Alias


text:
using CUDAnativelib
using Base.Test

include("curand.jl")


text:
import Base: \
import LinearAlgebra: lu

using Base: OneTo
using LinearAlgebra: givensAlgorithm

"""
    rotate_right!(R, from, to, Q) → nothing

Rotates the diagonal elements between R[from:to,from:to]
such that the eigenvalue at R[to,to] → R[from,from]
and all the others move R[i,i] → R[i+1,i+1].

Assumes that `from` and `to` refer to the first index of the 2×2 block in the
case of complex conjugate pairs of eigenvalues.
"""
function rotate_right!(R, from::Integer, to::Integer, Q = NotWanted())
    i = to

    while i > from
        # Let's just do this check always; sometimes a block might split.
        curr_11 = is_start_of_11_block(R, i)
        prev_11 = is_end_of_11_block(R, i-1)
        j = prev_11 ? i - 1 : i - 2
        swap!(R, j, prev_11, curr_11, Q)
        i = j
    end

    nothing
end

"""
    rotate_left!(R, from, to, Q) → nothing

Rotates the diagonal elements between R[from:to,from:to]
such that the eigenvalue at R[from,from] → R[to,to]
and all the others move R[i,i] → R[i-1,i-1].

Assumes that `from` and `to` refer to the first index of the 2×2 block in the
case of complex conjugate pairs of eigenvalues.
"""
function rotate_left!(R, from::Integer, to::Integer, Q = NotWanted())
    i = from

    while true
        # Let's just do this check always; sometimes a block might split.
        curr_11 = is_start_of_11_block(R, i)
        j = curr_11 ? i + 1 : i + 2
        j > to && break
        next_11 = is_start_of_11_block(R, j)

        swap!(R, i, curr_11, next_11, Q)
        i = next_11 ? i + 1 : i + 2
    end

    nothing
end

struct CompletelyPivotedLU{T,N,TA<:SMatrix{N,N,T},TP}
    A::TA
    p::TP
    q::TP
    singular::Bool
end

struct CompletePivoting end

"""
    lu(A, CompletePivoting) → CompletelyPivotedLU(B, p, q)

Computes the LU factorization of A using complete pivoting.
"""
function lu(A::SMatrix{N,N,T}, ::Type{CompletePivoting}) where {N,T}
    # A, p and q should allocate, but escape analysis will eliminate this!
    if isbitstype(T)
        A = MMatrix(A)
    else
        A = SizedMatrix(A)
    end
    q = @MVector fill(N, N)
    p = @MVector fill(N, N)
    singular = false

    # Maybe I should consider doing this recursively.
    for k = OneTo(N - 1)
        # Find max value in sub-part.
        m, n, maxval = 1, 1, zero(real(T))
        for j = k:N, i = k:N
            if abs(A[i,j]) > maxval
                m, n, maxval = i, j, abs(A[i,j])
            end
        end

        # Store the p and q
        p[k] = m
        q[k] = n

        # Swap row and col
        for j = k:N
            A[k, j], A[m, j] = A[m, j], A[k, j]
        end

        for j = k:N
            A[j, k], A[j, n] = A[j, n], A[j, k]
        end

        Akk = A[k,k]

        # It has actually happened :(
        if iszero(Akk)
            singular = true
            break
        end

        for i = k+1:N
            A[i, k] /= Akk
        end

        for j = k+1:N
            Akj = A[k,j]

            for i = k+1:N
                A[i,j] -= A[i,k] * Akj
            end
        end
    end

    if iszero(A[N,N])
        singular = true
    end

    # Back to immutable land!
    return CompletelyPivotedLU(SMatrix(A), SVector(p), SVector(q), singular)
end

function (\)(LU::CompletelyPivotedLU{T,N}, b::SVector{N}) where {T,N}
    if isbitstype(T)
        x = MVector(b)
    else
        x = SizedVector(b)
    end

    # x ← L \ (P * b)
    for i = OneTo(N)
        x[i], x[LU.p[i]] = x[LU.p[i]], x[i]
        for j = i+1:N
            x[j] -= LU.A[j, i] * x[i]
        end
    end

    # x ← Q * (U \ x)
    for i = N:-1:1
        for j = N:-1:i+1
            x[i] -= LU.A[i, j] * x[j]
        end
        x[i] /= LU.A[i,i]
        x[i], x[LU.q[i]] = x[LU.q[i]], x[i]
    end

    # Go back to immutable land!
    SVector(x)
end

@inline sylvsystem(A::SMatrix{1,1,T}, B::SMatrix{2,2,T}) where {T} =
    @SMatrix [A[1,1]-B[1,1] -B[2,1]       ;
              -B[1,2]        A[1,1]-B[2,2]]

@inline sylvsystem(A::SMatrix{2,2,T}, B::SMatrix{1,1,T}) where {T} =
    @SMatrix [A[1,1]-B[1,1] A[1,2]       ;
              A[2,1]        A[2,2]-B[1,1]]

@inline sylvsystem(A::SMatrix{2,2,T}, B::SMatrix{2,2,T}) where {T} =
    @SMatrix [A[1,1]-B[1,1] A[1,2]        -B[2,1]       T(0)         ;
              A[2,1]        A[2,2]-B[1,1] T(0)          -B[2,1]      ;
              -B[1,2]       T(0)          A[1,1]-B[2,2] A[1,2]       ;
              T(0)          -B[1,2]       A[2,1]        A[2,2]-B[2,2]]

"""
    sylv(A, B, C) → X, singular

Solve A * X - X * B = C for X, where A and B are 1×1 or 2×2 matrices.

It works by recasting the Sylvester equation to a linear system
(I ⊗ A + Bᵀ ⊗ I) vec(X) = vec(C) of size 2 or 4, which is then solved by
Gaussian elimination with complete pivoting.

If the eigenvalues of A and B are equal, then `singular = true`.
"""
@inline function sylv(A::SMatrix{N,N,T}, B::SMatrix{M,M,T}, C::SMatrix{N,M,T}) where {T,N,M}
    fact = lu(sylvsystem(A, B), CompletePivoting)
    rhs = SVector{N*M,T}(C)
    SMatrix{N,M,T}(fact \ rhs), fact.singular
end

"""
    swap22_rotations(X) → c₁, s₁, c₂, s₂, c₃, s₃, c₄, s₄

Construct two sets of two Given's rotations that transform
```
  -x₁₁ -x₁₂
  -x₂₁ -x₂₂
   1    .
   .    1
```
to upper triangular form:
```
  x x  ⇒  x x  ⇒  * *  ⇒  x x  ⇒  x x
  x x  ⇒  * *  ⇒  . *  ⇒  . x  ⇒  . *
  1 .  ⇒  . *  ⇒  . x  ⇒  . *  ⇒  . .
  . 1  ⇒  . 1  ⇒  . 1  ⇒  . .  ⇒  . .
```
"""
function swap22_rotations(X::SMatrix{2,2,T}) where {T}
    # Upper triangulize first column of X
    c₁, s₁, nrm₁ = givensAlgorithm(-X[2,1], T(1))
    c₂, s₂, nrm₂ = givensAlgorithm(-X[1,1], nrm₁)

    # Apply the Givens rotations to the second column of X
    X22 = c₁ * -X[2,2] # + s₁ * 0
    X32 = -s₁' * -X[2,2] # + c₁ * 0
    # X12 = c₂ * -X[1,2] + s₂ * X22 # aint gonna need it!
    X22 = -s₂' * -X[1,2] + c₂ * X22

    # Upper triangularize the second column of X
    c₃, s₃, nrm₃ = givensAlgorithm(X32, T(1))
    c₄, s₄, nrm₄ = givensAlgorithm(X22, nrm₃)

    return c₁, s₁, c₂, s₂, c₃, s₃, c₄, s₄
end

"""
    swap12_rotations(X) → c₁, s₁, c₂, s₂

Construct two Given's rotations that transform
```
  -x₁₁ -x₁₂
   1    .
   .    1
```
to upper triangular form:
```
  x x  ⇒  * *  ⇒  x x
  1 .  ⇒  . *  ⇒  . *
  . 1  ⇒  . 1  ⇒  . .
```
"""
function swap12_rotations(X::SMatrix{1,2,T}) where {T}
    # Upper triangulize first column of X
    c₁, s₁, nrm₁₁ = givensAlgorithm(-X[1,1], T(1))

    # Apply the Givens rotations to the second column of X
    # X12 = c₁₁ * -X[1,2] # + s₁₁ * 0 # ain't gonna need it!
    X22 = -s₁' * -X[1,2] # + c₁₁ * 0

    # Upper triangularize the second column of X
    c₂, s₂, nrm₁₂ = givensAlgorithm(X22, T(1))

    return c₁, s₁, c₂, s₂
end

"""
    swap21_rotations(X) → c₁, s₁, c₂, s₂

Construct two Given's rotations that transform
```
  -x₁₁
  -x₂₁
   1
```
to upper triangular form:
```
  x  ⇒  x  ⇒  x
  x  ⇒  *  ⇒  .
  1  ⇒  .  ⇒  .
```
"""
function swap21_rotations(X::SMatrix{2,1,T}) where {T}
    c₁, s₁, nrm₁ = givensAlgorithm(-X[2,1], T(1))
    c₂, s₂, nrm₂ = givensAlgorithm(-X[1,1], nrm₁)
    return c₁, s₁, c₂, s₂
end

"""
    swap22!(R, i) → R

Swap a 2×2 block with a 2×2 block in R[i:i+4,i:i+4] via unitary transformations.
```
    i
    ↓
i → x x * *  ⇒  y y * *
    x x * *  ⇒  y y * *
    . . y y  ⇒  . . x x
    . . y y  ⇒  . . x x
```
Swapping means similarity transformation.
"""
function swap22!(R::AbstractMatrix{T}, i::Integer, Q = NotWanted()) where {T}
    m, n = size(R)

    @inbounds begin

        # Copy the upper triangular blocks + connection between them.
        A = @SMatrix [R[i+0,i+0] R[i+0,i+1]; R[i+1,i+0] R[i+1,i+1]]
        B = @SMatrix [R[i+2,i+2] R[i+2,i+3]; R[i+3,i+2] R[i+3,i+3]]
        C = @SMatrix [R[i+0,i+2] R[i+0,i+3]; R[i+1,i+2] R[i+1,i+3]]

        # A * X - X * B = C
        X, singular = sylv(A, B, C)

        # No need to swap if eigenvalues are indistinguishable
        singular && return R

        # Rotations that upper triangularize X
        c₁,s₁, c₂,s₂, c₃,s₃, c₄,s₄ = swap22_rotations(X)
        G₁ = Rotation3(c₁, s₁, c₂, s₂, i+0)
        G₂ = Rotation3(c₃, s₃, c₄, s₄, i+1)

        # Apply to R
        lmul!(G₁, R, i, n)
        rmul!(R, G₁, 1, i+3)
        lmul!(G₂, R, i, n)
        rmul!(R, G₂, 1, i+3)

        # Zero out things.
        R[i+2,i+0] = zero(T)
        R[i+3,i+0] = zero(T)
        R[i+2,i+1] = zero(T)
        R[i+3,i+1] = zero(T)

        # Accumulate
        rmul!(Q, G₁)
        rmul!(Q, G₂)
    end

    R
end

"""
    swap21!(R, i) → R

Swap a 1×1 block with a 2×2 block in R[i:i+3,i:i+3] via unitary transformations
```
    i
    ↓
i → x x *  ⇒  y * *
    x x *  ⇒  . x x
    . . y  ⇒  . x x
```
Swapping means similarity transformation.
"""
function swap21!(R::AbstractMatrix{T}, i::Integer, Q = NotWanted()) where {T}
    m, n = size(R)

    @inbounds begin

        # Copy the upper triangular blocks + connection between them.
        A = @SMatrix [R[i+0,i+0] R[i+0,i+1]; R[i+1,i+0] R[i+1,i+1]]
        B = @SMatrix [R[i+2,i+2]]
        C = @SMatrix [R[i+0,i+2]; R[i+1,i+2]]

        # A * X - X * B = C
        X, singular = sylv(A, B, C)

        # No need to swap if eigenvalues are indistinguishable
        singular && return R

        # Rotations that upper triangularize X
        c₁,s₁, c₂,s₂ = swap21_rotations(X)
        G₁ = Rotation3(c₁, s₁, c₂, s₂, i)

        # Apply rotations
        lmul!(G₁, R, i, n)
        rmul!(R, G₁, 1, i+2)

        # Zero out things.
        R[i+1,i+0] = zero(T)
        R[i+2,i+0] = zero(T)

        # Accumulate
        rmul!(Q, G₁)
    end

    R
end

"""
    swap12!(R, i) → R

Swap a 1×1 block with a 2×2 block in R[i:i+3,i:i+3] via unitary transformations
```
    i
    ↓
i → y * *  ⇒  x x *
    . x x  ⇒  x x *
    . x x  ⇒  . . y
```
Swapping means similarity transformation.

TODO: we can optimize this swap to use just one Rotation3 rather than two
Rotation2's.
"""
function swap12!(R::AbstractMatrix{T}, i::Integer, Q = NotWanted()) where {T}
    m, n = size(R)

    @inbounds begin

        # Copy the upper triangular blocks + connection between them.
        A = @SMatrix [R[i+0,i+0]]
        B = @SMatrix [R[i+1,i+1] R[i+1,i+2]; R[i+2,i+1] R[i+2,i+2]]
        C = @SMatrix [R[i+0,i+1] R[i+0,i+2]]

        # A * X - X * B = C
        X, singular = sylv(A, B, C)

        # No need to swap if eigenvalues are indistinguishable
        singular && return R

        # Rotations that upper triangularize X
        c₁,s₁, c₂,s₂ = swap12_rotations(X)
        G₁ = Rotation2(c₁, s₁, i+0)
        G₂ = Rotation2(c₂, s₂, i+1)

        # Apply rotations
        lmul!(G₁, R, i, n)
        rmul!(R, G₁, 1, i+2)
        lmul!(G₂, R, i, n)
        rmul!(R, G₂, 1, i+2)

        # Zero out things.
        R[i+2,i+0] = zero(T)
        R[i+2,i+1] = zero(T)

        # Accumulate
        rmul!(Q, G₁)
        rmul!(Q, G₂)
    end

    R
end

"""
    swap11!(R, i) → R

Swap a 1×1 block with a 1×1 block in R[i:i+1,i:i+1] via unitary transformations
    i
    ↓
i → x *  ⇒  y *
    . y  ⇒  . x
"""
function swap11!(R::AbstractMatrix, i::Integer, Q = NotWanted())
    m, n = size(R)

    @inbounds begin
        R₁₁ = R[i+0,i+0]
        R₁₂ = R[i+0,i+1]
        R₂₂ = R[i+1,i+1]

        # Turns out the Sylvester equation is not so hard to solve in this case.
        G, = get_rotation(R₁₂, R₂₂ - R₁₁, i)

        # Miniscule optimization by not touching R[i:i+1,i:i+1]
        lmul!(G, R, i+2, n)
        rmul!(R, G, 1, i-1)
        R[i+0,i+0] = R₂₂
        R[i+1,i+1] = R₁₁

        # Accumulate
        rmul!(Q, G)
    end

    R
end

"""
    swap!(R, i, curr_11, next_11, Q) → nothing

Swap the two consecutive blocks of the Schur form starting at index i.
"""
function swap!(R::AbstractMatrix, i::Integer, curr_11::Bool, next_11::Bool, Q = NotWanted())
    if curr_11
        if next_11
            swap11!(R, i, Q)
        else
            swap12!(R, i, Q)
        end
    else
        if next_11
            swap21!(R, i, Q)
        else
            swap22!(R, i, Q)
        end
    end
end

@inline is_start_of_11_block(R, i) = i == size(R, 2) || @inbounds(iszero(R[i+1,i]))
@inline is_end_of_11_block(R, i) = i == 1 || @inbounds(iszero(R[i,i-1]))


text:
## -----------  Notation
#=
The notation is defined in a `numbers` spreadsheet. Everything else is generated from this.

Flow when modifying notation:
1. Edit the `numbers` file.
2. Export to `notation_table.csv`.
3. `write_notation_preamble()`
4. `write_notation_summary()`

Typeset `model_notation.lyx`

Since the notation material is also used in the paper (written by Oksana), all of these files live in Dropbox.
=#

const GpaLabel = "AFQT";

# This holds a copy of the symbol table, so it need not be reloaded all the time.
const SymTable = SymbolTable();

"""
	$(SIGNATURES)

Returns `SymbolTable` with notation. If `SymTable` is populated, returns it. Otherwise reads a CSV file.

This will be used by `write_notation_preamble`.

Todo: Allow for recursive symbol definitions. This can already be done, but one needs to ensure that dependencies are in groups that come earlier than uses.
"""
function symbol_table(; forceReload :: Bool = false)
    if forceReload  ||  isempty(SymTable)
        reload_symbol_table();
    end
    return SymTable
end

function reload_symbol_table()
    LatexLH.erase!(SymTable);
    copy_symbol_table_from_dropbox();
    load_from_csv!(SymTable, sym_table_path());
end

sym_table_path(computer = :current) =
    joinpath(notation_dir(computer), "notation_table.csv");

# Symbol table is hand-edited in a Dropbox dir. Needs to imported periodically.
function copy_symbol_table_from_dropbox()
    if !is_remote(get_computer(:current))
        tgDir = notation_dir();
        isdir(tgDir)  ||  mkpath(tgDir);
        srcPath = joinpath(dropbox_dir(), "lutz", "notation", "notation_table.csv");
        @assert isfile(srcPath)  "Not found: $srcPath";
        cp(srcPath, sym_table_path(); force = true);
    end
end

# Latex symbol by name
lsymbol(name :: Symbol) = latex(symbol_table(), name);
ldescription(name :: Symbol) = LatexLH.description(symbol_table(), name);

"""
	$(SIGNATURES)

Make an entry for a symbol. Returns a vector with (description, newcommand name, numerical value).

# Example
```
symbol_entry(st, :workTime, 1.2) == ("Work time", "workTimeV", 1.2)
```
"""
function symbol_entry(st :: SymbolTable, sName :: Symbol, value)
    si = st[sName];
    return LatexLH.description(si), "$(si.name)V", value
end

"""
	$(SIGNATURES)

Write preamble with newcommand statements that define notation for the paper.
Writes to a local directory so that code can run on the cluster.
Calls `symbol_table` to make code revisions easier.
"""
function write_notation_preamble(; fPath = notation_preamble_path(),
    copyToDropbox = false)
    st = symbol_table(forceReload = true);
    # fPath = notation_preamble_path();
    println("Writing notation preamble to ", fpath_to_show(fPath));
    make_dir(fPath);
    open(fPath, "w") do io
        write_preamble(io, st);
    end
    if copyToDropbox  &&  !is_remote(:current)
        tgPath = joinpath(notation_copy_dir(), "notation_preamble.tex");
        cp(fPath, tgPath; force = true);
    end
    println("Done");
    return fPath
end

notation_preamble_path() = joinpath(notation_dir(), "notation_preamble.tex");

"""
	$(SIGNATURES)

Write a tex file that lists notation. To be included in another document as a Section.
Writes to a local directory so that code can run on the cluster.
"""
function write_notation_summary(; fPath = notation_summary_path(),
    copyToDropbox = false)
    println("Writing notation summary to ", fpath_to_show(fPath));
    make_dir(fPath);
    open(fPath, "w") do io
        write_notation_tex(io, symbol_table(forceReload = true));
    end
    if copyToDropbox  &&  !is_remote(:current)
        tgPath = joinpath(notation_copy_dir(), "notation_summary.tex");
        cp(fPath, tgPath; force = true);
    end
    println("Done");
end

notation_summary_path() = joinpath(notation_dir(), "notation_summary.tex");

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


text:
function eval_definition(fn, init_val, lb, ub)
    return (
        eval(Meta.parse(fn)),
        eval(Meta.parse(init_val)),
        eval(Meta.parse(lb)),
        eval(Meta.parse(ub))
    )
end

function parse_contents(contents, filename)
    content_type, content_string = split(contents, ',')
    decoded = base64decode(content_string)
    df = DataFrame()
    try
      if occursin("csv", filename)
        str = String(decoded)
        df =  CSV.read(IOBuffer(str), DataFrame)
      end
    catch e
      print(e)
      return html_div([
          "There was an error processing this file."
      ])
    end
     #shuffle data
     n_row = nrow(df)
     shuffle_data = df[shuffle(1:n_row)[:], :]
     shuffle_data[!,"id"] = 1:n_row
    return shuffle_data
  end

  function load_default(val)
    if val == "1"
      df = CSV.read(download("https://raw.githubusercontent.com/efmanu/TuringDashApp.jl/master/TuringDash/datasets/data1.csv"), DataFrame)
    elseif val =="2"
      df = CSV.read(download("https://raw.githubusercontent.com/efmanu/TuringDashApp.jl/master/TuringDash/datasets/data2.csv"), DataFrame)
    else
      df = CSV.read(download("https://raw.githubusercontent.com/efmanu/TuringDashApp.jl/master/TuringDash/datasets/data3.csv"), DataFrame)
    end
    return df
  end


text:
@params struct IpoptOptions
    nt::NamedTuple
end
function IpoptOptions(;
    first_order = true,
    hessian_approximation = first_order ? "limited-memory" : "exact",
    kwargs...,
)
    h = hessian_approximation
    return IpoptOptions((;hessian_approximation = h, kwargs...))
end

@params mutable struct IpoptWorkspace <: Workspace
    model::Model
    problem::Ipopt.IpoptProblem
    x0::AbstractVector
    options::IpoptOptions
    counter::Base.RefValue{Int}
end
function IpoptWorkspace(
    model::Model, x0::AbstractVector = getinit(model);
    options = IpoptOptions(), kwargs...,
)
    problem, counter = getipopt_problem(
        model, x0,
        options.nt.hessian_approximation == "limited-memory",
    )
    return IpoptWorkspace(model, problem, x0, options, counter)
end
@params struct IpoptResult
    minimizer
    minimum
    problem
    status
    fcalls::Int
end

function optimize!(workspace::IpoptWorkspace)
    @unpack problem, options, counter = workspace
    counter[] = 0
    foreach(keys(options.nt)) do k
        v = options.nt[k]
        Ipopt.addOption(problem, string(k), v)
    end
    solvestat = Ipopt.solveProblem(problem)
    return IpoptResult(
        copy(problem.x), problem.obj_val,
        problem, solvestat, counter[]
    )
end

struct IpoptAlg{O}
    options::O
end
IpoptAlg(; kwargs...) = IpoptAlg(kwargs)

function Workspace(model::AbstractModel, optimizer::IpoptAlg, args...; kwargs...,)
    return IpoptWorkspace(model, args...; kwargs...)
end

# Implement these for sparse matrices
function fill_indices!(rows, cols, J0::Matrix; offset = 0)
    nconstr, nvars = size(J0)
    for j in 1:nvars
        cols[offset + 1 : offset + nconstr] .= j
        rows[offset + 1 : offset + nconstr] .= 1:nconstr
        offset += nconstr
    end
    return rows, cols
end
function fill_indices!(rows, cols, HL::LowerTriangular{<:Real, <:Matrix}; offset = 0)
    nvars = size(HL, 1)
    for j in 1:nvars
        cols[offset + 1 : offset + nvars - j + 1] .= j
        rows[offset + 1 : offset + nvars - j + 1] .= j:nvars
        offset += nvars - j + 1
    end
    return rows, cols
end
function fill_indices!(rows, cols, HL::SparseMatrixCSC; offset = 0)
    for col in 1:length(HL.colptr)-1
        indices = HL.colptr[col]:HL.colptr[col+1]-1
        nvars = length(indices)
        cols[offset + 1 : offset + nvars] .= col
        rows[offset + 1 : offset + nvars] = HL.rowval[indices]
        offset += nvars
    end
    return rows, cols
end

function add_values!(values, J::Matrix; offset = 0)
    nvars = length(J)
    values[offset+1:offset+nvars] .+= vec(J)
    return values
end
function add_values!(values, HL::LowerTriangular{<:Real, <:Matrix}; factor = 1, offset = 0)
    nvars = size(HL, 1)
    for j in 1:nvars
        values[offset + 1 : offset + nvars - j + 1] .+= HL[j:nvars, j] .* factor
        offset += nvars - j + 1
    end
    return values
end
function add_values!(values, HL::SparseMatrixCSC; factor = 1, offset = 0)
    nvars = length(HL.nzval)
    values[offset+1:offset+nvars] .= HL.nzval .* factor
    return values
end

nvalues(::Nothing) = 0
nvalues(J::Matrix) = length(J)
nvalues(H::LowerTriangular{<:Real, <:Matrix}) = (size(H, 1) + 1) * size(H, 1) ÷ 2
nvalues(H::SparseMatrixCSC) = length(H.nzval)

_dot(f, x, y) = dot(f(x), y)
_dot(::Nothing, ::Any, ::Any) = 0.0

function getipopt_problem(model::Model, x0::AbstractVector, first_order::Bool)
    eq = if length(model.eq_constraints.fs) == 0
        nothing
    else
        model.eq_constraints
    end
    ineq = if length(model.ineq_constraints.fs) == 0
        nothing
    else
        model.ineq_constraints
    end
    obj = CountingFunction(getobjective(model))
    return getipopt_problem(
        obj,
        ineq,
        eq,
        x0,
        getmin(model),
        getmax(model),
        first_order,
    ), obj.counter
end
function getipopt_problem(obj, ineq_constr, eq_constr, x0, xlb, xub, first_order)
    nvars = 0
    if ineq_constr !== nothing
        ineqJ0 = Zygote.jacobian(ineq_constr, x0)[1]
        ineq_nconstr, nvars = size(ineqJ0)
        Joffset = nvalues(ineqJ0)
    else
        ineqJ0 = nothing
        ineq_nconstr = 0
        Joffset = 0
    end
    if eq_constr !== nothing
        eqJ0 = Zygote.jacobian(eq_constr, x0)[1]
        eq_nconstr, nvars = size(eqJ0)
    else
        eqJ0 = nothing
        eq_nconstr = 0
    end
    @assert nvars > 0
    lag(factor, y) = x -> begin
        factor * obj(x) +
            _dot(ineq_constr, x, @view(y[1:ineq_nconstr])) +
            _dot(eq_constr, x, @view(y[ineq_nconstr+1:end]))
    end
    clb = [fill(-Inf, ineq_nconstr); zeros(eq_nconstr)]
    cub = zeros(ineq_nconstr + eq_nconstr)

    function eval_g(x::Vector{Float64}, g::Vector{Float64})
        if ineq_constr !== nothing
            g[1:ineq_nconstr] .= ineq_constr(x)
        end
        if eq_constr !== nothing
            g[ineq_nconstr+1:end] .= eq_constr(x)
        end
        return g
    end
    function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
        grad_f .= Zygote.gradient(obj, x)[1]
    end
    function eval_jac_g(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, values::Vector{Float64})
        if mode == :Structure
            ineqJ0 === nothing || fill_indices!(rows, cols, ineqJ0)
            eqJ0 === nothing || fill_indices!(rows, cols, eqJ0, offset = Joffset)
        else
            values .= 0
            if ineq_constr !== nothing
                ineqJ = Zygote.jacobian(ineq_constr, x)[1]
                add_values!(values, ineqJ)
            end
            if eq_constr !== nothing
                eqJ = Zygote.jacobian(eq_constr, x)[1]
                add_values!(values, eqJ, offset = Joffset)
            end
        end
    end

    if first_order
        eval_h = (x...) -> 0.0
        Hnvalues = 0
    else
        HL0 = LowerTriangular(
            Zygote.hessian(
                lag(1.0, ones(ineq_nconstr + eq_nconstr)),
                x0,
            ),
        )
        eval_h = function (x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, obj_factor::Float64, lambda::Vector{Float64}, values::Vector{Float64})
            if mode == :Structure
                fill_indices!(rows, cols, HL0)
            else
                HL = LowerTriangular(
                    Zygote.hessian(lag(obj_factor, lambda), x),
                )
                values .= 0
                add_values!(values, HL)
            end
        end
        Hnvalues = nvalues(HL0)
    end
    prob = Ipopt.createProblem(
        nvars, xlb, xub, ineq_nconstr + eq_nconstr, clb, cub,
        nvalues(ineqJ0) + nvalues(eqJ0), Hnvalues, obj,
        eval_g, eval_grad_f, eval_jac_g, eval_h,
    )
    prob.x = x0
    return prob
end