text:
# Determine MATLAB library path and provide facilities to load libraries with
# this path

function matlab_home_path()
    matlab_home = get(ENV, "MATLAB_HOME", "")
    if isempty(matlab_home)
        if islinux()
            matlab_home = dirname(dirname(realpath(chomp(readstring(`which matlab`)))))
        elseif isapple()
            default_dir = "/Applications"
            if isdir(default_dir)
                dirs = readdir(default_dir)
                filter!(app -> ismatch(r"^MATLAB_R[0-9]+[ab]\.app$", app), dirs)
                if !isempty(dirs)
                    matlab_home = joinpath(default_dir, maximum(dirs))
                end
            end
        elseif iswindows()
            default_dir = Sys.WORD_SIZE == 32 ? "C:\\Program Files (x86)\\MATLAB" : "C:\\Program Files\\MATLAB"
            if isdir(default_dir)
                dirs = readdir(default_dir)
                filter!(dir -> ismatch(r"^R[0-9]+[ab]$", dir), dirs)
                if !isempty(dirs)
                    matlab_home = joinpath(default_dir, maximum(dirs))
                end
            end
        end
    end
    if isempty(matlab_home)
        error("The MATLAB path could not be found. Set the MATLAB_HOME environmental variable to specify the MATLAB path.")
    end
    return matlab_home
end

function matlab_lib_path()
    # get path to MATLAB libraries
    matlab_home = matlab_home_path()
    matlab_lib_dir = if islinux()
        Sys.WORD_SIZE == 32 ? "glnx86" : "glnxa64"
    elseif isapple()
        Sys.WORD_SIZE == 32 ? "maci" : "maci64"
    elseif iswindows()
        Sys.WORD_SIZE == 32 ? "win32" : "win64"
    end
    matlab_lib_path = joinpath(matlab_home, "bin", matlab_lib_dir)
    if !isdir(matlab_lib_path)
        error("The MATLAB library path could not be found.")
    end
    return matlab_lib_path
end

function matlab_startcmd()
    matlab_home = matlab_home_path()
    if !iswindows()
        default_startcmd = joinpath(matlab_home, "bin", "matlab")
        if !isfile(default_startcmd)
            error("The MATLAB path is invalid. Set the MATLAB_HOME evironmental variable to the MATLAB root.")
        end
        default_startcmd = "exec $(Base.shell_escape(default_startcmd))"
    elseif iswindows()
        default_startcmd = joinpath(matlab_home, "bin", (Sys.WORD_SIZE == 32 ? "win32" : "win64"), "MATLAB.exe")
        if !isfile(default_startcmd)
            error("The MATLAB path is invalid. Set the MATLAB_HOME evironmental variable to the MATLAB root.")
        end
    end
    return default_startcmd
end

# helper library access function

engfunc(fun::Symbol) = Libdl.dlsym(libeng[], fun)
mxfunc(fun::Symbol)  = Libdl.dlsym(libmx[], fun)
matfunc(fun::Symbol) = Libdl.dlsym(libmat[], fun)


text:
module ValueSymbols

using Compat
import Compat: String

export ValueSymbol
immutable ValueSymbol
    ptr::UInt
    ValueSymbol(sym::Symbol) = new(Ptr{Cchar}(Cstring(sym)))
end

function Base.convert(::Type{Symbol}, vsym::ValueSymbol)
    Symbol(unsafe_wrap(String, Ptr{Cchar}(vsym.ptr)))
end

Base.promote_rule(::Type{ValueSymbol}, ::Type{Symbol}) = ValueSymbol

# Output
Base.show(io::IO, vsym::ValueSymbol) = show(io, convert(Symbol, vsym))

# Equality is based on pointer comparison
Base. ==(vsym1::ValueSymbol, vsym2::ValueSymbol) = vsym1.ptr == vsym2.ptr
Base. ==(vsym1::ValueSymbol, sym2::Symbol) = vsym1 == ValueSymbol(sym2)
Base. ==(sym1::Symbol, vsym2::ValueSymbol) = ValueSymbol(sym1) == vsym2.ptr
function Base.isless(vsym1::ValueSymbol, vsym2::ValueSymbol)
    unsafe_wrap(String, Ptr{Cchar}(vsym1.ptr)) <
        unsafe_wrap(String, Ptr{Cchar}(vsym2.ptr))
end
Base.isless(vsym1::ValueSymbol, sym2::Symbol) = vsym1 < ValueSymbol(sym2)
Base.isless(sym1::Symbol, vsym2::ValueSymbol) = ValueSymbol(sym1) < vsym2

typealias SerType begin
    VERSION < v"0.5.0-dev+4755" ? SerializationState : AbstractSerializer
end

function Base.serialize(ser::SerType, vsym::ValueSymbol)
    Base.serialize_type(ser, ValueSymbol)
    write(ser.io, Symbol(vsym))
end

function Base.deserialize(ser::SerType, ::Type{ValueSymbol})
    ValueSymbol(read(ser.io, Symbol))
end

end

ext: jl


text:
module FieldDicts

export FieldDict

"""
    FieldDict{V}(x) <: AbstractDict{Symbol,V}

Wraps `x` and provides access to its fields through a dictionary interface. The resulting
key value pairs correspond to `x`'s field names and respective values.

# Examples

```jldoctest fielddict_docstring
julia> using FieldDicts

julia> mutable struct Foo
           x::Int
           y::Float64
       end

julia> x = Foo(1, 2);

julia> d = FieldDict(x)
FieldDict{Real, Foo} with 2 entries:
  :x => 1
  :y => 2.0

```

The keys and properties are the same as the underlying structures field names.
```jldoctest fielddict_docstring
julia> keys(d) == propertynames(d) == (:x, :y)
true

```

The values are similarly accessible through the dictionary interface.
```jldoctest fielddict_docstring
julia> collect(values(d)) == [1, 2]
true

```

These fields can be accessed via traditional dictionary-like access or the dot-property notation.
```jldoctest fielddict_docstring
julia> d[:x] = 1;

julia> d.x == d[:x] == 1
true

julia> get(d, :y, 3)
2.0

julia> get(d, :z, 3)
3

```

"""
struct FieldDict{V,P} <: AbstractDict{Symbol,V}
    parent::P

    FieldDict{V,P}(p) where {V,P} = new{V,P}(p)
    FieldDict{V}(p::P) where {V,P} = FieldDict{V,P}(p)
    FieldDict(p::P) where {P} = FieldDict{_struct_eltype(P)}(p)
end

Base.length(x::FieldDict{V,P}) where {V,P} = fieldcount(P)

_struct_eltype(::Type{T}) where {T} = __struct_eltype(T, Val(fieldcount(T)))
@generated function __struct_eltype(T::DataType, ::Val{N}) where {N}
    if N === 0
        return :Any
    elseif N === 1
        return :(fieldtype(T, 1))
    else
        out = :(fieldtype(T, 1))
        for i in 2:N
            out = :(typejoin($out, fieldtype(T, $i)))
        end
        return Expr(:block, Expr(:meta, :inline), out)
    end
end

@inline function Base.iterate(x::FieldDict{V,P}) where {V,P}
    if fieldcount(P) === 0
        return nothing
    else
        return (Pair{Symbol,V}(fieldname(P, 1), getfield(getfield(x, :parent), 1)), 2)
    end
end
@inline function Base.iterate(x::FieldDict{V,P}, state::Int) where {V,P}
    if fieldcount(P) < state
        return nothing
    else
        return (Pair{Symbol,V}(fieldname(P, state), getfield(getfield(x, :parent), state)), state + 1)
    end
end

const FieldValues{V,P} = Base.ValueIterator{FieldDict{V,P}}
@inline function Base.iterate(x::FieldValues{V,P}) where {V,P}
    if fieldcount(P) === 0
        return nothing
    else
        return (getfield(getfield(getfield(x, :dict), :parent), 1), 2)
    end
end
@inline function Base.iterate(x::FieldValues{V,P}, state::Int) where {V,P}
    if fieldcount(P) < state
        return nothing
    else
        return (getfield(getfield(getfield(x, :dict), :parent), state), state + 1)
    end
end

Base.keys(x::FieldDict{V,P}) where {V,P} = fieldnames(P)

Base.propertynames(x::FieldDict) = keys(x)

@inline Base.getproperty(x::FieldDict, s::Symbol) = getfield(getfield(x, :parent), s)

@inline Base.setproperty!(x::FieldDict, s::Symbol, v) = setfield!(getfield(x, :parent), s, v)

@inline Base.getindex(x::FieldDict, s::Symbol) = getproperty(x, s)

@inline Base.setindex!(x::FieldDict, v, s::Symbol) = setproperty!(x, s, v)

@inline function Base.get(x::FieldDict{V,P}, s::Symbol, default) where {V,P}
    p = getfield(x, 1)
    i = Base.fieldindex(P, s, false)
    if isdefined(p, i)
        return getfield(p, i)
    else
        return default
    end
end
@inline function Base.get!(x::FieldDict{V,P}, s::Symbol, default) where {V,P}
    p = getfield(x, 1)
    i = Base.fieldindex(P, s, false)
    if isdefined(p, i)
        return getfield(p, i)
    else
        setfield!(p, i, default)
        return default
    end
end
@inline function Base.get(f::Union{Function,Type}, x::FieldDict{V,P}, s::Symbol) where {V,P}
    p = getfield(x, 1)
    i = Base.fieldindex(P, s, false)
    if isdefined(p, i)
        return getfield(p, i)
    else
        return f()
    end
end
@inline function Base.get!(f::Union{Function,Type}, x::FieldDict{V,P}, s::Symbol) where {V,P}
    p = getfield(x, 1)
    i = Base.fieldindex(P, s, false)
    if isdefined(p, i)
        return getfield(p, i)
    else
        default = f()
        setfield!(p, i, default)
        return default
    end
end

end


text:
struct ObjectiveIndex{ F <:MOI.AbstractFunction }
	value :: Int64
end

# TODO: should `ObjectiveIndex` be a subtype of MOI.AbstractFunction`?
# this would allow for constraints involving the objectives, i.e.
# min_x f, subject to f(x) ≤ 0 ?


text:
#    Problem : GROUP A
#    *********
#    TA simple nonlinear problem given as an example in Section 1.2.4 of
#    the LANCELOT Manual.
#    The problem is non convex and has several local minima.
#
#    Original SIF Source:
#    A.R. Conn, N. Gould and Ph.L. Toint,
#    "LANCELOT, A Fortran Package for Large-Scale Nonlinear Optimization
#    (Release A)"
#    Springer Verlag, 1992.
#
#    Implementation translated from Source:
#    http://eprints.tsu.ge/234/14/Tests%20collection-K-F.pdf
#
#    EG2.SIF classification OUR2-AN-1000-0
#
#    Number of variables (at least 2)
#
# Daniel Henderson, 08/2021

f = (x) -> begin
    fx = 0.0
    α = x[1] - 1.0
    for i in firstindex(x):lastindex(x)-1
        γ = α + x[i]^2
        fx += sin(γ)
    end
    return fx
end

g! = (g, x) -> begin
    α = x[1] - 1
    xn = x[lastindex(x)]
    g[lastindex(x)] = cos(xn^2)xn
    for i in firstindex(x):lastindex(x)-1
        γ = α + x[i]^2
        g[1] += cos(γ)
        g[i] += 2.0cos(γ)x[i]
    end
    return g
end

fg! = (g, x) -> begin
    fx = 0.0
    α = x[1] - 1.0
    xn = x[lastindex(x)]
    g[lastindex(x)] = cos(xn^2)xn
    for i in firstindex(x):lastindex(x)-1
        γ = α + x[i]^2
        fx += sin(γ)
        g[1] += cos(γ)
        g[i] += 2.0cos(γ)x[i]
    end
    return fx, g
end

init = (n::Int=1000) -> begin
    n < 2 && @warn("EG2: number of variables must be ≥ 2")
    n = max(n, 2)

    return n, zeros(n)
end

#TestSet["EG2"] = UncProgram("EG2", f, g!, fg!, init)


text:
using SCEDC
using Test

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


text:
# load packages
using CSV;
using DataFrames;

# setup filepaths
path_source = string(@__DIR__,"\\..\\source");
path_dev = string(@__DIR__,"\\..\\dev");

# merge files
df_data_appended = CSV.read(joinpath(path_dev,"19-21VAERSDATA.csv"), DataFrame);
df_vax_appended = CSV.read(joinpath(path_dev,"19-21VAERSVAX.csv"), DataFrame);
df_symp_appended = CSV.read(joinpath(path_dev,"19-21VAERSSYMPTOMS_dev.csv"), DataFrame);

# rename colnames of df_symp_appended
colnames = ["VAERS_ID", "SYMPTOMS"];
rename!(df_symp_appended, colnames)

# merge (left-join) VAERS data, vaccine files
print(size(df_data_appended));
print(size(df_vax_appended));
df_merged = leftjoin(df_data_appended, df_vax_appended, on = :VAERS_ID)
print(size(df_merged));

# merge (left-join) VAERS data + vax, symptoms files
print(size(df_merged));
print(size(df_symp_appended));
df_merged = leftjoin(df_merged, df_symp_appended, on = :VAERS_ID)
print(size(df_merged));

# write to dev
CSV.write(joinpath(path_dev,"19-21VAERSCOMB.csv"), df_merged);


text:
using Metatheory
using Metatheory.Rewriters
using Test
using TermInterface
# using SymbolicUtils

array_theory = @theory x y f g M N begin
    #map(f,x)[n:m] = map(f,x[n:m]) # but does NOT commute with filter
    map(f,fill(x,N))            == fill(apply(f,x), N) # hmm
    # cumsum(fill(x,N))           == collect(x:x:(N*x))
    fill(x,N)[y]                --> x
    length(fill(x,N))           --> N
    reverse(reverse(x))         --> x
    sum(fill(x,N))              --> x * N
    map(f,reverse(x))           == reverse(map(f, x))
    filter(f,reverse(x))        == reverse(filter(f,x))
    reverse(fill(x,N))          == fill(x,N)
    filter(f, fill(x,N))        == (if apply(f, x); fill(x,N) else fill(x,0) end)
    filter(f, filter(g, x))     == filter(fand(f,g), x) # using functional &&
    cat(fill(x,N),fill(x,M))    == fill(x,N + M)
    cat(map(f,x), map(f,y))     == map(f, cat(x,y))
    map(f, cat(x,y))            == cat(map(f,x), map(f,y))
    map(f,map(g,x))             == map(f ∘ g, x)
    reverse( cat(x,y) )         == cat(reverse(y), reverse(x))
    map(f,x)[y]                 == apply(f,x[y])
    apply(f ∘ g, x)             == apply(f, apply(g, x))

    reduce(g, map(f, x))        == mapreduce(f, g, x)
    foldl(g, map(f, x))         == mapfoldl(f, g, x)
    foldr(g, map(f, x))         == mapfoldr(f, g, x)
end

asymptot_t = @theory x y z n m f g begin
    (length(filter(f, x)) <= length(x)) => true
    length(cat(x, y)) --> length(x) + length(y)
    length(map(f, x)) => length(map)
    length(x::UnitRange) => length(x)
end

fold_theory = @theory x y z begin
    x::Number * y::Number => x*y
    x::Number + y::Number => x+y
    x::Number / y::Number => x/y
    x::Number - y::Number => x/y
    # etc...
end

# Simplify expressions like :(d->3:size(A,d)-3) given an explicit value for d
import Base.Cartesian: inlineanonymous

tryinlineanonymous(x) = nothing
function tryinlineanonymous(ex::Expr)
    exprhead(ex) != :call && return nothing
    f = operation(ex)
    (!(f isa Expr) || exprhead(f) !== :->) && return nothing
    arg = arguments(ex)[1]
    println(arg)
    try
        return inlineanonymous(f, arg)
    catch e
        return nothing
    end
end

normalize_theory = @theory x y z f g begin
    fand(f, g)  => Expr(:->, :x, :(($f)(x) && ($g)(x)))
    apply(f, x) => Expr(:call, f, x)
end

params = SaturationParams()

function stream_optimize(ex)
    g = EGraph(ex)
    rep = saturate!(g, array_theory, params)
    @info rep
    ex = extract!(g, astsize) # TODO cost fun with asymptotic complexity
    ex = Fixpoint(Postwalk(Chain([tryinlineanonymous, normalize_theory..., fold_theory...])))(ex)
    return ex
end

build_fun(ex) = eval(:(()->$ex))

@testset "Stream Fusion" begin
    ex = :( map(x -> 7 * x, fill(3,4)))
    opt = stream_optimize(ex)
    @test opt == :(fill(21, 4))

    ex = :( map(x -> 7 * x, fill(3,4) )[1])
    opt = stream_optimize(ex)
    @test opt == 21
end

# ['a','1','2','3','4']
ex = :(filter(ispow2, filter(iseven, reverse(reverse(fill(4, 100))))))
opt = stream_optimize(ex)

ex = :( map(x -> 7 * x, reverse(reverse(fill(13,40))) ))
opt = stream_optimize(ex)
opt = stream_optimize(opt)

macro stream_optimize(ex)
    stream_optimize(ex)
end


text:
# Calibrate convective adjustment closure parameters to LESbrary 2-day "free_convection" simulation

using OceanBoundaryLayerParameterEstimation
using ParameterEstimocean
using ParameterEstimocean.Transformations: Transformation
using LinearAlgebra, CairoMakie, DataDeps
using Oceananigans
using Oceananigans.Units
using Plots
using PyPlot
using PlotlyJS

using ParameterEstimocean.Parameters: unconstrained_prior

architecture = GPU()
Nensemble = 25
forward_map_description = ""
Δt = 10.0
iterations = 10

###
### Build an observation from "free convection" LESbrary simulation
###

data_path = datadep"two_day_suite_4m/free_convection_instantaneous_statistics.jld2" # Nz = 64
data_path_highres = datadep"two_day_suite_2m/free_convection_instantaneous_statistics.jld2" # Nz = 128

transformation = (b = Transformation(normalization=ZScore()),)

observation, observation_highres = SyntheticObservations.([data_path, data_path_highres];
                                    field_names=(:b,),
                                    times=[3hours, 12hours, 48hours],
                                    transformation,
                                    architecture,
                                    regrid=(1, 1, 32)
                                    )

observations = [observation]
output_map = ConcatenatedOutputMap()
y = [observation_map(output_map, observation)...]

Nobs = length(y)

# noise_cov_name = "noise_covariance_128_64_resolutions_to_32"
# Γy = estimate_η_covariance(output_map, [observation, observation_highres]) .+ Matrix(1e-10 * I, Nobs, Nobs)

noise_cov_name = "noise_covariance_01"
Γy = Matrix(0.01 * I, Nobs, Nobs)

directory = "calibrate_convadj_to_lesbrary/$(iterations)_iters_$(Nensemble)_particles_$(forward_map_description)/$(noise_cov_name)/"
isdir(directory) || mkpath(directory)

begin
    f = CairoMakie.Figure(resolution = (500, 600), fontsize = 24)
    ax1 = Axis(f[1, 1], title = "y")
    ax2 = Axis(f[1, 2], title = "diag(Γy)")

    lines!(ax1, y, 1:length(y), color=:purple, linewidth=4)

    noise_var = diag(Γy)
    lines!(ax2, noise_var, 1:length(noise_var), color=:purple, linewidth=4)

    save(joinpath(directory, "obs_and_noise_var.png"), f)
end

f = CairoMakie.heatmap(Γy)
save(joinpath(directory, "noise_cov_heatmap.png"), f)

###
### Build Inverse Problem
###

priors = (convective_κz = ScaledLogitNormal(bounds=(0.0, 2.0)),
          background_κz = ScaledLogitNormal(bounds=(0.0, 2.5e-3)))
free_parameters = FreeParameters(priors)

closure = ConvectiveAdjustmentVerticalDiffusivity()

simulation = lesbrary_ensemble_simulation(observations; Nensemble, architecture, closure, Δt)
calibration = InverseProblem(observations, simulation, free_parameters; output_map = output_map)

function build_inverse_problem(Nensemble)
    simulation = lesbrary_ensemble_simulation(observations; Nensemble, architecture, closure, Δt)
    calibration = InverseProblem(observations, simulation, free_parameters; output_map)
    return calibration
end

calibration = build_inverse_problem(Nensemble)

###
### Run EKI
###

unconstrained_prior1 = unconstrained_prior(priors[:convective_κz])
unconstrained_prior2 = unconstrained_prior(priors[:background_κz])

# n = 10
# uniform_unit_samples = collect(range(0.0; stop=1.0, length=n+2)[2:n+1])
# samples1 = quantile(unconstrained_prior1, uniform_unit_samples)
# samples2 = quantile(unconstrained_prior2, uniform_unit_samples)
# initial_ensemble = hcat([[s1, s2] for s1 in samples1, s2 in samples2]...)

eki = EnsembleKalmanInversion(calibration; noise_covariance = Γy)
iterate!(eki; iterations = iterations)

###
### Summary Plots
###

plot_parameter_convergence!(eki, directory)
plot_pairwise_ensembles!(eki, directory)
plot_error_convergence!(eki, directory)

visualize!(calibration, eki.iteration_summaries[end].ensemble_mean;
    field_names = (:b,),
    directory = directory,
    filename = "realizations.png"
)

###
### Visualize loss landscape
###

name = "Loss landscape"

# pvalues = Dict(
#     :convective_κz => collect(0.1:0.05:2.0),
#     :background_κz => collect(0.5e-4:0.5e-4:1e-3),
# )

pvalues = Dict(
    :convective_κz => collect(0.1:0.05:2.1),
    :background_κz => collect(0.5e-4:1e-4:2.5e-3),
)

p1 = pvalues[:convective_κz]
p2 = pvalues[:background_κz]

ni = length(p1)
nj = length(p2)

params = hcat([[p1[i], p2[j]] for i = 1:ni, j = 1:nj]...)
xc = params[1, :]
yc = params[2, :]

# build an `InverseProblem` that can accommodate `ni*nj` ensemble members
big_calibration = build_inverse_problem(ni*nj)

y = observation_map(big_calibration)

using FileIO

file = "calibrate_convadj_to_lesbrary/loss_landscape_$(forward_map_description).jld2"

# @time G = forward_map(big_calibration, params)
G = load(file)["G"]

θglobalmin = NamedTuple((:convective_κz => 0.275, :background_κz => 0.000275))
visualize!(calibration, θglobalmin;
    field_names = (:b,),
    directory = directory,
    filename = "realizations_θglobalmin.png"
)

using ParameterEstimocean.EnsembleKalmanInversions: eki_objective
Φs = [eki_objective(eki, params[:,j], G[:,j]; constrained = true) for j in 1:size(G, 2)]

save(file, Dict("G" => G,
                noise_cov_name*"/Φ2" => getindex.(Φs, 2),
                noise_cov_name*"/Φ1" => getindex.(Φs, 1)))

function plot_contour(eki, xc, yc, zc, name, directory; zlabel="MSE loss", title="EKI Particle Traversal Over Loss Landscape", plot_minimizer=false, plot_scatters=true)

    # 2D contour plot with EKI particles superimposed
    begin
        f = CairoMakie.Figure()
        ax1 = Axis(f[1, 1],
            title = title,
            xlabel = "convective_κz",
            ylabel = "background_κz")

        co = CairoMakie.contourf!(ax1, xc, yc, zc, levels = 50, colormap = :default)

        cvt(iter) = hcat(collect.(eki.iteration_summaries[iter].parameters)...)
        diffc = cvt(2) .- cvt(1)
        diff_mag = mapslices(norm, diffc, dims = 1)
        us = diffc[1, :]
        vs = diffc[2, :]
        xs = cvt(1)[1, :]
        ys = cvt(1)[2, :]

        plot_scatters && arrows!(xs, ys, us, vs, arrowsize = 10, lengthscale = 0.3,
            arrowcolor = :yellow, linecolor = :yellow)

        legend_labels = Vector{String}([])
        scatters = []

        if plot_minimizer

            # Ignore all NaNs
            not_nan_indices = findall(.!isnan.(zc))
            xc_no_nans = xc[not_nan_indices]
            yc_no_nans = yc[not_nan_indices]
            zc_no_nans = zc[not_nan_indices]

            am = argmin(zc_no_nans)
            minimizing_params = [xc_no_nans[am] yc_no_nans[am]]
            push!(scatters, CairoMakie.scatter!(ax1, minimizing_params, marker = :x, markersize = 30))
            push!(legend_labels, "Global minimum")
        end

        if plot_scatters
            for (i, iteration) in enumerate([0, 1, iterations])
                ensemble = eki.iteration_summaries[iteration].parameters
                ensemble = [[particle[:convective_κz], particle[:background_κz]] for particle in ensemble]
                ensemble = transpose(hcat(ensemble...)) # N_ensemble x 2
                push!(scatters, CairoMakie.scatter!(ax1, ensemble))
            end
            legend_labels = vcat(legend_labels, ["Initial ensemble", "Iteration 1", "Iteration $(iterations)"])
        end

        if plot_minimizer || plot_scatters
            Legend(f[1, 2], scatters, legend_labels, position = :lb)
        end

        save(joinpath(directory, "loss_contour_$(name).png"), f)
    end

    # 3D loss landscape
    begin
        f = CairoMakie.Figure()
        ax1 = Axis3(f[1, 1],
            title = title,
            xlabel = "convective_κz",
            ylabel = "background_κz",
            zlabel = zlabel,
            shading = true,
            shininess = 32.0,
            transparency = true
        )

        CairoMakie.surface!(ax1, xc, yc, zc)
        save(joinpath(directory, "loss_landscape_$(name).png"), f)

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

        # layout = PlotlyJS.Layout(
        #     title=title,
        #     xtitle="convective_κz",
        #     ytitle="background_κz",
        #     autosize=false,
        #     scene_camera_eye=attr(x=-1.25, y=-1.25, z=1.25),
        #     width=250, height=250
        # )
        # p = PlotlyJS.plot(PlotlyJS.surface(
        #     x=collect(1:ni),
        #     y=collect(1:nj),
        #     z=reshape(zc, (ni,nj)),
        #     opacity=0.9,
        #     contours_z=attr(
        #         show=true,
        #         usecolormap=true,
        #         highlightcolor="limegreen",
        #         project_z=true
        #     )
        # ), layout)
        # PlotlyJS.savefig(p, joinpath(directory, "v3_loss_landscape_$(name).png"), width=1000, height=1000)

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

        # xgrid = reshape(xc, (ni,nj))
        # ygrid = reshape(yc, (ni,nj))
        # zgrid = reshape(zc, (ni,nj))
        # fig = PyPlot.figure("pyplot_surfaceplot",figsize=(10,10))
        # ax = fig.add_subplot(2,1,1,projection="3d")
        # ax.plot_surface(xgrid, ygrid, zgrid, rstride=2, cstride=2, edgecolors="k", cmap="autumn_r", alpha=0.8, linewidth=0.0)
        # PyPlot.xlabel("convective_κz")
        # PyPlot.ylabel("background_κz")
        # PyPlot.zlabel(zlabel)
        # PyPlot.title("Surface Plot")

        # ax.contourf(xgrid, ygrid, zgrid, 10, linewidth=3, cmap="autumn_r", linestyles="solid", offset=-1)

        # PyPlot.savefig(joinpath(directory, "v3_loss_landscape_$(name).png"))
        # PyPlot.close(fig)

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

        # ax.contourf(xgrid, ygrid, zgrid, 10, linewidth=3, cstride=1, rstride=1, cmap="autumn_r", linestyles="solid", offset=-1)
        # ax.contourf(xgrid, ygrid, zgrid, 10, linewidth=3, colors="k", linestyles="solid")

        # PyPlot.subplot(212)
        # ax = fig.add_subplot(2,1,2)
        # cp = PyPlot.contour(xgrid, ygrid, zgrid, colors="black")
        # ax.clabel(cp, inline=1, fontsize=10)
        # PyPlot.xlabel("convective_κz")
        # PyPlot.ylabel("background_κz")
        # PyPlot.title("Contour Plot")
        # PyPlot.tight_layout()

    end
end

G = load(file)["G"]
zc = [mapslices(norm, G .- y, dims = 1)...]

Φ1 = load(file)[noise_cov_name*"/Φ1"]
Φ2 = load(file)[noise_cov_name*"/Φ2"]

plot_contour(eki, xc, yc, zc, "norm", directory; zlabel = "|G(θ) - y|", title="Mean Square Error", plot_minimizer=true, plot_scatters=false)
plot_contour(eki, xc, yc, Φ1, "Φ1", directory; zlabel = "Φ1", plot_scatters=false, title="Φ1")
plot_contour(eki, xc, yc, Φ2, "Φ2", directory; zlabel = "Φ2", plot_scatters=false, title="Φ2")
plot_contour(eki, xc, yc, Φ1 .+ Φ2, "Φ", directory; zlabel = "Φ", plot_minimizer=true, plot_scatters=false, title="EKI Objective Function, Φ")

particles_overlaid_directory = joinpath(directory, "with_particles_overlaid")
isdir(particles_overlaid_directory) || mkpath(particles_overlaid_directory)
plot_contour(eki, xc, yc, zc, "norm", particles_overlaid_directory; zlabel = "|G(θ) - y|", title="Mean Square Error", plot_minimizer=true, plot_scatters=true)
plot_contour(eki, xc, yc, Φ1, "Φ1", particles_overlaid_directory; zlabel = "Φ1", plot_scatters=true, title="Φ1")
plot_contour(eki, xc, yc, Φ2, "Φ2", particles_overlaid_directory; zlabel = "Φ2", plot_scatters=true, title="Φ2")
plot_contour(eki, xc, yc, Φ1 .+ Φ2, "Φ", particles_overlaid_directory; zlabel = "Φ", plot_minimizer=true, plot_scatters=true, title="EKI Objective Function, Φ")


text:
using Test
using TinyRPC

count = 0

function foo(x)
    global count
    count += x
    x
end

function bar(; extra="")
    return "BAR$extra"
end

echo(x) = x
echo(args...) = args

function delay()
    @info("delay start")
    sleep(2)
    @info("delay end")
    7
end

@testset "TinyRPC" begin

    port = 2000 + rand(UInt8)

    server, clients = TinyRPC.listen(;port=port)

    io = TinyRPC.connect("localhost"; port=port)

    r = TinyRPC.@remote io vcat([1, 2, 3], 4, 5, 6)
    @test r == [1, 2, 3, 4, 5, 6]

    r = TinyRPC.remote_call(io, :vcat, [1, 2, 3], 4, 5, 6)
    @test r == [1, 2, 3, 4, 5, 6]

    @show clients
    @show clients[1]
    r = TinyRPC.@remote clients[1] vcat([1, 2, 3], 4, 5, 6)
    @test r == [1, 2, 3, 4, 5, 6]

    r = TinyRPC.remote_call(clients[1], :vcat, [1, 2, 3], 4, 5, 6)
    @test r == [1, 2, 3, 4, 5, 6]

    global count
    count = 0
    r = TinyRPC.@remote io foo(7)
    @test r == 7
    @test count == 7

    r = TinyRPC.remote_call(io, :foo, 7)
    @test r == 7
    @test count == 14

    r = TinyRPC.@remote clients[1] foo(2)
    @test r == 2
    @test count == 16

    @test TinyRPC.remote_call(io, :bar) == "BAR"

    @test TinyRPC.remote_call(io, :echo, 7) == 7

    @test TinyRPC.remote_call(io, :echo, 1, 2, 3) == (1, 2, 3)

    @test TinyRPC.remote_call(io, :echo, :FOO) == :FOO

    @test TinyRPC.remote_call(io, :echo, "BAR") == "BAR"

    the_err = nothing
    @async try
        TinyRPC.remote_call(io, :delay)
    catch err
        @info "remote_call error", err
        the_err = err
    end
    sleep(0.2)
    close(io.io)
    sleep(0.2)
    @test the_err isa EOFError

    @test TinyRPC.remote_call(io, :echo, "BAR") == "BAR"
    @test TinyRPC.remote_call(io, :echo, 1, 2, 3) == (1, 2, 3)

    @test TinyRPC.remote_call(io, :bar) == "BAR"
    @test TinyRPC.remote_call(io, :bar; extra="_MORE") == "BAR_MORE"
end

module AMod

    using TinyRPC
    using Test

    module TheMod
        export mod_foo, mod_bar
        mod_foo(x) = "foo $x"
        mod_bar() = "bar"
    end

    port = 2001 + rand(UInt8)
    server, clients = TinyRPC.listen(;port=port, mod=TheMod)
    io = TinyRPC.connect("localhost"; port=port)

    TinyRPC.@remote_using io :TheMod

    @testset "TinyRPC remote_using" begin

        @test mod_bar() == "bar"
        @test mod_foo(8) == "foo 8"
        @test mod_foo("bar") == "foo bar"
    end
end


text:
#FileIO2 tools to manipulate text files
#-------------------------------------------------------------------------------

#==Implements a text reader.

NOTE: Uses inefficient high-level tools (readline/split) to provide functionnality.
==#

#==Constants
===============================================================================#
const splitter_default = [' ', '\t']
const splitter_csv = [' ', '\t', ',']

#==Main data structures
===============================================================================#
mutable struct TextReader{T<:TextFormat} <: AbstractReader{T}
	s::IO
	splitter
	linebuf::Vector
end

#Define default TextReader - when text format not (fully) specified:
(::Type{TextReader})(s::IO, splitter, linebuf::Vector) =
	TextReader{TextFormat{UnknownTextEncoding}}(s, splitter, linebuf)
(::Type{TextReader{TextFormat}})(s::IO, splitter, linebuf::Vector) =
	TextReader{TextFormat{UnknownTextEncoding}}(s, splitter, linebuf)

#Default splitter value:
(RT::Type{T})(s::IO, splitter = splitter_default) where T<:TextReader =
	RT(s, splitter, String[])

#==Helper functions
===============================================================================#
#Advance TextReader linebuf, if necessary:
#TODO: Re-implement with more efficient, lower level functions
function refreshbuf(r::TextReader)
	while length(r.linebuf) < 1
		if eof(r.s); throw(EOFError()); end
		linebuf = split(chomp(readline(r.s)), r.splitter)

		for token in linebuf
			if token != ""
				push!(r.linebuf, token)
			end
		end
	end
end

#==Open/read/close functions
===============================================================================#

open(RT::Type{T}, path::String, splitter = splitter_default) where T<:TextReader =
	RT(open(path, "r"), splitter)
open(RT::Type{TextReader{T}}, path::String, splitter = splitter_csv) where T<:CSVFormat =
	RT(open(path, "r"), splitter)

#Read in entire text file as string
function read(RT::Type{T}, path::String) where T<:TextReader
	open(RT, path) do reader
		return read(reader.s, String)
	end
end

#Read in next token as String:
function read(r::TextReader, ::Type{DT}) where DT<:AbstractString
	refreshbuf(r)
	v = popfirst!(r.linebuf)
	try
		#Read ahead if necessary (Update state for eof()):
		refreshbuf(r)
	catch
	end
	return v
end

#Read in next token & interpret as of type DT:
function read(r::TextReader, ::Type{DT}) where DT
	return parse(DT, read(r, String))
end

#Read in next token & interpret as most appropriate type:
function read(r::TextReader, ::Type{Any})
	return parse(read(r, String))
end

close(r::TextReader) = close(r.s)

#==Support functions
===============================================================================#
eof(r::TextReader) = (eof(r.s) && length(r.linebuf) < 1)

function Base.readline(r::TextReader)
	linebuf = []
	return readline(r.s)
end

function Base.read(r::TextReader, String)
	linebuf = []
	return read(r.s, String)
end

#Last line


text:
# Scalar product and norm: only valid for EuclideanSpace
Base.dot{S<:EuclideanSpace}(t1::Tensor{S},t2::Tensor{S})= (space(t1)==space(t2) ? dot(vec(t1),vec(t2)) : throw(SpaceError()))
Base.vecnorm{S<:EuclideanSpace}(t::Tensor{S})=vecnorm(t.data) # frobenius norm

# # Factorizations:
# #-----------------
# for (S,TT) in ((CartesianSpace,CartesianTensor),(ComplexSpace,ComplexTensor))
#     @eval function svd!(t::$TT,n::Int,::NoTruncation)
#         # Perform singular value decomposition corresponding to bipartion of the
#         # tensor indices into the left indices 1:n and remaining right indices,
#         # thereby destroying the original tensor.
#         N=numind(t)
#         spacet=space(t)
#         leftspace=spacet[1:n]
#         rightspace=spacet[n+1:N]
#         leftdim=dim(leftspace)
#         rightdim=dim(rightspace)
#         data=reshape(t.data,(leftdim,rightdim))
#         F=svdfact!(data)
#         newdim=length(F[:S])
#         newspace=$S(newdim)
#         U=tensor(F[:U],leftspace*newspace')
#         S=tensor(diagm(F[:S]),newspace*newspace')
#         V=tensor(F[:Vt],newspace*rightspace)
#         return U,S,V,abs(zero(eltype(t)))
#     end
#
#     @eval function svd!(t::$TT,n::Int,trunc::Union(TruncationDimension,TruncationSpace))
#         # Truncate rank corresponding to bipartition into left indices 1:n
#         # and remain right indices, based on singular value decomposition,
#         # thereby destroying the original tensor.
#         # Truncation parameters are given as keyword arguments: trunctol should
#         # always be one of the possible arguments for specifying truncation, but
#         # truncdim can be replaced with different parameters for other types of tensors.
#
#         N=numind(t)
#         spacet=space(t)
#         leftspace=spacet[1:n]
#         rightspace=spacet[n+1:N]
#         leftdim=dim(leftspace)
#         rightdim=dim(rightspace)
#         dim(trunc) >= min(leftdim,rightdim) && return svd!(t,n)
#         data=reshape(t.data,(leftdim,rightdim))
#         F=svdfact!(data)
#         sing=F[:S]
#
#         # compute truncation dimension
#         if eps(trunc)!=0
#             sing=F[:S]
#             normsing=vecnorm(sing)
#             truncdim=0
#             while truncdim<length(sing) && vecnorm(sing[(truncdim+1):end])>eps(trunc)*normsing
#                 truncdim+=1
#             end
#             truncdim=min(truncdim,dim(trunc))
#         else
#             truncdim=dim(trunc)
#         end
#         newspace=$S(truncdim)
#
#         # truncate
#         truncerr=vecnorm(sing[(truncdim+1):end])/vecnorm(sing)
#         U=tensor(F[:U][:,1:truncdim],leftspace*newspace')
#         S=tensor(diagm(sing[1:truncdim]),newspace*newspace')
#         V=tensor(F[:Vt][1:truncdim,:],newspace*rightspace)
#         return U,S,V,truncerr
#     end
#
#     @eval function svd!(t::$TT,n::Int,trunc::TruncationError)
#         # Truncate rank corresponding to bipartition into left indices 1:n
#         # and remain right indices, based on singular value decomposition,
#         # thereby destroying the original tensor.
#         # Truncation parameters are given as keyword arguments: trunctol should
#         # always be one of the possible arguments for specifying truncation, but
#         # truncdim can be replaced with different parameters for other types of tensors.
#
#         N=numind(t)
#         spacet=space(t)
#         leftspace=spacet[1:n]
#         rightspace=spacet[n+1:N]
#         leftdim=dim(leftspace)
#         rightdim=dim(rightspace)
#         data=reshape(t.data,(leftdim,rightdim))
#         F=svdfact!(data)
#
#         # find truncdim based on trunctolinfo
#         sing=F[:S]
#         normsing=vecnorm(sing)
#         truncdim=0
#         while truncdim<length(sing) && vecnorm(sing[(truncdim+1):end])>eps(trunc)*normsing
#             truncdim+=1
#         end
#         newspace=$S(truncdim)
#
#         # truncate
#         truncerr=vecnorm(sing[(truncdim+1):end])/normsing
#         U=tensor(F[:U][:,1:truncdim],leftspace*newspace')
#         S=tensor(diagm(sing[1:truncdim]),newspace*newspace')
#         V=tensor(F[:Vt][1:truncdim,:],newspace*rightspace)
#         return U,S,V,truncerr
#     end
#
#     @eval function leftorth!(t::$TT,n::Int,::NoTruncation)
#         # Create orthogonal basis U for indices 1:n, and remainder C for right
#         # indices, thereby destroying the original tensor.
#         # Decomposition should be unique, such that it always returns the same
#         # result for the same input tensor t. UC = QR is fastest but only unique
#         # after correcting for phases.
#         N=numind(t)
#         spacet=space(t)
#         leftspace=spacet[1:n]
#         rightspace=spacet[n+1:N]
#         leftdim=dim(leftspace)
#         rightdim=dim(rightspace)
#         data=reshape(t.data,(leftdim,rightdim))
#         if leftdim>rightdim
#             newdim=rightdim
#             tau=zeros(eltype(data),(newdim,))
#             Base.LinAlg.LAPACK.geqrf!(data,tau)
#
#             C=zeros(eltype(t),(newdim,newdim))
#             for j in 1:newdim
#                 for i in 1:j
#                     @inbounds C[i,j]=data[i,j]
#                 end
#             end
#             Base.LinAlg.LAPACK.orgqr!(data,tau)
#             U=data
#             for i=1:newdim
#                 tau[i]=sign(C[i,i])
#             end
#             scale!(U,tau)
#             scale!(conj!(tau),C)
#         else
#             newdim=leftdim
#             C=data
#             U=eye(eltype(data),newdim)
#         end
#
#         newspace=$S(newdim)
#         return tensor(U,leftspace*newspace'), tensor(C,newspace*rightspace), abs(zero(eltype(t)))
#     end
#
#     @eval function leftorth!(t::$TT,n::Int,trunc::Union(TruncationDimension,TruncationSpace))
#         # Truncate rank corresponding to bipartition into left indices 1:n
#         # and remain right indices, based on singular value decomposition,
#         # thereby destroying the original tensor.
#         # Truncation parameters are given as keyword arguments: trunctol should
#         # always be one of the possible arguments for specifying truncation, but
#         # truncdim can be replaced with different parameters for other types of tensors.
#
#         N=numind(t)
#         spacet=space(t)
#         leftspace=spacet[1:n]
#         rightspace=spacet[n+1:N]
#         leftdim=dim(leftspace)
#         rightdim=dim(rightspace)
#         dim(trunc) >= min(leftdim,rightdim) && return leftorth!(t,n)
#         data=reshape(t.data,(leftdim,rightdim))
#         F=svdfact!(data)
#         sing=F[:S]
#
#         # compute truncation dimension
#         if eps(trunc)!=0
#             sing=F[:S]
#             normsing=vecnorm(sing)
#             truncdim=0
#             while truncdim<length(sing) && vecnorm(sing[(truncdim+1):end])>eps(trunc)*normsing
#                 truncdim+=1
#             end
#             truncdim=min(truncdim,dim(trunc))
#         else
#             truncdim=dim(trunc)
#         end
#         newspace=$S(truncdim)
#
#         # truncate
#         truncerr=vecnorm(sing[(truncdim+1):end])/vecnorm(sing)
#         U=tensor(F[:U][:,1:truncdim],leftspace*newspace')
#         C=tensor(scale!(sing[1:truncdim],F[:Vt][1:truncdim,:]),newspace*rightspace)
#         return U,C,truncerr
#     end
#
#     @eval function leftorth!(t::$TT,n::Int,trunc::TruncationError)
#         # Truncate rank corresponding to bipartition into left indices 1:n
#         # and remain right indices, based on singular value decomposition,
#         # thereby destroying the original tensor.
#         # Truncation parameters are given as keyword arguments: trunctol should
#         # always be one of the possible arguments for specifying truncation, but
#         # truncdim can be replaced with different parameters for other types of tensors.
#
#         N=numind(t)
#         spacet=space(t)
#         leftspace=spacet[1:n]
#         rightspace=spacet[n+1:N]
#         leftdim=dim(leftspace)
#         rightdim=dim(rightspace)
#         data=reshape(t.data,(leftdim,rightdim))
#         F=svdfact!(data)
#
#         # compute truncation dimension
#         sing=F[:S]
#         normsing=vecnorm(sing)
#         truncdim=0
#         while truncdim<length(sing) && vecnorm(sing[(truncdim+1):end])>eps(trunc)*normsing
#             truncdim+=1
#         end
#         newspace=$S(truncdim)
#
#         # truncate
#         truncerr=vecnorm(sing[(truncdim+1):end])/vecnorm(sing)
#         U=tensor(F[:U][:,1:truncdim],leftspace*newspace')
#         C=tensor(scale!(sing[1:truncdim],F[:Vt][1:truncdim,:]),newspace*rightspace)
#
#         return U,C,truncerr
#     end
#
#     @eval function rightorth!(t::$TT,n::Int,::NoTruncation)
#         # Create orthogonal basis U for right indices, and remainder C for left
#         # indices. Decomposition should be unique, such that it always returns the
#         # same result for the same input tensor t. CU = LQ is typically fastest but only
#         # unique after correcting for phases.
#         N=numind(t)
#         spacet=space(t)
#         leftspace=spacet[1:n]
#         rightspace=spacet[n+1:N]
#         leftdim=dim(leftspace)
#         rightdim=dim(rightspace)
#         data=reshape(t.data,(leftdim,rightdim))
#         if leftdim<rightdim
#             newdim=leftdim
#             tau=zeros(eltype(data),(newdim,))
#             datat=transpose(data)
#             Base.LinAlg.LAPACK.geqrf!(datat,tau)
#
#             C=zeros(eltype(t),(newdim,newdim))
#             for j in 1:newdim
#                 for i in 1:j
#                     @inbounds C[j,i]=datat[i,j]
#                 end
#             end
#             Base.LinAlg.LAPACK.orgqr!(datat,tau)
#             Base.transpose!(data,datat)
#             U=data
#
#             for i=1:newdim
#                 tau[i]=sign(C[i,i])
#             end
#             scale!(tau,U)
#             scale!(C,conj!(tau))
#         else
#             newdim=rightdim
#             C=data
#             U=eye(eltype(data),newdim)
#         end
#
#         newspace=$S(newdim)
#         return tensor(C,leftspace ⊗ dual(newspace)), tensor(U,newspace ⊗ rightspace), abs(zero(eltype(t)))
#     end
#
#     @eval function rightorth!(t::$TT,n::Int,trunc::Union(TruncationDimension,TruncationSpace))
#         # Truncate rank corresponding to bipartition into left indices 1:n
#         # and remain right indices, based on singular value decomposition,
#         # thereby destroying the original tensor.
#         # Truncation parameters are given as keyword arguments: trunctol should
#         # always be one of the possible arguments for specifying truncation, but
#         # truncdim can be replaced with different parameters for other types of tensors.
#
#         N=numind(t)
#         spacet=space(t)
#         leftspace=spacet[1:n]
#         rightspace=spacet[n+1:N]
#         leftdim=dim(leftspace)
#         rightdim=dim(rightspace)
#         dim(trunc) >= min(leftdim,rightdim) && return rightorth!(t,n)
#         data=reshape(t.data,(leftdim,rightdim))
#         F=svdfact!(data)
#         sing=F[:S]
#
#         # compute truncation dimension
#         if eps(trunc)!=0
#             sing=F[:S]
#             normsing=vecnorm(sing)
#             truncdim=0
#             while truncdim<length(sing) && vecnorm(sing[(truncdim+1):end])>eps(trunc)*normsing
#                 truncdim+=1
#             end
#             truncdim=min(truncdim,dim(trunc))
#         else
#             truncdim=dim(trunc)
#         end
#         newspace=$S(truncdim)
#
#         # truncate
#         truncerr=vecnorm(sing[(truncdim+1):end])/vecnorm(sing)
#         C=tensor(scale!(F[:U][:,1:truncdim],sing[1:truncdim]),leftspace*newspace')
#         U=tensor(F[:Vt][1:truncdim,:],newspace*rightspace)
#         return C,U,truncerr
#     end
#
#     @eval function rightorth!(t::$TT,n::Int,trunc::TruncationError)
#         # Truncate rank corresponding to bipartition into left indices 1:n
#         # and remain right indices, based on singular value decomposition,
#         # thereby destroying the original tensor.
#         # Truncation parameters are given as keyword arguments: trunctol should
#         # always be one of the possible arguments for specifying truncation, but
#         # truncdim can be replaced with different parameters for other types of tensors.
#
#         N=numind(t)
#         spacet=space(t)
#         leftspace=spacet[1:n]
#         rightspace=spacet[n+1:N]
#         leftdim=dim(leftspace)
#         rightdim=dim(rightspace)
#         data=reshape(t.data,(leftdim,rightdim))
#         F=svdfact!(data)
#
#         # compute truncation dimension
#         sing=F[:S]
#         normsing=vecnorm(sing)
#         truncdim=0
#         while truncdim<length(sing) && vecnorm(sing[(truncdim+1):end])>eps(trunc)*normsing
#             truncdim+=1
#         end
#         newspace=$S(truncdim)
#
#         # truncate
#         truncerr=vecnorm(sing[(truncdim+1):end])/vecnorm(sing)
#         C=tensor(scale!(F[:U][:,1:truncdim],sing[1:truncdim]),leftspace*newspace')
#         U=tensor(F[:Vt][1:truncdim,:],newspace*rightspace)
#         return C,U,truncerr
#     end
# end
#
# # Methods below are only implemented for CartesianMatrix or ComplexMatrix:
# #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# typealias ComplexMatrix{T} ComplexTensor{T,2}
# typealias CartesianMatrix{T} CartesianTensor{T,2}
#
# function Base.pinv(t::Union(ComplexMatrix,CartesianMatrix))
#     # Compute pseudo-inverse
#     spacet=space(t)
#     data=copy(t.data)
#     leftdim=dim(spacet[1])
#     rightdim=dim(spacet[2])
#
#     F=svdfact!(data)
#     Sinv=F[:S]
#     for k=1:length(Sinv)
#         if Sinv[k]>eps(Sinv[1])*max(leftdim,rightdim)
#             Sinv[k]=one(Sinv[k])/Sinv[k]
#         end
#     end
#     data=F[:V]*scale(F[:S],F[:U]')
#     return tensor(data,spacet')
# end
#
# function Base.eig(t::Union(ComplexMatrix,CartesianMatrix))
#     # Compute eigenvalue decomposition.
#     spacet=space(t)
#     spacet[1] == spacet[2]' || throw(SpaceError("eigenvalue factorization only exists if left and right index space are dual"))
#     data=copy(t.data)
#
#     F=eigfact!(data)
#
#     Lambda=tensor(diagm(F[:values]),spacet)
#     V=tensor(F[:vectors],spacet)
#     return Lambda, V
# end
#
# function Base.inv(t::Union(ComplexMatrix,CartesianMatrix))
#     # Compute inverse.
#     spacet=space(t)
#     spacet[1] == spacet[2]' || throw(SpaceError("inverse only exists if left and right index space are dual"))
#
#     return tensor(inv(t.data),spacet)
# end


text:
function plot(atoms::Vector{T}) where T <: Molecules.Atom

    out = 0
    for i = eachindex(atoms)
        center = Point3f0(atoms[i].xyz)
        symbol = Molecules.symbol(atoms[i])
        r = symbol == "H" ? 0.3 : 0.5
        c = atom_color[symbol]

        sphere = Sphere(center, r)

        if i == 1
            out = mesh(sphere, show_axis=false, color=c)
        else
            mesh!(sphere, color=c)
        end
    end

    # Add chemical bonds for close atoms
    for i in eachindex(atoms)
        A1 = atoms[i]
        r1 = A1.xyz
        for j in (i+1):length(atoms)
            A2 = atoms[j]
            r2 = A2.xyz
            d = √sum((r1 .- r2).^2)
            if d < 1.55
                xyz = [Point3f0(A1.xyz)]
                uvw = [Point3f0(A2.xyz - A1.xyz)]
                arrows!(xyz, uvw, linewidth=0.1, arrowsize=0.0, linecolor=stegeman)
            end
        end
    end

    return out
end


text:
mutable struct Item
    index::Int64
    weight::Int64
    profit::Int64
    solution::Bool
end

struct Interval
    first::Int64
    last::Int64
    weight::Int64
end

mutable struct GlobalData
    items::Vector{Item}
    capacity::Int64

    br::Int64
    wsb::Int64
    psb::Int64
    z::Int64
    fsort::Int64
    lsort::Int64

    stack1::Vector{Interval}
    stack2::Vector{Interval}

    estack::Vector{Int64}

    GlobalData() = new([], 0, 0, 0, 0, 0, 0, 0, [], [], [])
end

function solveKnapExpCore(data::KnapData)
    n = length(data.items)
    if n == 0 return 0, Int64[] end

    gd = GlobalData()
    gd.items = [ Item(i, data.items[i].weight, data.items[i].profit, false) for i in 1:n ]
    gd.capacity = data.capacity

    interval = partsort(gd, 1, n, 0)
    gd.fsort = interval.first
    gd.lsort = interval.last

    gd.z = heuristic(gd, 1, n)

    elebranch(gd, 0, gd.wsb - gd.capacity, gd.br - 1, gd.br)

    return gd.z + gd.psb, definesolution(gd)
end

function det(a::Matrix{Int64})
    return a[1, 1] * a[2, 2] - a[1, 2] * a[2, 1]
end

function partsort(gd::GlobalData, f::Int64, l::Int64, ws::Int64)
    d = l - f + 1
    if d > 1
        m = f + d ÷ 2

        if det([ gd.items[f].profit gd.items[f].weight; gd.items[m].profit gd.items[m].weight ]) < 0
            gd.items[f], gd.items[m] = gd.items[m], gd.items[f]
        end

        if d > 2
            if det([ gd.items[m].profit gd.items[m].weight; gd.items[l].profit gd.items[l].weight ]) < 0
                gd.items[m], gd.items[l] = gd.items[l], gd.items[m]

                if det([ gd.items[f].profit gd.items[f].weight; gd.items[m].profit gd.items[m].weight ]) < 0
                    gd.items[f], gd.items[m] = gd.items[m], gd.items[f]
                end
            end
        end
    end

    if d <= 3
        return Interval(f, l, ws)
    else
        mp = gd.items[m].profit
        mw = gd.items[m].weight
        i = f
        j = l
        wi = ws
        while true
            while true
                wi += gd.items[i].weight
                i += 1
                if det([ gd.items[i].profit gd.items[i].weight; mp mw ]) <= 0 break end
            end
            while true
                j -= 1
                if det([ gd.items[j].profit gd.items[j].weight; mp mw ]) >= 0 break end
            end
            if i > j break end
            gd.items[i], gd.items[j] = gd.items[j], gd.items[i]
        end

        if wi > gd.capacity
            push!(gd.stack2, Interval(i, l, wi))
            return partsort(gd, f, i - 1, ws)
        else
            push!(gd.stack1, Interval(f, i - 1, ws))
            return partsort(gd, i, l, wi)
        end
    end
end

function heuristic(gd::GlobalData, f::Int64, l::Int64)
    ps = 0
    ws = gd.capacity

    i = f
    while gd.items[i].weight <= ws
        ws -= gd.items[i].weight
        ps += gd.items[i].profit
        gd.items[i].solution = true
        i += 1
    end

    gd.br = i
    gd.wsb = gd.capacity - ws
    gd.psb = ps

    dz = (gd.capacity - gd.wsb) * gd.items[gd.br].profit ÷ gd.items[gd.br].weight

    empty!(gd.estack)
    gd.z = 0
    if gd.z == dz return gd.z end

    r = gd.capacity - gd.wsb
    for i in gd.br:l
        if (gd.items[i].weight <= r) && (gd.items[i].profit > gd.z)
            empty!(gd.estack)
            push!(gd.estack, gd.items[i].index)
            gd.z = gd.items[i].profit
            if gd.z == dz return gd.z end
        end
    end

    r = gd.wsb + gd.items[gd.br].weight - gd.capacity
    pb = gd.items[gd.br].profit
    for i in gd.br - 1:-1:f
        if (gd.items[i].weight >= r) && (pb - gd.items[i].profit > gd.z)
            empty!(gd.estack)
            push!(gd.estack, gd.items[gd.br].index)
            push!(gd.estack, gd.items[i].index)
            gd.z = pb - gd.items[i].profit
            if gd.z == dz return gd.z end
        end
    end

    return gd.z
end

function reduce(gd::GlobalData, first::Int64, last::Int64)
    pb = gd.items[gd.br].profit
    wb = gd.items[gd.br].weight

    q = det([ gd.z + 1 gd.capacity - gd.wsb; pb wb ])
    i = first
    j = last

    if i <= gd.br
        k = gd.fsort - 1
        while i <= j
            if det([ -gd.items[j].profit -gd.items[j].weight; pb wb ]) < q
                gd.items[i], gd.items[j] = gd.items[j], gd.items[i]
                i += 1
            else
                gd.items[j], gd.items[k] = gd.items[k], gd.items[j]
                j -= 1
                k -= 1
            end
        end
        if k == gd.fsort - 1
            gd.items[first], gd.items[k] = gd.items[k], gd.items[first]
            k -= 1
        end
        last = gd.fsort - 1
        first = k + 1
    else
        k = gd.lsort + 1
        while i <= j
            if det([ gd.items[i].profit gd.items[i].weight; pb wb ]) < q
                gd.items[i], gd.items[j] = gd.items[j], gd.items[i]
                j -= 1
            else
                gd.items[i], gd.items[k] = gd.items[k], gd.items[i]
                i += 1
                k += 1
            end
        end
        if k == gd.lsort + 1
            gd.items[last], gd.items[k] = gd.items[k], gd.items[last]
            k += 1
        end
        first = gd.lsort + 1
        last = k - 1
    end
    return first, last
end

function sorti(gd::GlobalData, stack::Vector{Interval})
    if isempty(stack) return false end

    interval = pop!(stack)
    first, last = reduce(gd, interval.first, interval.last)
    interval = partsort(gd, first, last, interval.weight)
    if interval.first < gd.fsort gd.fsort = interval.first end
    if interval.last > gd.lsort gd.lsort = interval.last end

    return true
end

function elebranch(gd::GlobalData, ps::Int64, ws::Int64, s::Int64, t::Int64)
    improved = false

    if ws <= 0
        if ps > gd.z
            improved = true
            gd.z = ps
            empty!(gd.estack)
        end

        while true
            if t > gd.lsort && !sorti(gd, gd.stack2) break end
            if det([ ps - (gd.z + 1) ws; gd.items[t].profit gd.items[t].weight ]) < 0 break end
            if elebranch(gd, ps + gd.items[t].profit, ws + gd.items[t].weight, s, t + 1)
                improved = true
                push!(gd.estack, gd.items[t].index)
            end
            t += 1
        end
    else
        while true
            if s < gd.fsort && !sorti(gd, gd.stack1) break end
            if det([ ps - (gd.z + 1) ws; gd.items[s].profit gd.items[s].weight ]) < 0 break end
            if elebranch(gd, ps - gd.items[s].profit, ws - gd.items[s].weight, s - 1, t)
                improved = true
                push!(gd.estack, gd.items[s].index)
            end
            s -= 1
        end
    end

    return improved
end

function definesolution(gd::GlobalData)
    sort!(gd.items, by = x -> x.index)
    for index in gd.estack
        gd.items[index].solution = !gd.items[index].solution
    end

    solution = Int64[]
    for item in gd.items
        if item.solution
            push!(solution, item.index)
        end
    end
    return solution
end


text:
using DrWatson
@quickactivate "StatReth"

# %%
using CSV
# using DataFrames
using StatsBase
using Distributions
using Plots
using StatsPlots
using Optim
using Turing

# include(srcdir("quap.jl"))

# %%
d = CSV.read(datadir("exp_raw/Howell_1.csv"), copycols = true)
d2 = filter(row -> row.age >= 18, d)
d2 = d[d.age .>= 18, :]

# %%
@model height(heights) = begin
    μ ~ Normal(178, 20)
    # σ ~ Uniform(0, 50)
    # heights .~ Normal(μ, σ)
    heights .~ Normal(μ, 5.0)
end

@model line(x, y) = begin
    #priors
    alpha ~ Normal(178.0, 100.0)
    beta ~ Normal(0.0, 10.0)
    s ~ Uniform(0, 50)

    #model
    mu = alpha .+ beta*x
    y .~ Normal.(mu, s)
end

m = height(d2.height)
m2 = line(d2.weight, d2.height)

# %%
sample(m, NUTS(0.65), 1000)

opt = optimize(m, MAP(), [178.0, 25.0])
coeftable(opt)
informationmatrix(opt)

using LinearAlgebra
function quap(model::Turing.Model, args...; kwargs...)
    opt = optimize(model, MAP(), args...; kwargs...)

    map = opt.values.array
    var_cov_matrix = informationmatrix(opt)
    sym_var_cov_matrix = Symmetric(var_cov_matrix)  # lest MvNormal complains, loudly
    converged = Optim.converged(opt.optim_result)

    distr = if length(map) == 1
        Normal(map[1], √sym_var_cov_matrix[1])  # Normal expects stddev
    else
        @show sym_var_cov_matrix
        MvNormal(map, sym_var_cov_matrix)       # MvNormal expects variance matrix
    end

    (coef = map, vcov = sym_var_cov_matrix, converged = converged, distr = distr)
end

r = quap(m)
r2 = quap(m2, Optim.NelderMead())

rand(MvNormal([0.0], [0.5]), 10000) |> histogram

rand(MvNormal([0.0], [0.5]), 10000)' |> histogram
rand(Normal(0.0, 0.5), 10000) |> histogram

MvNormal([0.0], [0.5])

I now also return a normal distribution to sample from since the 1D case had the `Normal(MAP[1], vcov[1])` pitfall which would incorrectly use the variance instead of the std dev.


text:
for n in [2, 4, 16]
    H = matrixdepot("hadamard", Int, n)
    @test H*H' == n*eye(n)
end

try
    matrixdepot("hadamard", 0)
catch ArgumentError
    println("hadamard: arguementerror")
end

println("'hadamard' passed test...")


text:
#=
    `MVNSampler` is used to draw from multivariate normal distribution
=#

import Base: ==
import LinearAlgebra: BlasReal
import Random

struct MVNSampler{TM<:Real,TS<:Real,TQ<:BlasReal}
    mu::Vector{TM}
    Sigma::Matrix{TS}
    Q::Matrix{TQ}
end

function MVNSampler(mu::Vector{TM}, Sigma::Matrix{TS}) where {TM<:Real,TS<:Real}
    ATOL1, RTOL1 = 1e-8, 1e-8
    ATOL2, RTOL2 = 1e-8, 1e-14

    n = length(mu)

    if size(Sigma) != (n, n) # Check Sigma is n x n
        throw(ArgumentError(
            "Sigma must be 2 dimensional and square matrix of same length to mu"
        ))
    end

    issymmetric(Sigma) || throw(ArgumentError("Sigma must be symmetric"))

    C = cholesky(Symmetric(Sigma, :L), Val(true), check=false)
    A = C.factors
    r = C.rank
    p = invperm(C.piv)

    if r == n  # Positive definite
        Q = tril!(A)[p, p]
        return MVNSampler(mu, Sigma, Q)
    end

    non_PSD_msg = "Sigma must be positive semidefinite"

    for i in r+1:n
        A[i, i] >= -ATOL1 - RTOL1 * A[1, 1] ||
            throw(ArgumentError(non_PSD_msg))
    end

    tril!(view(A, :, 1:r))
    A[:, r+1:end] .= 0
    Q = A[p, p]
    isapprox(Q*Q', Sigma; rtol=RTOL2, atol=ATOL2) ||
        throw(ArgumentError(non_PSD_msg))

    return MVNSampler(mu, Sigma, Q)
end

# methods with the optional rng argument first
Random.rand(rng::AbstractRNG, d::MVNSampler) =
    d.mu + d.Q * randn(rng, length(d.mu))
Random.rand(rng::AbstractRNG, d::MVNSampler, n::Integer) =
    d.mu .+ d.Q * randn(rng, (length(d.mu), n))

# methods to draw from `MVNSampler`
Random.rand(d::MVNSampler) = rand(Random.GLOBAL_RNG, d)
Random.rand(d::MVNSampler, n::Integer) = rand(Random.GLOBAL_RNG, d, n)

==(f1::MVNSampler, f2::MVNSampler) =
    (f1.mu == f2.mu) && (f1.Sigma == f2.Sigma) && (f1.Q == f2.Q)


text:
@testset "Degree Independent Set" begin

    g3 = StarGraph(5)
    for g in testgraphs(g3)
        c = @inferred(independent_set(g, DegreeIndependentSet()))
        @test sort(c) == [2, 3, 4, 5]
    end

    g4 = CompleteGraph(5)
    for g in testgraphs(g4)
        c = @inferred(independent_set(g, DegreeIndependentSet()))
        @test length(c)== 1 #Exactly one vertex
    end

    #PathGraph(5) with additional edge 2-5
    g5 = Graph([0 1 0 0 0; 1 0 1 0 1; 0 1 0 1 0; 0 0 1 0 1; 0 1 0 1 0])
    for g in testgraphs(g5)
        c = @inferred(independent_set(g, DegreeIndependentSet()))
        @test sort(c) == [1, 3, 5]
    end
end


text:
@testset "Launch Configuration" begin
    kernel(x) = nothing

    device = AMDGPUnative.default_device()
    queue = AMDGPUnative.default_queue(device)

    # Group/grid size selection and aliases
    eval(:(@roc groupsize=2 $kernel(1)))
    eval(:(@roc groupsize=2 gridsize=4 $kernel(1)))
    eval(:(@roc gridsize=2 $kernel(1)))
    eval(:(@roc threads=2 $kernel(1)))
    eval(:(@roc blocks=4 $kernel(1)))
    eval(:(@roc threads=2 gridsize=4 $kernel(1)))

    # Device/queue selection and aliases
    eval(:(@roc device=$device $kernel(1)))
    eval(:(@roc device=$device queue=$queue $kernel(1)))
    eval(:(@roc queue=$queue $kernel(1)))
    eval(:(@roc agent=$device $kernel(1)))
    eval(:(@roc stream=$queue $kernel(1)))
    eval(:(@roc agent=$device queue=$queue $kernel(1)))
end


text:
using EconModels

# Test promote rules and converts

test1 = Setting(:test1, 22.0) # short constructor
vint = Setting(:data_vintage, "REF", true, "vint", "Data vintage") # full constructor

@test promote_rule(Setting{Float64}, Float16) == Float64
@test promote_rule(Setting{Bool}, Bool) == Bool
@test promote_rule(Setting{ASCIIString}, AbstractString) == UTF8String
@test convert(Int64, test1) == 22
@test convert(Float32, test1) == Float32(22.)
@test convert(ASCIIString, vint) == "REF"


text:
###############################################################################
#
# An example of writing cell comments to a worksheet using XlsxWriter.
#
# Each of the worksheets demonstrates different features of cell comments.
# Original Python Copyright 2013-2016, John McNamara, jmcnamara@cpan.org
# https://github.com/jmcnamara/XlsxWriter

using XlsxWriter

function test()

	wb = Workbook("comments2.xlsx")

	ws1 = add_worksheet!(wb)
	ws2 = add_worksheet!(wb)
	ws3 = add_worksheet!(wb)
	ws4 = add_worksheet!(wb)
	ws5 = add_worksheet!(wb)
	ws6 = add_worksheet!(wb)
	ws7 = add_worksheet!(wb)
	ws8 = add_worksheet!(wb)

	text_wrap = add_format!(wb, Dict("text_wrap"=> 1, "valign"=> "top"))

	###############################################################################
	#
	# Example 1. Demonstrates a simple cell comments without formatting.
	#            comments.
	#

	# Set up some formatting.
	set_column!(ws1, "C:C", 25)
	set_row!(ws1, 2, height=50)
	set_row!(ws1, 5, height=50)

	# Simple ASCII string.
	cell_text = "Hold the mouse over this cell to see the comment."

	comment = "This is a comment."

	write!(ws1, "C3", cell_text, fmt=text_wrap)
	write_comment!(ws1, "C3", comment)

	###############################################################################
	#
	# Example 2. Demonstrates visible and hidden comments.
	#

	# Set up some formatting.
	set_column!(ws2, "C:C", 25)
	set_row!(ws2, 2, height=50)
	set_row!(ws2, 5, height=50)

	cell_text = "This cell comment is visible."
	comment = "Hello."

	write!(ws2, "C3", cell_text, fmt=text_wrap)
	write_comment!(ws2, "C3", comment, Dict("visible"=> true))

	cell_text = "This cell comment isn\"t visible (the default)."

	write!(ws2, "C6", cell_text, fmt=text_wrap)
	write_comment!(ws2, "C6", comment)

	###############################################################################
	#
	# Example 3. Demonstrates visible and hidden comments set at the worksheet
	#            level.
	#

	# Set up some formatting.
	set_column!(ws3, "C:C", 25)
	set_row!(ws3, 2, height=50)
	set_row!(ws3, 5, height=50)
	set_row!(ws3, 8, height=50)

	# Make all comments on the worksheet visible.
	show_comments!(ws3)

	cell_text = "This cell comment is visible, explicitly."
	comment = "Hello."

	write!(ws3, "C3", cell_text, fmt=text_wrap)
	write_comment!(ws3, "C3", comment, Dict("visible"=> 1))

	cell_text = "This cell comment is also visible because of show_comments()."

	write!(ws3, "C6", cell_text, fmt=text_wrap)
	write_comment!(ws3, "C6", comment)

	cell_text = "However, we can still override it locally."

	write!(ws3, "C9", cell_text, fmt=text_wrap)
	write_comment!(ws3, "C9", comment, Dict("visible"=> false))

	###############################################################################
	#
	# Example 4. Demonstrates changes to the comment box dimensions.
	#

	# Set up some formatting.
	set_column!(ws4, "C:C", 25)
	set_row!(ws4, 2, height=50)
	set_row!(ws4, 5, height=50)
	set_row!(ws4, 8, height=50)
	set_row!(ws4, 15, height=50)

	show_comments!(ws4)

	cell_text = "This cell comment is default size."
	comment = "Hello."

	write!(ws4, "C3", cell_text, fmt=text_wrap)
	write_comment!(ws4, "C3", comment)

	cell_text = "This cell comment is twice as wide."

	write!(ws4, "C6", cell_text, fmt=text_wrap)
	write_comment!(ws4, "C6", comment, Dict("x_scale"=> 2))

	cell_text = "This cell comment is twice as high."

	write!(ws4, "C9", cell_text, fmt=text_wrap)
	write_comment!(ws4, "C9", comment, Dict("y_scale"=> 2))

	cell_text = "This cell comment is scaled in both directions."

	write!(ws4, "C16", cell_text, fmt=text_wrap)
	write_comment!(ws4, "C16", comment, Dict("x_scale"=> 1.2, "y_scale"=> 0.8))

	cell_text = "This cell comment has width and height specified in pixels."

	write!(ws4, "C19", cell_text, fmt=text_wrap)
	write_comment!(ws4, "C19", comment, Dict("width"=> 200, "height"=> 20))

	###############################################################################
	#
	# Example 5. Demonstrates changes to the cell comment position.
	#
	set_column!(ws5, "C:C", 25)
	set_row!(ws5, 2, height=50)
	set_row!(ws5, 5, height=50)
	set_row!(ws5, 8, height=50)
	set_row!(ws5, 11, height=50)

	show_comments!(ws5, )

	cell_text = "This cell comment is in the default position."
	comment = "Hello."

	write!(ws5, "C3", cell_text, fmt=text_wrap)
	write_comment!(ws5, "C3", comment)

	cell_text = "This cell comment has been moved to another cell."

	write!(ws5, "C6", cell_text, fmt=text_wrap)
	write_comment!(ws5, "C6", comment, Dict("start_cell"=> "E4"))

	cell_text = "This cell comment has been moved to another cell."

	write!(ws5, "C9", cell_text, fmt=text_wrap)
	write_comment!(ws5, "C9", comment, Dict("start_row"=> 8, "start_col"=> 4))

	cell_text = "This cell comment has been shifted within its default cell."

	write!(ws5, "C12", cell_text, fmt=text_wrap)
	write_comment!(ws5, "C12", comment, Dict("x_offset"=> 30, "y_offset"=> 12))

	###############################################################################
	#
	# Example 6. Demonstrates changes to the comment background color.
	#
	set_column!(ws6, "C:C", 25)
	set_row!(ws6, 2, height=50)
	set_row!(ws6, 5, height=50)
	set_row!(ws6, 8, height=50)

	show_comments!(ws6, )

	cell_text = "This cell comment has a different color."
	comment = "Hello."

	write!(ws6 ,"C3", cell_text, fmt=text_wrap)
	write_comment!(ws6, "C3", comment, Dict("color"=> "green"))

	cell_text = "This cell comment has the default color."

	write!(ws6, "C6", cell_text, fmt=text_wrap)
	write_comment!(ws6, "C6", comment)

	cell_text = "This cell comment has a different color."

	write!(ws6, "C9", cell_text, fmt=text_wrap)
	write_comment!(ws6, "C9", comment, Dict("color"=> "#CCFFCC"))

	###############################################################################
	#
	# Example 7. Demonstrates how to set the cell comment author.
	#
	set_column!(ws7, "C:C", 30)
	set_row!(ws7, 2, height=50)
	set_row!(ws7, 5, height=50)
	set_row!(ws7, 8, height=50)

	author = ""
	cell = "C3"

	cell_text = "Move the mouse over this cell and you will see 'Cell commented by (blank)' in the status bar at the bottom"

	comment = "Hello."

	write!(ws7, cell, cell_text, fmt=text_wrap)
	write_comment!(ws7, cell, comment)

	author = "Julia"
	cell = "C6"
	cell_text = "Move the mouse over this cell and you will see 'Cell commented by Julia' in the status bar at the bottom"

	write!(ws7, cell, cell_text, fmt=text_wrap)
	write_comment!(ws7, cell, comment, Dict("author"=> author))

	###############################################################################
	#
	# Example 8. Demonstrates the need to explicitly set the row height.
	#
	# Set up some formatting.
	set_column!(ws8, "C:C", 25)
	set_row!(ws8, 2, height=80)

	show_comments!(ws8)

	cell_text = "The height of this row has been adjusted explicitly using set_row(). The size of the comment box is adjusted accordingly by XlsxWriter."

	comment = "Hello."

	write!(ws8, "C3", cell_text, fmt=text_wrap)
	write_comment!(ws8, "C3", comment)

	cell_text = "The height of this row has been adjusted by Excel due to the text wrap property being set. Unfortunately this means that the height of the row is unknown to XlsxWriter at run time and thus the comment box is stretched as well.\n\nUse set_row() to specify the row height explicitly to avoid this problem."

	write!(ws8, "C6", cell_text, fmt=text_wrap)
	write_comment!(ws8, "C6", comment)

	close(wb)
	isfile("comments2.xlsx")
end

test()


text:
function remove_tasks(gamma::Float64, LS::LoadingSequence, CTS::Constants)
    new_LS = deepcopy(LS)
    mod_gamma = Int(floor(CTS.P*gamma))
    if new_LS.order[mod_gamma+1].task.p == 0 && new_LS.order[mod_gamma+1].task.c == 0
        mod_gamma -= 1
    end
    for t = 1:mod_gamma
        if new_LS.order[end].task.p != 0 && new_LS.order[end].task.c != 0
            new_LS.len -= 1
            new_LS.tasks_left += 1
        end
        deleteat!(new_LS.filled_pos, findall(x->x==new_LS.order[length(new_LS.order)].task.p, new_LS.filled_pos))
        deleteat!(new_LS.loaded_cont, findall(x->x==new_LS.order[length(new_LS.order)].task.c, new_LS.loaded_cont))
        deleteat!(new_LS.order, length(new_LS.order))
    end
    return(new_LS)
end

function last_qc_moves(LS::LoadingSequence, CTS::Constants)
    QC_MOVES=Dict{Int, Array{NamedTuple{(:bay, :end_time, :start_time, :time),Tuple{Int64,Int64,Int64,Int64}}, 1}}()
    min_q = (q=0, end_time=CTS.H, start_time=CTS.H, time=CTS.H)

    flag = Array{Int,1}()
    while Set(collect(1:CTS.Q)) != Set(flag)
        for t in reverse(LS.order)
            for q = 1:CTS.Q
                if t.qc == q && !(q in flag)
                    push!(flag, q)
                    if t.start_time + 2*t.task.t < min_q.end_time
                        min_q = (q=q, end_time=t.start_time + 2*t.task.t, start_time=t.start_time, time = t.task.t)
                        break
                    end
                end
            end
        end
    end

    flag = Array{Int,1}()
    while Set(collect(1:CTS.Q)) != Set(flag)
        for q = 1:CTS.Q
            l=Array{NamedTuple{(:bay, :end_time, :start_time, :time),Tuple{Int64,Int64,Int64,Int64}}, 1}()
            for t in reverse(LS.order)
                if t.qc == q
                    if t.start_time > min_q.end_time
                        push!(l, (bay=t.task.b, end_time=t.start_time + 2*t.task.t, start_time=t.start_time, time=t.task.t))
                    else
                        push!(l, (bay=t.task.b, end_time=t.start_time + 2*t.task.t, start_time=t.start_time, time=t.task.t))
                        push!(flag, q)
                        break
                    end
                end
            end
            QC_MOVES[q]=reverse(l)
        end
    end
    return(min_q, QC_MOVES)
end

function get_current_state(LS::LoadingSequence, CTS::Constants)
    TIME = init_timer(CTS)
    TIME.available_cranes = Array{Int, 1}()
    QC = init_qc(CTS)
    min_q, QC_MOVES = last_qc_moves(LS, CTS)
    QC_TRAVEL = get_qc_travel(QC_MOVES, CTS)

    TIME.period = min_q.end_time
    TIME.available_cranes = Array{Int, 1}()

    for q = 1:CTS.Q
        QC[q].task_buffer = Array{LTask, 1}()
        #idle crane
        if QC_MOVES[q][end].end_time == TIME.period && length(QC_MOVES[q]) == 1
            QC[q].status = "idle"
            QC[q].time_left = 0
            QC[q].current_bay = QC_MOVES[q][end].bay
            QC[q].next_bay = QC_MOVES[q][end].bay
            push!(TIME.available_cranes, q)
        #loading crane
        elseif TIME.period >= QC_MOVES[q][end].start_time
            QC[q].status = "loading"
            QC[q].time_left = QC_MOVES[q][end].end_time - TIME.period
            QC[q].current_bay = QC_MOVES[q][end].bay
            QC[q].next_bay = QC_MOVES[q][end].bay
        #moving or waiting to move/load crane
        elseif QC_MOVES[q][end-1].end_time <= TIME.period && TIME.period < QC_MOVES[q][end].start_time
            QC[q].next_bay = QC_MOVES[q][end].bay
            push!(QC[q].task_buffer, LTask(0, QC_MOVES[q][end].bay, 0, QC_MOVES[q][end].time))
            if QC_MOVES[q][end-1].bay == QC_MOVES[q][end].bay
                QC[q].status = "waiting to load"
                QC[q].time_left = QC_MOVES[q][end].start_time - TIME.period
                QC[q].current_bay = QC_MOVES[q][end].bay
            elseif QC_TRAVEL[q][1].start_time <= TIME.period
                QC[q].status = "moving"
                QC[q].time_left = QC_TRAVEL[q][1].end_time - TIME.period
                QC[q].current_bay = QC_MOVES[q][end].bay - QC[q].time_left/CTS.tt
            else
                QC[q].status = "waiting to move"
                QC[q].time_left = QC_TRAVEL[q][1].start_time - TIME.period
                QC[q].current_bay = QC_MOVES[q][end-1].bay
            end
        #still loading the previous task
        elseif  QC_MOVES[q][end-1].start_time <= TIME.period
            if length(QC_MOVES[q]) == 2
                QC[q].status = "loading"
                QC[q].time_left = QC_MOVES[q][end-1].end_time - TIME.period
                QC[q].current_bay = QC_MOVES[q][end-1].bay
                QC[q].next_bay = QC_MOVES[q][end].bay
                push!(QC[q].task_buffer, LTask(0, QC_MOVES[q][end].bay, 0, QC_MOVES[q][end].time))
            end
        else
            for i = 1:length(QC_MOVES[q])-1
                #moving or waiting to move/load crane
                if QC_MOVES[q][i].end_time <= TIME.period && TIME.period < QC_MOVES[q][i+1].start_time
                    QC[q].next_bay = QC_MOVES[q][i+1].bay
                    for i = 1:length(QC_MOVES[q])-1
                        push!(QC[q].task_buffer, LTask(0, QC_MOVES[q][i+1].bay, 0, QC_MOVES[q][i+1].time))
                    end
                    if QC_MOVES[q][i].bay == QC_MOVES[q][i+1].bay
                        QC[q].status = "waiting to load"
                        QC[q].time_left = QC_MOVES[q][i+1].start_time - TIME.period
                        QC[q].current_bay = QC_MOVES[q][i+1].bay
                    elseif QC_TRAVEL[q][1].start_time <= TIME.period
                        QC[q].status = "moving"
                        QC[q].time_left = QC_TRAVEL[q][1].end_time - TIME.period
                        QC[q].current_bay = QC_MOVES[q][i+1].bay - QC[q].time_left/CTS.tt
                    else
                        QC[q].status = "waiting to move"
                        QC[q].time_left = QC_TRAVEL[q][1].start_time - TIME.period
                        QC[q].current_bay = QC_MOVES[q][i].bay
                    end
                #still loading the previous task
                elseif  QC_MOVES[q][i].start_time <= TIME.period
                    QC[q].status = "loading"
                    QC[q].time_left = QC_MOVES[q][i].end_time - TIME.period
                    QC[q].current_bay = QC_MOVES[q][i].bay
                    QC[q].next_bay = QC_MOVES[q][i].bay
                    for i = 1:length(QC_MOVES[q])-1
                        push!(QC[q].task_buffer, LTask(0, QC_MOVES[q][i+1].bay, 0, QC_MOVES[q][i+1].time))
                    end
                end
            end
        end
    end
    return(true, TIME, QC)
end

function get_single_moves(QC_MOVES::Dict{Int, Array{NamedTuple{(:bay, :end_time, :start_time, :time),Tuple{Int64,Int64,Int64,Int64}}, 1}}, CTS::Constants)
    single_moves = Array{NamedTuple{(:qc, :bay, :start_time, :time),Tuple{Int64, Int64, Int64, Int64}},1}()
    check = false
    for q = 1:CTS.Q
        for t = 2:length(QC_MOVES[q])-1
            if QC_MOVES[q][t-1].bay < QC_MOVES[q][t].bay && QC_MOVES[q][t].bay > QC_MOVES[q][t+1].bay
                if QC_MOVES[q][t].time != 0
                    push!(single_moves, (qc=q, bay=QC_MOVES[q][t].bay, start_time = QC_MOVES[q][t].start_time, time = QC_MOVES[q][t].time))
                    check = true
                end
            elseif QC_MOVES[q][t-1].bay > QC_MOVES[q][t].bay && QC_MOVES[q][t].bay < QC_MOVES[q][t+1].bay
                if QC_MOVES[q][t].time != 0
                    push!(single_moves, (qc=q, bay=QC_MOVES[q][t].bay, start_time = QC_MOVES[q][t].start_time, time = QC_MOVES[q][t].time))
                    check = true
                end
            end
        end

        if QC_MOVES[q][end-1].bay != QC_MOVES[q][end-2].bay && QC_MOVES[q][end-1].time != 0
            push!(single_moves, (qc=q, bay=QC_MOVES[q][end-1].bay, start_time = QC_MOVES[q][end-1].start_time, time = QC_MOVES[q][end-1].time))
            check = true
        end
        if QC_MOVES[q][1].bay != init_bay(q,CTS) && QC_MOVES[q][1].bay != QC_MOVES[q][2].bay && QC_MOVES[q][2].time != 0
            push!(single_moves, (qc=q, bay=QC_MOVES[q][1].bay, start_time = QC_MOVES[q][1].start_time, time = QC_MOVES[q][1].time))
            check = true
        end
    end

    if check == false
        return(false)
    else
        return(single_moves)
    end
end

function get_single_task(single::NamedTuple{(:qc, :bay, :start_time, :time),Tuple{Int64,Int64,Int64,Int64}}, LS::LoadingSequence, CTS::Constants)
    single_task = (task=LTask(0,0,0,0), start_time=0, qc=0)
    left_flag = false
    right_flag = false
    l = Array{NamedTuple{(:task, :start_time, :qc),Tuple{LTask,Int64,Int64}},1}()
    r = Array{NamedTuple{(:task, :start_time, :qc),Tuple{LTask,Int64,Int64}},1}()
    for t in LS.order
        #single taks
        if t.start_time == single.start_time && t.qc == single.qc && t.task.b == single.bay && t.task.t == single.time
                single_task = t
        end
        if t.start_time < single.start_time && t.task.b == single.bay && t.qc == single.qc
            if left_flag == false
                left_flag = true
            elseif l[end].start_time + 2*l[end].task.t == t.start_time
                deleteat!(l, length(l))
            end
            push!(l, t)
        elseif single.start_time < t.start_time && t.task.b == single.bay && t.qc == single.qc
            if right_flag == false
                right_flag = true
            elseif r[end].start_time + 2*r[end].task.t == t.start_time
                deleteat!(r, length(r))
            end
            push!(r, t)
        end
    end

    return(single_task, left_flag, right_flag, l, r)
end

function trim_left(single::NamedTuple{(:task, :start_time, :qc),Tuple{LTask,Int64,Int64}}, left::NamedTuple{(:task, :start_time, :qc),Tuple{LTask,Int64,Int64}}, LS::LoadingSequence, CTS::Constants)
    LS1 = deepcopy(LS)
    LS2 = deepcopy(LS)
    LS3 = deepcopy(LS)

    #delete single_task
    deleteat!(LS1.filled_pos, findall(x->x==single.task.p, LS1.filled_pos))
    deleteat!(LS1.loaded_cont, findall(x->x==single.task.c, LS1.loaded_cont))
    deleteat!(LS1.order, findall(x->x==single, LS1.order))
    deleteat!(LS2.filled_pos, findall(x->x==single.task.p, LS2.filled_pos))
    deleteat!(LS2.loaded_cont, findall(x->x==single.task.c, LS2.loaded_cont))
    deleteat!(LS2.order, findall(x->x==single, LS2.order))
    deleteat!(LS3.filled_pos, findall(x->x==single.task.p, LS3.filled_pos))
    deleteat!(LS3.loaded_cont, findall(x->x==single.task.c, LS3.loaded_cont))
    deleteat!(LS3.order, findall(x->x==single, LS3.order))

    #trim LS
    for t in LS.order
        if  left.start_time < t.start_time
            deleteat!(LS1.filled_pos, findall(x->x==t.task.p, LS1.filled_pos))
            deleteat!(LS1.loaded_cont, findall(x->x==t.task.c, LS1.loaded_cont))
            deleteat!(LS1.order, findall(x->x==t, LS1.order))
            if single.start_time <= t.start_time
                deleteat!(LS2.filled_pos, findall(x->x==t.task.p, LS2.filled_pos))
                deleteat!(LS2.loaded_cont, findall(x->x==t.task.c, LS2.loaded_cont))
                deleteat!(LS2.order, findall(x->x==t, LS2.order))
            elseif t.start_time <= single.start_time + 2*single.task.t
                deleteat!(LS3.filled_pos, findall(x->x==t.task.p, LS3.filled_pos))
                deleteat!(LS3.loaded_cont, findall(x->x==t.task.c, LS3.loaded_cont))
                deleteat!(LS3.order, findall(x->x==t, LS3.order))
            end
        elseif t.start_time <= left.start_time
            deleteat!(LS2.filled_pos, findall(x->x==t.task.p, LS2.filled_pos))
            deleteat!(LS2.loaded_cont, findall(x->x==t.task.c, LS2.loaded_cont))
            deleteat!(LS2.order, findall(x->x==t, LS2.order))
            deleteat!(LS3.filled_pos, findall(x->x==t.task.p, LS3.filled_pos))
            deleteat!(LS3.loaded_cont, findall(x->x==t.task.c, LS3.loaded_cont))
            deleteat!(LS3.order, findall(x->x==t, LS3.order))
        end
    end

    return(LS1, LS2, LS3)
end

function trim_right(single::NamedTuple{(:task, :start_time, :qc),Tuple{LTask,Int64,Int64}}, right::NamedTuple{(:task, :start_time, :qc),Tuple{LTask,Int64,Int64}}, LS::LoadingSequence, CTS::Constants)
    LS1 = deepcopy(LS)
    LS2 = deepcopy(LS)
    LS3 = deepcopy(LS)

    #delete single_task
    deleteat!(LS1.filled_pos, findall(x->x==single.task.p, LS1.filled_pos))
    deleteat!(LS1.loaded_cont, findall(x->x==single.task.c, LS1.loaded_cont))
    deleteat!(LS1.order, findall(x->x==single, LS1.order))
    deleteat!(LS2.filled_pos, findall(x->x==single.task.p, LS2.filled_pos))
    deleteat!(LS2.loaded_cont, findall(x->x==single.task.c, LS2.loaded_cont))
    deleteat!(LS2.order, findall(x->x==single, LS2.order))
    deleteat!(LS3.filled_pos, findall(x->x==single.task.p, LS3.filled_pos))
    deleteat!(LS3.loaded_cont, findall(x->x==single.task.c, LS3.loaded_cont))
    deleteat!(LS3.order, findall(x->x==single, LS3.order))

    #trim LS
    for t in LS.order
        if  t.start_time >= single.start_time
            deleteat!(LS1.filled_pos, findall(x->x==t.task.p, LS1.filled_pos))
            deleteat!(LS1.loaded_cont, findall(x->x==t.task.c, LS1.loaded_cont))
            deleteat!(LS1.order, findall(x->x==t, LS1.order))
            if t.start_time >= right.start_time + 2*right.task.t
                deleteat!(LS2.filled_pos, findall(x->x==t.task.p, LS2.filled_pos))
                deleteat!(LS2.loaded_cont, findall(x->x==t.task.c, LS2.loaded_cont))
                deleteat!(LS2.order, findall(x->x==t, LS2.order))
            elseif t.start_time <= right.start_time
                deleteat!(LS3.filled_pos, findall(x->x==t.task.p, LS3.filled_pos))
                deleteat!(LS3.loaded_cont, findall(x->x==t.task.c, LS3.loaded_cont))
                deleteat!(LS3.order, findall(x->x==t, LS3.order))
            end
        elseif t.start_time < single.start_time
            deleteat!(LS2.filled_pos, findall(x->x==t.task.p, LS2.filled_pos))
            deleteat!(LS2.loaded_cont, findall(x->x==t.task.c, LS2.loaded_cont))
            deleteat!(LS2.order, findall(x->x==t, LS2.order))
            deleteat!(LS3.filled_pos, findall(x->x==t.task.p, LS3.filled_pos))
            deleteat!(LS3.loaded_cont, findall(x->x==t.task.c, LS3.loaded_cont))
            deleteat!(LS3.order, findall(x->x==t, LS3.order))
        end
    end
    return(LS1, LS2, LS3)
end

function merge_left(single::NamedTuple{(:task, :start_time, :qc),Tuple{LTask,Int64,Int64}}, LS1::LoadingSequence, LS2::LoadingSequence, LS3::LoadingSequence, CTS::Constants)
    #merge single task
    push!(LS1.order, (task=single.task, start_time = get_qc_last_time(single.qc, LS1) + travel_time(get_qc_last_bay(single.qc, LS1), single.task.b, CTS), qc=single.qc))
    push!(LS1.filled_pos, single.task.p)
    push!(LS1.loaded_cont, single.task.c)
    #update and merge LS2
    for t in LS2.order
        if t.qc == single.qc
            nt = (task=t.task, start_time = t.start_time + 2*single.task.t, qc=single.qc)
            push!(LS1.order, nt)
        else
            push!(LS1.order, t)
        end
        push!(LS1.filled_pos, t.task.p)
        push!(LS1.loaded_cont, t.task.c)
    end
    #update and merge LS3
    if length(LS3.order) > 0
    dif = 2*travel_time(get_qc_last_bay(single.qc, LS2), single.task.b, CTS)
        for t in LS3.order
            if t.qc == single.qc
                nt = (task=t.task, start_time = t.start_time - dif, qc=single.qc)
                push!(LS1.order, nt)
            else
                push!(LS1.order, t)
            end
            push!(LS1.filled_pos, t.task.p)
            push!(LS1.loaded_cont, t.task.c)
        end
    end
    return(LS1)
end

function merge_right(single::NamedTuple{(:task, :start_time, :qc),Tuple{LTask,Int64,Int64}}, LS1::LoadingSequence, LS2::LoadingSequence, LS3::LoadingSequence, CTS::Constants)
    #update and merge LS2
    if length(LS1.order) > 0
        dif = 2*travel_time(get_qc_last_bay(single.qc, LS1), single.task.b, CTS) + 2*single.task.t
    else
        dif = 2*travel_time(init_bay(single.qc, CTS), single.task.b, CTS) + 2*single.task.t
    end
    for t in LS2.order
        if t.qc == single.qc
            nt = (task=t.task, start_time = t.start_time - dif, qc=single.qc)
            push!(LS1.order, nt)
        else
            push!(LS1.order, t)
        end
        push!(LS1.filled_pos, t.task.p)
        push!(LS1.loaded_cont, t.task.c)
    end

    #merge single task
    push!(LS1.order, (task=single.task, start_time = get_qc_last_time(single.qc, LS1) + travel_time(get_qc_last_bay(single.qc, LS1), single.task.b, CTS), qc=single.qc))
    push!(LS1.filled_pos, single.task.p)
    push!(LS1.loaded_cont, single.task.c)

    #update and merge LS3
    dif = dif - 2*single.task.t
    for t in LS3.order
        if t.qc == single.qc
            nt = (task=t.task, start_time = t.start_time - dif, qc=single.qc)
            push!(LS1.order, nt)
        else
            push!(LS1.order, t)
        end
        push!(LS1.filled_pos, t.task.p)
        push!(LS1.loaded_cont, t.task.c)
    end
    return(LS1)
end

function remove_single_moves(prec::Dict{Int, Array}, makespan::Int, LS::LoadingSequence, CTS::Constants)
    QC_MOVES = get_qc_moves(LS, CTS)
    single_moves = get_single_moves(QC_MOVES, CTS)
    if single_moves == false
        return(LS)
    else
        for single in single_moves
            if single.qc == LS.order[end].qc
                single_task, left_flag, right_flag, l, r = get_single_task(single, LS, CTS)
                if left_flag == true
                    for l_task in l
                        LS1, LS2, LS3 = trim_left(single_task, l_task, LS, CTS)
                        new_LS = merge_left(single_task, LS1, LS2, LS3, CTS)
                        if check_solution(prec, new_LS, CTS) == true
                            if total_makespan(new_LS, CTS) < makespan
                                return(new_LS)
                            end
                        end
                    end
                end
                if right_flag == true
                    for r_task in r
                        LS1, LS2, LS3 = trim_right(single_task, r_task, LS, CTS)
                        new_LS = merge_right(single_task, LS1, LS2, LS3, CTS)
                        if check_solution(prec, new_LS, CTS) == true
                            if total_makespan(new_LS, CTS) < makespan
                                return(new_LS)
                            end
                        end
                    end
                end
            end
        end
        return(LS)
    end
end

function remove_useless_travel(LS::LoadingSequence)
    if LS.order[end].task.t == 0
        deleteat!(LS.order, length(LS.order))
    end
end

function get_final_state(LS::LoadingSequence, CTS::Constants)
    TIME = init_timer(CTS)
    TIME.available_cranes = Array{Int, 1}()
    QC = init_qc(CTS)
    min_q, moves, QC_MOVES = ls_qc_moves(best_LS, CTS)
    QC_TRAVEL = get_qc_travel(QC_MOVES, CTS)

    TIME.period = min_q.end_time
    TIME.available_cranes = Array{Int, 1}()

    for q = 1:CTS.Q
        QC[q].task_buffer = Array{LTask, 1}()
        #idle crane
        if QC_MOVES[q][1].end_time == TIME.period
            QC[q].status = "idle"
            QC[q].time_left = 0
            QC[q].current_bay = QC_MOVES[q][1].bay
            QC[q].next_bay = QC_MOVES[q][1].bay
            push!(TIME.available_cranes, q)
        #loading crane
        elseif QC_MOVES[q][1].start_time <= TIME.period && TIME.period < QC_MOVES[q][1].end_time
            QC[q].status = "loading"
            QC[q].time_left = QC_MOVES[q][1].end_time - TIME.period
            QC[q].current_bay = QC_MOVES[q][1].bay
            QC[q].next_bay = QC_MOVES[q][1].bay
        else
            QC[q].status = "idle"
            QC[q].time_left = 0
            QC[q].current_bay = QC_MOVES[q][1].bay
            QC[q].next_bay = QC_MOVES[q][1].bay
            push!(TIME.available_cranes, q)
        end
    end
    return(min_q.q, moves, TIME, QC)
end

function remove_final_tasks(moves::Int, LS::LoadingSequence, CTS::Constants)
    new_LS = deepcopy(LS)

    if moves > 2
        while new_LS.len >= LS.len - (moves - 2)
            if new_LS.order[end].task.p != 0 && new_LS.order[end].task.c != 0
                new_LS.len -= 1
                new_LS.tasks_left += 1
            end
            deleteat!(new_LS.filled_pos, findall(x->x==new_LS.order[length(new_LS.order)].task.p, new_LS.filled_pos))
            deleteat!(new_LS.loaded_cont, findall(x->x==new_LS.order[length(new_LS.order)].task.c, new_LS.loaded_cont))
            deleteat!(new_LS.order, length(new_LS.order))
        end
        if new_LS.order[end].task.p == 0 && new_LS.order[end].task.c == 0
            deleteat!(new_LS.order, length(new_LS.order))
        end
    end
    return(new_LS)
end

function ls_qc_moves(LS::LoadingSequence, CTS::Constants)
    QC_MOVES=Dict{Int, Array{NamedTuple{(:bay, :end_time, :start_time, :time),Tuple{Int64,Int64,Int64,Int64}}, 1}}()
    min_q = (q=0, end_time=CTS.H, start_time=CTS.H, time=CTS.H)
    moves = 0
    flag = Array{Int,1}()
    while Set(collect(1:CTS.Q)) != Set(flag)
        for t in reverse(LS.order)
            for q = 1:CTS.Q
                if t.qc == q && !(q in flag) && t.task.t != 0
                    push!(flag, q)
                    if t.start_time + 2*t.task.t < min_q.end_time
                        min_q = (q=q, end_time=t.start_time + 2*t.task.t, start_time=t.start_time, time = t.task.t)
                        break
                    end
                end
            end
        end
    end

    flag = Array{Int,1}()
    while Set(collect(1:CTS.Q)) != Set(flag)
        for q = 1:CTS.Q
            l=Array{NamedTuple{(:bay, :end_time, :start_time, :time),Tuple{Int64,Int64,Int64,Int64}}, 1}()
            for t in reverse(LS.order)
                if t.qc == q
                    if t.start_time > min_q.end_time
                        push!(l, (bay=t.task.b, end_time=t.start_time + 2*t.task.t, start_time=t.start_time, time=t.task.t))
                        if t.task.t != 0
                            moves += 1
                        end
                    elseif t.start_time < min_q.start_time && t.start_time + 2*t.task.t < min_q.start_time
                        push!(l, (bay=t.task.b, end_time=t.start_time, start_time=t.start_time, time=0))
                        push!(flag, q)
                        break
                    else
                        push!(l, (bay=t.task.b, end_time=t.start_time + 2*t.task.t, start_time=t.start_time, time=t.task.t))
                        push!(flag, q)
                        if t.task.t != 0
                            moves += 1
                        end
                        break
                    end
                end
            end
            QC_MOVES[q]=reverse(l)
        end
    end
    return(min_q, moves, QC_MOVES)
end

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
function matrix_cone(S::Type{<:MOI.AbstractSymmetricMatrixSetTriangle},
                     side_dimension)
    if iszero(side_dimension)
        # Some solvers such as Mosek does not support 0-dimensional PSD cone
        return EmptyCone()
    elseif isone(side_dimension)
        # PSD constraints on 1x1 matrices are equivalent to the nonnegativity
        # of the only entry.
        return MOI.Nonnegatives(1)
    elseif side_dimension == 2
        # PSD constraints on 2x2 matrices are SOC representable.
        return PositiveSemidefinite2x2ConeTriangle()
    else
        # PSD constraints on nxn matrices with n ≥ 3 is not SOC representable,
        # see [F18].
        #
        # [F18] Fawzi, Hamza
        # On representing the positive semidefinite cone using the second-order
        # cone.
        # Mathematical Programming (2018): 1-10.
        return S(side_dimension)
    end
end

"""
    struct DiagonallyDominantConeTriangle <: MOI.AbstractSymmetricMatrixSetTriangle
        side_dimension::Int
    end

See Definition 4 of [AM17] for a precise definition of the last two items.

[AM17] Ahmadi, A. A. & Majumdar, A.
*DSOS and SDSOS Optimization: More Tractable Alternatives to Sum of Squares and Semidefinite Optimization*
ArXiv e-prints, **2017**.
"""
struct DiagonallyDominantConeTriangle <: MOI.AbstractSymmetricMatrixSetTriangle
    side_dimension::Int
end

function matrix_cone(S::Type{DiagonallyDominantConeTriangle},
                     side_dimension)
    if iszero(side_dimension)
        return EmptyCone()
    elseif isone(side_dimension)
        return MOI.Nonnegatives(1)
    else
        # With `side_dimension` = 2, we want to avoid using SOC and only use LP
        return S(side_dimension)
    end
end

"""
    struct ScaledDiagonallyDominantConeTriangle <: MOI.AbstractSymmetricMatrixSetTriangle
        side_dimension::Int
    end

See Definition 4 of [AM17] for a precise definition of the last two items.

[AM17] Ahmadi, A. A. & Majumdar, A.
*DSOS and SDSOS Optimization: More Tractable Alternatives to Sum of Squares and Semidefinite Optimization*
ArXiv e-prints, **2017**.
"""
struct ScaledDiagonallyDominantConeTriangle <: MOI.AbstractSymmetricMatrixSetTriangle
    side_dimension::Int
end

function matrix_cone(S::Type{ScaledDiagonallyDominantConeTriangle},
                     side_dimension)
    if iszero(side_dimension)
        return EmptyCone()
    elseif isone(side_dimension)
        return MOI.Nonnegatives(side_dimension)
    elseif side_dimension == 2
        return PositiveSemidefinite2x2ConeTriangle()
    else
        return S(side_dimension)
    end
end

# isbits types, nothing to copy
function Base.copy(set::Union{DiagonallyDominantConeTriangle,
                              ScaledDiagonallyDominantConeTriangle})
    return set
end


text:
# included in module Wut

struct CoordinateEncoder

    w::Number            # width of the contiguous 1 bits, must be odd
    n::Number            # output bit width (must be greater than w)

    function CoordinateEncoder(;w=0,n=0)
        res=(maxval-minval)/n
        new(w,n)

end

end
export CoordinateEncoder

function encode(e::CoordinateEncoder, n::Number)
    encodeIntoArray(e, n, BitPat(e.n,n))
end
export encode

function getWidth(e::CoordinateEncoder)
    e.n
end
export getWidth

function getDescription(e::CoordinateEncoder)
   e.name
end
export getDescription

function getBucketIndices(e::CoordinateEncoder)

end
export getBucketIndices

function encodeIntoArray(e::CoordinateEncoder,n::Number,b::BitPat;learn=true)

    b
end
export encodeIntoArray

function decode(e::CoordinateEncoder)

end
export decode

function getBucketValues(e::CoordinateEncoder)

end
export getBucketValues

function getBucketInfo(e::CoordinateEncoder)

end
export getBucketInfo

function topDownCompute(e::CoordinateEncoder)

end
export topDownCompute

function closenessScores(e::CoordinateEncoder)

end
export closenessScores

# end # module