text:
struct FDPoisson2DProblem{TF, TBCT, TBCB, TBCL, TBCR}
    a # (a,b)×(c,d) define the domain
    b
    c
    d
    Δx # mesh spacing in x
    nx # number of points in the interior, excluding a and b
    Δy # mesh spacing in y
    ny # number of points in the interior, excluding c and d
    x # interior mesh
    y # interior mesh
    f::TF # value of f at the interior nodes
    A # the matrix
    rhs # right hand side of our system
    top_bc::TBCT
    bottom_bc::TBCB
    left_bc::TBCL
    right_bc::TBCR
end

struct DirichletBC
    u_bc
end

struct NeumannBC
    v_bc
end

struct PeriodicBC
end


text:
#=
    Code for interacting with a Graph structure by:
        - Adding nodes
        - Adding edges
        - Checking if a given edge exists
        - Retrieving neighbors of a given node.
        - Getting formatted edgelist of graph.
        - Creating a copy of a graph.
=#

"""
    Create a new Vertex in the given Graph.

    Args:
        graph:
            'Graph' structure for which we want to add a node.
"""
function add_node(graph::Graph)
    index = graph.N+1
    new_node = Vertex(index)
    push!(graph.vertices, new_node)
    graph.N += 1
end

"""
    Create an edge between vertex 'src' and vertex 'tgt'.
    Edges are created as directed edges (src->tgt) whether
    the network is directed or not.

    The direction of the network defines the behavior of
    other functions.

    Args:
        src:
            Source vertex of the created edge.
        tgt:
            Target vertex of the created edge.
        graph:
            'Graph' structure where the vertices and edges are placed.
"""
function add_edge(src, tgt, graph::Graph)
    is_directed = graph.is_directed
    node_i = graph.vertices[src]
    node_j = graph.vertices[tgt]

    # -- Create the edge object --
    new_edge = Edge(node_i.index, node_j.index)
    push!(node_i.edges_source, new_edge)
    push!(node_j.edges_target, new_edge)
    push!(graph.edges, new_edge)
    new_edge.index = length(graph.edges)
    graph.M += 1
end

"""
    Check if the edge(src,tgt) already exists in the network.
    If the network is directed, then we check for src->tgt, while
    for the undirected, src<->tgt.

    Args:
        src:
            Source vertex of the searched edge.
        tgt:
            Target vertex of the searched edge.
        graph:
        'Graph' structure where the vertices and edges are placed.
"""
function is_edge(src, tgt, graph::Graph)
    node_i = graph.vertices[src]
    node_j = graph.vertices[tgt]
    is_directed = graph.is_directed

    if is_directed
        out_edges = node_i.edges_source
        for cur_edge in out_edges
            if cur_edge.target == node_j.index
                return true
            end
        end
    else
        out_edges = node_i.edges_source
        in_edges = node_i.edges_target
        for cur_edge in out_edges
            if cur_edge.target == node_j.index
                return true
            end
        end
        for cur_edge in in_edges
            if cur_edge.source == node_j.index
                return true
            end
        end
    end
    return false
end

"""
    If the network is directed, returns outcoming neighbors.
    Otherwise, returns all neighbors.

    Args:
        v_index:
            Vertex for which we want to retrieve the neighbors.
        graph:
            'Graph' structure where the vertices and edges are placed.
"""
function get_all_neighbors_aware(v_index::Int, graph::Graph)
    node_v = graph.vertices[v_index]
    neighbors = Int[]
    out_edges_v = node_v.edges_source
    for cur_edge in out_edges_v
        append!(neighbors, [cur_edge.target])
    end
    if graph.is_directed
        return collect(Int, Set(neighbors))
    end
    in_edges_v = node_v.edges_target
    for cur_edge in in_edges_v
        append!(neighbors, [cur_edge.source])
    end
    return collect(Int, Set(neighbors))
end

# -----> COPY FUNCTIONS FOR GRAPH <----- #

"""
    Returns a (M,2) matrix representing the edgelist of the
    parsed 'graph'.

    Args:
        graph:
            'Graph' structure where the vertices and edges are placed.
"""
function get_edgelist(graph::Graph)
    M = graph.M
    edgelist = zeros(Int, (M,2))
    for (j, edge) in enumerate(graph.edges)
        edgelist[j,1] = edge.source
        edgelist[j,2] = edge.target
    end
    return edgelist
end

function copy_graph(graph::Graph)
    edges = get_edgelist(graph)
    new_graph = graph_from_edgelist(edges, graph.is_directed)

    new_graph.int_vproperties = copy(graph.int_vproperties)
    new_graph.float_vproperties = copy(graph.float_vproperties)
    new_graph.string_vproperties = copy(graph.string_vproperties)
    new_graph.int_eproperties = copy(graph.int_eproperties)
    new_graph.float_eproperties = copy(graph.float_eproperties)
    new_graph.string_eproperties = copy(graph.string_eproperties)
    return new_graph
end

ext: jl


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")

p = [0.0, 0.0]

mergemapper = CutCellDG.MergeMapper()

south = CutCellDG.cell_map_to_south()
@test allapprox(south(p), [0.0, -2.0])
@test allapprox(mergemapper[1](p), [0.0, -2.0])

southeast = CutCellDG.cell_map_to_south_east()
@test allapprox(southeast(p), [2.0, -2.0])
@test allapprox(mergemapper[5](p), [2.0, -2.0])

east = CutCellDG.cell_map_to_east()
@test allapprox(east(p), [2.0, 0.0])
@test allapprox(mergemapper[2](p), [2.0, 0.0])

northeast = CutCellDG.cell_map_to_north_east()
@test allapprox(northeast(p), [2.0, 2.0])
@test allapprox(mergemapper[6](p), [2.0, 2.0])

north = CutCellDG.cell_map_to_north()
@test allapprox(north(p), [0.0, 2.0])
@test allapprox(mergemapper[3](p), [0.0, 2.0])

northwest = CutCellDG.cell_map_to_north_west()
@test allapprox(northwest(p), [-2.0, 2.0])
@test allapprox(mergemapper[7](p), [-2.0, 2.0])

west = CutCellDG.cell_map_to_west()
@test allapprox(west(p), [-2.0, 0.0])
@test allapprox(mergemapper[4](p), [-2.0, 0.0])

southwest = CutCellDG.cell_map_to_south_west()
@test allapprox(southwest(p), [-2.0, -2.0])
@test allapprox(mergemapper[8](p), [-2.0, -2.0])

polyorder = 1
numqp = 2

elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
levelsetbasis = HermiteTensorProductBasis(2)

x0 = [0.0, 0.0]
meshwidths = [2.0, 1.0]
nelements = [2, 1]
refpoints = interpolation_points(elasticitybasis)
dim,nf = size(interpolation_points(levelsetbasis))
dgmesh = CutCellDG.DGMesh(x0, meshwidths, nelements, refpoints)
cgmesh = CutCellDG.CGMesh(x0, meshwidths, nelements, nf)

normal = [1.0, 0.0]
xI = [1.1, 0.0]
levelset = CutCellDG.LevelSet(
    x -> plane_distance_function(x, normal, xI),
    cgmesh,
    levelsetbasis,
)

cutmesh = CutCellDG.CutMesh(dgmesh, levelset)
cellquads = CutCellDG.CellQuadratures(cutmesh, levelset, numqp)
interfacequads = CutCellDG.InterfaceQuadratures(cutmesh, levelset, numqp)
facequads = CutCellDG.FaceQuadratures(cutmesh, levelset, numqp)

mergedmesh = CutCellDG.MergedMesh!(cutmesh, cellquads, facequads, interfacequads)
@test CutCellDG.number_of_nodes(mergedmesh) == 8
@test CutCellDG.has_merged_cells(mergedmesh)
@test CutCellDG.merge_direction(mergedmesh,+1,1) == 0
@test CutCellDG.merge_direction(mergedmesh,+1,2) == 0
@test CutCellDG.merge_direction(mergedmesh,-1,1) == 0
@test CutCellDG.merge_direction(mergedmesh,-1,2) == 2


text:
# Creates the sampling function needed for ApproximateGPs.pathwise_sample
function create_prior_sample_function(num_basis_features=100)
    function prior_sample_function(rng, prior, input_dims, num_samples)
        ϕ = sample_basis(rng, prior.kernel, num_basis_features, input_dims)
        w = rand(rng, Normal(), (num_basis_features, num_samples))
        return ϕ, w
    end
    return prior_sample_function
end

## Now in ApproximateGPs
# struct PosteriorSample
#     prior  # prior(x) returns an approx prior sample
#     update  # update(x) computes the pathwise update
# end

# (f::PosteriorSample)(x) = f([x])
# (f::PosteriorSample)(x::AbstractVector) = f.prior(x) + f.update(x)

# function sample_posterior_functions(f::ApproxPosteriorGP{SVGP}, n_features=10, n_samples=1)
#     dims = length(f.data.u[1])

#     # 1. sample a prior basis
#     ϕ = prior_basis(f.prior.kernel, n_features, dims)

#     # 2. sample a prior function
#     prior, w = sample_prior_functions(ϕ, n_features, n_samples)

#     # 3. pathwise update
#     z = f.data.u
#     Pw = w'hcat(ϕ.(z)...)
#     q = f.approx.q
#     u = rand(q, n_samples)
#     v = f.data.Kuu \ (u - Pw')
#     update(x) = v'cov(f.prior, z, x)

#     return PosteriorSample(prior, update)
# end


text:
import Base: show, getindex

abstract type AggregationFunction end

struct Aggregation{N} <: AggregationFunction
    fs::NTuple{N, AggregationFunction}
    Aggregation(fs::Vararg{AggregationFunction, N}) where N = new{N}(fs)
    Aggregation(fs::NTuple{N, AggregationFunction}) where N = new{N}(fs)
end

Flux.@functor Aggregation

(a::Aggregation)(x::Union{AbstractArray, Missing}, args...) = vcat([f(x, args...) for f in a.fs]...)
(a::AggregationFunction)(x::ArrayNode, args...) = mapdata(x -> a(x, args...), x)

Base.getindex(a::Aggregation, i) = a.fs[i]

const AggregationWeights = Union{Nothing,
                                 AbstractVector{T} where T <: Real,
                                 AbstractMatrix{T} where T <: Real}
const MaybeMatrix = Union{Missing,
                          AbstractMatrix{T} where T <: Real}

bagnorm(w::Nothing, b) = length(b)
bagnorm(w::AbstractVector, b) = @views sum(w[b])
bagnorm(w::AbstractMatrix, b) = @views vec(sum(w[:, b], dims=2))

weight(w::Nothing, _, _) = 1
weight(w::AbstractVector, _, j) = w[j]
weight(w::AbstractMatrix, i, j) = w[i, j]

weightsum(ws::Real, _) = ws
weightsum(ws::AbstractVector, i) = ws[i]

# more stable definitions for r_map and p_map
Zygote.@adjoint softplus(x) = softplus.(x), Δ -> (Δ .* σ.(x),)

include("segmented_sum.jl")
include("segmented_mean.jl")
include("segmented_max.jl")
include("segmented_pnorm.jl")
include("segmented_lse.jl")
# include("transformer.jl")

Base.show(io::IO, ::MIME"text/plain", a::T) where T <: AggregationFunction = print(io, "$(T.name)($(length(a.C)))")
Base.show(io::IO, m::MIME"text/plain", a::Aggregation{N}) where N = print(io, "⟨" * join(repr(m, f) for f in a.fs ", ") * "⟩")

export SegmentedSum, SegmentedMean, SegmentedMax, SegmentedPNorm, SegmentedLSE

const names = ["Sum", "Mean", "Max", "PNorm", "LSE"]
for idxs in powerset(collect(1:length(names)))
    length(idxs) > 1 || continue
    for p in permutations(idxs)
        s = Symbol("Segmented", names[p]...)
        @eval function $s(d::Int)
            Aggregation($((Expr(:call, Symbol("Segmented" * n), :d)
                           for n in names[p])...))
        end
        @eval function $s(D::Vararg{Int, $(length(p))})
            Aggregation($((Expr(:call, Symbol("Segmented" * n), :(D[$i]))
                           for (i,n) in enumerate(names[p]))...))
        end
        @eval export $s
    end
end


text:
export connectnodes

"""
    connectnodes(nodes::Array{Node,1}, [parameters::SegmentParameters])

Returns an array of `Segment`s connecting `Node`s.
"""
function connectnodes(nodes::Array{Node,1},sp::SegmentParameters = SegmentParameters())
  segments = Array{Segment}(undef,length(nodes)-1)
  for i in eachindex(segments)
    segments[i] = Segment(nodes[i],nodes[i+1],sp)
  end
  return segments
end


text:
using Pkg
Pkg.activate(@__DIR__)
CI = get(ENV, "CI", nothing) == "true"
ENV["PYTHON"] = ""
Pkg.build("PyCall")
using Documenter, ClimateTools

makedocs(sitename = "ClimateTools.jl",
    doctest = false,
    format = Documenter.HTML(
    prettyurls = CI,
    ),
    pages = [
       "index.md",
       "installation.md",
       "datasets.md",
       "indices.md",
       "interpolation.md",
       "biascorrection.md",
       "maps.md",
       "interface.md",
       "examples.md",
       "functions.md",
       ]

)

if CI
    deploydocs(
    repo   = "github.com/JuliaClimate/ClimateTools.jl.git",
    target = "build"
    )
end


text:
# using LinearAlgebra, Plots, JSON
# include("../src/SFFM.jl")
include("METools.jl")
T = [-2.0 2.0; 1.0 -1.0]#[-2.0 2.0 0; 1.0 -2.0 1; 1 1 -2]
C = [1.0; -2.0]#; -1]
fn(x) = [ones(size(x)) ones(size(x))]# ones(size(x))]
model = SFFM.Model(;T=T,C=C,r=(r=fn,R=fn),Bounds=[0 10;-Inf Inf])
N₋ = sum(C.<=0)
N₊ = sum(C.>=0)
NPhases = length(C)

t = 4.0
τ = SFFM.FixedTime(T=t)
NSim = 100_000
sims = SFFM.SimSFM(model=model,StoppingTime=τ,InitCondition=(φ=2*ones(Int,NSim),X=zeros(NSim)))

function MakeGlobalApprox(;NCells = 3,up, down,T,C,bkwd=false,D=[],plusI = false)
    αup,Qup,qup = up
    αdown,Qdown,qdown = down
    N₋ = sum(C.<=0)
    N₊ = sum(C.>=0)
    NPhases = length(C)
    NBases = length(αup)
    Q = zeros(NCells*NBases*NPhases,NCells*NBases*NPhases)
    for n in 1:NCells
        for i in 1:NPhases
            idx = (1:NBases) .+ (n-1)*(NBases*NPhases) .+ (i-1)*NBases
            if C[i]>0
                Q[idx,idx] = Qup*abs(C[i])
                if n<NCells
                    Q[idx,idx .+ NBases*NPhases] = qup*αup*abs(C[i])
                end
            elseif C[i]<0
                Q[idx,idx] = Qdown*abs(C[i])
                if n>1
                    Q[idx,idx .- NBases*NPhases] = qdown*αdown*abs(C[i])
                end
            end
        end
    end
    T₋₋ = T[C.<=0,C.<=0]
    T₊₋ = T[C.>=0,:].*((C.<0)')
    T₋₊ = T[C.<=0,:].*((C.>0)')
    T₊₊ = T[C.>=0,C.>=0]

    inLower = [kron(diagm(abs.(C).*(C.<0)),qdown)[:,C.<0]; zeros((NCells-1)*NPhases*NBases,N₋)]
    outLower = [kron(1,kron(T₋₊,αup)) zeros(N₋,N₊+(NCells-1)*NPhases*NBases)]
    inUpper = [zeros((NCells-1)*NPhases*NBases,N₊);kron(diagm(abs.(C).*(C.>0)),qup)[:,C.>0]]
    outUpper = [zeros(N₊,N₋+(NCells-1)*NPhases*NBases) kron(1,kron(T₊₋,αdown))]

    Tdiag = diagm(diag(T))
    Toff = T-diagm(diag(T))

    G = kron(I(NCells),kron(Tdiag,I(NPhases*NBases÷2))) +
        kron(I(NCells),kron(kron(Toff,I(2)[end:-1:1,:]),I(NBases÷2)))
        display(size(G))
        display(size(Q))
    B = [
        T₋₋ outLower;
        inLower G+Q inUpper;
        outUpper T₊₊;
    ]

    return Q, B
end

let
    globalerrME = []
    globalerrDG = []
    BasesVec = 1:2:21
    for NBases in BasesVec
        Δ = 1
        # NBases = 2
        Nodes = collect(0:Δ:10)
        NCells = length(Nodes)-1
        Erlang = MakeErlang(NBases,mean=Δ)#MakeME(CMEParams[NBases], mean = Δ)#
        AugQ = [Erlang.Q 0*Erlang.Q;
            0*Erlang.Q Erlang.Q']
        # AugQ = [Erlang.Q 0*Erlang.Q;
        #     0*Erlang.Q (Erlang.Q-diagm(diag(Erlang.Q)))']
        # AugQ = [Erlang.Q 0*Erlang.Q;
        #     0*Erlang.Q -Erlang.Q]
        # AugQ[NBases+1,NBases+1] = 0
        rev = reversal(Erlang)
        AugQ = [Erlang.Q 0*Erlang.Q;
            0*Erlang.Q rev.Q+rev.q*rev.α]
        # Augα = [Erlang.α 0*Erlang.α*exp(Erlang.Q*Δ)./sum(Erlang.α*exp(Erlang.Q*Δ))]
        Augα = [Erlang.α rev.α]
        Augmented = (
            α = Augα,
            Q = AugQ,
            q = [-sum(Erlang.Q,dims=2); zeros(size(Erlang.Q,1))]
            )

        Q, B = MakeGlobalApprox(
            NCells = NCells,
            up = Augmented,
            down = Augmented,
            T = T,
            C = C,
        )
        if NBases<6
            display(Erlang.Q)
            display(B)
        end
        # display(Q)

        DGMesh = SFFM.MakeMesh(model=model,NBases=1,Nodes=collect(Nodes[1]:Δ/NBases:Nodes[end]),Basis="lagrange")
        All = SFFM.MakeAll(model=model,mesh=DGMesh)

        initDist = zeros(1,size(All.B.B,1))
        initDist[1] = 1

        temp = initDist*exp(Matrix(All.B.B)*t)#SFFM.EulerDG(D=All.B.B,y=t,x0=initDist)#
        DGdist_t = SFFM.Coeffs2Dist(model=model,mesh=DGMesh,Coeffs=temp,type="probability")

        initDist = zeros(1,size(B,1))
        initDist[1] = 1
        dist_t = initDist*exp(B*t)#SFFM.EulerDG(D=B,y=t,x0=initDist)#
        display(dist_t)
        display(DGdist_t.distribution)
        pm_t = dist_t[[1:N₋;(end-N₊+1):end]]
        manipulate(dist_t) = begin
            temp = zeros(NCells*NPhases)
            for n in 1:NBases
                temp += dist_t[(N₋+n):(2*NBases):(end-N₊)]
            end
            dist_t = reshape(temp,NPhases,NCells)
            return dist_t
        end
        dist_t = manipulate(dist_t)

        mesh = SFFM.MakeMesh(model=model,NBases=NBases,Nodes=Nodes,Basis="lagrange")
        simDist = SFFM.Sims2Dist(model=model,mesh=mesh,sims=sims,type="probability")

        localerrME = sum(abs.(pm_t-simDist.pm))
        localerrDG = sum(abs.(DGdist_t.pm-simDist.pm))

        q = plot(layout = (2,1))
        for n in 1:NCells
            x = simDist.x[n]
            for i in 1:NPhases
                if n==1
                    label1 = "DG"
                    label2 = "ME"
                else
                    label1 = false
                    label2 = false
                end
                yvalsME = [sum(dist_t[i,n])]
                yvalsDG = [sum(DGdist_t.distribution[:,(1:NBases).+NBases*(n-1),i])]
                scatter!([x],yvalsDG,label=label1,subplot=i,color=:blue,markershape=:rtriangle)
                scatter!([x],yvalsME,label=label2,subplot=i,color=:black,markershape=:ltriangle)
                localerrME += abs(sum(yvalsME-simDist.distribution[:,n,i]))
                localerrDG += abs(sum(yvalsDG-simDist.distribution[:,n,i]))
            end
        end
        display(q)
        push!(globalerrME,localerrME)
        push!(globalerrDG,localerrDG)
        # for i in 1:NPhases
        #     scatter!(simDist.x,simDist.distribution[:,:,i][:],subplot=i,label="sim")
        # end
        # display(plot!())
    end
    p = plot(BasesVec,log.(globalerrME),label=false,linestyle=:dash,colour=1)
    plot!(p,BasesVec,log.(globalerrDG),label=false,linestyle=:solid,colour=3)
    display(p)
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
module AlignUtils

using Distributed
using DataFrames
using Gadfly
using Colors
using JLD
# using Levenshtein
using ..Levenshtein # TODO: use the above line when the Levenshtein package is fixed.
using SynapseClient
using ..BamReader
using ..DISSEQT
using ..SynapseTools

export
	Sample,
	find_samples,
	find_aligned,
	assign_reference!,
	align_sample!,
	align_samples!,
	askforconfirmation,
	makecleanfolder,
	getreferenceinfo,
	reference_sanity_check,
	computecodonfrequencies,
	computenucleotidefrequencies,
	uploadaligned,
	uploadswarms,
	coverageplots

include("log.jl")
include("sample.jl")
include("synapseutils.jl")
include("align.jl")
include("misc.jl")
include("consensus.jl")
include("codonfrequency.jl")
include("nucleotidefrequency.jl")
include("readcoverage.jl")

end


text:
using PseudoArcLengthContinuation
const PALC  =  PseudoArcLengthContinuation

function test_newton()
	println("--> Test Newton method")
	N = 100
	x0 = ones(N) .+ rand(N) * 0.1
	F(x) = x.^3 .- 1.0
	Jac(x) = diagm(0 => 3.0 * x.^2)

	opts = PALC.NewtonPar(verbose = false, maxIter = 8)
	sol, hist, flag, _ = @time PALC.newton(F, Jac, x0, opts)
	sol, hist, flag, _ = @time PALC.newton(F, Jac, x0, opts, normN = x->norm(x,Inf64))

end

test_newton()


text:
using GeneralizedGenerated
using Random: GLOBAL_RNG
using TupleVectors: chainvec

export rand
EmptyNTtype = NamedTuple{(),Tuple{}} where T<:Tuple

function Base.rand(rng::AbstractRNG, d::ConditionalModel, N::Int)
    r = chainvec(rand(rng, d), N)
    for j in 2:N
        @inbounds r[j] = rand(rng, d)
    end
    return r
end

Base.rand(d::ConditionalModel, N::Int) = rand(GLOBAL_RNG, d, N)

@inline function Base.rand(rng::AbstractRNG, c::ConditionalModel)
    m = Model(c)
    return _rand(getmoduletypencoding(m), m, argvals(c))(rng)
end

@inline function Base.rand(m::ConditionalModel)
    rand(GLOBAL_RNG, m)
end

@inline function Base.rand(rng::AbstractRNG, m::Model)
    return _rand(getmoduletypencoding(m), m, NamedTuple())(rng)
end

rand(m::Model) = rand(GLOBAL_RNG, m)

sourceRand(m::Model) = sourceRand()(m)
sourceRand(jd::ConditionalModel) = sourceRand(jd.model)

export sourceRand
function sourceRand()
    function(_m::Model)
        proc(_m, st::Assign)  = :($(st.x) = $(st.rhs))
        proc(_m, st::Sample)  = :($(st.x) = rand(_rng, $(st.rhs)))
        proc(_m, st::Return)  = :(return $(st.rhs))
        proc(_m, st::LineNumber) = nothing

        vals = map(x -> Expr(:(=), x,x),parameters(_m))

        wrap(kernel) = @q begin
            _rng -> begin
                $kernel
                $(Expr(:tuple, vals...))
            end
        end

        buildSource(_m, proc, wrap) |> MacroTools.flatten
    end
end

@gg function _rand(M::Type{<:TypeLevel}, _m::Model, _args)
    body = type2model(_m) |> sourceRand() |> loadvals(_args, NamedTuple())
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end

@gg function _rand(M::Type{<:TypeLevel}, _m::Model, _args::NamedTuple{()})
    body = type2model(_m) |> sourceRand()
    @under_global from_type(_unwrap_type(M)) @q let M
        $body
    end
end


text:
mutable struct Table
    rowheights::Array{T, 1} where T <: Real
    colwidths::Array{T, 1} where T <: Real
    nrows::Int
    ncols::Int
    currentrow::Int
    currentcol::Int
    leftmargin::Real
    topmargin::Real
    center::Point
end

"""
    t = Table(nrows, ncols)
    t = Table(nrows, ncols, colwidth, rowheight)
    t = Table(rowheights, columnwidths)

Tables are centered at `O`, but you can supply a point after the specifications.

    t = Table(nrows, ncols, centerpoint)
    t = Table(nrows, ncols, colwidth, rowheight, centerpoint)
    t = Table(rowheights, columnwidths, centerpoint)

Examples

Simple tables

    t = Table(4, 3) # 4 rows and 3 cols, default is 100w, 50 h
    t = Table(4, 3, 80, 30)   # 4 rows of 30pts high, 3 cols of 80pts wide
    t = Table(4, 3, (80, 30)) # same
    t = Table((4, 3), (80, 30)) # same

Specify row heights and column widths instead of quantities:

    t = Table([60, 40, 100], 50) # 3 different height rows, 1 column 50 wide
    t = Table([60, 40, 100], [100, 60, 40]) # 3 rows, 3 columns
    t = Table(fill(30, (10)), [50, 50, 50]) # 10 rows 30 high, 3 columns 10 wide
    t = Table(50, [60, 60, 60]) # just 1 row (50 high), 3 columns 60 wide
    t = Table([50], [50]) # just 1 row, 1 column, both 50 units wide
    t = Table(50, 50, 10, 5) # 50 rows, 50 columns, 10 units wide, 5 units high
    t = Table([6, 11, 16, 21, 26, 31, 36, 41, 46], [6, 11, 16, 21, 26, 31, 36, 41, 46])
    t = Table(15:5:55, vcat(5:2:15, 15:-2:5))
     #  table has 108 cells, with:
     #  row heights: 15 20 25 30 35 40 45 50 55
     #  col widths:  5 7 9 11 13 15 15 13 11 9 7 5
    t = Table(vcat(5:10:60, 60:-10:5), vcat(5:10:60, 60:-10:5))
    t = Table(vcat(5:10:60, 60:-10:5), 50) # 1 column 50 units wide
    t = Table(vcat(5:10:60, 60:-10:5), 1:5:50)

A Table is an iterator that, for each iteration, returns a tuple of:

- the `x`/`y` point of the center of cells arranged in rows and columns (relative to current 0/0)

- the number of the cell (left to right, then top to bottom)

`nrows`/`ncols` are the number of rows and columns required.

It's sometimes useful to know which row and column you're currently on while iterating:

```
t.currentrow
t.currentcol
```

and row heights and column widths are available in:

```
t.rowheights
t.colwidths
```

`box(t::Table, r, c)` can be used to fill table cells:

```
@svg begin
    for (pt, n) in (t = Table(8, 3, 30, 15))
        randomhue()
        box(t, t.currentrow, t.currentcol, :fill)
        sethue("white")
        text(string(n), pt)
    end
end
```

or without iteration, using cellnumber:

```
@svg begin
    t = Table(8, 3, 30, 15)
    for n in eachindex(t)
        randomhue()
        box(t, n, :fill)
        sethue("white")
        text(string(n), t[n])
    end
end
```

To use a Table to make grid points:

```
julia> first.(collect(Table(10, 6)))
60-element Array{Luxor.Point,1}:
 Luxor.Point(-10.0, -18.0)
 Luxor.Point(-6.0, -18.0)
 Luxor.Point(-2.0, -18.0)
 ⋮
 Luxor.Point(2.0, 18.0)
 Luxor.Point(6.0, 18.0)
 Luxor.Point(10.0, 18.0)
```

which returns an array of points that are the center points of the cells in the table.
"""
function Table(nrows::Int, ncols::Int, cellwidth, cellheight, center=O)
    currentrow = 1
    currentcol = 1
    rowheights = fill(cellheight, nrows)
    colwidths  = fill(cellwidth, ncols)
    leftmargin = center.x - sum(colwidths)/2
    topmargin  = center.y - sum(rowheights)/2
    Table(rowheights, colwidths, nrows, ncols, currentrow, currentcol, leftmargin, topmargin, center)
end

# simple: nrows, ncols, (width, height)
Table(nrows::Int, ncols::Int, wh::NTuple{2, T}, center=O) where T <: Real =
    Table(nrows, ncols, wh[1], wh[2], center)

# simple: (nrows, ncols)
Table(rc::NTuple{2, T}, center=O) where T <: Int =
    Table(rc[1], rc[2], center)

# simple: (nrows, ncols), (width, height)
Table(rc::NTuple{2, T1}, wh::NTuple{2, T2}, center=O) where T1 <: Int where T2 <: Real =
    Table(rc[1], rc[2], wh[1], wh[2], center)

# simple: nrows, ncols using default of 100w 50 h
function Table(nrows::Int, ncols::Int, center=O)
    cellwidth  = 100
    cellheight = 50
    Table(nrows, ncols, cellwidth, cellheight, center)
end

# row heights in array, column widths in array
function Table(rowheights::Array{T1, 1}, colwidths::Array{T2, 1}, center=O) where T1 <: Real where T2 <: Real
    rowheights = collect(Iterators.flatten(rowheights))
    colwidths =  collect(Iterators.flatten(colwidths))
    nrows = length(rowheights)
    ncols = length(colwidths)
    currentrow = 1
    currentcol = 1
    leftmargin = center.x - sum(colwidths)/2
    topmargin  = center.y - sum(rowheights)/2
    Table(rowheights, colwidths, nrows, ncols, currentrow, currentcol, leftmargin, topmargin, center)
end

# row heights in array, a single column width
Table(rowheights::Array{T1, 1}, colwidth::T2, center=O) where T1 <: Real where T2 <: Real =
    Table(rowheights, [colwidth], center)

# a single row height, column widths in an array
Table(row_height::T1, colwidths::Array{T2, 1}, center=O) where T1 <: Real where T2 <: Real =
    Table([row_height], colwidths, center)

# a range of row heights, a range of column widths
Table(rowheights::AbstractRange{T1}, colwidths::AbstractRange{T2}, center=O) where T1 <: Real where T2 <: Real =
    Table(collect(rowheights), collect(colwidths), center)

# an array of row heights, a range of column widths
Table(rowheights::Array{T1, 1}, colwidths::AbstractRange{T2}, center=O) where T1 <: Real where T2 <: Real =
    Table(rowheights, collect(colwidths), center)

# a range of row heights, an array of column widths
Table(rowheights::AbstractRange{T1}, colwidths::Array{T2, 1}, center=O) where T1 <: Real where T2 <: Real =
    Table(collect(rowheights), colwidths, center)

# interfaces

function Base.iterate(t::Table)
    x = t.leftmargin + (t.colwidths[1]/2)
    y = t.topmargin  + (t.rowheights[1]/2)
    cellnumber = 2
    t.currentrow = min(div(cellnumber - 1, t.ncols) + 1, 1)
    t.currentcol = mod1(cellnumber - 1, t.ncols)
    x1 = t.leftmargin + sum(t.colwidths[1:t.currentcol - 1]) + t.colwidths[t.currentcol]/2
    y1 = t.topmargin  + sum(t.rowheights[1:t.currentrow - 1]) + t.rowheights[t.currentrow]/2
    nextpoint = Point(x1, y1)
    return ((Point(x, y), 1), (nextpoint, 2))
end

function Base.iterate(t::Table, state)
    if state[2] > t.nrows * t.ncols
        return
    end
    # state[1] is the Point
    x = state[1].x
    y = state[1].y
    # state[2] is the cellnumber
    cellnumber = state[2]
    # next pos
    t.currentrow = div(cellnumber - 1, t.ncols) + 1
    t.currentcol = mod1(cellnumber, t.ncols)
    x1 = t.leftmargin + sum(t.colwidths[1:t.currentcol - 1]) + t.colwidths[t.currentcol]/2
    y1 = t.topmargin  + sum(t.rowheights[1:t.currentrow - 1]) + t.rowheights[t.currentrow]/2
    nextpoint = Point(x1, y1)
    return ((nextpoint, cellnumber), (nextpoint, cellnumber + 1))
end

function Base.size(t::Table)
    return (t.nrows, t.ncols)
end

function Base.size(t::Table, n::Int)
    return n == 1 ? t.nrows :
           n == 2 ? t.ncols : 1
end

function Base.length(t::Table)
    t.nrows * t.ncols
end

Base.lastindex(t::Table) = length(t)

Base.eltype(::Type{Table}) = Tuple

function Base.getindex(t::Table, i::Int)
    1 <= i <= t.ncols * t.nrows || throw(BoundsError(t, i))
    currentrow = div(i - 1, t.ncols) + 1
    currentcol = mod1(i, t.ncols)
    x1 = t.leftmargin + sum(t.colwidths[1:currentcol - 1]) + t.colwidths[currentcol]/2
    y1 = t.topmargin  + sum(t.rowheights[1:currentrow - 1]) + t.rowheights[currentrow]/2
    return  Point(x1, y1)
end

# t[r, c] -> Luxor.Point(0.0, -192.5)
function Base.getindex(t::Table, r::Int, c::Int)
    r <= t.nrows || throw(BoundsError(t, r))
    c <= t.ncols || throw(BoundsError(t, c))
    x1 = t.leftmargin + sum(t.colwidths[1:c - 1]) + t.colwidths[c]/2
    y1 = t.topmargin  + sum(t.rowheights[1:r - 1]) + t.rowheights[r]/2
    return  Point(x1, y1)
end

Base.getindex(t::Table, I) = [t[i] for i in I]

# get row:    t[1, :]
Base.getindex(t::Table, r::T, ::Colon) where T <: Integer = [t[r, n] for n in 1:t.ncols]

# get column: t[:, 3]
Base.getindex(t::Table, ::Colon, c::T) where T <: Integer = [t[n, c] for n in 1:t.nrows]

Base.eachindex(t::Table) = 1:length(t)

# box extensions
"""
    box(t::Table, r::Integer, c::Integer, action::Symbol)
    box(t::Table, r::Integer, c::Integer; action=:none)

Draw a box in table `t` at row `r` and column `c`.
"""
function box(t::Table, r::Integer, c::Integer, action::Symbol; kwargs...)
    cellw, cellh = t.colwidths[c], t.rowheights[r]
    box(t[r, c], cellw, cellh, action; kwargs...)
end

box(t::Table, r::Integer, c::Integer; action=:none, reversepath=false, vertices=false) =
    box(t, r, c, action, reversepath=reversepath, vertices=vertices)

"""
    box(t::Table, cellnumber::Int, action::Symbol=:none; vertices=false)
    box(t::Table, cellnumber::Int; action=:none, vertices=false)

Make box around cell `cellnumber` in table `t`.
"""
function box(t::Table, cellnumber::Int, action::Symbol; vertices=false)
    r = div(cellnumber - 1, t.ncols) + 1
    c = mod1(cellnumber, t.ncols)

    cellw, cellh = t.colwidths[c], t.rowheights[r]
    box(t[r, c], cellw, cellh, action; vertices=vertices)
end

box(t::Table, cellnumber::Int; action=:none, vertices=false) = box(t, cellnumber, action, vertices=vertices)

"""
    highlightcells(t::Table, cellnumbers, action::Symbol=:stroke;
            color::Colorant=colorant"red",
            offset = 0)

Highlight (draw or fill) one or more cells of table `t`. `cellnumbers` is a range,
array, or an array of row/column tuples.

    highlightcells(t, 1:10, :fill, color=colorant"blue")
    highlightcells(t, vcat(1:5, 150), :stroke, color=colorant"magenta")
    highlightcells(t, [(4, 5), (3, 6)])
"""
function highlightcells(t::Table, cellnumbers, action::Symbol=:stroke;
        color::Colorant=colorant"red",
        offset = 0)
    sethue(color)
    for cell in cellnumbers
        if isa(cell, Tuple)
            row, col = cell
            box(t[row, col], t.colwidths[col] + offset, t.rowheights[row] + offset, action)
        else
            ci = CartesianIndices((t.ncols, t.nrows))[cell]
            col = ci.I[1]
            row = ci.I[2]
            box(t[cell], t.colwidths[col] + offset, t.rowheights[row] + offset, action)
        end
    end
end


text:
using LimitOrderBook
using Profile
using Printf
using DataStructures

function generate_orderbook(max_price)
    ob = OrderBook(UInt128(max_price), UInt128(1))
end

#Profile.clear()
#@profile generate_orderbook(UInt128(100000))
#Profile.print(format=:flat)

for max_price in [1e2, 1e6, 1e9]
    @printf "Creates orderbook with max price:%.0E \n" max_price
    @time generate_orderbook(max_price)
end

function send_passive_limit_orders(iterations)
    max_price = UInt128(100000)
    min_price = UInt128(1)
    ob = OrderBook(max_price, min_price)
    t = time()
    for i in 1:iterations
        for price in min_price + 1:max_price - 1
            limitorder!(ob, BUY, price, UInt128(1), UInt128(1))
        end
    end
    @printf "Passive limit orders per second:%.2E \n" (iterations*max_price)/(time() - t)
end

println("1 iterations")
@time send_passive_limit_orders(1)

println("10 iterations")
@time send_passive_limit_orders(10)

#Profile.clear()
#@profile send_passive_limit_orders()
#Profile.print(format=:flat)

function cancel_orders(iterations)
    max_price = UInt128(100000)
    min_price = UInt128(1)
    ob = OrderBook(max_price, min_price)
    order_ids = Deque{UInt128}()
    for i in 1:iterations
        for price in min_price + 1:max_price - 1
            trades, order_id = limitorder!(ob, BUY, price, UInt128(1), UInt128(1))
            push!(order_ids, order_id)
        end
    end

    n_orders = order_ids.len
    t = time()
    for order_id in order_ids
        cancel!(ob, order_id)
    end
    @printf "Canceled %.0E orders per second \n" n_orders/(time() - t)

end

println("Cancels orders: 1 iteration")
@time cancel_orders(1)

println("Cancels orders: 10 iteration")
@time cancel_orders(10)

function match_orders(iterations)
    max_price = UInt128(100000)
    min_price = UInt128(1)
    ob = OrderBook(max_price, min_price)
    order_ids = Deque{UInt128}()
    size = UInt128(1)
    n_orders = UInt128(1)
    for i in 1:iterations
        for price in min_price + 1:max_price - 1
            trades, order_id = limitorder!(ob, BUY, price, size, UInt128(1))
            n_orders = n_orders + 1
        end
    end

    # starts selling
    t = time()
    for i in 1:n_orders - 2
        limitorder!(ob, SELL, min_price, size, UInt128(1))
    end
    @printf "Matched %.3E orders per second \n" n_orders/(time() - t)
end

println("Matches orders: 1 iteration")
@time match_orders(1)

println("Matches orders: 10 iteration")
@time match_orders(10)

function match_large_order(iterations)
    max_price = UInt128(100000)
    min_price = UInt128(1)
    ob = OrderBook(max_price, min_price)
    order_ids = Deque{UInt128}()
    size = UInt128(1)
    n_orders = UInt128(1)
    for i in 1:iterations
        for price in min_price + 1:max_price - 1
            trades, order_id = limitorder!(ob, BUY, price, size, UInt128(1))
            n_orders = n_orders + 1
        end
    end
    t = time()
    limitorder!(ob, SELL, min_price, size*(n_orders - 2), UInt128(1))
    @printf "Matched %.0E orders per second with one order \n" n_orders/(time() - t)
end

println("Matches orders with one large order: 1 iteration")
@time match_large_order(1)

println("Matches orders with one large order: 10 iteration")
@time match_large_order(10)


text:
using Gadfly

function testMultipleStarts(tsize::Int=300, npoints::Int=20)
    SpeedImage = ones(tsize, tsize)
    SourcePoint = rand(2, npoints) .* tsize .+ 1
    return FastMarching.msfm(SpeedImage, SourcePoint, true, true)
end  # function testMultipleStarts

t1 = testMultipleStarts()
draw(PNG("multiplestart.png"), plot(; z=t1, Geom.contour))


text:
@testset "tikhonov" begin
	include("tikhonovRegTest.jl")
end


text:
include("count_table.jl")
include("nri.jl")
include("overlap.jl")
include("sbm.jl")
include("syn_recovery.jl")


text:
using ParametricOptInterface
using BenchmarkTools
using JuMP
using LinearAlgebra
const POI = ParametricOptInterface
import Pkg

function moi_add_variables(N::Int)
    model = direct_model(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            MOI.Utilities.AUTOMATIC,
        ),
    )
    @variable(model, x[i = 1:N])
    return nothing
end

function poi_add_variables(N::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:N])
    return nothing
end

function poi_add_parameters(N::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:N] in POI.Parameter(1))
    return nothing
end

function poi_add_parameters_and_variables(N::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:N] in POI.Parameter(1))
    return nothing
end

function poi_add_parameters_and_variables_alternating(N::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    for i in 1:Int(N / 2)
        @variable(model)
        @variable(model, set = POI.Parameter(1))
    end
    return nothing
end

function moi_add_saf_ctr(N::Int, M::Int)
    model = direct_model(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            MOI.Utilities.AUTOMATIC,
        ),
    )
    @variable(model, x[i = 1:N])
    @constraint(model, cons[i = 1:M], sum(x) >= 1)
    return nothing
end

function poi_add_saf_ctr(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:N])
    @constraint(model, cons[i = 1:M], sum(x) >= 1)
    return nothing
end

function poi_add_saf_variables_and_parameters_ctr(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:Int(N / 2)])
    @variable(model, p[i = 1:Int(N / 2)] in POI.Parameter.(0))
    @constraint(model, con[i = 1:M], sum(x) + sum(p) >= 1)
    return nothing
end

function poi_add_saf_variables_and_parameters_ctr_parameter_update(
    N::Int,
    M::Int,
)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:Int(N / 2)])
    @variable(model, p[i = 1:Int(N / 2)] in POI.Parameter.(0))
    @constraint(model, con[i = 1:M], sum(x) + sum(p) >= 1)
    MOI.set.(model, POI.ParameterValue(), p, 0.5)
    POI.update_parameters!(backend(model))
    return nothing
end

function moi_add_sqf_variables_ctr(N::Int, M::Int)
    model = direct_model(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            MOI.Utilities.AUTOMATIC,
        ),
    )
    @variable(model, x[i = 1:N])
    @constraint(model, con[i = 1:M], dot(x , x) >= 1)
    return nothing
end

function poi_add_sqf_variables_ctr(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:N])
    @constraint(model, con[i = 1:M], dot(x, x) >= 1)
    return nothing
end

function poi_add_sqf_variables_parameters_ctr(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:Int(N / 2)])
    @variable(model, p[i = 1:Int(N / 2)] in POI.Parameter.(1))
    @constraint(model, con[i = 1:M], dot(x , p) >= 1)
    return nothing
end

function poi_add_sqf_variables_parameters_ctr_parameter_update(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:Int(N / 2)])
    @variable(model, p[i = 1:Int(N / 2)] in POI.Parameter.(1))
    @constraint(model, con[i = 1:M], dot(x , p) >= 1)
    MOI.set.(model, POI.ParameterValue(), p, 0.5)
    POI.update_parameters!(backend(model))
    return nothing
end

function poi_add_sqf_parameters_parameters_ctr(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:Int(N / 2)])
    @variable(model, p[i = 1:Int(N / 2)] in POI.Parameter.(1))
    @constraint(model, con[i = 1:M], dot(p , p) >= 1)
    return nothing
end

function poi_add_sqf_parameters_parameters_ctr_parameter_update(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:Int(N / 2)])
    @variable(model, p[i = 1:Int(N / 2)] in POI.Parameter.(1))
    @constraint(model, con[i = 1:M], dot(p , p) >= 1)
    MOI.set.(model, POI.ParameterValue(), p, 0.5)
    POI.update_parameters!(backend(model))
    return nothing
end

function moi_add_saf_obj(N::Int, M::Int)
    model = direct_model(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            MOI.Utilities.AUTOMATIC,
        ),
    )
    @variable(model, x[i = 1:N])
    @objective(model, Min, sum(x))
    return nothing
end

function poi_add_saf_obj(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:N])
    for _ in 1:M
        @objective(model, Min, sum(x))
    end
    return nothing
end

function poi_add_saf_variables_and_parameters_obj(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:Int(N / 2)])
    @variable(model, p[i = 1:Int(N / 2)] in POI.Parameter.(1))
    for _ in 1:M
        @objective(model, Min, sum(x) + sum(p))
    end
    return nothing
end

function poi_add_saf_variables_and_parameters_obj_parameter_update(
    N::Int,
    M::Int,
)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:Int(N / 2)])
    @variable(model, p[i = 1:Int(N / 2)] in POI.Parameter.(1))
    for _ in 1:M
        @objective(model, Min, sum(x) + sum(p))
    end
    for _ in 1:M
        MOI.set.(model, POI.ParameterValue(), p, 0.5)
        POI.update_parameters!(backend(model))
    end
    return nothing
end

function moi_add_sqf_variables_obj(N::Int, M::Int)
    model = direct_model(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            MOI.Utilities.AUTOMATIC,
        ),
    )
    @variable(model, x[i = 1:N])
    for _ in 1:M
        @objective(model, Min, dot(x , x))
    end
    return nothing
end

function poi_add_sqf_variables_obj(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:N])
    for _ in 1:M
        @objective(model, Min, dot(x , x))
    end
    return nothing
end

function poi_add_sqf_variables_parameters_obj(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:Int(N / 2)])
    @variable(model, p[i = 1:Int(N / 2)] in POI.Parameter.(1))
    for _ in 1:M
        @objective(model, Min, dot(x , p))
    end
    return nothing
end

function poi_add_sqf_variables_parameters_obj_parameter_update(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:Int(N / 2)])
    @variable(model, p[i = 1:Int(N / 2)] in POI.Parameter.(1))
    for _ in 1:M
        @objective(model, Min, dot(x , p))
    end
    for _ in 1:M
        MOI.set.(model, POI.ParameterValue(), p, 0.5)
        POI.update_parameters!(backend(model))
    end
    return nothing
end

function poi_add_sqf_parameters_parameters_obj(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:Int(N / 2)])
    @variable(model, p[i = 1:Int(N / 2)] in POI.Parameter.(1))
    for _ in 1:M
        @objective(model, Min, dot(p , p))
    end
    return nothing
end

function poi_add_sqf_parameters_parameters_obj_parameter_update(N::Int, M::Int)
    model = direct_model(
        POI.Optimizer(
            MOI.Utilities.CachingOptimizer(
                MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
                MOI.Utilities.AUTOMATIC,
            ),
        ),
    )
    @variable(model, x[i = 1:Int(N / 2)])
    @variable(model, p[i = 1:Int(N / 2)] in POI.Parameter.(1))
    for _ in 1:M
        @objective(model, Min, dot(p , p))
    end
    for _ in 1:M
        MOI.set.(model, POI.ParameterValue(), p, 0.5)
        POI.update_parameters!(backend(model))
    end
    return nothing
end

function run_benchmarks(N::Int, M::Int)
    println("Pkg status:")
    Pkg.status()
    println("")
    GC.gc()
    println("variables on a MOIU.Model.")
    @btime moi_add_variables($N)
    GC.gc()
    println("variables on a POI.Optimizer.")
    @btime poi_add_variables($N)
    GC.gc()
    println("parameters on a POI.Optimizer.")
    @btime poi_add_parameters($N)
    GC.gc()
    println("parameters and variables on a POI.Optimizer.")
    @btime poi_add_parameters_and_variables($N)
    GC.gc()
    println("alternating parameters and variables on a POI.Optimizer.")
    @btime poi_add_parameters_and_variables_alternating($N)
    GC.gc()
    println("SAF constraint with variables on a MOIU.Model.")
    @btime moi_add_saf_ctr($N, $M)
    GC.gc()
    println("SAF constraint with variables on a POI.Optimizer.")
    @btime poi_add_saf_ctr($N, $M)
    GC.gc()
    println("SAF constraint with variables and parameters on a POI.Optimizer.")
    @btime poi_add_saf_variables_and_parameters_ctr($N, $M)
    GC.gc()
    println("SQF constraint with variables on a MOIU.Model{Float64}.")
    @btime moi_add_sqf_variables_ctr($N, $M)
    GC.gc()
    println("SQF constraint with variables on a POI.Optimizer.")
    @btime poi_add_sqf_variables_ctr($N, $M)
    GC.gc()
    println("SQF constraint with product of variables and parameters on a POI.Optimizer.")
    @btime poi_add_sqf_variables_parameters_ctr($N, $M)
    GC.gc()
    println("SQF constraint with product of parameters on a POI.Optimizer.")
    @btime poi_add_sqf_parameters_parameters_ctr($N, $M)
    GC.gc()
    println("SAF objective with variables on a MOIU.Model.")
    @btime moi_add_saf_obj($N, $M)
    GC.gc()
    println("SAF objective with variables on a POI.Optimizer.")
    @btime poi_add_saf_obj($N, $M)
    GC.gc()
    println("SAF objective with variables and parameters on a POI.Optimizer.")
    @btime poi_add_saf_variables_and_parameters_obj($N, $M)
    GC.gc()
    println("SQF objective with variables on a MOIU.Model.")
    @btime moi_add_sqf_variables_obj($N, $M)
    GC.gc()
    println("SQF objective with variables on a POI.Optimizer.")
    @btime poi_add_sqf_variables_obj($N, $M)
    GC.gc()
    println("SQF objective with product of variables and parameters on a POI.Optimizer.")
    @btime poi_add_sqf_variables_parameters_obj($N, $M)
    GC.gc()
    println("SQF objective with product of parameters on a POI.Optimizer.")
    @btime poi_add_sqf_parameters_parameters_obj($N, $M)
    GC.gc()
    println("Update parameters in SAF constraint with variables and parameters on a POI.Optimizer.")
    @btime poi_add_saf_variables_and_parameters_ctr_parameter_update($N, $M)
    GC.gc()
    println("Update parameters in SAF objective with variables and parameters on a POI.Optimizer.")
    @btime poi_add_saf_variables_and_parameters_obj_parameter_update($N, $M)
    GC.gc()
    println("Update parameters in SQF constraint with product of variables and parameters on a POI.Optimizer.")
    @btime poi_add_sqf_variables_parameters_ctr_parameter_update($N, $M)
    GC.gc()
    println("Update parameters in SQF constraint with product of parameters on a POI.Optimizer.")
    @btime poi_add_sqf_parameters_parameters_ctr_parameter_update($N, $M)
    GC.gc()
    println("Update parameters in SQF objective with product of variables and parameters on a POI.Optimizer.")
    @btime poi_add_sqf_variables_parameters_obj_parameter_update($N, $M)
    GC.gc()
    println("Update parameters in SQF objective with product of parameters on a POI.Optimizer.")
    @btime poi_add_sqf_parameters_parameters_obj_parameter_update($N, $M)
    return nothing
end

N = 10_000
M = 100
run_benchmarks(N, M)


text:
include("Include.jl")

# extra:
using Flux
using Flux: @epochs
using BSON: @load

# load training set -
full_training_data_frame = load_training_data()

# filter -
has_TM_flag = 0
experimental_data_table = filter([:visitid, :TM] => (x, y) -> ((x == 2 || x == 3) && y == has_TM_flag), full_training_data_frame)

# what is the size of the system?
(P,D) = size(experimental_data_table)

# initialize some space -
outdim = 5
validation_data_array = Array{Float32,2}(undef, P, outdim)
simulated_data_array = Array{Array{Float32,1},1}()

for i ∈ 1:P

    # we used all data *except* i to train the model -
    validation_input_data = convert.(Float32, Vector(experimental_data_table[i, 3:14]))
    validation_output_data = convert.(Float32, Vector(experimental_data_table[i, 15:19]))

    # grab this output data -
    for k ∈ 1:outdim
        validation_data_array[i,k] = validation_output_data[k]
    end

    for j = 1:P

        # load a model -
        model_name = "deep_coag_model-L$(j)O-TF-TM-$(has_TM_flag)-V2-V3.bson"
        model_file_path = joinpath(_PATH_TO_MODELS, model_name)
        @load model_file_path deep_coag_model

        # sample -
        simulated_output_data = sample(deep_coag_model, validation_input_data)

        # grab -
        push!(simulated_data_array,abs.(simulated_output_data))
    end
end

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
module Service

using ..Model, ..Mapper

function createPaper(obj)
    # obj = json dict
    @assert haskey(obj, :title) && !isempty(obj.title)
    @assert haskey(obj, :abst) && !isempty(obj.abst)
    @assert haskey(obj, :authors) && !isempty(obj.authors)
    @assert haskey(obj, :publishdate) && !isempty(obj.publishdate)

    paper = Paper(obj.title, obj.abst, obj.authors, obj.publishdate)
    Mapper.store!(paper)
    return paper
end

getPaper(id) = Mapper.get(id)

function updatePaper(id, updated)
    paper=Mapper.get(id)
    paper.title = updated.title
    paper.abst = updated.abst
    paper.authors = updated.authors
    paper.publishdate = updated.publishdate
    Mapper.store!(paper)
    return paper
end

function deletePaper(id)
    Mapper.delete(id)
    return
end

function pickPaperToRead()
    papers = Mapper.getAllPapers()
    leastTimesPicked = minimum(x -> x.timespicked, papers)
    leastPickedPapers = filter(x -> x.timespicked == leastTimesPicked, papers)
    pickedPaper = rand(leastPickedPapers)
    pickedPaper.timespicked += 1
    Mapper.store!(pickedPaper)
    return pickedPaper
end

end


text:
#Mapping functions
##################
# function index2d(k,nx)
#     j = floor(Int,(k-1)/nx)+1
#     # TODO: add consistency check with ny?
#     return k-(j-1)*nx,j
# end

function index1d(i,j,nx)
    # TODO: add consistency check with ny?
    return i+(j-1)*nx
end

"""
    compute_equilibrium!(...)

* rel_error is the tolerance for global temperature equilibrium (default is 2e-5).
* max_years is the maximum number of annual cycles to be computed when searching for equilibrium
* update_heat_capacity : manage the monthly update function for the heat capacity
* update_solar_forcing : manage the monthly update function for the solar forcing
"""
function compute_equilibrium!(discretization; max_years=100, rel_error=2e-5, verbose=true, update_heat_capacity = false,update_solar_forcing = false)
    @unpack mesh, model, num_steps_year, annual_temperature, rhs, last_rhs  = discretization
    @unpack nx, dof = mesh

    average_temperature = average_temperature_old = area_weighted_average(view(annual_temperature, :, num_steps_year), mesh)

    if verbose
      println("year","  ","Average Temperature")
      println(0,"  ",average_temperature_old)
    end

    for year in 1:max_years
        average_temperature = 0.0
        for time_step in 1:num_steps_year

            update_monthly_params!(model, discretization.lu_decomposition, mesh, num_steps_year, time_step, update_heat_capacity, update_solar_forcing)

            old_time_step = (time_step == 1) ? num_steps_year : time_step - 1
            update_rhs!(rhs, mesh, num_steps_year, time_step, view(annual_temperature, :, old_time_step), model, last_rhs)

            # Use in-place operation `ldiv!` instead of `\` to avoid allocations
            ldiv!(view(annual_temperature, :, time_step), discretization.lu_decomposition, rhs)

            annual_temperature[1:nx, time_step] .= annual_temperature[1, time_step]
            annual_temperature[dof-nx+1:dof, time_step] .= annual_temperature[dof, time_step]

            average_temperature += area_weighted_average(view(annual_temperature, :, time_step), mesh)
        end
        average_temperature = average_temperature/num_steps_year
        if verbose
          println(year,"  ",average_temperature)
        end

        if (abs(average_temperature-average_temperature_old)<rel_error)
            if verbose
              println("EQUILIBRIUM REACHED!")
            end
            break
        end

        average_temperature_old = average_temperature
    end

    return average_temperature
end

"""
    compute_evolution!(...)

Compute the evolution of the mean temperature with varying CO2 levels.
"""
function compute_evolution!(discretization, co2_concentration_at_step, year_start, year_end; verbose=true, update_heat_capacity = false,update_solar_forcing = false)
    @unpack mesh, model, num_steps_year, annual_temperature, rhs, last_rhs  = discretization
    @unpack nx, dof = mesh

    if verbose
      println("year","  ","Average Temperature")
    end

    max_years = year_end - year_start + 1

    # Allocate output arrays
    year_at_step = zeros(Float64,max_years*num_steps_year+1)
    mean_temperature_at_step = zeros(Float64,max_years*num_steps_year+1)
    year_array = zeros(Float64,max_years)
    mean_temperature_yearly = zeros(Float64,max_years)

    # Fill year arrays
    year_at_step[1] = year_start + 0.21644  # The simulation starts at the vernal equinox
    for year in 2:size(year_at_step,1)
        year_at_step[year] = year_at_step[year-1] + 1.0/num_steps_year
    end
    year_array[1] = year_start + 0.71644    # The simulation starts at the vernal equinox
    for year in 2:max_years
        year_array[year] = year_array[year-1] + 1.0
    end

    # Initialize the temperature array
    mean_temperature_at_step[1] = area_weighted_average(view(annual_temperature, :, num_steps_year), mesh)
    step = 1

    for year in 1:max_years
        average_temperature = 0.0
        for time_step in 1:num_steps_year
            set_co2_concentration!(model, co2_concentration_at_step[step])
            update_monthly_params!(model, discretization.lu_decomposition, mesh, num_steps_year, time_step, update_heat_capacity, update_solar_forcing)
            old_time_step = (time_step == 1) ? num_steps_year : time_step - 1
            update_rhs!(rhs, mesh, num_steps_year, time_step, view(annual_temperature, :, old_time_step), model, last_rhs)

            # Use in-place operation `ldiv!` instead of `\` to avoid allocations
            ldiv!(view(annual_temperature, :, time_step), discretization.lu_decomposition, rhs)

            annual_temperature[1:nx, time_step] .= annual_temperature[1, time_step]
            annual_temperature[dof-nx+1:dof, time_step] .= annual_temperature[dof, time_step]

            # Compute mean temperature
            step += 1
            mean_temperature_at_step[step] = area_weighted_average(view(annual_temperature, :, time_step), mesh)
            average_temperature += mean_temperature_at_step[step]
        end
        average_temperature = average_temperature/num_steps_year
        mean_temperature_yearly[year] = average_temperature
        if verbose
          println(year_start+year-1,"  ",average_temperature)
        end
    end

    return (; year_at_step, mean_temperature_at_step, year_array, mean_temperature_yearly)
end

"""
Computes mean temperature in the globe at a specific time
"""
function area_weighted_average(vector,mesh)
    @unpack nx,ny,area,dof = mesh

    average_vector = 0.0

    # Contribution of inner points
    for j=2:ny-1
        for i=1:nx
            row_idx = index1d(i,j,nx)
            average_vector += area[j] * vector[row_idx]
        end
    end

    # Poles
    average_vector += area[1] * (vector[1] + vector[dof])

    return average_vector
end

"""
compute_matrix(mesh,num_steps_year,model)
Compute the system matrix for the linear solver
"""
function compute_matrix(mesh,num_steps_year,model)
    @unpack nx,ny,dof,h,geom,csc2,cot,area = mesh
    @unpack diffusion_coeff,heat_capacity,radiative_cooling_feedback = model
    matrix = spzeros(Float64,dof,dof)
    sh2 = 1 / h^2
    # Inner DOFs (c coefficients are divided by h^2)
    for j=2:ny-1
        for i=1:nx

            # Compute coefficients
            c0 = 2 * num_steps_year + (2 * sh2 * diffusion_coeff[i,j] * (1 + csc2[j]) + radiative_cooling_feedback) / heat_capacity[i,j,2]

            # Get diffusion coefficient
            if (i == 1) # Periodic BC
                d_phi = (diffusion_coeff[2,j] - diffusion_coeff[nx,j]) / (2 * h)
            elseif (i == nx) # Periodic BC
                d_phi = (diffusion_coeff[1,j] - diffusion_coeff[nx-1,j]) / (2 * h)
            else # Inner DOFs
                d_phi = (diffusion_coeff[i+1,j] - diffusion_coeff[i-1,j]) / (2 * h)
            end

            c1 = sh2 * csc2[j] * (diffusion_coeff[i,j] - 0.5 * h * d_phi)
            c3 = sh2 * csc2[j] * (diffusion_coeff[i,j] + 0.5 * h * d_phi)

            d_theta= (diffusion_coeff[i,j+1] - diffusion_coeff[i,j-1]) / (2 * h)
            c2 = sh2 * (diffusion_coeff[i,j] - 0.5 * h * (diffusion_coeff[i,j] * cot[j] + d_theta))
            c4 = sh2 * (diffusion_coeff[i,j] + 0.5 * h * (diffusion_coeff[i,j] * cot[j] + d_theta))

            # Fill matrix A
            row_idx = index1d(i,j,nx)
            col_idx_c0 = row_idx

            if (i==1) # Periodic BC
                col_idx_c1 = index1d(nx,j,nx)
                col_idx_c3 = row_idx + 1
            elseif (i==nx) # Periodic BC
                col_idx_c1 = row_idx - 1
                col_idx_c3 = index1d(1,j,nx)
            else
                col_idx_c1 = row_idx - 1
                col_idx_c3 = row_idx + 1
            end
            col_idx_c2 = row_idx - nx
            col_idx_c4 = row_idx + nx

            matrix[row_idx,col_idx_c0] = c0
            matrix[row_idx,col_idx_c1] = -c1 / heat_capacity[i,j,2]
            matrix[row_idx,col_idx_c2] = -c2 / heat_capacity[i,j,2]
            matrix[row_idx,col_idx_c3] = -c3 / heat_capacity[i,j,2]
            matrix[row_idx,col_idx_c4] = -c4 / heat_capacity[i,j,2]

            # Get rhs
            # rhs[row_idx] = -(c0 * Temp[col_idx_c0] - c1 * Temp[col_idx_c0] - c2 * Temp[col_idx_c2] - c3 * Temp[col_idx_c3] - c4 * Temp[col_idx_c4]) + 4 * heat_capacity[i,j]
        end
    end

    # Poles
    total_area = area[1] + area[2]

    diff_total_np = zeros(Float64,nx)
    diff_total_sp = zeros(Float64,nx)
    for i=1:nx
        diff_total_np[i] = (area[1] * diffusion_coeff[1,1]  + area[2] * diffusion_coeff[i,2]) / total_area
        diff_total_sp[i] = (area[1] * diffusion_coeff[1,ny] + area[2] * diffusion_coeff[i,ny-1]) / total_area
    end
    diff_total_sum_np = sum(diff_total_np)
    diff_total_sum_sp = sum(diff_total_sp)

    gc_np = 2 * num_steps_year + (geom * diff_total_sum_np + radiative_cooling_feedback) / heat_capacity[1,1,2]    # north pole
    gc_sp = 2 * num_steps_year + (geom * diff_total_sum_sp + radiative_cooling_feedback) / heat_capacity[1,ny,2]    # south pole

    for i=1:nx
        #North pole
        row_idx = i
        matrix[row_idx,row_idx] = gc_np
        for ii=1:nx
            matrix[row_idx,nx+ii] = - geom * diff_total_np[ii] / heat_capacity[1,1,2]
        end

        #South pole
        row_idx = dof + 1 - i
        matrix[row_idx,row_idx] = gc_sp
        for ii=1:nx
            matrix[row_idx,dof-2*nx+ii] = - geom * diff_total_sp[ii] / heat_capacity[1,ny,2]
        end
    end
    return matrix
end

function update_rhs!(rhs, mesh, num_steps_year, time_step, temperature, model,  last_rhs)
    @unpack nx,ny = mesh
    @unpack heat_capacity, solar_forcing, radiative_cooling_co2 = model
    for j=1:ny
        for i=1:nx
            row_idx = index1d(i,j,nx)

            rhs[row_idx] = 4 * temperature[row_idx] * num_steps_year  - last_rhs[row_idx] + solar_forcing[i,j,time_step] / heat_capacity[i,j,2] - (radiative_cooling_co2[2] / heat_capacity[i,j,2] + radiative_cooling_co2[1] / heat_capacity[i,j,1])
            if (time_step == 1)
                rhs[row_idx] += solar_forcing[i,j,num_steps_year] / heat_capacity[i,j,1]
            else
                rhs[row_idx] += solar_forcing[i,j,time_step-1] / heat_capacity[i,j,1]
            end
        end
    end
    last_rhs .= rhs
end

"""
 update_monthly_params!(...)

* The routine updates monthly all important parameters of the model,
  which are necessary for the passed update settings.
  When updating the heat capacity, the LU decomposition must take place again.
  Update start is the first week of April.

"""
function update_monthly_params!(model, lu_decomposition, mesh, num_steps_year, time_step, update_heat_capacity, update_solar_forcing)

    if mod(time_step,4) == 2 && time_step in (2:48)
        if update_heat_capacity == true && update_solar_forcing == true
         update_model(model, mesh,time_step, true,true,false,true,true)
         lu_decomposition = compute_lu_matrices(mesh, model, num_steps_year)
        elseif update_heat_capacity == true && update_solar_forcing == false
            update_model(model, mesh, time_step, true,false,false,true,false)
            lu_decomposition = compute_lu_matrices(mesh, model, num_steps_year)
        elseif update_heat_capacity == false && update_solar_forcing == true
            update_model(model, mesh, time_step, true,true,false,false,true)
        elseif update_heat_capacity == false && update_solar_forcing == false
            update_model(model, mesh, time_step, false,false,false,false,false)
        end
    else
        if update_heat_capacity
            model.heat_capacity[:,:,1] = model.heat_capacity[:,:,2]
        end
    end
end


text:
using RandomMatrix, LinearAlgebra
using Test

@test randDiagonal(2) !== nothing
@test randTriangular(1:3,3,upper=false,Diag=false) !== nothing
@test !isreal(randMatrix(10)|>eigvals)

@test isreal(randHermitian(3)|>eigvals)
@test randHermitian(ComplexNormal(im,2),3,diag=Elliptic(0.1,c=im,R=9)) !==nothing

@test pdf(MarchenkoPastur(rand()),0) == 0


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ ea747ce8-5f4f-11ec-3fff-191d0010c672
using DelimitedFiles

# ╔═╡ c2b811e0-ed4f-4b5c-9aad-79b02298fa2c
using Statistics

# ╔═╡ 7ec3a4cd-95d5-42c3-878a-bb29228fae89
d = readdlm("inputs/day_03.txt", Int64)

# ╔═╡ 5047baba-22aa-4fe3-8cec-ab8b1edbfb14

# ╔═╡ 19e8e799-59f0-47e0-8e72-69fa572bd082
# each column is a binary number
raw_numbers = hcat(reverse.(digits.(d, pad=12))...)

# ╔═╡ 30a18fb3-55f7-4404-85a4-d8e866223bd5
most_common(A) = mean(A) > 0.5 ? 1 : 0

# ╔═╡ 9f7c1bf1-2984-4bb5-8c2a-e68a830d0756
least_common(A) = mean(A) > 0.5 ? 0 : 1

# ╔═╡ 0e10b830-d1ed-4fb5-98d8-04f7a0afe1af
γ = parse(Int, string(mapslices(most_common, raw_numbers; dims=2)'...), base=2)

# ╔═╡ e1af0c80-14d7-4b4b-929d-b44b119f3bf9
ϵ = parse(Int, string(mapslices(least_common, raw_numbers; dims=2)'...), base=2)

# ╔═╡ a15c4c54-1aa1-40bf-a4a8-458b5f26bab1
γ * ϵ

# ╔═╡ b2bc6205-66e1-43f6-b642-d0eac97f8d6e
md"""
# Part 2
"""

# ╔═╡ 536c91f5-2962-4d81-afa1-b5dd659c1a7f

# ╔═╡ b3644b64-fa51-4516-a3c6-4b51a26332fd

# ╔═╡ c4bf6648-cc81-4b80-b790-fe6542505151

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0"
manifest_format = "2.0"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
"""

# ╔═╡ Cell order:
# ╠═ea747ce8-5f4f-11ec-3fff-191d0010c672
# ╠═c2b811e0-ed4f-4b5c-9aad-79b02298fa2c
# ╠═7ec3a4cd-95d5-42c3-878a-bb29228fae89
# ╠═5047baba-22aa-4fe3-8cec-ab8b1edbfb14
# ╠═19e8e799-59f0-47e0-8e72-69fa572bd082
# ╠═30a18fb3-55f7-4404-85a4-d8e866223bd5
# ╠═9f7c1bf1-2984-4bb5-8c2a-e68a830d0756
# ╠═0e10b830-d1ed-4fb5-98d8-04f7a0afe1af
# ╠═e1af0c80-14d7-4b4b-929d-b44b119f3bf9
# ╠═a15c4c54-1aa1-40bf-a4a8-458b5f26bab1
# ╠═b2bc6205-66e1-43f6-b642-d0eac97f8d6e
# ╠═536c91f5-2962-4d81-afa1-b5dd659c1a7f
# ╠═b3644b64-fa51-4516-a3c6-4b51a26332fd
# ╠═c4bf6648-cc81-4b80-b790-fe6542505151
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using Adapt

using Oceananigans.Architectures: AbstractArchitecture
using Oceananigans.Grids: AbstractGrid
using Oceananigans.BoundaryConditions
using Oceananigans.BoundaryConditions: AbstractBoundaryConditionClassification

using KernelAbstractions: NoneEvent

import Base: show

import Oceananigans.BoundaryConditions: bc_str, print_condition

import Oceananigans.BoundaryConditions:
    apply_x_bcs!,
    apply_y_bcs!,
    apply_x_east_bc!,
    apply_x_west_bc!,
    apply_y_south_bc!,
    apply_y_north_bc!

struct CubedSphereExchange <: AbstractBoundaryConditionClassification end

const CubedSphereExchangeBC = BoundaryCondition{<:CubedSphereExchange}

bc_str(::CubedSphereExchangeBC) ="CubedSphereExchange"

CubedSphereExchangeBoundaryCondition(val; kwargs...) = BoundaryCondition(CubedSphereExchange, val; kwargs...)

struct CubedSphereExchangeInformation{F, S}
    from_face :: F
      to_face :: F
    from_side :: S
      to_side :: S
end

CubedSphereExchangeInformation(; from_face, to_face, from_side, to_side) =
    CubedSphereExchangeInformation(from_face, to_face, from_side, to_side)

Base.show(io::IO, ex::CubedSphereExchangeInformation) =
    print(io, "CubedSphereExchangeInformation: (from: face $(ex.from_face) $(ex.from_side) side, to: face $(ex.to_face) $(ex.to_side) side)")

print_condition(info::CubedSphereExchangeInformation) =
    "(from: face $(info.from_face) $(info.from_side) side, to: face $(info.to_face) $(info.to_side) side)"

function inject_cubed_sphere_exchange_boundary_conditions(field_bcs, face_number, face_connectivity)

    west_exchange_info = CubedSphereExchangeInformation(
        from_face = face_number,
        from_side = :west,
          to_face = face_connectivity[face_number].west.face,
          to_side = face_connectivity[face_number].west.side
    )

    east_exchange_info = CubedSphereExchangeInformation(
        from_face = face_number,
        from_side = :east,
          to_face = face_connectivity[face_number].east.face,
          to_side = face_connectivity[face_number].east.side
    )

    south_exchange_info = CubedSphereExchangeInformation(
        from_face = face_number,
        from_side = :south,
          to_face = face_connectivity[face_number].south.face,
          to_side = face_connectivity[face_number].south.side
    )

    north_exchange_info = CubedSphereExchangeInformation(
        from_face = face_number,
        from_side = :north,
          to_face = face_connectivity[face_number].north.face,
          to_side = face_connectivity[face_number].north.side
    )

    west_exchange_bc = CubedSphereExchangeBoundaryCondition(west_exchange_info)
    east_exchange_bc = CubedSphereExchangeBoundaryCondition(east_exchange_info)
    south_exchange_bc = CubedSphereExchangeBoundaryCondition(south_exchange_info)
    north_exchange_bc = CubedSphereExchangeBoundaryCondition(north_exchange_info)

    return FieldBoundaryConditions(west_exchange_bc,
                                   east_exchange_bc,
                                   south_exchange_bc,
                                   north_exchange_bc,
                                   field_bcs.bottom,
                                   field_bcs.top,
                                   field_bcs.immersed)
end

Adapt.adapt_structure(to, ::CubedSphereExchangeInformation) = nothing
Adapt.adapt_structure(to, ::CubedSphereExchangeBC) = nothing

# Don't "apply fluxes" across CubedSphere boundaries
@inline apply_x_east_bc!(  Gc, loc, ::CubedSphereExchangeBC, args...) = nothing
@inline apply_x_west_bc!(  Gc, loc, ::CubedSphereExchangeBC, args...) = nothing
@inline apply_y_north_bc!( Gc, loc, ::CubedSphereExchangeBC, args...) = nothing
@inline apply_y_south_bc!( Gc, loc, ::CubedSphereExchangeBC, args...) = nothing

apply_x_bcs!(Gc, ::AbstractGrid, c, ::CubedSphereExchangeBC, ::CubedSphereExchangeBC, ::AbstractArchitecture, args...) = NoneEvent()
apply_y_bcs!(Gc, ::AbstractGrid, c, ::CubedSphereExchangeBC, ::CubedSphereExchangeBC, ::AbstractArchitecture, args...) = NoneEvent()


text:
############################################################################################
function init(
    ::Type{NUTS},
    config::ConfigNUTS,
    objective::Objective,
    difftune::AbstractDifferentiableTune,
)
    return NUTS(
        ModelWrappers.log_density_and_gradient(objective, difftune), difftune, config.energy
    )
end

function infer(
    _rng::Random.AbstractRNG,
    diagnostics::Type{AbstractDiagnostics},
    kernel::NUTS,
    model::ModelWrapper,
    data::D,
) where {D}
    return DiagnosticsNUTS{model.info.reconstruct.default.output}
end

############################################################################################
# Export
export init, infer


text:
function findtwin(edges,e)
    #edges = decompose(Face{2,Int},faces)
    for j in 1:length(edges)
        if edges[j][2]==e[1] && edges[j][1]==e[2]
            return j
        end
    end
end

"""
    EdgeDS(faces::PlainDS)

Constructs and returns edge based datastructure `EdgeDS` from plain face based datastructure PlainDS. Half edge based datastructure `EdgeDS` stores list of edges consisting of base vertex, next edge index and twin edge index.
"""
struct EdgeDS
    edges::Array{Face{3,Int},1} ### It is equal to PlainDS. That is why we should introduce a new type.

    # row is data assciated with edge. It contains a base vertex, next edge and the twin edge
    #neighs::Array{Face{3,Int},1} ### For now keeping simple
    #vfaces::Array{Int,1}
    function EdgeDS(faces::PlainDS)

        edges = decompose(Face{2,Int},faces)

        rows = Face{3,Int}[]

        for ei in edges
            basev = ei[1]
            twin = -1
            for j in 1:length(edges)
                if edges[j][2]==ei[1] && edges[j][1]==ei[2]
                    twin = j
                    break
                end
            end

            ### Now I could find triangle with specific edge
            ### Orientation is important

            ti = findtriangle(faces,ei)
            face = faces[ti]
            nedge = getedge(face,ei)

            next = -1
            for j in 1:length(edges)
                if edges[j]==nedge
                    next = j
                    break
                end
            end

            # @show basev, next, twin
            push!(rows,Face(basev,next,twin))
        end

        return new(rows)
    end

end

function findtriangle(faces::PlainDS,edge)

    equal(v1,v2) = v1==edge[1] && v2==edge[2]

    for i in 1:length(faces)
        if in(edge[1],faces[i]) & in(edge[2],faces[i])
            v1,v2,v3 = faces[i]
            if equal(v1,v2) || equal(v2,v3) || equal(v3,v1)
                return i
            end
        end
    end

end

# Returns an edge which we need to look up
function getedge(t,edge)
    v1,v2,v3=t
    if edge==Face(v1,v2)
        return Face(v2,v3)
    elseif edge==Face(v2,v3)
        return Face(v3,v1)
    else
        return Face(v1,v2)
    end
end

### A little bit of thought

### the state could consist of

### This one would construct the beginign vertex
### So aftger that one is able to do VertexRing(v,t::EdgeDS)

function VertexRing(v::Int,t::EdgeDS)
    edges = t.edges

    for ei in edges
        if v==ei[1]
            next = ei[2]
            return VertexRing(v,next,t)
        end
    end
end

function Base.iterate(iter::VertexRing{EdgeDS})
    e = iter.start
    v = iter.t.edges[e][1]
    return (v,e)
end

function Base.iterate(iter::VertexRing{EdgeDS},e)
    edges = iter.t.edges

    ea = edges[e][2]
    eb = edges[ea][3] # the twin
    ec = edges[eb][2]

    if ec==iter.start
        return nothing
    else
        return (edges[ec][1],ec)
    end
end

function EdgeRing(v::Int,t::EdgeDS)
    iter = VertexRing(v,t)
    return PairIterator(iter)
end

### One would need to mark used rows one by one and loop over the rows. Implementation is rather easy. Also important for exporting.
Faces(t::EdgeDS) = error("Not implemented")

### Selects only edges if their next index is larger to avoid dublications. For simplicity let's not write it as an iterator.
Edges(t::EdgeDS) = error("Not implemented")


text:
using Parameters, Plots, JuMP, LaTeXStrings

const block_colors = cgrad(:inferno)

# Function for drawing rectangles
rect(x, y, w, h) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])

"""Creates a planogram which visualizes the product placement on the shelves."""
function plot_planogram(products, shelves, blocks, P_b, H_s, H_p, W_p, W_s,
                        SK_p, n_ps, o_s, b_bs, x_bs, z_bs;
                        draw_block_allocations=true)
    # Initialize the plot
    plt = plot(
        legend=:none,
        background=:lightgray,
        size=(780, 400)
    )

    # Cumulative shelf heights
    y_s = vcat([0], cumsum([H_s[s] for s in shelves]))

    # Draw products
    for b in blocks, s in shelves
        x = x_bs[b, s]
        for p in P_b[b]
            stack = max(min(div(H_s[s], H_p[p]), SK_p[p]), 1)
            for i in 1:n_ps[p, s]
                y = y_s[s-shelves[1]+1]
                for j in 1:stack
                    plot!(plt, rect(x, y, W_p[p], H_p[p]),
                          color=block_colors[b/length(blocks)])
                    y += H_p[p]
                end
                x += W_p[p]
            end
        end
    end

    # Draw shelves
    for s in shelves
        y = y_s[s-shelves[1]+1]
        plot!(plt, [0, W_s[s]], [y, y], color=:black)
    end
    plot!(plt, [0, W_s[shelves[end]]], [y_s[end], y_s[end]],
          color=:black, linestyle=:dash)

    # Draw block allocations
    if draw_block_allocations
        for b in blocks, s in shelves
            if iszero(z_bs[b, s])
                continue
            end
            color = block_colors[b/length(blocks)]
            y = y_s[s-shelves[1]+1]
            scatter!(plt, [x_bs[b, s]], [y],
                     color=color, markerstrokewidth=0, markersize = 2.5)
            plot!(plt, [x_bs[b, s], x_bs[b, s] + b_bs[b, s]], [y, y],
                  color=color, linewidth=3)
        end
    end

    return plt
end

"""Creates a planogram which visualizes the product placement on the shelves without blocks."""
function plot_planogram_no_blocks(products, shelves, blocks, P_b, H_s, H_p,
                                  W_p, W_s, SK_p, n_ps, o_s)
    # Initialize the plot
    plt = plot(
        legend=:none,
        background=:lightgray,
        size=(780, 400)
    )

    # Cumulative shelf heights
    y_s = vcat([0], cumsum([H_s[s] for s in shelves]))

    # Draw products
    for s in shelves
        x = 0
        for b in blocks
            for p in P_b[b]
                stack = max(min(div(H_s[s], H_p[p]), SK_p[p]), 1)
                for i in 1:n_ps[p, s]
                    y = y_s[s-shelves[1]+1]
                    for j in 1:stack
                        plot!(plt, rect(x, y, W_p[p], H_p[p]),
                              color=block_colors[b/length(blocks)])
                        y += H_p[p]
                    end
                    x += W_p[p]
                end
            end
        end
    end

    # Draw shelves
    for s in shelves
        y = y_s[s-shelves[1]+1]
        plot!(plt, [0, W_s[s]], [y, y], color=:black)
    end
    plot!(plt, [0, W_s[shelves[end]]], [y_s[end], y_s[end]],
          color=:black, linestyle=:dash)

    return plt
end

"""Create a planogram for each module."""
function plot_planograms(parameters::Params, variables::Variables)
    @unpack products, shelves, blocks, S_m, P_b, H_s, H_p, W_p, W_s, SK_p =
            parameters
    @unpack n_ps, o_s, b_bs, x_bs, z_bs = variables
    return [plot_planogram(products, shelves′, blocks, P_b, H_s, H_p, W_p, W_s,
                           SK_p, n_ps, o_s, b_bs, x_bs, z_bs)
                           for shelves′ in S_m]
end

"""Create a planogram for each module without blocks."""
function plot_planograms_no_blocks(parameters::Params, variables::Variables)
    @unpack products, shelves, blocks, S_m, P_b, H_s, H_p, W_p, W_s, SK_p =
            parameters
    @unpack n_ps, o_s = variables
    return [plot_planogram_no_blocks(
                products, shelves′, blocks, P_b, H_s, H_p, W_p, W_s, SK_p,
                n_ps, o_s) for shelves′ in S_m]
end

"""Creates a barchart of number of product facings per product."""
function plot_product_facings(parameters::Params, variables::Variables)
    @unpack products, shelves, blocks, P_b, N_p_max = parameters
    @unpack n_ps = variables

    colors = [block_colors[b/length(blocks)] for b in blocks for p in P_b[b]]

    # Plot maximum number of facings.
    plt = bar(
        N_p_max,
        xlabel=L"$p$",
        ylabel=L"$n_{p,s}$",
        linewidth=0,
        color=colors,
        background=:lightgray,
        legend=:none,
        alpha=0.3,
        tickfontsize=6,
        xticks=vcat([1], [last(P_b[b]) for b in blocks]),
        size=(780, 400)
    )

    # Plot number of facings placed on to the shelves.
    bar!(
        plt,
        [sum(n_ps[p, s] for s in shelves) for p in products],
        color=colors,
        linewidth=0,
        legend=:none,
        background=:lightgray
    )

    return plt
end

"""Bar chart of demand and sales per product."""
function plot_demand_and_sales(parameters::Params, variables::Variables)
    @unpack blocks, P_b, D_p = parameters
    @unpack s_p = variables

    # FIXME?
    colors = [block_colors[b/length(blocks)] for b in blocks for p in P_b[b]]

    # Plot demand of products
    plt = bar(
        D_p,
        alpha=0.3,
        linewidth=0,
        background=:lightgray,
        xlabel=L"$p$",
        ylabel=L"$D_p$ and $s_p$",
        color=colors,
        tickfontsize=6,
        xticks=vcat([1], [last(P_b[b]) for b in blocks]),
        legend=:none,
        size=(780, 400)
    )

    # Plot products sold
    bar!(
        plt,
        s_p,
        alpha=1.0,
        linewidth=0,
        color=colors
    )
    return plt
end

"""Percentage of demand satisfied by sales per block."""
function plot_demand_sales_percentage(parameters::Params, variables::Variables)
    @unpack blocks, P_b, D_p = parameters
    @unpack s_p = variables

    pr = [sum(s_p[p] for p in P_b[b])/sum(D_p[p] for p in P_b[b])
          for b in blocks]
    bar(
        pr,
        color=[block_colors[b/length(blocks)] for b in blocks],
        ylims=(0, 1),
        xlabel=L"$b$",
        xticks=1:1:size(blocks, 1),
        legend=:none,
        background=:lightgray,
    )
end

"""Plot the total amount of allocated facings per product per block."""
function plot_allocation_amount(parameters::Params, variables::Variables)
    @unpack shelves, blocks, P_b = parameters
    @unpack n_ps = variables

    pr = [sum(n_ps[p, s] for s in shelves for p in P_b[b]) for b in blocks]
    plt = bar(
        pr,
        color=[block_colors[b/length(blocks)] for b in blocks],
        xlabel=L"$b$",
        # ylabel=L"",
        xticks=1:1:size(blocks, 1),
        legend=:none,
        background=:lightgray
    )
    return plt
end

"""Plot the percentage of allocated facings of maximum facings per block."""
function plot_allocation_percentage(parameters::Params, variables::Variables, optimizer)
    @unpack products, shelves, blocks, modules, P_b, S_m, N_p_min, N_p_max,
            G_p, R_p, D_p, L_p, W_p, H_p, M_p, SK_p, M_s_min, M_s_max, W_s,
            H_s, L_s, P_ps, SL, w1, w2, w3 = parameters
    @unpack n_ps = variables

    pr = [sum(n_ps[p, s] for p in P_b[b] for s in shelves) for b in blocks]
    pr_max = []
    for b in blocks
        specs = Specs(height_placement=false, blocking=false)
        parameters2 = Params(
            P_b[b], shelves, Integer[], Integer[], [Integer[]], [Integer[]],
            N_p_min, N_p_max, G_p, R_p, D_p, L_p, W_p, H_p, M_p, SK_p, M_s_min,
            M_s_max, W_s, H_s, L_s, P_ps, SL, w1, w2, w3)
        model = ShelfSpaceAllocationModel(parameters2, specs)
        optimize!(model, optimizer)
        n_ps_max = value.(model[:n_ps])
        push!(pr_max, sum(n_ps_max[p, s] for p in P_b[b] for s in shelves))
    end

    plt = bar(
        pr ./ pr_max,
        color=[block_colors[b/length(blocks)] for b in blocks],
        ylims=(0, 1),
        xlabel=L"$b$",
        # ylabel=L"",
        xticks=1:1:size(blocks, 1),
        legend=:none,
        background=:lightgray
    )
    return plt
end


text:
using Statistics

struct GaussNum{T<:Real} <: Real
    μ::T
    σ::T
end
Statistics.mean(x::GaussNum) = x.μ
Statistics.var(x::GaussNum) = x.σ^2
Statistics.std(x::GaussNum) = x.σ
GaussNum(x,y) = GaussNum(promote(x,y)...)
±(x,y) = GaussNum(x,y)
Base.convert(::Type{T}, x::T) where T<:GaussNum = x
Base.convert(::Type{GaussNum{T}}, x::Number) where T = GaussNum(x,zero(T))
Base.promote_rule(::Type{GaussNum{T}}, ::Type{S}) where {T,S} = GaussNum{T}
Base.promote_rule(::Type{GaussNum{T}}, ::Type{GaussNum{T}}) where T = GaussNum{T}

isuncertain(x::GaussNum) = x.σ!=0
isuncertain(x::Number) = false

# TODO: add AD/arithmetic funcs

using Plots
@recipe function plot(ts::AbstractVector, xs::AbstractVector{<:GaussNum})
    # you can set a default value for an attribute with `-->`
    # and force an argument with `:=`
    μs = [x.μ for x in xs]
    σs = [x.σ for x in xs]
    @series begin
        :seriestype := :path
        # ignore series in legend and color cycling
        primary := false
        linecolor := nothing
        fillcolor := :gray
        fillalpha := 0.5
        fillrange := μs .- σs
        # ensure no markers are shown for the error band
        markershape := :none
        # return series data
        ts, μs .+ σs
    end
    ts, μs
end


text:
"""
get_filtermatrix(tin,tout)
# Arguments
- `tin`: input times
- `tout`: output times
# Output
- `Eout2in`: matrix that maps from tout to tin
- `Fin2out`: matrix that maps from tin to tout
"""
function get_filtermatrix(tin,tout)
    # function get_matrixfilter(tin,tout)

    ntin  = length(tin)
    ntout = length(tout)

    # assume evenly spaced in time
    Δtout = tout[2]-tout[1]

    # find all hi-res values bounded by tiepoint 1 and 2
    inrange = x -> ( tout[1] < x <  tout[2]) # found on plot
    iblock = findall(inrange,collect(tin))

    nblock = length(iblock)
    Eblock = zeros(nblock,2)
    for ii = 1:nblock
        Eblock[ii,1] = -(tin[iblock[ii]] - tout[2])/Δtout
        Eblock[ii,2] = (tin[iblock[ii]] - tout[1])/Δtout
    end

    # put the blocks together in a big sparse matrix.
    Eout2in = zeros(ntin+56,ntout) # make ntin too long and then trim later

    # for a real function, would want to fix this.
    ival = 3 # first two 6-hourly values off range
    jval = 1

    for ii = 1:ntout-1
        Eout2in[ival:ival+nblock-1,jval:jval+1] = Eblock
        ival += nblock
        jval += 1
    end

    # some trimming must be done.
    Eout2in = Eout2in[1:ntin,:]
    ETE = Eout2in'*Eout2in
    Fin2out = ETE\(Eout2in')

    Eout2in = convert(Array{Float32,2},Eout2in)
    Fin2out = convert(Array{Float32,2},Fin2out)

    return Eout2in, Fin2out
end

"""
    matrixfilter(F,froot,years,γ)
    writing it in a funny way to save computation
    issue with timeseries being read in different files
# Arguments
- `F`: filter in matrix form
- `froot`: filename root
- `years`: iterator for multiple files
- `γ`: GCM grid (meshArray type)
"""
function matrixfilter(F,froot,years,γ)

    # initialize known variables

    # use size of F to determine size of input/output
    nout = size(F,1)
    nin  = size(F,2)

    nyr = length(years)

    fname = froot*string(years[1])
    println("initialize θout")
    field = read_bin(fname,Float32,γ)

    # initialize θh
    # θout = 0.0f0 .*field[:,1:nout] # prone to have NaNs

    θout = MeshArray(γ,Float32,nout) # some nans here
    fill!(θout,0.0)
    #for tt= 1:nout
        # initialize a sub-optimal way
    #    θout[:,tt]=γ.read(tmp1,MeshArray(γ,Float32))
    #end
    println("initialization finished")

    istart = 1
    for tt = 1:nyr
        fname = froot*string(years[tt])
        println(fname)
        field = read_bin(fname,Float32,γ)

        # may need to keep track of indices
        nt = size(field,2)
        iend = istart+nt-1
        @time annualcontrib!(θout,F[:,istart:iend],field)
        istart = iend+1
    end

    return θout
end

"""
    function matrixfilter_annualcontrib!(θout,F,θin)
# Arguments
- `θout`: filtered values
- `F`: filter function as a matrix
- `θin`: input values at all times t (e.g., θ(t))
"""
function annualcontrib!(θout,F,θin)
    nt = size(θin,2)
    for tt = 1:nt
         # igood = findall(x -> abs(x)>1e-8,F[:,tt])
         # if !isempty(igood)
         #     matmul2!(θout[:,igood],F[igood,tt],θin[:,tt])
         # end
        columnscale!(θout,F[:,tt],θin[:,tt])
    end
end

"""
    function columnscale!(product,M,flux)
- `product::MeshArrays.gcmarray{Float32,2,Array{Float32,2}}`
- `M::Array{Float32,1}`
- `flux::MeshArrays.gcmarray{Float32,1,Array{Float32,2}}`
"""
function columnscale!(product::MeshArrays.gcmarray{Float32,2,Array{Float32,2}},M::Array{Float32,1},flux::MeshArrays.gcmarray{Float32,1,Array{Float32,2}})

    nM = size(M,1)  # matrix size M
    nQ  = size(flux,1) # repeat matmul Q times

     for qq = 1:nQ # inner product over faces
         for mm = 1:nM
             if abs(M[mm]) > 1e-8
                 # igood = findall(x -> abs(x)>1e-8,F[:,tt])
                 # if !isempty(igood)
                 #     matmul2!(θout[:,igood],F[igood,tt],θin[:,tt])
                 # end
                 product[qq,mm] += flux[qq] * M[mm]
             end
        end
    end
end

"""
    function columnscale!(product,M,flux)
- `product::MeshArrays.gcmarray{Float32,2,Array{Float32,2}}`
- `M::Array{Float32,1}`
- `flux::MeshArrays.gcmarray{Float32,2,Array{Float32,2}}`
"""
function columnscale!(product::MeshArrays.gcmarray{Float32,2,Array{Float32,2}},M::Array{Float32,2},flux::MeshArrays.gcmarray{Float32,2,Array{Float32,2}})

    nM = size(M,1)  # matrix size M
    nQ  = size(flux,1) # repeat matmul Q times
    nN  = size(M,2) # repeat matmul Q times

    for mm = 1:nM
        for qq = 1:nQ # inner product over faces
            for nn = 1:nN # inner product over faces
                product[qq,mm] += flux[qq,nn] * M[mm,nn]
            end
        end
    end
end

"""
    matrixspray(F,rmfield,frootin,frootout,years,γ)

    writing it in a funny way to save computation
    issue with timeseries being read in different files
# Arguments
- `F`: spray operator in matrix form
- `rmfield`: field to be removed from files
- `frootin`: filename root of input field
- `frootout`: filename root of output field
- `years`: iterator for multiple files
- `γ`: GCM grid (meshArray type)
"""
function matrixspray(F,rmfield,frootin,frootout,years,γ)

    # initialize known variables

    # use size of F to determine size of input/output
    nout = size(F,1)
    nin  = size(F,2)

    nyr = length(years)

    fnamein = frootin*string(years[1])
    println("initialize θout")
    field = read_bin(fnamein,Float32,γ)

    istart = 1
    nseries = []
    for tt = 1:nyr
        fnamein = frootin*string(years[tt])
        println("reading file "*fnamein)
        field = read_bin(fnamein,Float32,γ)

        # may need to keep track of indices
        push!(nseries,size(field,2))
        iend = istart+nseries[end]-1
        @time annualcontrib!(field,F[istart:iend,:],rmfield)
        istart = iend+1

        fnameout = frootout*string(years[tt])
        println("saved file "*fnameout)
        write(fnameout,field)
    end
end


text:
using EchoviewEcs

using Test

filename = joinpath(dirname(@__FILE__),
                  "data/sample.ecs")

calibrations = load(filename)

@test length(calibrations) == 4
@test calibrations[1]["Frequency"] == 18.0
@test calibrations[2]["Frequency"] == 38.0
@test calibrations[3]["Frequency"] == 120.0
@test calibrations[4]["Frequency"] == 200.0
@test calibrations[4]["EK60SaCorrection"] == -1.54f0

@test calibrations[1]["SoundSpeed"] == 1520f0 # Specified in SourceCal
@test calibrations[2]["SoundSpeed"] == 1541.21f0 # Inheritted from FileSet


text:
using VertexEliminationOrder
using Test
using Graphs
using SparseArrays

# test augment_flow!
G = smallgraph("house")
n = nv(G)
cap = SparseMatrixCSC{Int64, Int64}(adjacency_matrix(G))
flow = spzeros(Int64, n, n)
expected_flow = spzeros(Int64, n, n)
expected_flow[5, 4] = 1
expected_flow[4, 5] = -1
expected_flow[4, 2] = 1
expected_flow[2, 4] = -1
expected_flow[2, 1] = 1
expected_flow[1, 2] = -1
@info "Testing augment_flow! with house graph"
@test augment_flow!(flow, cap, G, 5, 1) == 1
@test flow == expected_flow
expected_flow[5, 3] = 1
expected_flow[3, 5] = -1
expected_flow[3, 1] = 1
expected_flow[1, 3] = -1
@test augment_flow!(flow, cap, G, 5, 1) == 1
@test flow == expected_flow
# no more augmenting path available
@test augment_flow!(flow, cap, G, 5, 1) == 0

# undo the last flow augmentation
flow[5, 3] = 0
flow[3, 5] = 0
flow[3, 1] = 0
flow[1, 3] = 0

# test forward growing
SR = falses(5); SR[5] = true
expected_SR = falses(5); expected_SR[[1, 3, 4, 5]] .= true
forward_grow!(SR, G, flow, cap)
@test SR == expected_SR

# test backward growing
TR = falses(5); TR[1] = true
expected_TR = falses(5); expected_TR[[1, 3, 4, 5]] .= true
@info "Testing forward_grow! (backward) with house graph"
forward_grow!(TR, G, flow, cap)
@test TR == expected_TR

# test piercing_node
dist = floyd_warshall_shortest_paths(G).dists
cut = [4=>2, 4=>3, 5=>3]
SR=falses(5); SR[[4, 5]] .= true
TR=falses(5); SR[[1, 3]] .= true
@info "Testing piercing_node with house graph"
@test piercing_node(cut, SR, TR, true, dist) == 2

#test flowcutter
@info flowcutter(G, 5, 1)


text:
function make_copy(device::PSY.Component, name::AbstractString)
    new_device = deepcopy(device)
    IS.clear_forecasts!(new_device)
    set_device_name!(new_device, name)
    set_uuid!(new_device)
    IS.copy_forecasts!(device, new_device)
    IS.set_time_series_storage!(new_device, nothing)
    return new_device
end

function set_uuid!(device::PSY.Component)
    device.internal.uuid = UUIDs.uuid4()
    return
end

function set_device_name!(device::PSY.Component, name::AbstractString)
    device.name = name
    return
end

function add_copies!(
    sys::PSY.System,
    device_list::Dict{Type{<:PSY.Component}, Array{AbstractString}},
)
    buses = collect(PSY.get_components(PSY.Bus, sys))
    for (type, device_names) in device_list
        if !isempty(device_names)
            for (ix, name) in enumerate(device_names)
                bus = buses[rand(1:3)]
                device = PSY.get_components_by_name(supertype(type), sys, name)[1]
                new_device = make_copy(device, name * "_copy_$(ix)")
                PSY.set_bus!(new_device, bus)
                services = PSY.get_services(new_device)
                _services = deepcopy(PSY.get_services(new_device))
                while !isempty(services)
                    PSY._remove_service!(new_device, services[1])
                end
                PSY.add_component!(sys, new_device)
                for service in _services
                    PSY.add_service!(new_device, service)
                end

            end
        end
    end
    return
end

function ConfigDataContainer()
    return Dict{Type{<:PSY.Component}, Array{AbstractString}}(
        PSY.ThermalStandard => Vector{String}(),
        PSY.RenewableFix => Vector{String}(),
        PSY.RenewableDispatch => Vector{String}(),
        PSY.GenericBattery => Vector{String}(),
        PSY.HydroDispatch => Vector{String}(),
        PSY.HydroEnergyReservoir => Vector{String}(),
    )
end

function update_config(future_config_dict, future, sample_count)
    for config in future
        future_config_dict[sample_count] = ConfigDataContainer()
        future_config_dict[sample_count][PSY.ThermalStandard] = config[1]
        future_config_dict[sample_count][PSY.HydroEnergyReservoir] = config[2]
        future_config_dict[sample_count][PSY.RenewableDispatch] = config[3]
        future_config_dict[sample_count][PSY.GenericBattery] = config[4]
        sample_count += 1
    end
    return future_config_dict, (sample_count - 1)
end