text:
using LibBSON, Mongo, Base.Test

@testset "Mongo" begin
    client = MongoClient()
    collection = MongoCollection(client, "foo", "bar")
    oid = BSONOID()

    @testset "variables" begin
      @test collection.client == client
      @test collection.db == "foo"
      @test collection.name == "bar"
    end

    @testset "insert" begin
        insert(collection, ("_id" => oid, "hello" => "before"))
        @test count(collection, ("_id" => oid)) == 1
        @test count(collection) == 1
        for item in find(collection, ("_id" => oid), ("_id" => false, "hello" => true))
            @test dict(item) == Dict("hello" => "before")
        end
    end

    @testset "update" begin
        update(
            collection,
            ("_id" => oid),
            set("hello" => "after")
            )

        @test count(collection, ("_id" => oid)) == 1
        for item in find(collection, ("_id" => oid), ("_id" => false, "hello" => true))
            @test dict(item) == Dict("hello" => "after")
        end
    end

    @testset "command_simple" begin
        reply = command_simple(
            client,
            "foo",
            Dict(
               "count" => "bar",
               "query" => Dict("_id" => oid))
            )
        @test reply["n"] == 1
    end

    @testset "delete" begin
        delete(
            collection,
            ("_id" => oid)
            )
        @test count(collection, ("_id" => oid)) == 0
        @test count(collection) == 0
    end
end

@testset "Mongo: bad host/port" begin
    client = MongoClient("bad-host-name", 9999)
    collection = MongoCollection(client, "foo", "bar")
    @test_throws ErrorException insert(collection, ("foo" => "bar"))
end

@testset "Query building helpers" begin
    client = MongoClient()
    ppl = MongoCollection(client, "foo", "ppl")
    person(name, age) = insert(ppl, ("name" => name, "age" => age))
    person("Tim", 25)
    person("Jason", 21)
    person("Jim", 87)

    @testset "orderby" begin
        @test first(find(ppl, (query(), orderby("age" => -1))))["name"] == "Jim"
        @test first(find(ppl, (query(), orderby("age" => 1))))["name"] == "Jason"
    end

    @testset "gt and lt" begin
        @test first(find(ppl, query("age" => lt(25))))["name"] == "Jason"
        @test first(find(ppl, query("age" => gt(50))))["name"] == "Jim"
    end

    @testset "in and nin" begin
        @test first(find(ppl, query("age" => in([21]))))["name"] == "Jason"
        @test first(find(ppl, query("age" => nin([21,25]))))["name"] == "Jim"
    end

    @testset "eq and ne" begin
        @test first(find(ppl, query("age" => eq(21))))["name"] == "Jason"
        @test first(find(ppl, query("age" => ne(87))))["name"] != "Jim"
    end

    @testset "update with operator" begin
        update(ppl, ("age" => 87), set("age" => 88))
        @test first(find(ppl, query("name" => "Jim")))["age"] == 88
    end

    delete(ppl, ())
end


text:
# + is not defined between SingleVariable
function MP.polynomial(p::GramMatrix{MOI.SingleVariable})
    Q = convert(Vector{MOI.ScalarAffineFunction{Float64}}, p.Q.Q)
    MP.polynomial(GramMatrix(SymMatrix(Q, p.Q.n), p.x))
end
function MP.polynomial(p::GramMatrix{F}) where {F <: MOI.AbstractFunction}
    MP.polynomial(p, MOIU.promote_operation(+, Float64, F, F))
end

function primal_value(model, p::GramMatrix{MOI.SingleVariable})
    # TODO [perf] use MOI typed mapped array
    Q = MOI.get(model, MOI.VariablePrimal(),
                MOI.VariableIndex[sv.variable for sv in p.Q.Q])
    return GramMatrix(SymMatrix(Q, p.Q.n), p.x)
end

ext: jl


text:
module ReferenceTests

using Test
using ImageCore
using Distances
using FileIO
using ImageInTerminal
using SHA
using DeepDiffs
using Random

export
    @withcolor,
    @io2str,
    @test_reference,
    psnr_equality

include("utils.jl")
include("test_reference.jl")
include("fileio.jl")
include("equality_metrics.jl")
include("render.jl")

if Base.VERSION >= v"1.4.2" && ccall(:jl_generating_output, Cint, ()) == 1
    precompile(Tuple{Core.kwftype(typeof(test_reference)),NamedTuple{(:by,), Tuple{typeof(isequal)}},typeof(test_reference),String,Matrix{RGB{Float64}}})  || @warn "Failure to precompile ReferenceTests.test_reference(String,Matrix{RGB})"
    precompile(Tuple{Core.kwftype(typeof(test_reference)),NamedTuple{(:by,), Tuple{typeof(isequal)}},typeof(test_reference),String,Matrix{RGBA{Float64}}}) || @warn "Failure to precompile ReferenceTests.test_reference(String,Matrix{RGBA})"
    if isdefined(FileIO, :action)
        precompile(Tuple{typeof(test_reference),File{DataFormat{:PNG},String},Matrix{RGB{Float64}},Function,Nothing})  || @warn "Failure to precompile ReferenceTests.test_reference(File{format\"PNG\",String},Matrix{RGB})"
        precompile(Tuple{typeof(test_reference),File{DataFormat{:PNG},String},Matrix{RGBA{Float64}},Function,Nothing}) || @warn "Failure to precompile ReferenceTests.test_reference(File{format\"PNG\",String},Matrix{RGBA})"
    else
        precompile(Tuple{typeof(test_reference),File{DataFormat{:PNG}},Matrix{RGB{Float64}},Function,Nothing})  || @warn "Failure to precompile ReferenceTests.test_reference(File{format\"PNG\"},Matrix{RGB})"
        precompile(Tuple{typeof(test_reference),File{DataFormat{:PNG}},Matrix{RGBA{Float64}},Function,Nothing}) || @warn "Failure to precompile ReferenceTests.test_reference(File{format\"PNG\"},Matrix{RGBA})"
    end
end

end # module


text:
# equal probability constructor
function MPS{T}(lx::Int, d::Int=2;
                noise::Float64=0.0,
                seed=1911) where {T<:Number}

    rng = MersenneTwister(seed);
    As = [ sqrt(1/d) * ones(T, 1, d, 1) + noise * randn(rng, T, 1, d, 1)
           for i=1:lx ]
    mps = MPS{T}(As, 0)
    center_at!(mps, 1)
    mps
end

# constructor with intial configuration vector
function MPS{T}(lx::Int, d::Int,
                initconf::Vector{Int}) where {T<:Number}
    all((0 .<= initconf) .& (initconf .< d)) ||
        error("conf vector out of bounds!")

    mps = MPS{T}([zeros(T, 1, d, 1) for i=1:lx ], 0)
    for site=1:lx
        mps.As[site][1, initconf[site]+1, 1] = 1.
    end
    center_at!(mps, lx)
    mps
end

# constructor with intial configuration vector and noise
function MPS{T}(lx::Int, d::Int,
                initconf::Vector{Int},
                maxdim::Int,
                noise::Float64;
                seed=1911) where {T<:Number}
    all((0 .<= initconf) .& (initconf .< d)) ||
        error("conf vector out of bounds!")

    dims = ones(Int, lx+1)
    ### this can be wrapped into a function itself
    for l = 2:lx
        dims[l] = min(maxdim, dims[l-1]*d)
    end
    for l = lx:-1:1
        dims[l] = min(dims[l], dims[l+1]*d)
    end
    ### till here!

    rng = MersenneTwister(seed);
    mps = MPS{T}([ noise * randn(rng, T, dims[i], d, dims[i+1]) for i=1:lx ], 0)
    for site=1:lx
        mps.As[site][1, initconf[site]+1, 1] = 1.
    end
    center_at!(mps, lx)
    mps
end

#constructor from a ketstate given in Ising basis
function MPS(lx::Int, d::Int,
             ketstate::Vector{T};
             svtruncation::Bool=false) where {T<:Number}

    @assert length(ketstate) == d^lx
    mps = MPS{T}()
    #As = Array{T,3}[]

    dims = zeros(Int64, lx+1)
    dims[1] = 1
    rdim = d^(lx-1)

    # state_matrix is now : dims[1]*d x rdim
    state_matrix = transpose(reshape(ketstate, rdim, dims[1]*d))

    for link=2:lx-1
        fact = svd(state_matrix, full=false)

        if svtruncation
            S, n, ratio =  truncate(fact.S, threshold=1.e-15)
        else
            S, n, ratio = fact.S, length(fact.S), 1.
        end
        dims[link] = n

        # U is : dims[link-1] * d, dims[link]
        U = fact.U[:,1:n]
        push!(mps, reshape(U, dims[link-1], d, dims[link]))

        # S*Vt is : dims[link] x rdim where rdim is the new rdim x d
        rdim = div(rdim , d)
        state_matrix = reshape(permutedims(
            reshape(Diagonal(S)*fact.Vt[1:n,:], dims[link], rdim, d),
            [1,3,2]), dims[link]*d, rdim)
        # final results is dims[link]*d x rdim
    end

    fact = svd(state_matrix, full=false)

    if svtruncation
        S, n, ratio =  truncate(fact.S)
    else
        S, n, ratio = fact.S, length(fact.S), 1.
    end

    dims[lx] = n

    U = fact.U[:,1:n]
    push!(mps, reshape(U, dims[lx-1], d, dims[lx]))

    dims[lx+1] = 1
    push!(mps, permutedims(reshape(Diagonal(S) * fact.Vt[1:n,:],
                                   dims[lx], dims[lx+1], d), [1,3,2]))

    return mps
end

"""
    mps2ketstate(mps)

make a ketstate from a `mps` by multiplication the matrices
corresponding to each Ising configuration. Note that this function is
very expensive and is only meant for testing.

"""
function mps2ketstate(mps::MPS{T}) where {T<:Number}

    lx = length(mps)
    ds = [dim(sitespace(mps, l)) for l=1:lx]

    ketstate = zeros(T, prod(ds))

    for ketindex::Int = 0:prod(ds)-1
        resolve::Int = ketindex
        amplitude = T[1]
        ## TODO: find a better way!
        for i::Int64 = lx:-1:1
            d = ds[i]
            A = mps.As[i]
            amplitude = reshape(A[:,(resolve % d) + 1,:],
                                size(A, 1), size(A, 3)) * amplitude
            resolve = div(resolve, d)
        end
        ketstate[ketindex+1] = amplitude[1]
    end
    ketstate
end

# function randmps(T::Type{<:RLorCX}, rng::AbstractRNG, lx::Int, d::Int, maxdim::Int)
#     dims = Vector{Int}(undef, lx+1)
#     dims[1] = 1
#     dims[lx+1] = 1

#     ### this can be wrapped into a function itself
#     for l = 2:lx
#         dims[l] = min(maxdim, dims[l-1]*d)
#     end
#     for l = lx:-1:1
#         dims[l] = min(dims[l], dims[l+1]*d)
#     end
#     ### till here!

#     matrices = Vector{Array{T, 3}}(undef, lx)
#     for l = 1:lx
#         Dr = dims[l+1]
#         Dl = dims[l]
#         matrices[l] = reshape(randisometry(T,rng, Dl, d*Dr), Dl, d, Dr)
#     end
#     return MatrixProductState{T}(lx, d, dims, matrices, lx)
# end

# ### conversions
# ###############

# function convert(::Type{MatrixProductState{ComplexF64}},
#                  mps::MatrixProductState{Float64})
#     MatrixProductState{ComplexF64}(mps.lx, mps.d, mps.dims,
#                                    convert(Vector{Array{ComplexF64,3}}, mps.matrices),
#                                    mps.center)
# end

# ### TOOLS
# #########

# function truncate(svector::Vector{Float64};
#                   maxdim::Int64=length(svector),
#                   threshold::Float64=1.e-15)

#     ## NOTE: S is assumed be all non-negative and sorted in
#     ## descending order
#     n = min(maxdim, sum(svector .> svector[1]*threshold))
#     return svector[1:n], n, sum(svector[1:n])/sum(svector)
# end

# function normalize!(mps::MatrixProductState{T}) where{T<:RLorCX}
#     A = mps.matrices[mps.center]
#     if mps.center > div(mps.lx, 2)
#         U,S,Vt = svd(reshape(A, prod(size(A)[1:2]), size(A,3)))
#     else
#         U,S,Vt = svd(reshape(A, size(A,1), prod(size(A)[2:3])))
#     end
#     normalize!(S)
#     mps.matrices[mps.center] = reshape(U * S * Vt, size(A))
#     S
# end

# """
#     canonicalize_at!(matrices, center)

# canonicalize at `center` from the both ends of an MPS. This function
# is mainly used during the initialization so one should not truncate
# during the canonicalization process!

# """
# function canonicalize_at!(matrices::Vector{Array{T, 3}},
#                           center::Int64) where {T<:RLorCX}
#     lx = length(matrices)
#     @assert center > 0 && center < lx + 1

#     for site=1:center-1
#         isometrize_push_right!(matrices, site, svtruncation=false)
#     end

#     for site=lx:-1:center+1
#         isometrize_push_left!(matrices, site, svtruncation=false)
#     end
# end

# """
#    isometrize_push_right!(matrices, site)

# perform a right step of the canonicalization procedure of a vector of
# matrices at `site`. This is to perform an SVD on the matrices at
# `site`, and then multiply the singluar values and right unitary
# matricx ``V^{†}`` to matrices of `site+1`. This ensures the matrices
# of `site` are left isometric, because U^{†}U = I.

# """
# function isometrize_push_right!(matrices::Vector{Array{T, 3}},
#                                 site::Int64;
#                                 svtruncation::Bool=false) where {T<:RLorCX}
#     lx = length(matrices)
#     if site < lx
#         a = matrices[site]
#         dims = size(a)
#         fact = svd(reshape(a, dims[1]*dims[2], dims[3]), full=false)

#         if svtruncation
#             S, n, ratio =  truncate(fact.S, threshold=1.e-15)
#             U = fact.U[:,1:n]
#             Vt = fact.Vt[1:n,:]
#         else
#             S, n, ratio = fact.S, length(fact.S), 1.
#             U = fact.U
#             Vt = fact.Vt
#         end

#         matrices[site] = reshape(U, dims[1], dims[2], n)

#         @tensor matrices[site+1][i,d,j] := (Diagonal(S) * Vt)[i, k] *
#             matrices[site+1][k,d,j]
#     end
#     return nothing
# end

# """
#     isometrize_push_left!(matrices, site)

# perform a left step of the canonicalization procedure of a vector of
# matrices at `site`. This is perform an SVD on the matrices at `site`,
# and then multiply the left unitary matricx ``U`` and singular values
# to matrices of `site-1`. This ensures the matrices of `site` are
# right isometric, because VV^{†} = I.

# """
# function isometrize_push_left!(matrices::Vector{Array{T, 3}},
#                                site::Int64;
#                                svtruncation::Bool=false) where {T<:RLorCX}
#     if site > 0
#         a = matrices[site]
#         dims = size(a)
#         fact = svd(reshape(a, dims[1], dims[2] * dims[3]), full=false)

#         if svtruncation
#             S, n, ratio =  truncate(fact.S)
#             U = fact.U[:,1:n]
#             Vt = fact.Vt[1:n,:]
#         else
#             S, n, ratio = fact.S, length(fact.S), 1.
#             U = fact.U
#             Vt = fact.Vt
#         end

#         #matrices[site] = permutedims(reshape(Vt, n, dims[3], dims[2]), [1,3,2])
#         matrices[site] = reshape(Vt, n, dims[2], dims[3])

#         @tensor matrices[site-1][i,d,j] := matrices[site-1][i,d,k] *
#             (U * Diagonal(S))[k,j]
#     end
#     nothing
# end
# """
#     correctdims!(mps)

# correctly store the bond dimensions of the MPS. This is needed when
# some function manipulate the matrices of the MPS but don't see/fix the
# dims field.

# """
# function correctdims!(mps::MatrixProductState{T}) where {T<:RLorCX}
#     dims = Int64[]
#     push!(dims,size(mps.matrices[1])[1])
#     for n=1:mps.lx
#         push!(dims, size(mps.matrices[n])[3])
#     end
#     mps.dims = dims
#     nothing
# end

# """
#     move_center!(mps, new_center)

# move center of mps to a new location at `new_center`.

# """
# function move_center!(mps::MatrixProductState{T},
#                       new_center::Int64;
#                       svnormalize::Bool=false) where {T<:RLorCX}
#     lx= mps.lx
#     @assert new_center > 0 && new_center <= lx

#     center = mps.center
#     if new_center > center
#         for p=center:new_center-1
#             isometrize_push_right!(mps.matrices, p)
#         end
#     elseif new_center < center
#         for p=center:-1:new_center+1
#             isometrize_push_left!(mps.matrices, p)
#         end
#     end
#     correctdims!(mps)
#     mps.center = new_center
#     nothing
# end

## TODO: combine these with the symmetric cases!
# specific measure functions
function measure(mps::MPS{T},
                 op::Matrix{T}, l::Int) where {T}
    lx = length(mps)
    @assert (l <= lx) && (l > 0)
    @assert sitespace(mps, l) == space(op, 2) == dual(space(op, 1))

    center_at!(mps, l)
    A = mps.As[l]
    @tensor v = scalar((A[l,d',r] * conj(A)[l,d,r]) * op[d,d'])
    Dict{Int, T}(l=>v)
end

function measure(mps::MPS{T},
                 op::Matrix{T}) where {T}
    Vd = space(op, 1)
    @assert space(op, 2) == Vd

    lx = length(mps)
    values = Dict{Int, T}()
    center_at!(mps, 1)
    L = ones(T, 1, 1)

    for x = 1:lx
        A = mps.As[x]
        @tensor v = scalar((L[lu, ld] * A[lu, d', r]) * op[d, d'] * conj(A)[ld, d, r])
        values[x] = v
        @tensor L[ru, rd] := (L[lu, ld] * A[lu, d, ru]) * conj(A)[ld, d, rd]
    end
    values
end
function measure(mps::MPS{T1},
                 op::Matrix{T2}) where {T1,T2}
    measure(mps, convert(Matrix{T1}, op))
end

function measure(mps::MPS{T},
                 op1::Matrix{T},
                 op2::Matrix{T},
                 x1::Int,
                 x2::Int) where {T}
    lx = length(mps)
    @assert (0 < x1 < x2) && (x2 <= lx)

    Vd1 = space(op1, 1)
    @assert space(op1, 2) == Vd1
    Vd2 = space(op2, 1)
    @assert space(op2, 2) == Vd2

    values = Dict{NTuple{2, Int}, T}()
    center_at!(mps, x1)
    A = mps.As[x1]

    @tensor L[ru,rd] := A[l,d',ru] * op1[d,d'] * conj(A)[l,d,rd]
    for x=x1+1:x2-1
        A = mps.As[site]
        @tensor L[ru,rd] := L[lu, ld] * A[lu,d,ru] * conj(A)[ld,d,rd]
    end
    A = mps.As[x2]
    @tensor v = scalar((L[lu, ld] * A[lu,d',r] * conj(A)[ld,d,r]) * op2[d,d'])
    values[(x1, x2)] = v
end

function measure(mps::MPS{T},
                 op1::Matrix{T},
                 op2::Matrix{T}) where {T}
    lx = length(mps)
    Vd1 = space(op1, 1)
    @assert space(op1, 2) == dual(Vd1)
    Vd2 = space(op2, 1)
    @assert space(op2, 2) == dual(Vd2)

    values = Dict{NTuple{2,Int}, T}()
    for x1=1:lx-1
        center_at!(mps, x1)
        A = mps.As[x1]
        @tensor L[ru,rd] := A[l,d',ru] * op1[d,d'] * conj(A)[l,d,rd]
        for x=x1+1:lx-1
            A = mps.As[x]
            @tensor v = scalar((L[lu, ld] * A[lu,d',r] * conj(A)[ld,d,r]) * op2[d,d'])
            values[(x1, x)] = v
            @tensor L[ru,rd] := L[lu, ld] * A[lu,d,ru] * conj(A)[ld,d,rd]
        end
        A = mps.As[lx]
        @tensor v = scalar((L[lu, ld] * A[lu,d',r] * conj(A)[ld,d,r]) * op2[d,d'])
        values[x1, lx] = v
    end
    values
end

function measure(mps::MPS{T1},
                 op1::Matrix{T2},
                 op2::Matrix{T3}) where {T1,T2,T3}
    measure(mps, convert(Matrix{T1}, op1), convert(Matrix{T1}, op2))
end

function measure(mps::MPS{T},
                 mpo::MPO{T}) where {T}

    @assert length(mpo) == length(mps)
    values = Dict{Int, T}()
    L = ones(T, 1, 1, 1)
    for x=1:length(mps)
        A = mps.As[x]
        W = mpo.Ws[x]
        @tensor L[ru,rm,rd] := (L[lu,lm,ld] * A[lu,d',ru]) *
            W[lm, d, rm, d'] * conj(A)[ld,d,rd]
    end
    values[1] = L[1,1,1]
    values
end

function apply!(mps         :: MPS{T},
                op          :: Array{T,4},
                x           :: Int;
                maxdim      :: Int=bonddim(mps, x+1),
                pushto      :: Symbol=:R,
                svnormalize :: Bool=false) where {T}

    @boundscheck 0 < x < length(mps) || throw()

    if mps.center < x
        center_at!(mps, x)
    elseif mps.center > x+1
        center_at!(mps, x+1)
    end

    A1 = mps.As[x]
    A2 = mps.As[x+1]

    @tensor R[a,i,j,b] := op[i,j,k,l] * (A1[a,k,c] * A2[c,l,b])

    u,s,v = svdtrunc(reshape(R, prod(size(R)[1:2]), prod(size(R)[3:4])),
                     maxdim=maxdim, tol=1.e-14)

    svnormalize && normalize!(s)

    _space = space(R)
    if (pushto == :R)
        mps.As[x] = splitleg(u, 1, _space[1:2])
        mps.As[x+1] = splitleg(s*v, 2, _space[3:4])
        mps.center = x+1
    elseif (pushto == :L)
        mps.As[x] = splitleg(u*s, 1, _space[1:2])
        mps.As[x+1] = splitleg(v, 2, _space[3:4])
        mps.center = x
    else
        error("invalid push_to :", pushto)
    end
    mps
end

function apply!(mps         :: MPS{T1},
                op          :: Array{T2,4},
                x           :: Int;
                maxdim      :: Int=bonddim(mps, x+1),
                pushto      :: Symbol=:R,
                svnormalize :: Bool=false) where {T1,T2}

    T = promote_type(T1, T2)
    apply!(convert(MPS{T}, mps),
           convert(Array{T,4}, op),
           x, maxdim=maxdim, pushto=pushto, svnormalize=svnormalize)
end

function apply!(mps      :: MPS{T},
                op       :: Array{T, 2},
                x        :: Int) where {T}

    @boundscheck 0 < x < length(mps) || throw()

    center_at!(mps, x)

    A = mps.As[x]
    @tensor A[l,o,r] := A[l,o',r] * op[o,o']
    mps.As[x] = A
    mps
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
"""
    UnsupportedError <: Exception

Abstract type for error thrown when an element is not supported by the model.
"""
abstract type UnsupportedError <: Exception end

"""
    element_name(err::UnsupportedError)

Return the name of the element that is not supported
"""
function element_name end

function Base.showerror(io::IO, err::UnsupportedError)
    print(io, typeof(err), ": ", element_name(err),
          " is not supported by the the model")
    m = message(err)
    if Base.isempty(m)
        print(io, ".")
    else
        print(io, ": ", m)
    end
end

"""
    NotAllowedError <: Exception

Abstract type for error thrown when an operation is supported but cannot be
applied in the current state of the model.
"""
abstract type NotAllowedError <: Exception end

"""
    operation_name(err::NotAllowedError)

Return the name of the operation throwing the error in a gerund (i.e. -ing
form).
"""
function operation_name end

function Base.showerror(io::IO, err::NotAllowedError)
    print(io, typeof(err), ": ", operation_name(err),
          " cannot be performed")
    m = message(err)
    if Base.isempty(m)
        print(io, ".")
    else
        print(io, ": ", m)
    end
    print(io, " You may want to use a `CachingOptimizer` in `Automatic` mode",
          " or you may need to call `resetoptimizer!` before doing this",
          " operation if the `CachingOptimizer` is in `Manual` mode.")
end

"""
    message(err::Union{UnsupportedError, NotAllowedError})

Return a `String` containing a human-friendly explanation why the operation
is not supported/cannot be performed. It is printed in the error message if it
is not empty. By convention, it should be stored in the `message` field, it is
it the case, the `message` method does not have to be implemented.
"""
message(err::Union{UnsupportedError, NotAllowedError}) = err.message


text:
using UnicodePlots: barplot
using Base: show, print, popfirst!

"""
    AcctMap{Sz,Px,Oid,Aid}

Collection of open orders by account.

`(Sz,Px,Oid,Aid)` characterize the type of Order present in the `AcctMap`.
See documentation on [`Order`](@ref) for more information on the meaning of types.

The account map is implemented as a `Dict` containing `AVLTree`s.
    AcctMap{Sz,Px,Oid,Aid} = Dict{Aid,AVLTree{Oid,Order{Sz,Px,Oid,Aid}}}
The outer key is the account id, mapping to an `AVLTree` of `Order`s keyed by order id.
"""
AcctMap{Sz<:Real,Px<:Real,Oid<:Integer,Aid<:Integer} = Dict{
    Aid,AVLTree{Oid,Order{Sz,Px,Oid,Aid}}
}

# Add functions for adding and removing orders from account map
@inline function _add_order_acct_map!(
    acct_map::AcctMap{Sz,Px,Oid,Aid}, acct::Aid, order::Order{Sz,Px,Oid,Aid}
) where {Sz,Px,Oid,Aid}
    if !haskey(acct_map, acct)
        acct_map[acct] = AVLTree{Oid,Order{Sz,Px,Oid,Aid}}() # If account isn't registered, register
    end
    insert!(acct_map[acct], order.orderid, order) # Add order to account map
    return nothing
end

@inline function _delete_order_acct_map!(
    acct_map::AcctMap{Sz,Px,Oid,Aid}, acct::Aid, orderid::Oid
) where {Sz,Px,Oid,Aid}
    haskey(acct_map, acct) && delete!(acct_map[acct], orderid)
    return nothing
end

"""
    OrderBook{Sz,Px,Oid,Aid}

An `OrderBook` is a data structure containing __limit orders__ represented as objects of type `Order{Sz,Px,Oid,Aid}`.

See documentation on [`Order`](@ref) for more information on this type.

How to use `Orderbook`:
 - Initialize an empty limit order book as `OrderBook{Sz,Px,Oid,Aid}()`
 - __Submit__ or __cancel__ limit orders with [`submit_limit_order!`](@ref) and [`cancel_order!`](@ref).
 - Submit __market orders__ with [`submit_market_order!`](@ref)
 - Retrieve order book state information with `print` or `show` methods, as well as [`book_depth_info`](@ref), [`best_bid_ask`](@ref), [`volume_bid_ask`](@ref), [`n_orders_bid_ask`](@ref) and [`get_acct`](@ref)
 - Write book state to `csv` file with [`write_csv`](@ref).

"""
mutable struct OrderBook{Sz<:Real,Px<:Real,Oid<:Integer,Aid<:Integer}
    bid_orders::OneSidedBook{Sz,Px,Oid,Aid} # bid orders
    ask_orders::OneSidedBook{Sz,Px,Oid,Aid} # ask orders
    acct_map::AcctMap{Sz,Px,Oid,Aid} # Map from acct_id::Aid to AVLTree{order_id::Oid,Order{Sz,Px,Oid,Aid}}
    flags::Dict{Symbol,Any} # container for additional order book logic flags (not yet implemented)
    function OrderBook{Sz,Px,Oid,Aid}() where {Sz,Px,Oid,Aid}
        return new{Sz,Px,Oid,Aid}(
            OneSidedBook{Sz,Px,Oid,Aid}(; is_bid_side=true),
            OneSidedBook{Sz,Px,Oid,Aid}(; is_bid_side=false),
            AcctMap{Sz,Px,Oid,Aid}(),
            Dict{Symbol,Any}(:PlotTickMax => 5),
        )
    end
end

"""
    clear_book!(ob::OrderBook,n_keep::Int64=10)

Remove all orders beyond `n_keep ≥ 0` from the best bid and best ask.
When `n_keep==0`, all orders are cleared.

"""
function clear_book!(ob::OrderBook{Sz,Px,Oid,Aid}, n_keep::Int64=10) where {Sz,Px,Oid,Aid}
    (n_keep < 0) && error("$n_keep should be non-negative")
    # clear the bids
    cleared_bids = Vector{Order{Sz,Px,Oid,Aid}}()
    bids_to_clear = [abs(k) for (k, v) in ob.bid_orders.book]
    bids_to_clear = x -> last(x, max(length(x) - n_keep, 0))(bids_to_clear) # clear all but last n
    for px in bids_to_clear
        cleared_queue = _popat_queue!(ob.bid_orders, px)
        append!(cleared_bids, cleared_queue.queue)
    end
    # clear the asks
    cleared_asks = Vector{Order{Sz,Px,Oid,Aid}}()
    asks_to_clear = [abs(k) for (k, v) in ob.ask_orders.book]
    asks_to_clear = x -> first(x, max(length(x) - n_keep, 0))(asks_to_clear) # clear all but last n
    for px in asks_to_clear
        cleared_queue = _popat_queue!(ob.ask_orders, px)
        append!(cleared_asks, cleared_queue.queue)
    end
    return cleared_bids, cleared_asks
end

## Order book statistics functions

@inline function _sidebook_stats(sb::OneSidedBook, max_depth)
    # get book statistics until fixed price depth
    raw_info_list = [
        (pq.price[], pq.total_volume[], pq.num_orders[]) for
        (pk, pq) in Base.Iterators.take(sb.book, max_depth)
    ]
    # Compile and return as dict of vectors
    return Dict(
        :side => isbidbook(sb) ? :BID : :ASK,
        :price => [x[1] for x in raw_info_list],
        :volume => [x[2] for x in raw_info_list],
        :orders => [x[3] for x in raw_info_list],
    )
end

"""
    book_depth_info(ob::OrderBook, max_depth=5)

Returns prices, volumes and order counts at bid and ask in `ob::OrderBook`
until fixed depth `max_depth` as a nested `Dict`.
"""
function book_depth_info(ob::OrderBook, max_depth=5)
    return Dict(
        :BID => _sidebook_stats(ob.bid_orders, max_depth),
        :ASK => _sidebook_stats(ob.ask_orders, max_depth),
    )
end

"""
    best_bid_ask(ob::OrderBook)

Return best bid/ask prices in order book as a `Tuple`
"""
function best_bid_ask(ob::OrderBook{Sz,Px,Oid,Aid}) where {Sz,Px,Oid,Aid}
    return (
        ob.bid_orders.best_price::Union{Nothing,Px},
        ob.ask_orders.best_price::Union{Nothing,Px},
    )
end

"""
    volume_bid_ask(ob::OrderBook)

Return total bid and ask volume from order book as a `Tuple`.
"""
function volume_bid_ask(ob::OrderBook{Sz,Px,Oid,Aid}) where {Sz,Px,Oid,Aid}
    return (
        ob.bid_orders.total_volume::Union{Nothing,Sz},
        ob.ask_orders.total_volume::Union{Nothing,Sz},
    )
end

"""
    n_orders_bid_ask(ob::OrderBook)

Return total number of orders on each side of order book, returned as a `Tuple`
"""
function n_orders_bid_ask(ob::OrderBook)
    return (
        ob.bid_orders.num_orders::Union{Nothing,Int32},
        ob.ask_orders.num_orders::Union{Nothing,Int32},
    )
end

# Order Book utility functions

"""
    order_types(::Order{Sz,Px,Oid,Aid})
    order_types(::OneSidedBook{Sz,Px,Oid,Aid})
    order_types(::OrderBook{Sz,Px,Oid,Aid})

Returns `(Sz,Px,Oid,Aid)`.
"""
order_types(::Order{Sz,Px,Oid,Aid}) where {Sz,Px,Oid,Aid} = Sz, Px, Oid, Aid
order_types(::OneSidedBook{Sz,Px,Oid,Aid}) where {Sz,Px,Oid,Aid} = Sz, Px, Oid, Aid
order_types(::OrderBook{Sz,Px,Oid,Aid}) where {Sz,Px,Oid,Aid} = Sz, Px, Oid, Aid

"""
    get_acct(ob::OrderBook{Sz,Px,Oid,Aid},acct_id::Aid)

Return all open orders assigned to account `acct_id`
"""
function get_acct(ob::OrderBook{Sz,Px,Oid,Aid}, acct_id::Aid) where {Sz,Px,Oid,Aid}
    return get(ob.acct_map, acct_id, nothing)
end

"""
    ask_orders(sb::OneSidedBook)

Return an iterator over all __ask orders__ by price/time priority.
"""
ask_orders(ob::OrderBook) = Iterators.flatten(q for (k,q) in ob.ask_orders.book)

"""
    bid_orders(sb::OneSidedBook)

Return an iterator over all __bid orders__ by price/time priority.
"""
bid_orders(ob::OrderBook) = Iterators.flatten(q for (k,q) in ob.bid_orders.book)

"""
    write_csv(
        io::IO,
        ob::OrderBook;
        row_formatter = _order_to_csv,
        header = "TRD,ID,SIDE,SIZE,PX,ACCT",
        )

Write OrderBook `ob` to an IO stream into `csv` format where each row corresponds to an order
The formatting for each row is given by the function argument `row_formatter(::Order)::String`.
The `csv` header can be provided as an argument where setting it to `nothing` writes no header.
"""
function write_csv(
    io::IO, ob::OrderBook; row_formatter=_order_to_csv, header="TRD,ID,SIDE,SIZE,PX,ACCT"
)
    cnt = 0
    !is_nothing(write_header) && (cnt += write(io, header, '\n'))
    # write all of the bids
    for (pk, pq) in ob.bid_orders.book
        for ord in pq.queue
            cnt += write(io, _order_to_csv(ord), '\n')
        end
    end
    # write all of the asks
    for (pk, pq) in ob.ask_orders.book
        for ord in pq.queue
            cnt += write(io, _order_to_csv(ord) * '\n')
        end
    end
    return cnt
end

@inline function _print_book_barplot(
    io::IO, ob::OrderBook{Sz,Px,Oid,Aid}
) where {Sz,Px,Oid,Aid}
    # Get book info
    max_depth = ob.flags[:PlotTickMax]
    sb_info = LimitOrderBook.book_depth_info(ob, max_depth)
    all_prices = [sb_info[:BID][:price]; sb_info[:ASK][:price]]

    println(io, "\n Order Book histogram (within $max_depth ticks of center):\n")

    if !isempty(sb_info[:BID][:volume])
        # Get max price str length
        max_len = maximum(length(string.(all_prices)))

        bid_plt = barplot(
            lpad.(string.(reverse(sb_info[:BID][:price])), max_len, " "),
            reverse(sb_info[:BID][:volume]);
            color=:red,
            ylabel=":BID",
            border=:none,
            padding=0,
        )

        println(io, bid_plt, '\n')
    else
        print(io, "\n    :BID   <empty>\n")
    end

    if !isempty(sb_info[:ASK][:volume])
        # Get max price str length
        max_len = maximum(length(string.(all_prices)))

        ask_plt = barplot(
            lpad.(string.(sb_info[:ASK][:price]), max_len, " "),
            sb_info[:ASK][:volume];
            ylabel=":ASK",
            border=:none,
            padding=0,
        )

        println(io, ask_plt)
    else
        print(io, "\n    :ASK   <empty>\n")
    end

    return nothing
end

function _print_book_info(io::IO, ob::OrderBook{Sz,Px,Oid,Aid}) where {Sz,Px,Oid,Aid}
    return print(
        io,
        "OrderBook{Sz=$Sz,Px=$Px,Oid=$Oid,Aid=$Aid} with properties:\n",
        "  ⋄ best bid/ask price: $(best_bid_ask(ob))\n",
        "  ⋄ total bid/ask volume: $(volume_bid_ask(ob))\n",
        "  ⋄ total bid/ask orders: $(n_orders_bid_ask(ob))\n",
        "  ⋄ flags = $([ k => v for (k,v) in ob.flags])",
    )
end

Base.print(io::IO, ob::OrderBook) = _print_book_info(io, ob)

function Base.show(io::IO, ::MIME"text/plain", ob::OrderBook)
    println(io, ob)
    _print_book_barplot(io, ob)
    return nothing
end

function Base.show(io::IO, ob::OrderBook)
    println(io, ob)
    return nothing
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -4.572700962608518 <= q <= 2.439967378874661)

                     add_NL_constraint(m, :(swish(-0.9934015186054679 + 0.8375762785867722*$(x[1]) + 0.668300383494036*$(x[2])) + swish(0.08494880908446678 + -0.07388095057292743*$(x[1]) + -0.9863741327473892*$(x[2])) + swish(-0.1579140823459273 + 0.7305802491259836*$(x[1]) + -0.2096221762144812*$(x[2])) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
export deepmap_parser

"""
    deepmap_parser(f::Function[, mem::AbstractDict=IdDict()], x::CombinedParser,a...;kw...)

Perform a deep transformation of a `x`.

Default method
1. Returns cached result if `haskey(x, mem)` to avoid infinite recursion.
2. construct deep transformation `dt = _deepmap_parser(f, mem, x, a...; kw...)`
3. cache and return `f(dt, a...; kw...)`

Used for [`log_names`](@ref).

## For a new `CombinedParser`,
define either `deepmap_parser` or `_deepmap_parser`.

## For a parser transformation `f`,
define either custom
- `deepmap_parser(::typeof(f),...)` (see example implementation [`substitute`](@ref))
- construction method `_deepmap_parser(::typeof(f),...)`  (see example implementation [`caseless`](@ref))
- leaf method `f` (see example implementation [`deepmap`](@ref))
"""
deepmap_parser(f,x::CombinedParser, a...;kw...) =
    deepmap_parser(f,IdDict(),x,a...;kw...)

deepmap_parser(f,mem::AbstractDict, x, a...; kw...) =
    get!(mem,x) do
        dt = _deepmap_parser(f, mem, x, a...; kw...)
        f(dt, a...; kw...)
    end

_deepmap_parser(f, mem::AbstractDict, x::LeafParser,a...;kw...) = x

_deepmap_parser(f, mem::AbstractDict, x::Union{AtStart,AtEnd,Always,Never}, a...; kw...) = x

_deepmap_parser(f,mem::AbstractDict, x::ConstantParser,a...;kw...) = x

"""
    deepmap_parser(f,mem::AbstractDict,x,a...;kw...)

Perform a deep transformation of a CombinedParser.

!!! note
    For a custom parser `P<:CombinedParser` with sub-parsers, provide a method
    ```julia
    CombinedParsers._deepmap_parser(f,mem::AbstractDict,x::P,a...;kw...) =
         ## construct replacement, e.g. if P <: WrappedParser
         P(deepmap_parser(f,mem,x.parser,a...;kw...))
    ```
"""
_deepmap_parser(f,mem::AbstractDict,x::CombinedParser,a...;kw...) =
    error("""
    For a custom parser `$(typeof(x))` with sub-parsers, provide a method
    ```julia
    CombinedParsers._deepmap_parser(f,mem::AbstractDict,x::$(typeof(x)),a...;kw...) =
          ## construct replacement, e.g. if P <: WrappedParser
          $(typeof(x))(deepmap_parser(f,mem,x.parser,a...;kw...))
    ```
""")

_deepmap_parser(f,mem::AbstractDict,x::Atomic,a...;kw...) =
    Atomic(deepmap_parser(f,mem,x.parser,a...;kw...))

_deepmap_parser(f,mem::AbstractDict,x::Lazy,a...;kw...) =
    Lazy(deepmap_parser(f,mem,x.parser,a...;kw...))

_deepmap_parser(f,mem::AbstractDict,x::Sequence,a...;kw...) =
    Sequence( ( deepmap_parser(f,mem,p,a...;kw...)
                for p in x.parts)... )

_deepmap_parser(f,mem::AbstractDict,x::Optional,a...;kw...) =
    Optional(deepmap_parser(f,mem,x.parser,a...;kw...);
             default=x.default)

_deepmap_parser(f,mem::AbstractDict,x::NamedParser,a...;kw...) =
    NamedParser(x.name,deepmap_parser(f,mem,x.parser,a...;kw...))

_deepmap_parser(f,mem::AbstractDict,x::Repeat,a...;kw...) =
    Repeat(x.range,
           deepmap_parser(f,mem,x.parser,a...))

_deepmap_parser(f,mem::AbstractDict,x::MappedSequenceParser,a...;kw...) =
    MappedSequenceParser(x.f,deepmap_parser(f,mem,x.parser,a...;kw...))

_deepmap_parser(f,mem::AbstractDict,x::SideeffectParser,a...;kw...) =
    SideeffectParser(
        x.effect,
        deepmap_parser(f,mem,x.parser,a...;kw...),
        x.args...)

_deepmap_parser(f,mem::AbstractDict,x::FlatMap,a...;kw...) =
    FlatMap{result_type(x)}(
        x.right,# v -> deepmap_parser(f,mem,x.right(v),a...;kw...),
        deepmap_parser(f,mem,x.left,a...;kw...))

@inline deepmap_parser(f,mem::AbstractDict,x::Either,a...;kw...) =
    get!(mem,x) do
        f(deepmap_either(f,mem,x,a...;kw...),a...;kw...)
    end

"""
    strip_either1(x::CombinedParser)

Replace all `Either` parsers with one option with that option.

Used in 2-stage [`substitute`](@ref) (stage 1: collect for recursion, stage 2: simplify).
"""
strip_either1(x::CombinedParser) = deepmap_parser(_strip_either1, x)
_strip_either1(x::CombinedParser) = x
deepmap_parser(::typeof(_strip_either1),mem::AbstractDict,x::Either) =
    if length(x.options) == 1
        deepmap_parser(_strip_either1,mem,first(x.options))
    else
        deepmap_either(_strip_either1,mem,x)
    end

deepmap_either(f,mem::AbstractDict,x::Either{<:Tuple},a...;kw...) =
    Either((deepmap_parser(f,mem,p,a...;kw...) for p in x.options)... )

function deepmap_either(f,mem::AbstractDict,x::Either{<:Vector},a...;kw...)
    mem[x] = r = Either{result_type(x)}(Any[])
    for p in x.options
        push!(r,deepmap_parser(f,mem,p,a...;kw...))
    end
    r
end

export substitute

struct Substitution<:CombinedParser{Nothing,Nothing}
    name::Symbol
end
"""
    substitute(name::Symbol)

Define a parser substitution.

    substitute(parser::CombinedParser)

Apply parser substitution, respecting scope in the defined tree:

- Parser variables are defined within scope of `Either`s, for all its `NamedParser` options.
- `Substitution` parsers are replaced with parser variables.
- [`strip_either1`](@ref) is used to simplify in a second phase.

!!! note
    Substitution implementation is experimental pending feedback.

    todo: scope NamedParser objects in WrappedParser, Sequence, etc.?

```jldoctest
julia> Either(:a => !Either(
                 :b => "X",
                 :d => substitute(:b),
                 substitute(:c)),
              :b => "b",
              :c => substitute(:b)
              ) |> substitute
|🗄 Either
├─ |🗄 Either |> ! |> with_name(:a)
│  ├─ X  |> with_name(:b)
│  ├─ X  |> with_name(:b) |> with_name(:d)
│  └─ b  |> with_name(:b) |> with_name(:c)
├─ b  |> with_name(:b)
└─ b  |> with_name(:b) |> with_name(:c)
::SubString{String}
```

# Example
With `substitute` you can write recursive parsers in a style inspired by (E)BNF.
[`CombinedParsers.BNF.ebnf`](@ref) uses `substitute`.

```jldoctest
julia> def = Either(:integer => !Either("0", Sequence(Optional("-"), substitute(:natural_number))),
                    :natural_number => !Sequence(substitute(:nonzero_digit), Repeat(substitute(:digit))),
                    :nonzero_digit => re"[1-9]",
                    :digit => Either("0", substitute(:nonzero_digit)))
|🗄 Either
├─ |🗄 Either |> ! |> with_name(:integer)
│  ├─ 0
│  └─ 🗄 Sequence
│     ├─ \\-? |
│     └─  natural_number call substitute!
├─ 🗄 Sequence |> ! |> with_name(:natural_number)
│  ├─  nonzero_digit call substitute!
│  └─ * digit call substitute! |> Repeat
├─ [1-9] ValueIn |> with_name(:nonzero_digit)
└─ |🗄 Either |> with_name(:digit)
   ├─ 0
   └─  nonzero_digit call substitute!
::Union{Nothing, Char, SubString{String}}

julia> substitute(def)
|🗄 Either
├─ |🗄 Either |> ! |> with_name(:integer)
│  ├─ 0
│  └─ 🗄 Sequence
│     ├─ \\-? |
│     └─ 🗄 Sequence |> ! |> with_name(:natural_number) # branches hidden
├─ 🗄 Sequence |> ! |> with_name(:natural_number)
│  ├─ [1-9] ValueIn |> with_name(:nonzero_digit)
│  └─ |🗄* Either |> with_name(:digit) |> Repeat
│     ├─ 0
│     └─ [1-9] ValueIn |> with_name(:nonzero_digit)
├─ [1-9] ValueIn |> with_name(:nonzero_digit)
└─ |🗄 Either |> with_name(:digit)
   ├─ 0
   └─ [1-9] ValueIn |> with_name(:nonzero_digit)
::Union{Char, SubString{String}}
```
"""
substitute(name::Symbol) =
    Substitution(name)
substitute(name::AbstractString) =
    Substitution(Symbol(name))
CombinedParsers.iterate_state(parser::Substitution, a...) = error(" call substitute")
function CombinedParsers.print_constructor(io::IO, x::Substitution)
    printstyled(io, x.name, color=:red)
    print(io, " call substitute!")
end
CombinedParsers._deepmap_parser(f,mem::AbstractDict,x::Substitution,a...;kw...) = x

substitute(x::CombinedParser) =
    strip_either1(deepmap_parser(_substitute, x, Dict{Symbol,CombinedParser}()))

_substitute(parser, assignments) = parser
function _substitute(parser::Substitution, assignments::AbstractDict)
    get(assignments,parser.name) do
        error("parser $(parser.name) is not defined")
    end
end

# set assignments
function deepmap_parser(::typeof(_substitute), mem::AbstractDict, x::Either, assignments)
    _assignments = copy(assignments)
    for o in x.options
        while o isa WrappedParser
            if o isa NamedParser
                _assignments[o.name] = Either(Any[])
            end
            o = o.parser
        end
    end
    for o in x.options
        while o isa WrappedParser
            if o isa NamedParser
                push!(_assignments[o.name], deepmap_parser(_substitute, mem, o, _assignments))
            end
            o = o.parser
        end
    end
    # deepmap_parser(_strip_either1, deepmap_either(_substitute, mem, x, _assignments))
    deepmap_either(_substitute, mem, x, _assignments)
end

# skip mem lookup (accept no assignments outside call/nesting parser stack!)
function deepmap_parser(::typeof(_substitute), mem::AbstractDict, x::Substitution, assignments)
    _substitute(x, assignments)
end


text:
# ---------------------------------------------------------------------------- #
#
#   writeTecplot.jl
#
#   Several functions to write Tecplot output
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    writeTecplot( flname::String, prob::Elas, mesh::Mesh2D,
                  uh::Array{Float64,3}, σh::Array{Float64,3};
                  σmises = Array{Float64,3}(0,0,0) )

Outputs tecplot data file to `flname` for 2D structural elasticity problems.
"""
function writeTecplot( flname::String, prob::Elas, mesh::Mesh2D, uh::Array{Float64,3},
  σh::Array{Float64,3}; σmises = Array{Float64,3}(0,0,0) )

  nelem   = size( mesh.nodes, 3 )
  nnodes  = size( mesh.nodes, 1 )
  ntriang = size( mesh.tloc,  1 ) # local number of triangles

  nnodesTot = nnodes * nelem

  nunkTot = mesh.dim + mesh.dim + mesh.dim^2 + ( size(σmises,1) > 0 )

  # Open file
  fid = open( flname, "w" )

  ### Write header
  @printf( fid, "TITLE = \"%s\"\n", prob.name )
  @printf( fid, "VARIABLES = \"X\", \"Y\"\n" )
  for jj in 1:mesh.dim
    @printf( fid, ", \"U%i\"",jj )
  end
  for ii in 1:mesh.dim, jj in 1:mesh.dim
    @printf( fid, ", \"<greek>s</greek>%i%i\"", ii, jj )
  end
  if size( σmises, 1 ) > 0
    @printf( fid, ", \"<greek>s</greek>mises\"" )
  end
  @printf( fid, "\n" )
  @printf( fid, "ZONE, DATAPACKING=POINT, NODES=%i, ELEMENTS=%i, ZONETYPE=FETRIANGLE\n",
    nnodesTot, nelem * ntriang )

  ### Write values
  for pp in 1:nelem

    wrtval = fill( 0.0, nnodes, nunkTot )

    if size( σmises, 1 ) == 0
        wrtval = [mesh.nodes[:,:,pp] uh[:,:,pp] σh[:,:,pp] ]
    else
        wrtval = [mesh.nodes[:,:,pp] uh[:,:,pp] σh[:,:,pp] σmises[:,:,pp] ]
    end

    for jj in 1:nnodes
      for ii in 1:nunkTot
        @printf(fid, "%16.15e\t", wrtval[jj,ii] )
      end
      @printf(fid,"\n")
    end
  end

  ### Write connectivity
  tnode = 0
  for pp in 1:nelem

    for jj in 1:ntriang
      for kk in 1:(mesh.dim+1)
        @printf(fid, "%i\t", tnode + mesh.tloc[jj,kk] )
      end
      @printf(fid,"\n")
    end

    tnode = tnode + nnodes

  end

  close(fid)

end

"""
    writeTecplot( flname::String, prob::Elas, mesh::Mesh3D,
                  uh::Array{Float64,3}, σh::Array{Float64,3};
                  σmises = Array{Float64,3}(0,0,0) )

Outputs tecplot data file to `flname` for 3D structural elasticity problems.
"""
function writeTecplot( flname::String, prob::Elas, mesh::Mesh3D, uh::Array{Float64,3},
  σh::Array{Float64,3}; σmises = Array{Float64,3}(0,0,0) )

  nelem   = size( mesh.nodes, 3 )
  nnodes  = size( mesh.nodes, 1 )
  ntets   = size( mesh.tloc,  1 ) # local number of tetrahedrons

  nnodesTot = nnodes * nelem

  nunkTot = mesh.dim + mesh.dim + mesh.dim^2 + ( size(σmises,1) > 0 )

  # Open file
  fid = open( flname, "w" )

  ### Write header
  @printf( fid, "TITLE = \"%s\"\n", prob.name )
  @printf( fid, "VARIABLES = \"X\", \"Y\", \"Z\"" )
  for jj in 1:mesh.dim
    @printf( fid, ", \"U%i\"",jj )
  end
  for ii in 1:mesh.dim, jj in 1:mesh.dim
    @printf( fid, ", \"<greek>s</greek>%i%i\"", ii, jj )
  end
  if size( σmises, 1 ) > 0
    @printf( fid, ", \"<greek>σ</greek>mises\"" )
  end
  @printf( fid, "\n" )
  @printf( fid, "ZONE, DATAPACKING=POINT, NODES=%i, ELEMENTS=%i, ZONETYPE=FETETRAHEDRON\n",
    nnodesTot, nelem * ntets )

  ### Write values
  for pp in 1:nelem

    wrtval = fill( 0.0, nnodes, nunkTot )

    if size( σmises, 1 ) == 0
        wrtval = [mesh.nodes[:,:,pp] uh[:,:,pp] σh[:,:,pp] ]
    else
        wrtval = [mesh.nodes[:,:,pp] uh[:,:,pp] σh[:,:,pp] σmises[:,:,pp] ]
    end

    for jj in 1:nnodes
      for ii in 1:nunkTot
        @printf(fid, "%16.15e\t", wrtval[jj,ii] )
      end
      @printf(fid,"\n")
    end
  end

  ### Write connectivity
  tnode = 0::Int64
  for pp in 1:nelem

    for jj in 1:ntets
      for kk in 1:(mesh.dim+1)
        @printf(fid, "%i\t", tnode + mesh.tloc[jj,kk] )
      end
      @printf(fid,"\n")
    end

    tnode = tnode + nnodes

  end

  close(fid)

end

"""
    writeTecplotP1( flname::String, prob::Elas, mesh::Mesh3D,
                  uh::Array{Float64,3}, σh::Array{Float64,3};
                  σmises = Array{Float64,3}(0,0,0) )

Outputs tecplot data file to `flname` for 3D structural elasticity problems.
"""
function writeTecplotP1( flname::String, prob::Elas, mesh::Mesh3D, uh::Array{Float64,3},
                         σh::Array{Float64,3}, σmises::Array{Float64,3} )

  nelem   = size( mesh.nodes, 3 )
  nnodes  = size( mesh.nodes, 1 )

  nnodesTot = mesh.n

  nunkTot = mesh.dim + mesh.dim + mesh.dim^2 + 1

  # Open file
  fid = open( flname, "w" )

  ### Write header
  @printf( fid, "TITLE = \"%s\"\n", prob.name )
  @printf( fid, "VARIABLES = \"X\", \"Y\", \"Z\"" )
  for jj in 1:mesh.dim
    @printf( fid, ", \"U%i\"",jj )
  end
  for ii in 1:mesh.dim, jj in 1:mesh.dim
    @printf( fid, ", \"<greek>s</greek>%i%i\"", ii, jj )
  end
  @printf( fid, ", \"<greek>s</greek>mises\"" )
  @printf( fid, "\n" )
  @printf( fid, "ZONE, DATAPACKING=POINT, NODES=%i, ELEMENTS=%i, ZONETYPE=FETETRAHEDRON\n",
    nnodesTot, nelem )

  ### Interpolate values on nodes
  nelpoin = fill( 0,  mesh.n )
  values  = fill( 0., mesh.n, nunkTot )
  values[:, 1:mesh.dim] = mesh.p

  for el in 1:nelem, jj in 1:(mesh.dim+1)

      p0 = mesh.t[el,jj]
      values[ p0, mesh.dim+(1:mesh.dim)   ] += uh[jj,:,el]
      values[ p0, 2*mesh.dim+(1:mesh.dim^2) ] += σh[jj,:,el]
      values[ p0, end ] += σmises[jj,1,el]
      nelpoin[p0] += 1

  end

  for pp in 1:mesh.n

    values[pp,(mesh.dim+1):end] ./= nelpoin[pp]

    for ii in 1:nunkTot
      @printf(fid, "%16.15e\t", values[pp,ii] )
    end
    @printf(fid,"\n")

  end

  ### Write connectivity
  for ee in 1:nelem

    for kk in 1:(mesh.dim+1)
      @printf(fid, "%i\t", mesh.t[ee,kk] )
    end
    @printf(fid,"\n")

  end

  close(fid)

end


text:
## activate project environment
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

## load packages
using DelimitedFiles, LinearAlgebra, Statistics
using BenchmarkTools
using Trixi

## set up benchmark code
"""
    initial_condition_isentropic_vortex(x, t, equations)

The classical isentropic vortex test case of
- Chi-Wang Shu (1997)
  Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory
  Schemes for Hyperbolic Conservation Laws.
  [NASA/CR-97-206253](https://ntrs.nasa.gov/citations/19980007543)
"""
function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)
  ϱ0 = 1.0               # background density
  v0 = SVector(1.0, 1.0) # background velocity
  p0 = 10.0              # background pressure
  ε = 20.0               # vortex strength
  L = 10.0               # size of the domain per coordinate direction

  T0 = p0 / ϱ0           # background temperature
  γ = equations.gamma    # ideal gas constant

  x0 = v0 * t            # current center of the vortex
  dx = vortex_center.(x - x0, L)
  r2 = sum(abs2, dx)

  # perturbed primitive variables
  T = T0 - (γ - 1) * ε^2 / (8 * γ * π^2) * exp(1 - r2)
  v = v0 + ε / (2 * π) * exp(0.5 * (1 - r2)) * SVector(-dx[2], dx[1])
  ϱ = ϱ0 * (T / T0)^(1 / (γ - 1))
  p = ϱ * T

  return prim2cons(SVector(ϱ, v..., p), equations)
end

function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations3D)
  ϱ0 = 1.0               # background density
  v0 = SVector(1.0, 1.0, 0.0) # background velocity
  p0 = 10.0              # background pressure
  ε = 20.0               # vortex strength
  L = 10.0               # size of the domain per coordinate direction

  T0 = p0 / ϱ0           # background temperature
  γ = equations.gamma    # ideal gas constant

  x0 = v0 * t            # current center of the vortex
  dx = vortex_center.(x - x0, L)
  r2 = sum(abs2, dx)

  # perturbed primitive variables
  T = T0 - (γ - 1) * ε^2 / (8 * γ * π^2) * exp(1 - r2)
  v = v0 + ε / (2 * π) * exp(0.5 * (1 - r2)) * SVector(-dx[2], dx[1], 0.0)
  ϱ = ϱ0 * (T / T0)^(1 / (γ - 1))
  p = ϱ * T

  return prim2cons(SVector(ϱ, v..., p), equations)
end

vortex_center(x, L) = mod(x + L/2, L) - L/2

# 2D version
function overintegration_kernel!(du_overinte, du_fluxdiff, u_overinte, u_fluxdiff,
                                 lo2hi, tmp_lo2hi, hi2lo, tmp_hi2lo,
                                 element, mesh,
                                 nonconservative_terms, equations,
                                 dg, cache)

  Trixi.multiply_dimensionwise!(view(u_overinte, :, :, :, element),
                                lo2hi,
                                view(u_fluxdiff, :, :, :, element),
                                tmp_lo2hi)

  Trixi.weak_form_kernel!(du_overinte, u_overinte, element, mesh,
                          nonconservative_terms, equations, dg, cache)

  Trixi.multiply_dimensionwise!(view(du_fluxdiff, :, :, :, element),
                                hi2lo,
                                view(du_overinte, :, :, :, element),
                                tmp_hi2lo)

  return nothing
end

function run_benchmarks_2d(polydeg_fluxdiff, polydeg_overinte, numerical_flux)
  println("Benchmarking ", numerical_flux, " in 2D with polydeg=", polydeg_fluxdiff)

  # setup basic code
  equations = CompressibleEulerEquations2D(1.4)

  coordinates_min = (-5.0, -5.0)
  coordinates_max = ( 5.0,  5.0)
  mesh = TreeMesh(coordinates_min, coordinates_max,
                  initial_refinement_level=0, n_cells_max=100_000, periodicity=true)

  solver_fluxdiff = DGSEM(polydeg=polydeg_fluxdiff, surface_flux=numerical_flux,
                          volume_integral=VolumeIntegralFluxDifferencing(numerical_flux))

  solver_overinte = DGSEM(polydeg=polydeg_overinte, surface_flux=numerical_flux,
                          volume_integral=VolumeIntegralWeakForm())

  semi_fluxdiff = SemidiscretizationHyperbolic(mesh, equations,
                                              initial_condition_isentropic_vortex, solver_fluxdiff)

  semi_overinte = SemidiscretizationHyperbolic(mesh, equations,
                                              initial_condition_isentropic_vortex, solver_overinte)

  u_ode_fluxdiff  = compute_coefficients(0.0, semi_fluxdiff)
  du_ode_fluxdiff = similar(u_ode_fluxdiff)
  u_ode_overinte  = compute_coefficients(0.0, semi_overinte)
  du_ode_overinte = similar(u_ode_overinte)
  GC.@preserve u_ode_fluxdiff du_ode_fluxdiff u_ode_overinte du_ode_overinte begin
    u_fluxdiff  = Trixi.wrap_array(u_ode_fluxdiff,  semi_fluxdiff)
    du_fluxdiff = Trixi.wrap_array(du_ode_fluxdiff, semi_fluxdiff)
    u_overinte  = Trixi.wrap_array(u_ode_overinte,  semi_overinte)
    du_overinte = Trixi.wrap_array(du_ode_overinte, semi_overinte)
    element = 1

    # benchmark flux differencing volume terms
    result_fluxdiff = @benchmark Trixi.split_form_kernel!(
      $du_fluxdiff, $u_fluxdiff, $element, $mesh,
      $(Trixi.have_nonconservative_terms(equations)), $equations,
      $numerical_flux, $solver_fluxdiff, $(semi_fluxdiff.cache)) setup=(fill!($du_fluxdiff, zero(eltype($du_fluxdiff))))
    display(result_fluxdiff)

    # benchmark overintegration volume terms
    nodes_lo = solver_fluxdiff.basis.nodes
    nodes_hi = solver_overinte.basis.nodes
    lo2hi = Trixi.polynomial_interpolation_matrix(nodes_lo, nodes_hi)

    V_lo, _ = Trixi.vandermonde_legendre(nodes_lo)
    V_lo = [V_lo zeros(length(nodes_lo), length(nodes_hi) - length(nodes_lo))]
    V_hi, _ = Trixi.vandermonde_legendre(nodes_hi)
    d = zeros(length(nodes_hi)); d[1:length(nodes_lo)] .= 1
    proj = Diagonal(d)
    hi2lo = V_lo * proj / V_hi

    tmp_lo2hi = zeros(size(u_fluxdiff, 1), size(lo2hi)...)
    tmp_hi2lo = zeros(size(u_fluxdiff, 1), size(hi2lo)...)

    result_overinte = @benchmark overintegration_kernel!(
      $du_overinte, $du_fluxdiff, $u_overinte, $u_fluxdiff,
      $lo2hi, $tmp_lo2hi, $hi2lo, $tmp_hi2lo,
      $element, $mesh,
      $(Trixi.have_nonconservative_terms(equations)), $equations,
      $solver_overinte, $(semi_overinte.cache)) setup=(fill!($du_overinte, zero(eltype($du_overinte))))
    display(result_overinte)
  end
  println()

  return result_fluxdiff, result_overinte
end

# 3D version
function overintegration_kernel!(du_overinte, du_fluxdiff, u_overinte, u_fluxdiff,
                                 lo2hi, tmp1_lo2hi, tmp2_lo2hi,
                                 hi2lo, tmp1_hi2lo, tmp2_hi2lo,
                                 element, mesh,
                                 nonconservative_terms, equations,
                                 dg, cache)

  Trixi.multiply_dimensionwise!(view(u_overinte, :, :, :, :, element),
                                lo2hi,
                                view(u_fluxdiff, :, :, :, :, element),
                                tmp1_lo2hi, tmp2_lo2hi)

  Trixi.weak_form_kernel!(du_overinte, u_overinte, element, mesh,
                          nonconservative_terms, equations, dg, cache)

  Trixi.multiply_dimensionwise!(view(du_fluxdiff, :, :, :, :, element),
                                hi2lo,
                                view(du_overinte, :, :, :, :, element),
                                tmp1_hi2lo, tmp2_hi2lo)

  return nothing
end

function run_benchmarks_3d(polydeg_fluxdiff, polydeg_overinte, numerical_flux)
  println("Benchmarking ", numerical_flux, " in 3D with polydeg=", polydeg_fluxdiff)

  # setup basic code
  equations = CompressibleEulerEquations3D(1.4)

  coordinates_min = (-5.0, -5.0, -5.0)
  coordinates_max = ( 5.0,  5.0,  5.0)
  mesh = TreeMesh(coordinates_min, coordinates_max,
                  initial_refinement_level=0, n_cells_max=100_000, periodicity=true)

  solver_fluxdiff = DGSEM(polydeg=polydeg_fluxdiff, surface_flux=numerical_flux,
                          volume_integral=VolumeIntegralFluxDifferencing(numerical_flux))

  solver_overinte = DGSEM(polydeg=polydeg_overinte, surface_flux=numerical_flux,
                          volume_integral=VolumeIntegralWeakForm())

  semi_fluxdiff = SemidiscretizationHyperbolic(mesh, equations,
                                              initial_condition_isentropic_vortex, solver_fluxdiff)

  semi_overinte = SemidiscretizationHyperbolic(mesh, equations,
                                              initial_condition_isentropic_vortex, solver_overinte)

  u_ode_fluxdiff  = compute_coefficients(0.0, semi_fluxdiff)
  du_ode_fluxdiff = similar(u_ode_fluxdiff)
  u_ode_overinte  = compute_coefficients(0.0, semi_overinte)
  du_ode_overinte = similar(u_ode_overinte)
  GC.@preserve u_ode_fluxdiff du_ode_fluxdiff u_ode_overinte du_ode_overinte begin
    u_fluxdiff  = Trixi.wrap_array(u_ode_fluxdiff,  semi_fluxdiff)
    du_fluxdiff = Trixi.wrap_array(du_ode_fluxdiff, semi_fluxdiff)
    u_overinte  = Trixi.wrap_array(u_ode_overinte,  semi_overinte)
    du_overinte = Trixi.wrap_array(du_ode_overinte, semi_overinte)
    element = 1

    # benchmark flux differencing volume terms
    result_fluxdiff = @benchmark Trixi.split_form_kernel!(
      $du_fluxdiff, $u_fluxdiff, $element, $mesh,
      $(Trixi.have_nonconservative_terms(equations)), $equations,
      $numerical_flux, $solver_fluxdiff, $(semi_fluxdiff.cache)) setup=(fill!($du_fluxdiff, zero(eltype($du_fluxdiff))))
    display(result_fluxdiff)

    # benchmark overintegration volume terms
    nodes_lo = solver_fluxdiff.basis.nodes
    nodes_hi = solver_overinte.basis.nodes
    lo2hi = Trixi.polynomial_interpolation_matrix(nodes_lo, nodes_hi)

    V_lo, _ = Trixi.vandermonde_legendre(nodes_lo)
    V_lo = [V_lo zeros(length(nodes_lo), length(nodes_hi) - length(nodes_lo))]
    V_hi, _ = Trixi.vandermonde_legendre(nodes_hi)
    d = zeros(length(nodes_hi)); d[1:length(nodes_lo)] .= 1
    proj = Diagonal(d)
    hi2lo = V_lo * proj / V_hi

    tmp1_lo2hi = zeros(size(u_fluxdiff, 1), size(lo2hi, 1), size(lo2hi, 2), size(lo2hi, 2))
    tmp2_lo2hi = zeros(size(u_fluxdiff, 1), size(lo2hi, 1), size(lo2hi, 1), size(lo2hi, 2))
    tmp1_hi2lo = zeros(size(u_fluxdiff, 1), size(hi2lo, 1), size(hi2lo, 2), size(hi2lo, 2))
    tmp2_hi2lo = zeros(size(u_fluxdiff, 1), size(hi2lo, 1), size(hi2lo, 1), size(hi2lo, 2))

    result_overinte = @benchmark overintegration_kernel!(
      $du_overinte, $du_fluxdiff, $u_overinte, $u_fluxdiff,
      $lo2hi, $tmp1_lo2hi, $tmp2_lo2hi, $hi2lo, $tmp1_hi2lo, $tmp2_hi2lo,
      $element, $mesh,
      $(Trixi.have_nonconservative_terms(equations)), $equations,
      $solver_overinte, $(semi_overinte.cache)) setup=(fill!($du_overinte, zero(eltype($du_overinte))))
    display(result_overinte)
  end
  println()

  return result_fluxdiff, result_overinte
end

# general code
Base.ndims(::typeof(run_benchmarks_2d)) = 2
Base.ndims(::typeof(run_benchmarks_3d)) = 3

function run_benchmarks(run_benchmarks_inner, polydegs, numerical_flux)
  results = zeros(length(polydegs), 1)
  results .= polydegs
  header = "# Polydeg"

  # flux differencing
  println("#"^80)
  println("flux differencing\n")
  means = zeros(length(polydegs))
  stds  = zero(means)
  for (i, polydeg) in enumerate(polydegs)
    result, _ = run_benchmarks_inner(polydeg, polydeg, numerical_flux)
    means[i] = time(mean(result)) * 1.0e-9
    stds[i]  = time(std(result))  * 1.0e-9
  end
  results = hcat(results, means, stds)
  header = header * "\tFluxDifferencing-mean\tstd"

  # overintegration, polydeg + 1
  println("#"^80)
  println("overintegration, polydeg + 1\n")
  means = zeros(length(polydegs))
  stds  = zero(means)
  for (i, polydeg) in enumerate(polydegs)
    _, result = run_benchmarks_inner(polydeg, polydeg + 1, flux_central)
    means[i] = time(mean(result)) * 1.0e-9
    stds[i]  = time(std(result))  * 1.0e-9
  end
  results = hcat(results, means, stds)
  header = header * "\tpolydeg+1-mean\tstd"

  # overintegration, polydeg + 2
  println("#"^80)
  println("overintegration, polydeg + 2\n")
  means = zeros(length(polydegs))
  stds  = zero(means)
  for (i, polydeg) in enumerate(polydegs)
    _, result = run_benchmarks_inner(polydeg, polydeg + 2, flux_central)
    means[i] = time(mean(result)) * 1.0e-9
    stds[i]  = time(std(result))  * 1.0e-9
  end
  results = hcat(results, means, stds)
  header = header * "\tpolydeg+2-mean\tstd"

  # overintegration, 3 * polydeg ÷ 2
  println("#"^80)
  println("overintegration, 3 * polydeg ÷ 2\n")
  means = zeros(length(polydegs))
  stds  = zero(means)
  for (i, polydeg) in enumerate(polydegs)
    _, result = run_benchmarks_inner(polydeg, (3 * polydeg) ÷ 2, flux_central)
    means[i] = time(mean(result)) * 1.0e-9
    stds[i]  = time(std(result))  * 1.0e-9
  end
  results = hcat(results, means, stds)
  header = header * "\t3*polydeg÷2-mean\tstd"

  # overintegration, 2 * polydeg
  println("#"^80)
  println("overintegration, 2 * polydeg\n")
  means = zeros(length(polydegs))
  stds  = zero(means)
  for (i, polydeg) in enumerate(polydegs)
    _, result = run_benchmarks_inner(polydeg, 2 * polydeg, flux_central)
    means[i] = time(mean(result)) * 1.0e-9
    stds[i]  = time(std(result))  * 1.0e-9
  end
  results = hcat(results, means, stds)
  header = header * "\t2*polydeg-mean\tstd"

  # write results to a file
  open(joinpath(@__DIR__, "overintegration_$(ndims(run_benchmarks_inner))D_$(numerical_flux).dat"), "w") do io
    println(io, header)
    writedlm(io, results)
  end

  return results
end

function run_benchmarks_2d(polydegs, numerical_flux)
  run_benchmarks(run_benchmarks_2d, polydegs, numerical_flux)
end

function run_benchmarks_3d(polydegs, numerical_flux)
  run_benchmarks(run_benchmarks_3d, polydegs, numerical_flux)
end

## run benchmarks
polydegs = 3:15
run_benchmarks_2d(polydegs, flux_shima_etal)
run_benchmarks_2d(polydegs, flux_ranocha)

run_benchmarks_3d(polydegs, flux_shima_etal)
run_benchmarks_3d(polydegs, flux_ranocha)


text:
module FileIO

export DataFormat,
       File,
       Formatted,
       Stream,
       @format_str,
       add_format,
       del_format,
       add_loader,
       add_saver,
       filename,
       file_extension,
       load,
       loadstreaming,
       magic,
       query,
       save,
       savestreaming,
       skipmagic,
       stream,
       unknown,
       metadata

import Base.showerror
using Base: RefValue
using Pkg

include("types.jl")
include("registry_setup.jl")
include("query.jl")
include("error_handling.jl")
include("loadsave.jl")
include("mimesave.jl")
include("registry.jl")

"""
`FileIO` API (brief summary, see individual functions for more detail):

- `format"PNG"`: specifies a particular defined format
- `File{fmt}` and `Stream{fmt}`: types of objects that declare that a resource has a particular format `fmt`

- `load([filename|stream])`: read data in formatted file, inferring the format
- `load(File(format"PNG",filename))`: specify the format manually
- `loadstreaming([filename|stream])`: similar to `load`, except that it returns an object that can be read from
- `save(filename, data...)` for similar operations involving saving data
- `savestreaming([filename|stream])`: similar to `save`, except that it returns an object that can be written to

- `io = open(f::File, args...)` opens a file
- `io = stream(s::Stream)` returns the IOStream from the query object `s`

- `query([filename|stream])`: attempt to infer the format of `filename`
- `unknown(q)` returns true if a query can't be resolved
- `skipmagic(io, fmt)` sets the position of `io` to just after the magic bytes
- `magic(fmt)` returns the magic bytes for format `fmt`
- `info(fmt)` returns `(magic, extensions)` for format `fmt`

- `add_format(fmt, magic, extension)`: register a new format
- `add_loader(fmt, :Package)`: indicate that `Package` supports loading files of type `fmt`
- `add_saver(fmt, :Package)`: indicate that `Package` supports saving files of type `fmt`
"""
FileIO

end


text:
# A supported provider whose hosted repositories can be registered.
Base.@kwdef struct Provider{F <: GitForge.Forge}
    name::String
    client::F
    client_id::String
    client_secret::String
    auth_url::String
    token_url::String
    scope::String
    include_state::Bool = true
    token_type::Type = typeof(client.token)
end

function init_providers()
    if  haskey(CONFIG, "github")
        github = CONFIG["github"]
        PROVIDERS["github"] = Provider(;
            name="GitHub",
            client=GitHubAPI(;
                url=get(github, "api_url", GitHub.DEFAULT_URL),
                token=Token(github["token"]),
                has_rate_limits=!get(github, "disable_rate_limits", false),
            ),
            client_id=github["client_id"],
            client_secret=github["client_secret"],
            auth_url=get(github, "auth_url", "https://github.com/login/oauth/authorize"),
            token_url=get(github, "token_url", "https://github.com/login/oauth/access_token"),
            scope="public_repo",
        )
    end

    if haskey(CONFIG, "gitlab")
        gitlab = CONFIG["gitlab"]
        PROVIDERS["gitlab"] = Provider(;
            name="GitLab",
            client=GitLabAPI(;
                url=get(gitlab, "api_url", GitLab.DEFAULT_URL),
                token=PersonalAccessToken(gitlab["token"]),
                has_rate_limits=!get(gitlab, "disable_rate_limits", false),
            ),
            client_id=gitlab["client_id"],
            client_secret=gitlab["client_secret"],
            auth_url=get(gitlab, "auth_url", "https://gitlab.com/oauth/authorize"),
            token_url=get(gitlab, "token_url", "https://gitlab.com/oauth/token"),
            scope="read_user api",
            include_state=false,
            token_type=OAuth2Token,
        )
    end

    haskey(CONFIG, "extra_providers") && include(CONFIG["extra_providers"])
end


text:
"""
    ME( y, yhat )

Calculates Mean Error from vectors `Y` and `YHat`.
"""
function ME( y, yhat )
    return ( 1.0 / size(Y)[1] ) * sum( ( y - yhat ) )
end

"""
    MAE( y, yhat )

Calculates Mean Average Error from vectors `Y` and `YHat`
"""
function MAE( y, yhat )
    return ( 1.0 / size(Y)[1] ) * sum( abs.( ( y - yhat ) ) )
end

"""
    MAPE( y, yhat )

Calculates Mean Average Percent Error from vectors `Y` and `YHat`
"""
function MAPE( y, yhat )
    return ( 100.0 / size(Y)[1] ) * sum( abs.( ( y - yhat ) / y ) )
end

"""
    SSE( y, yhat )

Calculates Sum of Squared Errors from vectors `Y` and `YHat`
"""
function SSE( y, yhat )
    return sum( (yhat .- y) .^ 2  )
end

"""
    MSE( y, yhat )

Calculates Mean Squared Error from vectors `Y` and `YHat`
"""
function MSE( y, yhat )
    return SSE(y, yhat) / size(y)[1]
end

"""
    RMSE( y, yhat )

Calculates Root Mean Squared Error from vectors `Y` and `YHat`
"""
function RMSE( y, yhat )
    return sqrt( SSE(y, yhat) / size(y)[1] )
end

"""
    PercentRMSE( y, yhat; aspercent = true )

Calculates Percent Root Mean Squared Error from vectors `Y` and `YHat`
"""
function PercentRMSE( y, yhat; aspercent = true )
    mini = reduce(min, yhat)
    maxi = reduce(max, yhat)
    scalar = (aspercent) ? 100.0 : 1
    return scalar * RMSE(y,yhat) / (maxi - mini)
end

"""
    SSTotal( y, yhat )

Calculates Total Sum of Squared Deviations from vectors `Y` and `YHat`
"""
function SSTotal( y )
    return sum( ( y .- StatsBase.mean( y ) ) .^ 2 )
end

"""
    SSReg( y, yhat )

Calculates Sum of Squared Deviations due to Regression from vectors `Y` and `YHat`
"""
function SSReg( y, yhat )
    return sum( ( yhat .- StatsBase.mean( y ) ) .^ 2 )
end

"""
    SSRes( y, yhat )

Calculates Sum of Squared Residuals from vectors `Y` and `YHat`
"""
function SSRes( y, yhat )
    return sum( ( y - yhat ) .^ 2 )
end

"""
    RSquare( y, yhat )

Calculates R^2 from `Y` and `YHat`
"""
function RSquare( y, yhat )
    return 1.0 - ( SSRes(y, yhat) / SSTotal(y) )
end

"""
    FNorm( X )

Calculates the Froebinius norm of matrix X.
"""
FNorm( X ) = sqrt.( sum( X .^ 2 ) )

"""
    PearsonCorrelationCoefficient( y, yhat )

Calculates The Pearson Correlation Coefficient from vectors `Y` and `YHat`
"""
function PearsonCorrelationCoefficient(y, yhat)
    return StatsBase.cov( y, yhat ) / ( StatsBase.std( y ) * StatsBase.std( yhat )  )
end


text:
using SentinelArrays, Test

@testset "SentinelArrays" begin

x = SentinelVector{Float64}(undef, 10)
fill!(x, missing)
@test length(x) == 10
@test all(isequal.(x.data, x.sentinel))

@test x[1] === missing
# force recoding
sent = x.sentinel
x[1] = x.sentinel
@test x[1] === sent

@test size(x) == (10,)
x[1] = 3.14
@test x[1] === 3.14

resize!(x, length(x) + 1)
@test x[end] === missing

x = SentinelVector{Int64}(undef, 1)
x[1] = missing
@test x[1] === missing

x = SentinelVector{Union{Bool, Missing}}(undef, 1, missing, missing)
@test x[1] === missing

x = SentinelVector{String}(undef, 10)
@test x[1] === missing
x[1] = "hey"
@test x[1] == "hey"
x[1] = missing
@test x[1] === missing

@test all(x .=== copy(x))
@test length(empty!(x)) == 0
@test length(resize!(x, 10)) == 10
@test x[1] === missing

push!(x, missing)
@test x[end] === missing
push!(x, "ho")
@test x[end] == "ho"

deleteat!(x, length(x))
@test x[end] === missing
@test length(x) == 11

empty!(x)
append!(x, ["hey", "ho", missing])
@test length(x) == 3
@test isequal(x, ["hey", "ho", missing])

pushfirst!(x, missing)
@test x[1] === missing
prepend!(x, [missing, "first"])
@test x[1] === missing
@test x[2] == "first"

@test popfirst!(x) === missing
@test x[1] == "first"

insert!(x, 1, missing)
@test x[1] === missing
insert!(x, length(x) + 1, "pirate")
@test x[end] == "pirate"

@test splice!(x, length(x)) == "pirate"

t = [SentinelVector{Int64}(undef, 10), SentinelVector{Int64}(undef, 10), SentinelVector{Int64}(undef, 5)]
sent = t[1].sentinel
@test all(x->x.sentinel == sent, t)
SentinelArrays.newsentinel!(t...; force=false)
# make sure nothing got recoded
@test all(x->x.sentinel == sent, t)
SentinelArrays.newsentinel!(t...; force=true)
# make sure all got recoded
@test all(x->x.sentinel != sent, t)

# force recode of just one vector
sent = t[1].sentinel
t[2][1] = t[2].sentinel
@test !all(x->x.sentinel == sent, t)
SentinelArrays.newsentinel!(t...)
# make sure all got recoded to same
@test all(x->x.sentinel == t[1].sentinel, t)

A = SentinelArray(Int64[i for i = 1:10])
B = SentinelVector{Int64}(undef, 10)

C = vcat(A, B)
@test C[1:10] == collect(1:10)
@test all(C[11:20] .=== missing)

append!(A, B)
@test A[1:10] == collect(1:10)
@test all(A[11:20] .=== missing)

A = SentinelArray(Int64[i for i = 1:10])
B = SentinelVector{Int64}(undef, 10)
# force B to recode
B[1] = B.sentinel
B[1] = missing

C = vcat(A, B)
@test C[1:10] == collect(1:10)
@test all(C[11:20] .=== missing)

append!(A, B)
@test A[1:10] == collect(1:10)
@test all(A[11:20] .=== missing)

# make sure we bail when can't find a new automatic sentinel
A = SentinelArray([i for i = 0x00:0xff])
@test_throws SentinelCollisionError setindex!(A, A.sentinel, 1)

@test_throws ErrorException SentinelVector{Bool}(undef, 1)

t = SentinelVector{Tuple{Int32, Int32}}(undef, 1)
@test t.data[1] === t.sentinel

t = SentinelMatrix{Float64}(undef, (10, 10))
@test size(t) == (10, 10)

end # @testset


text:
struct WebSocketCommand
    syncId::String
    command::Symbol
    subCommand::Optional{Symbol}
    content::Any
end

WebSocketCommand(syncId, cmd::GeneralCommand) = WebSocketCommand(syncId, cmd.command, cmd.subCommand, cmd.content)

StructTypes.StructType(::Type{WebSocketCommand}) = StructTypes.Struct()
StructTypes.omitempties(::Type{WebSocketCommand}) = (:subCommand, :content)

mutable struct WebSocketAdapter <: ProtocolAdapter
    reservedSyncId::String
    output_channel::OutputChannel
    output_dict::Dict{String,OutputChannel}
    ready::Bool
    ready_cond::Condition
    socket::HTTP.WebSockets.WebSocket
    function WebSocketAdapter(reservedSyncId = "-1")
        new(reservedSyncId, OutputChannel(DEFAULT_BUFFER_SIZE), Dict{String,OutputChannel}(), false, Condition())
    end
end

function loop(adp::WebSocketAdapter, server, qq, verifyKey)
    url = "ws://$server/message?verifyKey=$verifyKey&qq=$qq"
    HTTP.WebSockets.open(url) do ws
        adp.socket = ws
        hello = JSON3.read(readavailable(ws))
        @assert hasproperty(hello, :data)
        set_ready(adp)
        put!(adp.output_channel, hello.data)
        while !eof(ws)
            with_catch_backtrace(adp.output_channel) do
                bytes = readavailable(ws)
                isempty(bytes) && eof(ws) && return
                @debug "Received data: $(String(copy(bytes)))"
                data = JSON3.read(bytes)::JSON3.Object
                dispatch(adp, data[:syncId], data[:data])
            end
        end
        # close all pending tasks
        close(adp.output_channel)
        foreach(close, values(adp.output_dict))
        @info "WebSocket adapter receiver quitted"
        @info "WebSocket adapter quitted"
    end
end

function with_entry(f, dict::AbstractDict, key, value)
    dict[key] = value
    try
        f()
    finally
        delete!(dict, key)
    end
end

function send(adp::WebSocketAdapter, cmd::GeneralCommand; session_key_position)
    wait_ready(adp)
    wscmd = WebSocketCommand(get_sync_id(), cmd)
    ret_ch = OutputChannel(0) do ch
        with_entry(adp.output_dict, wscmd.syncId, ch) do
            @debug "Command: $cmd"
            str = JSON3.write(wscmd)
            @debug "Sending command: $str"
            write(adp.socket, str)
            wait(ch)
        end
    end
    take!(ret_ch)
end

function dispatch(adp::WebSocketAdapter, syncId, data)
    ch = get(adp.output_dict, syncId, nothing)
    if isnothing(ch)
        msg = try_convert(EventOrMessage, data)
        @debug msg
        put!(adp.output_channel, msg)
    else
        put!(ch, data)
    end
end

function Base.close(adp::WebSocketAdapter)
    close(adp.socket)
end

function wait_ready(adp::WebSocketAdapter)
    while !adp.ready
        wait(adp.ready_cond)
    end
end

function set_ready(adp::WebSocketAdapter)
    adp.ready = true
    notify(adp.ready_cond)
end


text:
struct ImplicitDynamics{T,R,RZ,Rθ,M<:RoboDojo.Model{T},P<:RoboDojo.Policy{T},D<:RoboDojo.Disturbances{T},I} <: Model{T}
    n::Int
    m::Int
    d::Int
	eval_sim::Simulator{T,R,RZ,Rθ,M,P,D}
	grad_sim::Simulator{T,R,RZ,Rθ,M,P,D}
	q1::Vector{T}
	q2::Vector{T}
	v1::Vector{T}
	idx_q1::Vector{Int}
	idx_q2::Vector{Int}
	idx_u1::Vector{Int}
	info::I
end

function get_simulator(model, h, r_func, rz_func, rθ_func;
	T=1, r_tol=1.0e-8, κ_eval_tol=1.0e-4, nc=model.nc, nb=model.nc, diff_sol=true)

	sim = Simulator(model, T;
        h=h,
        residual=r_func,
        jacobian_z=rz_func,
        jacobian_θ=rθ_func,
        diff_sol=diff_sol,
        solver_opts=InteriorPointOptions(
            undercut=Inf,
            γ_reg=0.1,
            r_tol=r_tol,
            κ_tol=κ_eval_tol,
            max_ls=25,
            ϵ_min=0.25,
            diff_sol=diff_sol,
            verbose=false))

    # set trajectory sizes
	sim.traj.γ .= [zeros(nc) for t = 1:T]
	sim.traj.b .= [zeros(nb) for t = 1:T]

    sim.grad.∂γ1∂q1 .= [zeros(nc, model.nq) for t = 1:T]
	sim.grad.∂γ1∂q2 .= [zeros(nc, model.nq) for t = 1:T]
	sim.grad.∂γ1∂v1 .= [zeros(nc, model.nq) for t = 1:T]
	sim.grad.∂γ1∂u1 .= [zeros(nc, model.nu) for t = 1:T]
	sim.grad.∂b1∂q1 .= [zeros(nb, model.nq) for t = 1:T]
	sim.grad.∂b1∂q2 .= [zeros(nb, model.nq) for t = 1:T]
	sim.grad.∂b1∂v1 .= [zeros(nb, model.nq) for t = 1:T]
	sim.grad.∂b1∂u1 .= [zeros(nb, model.nu) for t = 1:T]

    return sim
end

function ImplicitDynamics(model, h, r_func, rz_func, rθ_func;
	T=1, r_tol=1.0e-8, κ_eval_tol=1.0e-6, κ_grad_tol=1.0e-6,
	no_impact=false, no_friction=false,
	n=(2 * model.nq), m=model.nu, d=model.nw, nc=model.nc, nb=model.nc,
	info=nothing)

	# set trajectory sizes
	no_impact && (nc = 0)
	no_friction && (nb = 0)

	eval_sim = get_simulator(model, h, r_func, rz_func, rθ_func;
			T=T, r_tol=r_tol, κ_eval_tol=κ_eval_tol, nc=nc, nb=nb, diff_sol=false)

	grad_sim = get_simulator(model, h, r_func, rz_func, rθ_func;
			T=T, r_tol=r_tol, κ_eval_tol=κ_grad_tol, nc=nc, nb=nb, diff_sol=true)

	q1 = zeros(model.nq)
	q2 = zeros(model.nq)
	v1 = zeros(model.nq)

	idx_q1 = collect(1:model.nq)
	idx_q2 = collect(model.nq .+ (1:model.nq))
	idx_u1 = collect(1:model.nu)

	ImplicitDynamics(n, m, d,
		eval_sim, grad_sim,
		q1, q2, v1,
		idx_q1, idx_q2, idx_u1, info)
end

function f(d, model::ImplicitDynamics, x, u, w)
	q1 = @views x[model.idx_q1]
	q2 = @views x[model.idx_q2]
	model.v1 .= q2
	model.v1 .-= q1
	model.v1 ./= model.eval_sim.h

	q3 = RoboDojo.step!(model.eval_sim, q2, model.v1, u, 1)

	d[model.idx_q1] .= q2
	d[model.idx_q2] .= q3

	return d
end

function fx(dx, model::ImplicitDynamics, x, u, w)
	q1 = @views x[model.idx_q1]
	q2 = @views x[model.idx_q2]
	model.v1 .= q2
	model.v1 .-= q1
	model.v1 ./= model.grad_sim.h

	RoboDojo.step!(model.grad_sim, q2, model.v1, u, 1)

	nq = model.grad_sim.model.nq
	for i = 1:nq
		dx[model.idx_q1[i], model.idx_q2[i]] = 1.0
	end

	dx[model.idx_q2, model.idx_q1] .= model.grad_sim.grad.∂q3∂q1[1]
	dx[model.idx_q2, model.idx_q2] .= model.grad_sim.grad.∂q3∂q2[1]

	return dx
end

function fu(du, model::ImplicitDynamics, x, u, w)
	q1 = @views x[model.idx_q1]
	q2 = @views x[model.idx_q2]
	model.v1 .= q2
	model.v1 .-= q1
	model.v1 ./= model.grad_sim.h

	RoboDojo.step!(model.grad_sim, q2, model.v1, u, 1)

	du[model.idx_q2, :] .= model.grad_sim.grad.∂q3∂u1[1]

	return du
end

function state_to_configuration(x::Vector{Vector{T}}) where T
	H = length(x)
	n = length(x[1])
	nq = convert(Int, floor(length(x[1]) / 2))
	q = Vector{T}[]

	for t = 1:H
		if t == 1
			push!(q, x[t][1:nq])
		end
		push!(q, x[t][nq .+ (1:nq)])
	end

	return q
end

# using BenchmarkTools
# using InteractiveUtils
# x = x̄[1]
# u = ū[1]
# w_ = w[1]
# d = zeros(nx)
# dx = zeros(nx, nx)
# du = zeros(nx, nu)

# f(d, im_dyn, x, u, w_)
# @benchmark f($d, $im_dyn, $x, $u, $w_)
# @code_warntype f(d, im_dyn, x, u, w_)

# fx(dx, im_dyn, x, u, w_)
# @benchmark fx($dx, $im_dyn, $x, $u, $w_)
# @code_warntype fx(dx, im_dyn, x, u, w_)

# fu(du, im_dyn, x, u, w_)
# @benchmark fu($du, $im_dyn, $x, $u, $w_)
# @code_warntype fu(du, im_dyn, x, u, w_)


text:
# Translation of minpack subroutine cstep
# Dianne O'Leary   July 1991
#
# Subroutine cstep
#
# The purpose of cstep is to compute a safeguarded step for
# a linesearch and to update an interval of uncertainty for
# a minimizer of the function.
#
# The parameter stx contains the step with the least function
# value. The parameter stp contains the current step. It is
# assumed that the derivative at stx is negative in the
# direction of the step. If bracketed is set true then a
# minimizer has been bracketed in an interval of uncertainty
# with endpoints stx and sty.
#
# The subroutine statement is
#
# subroutine cstep(stx, fx, dgx,
#                  sty, fy, dgy,
#                  stp, f, dp,
#                  bracketed, stpmin, stpmax, info)
#
# where
#
# stx, fx, and dgx are variables which specify the step,
#   the function, and the derivative at the best step obtained
#   so far. The derivative must be negative in the direction
#   of the step, that is, dgx and stp-stx must have opposite
#   signs. On output these parameters are updated appropriately
#
# sty, fy, and dgy are variables which specify the step,
#   the function, and the derivative at the other endpoint of
#   the interval of uncertainty. On output these parameters are
#   updated appropriately
#
# stp, f, and dp are variables which specify the step,
#   the function, and the derivative at the current step.
#   If bracketed is set true then on input stp must be
#   between stx and sty. On output stp is set to the new step
#
# bracketed is a logical variable which specifies if a minimizer
#   has been bracketed. If the minimizer has not been bracketed
#   then on input bracketed must be set false. If the minimizer
#   is bracketed then on output bracketed is set true
#
# stpmin and stpmax are input variables which specify lower
#   and upper bounds for the step
#
# info is an integer output variable set as follows:
#   If info = 1,2,3,4,5, then the step has been computed
#   according to one of the five cases below. Otherwise
#   info = 0, and this indicates improper input parameters
#
# Argonne National Laboratory. MINPACK Project. June 1983
# Jorge J. More', David J. Thuente

function cstep(stx::Real, fx::Real, dgx::Real,
               sty::Real, fy::Real, dgy::Real,
               stp::Real, f::Real, dg::Real,
               bracketed::Bool, stpmin::Real, stpmax::Real)

   info = 0

   #
   # Check the input parameters for error
   #

   if (bracketed && (stp <= min(stx, sty) || stp >= max(stx, sty))) ||
        dgx * (stp - stx) >= 0.0 || stpmax < stpmin
      throw(ArgumentError("Minimizer not bracketed"))
   end

   #
   # Determine if the derivatives have opposite sign
   #

   sgnd = dg * (dgx / abs(dgx))

   #
   # First case. A higher function value.
   # The minimum is bracketed. If the cubic step is closer
   # to stx than the quadratic step, the cubic step is taken,
   # else the average of the cubic and quadratic steps is taken
   #

   if f > fx
      info = 1
      bound = true
      theta = 3 * (fx - f) / (stp - stx) + dgx + dg
      s = max(theta, dgx, dg)
      gamma = s * sqrt((theta / s)^2 - (dgx / s) * (dg / s))
      if stp < stx
          gamma = -gamma
      end
      p = gamma - dgx + theta
      q = gamma - dgx + gamma + dg
      r = p / q
      stpc = stx + r * (stp - stx)
      stpq = stx + ((dgx / ((fx - f) / (stp - stx) + dgx)) / 2) * (stp - stx)
      if abs(stpc - stx) < abs(stpq - stx)
         stpf = stpc
      else
         stpf = stpc + (stpq - stpc) / 2
      end
      bracketed = true

   #
   # Second case. A lower function value and derivatives of
   # opposite sign. The minimum is bracketed. If the cubic
   # step is closer to stx than the quadratic (secant) step,
   # the cubic step is taken, else the quadratic step is taken
   #

   elseif sgnd < 0.0
      info = 2
      bound = false
      theta = 3 * (fx - f) / (stp - stx) + dgx + dg
      s = max(theta, dgx, dg)
      gamma = s * sqrt((theta / s)^2 - (dgx / s) * (dg / s))
      if stp > stx
         gamma = -gamma
      end
      p = gamma - dg + theta
      q = gamma - dg + gamma + dgx
      r = p / q
      stpc = stp + r * (stx - stp)
      stpq = stp + (dg / (dg - dgx)) * (stx - stp)
      if abs(stpc - stp) > abs(stpq - stp)
         stpf = stpc
      else
         stpf = stpq
      end
      bracketed = true

   #
   # Third case. A lower function value, derivatives of the
   # same sign, and the magnitude of the derivative decreases.
   # The cubic step is only used if the cubic tends to infinity
   # in the direction of the step or if the minimum of the cubic
   # is beyond stp. Otherwise the cubic step is defined to be
   # either stpmin or stpmax. The quadratic (secant) step is also
   # computed and if the minimum is bracketed then the the step
   # closest to stx is taken, else the step farthest away is taken
   #

   elseif abs(dg) < abs(dgx)
      info = 3
      bound = true
      theta = 3 * (fx - f) / (stp - stx) + dgx + dg
      s = max(theta, dgx, dg)
      #
      # The case gamma = 0 only arises if the cubic does not tend
      # to infinity in the direction of the step
      #
      gamma = s * sqrt(max(0.0, (theta / s)^2 - (dgx / s) * (dg / s)))
      if stp > stx
          gamma = -gamma
      end
      p = gamma - dg + theta
      q = gamma + dgx - dg + gamma
      r = p / q
      if r < 0.0 && gamma != 0.0
         stpc = stp + r * (stx - stp)
      elseif stp > stx
         stpc = stpmax
      else
         stpc = stpmin
      end
      stpq = stp + (dg / (dg - dgx)) * (stx - stp)
      if bracketed
         if abs(stp - stpc) < abs(stp - stpq)
            stpf = stpc
         else
            stpf = stpq
         end
      else
         if abs(stp - stpc) > abs(stp - stpq)
            stpf = stpc
         else
            stpf = stpq
         end
      end

   #
   # Fourth case. A lower function value, derivatives of the
   # same sign, and the magnitude of the derivative does
   # not decrease. If the minimum is not bracketed, the step
   # is either stpmin or stpmax, else the cubic step is taken
   #

   else
      info = 4
      bound = false
      if bracketed
         theta = 3 * (f - fy) / (sty - stp) + dgy + dg
         s = max(theta, dgy, dg)
         gamma = s * sqrt((theta / s)^2 - (dgy / s) * (dg / s))
         if stp > sty
             gamma = -gamma
         end
         p = gamma - dg + theta
         q = gamma - dg + gamma + dgy
         r = p / q
         stpc = stp + r * (sty - stp)
         stpf = stpc
      elseif stp > stx
         stpf = stpmax
      else
         stpf = stpmin
      end
   end

   #
   # Update the interval of uncertainty. This update does not
   # depend on the new step or the case analysis above
   #

   if f > fx
      sty = stp
      fy = f
      dgy = dg
   else
      if sgnd < 0.0
         sty = stx
         fy = fx
         dgy = dgx
      end
      stx = stp
      fx = f
      dgx = dg
   end

   #
   # Compute the new step and safeguard it
   #

   stpf = min(stpmax, stpf)
   stpf = max(stpmin, stpf)
   stp = stpf
   if bracketed && bound
      if sty > stx
         stp = min(stx + 0.66 * (sty - stx), stp)
      else
         stp = max(stx + 0.66 * (sty - stx), stp)
      end
   end

   return stx, fx, dgx, sty, fy, dgy, stp, f, dg, bracketed, info
end


text:
using CleanCMB
using Healpix
using Optim, LinearAlgebra
using Printf, CSV, DataFrames
using Mmap
using Statistics
using Tables
dir = "../data"
# %% Simulation parameters
nrz = 1500
ℓmin, ℓmax = 30, 260 # ℓ range for fitting
Alens = 1
Δℓ = 10 # multipole binning size
# %% Foreground cleaning parameters
showβ = false # Show the fitted foreground parameters (for the smoothed covariance)?
β0 = βs0, βd0, Td0 = [-3.0, 1.6, 19.6] # starting foreground parameters for minimisation of -log(likelihood) by `res = optimize(func, [βs0, βd0, Td0])`
ℓswitch = 30 # the multipole below which the foreground parameters are fitted for each band-power
smooth_FWHM = 0.5 # smoothing for the covariance matrix in units of degrees
σβ = σβs, σβd, σTd = [0.5, 0.5, 5] # Gaussian priors: -2*log(likelihood) = (x[1] - βs0)^2 / σβs^2 + (x[2] - βd0)^2 / σβd^2 + (x[3] - Td0)^2 / σTd^2
Σβ = Diagonal(σβ .^ 2)
lnlike_fgprior(x) = -0.5 * (x .- β0)' * (Σβ \ (x .- β0))
# %% Specification of the experiments
ν = [27, 39, 93, 145, 225, 280, 350, 410, 850] # in GHz
νused = [true, true, true, true, true, true, true, true, true] # Which frequencies to use for fitting?
#νused = [true, true, true, true, true, true, false, false, false]
nν = length(ν)
kν = findall(x -> x == true, νused)
FWHM = [91, 63, 30, 17, 11, 9, 0.58, 0.5, 0.23] # in arcmin
σ = FWHM * π / 10800 / √(8 * log(2)) # in radians
uKarcmin = [35, 21, 2.6, 3.3, 6.3, 16, 105, 372, 5.7e5] # in μK arcmin (for temperature; x√2 for pol)
lknee = [30, 30, 50, 50, 70, 100, 700, 700, 700]
αknee = [-2.4, -2.4, -2.5, -3, -3, -3, -1.4, -1.4, -1.4]
# polarisation beam, Eq.(5.8) of Ng & Liu, Int.J.Mod.Phys.D, 8, 61 (1999)
bPl(ℓ, σb) = ifelse(ℓ >= 2, exp(-(ℓ * (ℓ + 1) - 4) * σb^2 / 2), 0)
# %% Read in the hits map and calculate fsky need for the likelihood
nhitsfile = joinpath(dir, "nhits_SAT_r7.FITS")
nhits = readMapFromFITS(nhitsfile, 1, Float64)
nside = nhits.resolution.nside
nhits /= maximum(nhits)
fsky = mean(nhits)^2 / mean(nhits .^ 2)
# %% Read in binned theory power spectra
clt_th = CSV.read(joinpath(dir, "tensor_eebb_binned.csv"), DataFrame)
cls_th = CSV.read(joinpath(dir, "scalar_eebb_binned.csv"), DataFrame)
ell_eff = clt_th.leff
nbands = length(ell_eff)
rclass = 0.01
# %% Read in covariance matrices and perform parametric maxiimum likelihood fitting
io = open(joinpath(dir, "cov_bb_foreg.dat"), "r")
cov3 = Mmap.mmap(io, Array{Float64,3}, (nν, nν, nbands))
close(io)
cl1, cl2, cl3 = zeros(nbands, nrz), zeros(nbands, nrz), zeros(nbands, nrz)
for irz = 1:nrz
    ioin = open(joinpath(dir, @sprintf("cov_bb_total_irz%04d.dat", irz)), "r")
    cov1 = Mmap.mmap(ioin, Array{Float64,3}, (nν, nν, nbands))
    close(ioin)
    ioin = open(joinpath(dir, @sprintf("cov_bb_noise_irz%04d.dat", irz)), "r")
    cov2 = Mmap.mmap(ioin, Array{Float64,3}, (nν, nν, nbands))
    close(ioin)
    nij = zeros(nν, nν, nbands) # Noise covariance matrix for the ML analysis
    for iν = 1:nν
        ell_fact = ell_eff .* (ell_eff .+ 1) / 2π
        nij[iν, iν, :] =
            2 *
            (uKarcmin[iν] * π / 10800)^2 *
            (1 .+ (ell_eff / lknee[iν]) .^ αknee[iν]) .* ell_fact .*
            bPl.(ell_eff, σ[iν]) .^ -2
    end
    ## Use Parametric Maximum Likelihood method to obtain power spectra of clean maps of the CMB
    A(x) = [ones(length(kν)) synch.(ν[kν], βs = x[1]) dust1.(
        ν[kν],
        βd = x[2],
        Td = x[3],
    )] # Frequency response matrix
    # For ℓ > ℓswitch: Apply parametric maximum likelihood method using a smoothed covariance matrix.
    func_sum(x, cij) =
        -lnlike_fgprior(x) -
        fsky *
        Δℓ *
        sum(
            (2 * ell_eff[jb] + 1) *
            loglike_beta(nij[kν, kν, jb], A(x), cij[kν, kν, jb]) for
            jb = 1:nbands
        ) # -log(likelihood) to minimise by `optimize`
    # Smooth covariance matrices to `smooth_FWHM` resolution
    cij = zeros(nν, nν, nbands)
    σsmo = smooth_FWHM * π / 180 / sqrt(8 * log(2))
    for ib = 1:nbands
        cij[:, :, ib] = cov1[:, :, ib] * bPl(ell_eff[ib], σsmo)^2
    end
    res = optimize(x -> func_sum(x, cij), β0)
    if showβ
        println("Fitted parameters: B-mode")
        @show res
    end
    β = Optim.minimizer(res)
    B = [synch.(ν[kν], βs = β[1]) dust1.(ν[kν], βd = β[2], Td = β[3])] # Best-fitting frequency response of the FG
    w_smooth = milca_weights(cov1[kν, kν, :], ones(length(kν)), B)
    cl1[:, irz] = ilc_clean_cij(cov1[kν, kν, :], w_smooth)
    cl2[:, irz] = ilc_clean_cij(cov2[kν, kν, :], w_smooth)
    cl3[:, irz] = ilc_clean_cij(cov3[kν, kν, :], w_smooth)
    # For for ℓ ≤ ℓswitch: Apply parametric maximum likelihood method for each band-power.
    iib = findall(x -> x ≤ ℓswitch, ell_eff)
    if iib ≠ []
        for ib = 1:maximum(iib)
            func(x, cij) =
                -lnlike_fgprior(x) -
                fsky *
                (2 * ell_eff[ib] + 1) *
                loglike_beta(nij[kν, kν, ib], A(x), cij[kν, kν, ib])
            res = optimize(x -> func(x, cov1), β0)
            #@show res
            β = Optim.minimizer(res)
            B = [synch.(ν[kν], βs = β[1]) dust1.(ν[kν], βd = β[2], Td = β[3])]
            w = milca_weights(cov1[kν, kν, ib], ones(length(kν)), B)
            cl1[ib, irz] = ilc_clean_cij(cov1[kν, kν, ib], w)
            cl2[ib, irz] = ilc_clean_cij(cov2[kν, kν, ib], w)
            cl3[ib, irz] = ilc_clean_cij(cov3[kν, kν, ib], w)
        end
    end
end
m1, m2, m3 = mean(cl1, dims = 2), mean(cl2, dims = 2), mean(cl3, dims = 2)
v1 = var(cl1, dims = 2)

# %% Joint fit for the tensor-to-scalar ratio and the foreground amplitude
ii = findall(x -> x >= ℓmin && x <= ℓmax, ell_eff)
r, rw = zeros(nrz), zeros(nrz, 2)
y1 = clt_th.bb[ii] / rclass
y2 = m3[ii]
v = v1[ii]
Fij = [
    sum(y1 .* y1 ./ v) sum(y1 .* y2 ./ v)
    sum(y2 .* y1 ./ v) sum(y2 .* y2 ./ v)
] # 2x2 Fisher matrix for the tensor-to-scalar ratio and the foreground amplitude
Cij = inv(Fij) # Covariance matrix
for irz = 1:nrz
    x = cl1[ii, irz] .- m2[ii] .- Alens * cls_th.bb[ii]
    r[irz] = sum(x .* y1 ./ v) / sum(y1 .^ 2 ./ v)
    z = [sum(x .* y1 ./ v), sum(x .* y2 ./ v)]
    rw[irz, 1:2] = Fij \ z
end
println("Fitted ℓs: ", ell_eff[ii])
println("Used νs: ", ν[kν], " GHz")
println("Without foreground marginalisation:")
println("- r = ", mean(r), " ± ", std(r))
println("- Fisher error = ", 1 / √sum(y1 .^ 2 ./ v))
println("With foreground marginalisation:")
println("- r = ", mean(rw[:, 1]), " ± ", std(rw[:, 1]))
println("- Fisher error = ", √Cij[1, 1])
# println("- FG = ", mean(rw[:, 2]), " ± ", std(rw[:, 2]))
# println("- Fisher error = ", √Cij[2, 2])
# %% Save to the file
t = Tables.table([1:nrz r rw[:, 1]])
CSV.write("milca_results.csv", t, header = ["irz", "r_wo_FGmarg", "r_w_FGmarg"])


text:
"""
    perpendicular_vector(vec)

Compute a vector perpendicular to `vec` by switching the two elements with
largest absolute value, flipping the sign of the second largest, and setting the
remaining element to zero.
"""
function perpendicular_vector(vec::SVector{3})
    T = eltype(vec)

    # find indices of the two elements of vec with the largest absolute values:
    absvec = abs.(vec)
    ind1 = argmax(absvec) # index of largest element
    tmin = typemin(T)
    @inbounds absvec2 = @SVector [ifelse(i == ind1, tmin, absvec[i]) for i = 1 : 3] # set largest element to typemin(T)
    ind2 = argmax(absvec2) # index of second-largest element

    # perp[ind1] = -vec[ind2], perp[ind2] = vec[ind1], set remaining element to zero:
    @inbounds perpind1 = -vec[ind2]
    @inbounds perpind2 = vec[ind1]
    tzero = zero(T)
    perp = @SVector [ifelse(i == ind1, perpind1, ifelse(i == ind2, perpind2, tzero)) for i = 1 : 3]
end

@noinline length_error(v, len) =
    throw(DimensionMismatch("Expected length $len, got length $(length(v))"))

@inline function check_length(v, len)
    if length(v) != len
        length_error(v, len)
    end
end

function skew(v::AbstractVector)
    check_length(v, 3)
    @SMatrix [0   -v[3]  v[2];
              v[3] 0    -v[1];
             -v[2] v[1]  0]
end

function vee(S::AbstractMatrix)
    return @SVector [S[3,2], S[1,3], S[2,1]]
end

"""
The element type for a rotation matrix with a given angle type is composed of
trigonometric functions of that type.
"""
Base.@pure rot_eltype(angle_type) = typeof(sin(zero(angle_type)))


text:
#!/usr/bin/env julia
###############################################################################
# Params.

needed_quantity = 10

#TODO Get from the internet.
currencies_to_eur = Dict(
	"EUR" => 1,
	"GBP" => 1.1845,
	"USD" => 0.9576
)

simple_speed_grades = false

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

using Tables
using Octopart

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

families = read_families()

println("Getting offers from Octopart...")

@time offers = get_offers_for_FPGAs(
	families,
	needed_quantity,
	currencies_to_eur,
	simple_speed_grades
)

println("Done!")

# Save offers with params to XLS.
write_table("tmp/offers.xls", "offers", offers)

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

println("End")

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