text:
function analyzescenarios(Ψ::Simulation,
  τscenarios::Vector{Vector{Float64}},
  τscenarionames::Vector{Symbol};
  prefix::String = "",
  outpaper::String=OUT_PAPER,
  extension::String = "pdf")

  local x::UnitRange = 1:Ψ.nsims
  local y::Matrix{Float64}

  gr()
  nscenarios::Int = length(τscenarios)

  highproductionplots = Vector()
  lowproductionplots = Vector()
  investmentplots = Vector()
  capitalplots = Vector()
  valueplots = Vector()

  for i ∈ 1:nscenarios
    Random.seed!(11)
    updatefirmτ!(Ψ.Φ, τscenarios[i])
    Simulate!(Ψ)

    legend::Union{Bool, Symbol} = i==1 ? :topright : false
    xlab::String = i==nscenarios ? "period" : ""

    p = plot(x, Ψ.Yh, lab = "high state", legend=legend)
    plot!(p, x, Ψ.S , lab = "economic state (H/L)", line=:dot)
    #plot!(p, x, Ψ.Y)
    xaxis!(p, xlab)#"period")
    yaxis!(p, "production")
    plot!(p, title = "High State Allocations:  $(τscenarionames[i])")
    push!(highproductionplots, p)

    p = plot(x, Ψ.Yl, lab = "low state", legend=legend)
    plot!(p, x, Ψ.S , lab = "economic state (H/L)", line=:dot)
    #plot!(p, x, Ψ.Y)
    xaxis!(p, xlab)
    yaxis!(p, "production")
    plot!(p, title = "Low State Allocations:  $(τscenarionames[i])")
    push!(lowproductionplots, p)

    p = plot(x, Ψ.I, lab = "investment", legend=legend)
    plot!(p, x, Ψ.S , lab = "economic state (H/L)", line=:dot)
    #plot!(p, x, Ψ.Y)
    xaxis!(p, xlab)
    yaxis!(p, "investment")
    plot!(p, title = "Firm Investment:  $(τscenarionames[i])")
    push!(investmentplots, p)

    p = plot(x, Ψ.K, lab = "capital", legend=legend)
    plot!(p, x, Ψ.S , lab = "economic state (H/L)", line=:dot)
    #plot!(p, x, Ψ.Y)
    xaxis!(p, xlab)
    yaxis!(p, "capital")
    plot!(p, title = "Firm Capital:  $(τscenarionames[i])")
    push!(capitalplots, p)

    p = plot(x, Ψ.V, lab = "firm value", legend=legend)
    #plot!(twinx(p), x, Ψ.S , lab = "economic state (H/L)", line=:dot)
    plot!(p, x, Ψ.S , lab = "economic state (H/L)", line=:dot)
    #plot!(p, x, Ψ.Y)
    xaxis!(p, xlab)
    yaxis!(p, "firm value")
    plot!(p, title = "Firm Value:  $(τscenarionames[i])")
    push!(valueplots, p)
  end

  statespervar::Int = Int(round(sqrt(Ψ.Φ.N)))

  outproduction = plot(highproductionplots..., layout= (nscenarios,1), size=(1000, 900))
  savefig(outproduction, "$outpaper\\$(prefix)highscen_s$(statespervar)_n$(Ψ.nsims).$extension")

  outproduction = plot(lowproductionplots..., layout= (nscenarios,1), size=(1000, 900))
  savefig(outproduction, "$outpaper\\$(prefix)lowscen_s$(statespervar)_n$(Ψ.nsims).$extension")

  outproduction = plot(investmentplots..., layout= (nscenarios,1), size=(1000, 900))
  savefig(outproduction, "$outpaper\\$(prefix)investment_s$(statespervar)_n$(Ψ.nsims).$extension")

  outproduction = plot(capitalplots..., layout= (nscenarios,1), size=(1000, 900))
  savefig(outproduction, "$outpaper\\$(prefix)capital_s$(statespervar)_n$(Ψ.nsims).$extension")

  outproduction = plot(valueplots..., layout= (nscenarios,1), size=(1000, 900))
  savefig(outproduction, "$outpaper\\$(prefix)value_s$(statespervar)_n$(Ψ.nsims).$extension")
end

mutable struct τmeans
  τvals::Vector{Float64}
  highproductionmean::Vector{Float64}
  lowproductionmean::Vector{Float64}
  investmentmean::Vector{Float64}
  capitalmean::Vector{Float64}
  valuemean::Vector{Float64}
end

#this method does the tax analysis
function analyzetaxes(Ψ::Simulation,
  τvals::Vector{Float64};
  outpath = OUT_PATH, prefix::String = "")

  #these will hold the means
  nvals::Int = length(τvals)
  highproductionmean = Vector{Float64}(undef, nvals)
  lowproductionmean = Vector{Float64}(undef, nvals)
  investmentmean = Vector{Float64}(undef, nvals)
  capitalmean = Vector{Float64}(undef, nvals)
  valuemean = Vector{Float64}(undef, nvals)

  statespervar::Int = Int(round(sqrt(Ψ.Φ.N)))
  for i ∈ 1:nvals
    Random.seed!(11)
    updatefirmτ!(Ψ.Φ, [τvals[i], 0.0])
    Simulate!(Ψ)

    highproductionmean[i] = mean(Ψ.Yh)
    lowproductionmean[i] = mean(Ψ.Yl)
    investmentmean[i]  = mean(Ψ.I)
    capitalmean[i]  = mean(Ψ.K)
    valuemean[i]  = mean(Ψ.V)
  end

  writeto::String = "$outpath\\$(prefix)highstatetax_s$(statespervar)_n$(Ψ.nsims).jls"
  T::τmeans = τmeans(τvals, highproductionmean, lowproductionmean, investmentmean, capitalmean, valuemean)
  serialize(writeto, T)
  println("Wrote tax summary data to $writeto")

  for i ∈ 1:nvals
    Random.seed!(11)
    updatefirmτ!(Ψ.Φ, [0.0, τvals[i]])
    Simulate!(Ψ)

    highproductionmean[i] = mean(Ψ.Yh)
    lowproductionmean[i] = mean(Ψ.Yl)
    investmentmean[i]  = mean(Ψ.I)
    capitalmean[i]  = mean(Ψ.K)
    valuemean[i]  = mean(Ψ.V)
  end

  writeto = "$outpath\\$(prefix)lowstatetax_s$(statespervar)_n$(Ψ.nsims).jls"
  T= τmeans(τvals, highproductionmean, lowproductionmean, investmentmean, capitalmean, valuemean)
  serialize(writeto, T)
  println("Wrote tax summary data to $writeto")
end

function graphtaxanalysis(Ψ::Simulation;
  outpath = OUT_PATH,
  prefix::String = "",
  statespervar::Int = Int(round(sqrt(Ψ.Φ.N))),
  nsims::Int = Ψ.nsims,
  outpaper::String=OUT_PAPER,
  extension::String = "pdf"
  )

  gr()

  local plotholder::Vector = Vector()
  local plotnames::Vector{String} = Vector{String}()

  writefrom::String = "$outpath\\$(prefix)highstatetax_s$(statespervar)_n$(Ψ.nsims).jls"
  Thigh::τmeans = deserialize(writefrom)
  writefrom = "$outpath\\$(prefix)lowstatetax_s$(statespervar)_n$(Ψ.nsims).jls"
  Tlow = deserialize(writefrom)

  Ts::Vector{τmeans} = [Thigh, Tlow]
  taxstatenames::Vector{String} = ["high_state_taxes", "low_state_taxes"]

  for i ∈ 1:length(Ts)
    legend::Union{Bool, Symbol} = false
    xlab = "Tax Rates"

    p = plot(Ts[i].τvals, Ts[i].highproductionmean, lab = "high production alloc.", legend=:topright, size=(1000, 900))
    plot!(p, Ts[i].τvals, Ts[i].lowproductionmean , lab = "low production alloc.")
    xaxis!(p, xlab)#
    yaxis!(p, "mean production")
    plot!(p, title = "$(taxstatenames[i]): High and Low State Production")
    push!(plotholder, p)
    push!(plotnames, "production_$(taxstatenames[i])" )

    p = plot(Ts[i].τvals, Ts[i].capitalmean, lab = "capital", legend=legend, size=(1000, 900))
    xaxis!(p, xlab)#
    yaxis!(p, "mean capital")
    plot!(p, title = "$(taxstatenames[i]): Capital")
    push!(plotholder, p)
    push!(plotnames, "capital_$(taxstatenames[i])" )

    p = plot(Ts[i].τvals, Ts[i].investmentmean, lab = "investment", legend=legend, size=(1000, 900))
    xaxis!(p, xlab)#
    yaxis!(p, "mean investment")
    plot!(p, title = "$(taxstatenames[i]): Investment")
    push!(plotholder, p)
    push!(plotnames, "capital_$(taxstatenames[i])" )

    p = plot(Ts[i].τvals, Ts[i].valuemean, lab = "value", legend=legend, size=(1000, 900))
    xaxis!(p, xlab)#
    yaxis!(p, "mean value")
    plot!(p, title = "$(taxstatenames[i]): Value")
    push!(plotholder, p)
    push!(plotnames, "value_$(taxstatenames[i])" )

    #=p = plot(Ts[i].τvals, Ts[i].lowproductionmean, lab = "capital", legend=legend, size=(1000, 900))
    xaxis!(p, xlab)#
    yaxis!(p, "mean production")
    plot!(p, title = "$(taxstatenames[i]): Capital")
    push!(plotholder, p)
    push!(plotnames, "value_$(taxstatenames[i])" )=#
  end

  outplotholder = Vector()

  push!(outplotholder, plot(plotholder[1:4]...,  layout= (2,2)))
  push!(outplotholder, plot(plotholder[5:8]...,  layout= (2,2)))
  for i ∈ 1:length(outplotholder)
    savefig(outplotholder[i], "$outpaper\\$(prefix)$(taxstatenames[i])_s$(statespervar)_n$(Ψ.nsims).$extension")
  end

end

function graphscript!(Ψ::Simulation;
  τvals = Vector{Float64}(),
  τscenarios::Vector{Vector{Float64}} = Vector{Vector{Float64}}(),
  τscenarionames::Vector{Symbol} = Vector{Symbol}(),
  scenariographs::Bool = true,
  prefix::String = "", taxanalysis::Bool=true, taxgraphs::Bool= true)

  scenariographs && analyzescenarios(Ψ, τscenarios, τscenarionames,prefix=prefix)
  taxanalysis && analyzetaxes(Ψ::Simulation, τvals, prefix=prefix)
  taxgraphs && graphtaxanalysis(Ψ::Simulation, prefix=prefix)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module ACE1pack

# load and reexport JuLIP, ACE1, IPFitting
using Reexport
@reexport using JuLIP
@reexport using ACE1
@reexport using IPFitting
export JuLIP, ACE1, IPFitting

# Convenience Layer

include("fit.jl")

include("data.jl")

include("basis.jl")

include("solver.jl")

include("regularizer.jl")

include("read_params.jl")

end

ext: jl


text:
module IntegerMask

export Mask

import Base: ==, \

type Mask{N}
    indices::Vector{Int}
end

Mask{N}(indices::Int...) where N = Mask{N}(Int[indices...])

Base.copy{N}(m::Mask{N}) = Mask{N}(copy(m.indices))

Base.length(m::Mask) = length(m.indices)

=={N}(m1::Mask{N}, m2::Mask{N}) = m1.indices == m2.indices

function Base.getindex{N}(a::AbstractArray, m::Mask{N})
    @assert length(a) == N
    @inbounds result = a[m.indices]
    result
end

function Base.getindex{N}(m1::Mask{N}, m2::Mask{N})
    _indices = Int[]
    @inbounds for i=1:length(m2.indices)
        if m2.indices[i] ∈ m1.indices
            push!(_indices, i)
        end
    end
    Mask{length(m2.indices)}(_indices)
end

function Base.transpose{N}(m::Mask{N})
    L = length(m.indices)
    x = Vector{Int}(N - L)
    i_ = 1 # Position in the x vector
    j = 1 # Position in indices vector
    @inbounds for i=1:N
        if j > L || m.indices[j]!=i
            x[i_] = i
            i_ += 1
        else
            j += 1
        end
    end
    Mask{N}(x)
end

function Base.union{N}(m1::Mask{N}, m2::Mask{N})
    xvec = Vector{Int}()
    result = Mask{N}(xvec)
    ind1 = m1.indices
    ind2 = m2.indices
    i1 = 1
    i2 = 1
    N1 = length(ind1)
    N2 = length(ind2)
    @inbounds while true
        if i1 > N1
            for j=i2:N2
                push!(xvec, ind2[j])
            end
            return result
        elseif i2 > N2
            for j=i1:N1
                push!(xvec, ind1[j])
            end
            return result
        end
        x1 = ind1[i1]
        x2 = ind2[i2]
        if x1 == x2
            i1 += 1
            i2 += 1
            push!(xvec, x1)
        elseif x1 < x2
            i1 += 1
            push!(xvec, x1)
        else
            i2 += 1
            push!(xvec, x2)
        end
    end
end

function Base.intersect{N}(m1::Mask{N}, m2::Mask{N})
    xvec = Vector{Int}()
    result = Mask{N}(xvec)
    ind1 = m1.indices
    ind2 = m2.indices
    i1 = 1
    i2 = 1
    N1 = length(ind1)
    N2 = length(ind2)
    if i1 > N1 || i2 > N2
        return result
    end
    x1 = ind1[i1]
    x2 = ind2[i2]
    @inbounds while true
        if x1 == x2
            i1 += 1
            i2 += 1
            push!(xvec, x1)
            if i1 > N1 || i2 > N2
                return result
            end
            x1 = ind1[i1]
            x2 = ind2[i2]
        elseif x1 < x2
            i1 += 1
            if i1 > N1
                return result
            end
            x1 = ind1[i1]
        else
            i2 += 1
            if i2 > N2
                return result
            end
            x2 = ind2[i2]
        end
    end
end

function Base.setdiff{N}(m1::Mask{N}, m2::Mask{N})
    Mask{N}(setdiff(m1.indices, m2.indices))
end

\{N}(m1::Mask{N}, m2::Mask{N}) = setdiff(m1, m2)

Base.in{N}(i::Int, m2::Mask{N}) = in(i, m2.indices)

function Base.prod{N}(a::AbstractArray, m::Mask{N})
    @assert length(a) == N
    c = one(eltype(a))
    @inbounds for i in m.indices
        c *= a[i]
    end
    c
end

function Base.sum{N}(a::AbstractArray, m::Mask{N})
    @assert length(a) == N
    c = zero(eltype(a))
    @inbounds for i in m.indices
        c += a[i]
    end
    c
end

function Base.reduce{N}(f::Function, a::AbstractArray, m::Mask{N})
    @assert length(a) == N
    l = length(m)
    @assert l > 0
    c = a[m.indices[1]]
    for i in 2:l
        c = f(c, a[m.indices[i]])
    end
    c
end

end #module

max_issues_repo_name: bastikr/IntegerMask

max_forks_repo_name: bastikr/IntegerMask


text:
using Plots
using Distributed
function evaluate_decision_wrapper(p, decision, scenario)
    cost = 0.
    try cost = evaluate_decision(p, decision, scenario)
    catch e
        cost = Inf
    end
    return cost
end

function test_decision(p, timesteps)
    decision = optimal_decision(p)
    L = length(timesteps)
    cost_pos_flex = zeros(L)
    cost_neg_flex = zeros(L)
    potential_pos_flex = zeros(L)
    potential_neg_flex = zeros(L)
    @distributed for t in timesteps
        potential_pos_flex[t], cost_pos_flex[t] = find_f_max(p,t,1,decision,10.)
        potential_neg_flex[t], cost_neg_flex[t] = find_f_max(p,t,-1,decision,10.)
    end

    return cost_pos_flex, potential_pos_flex, cost_neg_flex, potential_neg_flex
end

function plot_flexibility(timesteps, cost_pos_flex, potential_pos_flex, cost_neg_flex, potential_neg_flex, obj_value)
    plt_cost = plot()
    plt_pot = plot()
    plot!(plt_cost, timesteps, (cost_pos_flex .- obj_value)./ potential_pos_flex, label = "price of positive flexibility")
    plot!(plt_cost, timesteps, (cost_neg_flex .- obj_value)./ potential_neg_flex, label = "price of negative flexibility")
    plot!(plt_pot, timesteps, potential_pos_flex, fillrange = 0, fillalpha = 0.35, label = "positive flexibility potential")
    plot!(plt_pot, timesteps, potential_neg_flex, fillrange = 0, fillalpha = 0.35, label = "negative flexibility potential")
    display(plot(plt_cost, plt_pot, layout = (2, 1)))
end

function flexibility_availability(plt, flex_potential; plot_options...)
        p_sorted = sort(unique(flex_potential))
        fraction = 1 .-[sum(abs.(flex_potential[:]).<=abs(f)) for f in p_sorted]./length(flex_potential)
        plot!(plt, p_sorted, fraction)
    display(plt_av)
end

function find_f_max(p,t,s,od,tol; maxiter = 100)
    a = 0.
    b = 10000.
    i = 0
    scen = @scenario t_xi = t s_xi = s F_xi = 0. probability = 1.
    cost = evaluate_decision_wrapper(p,od,scen)
    if cost == Inf
        return 0.,cost
    else
        while b-a>tol && i<=maxiter
            i+=1
            scen = @scenario t_xi = t s_xi = s F_xi = b probability = 1.
            cost = evaluate_decision_wrapper(p,od,scen)
            if cost == Inf
                b = (a+b)/2
            else
                a = b
                b *= 2
            end
        end
        if cost == Inf
            scen = @scenario t_xi = t s_xi = s F_xi = b-tol probability = 1.
            cost = evaluate_decision_wrapper(p,od,scen)
        end
        println(i)
        return s*b,cost
    end
end


text:
# For each Salsa call that is defined, there will be a unique instance of `AbstractKey` that
# identifies requests for the cached value (ie accessing inputs or calling derived
# functions). Instances of these key types are used to identify both _which method_ was
# called, and what the arguments to that call were.
#
# Dependencies between Salsa computations are represented via instances of these keys,
# which specify everything needed to rerun the computation: (key, function-args)
# The stored arguments do not include the Salsa Runtime object itself, since this is always
# present, and changes as computations are performed.
#
# Given:
#   @declare_input input_str(::Int, ::Int) :: String
#   @derived function foo(rt, x::Int, y::Any, z::Number) ... end
#   @derived function foo(rt, x,y,z) ... end
#
# Examples:
#  derived_foo(rt, 2,3) -> DerivedKey{typeof(derived_foo)}((2,3))
#  input_bar(rt, 2,3) -> InputKey{typeof(input_bar)}((2,3))
#  foo(rt,1,2,3)  -> DerivedKey{typeof(foo), Tuple{Int,Any,Number}}((1, 2, 3))
#  input_str(1,2) -> InputKey{typeof(input_str), Tuple{Int,Int}}((1, 2))
abstract type AbstractKey end

struct DerivedKey{F<:Function,TT<:Tuple} <: AbstractKey
    args::TT
end
DerivedKey{F}(args::TT) where {F<:Function,TT<:Tuple} = DerivedKey{F,TT}(args)

# NOTE: After several iterations, the InputKeys are now essentially identical to the
# DerivedKeys. They only differ to allow distinguishing them for dispatch. We might want to
# do some refactoring to share more code below.
struct InputKey{F<:Function,TT<:Tuple} <: AbstractKey
    args::TT
end
InputKey{F}(args::TT) where {F<:Function,TT<:Tuple} = InputKey{F,TT}(args)

# TODO: Probably don't need both this union and the abstract type. Just pick one.
# Convenience Union for sharing code.
const DependencyKey{F,TT} = Union{DerivedKey{F,TT}, InputKey{F,TT}}

# TODO(NHD): Use AutoHashEquals.jl here instead to autogenerate these.
# Note that floats should be compared for equality, not NaN-ness
function Base.:(==)(x1::DK, x2::DK) where DK <: DependencyKey
    return isequal(x1.args, x2.args)
end
function Base.isless(x1::DK, x2::DK) where DK <: DependencyKey
    return isless(x1.args, x2.args)
end

const INPUT_KEY_SEED = hash(codeunits(string(InputKey)))
const DERIVED_KEY_SEED = hash(codeunits(string(DerivedKey)))

function Base.hash(x::InputKey{F}, h::UInt) where {F}
    return hash(x.args, hash(codeunits(string(F)), h + INPUT_KEY_SEED))
end
function Base.hash(x::DerivedKey{F}, h::UInt) where {F}
    return hash(x.args, hash(codeunits(string(F)), h + DERIVED_KEY_SEED))
end

# Override `Base.show` to minimize redundant printing (skip module name).
# Don't print the TT type on the DependencyKey, since it's recovered by the fields.
function Base.show(io::IO, key::InputKey{F,TT}) where {F,TT}
    print(io, "InputKey{$F}($(key.args))")
end
function Base.show(io::IO, key::DerivedKey{F,TT}) where {F,TT}
    print(io, "DerivedKey{$F}($(key.args))")
end

# Pretty-print a DependencyKey for tracing and printing in DerivedFunctionExceptions:
# @input foo(::Runtime, 1::Int, 2::Int)
# @derived foo(::Runtime, 1::Int, 2::Any, 3::Number)
function _print_dep_key_as_call(
    io::IO,
    dependency::DependencyKey{F,TT}
) where {F,TT}
    call_args = dependency.args
    f = isdefined(F, :instance) ? nameof(F.instance) : nameof(F)
    argsexprs = [
        Expr(:(::), :Runtime),
        (Expr(:(::), call_args[i], fieldtype(TT, i)) for i = 1:length(call_args))...,
    ]
    # Display f in a copy/paste executable way, by exploiting Expr printing.
    f_str = string(:(($f,)))[2:end-2]  # (wraps f name in var"" if needed)
    print(io, "$f_str($(join(argsexprs, ", ")))")
end
function Base.print(io::IO, dependency::InputKey)
    print(io, "@input ")
    _print_dep_key_as_call(io, dependency)
end
function Base.print(io::IO, dependency::DerivedKey)
    print(io, "@derived ")
    _print_dep_key_as_call(io, dependency)
end


text:
export FreeflyerSE2
export init_traj_nothing, init_traj_straightline

mutable struct FreeflyerSE2 <: DynamicsModel
  x_dim   # state: r θ v ω
  u_dim
  clearance

  # Parameters that can be updated
  f::Vector
  A
  B
end

function FreeflyerSE2()
  x_dim = 6
  u_dim = 3
  clearance = 0.05
  FreeflyerSE2(x_dim, u_dim, clearance, [], [], [])
end

function SCPParam(model::FreeflyerSE2, fixed_final_time::Bool)
  convergence_threshold = 1.0e-2
  SCPParam(fixed_final_time, convergence_threshold)
end

function SCPParam_GuSTO(model::FreeflyerSE2)
  Δ0 = 3.
  ω0 = 1.
  ω_max = 1.0e10
  ε = 1.0e-2
  ρ0 = 0.1
  ρ1 = 0.3
  β_succ = 2.
  β_fail = 0.5
  γ_fail = 10.

  SCPParam_GuSTO(Δ0, ω0, ω_max, ε, ρ0, ρ1, β_succ, β_fail, γ_fail)
end

function SCPParam_Mao(model::FreeflyerSE2)
  ρ = [0.;0.25;0.9]
  Δ_u0 = 0.1
  λ = 1.
  α = 2.
  SCPParam_Mao(ρ, Δ_u0, λ, α)
end

function SCPParam_TrajOpt(model::FreeflyerSE2)
  μ0 = 1.
  s0 = 1.
  c = 10.
  τ_plus = 2.
  τ_minus = 0.5
  k = 5.
  ftol = 0.01
  xtol = 0.1
  ctol = 0.01
  max_penalty_iteration = 5
  max_convex_iteration = 5
  max_trust_iteration = 5

  SCPParam_TrajOpt(μ0, s0, c, τ_plus, τ_minus, k, ftol, xtol, ctol, max_penalty_iteration,max_convex_iteration,max_trust_iteration)
end

function cost_true(traj, traj_prev::Trajectory, OAP::A) where A <: OptAlgorithmProblem{Freeflyer{T}, FreeflyerSE2, E} where {T,E}
  u_dim = OAP.PD.model.u_dim
  U, N, dtp = traj.U, OAP.N, traj_prev.dt
  Jm = 0

  # Trapezoidal
  for k in 2:N
    Jm += sum(1/2*dtp*(U[j,k-1]^2 + U[j,k]^2) for j = 1:u_dim)
  end
  return Jm
end

function cost_true_convexified(traj, traj_prev::Trajectory, OAP::A) where A <: OptAlgorithmProblem{Freeflyer{T}, FreeflyerSE2, E} where {T,E}
  cost_true(traj, traj_prev, OAP)
end

#############################
# Trajectory Initializations
#############################
function init_traj_nothing(TOP::TrajectoryOptimizationProblem{Freeflyer{T}, FreeflyerSE2, E}) where {T,E}
  model, x_init, goal_set = TOP.PD.model, TOP.PD.x_init, TOP.PD.goal_set
  x_dim, u_dim, N, tf_guess = model.x_dim, model.u_dim, TOP.N, TOP.tf_guess

  goal_final = get_first_goal_at_time(goal_set, tf_guess)
  x_goal = goal_final.params.point

  X = repmat(0.5(x_init + x_goal),1,N)
  U = zeros(u_dim,N)
  Trajectory(X, U, tf_guess)
end

function init_traj_straightline(TOP::TrajectoryOptimizationProblem{Freeflyer{T}, FreeflyerSE2, E}) where {T,E}
  model, x_init, goal_set = TOP.PD.model, TOP.PD.x_init, TOP.PD.goal_set
  x_dim, u_dim, N, tf_guess = model.x_dim, model.u_dim, TOP.N, TOP.tf_guess

  # Set last state to the center of the goal defined latest in time
  t_goal_final = last(goal_set.goals)[1]
  x_goal = zeros(x_dim)
  for goal in values(inclusive(goal_set.goals, searchequalrange(goal_set.goals, t_goal_final)))
    x_goal[goal.ind_coordinates] = center(goal.params)
  end

  X = hcat(range(x_init, stop=x_goal, length=N)...)
  U = zeros(u_dim, N)
  Trajectory(X, U, tf_guess)
end

####################
# Constraint-adding
####################
function initialize_model_params!(SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}, traj_prev::Trajectory) where {T,E}
  N, robot, model = SCPP.N, SCPP.PD.robot, SCPP.PD.model
  x_dim, u_dim = model.x_dim, model.u_dim
  Xp, Up = traj_prev.X, traj_prev.U

  model.f, model.A, model.B = [], A_dyn(Xp[:,1],robot,model), B_dyn(Xp[:,1],robot,model)
  for k = 1:N
    push!(model.f, f_dyn(Xp[:,k],Up[:,k],robot,model))
  end

  SCPP.PD.robot.Jcollision = []
  for k = 1:N
    push!(SCPP.PD.robot.Jcollision, [Eye(2) zeros(2,1); zeros(1,3)])
  end
  update_model_params!(SCPP, traj_prev)
end

function update_model_params!(SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}, traj_prev::Trajectory) where {T,E}
  N, robot, model = SCPP.N, SCPP.PD.robot, SCPP.PD.model
  Xp, Up, f = traj_prev.X, traj_prev.U, model.f

  for k = 1:N
    update_f!(f[k], Xp[:,k], Up[:,k], robot, model)
  end

  xb, Jcollision = SCPP.PD.robot.xb, SCPP.PD.robot.Jcollision
  for k = 1:N
    thetap = traj_prev.X[3,k]
    Jcollision[k][1,3] = -xb[1]*sin(thetap)-xb[2]*cos(thetap)
    Jcollision[k][2,3] =  xb[1]*cos(thetap)-xb[2]*sin(thetap)
  end
end

macro scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)
  quote
    X, U, Tf = $(esc(traj)).X, $(esc(traj)).U, $(esc(traj)).Tf
    Xp, Up, Tfp, dtp = $(esc(traj_prev)).X, $(esc(traj_prev)).U, $(esc(traj_prev)).Tf, $(esc(traj_prev)).dt
    robot, model, WS, x_init, goal_set = $(esc(SCPP)).PD.robot, $(esc(SCPP)).PD.model, $(esc(SCPP)).WS, $(esc(SCPP)).PD.x_init, $(esc(SCPP)).PD.goal_set
    x_dim, u_dim, N, dh = model.x_dim, model.u_dim, $(esc(SCPP)).N, $(esc(SCPP)).dh
    X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh
  end
end

## Dynamics constraints
function dynamics_constraints(traj, traj_prev::Trajectory, SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}, k::Int) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)

  Xkp, Ukp, Xk = X[:,k-1], U[:,k-1], X[:,k]
  fpkp, Apkp, Bpkp, Xpkp, Upkp = get_f(k-1, model), get_A(k-1, model), get_B(k-1, model), Xp[:,k-1], Up[:,k-1]

  # Just Trapezoidal rule
  Uk = U[:,k]
  fpk, Apk, Bpk, Xpk, Upk = get_f(k, model), get_A(k, model), get_B(k, model), Xp[:,k], Up[:,k]
  return (Xkp-Xk) + 1/2*dtp.*(fpkp + Apkp*(Xkp-Xpkp) + Bpkp*(Ukp-Upkp) +
                              fpk  + Apk *(Xk-Xpk)   + Bpk *(Uk-Upk))

  # return A_dyn_discrete(X[:,k],dtp,robot,model)*X[:,k] + B_dyn_discrete(X[:,k],dtp,robot,model)*U[:,k] - X[:,k+1]
end

# Get current dynamics structures for a time step
get_f(k::Int, model::FreeflyerSE2) = model.f[k]
get_A(k::Int, model::FreeflyerSE2) = model.A
get_B(k::Int, model::FreeflyerSE2) = model.B

# Generate full dynamics structures for a time step
function f_dyn(x::Vector, u::Vector, robot::Robot, model::FreeflyerSE2)
  x_dim = model.x_dim
  f = zeros(x_dim)
  update_f!(f, x, u, robot, model)
  return f
end

function update_f!(f, x::Vector, u::Vector, robot::Robot, model::FreeflyerSE2)
  F, M = u[1:2], u[3]
  f[1:3] = x[4:6]
  f[4:5] = 1/robot.mass_ff*F
  f[6] = robot.J_ff_inv*M
end

function A_dyn(x::Vector, robot::Robot, model::FreeflyerSE2)
  kron([0 1; 0 0], Eye(3))
end

function B_dyn(x::Vector, robot::Robot, model::FreeflyerSE2)
  B = zeros(6,3)
  B[1:6,1:3] = kron([0;1], Eye(3))
  B[4:5,:] = 1/robot.mass_ff * B[4:5,:]
  B[6,:] = robot.J_ff_inv * B[6,:]
  return B
end

# Generate full discrete update version of dynamics matrices for a time step
# TODO(ambyld): Rename these? Initialize these once?
function A_dyn_discrete(x, dt, robot::Robot, model::FreeflyerSE2)
  kron([1 dt; 0 1], Eye(3))
end

function B_dyn_discrete(x, dt, robot::Robot, model::FreeflyerSE2)
  Jzz_inv = robot.J_ff_inv
  B = [0.5*dt^2*Eye(3); dt*Eye(3)]
  B[1:2,:] *= 1/robot.mass_ff  # translational force -> acceleration
  B[4:5,:] *= 1/robot.mass_ff
  B[3,:] *= Jzz_inv      # rotational moment -> acceleration
  B[6,:] *= Jzz_inv
  return B
end

## Convex state inequality constraints
function csi_translational_velocity_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}, k::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)
  return sum(X[3+j,k]^2 for j = 1:2)  - robot.hard_limit_vel^2
end

function csi_angular_velocity_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}, k::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)
  return X[6,k]^2 - robot.hard_limit_ω^2
end

## Convex control inequality constraints
function cci_translational_accel_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}, k::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)
  return 1/robot.mass_ff^2*sum(U[j,k]^2 for j = 1:2)  - robot.hard_limit_accel^2
end

function cci_angular_accel_bound(traj, traj_prev::Trajectory, SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}, k::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)
  Jzz_inv = robot.J_ff_inv
  return (Jzz_inv*U[3,k])^2 - robot.hard_limit_α^2
end

## Nonconvex state inequality constraints
function ncsi_body_obstacle_avoidance_constraints(traj, traj_prev::Trajectory, SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}, k::Int, i::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)

  rb_idx, env_idx = 1, i
  env_ = WS.btenvironment_keepout

  clearance = model.clearance

  r = [X[1:2,k];0]
  dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r,env_idx)

  # See Eq. 12b in "Convex optimization for proximity maneuvering of a spacecraft with a robotic manipulator"
  return clearance - dist
end

## Nonconvex state inequality constraints (convexified)
function ncsi_body_obstacle_avoidance_constraints_convexified(traj, traj_prev::Trajectory, SCPP::SCPProblem, k::Int, i::Int)
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)

  env_ = WS.btenvironment_keepout
  env_idx = i
  clearance = model.clearance

  # Base
  rb_idx = 1
  r0 = get_workspace_location(traj_prev, SCPP, k)
  dist, xbody, xobs = BulletCollision.distance(env_, rb_idx, r0, env_idx)
  r = get_workspace_location(traj, SCPP, k)

  if dist < SCPP.param.obstacle_toggle_distance
    # See Eq. 12b in "Convex optimization for proximity maneuvering of a spacecraft with a robotic manipulator"
    nhat = dist > 0 ?
      (xbody-xobs)./norm(xbody-xobs) :
      (xobs-xbody)./norm(xobs-xbody)

    return clearance - (dist + nhat'*(r-r0))
  else
    return 0.
  end
end

function ncsi_arm_obstacle_avoidance_constraints_convexified(traj, traj_prev::Trajectory, SCPP::SCPProblem, k::Int, i::Int)
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)

  env_ = WS.btenvironment_keepout
  env_idx = i
  clearance = model.clearance

  r0 = get_workspace_location(traj_prev, SCPP, k)

  # Arm
  rb_idx = 2
  xb = robot.xb
  thetap = Xp[3,k]
  R01p = [cos(thetap) -sin(thetap) 0.;
          sin(thetap)  cos(thetap) 0.;
          0. 0. 1.]
  r1 = R01p*xb + r0
  dist, xbody, xobs = BulletCollision.distance(env_, rb_idx, r1, env_idx)

  # return 0.   # toggle off

  if dist < SCPP.param.obstacle_toggle_distance
    nhat = dist > 0 ?
      (xbody-xobs)./norm(xbody-xobs) :
      (xobs-xbody)./norm(xobs-xbody)

    return clearance - (dist + nhat'*robot.Jcollision[k]*(X[1:3,k]-r0))
  else
    return 0.
  end
end

## State trust region inequality constraints
function stri_state_trust_region(traj, traj_prev::Trajectory, SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}, k::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)
  return sum((X[j,k]-Xp[j,k])^2 for j = 1:x_dim)
end

## Convex state inequality constraints
function ctri_control_trust_region(traj, traj_prev::Trajectory, SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}, k::Int) where {T,E}
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)
  return sum((U[j,k]-Up[j,k])^2 for j = 1:u_dim)
end

function get_workspace_location(traj, SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}, k::Int, i::Int=0) where {T,E}
  return [traj.X[1:2,k]; 0]
end

function SCPConstraints(SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}) where {T,E}
  model = SCPP.PD.model
  x_dim, u_dim, N = model.x_dim, model.u_dim, SCPP.N
  WS = SCPP.WS

  SCPC = SCPConstraints()

  ## Dynamics constraints
  add_constraint_category!(SCPC.dynamics, dynamics_constraints, :array, 2:N)

  ## State init constraints
  add_constraint_category!(SCPC.state_init_eq, sie_init_constraints, :scalar, 1, 1:x_dim)

  ## State boundary condition constraints
  goal_set, tf_guess = SCPP.PD.goal_set, SCPP.tf_guess
  for goal in values(inclusive(goal_set.goals, searchsortedfirst(goal_set.goals, tf_guess), searchsortedlast(goal_set.goals, tf_guess)))
    if typeof(goal.params) == PointGoal
      add_constraint_category!(SCPC.convex_state_boundary_condition_eq, csbce_goal_constraints, goal, :array)
    else
      add_constraint_category!(SCPC.convex_state_boundary_condition_ineq, csbci_goal_constraints, goal, :array)
    end
  end

  ## Convex state inequality constraints
  add_constraint_category!(SCPC.convex_state_ineq, csi_translational_velocity_bound, :scalar, 1:N)
  add_constraint_category!(SCPC.convex_state_ineq, csi_angular_velocity_bound, :scalar, 1:N)

  ## Nonconvex state equality constraints (convexified)
  nothing

  ## Nonconvex state inequality constraints
  env_ = WS.btenvironment_keepout
  N_obs = length(WS.btenvironment_keepout.convex_env_components)
  add_constraint_category!(SCPC.nonconvex_state_ineq, ncsi_body_obstacle_avoidance_constraints, :scalar, 1:N, 1:N_obs)

  ## Nonconvex state inequality constraints (convexified)
  add_constraint_category!(SCPC.nonconvex_state_convexified_ineq, ncsi_body_obstacle_avoidance_constraints_convexified, :scalar, 1:N, 1:N_obs)

  ## Convex control equality constraints
  nothing

  ## Convex control inequality constraints
  add_constraint_category!(SCPC.convex_control_ineq, cci_translational_accel_bound, :scalar, 1:N-1)
  add_constraint_category!(SCPC.convex_control_ineq, cci_angular_accel_bound, :scalar, 1:N-1)

  ## State trust region ineqality constraints
  add_constraint_category!(SCPC.state_trust_region_ineq, stri_state_trust_region, :scalar, 1:N)

  ## Constrol trust region inequality constraints
  add_constraint_category!(SCPC.control_trust_region_ineq, ctri_control_trust_region, :scalar, 1:N-1)

  return SCPC
end

function trust_region_ratio_gusto(traj, traj_prev::Trajectory, SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)
  fp, Ap = model.f, model.A
  num,den = 0, 0
  env_ = WS.btenvironment_keepout

  for k in 1:N-1
    linearized = fp[k] + Ap*(X[:,k]-Xp[:,k])
    num += norm(f_dyn(X[:,k],U[:,k],robot,model) - linearized)
    den += norm(linearized)
  end

  clearance = model.clearance

  # TODO(ambyld): Update for arm
  for k in 1:N
    r0 = get_workspace_location(traj_prev, SCPP, k)
    r = get_workspace_location(traj, SCPP, k)
    for (rb_idx, body_point) in enumerate(env_.convex_robot_components)
      for (env_idx,convex_env_component) in enumerate(env_.convex_env_components)
        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r0,env_idx)
        nhat = dist > 0 ?
          (xbody-xobs)./norm(xbody-xobs) :
          (xobs-xbody)./norm(xobs-xbody)
        linearized = clearance - (dist + nhat'*(r-r0))

        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r,env_idx)

        num += abs((clearance-dist) - linearized)
        den += abs(linearized)
      end
    end
  end
  return num/den
end

function trust_region_ratio_trajopt(traj, traj_prev::Trajectory, SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)
  fp = model.f
  num, den = 0, 0
  env_ = WS.btenvironment_keepout
  dt = traj.dt

  for k in 1:N-1
    ϕ_old = norm(fp[k] - (Xp[:,k]-Xp[:,k])/dtp, 1)
    ϕ_new = norm(f_dyn(X[:,k],U[:,k],robot,model) - (X[:,k]-X[:,k])/dt, 1)
    ϕ_hat_new = norm(dynamics_constraints(traj,traj_prev,SCPP,k,0), 1)
    num += (ϕ_old-ϕ_new)
    den += (ϕ_old-ϕ_hat_new)
  end

  clearance = model.clearance

  for k in 1:N
    r0 = get_workspace_location(traj_prev, SCPP, k)
    r = get_workspace_location(traj, SCPP, k)
    for (rb_idx,body_point) in enumerate(env_.convex_robot_components)
      for (env_idx,convex_env_component) in enumerate(env_.convex_env_components)
        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r0,env_idx)
        ϕ_old = clearance-dist

        dist,xbody,xobs = BulletCollision.distance(env_,rb_idx,r,env_idx)
        nhat = dist > 0 ?
          (xbody-xobs)./norm(xbody-xobs) :
          (xobs-xbody)./norm(xobs-xbody)
        ϕ_new = clearance-dist
        ϕ_hat_new = clearance - (dist + nhat'*(r-r0))

        num += (ϕ_old-ϕ_new)
        den += (ϕ_old-ϕ_hat_new)
      end
    end
  end

  return num/den
end

function trust_region_ratio_mao(traj, traj_prev::Trajectory, SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}) where {T,E}
  # Where i is the state index, and k is the timestep index
  X,U,Tf,Xp,Up,Tfp,dtp,robot,model,WS,x_init,goal_set,x_dim,u_dim,N,dh = @scp_shortcut_freeflyerSE2(traj, traj_prev, SCPP)
  num,den = 0, 0

  cost_true_prev = cost_true(traj_prev, traj_prev, SCPP)
  cost_true_new = cost_true(traj, traj, SCPP)
  for k in 1:N-1
    cost_true_prev += norm(dynamics_constraints(traj_prev, traj_prev, SCPP, k, 0),Inf)
    cost_true_new += norm(dynamics_constraints(traj, traj, SCPP, k, 0),Inf)
  end
  cost_linearized_new = cost_true(traj, traj, SCPP)
  return (cost_true_prev-cost_true_new)/(cost_true_prev-cost_linearized_new)
end

function get_dual_jump(SCPC::SCPConstraints, SCPP::SCPProblem{Freeflyer{T}, FreeflyerSE2, E}) where {T,E}
  init_constraint_category = SCPC.state_init_eq[:sie_init_constraints][1]
  -JuMP.dual.([init_constraint_category.con_reference[(i,)] for i = init_constraint_category.ind_other[1]])
end

##################
# Dynamics
##################
function A_SE2(rb::Freeflyer)
  kron([0 1; 0 0], Eye(3))
end

function B_SE2(rb::Freeflyer)
  B = zeros(6,3)
  B[1:6,1:3] = kron([0;1], Eye(3))
  B[4:5,:] = 1/rb.mass_ff * B[4:5,:]
  B[6,:] = rb.J_ff_inv * B[6,:]
  return B
end

function Ak_SE2(rb::Freeflyer, dt)
  kron([1 dt; 0 1], Eye(3))
end

function Bk_SE2(rb::Freeflyer, dt)
  Jzz_inv = robot.J_ff_inv
  Bk = [0.5*dt^2*Eye(3); dt*Eye(3)]
  Bk[1:2,:] *= 1/rb.mass_ff  # translational force -> acceleration
  Bk[4:5,:] *= 1/rb.mass_ff
  Bk[3,:] *= Jzz_inv      # rotational moment -> acceleration
  Bk[6,:] *= Jzz_inv
  return Bk
end

# continuous time dynamics
function f_SE2(X::Vector, U::Vector, rb::Freeflyer)
  F,M = U[1:2], U[3]
  Jzz_inv = robot.J_ff_inv
  xdot = zeros(length(X))
  xdot[1:3] = xdot[4:6]
  xdot[4:5] = 1/rb.mass_ff*F
  xdot[6] = Jzz_inv*M
  return xdot
end

max_stars_repo_name: StanfordASL/GuSTO

max_issues_repo_name: StanfordASL/GuSTO

max_forks_repo_name: StanfordASL/GuSTO


text:
# include the "exported" functions of the module
using Printf
using LinearAlgebra
using Statistics
using Random

include("byrow.jl")
include("dare.jl")
include("dlyap.jl")
include("simula01.jl")
include("ACQR_simula.jl")
include("ACQR_kfilter.jl")
include("ACQR_kfilter_s.jl")
include("ACQR_ksmoother.jl")
include("ACQR_ksmoother_s.jl")
include("ACQR_em_starting.jl")
include("ACQR_em.jl")
include("ACQR_em_s.jl")
include("ACQRS_sub.jl")
include("ACQR_em_run.jl")
include("ACQR_em_s_run.jl")


text:
using Documenter, PkgName

makedocs(
    sitename = "PkgName Documentation",
    modules = [PkgName]
)

max_stars_repo_name: ryapric/boilerplate

max_issues_repo_name: ryapric/boilerplate

max_forks_repo_name: ryapric/boilerplate


text:
# vectors

a = 1.2v1 + 1.89v3 + 1.2v4

# Homogeneous multivectors

A = 1.3v1 + 2.7v12
B = 0.7v3 + 1.1v123 + 3.05v
C = 0.7v4 + 0.1v2 + 1.59v23 + 1.1v124 + 3.05v1234

A₀ = 0.9v
A₁ = 1.1v1
A₂ = 1.96v23
A₃ = 1.2v123
A₄ = 1.0v1234

B₀ = 2.3v
B₁ = 0.8v2
B₂ = -1.8v13
B₃ = 2.5v234
B₄ = 1.7v1234

C₀ = 0.3v
C₁ = 0.4v4
C₂ = 0.92v24
C₃ = 1.0v134
C₄ = 1.1v1234

As = (A₀, A₁, A₂, A₃, A₄)
Bs = (B₀, B₁, B₂, B₃, B₄)
Cs = (C₀, C₁, C₂, C₃, C₄)

Zs = tuple(As..., Bs..., Cs...)

# Non-homogeneous multivectors

SA = sum(As)
SB = sum(Bs)
SC = sum(Cs)

NHB₀ = 2.3v
NHB₁ = 0.3v1 + 0.8v2
NHB₂ = 0.47v12 + 0.92v23 - 1.8v13
NHB₃ = 1.0v123 + 2.5v234
NHB₄ = 1.7v1234

NHC₀ = 0.3v
NHC₁ = 0.4v4 + 0.3v2
NHC₂ = 0.47v13 + 0.92v24 - 1.8v14
NHC₃ = 1.0v134 + 2.5v124
NHC₄ = 1.1v1234

NHBs = (NHB₀, NHB₁, NHB₂, NHB₃, NHB₄)
NHCs = (NHC₀, NHC₁, NHC₂, NHC₃, NHC₄)

NHZs = tuple(NHBs..., NHCs...)


text:
# --- PRIVATE METHODS BELOW HERE ----------------------------------------------------------------- #
# --- PRIVATE METHODS ABOVE HERE ----------------------------------------------------------------- #

# --- PUBLIC METHODS BELOW HERE ------------------------------------------------------------------ #
"""
    𝝙(contractSet::Set{VLAbstractAsset}, latticeModel::VLBinomialLattice, underlyingPrice::Float64;
        decisionLogic::Function=_american_decision_logic) -> VLResult

Magical description goes here!
"""
function 𝝙(contractSet::Set{VLAbstractAsset}, latticeModel::VLBinomialLattice, underlyingPrice::Float64;
    decisionLogic::Function=_american_decision_logic)::VLResult

    try

        # initialize -
        baseUnderlyingPrice = underlyingPrice
        𝛅 = 0.01 * baseUnderlyingPrice
        downPrice = (baseUnderlyingPrice - 𝛅)
        upPrice = (baseUnderlyingPrice + 𝛅)

        # compute the price for the base underlying value -
        up_price_tree = binomial_price(contractSet, latticeModel, upPrice; decisionLogic=decisionLogic) |> check
        down_price_tree = binomial_price(contractSet, latticeModel, downPrice; decisionLogic=decisionLogic) |> check

        # compute delta -
        delta_value = (up_price_tree[1,3] - down_price_tree[1,3]) / (2 * 𝛅)

        # return -
        return VLResult(delta_value)
    catch error
        return VLResult(error)
    end
end

"""
    ϴ(contractSet::Set{VLAbstractAsset}, latticeModel::VLBinomialLattice, underlyingPrice::Float64, eps::Float64;
        decisionLogic::Function=_american_decision_logic) -> VLResult

Magical description goes here!
"""
function ϴ(contractSet::Set{VLAbstractAsset}, latticeModel::VLBinomialLattice, underlyingPrice::Float64, eps::Float64;
    decisionLogic::Function=_american_decision_logic)::VLResult

    try

        baseUnderlyingPrice = underlyingPrice
        number_of_levels = latticeModel.numberOfLevels
        current_delta_T = latticeModel.𝝙t
        𝛅 = eps * current_delta_T

        # compute up and down price -
        theta_value = 0.0
        base_price_tree = binomial_price(contractSet, latticeModel, baseUnderlyingPrice; decisionLogic=decisionLogic) |> check
        latticeModel.𝝙t = 𝛅
        down_price_tree = binomial_price(contractSet, latticeModel, baseUnderlyingPrice; decisionLogic=decisionLogic) |> check

        # compute theta -
        theta_value = (down_price_tree[1,3] - base_price_tree[1,3])

        # return -
        return VLResult(theta_value)
    catch error
        return VLResult(error)
    end
end
# --- PUBLIC METHODS ABOVE HERE ------------------------------------------------------------------ #


text:
"""
Background for the example in the README
"""

loc_coords=InterpolationFunction(
    :(a0+a1*xi+a2*eta), # expression
    [:xi, :eta], # symbols for coordinates
    [:a0, :a1, :a2], # symbols for interpolation coefficients
    [ # degrees of freedom
        [0, 0] # values at corners (see explanation below)
    ];
    maximum_derivative=1 # we'll only really need first order derivatives
)

abs_coords=InterpolationFunction(
    :(a0+a1*x+a2*y),
    [:x, :y],
    [:a0, :a1, :a2],
    [
        [0, 0]
    ];
    maximum_derivative=1
)

dmn=Domain(
    loc_coords, # local coordinates InterpolationFunction struct
    [ # frontier of the domain in each coordinate of loc_coords
        ( # xi
            0.0, 1.0
        ),
        ( # eta
            0.0, :(1.0-xi) # linear triangular element
        )
    ],
    [ # points at which the interpolation from local to absolute coordinates
        [0.0, 0.0], # is imposed (triangle corners)
        [1.0, 0.0],
        [0.0, 1.0] # (xi, eta)
    ]
)

fin=Finel(
    [ # corners
        [1.0, 1.0],
        [2.0, 1.0],
        [2.0, 2.0]
    ],
    dmn # domain of integration
)

resd=Residual(
    [
        ( # first variable: weight function
            :w, abs_coords, [] # no derivatives calculated
        ),
        ( # second variable: u
            :u, abs_coords, [
                (:du!dx, [1, 0]) # first order derivative in x axis calculated
            ]
        )
    ],
    :(w*(du!dx-v)), # expression for the residual
    3; # Gauss-Legendre quadrature order
    extra_args=[:v] # v, a constant
)

resfunc=get_volume_residual_function(resd)

ws=[0.0, 1.0, 1.0] # format compatible with abs_coords InterpolationFunction,
us=[0.5, 1.0, 1.0] # as detailed before
v=-1.0

R=resfunc([ws, us], fin, v) # and any other extra_args symbols

# Test it!! :)
@assert isapprox(R, 0.5)

surface_resd=Residual(
    [
        ( # first variable: weight function
            :w, abs_coords, [] # no derivatives calculated
        ),
        ( # second variable: u
            :u, abs_coords, [
                (:du!dx, [1, 0]) # first order derivative in x axis calculated
            ]
        )
    ],
    :(w.*[du!dx, v]), # expression for the residual
    3; # Gauss-Legendre quadrature order, now at domain faces
    extra_args=[:v] # v, a constant
)

surface_resfunc=get_surface_residual_function(surface_resd)

R=surface_resfunc([ws, us], fin, v)

# test it!! :)
@assert isapprox(R, 0.25)


text:
import SimulateGalaxy
using Base.Test
include("galaxy_generation.jl")

#Test NFW Definition
nfwMR = [2.0, -5.3, 2.5, 0.16, 1.5, -9.0, 6.9, 0.086, 21.0, 1.5]
nfwMP = [2.4, -7.9, 1.1, 0.17, 3.0, 0.0, 8.2, 0.086, 21.0, 1.5]
nfwpMR = SimulateGalaxy.NFWParameters(nfwMR...)
nfwpMP = SimulateGalaxy.NFWParameters(nfwMP...)

@testset "Computed NFW Parameters" begin
    @testset "MR" begin
        @test nfwpMR.rs == nfwMR[10] / 2.16
        @test nfwpMR.Φs == (nfwMR[9] / 0.465) ^ 2
        @test nfwpMR.xlim == nfwMR[5] / nfwMR[10] * 2.16
        @test nfwpMR.Φlim == (nfwMR[9] / 0.465) ^ 2 * GravPotential(nfwMR[5] / nfwMR[10] * 2.16)
        @test nfwpMR.adjEc == nfwMR[4] * (nfwMR[9] / 0.465) ^ 2
        @test_approx_eq nfwpMR.adjJb nfwMR[8] * nfwMR[10] / 2.16 * sqrt((nfwMR[9] / 0.465) ^ 2)
    end
    @testset "MP" begin
        @test nfwpMP.rs == nfwMP[10] / 2.16
        @test nfwpMP.Φs == (nfwMP[9] / 0.465) ^ 2
        @test nfwpMP.xlim == nfwMP[5] / nfwMP[10] * 2.16
        @test nfwpMP.Φlim == (nfwMP[9] / 0.465) ^ 2 * GravPotential(nfwMP[5] / nfwMP[10] * 2.16)
        @test nfwpMP.adjEc == nfwMP[4] * (nfwMP[9] / 0.465) ^ 2
        @test_approx_eq nfwpMP.adjJb nfwMP[8] * nfwMP[10] / 2.16 * sqrt((nfwMP[9] / 0.465) ^ 2)
    end
end

#Test SFW Definition
sfwMR = [2.0, -5.3, 2.5, 0.16, 1.5, -9.0, 6.9, 0.086, 4229.2, .69444, 1.2, 3.05, 1.1]
sfwMP = [2.4, -7.9, 1.1, 0.17, 3.0, 0.0, 8.2, 0.0, 4229.2, .69444, 1.2, 3.05, 1.1]
sfwpMR = SimulateGalaxy.SFWParameters(sfwMR...)
sfwpMP = SimulateGalaxy.SFWParameters(sfwMP...)

@testset "Computed NFW Parameters" begin
    @testset "MR" begin
        @test sfwpMR.Φs == sfwMR[9] * sfwMR[10]^2
        @test sfwpMR.xlim == sfwMR[5] / sfwMR[10]
        @test sfwpMR.Φ0 == genphi0(10.0^-8, sfwMR[11], sfwMR[12], sfwMR[13], sfwMR[9], sfwMR[10])
        @test sfwpMR.Φlim == genphi(sfwMR[5] / sfwMR[10], sfwMR[11], sfwMR[12], sfwMR[13], sfwMR[9], sfwMR[10])
        @test sfwpMR.adjEc == sfwMR[4] * sfwMR[9] * sfwMR[10]^2
        @test_approx_eq sfwpMR.adjJb sfwMR[8] * sqrt(sfwMR[9]) * sfwMR[10]^2
    end
    @testset "MP" begin
        @test sfwpMP.Φs == sfwMP[9] * sfwMP[10]^2
        @test sfwpMP.xlim == sfwMP[5] / sfwMP[10]
        @test sfwpMP.Φ0 == genphi0(10.0^-8, sfwMP[11], sfwMP[12], sfwMP[13], sfwMP[9], sfwMP[10])
        @test sfwpMP.Φlim == genphi(sfwMP[5] / sfwMP[10], sfwMP[11], sfwMP[12], sfwMP[13], sfwMP[9], sfwMP[10])
        @test sfwpMP.adjEc == sfwMP[4] * sfwMP[9] * sfwMP[10]^2
        @test_approx_eq sfwpMP.adjJb sfwMP[8] * sqrt(sfwMP[9]) * sfwMP[10]^2
    end
end

#Test gravitational_potential
n = 5
xMR = rand(n) * nfwpMR.rlim / nfwpMR.rs
xMP = rand(n) * nfwpMP.rlim / nfwpMP.rs
@testset "gravitational_potential Tests" begin
    @testset "MR" begin
        for x in xMR
            @test SimulateGalaxy.gravitational_potential(x) == GravPotential(x)
            @test SimulateGalaxy.gravitational_potential(x, nfwpMR) == GravPotential(x) * (nfwMR[9] / 0.465) ^ 2
            @test SimulateGalaxy.gravitational_potential(x, sfwpMR) == genphi(x, sfwMR[11], sfwMR[12], sfwMR[13], sfwMR[9], sfwMR[10])
        end
    end
    @testset "MP" begin
        for x in xMP
            @test SimulateGalaxy.gravitational_potential(x) == GravPotential(x)
            @test SimulateGalaxy.gravitational_potential(x, nfwpMP) == GravPotential(x) * (nfwMP[9] / 0.465) ^ 2
            @test SimulateGalaxy.gravitational_potential(x, sfwpMP) == genphi(x, sfwMP[11], sfwMP[12], sfwMP[13], sfwMP[9], sfwMP[10])
        end
    end
end

#Test escape_velocity
@testset "escape_velocity Tests" begin
    @testset "MR" begin
        for x in xMR
            @test_approx_eq SimulateGalaxy.escape_velocity(x, nfwpMR) EscapeVelocity(x, nfwMR[5], nfwMR[10] / 2.16, (nfwMR[9] / 0.465) ^ 2)
            @test SimulateGalaxy.escape_velocity(x, sfwpMR) == escape_velocity(x, sfwMR[5], sfwMR[9], sfwMR[10], sfwMR[11], sfwMR[12], sfwMR[13])
        end
    end
    @testset "MP" begin
        for x in xMP
            @test_approx_eq SimulateGalaxy.escape_velocity(x, nfwpMP) EscapeVelocity(x, nfwMP[5], nfwMP[10] / 2.16, (nfwMP[9] / 0.465) ^ 2)
            @test SimulateGalaxy.escape_velocity(x, sfwpMP) == escape_velocity(x, sfwMP[5], sfwMP[9], sfwMP[10], sfwMP[11], sfwMP[12], sfwMP[13])
        end
    end
end

#Test NFW profile_density
vrMR = rand(n) * nfwpMR.vmax
vtMR = rand(n) * nfwpMR.vmax
vrMP = rand(n) * nfwpMP.vmax
vtMP = rand(n) * nfwpMP.vmax

@testset "NFW profile_density tests" begin
    @testset "MR" begin
        for (x, vr, vt) in zip(xMR, vrMR, vtMR)
            @test_approx_eq SimulateGalaxy.profile_density(x, vr, vt, nfwpMR) NFWDensity(x, vr, vt, nfwMR)
        end
    end
    @testset "MP" begin
        for (x, vr, vt) in zip(xMP, vrMP, vtMP)
            @test SimulateGalaxy.profile_density(x, vr, vt, nfwpMP) == NFWDensity(x, vr, vt, nfwMP)
        end
    end
end

#Test SFW profile_density
@testset "SFW profile_density tests" begin
    @testset "MR" begin
        for (x, vr, vt) in zip(xMR, vrMR, vtMR)
            @test_approx_eq SimulateGalaxy.profile_density(x, vr, vt, sfwpMR) SFW_density(x, vr, vt, sfwMR)
        end
    end
    @testset "MP" begin
        for (x, vr, vt) in zip(xMP, vrMP, vtMP)
            @test SimulateGalaxy.profile_density(x, vr, vt, sfwpMP) == SFW_density(x, vr, vt, sfwMP)
        end
    end
end

#Test simulate_galaxy NFW
@time SimulateGalaxy.simulate_galaxy(nfwpMR, 400, rate = true)
@time SimulateGalaxy.simulate_galaxy(nfwpMP, 400, rate = true)

#Test simulate_galaxy SFW
@time SimulateGalaxy.simulate_galaxy(sfwpMR, 400, rate = true)
@time SimulateGalaxy.simulate_galaxy(sfwpMP, 400, rate = true)

#Test simulate_galaxy NFW with metallicity parameters
@time SimulateGalaxy.simulate_galaxy(nfwpMR, -2.0,  .003, 400, rate = true)
@time SimulateGalaxy.simulate_galaxy(nfwpMP, -2.0,  .003, 400, rate = false)

#Test simulate_galaxy SFW with metallicity parameters
@time SimulateGalaxy.simulate_galaxy(sfwpMR, -2.0,  .003, 400, rate = true)
@time SimulateGalaxy.simulate_galaxy(sfwpMP, -2.0,  .003, 400, rate = false)

#=
p = nfwpMR
samplesize = 10
function nfwopt(x::Array{Float64, 1})
    return -(x[1] * x[1] * x[3] * SimulateGalaxy.profile_density(x[1], x[2], x[3], p))
end

minProb = Optim.optimize(nfwopt, [0.1, 1.0, 1.0], method = Optim.NelderMead())
fmax = max(0.0,  -1.1 * minProb.minimum)

vmax0 = SimulateGalaxy.escape_velocity(10.0^-8, p)
sampledr = Array{Float64}(samplesize)
sampledvr = Array{Float64}(samplesize)
sampledvt = Array{Float64}(samplesize)
if rate
    tested = 0
end
rate = false
accepted = 0
while accepted < samplesize
    x = (p.rlim  / p.rs) * rand()
    vr = vmax0 * rand()
    vt = vmax0 * rand()
    u = rand()
    sampleDensity = x * x * vt * SimulateGalaxy.profile_density(x, vr, vt, p)
    if rate
        tested += 1
    end
    if sampleDensity >= (u * fmax)
        accepted += 1
        sampledr[accepted] = x * p.rs
        sampledvr[accepted] =  vr
        sampledvt[accepted] = vt
    end
end
SimulateGalaxy.SphericalGalaxy(sampledr, sampledvr, sampledvt)
=#


text:
using StatsBase

input = joinpath(@__DIR__, "input")
lines = readlines(input)

function endpoints(line)
    components = map(n -> parse(Int, n), split(line, r",| -> "))
    return ((components[1], components[2]), (components[3], components[4]))
end

function orthogonals(coords)
    map(c -> begin
        (start, finish) = c
        if start[1] == finish[1]
            [(start[1], i) for i in minimum([start[2], finish[2]]):maximum([start[2], finish[2]])]
        elseif start[2] == finish[2]
            [(i, start[2]) for i in minimum([start[1], finish[1]]):maximum([start[1], finish[1]])]
        else
            nothing
        end
    end, coords)
end

function diagonals(coords)
    map(c -> begin
        (start, finish) = c
        if start[1] < finish[1] && start[2] < finish[2]
            return collect(zip(start[1]:finish[1], start[2]:finish[2]))
        elseif start[1] > finish[1] && start[2] > finish[2]
            return collect(zip(start[1]:-1:finish[1], start[2]:-1:finish[2]))
        elseif start[1] < finish[1] && start[2] > finish[2]
            return collect(zip(start[1]:finish[1], start[2]:-1:finish[2]))
        elseif start[1] > finish[1] && start[2] < finish[2]
            return collect(zip(start[1]:-1:finish[1], start[2]:finish[2]))
        else
            nothing
        end
    end, coords)
end

coords = map(endpoints, lines)
orthos = filter(!isnothing, orthogonals(coords))
dias = filter(!isnothing, diagonals(coords))

p1 = length(filter(p -> last(p) > 1, countmap(Iterators.flatten(orthos))))
p2 = length(filter(p -> last(p) > 1, countmap(Iterators.flatten(append!(orthos, dias)))))

println("-----------------------------------------------------------------------")
println("hydrothermal venture -- part one :: $p1")
println("hydrothermal venture -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 6005)
@assert(p2 == 23864)


text:
"""
mrf : code associated with dealing with Markov Random Fields, calculating pseudolikelihoods, reading inputs, etc.
		We try to accommodate both Gremlin and CCMPred outputs and we use a parameterization of w1/w2 matrices
		that is easier to code with, so some code to glue the formats together.
"""

module mrf

using JLD, LinearAlgebra

include("types.jl")
include("math.jl")
include("bio_seq.jl")
include("utils.jl")

function debug(x)
	println("DEBUG: $x")
end

function err(x)
	println("ERR: $x")
end

function get_gremlin_consensus(file_name, prot_name) #NOTE: extra dependency (optional_consensus).
	#I maintain some protein translations in a file because they align better with hmmalign or are wild-type versions of proteins.
	#The wild-type thing is relevant b/c you'd presumably be inserting most genes into broader wild-type contexts you'd like long-range interactions to be adapted towards.
	gen_con = true
	if "optional_consensus.fa" in readdir()
		gen_con = false
		reference_seqs = bio_seq.load_fasta("optional_consensus.fa")
		if prot_name in keys(reference_seqs)
			return reference_seqs[prot_name]
		else
			gen_con = true
		end
	end

	if gen_con
		#If I don't have a specific sequence in mind, just generate the consensus sequence (based on w1 matrix only) directly from MRF.
		lookup_nodes = Dict{Int64, Int64}() #based on number of weights, determine number of nodes in MRF.
		for i=1:1500 #assume we don't get proteins longer than this.
			lookup_nodes[((i*(i-1)*21*21)/2) + (i*20)] = i
		end
		w1, w2 = read_weights(file_name)
		nNodes = div(size(w1)[1], 21) #lookup_nodes[length(weights)]
		w1 = reshape(w1, (21, nNodes))'[1:end, 1:20] #this is sorta legacy code to get consensus.
		int_dict = Dict{Int64, String}(1 => "A", 2 => "R", 3 => "N", 4 => "D", 5 => "C", 6 => "Q", 7 => "E", 8 => "G", 9 => "H", 10 => "I", 11 => "L", 12 => "K", 13 => "M", 14 => "F", 15 => "P", 16 => "S", 17 => "T", 18 => "W", 19 => "Y", 20 => "V", 21 => "X")
		prot_con = ""
		for i in 1:nNodes
			prot_con *= int_dict[findmax(w1[i,1:end])[2]]
		end

		return prot_con
	end
end

#Split dense form of weights into w1/w2 variables.
function split_weights(w, nNodes, nStates, nEdges)
	w1 = reshape(w[1:(nStates - 1)*nNodes], nNodes, nStates-1)
	offset = (nStates - 1) * nNodes
	w2 = reshape(w[offset+1:offset+nEdges*nStates^2], nStates, nStates, nEdges)
	return w1, w2
end

function read_weights(jld_file)
	weights = load(replace(jld_file, ".csv"=>".jld")) #["weights"]
	#return weights
	w1 = weights["w1"]
	w2 = weights["w2"]
	return w1, w2
end

function read_csv_weights(csv_file) #Old version just read csv file, jld file is slightly more space efficient for storage.
	#For now reading weights from text file
	in_file = open(csv_file)
	line = readline(in_file)
	weights = Float32[]
	while !eof(in_file)
		push!(weights, parse(Float32, line))
		line = readline(in_file)
	end
	push!(weights, parse(Float32, line))
	close(in_file)
	return weights
end

function transform_weights(w1, w2, nNodes) #transforms strange vector of weights into more interpretable matrix of weights.
	new_w1 = hcat(w1, zeros(nNodes))
	new_w1 = reshape(new_w1', nNodes * 21, 1)
	new_w2 = zeros(Float32, nNodes * 21, nNodes * 21)
	edge_counter = 0
	for edge in collect(utils.ifilter(utils.two_lt, utils.big_vec(nNodes)))
		edge_counter += 1
		ei, ej = edge
		square_mat = w2[1:21, 1:21, edge_counter]
		for mi in 1:21
			for mj in 1:21
				new_w2[((ei - 1) * 21) + mj, ((ej - 1) * 21) + mi] = square_mat[mj, mi]
				new_w2[((ej - 1) * 21) + mi, ((ei - 1) * 21) + mj] = square_mat[mj, mi]
			end
		end
	end
	return new_w1, new_w2
end

function standardize_weights(w1, w2, nNodes)
	#This feels pretty deprecatable.
	new_w1, new_w2 = transform_weights(w1, w2, nNodes)
	return map(Float32, new_w1), map(Float32, new_w2)
end

function init_model(grem_file)
	lookup_nodes = Dict{Int64, Int64}()
	for i=1:1500 #assume we don't get proteins longer than this.
		lookup_nodes[((i*(i-1)*21*21)/2) + (i*20)] = i #Finally don't have to give numNodes. Works for even or odd.
	end
	w1, w2 = read_weights(grem_file)

	nNodes = div(size(w1)[1], 21)

	return types.GremModel(w1, w2, nNodes, 0) #not prot or base_logpot. That comes later.
end

function cpu_assess(prot_mat, w1, w2, nNodes, nProt) #removes cuda dependency at expense of speed.
	my_pv_w2 = prot_mat * w2
	my_ull = my_pv_w2 * prot_mat'
	my_pv_w1 = prot_mat * w1
	numerator = my_pv_w1 + diag(my_ull)
	partitions = reshape((w1' .+ my_pv_w2)', 21, nNodes * nProt)
	logZs = sum(reshape(math.mylogsumexp(partitions'), nNodes, nProt), dims=1)'
	all_scores = numerator - logZs
	return all_scores, my_ull, my_pv_w1, my_pv_w2
end

function psl(prot_seq, w1, w2, verbose = false)
	nNodes = length(prot_seq)
	prot_vec = bio_seq.convert_protein(prot_seq)
	pv_w2 = prot_vec * w2
	pv_w1 = prot_vec * w1
	ull = pv_w1 .+ (pv_w2 * prot_vec')
	ull = ull[1] * -1

	if verbose
		debug("The unnormalized LL is $ull")
	end

	new_idea = w1' + pv_w2
	new_idea = reshape(new_idea, 21, nNodes)
	new_logZ = math.mylogsumexp(new_idea')
	if verbose
		debug("\tThe logZ is $(new_logZ' * ones(nNodes, 1))")
	end
	psi_nll = ull .+ new_logZ' * ones(nNodes, 1)
	return psi_nll[1]
end

function basic_energy_calc(prot_seq, w1, w2)
	nNodes = length(prot_seq)
	prot_vec = bio_seq.convert_protein(prot_seq)
	energy = prot_vec * w1 .+ prot_vec * w2 * prot_vec'
	return energy[1] * -1
end

end


text:
#' @@name count
#'
#' @@description
#'
#' Count the number of missing values in every column of an
#' AbstractDataFrame.
#'
#' @@arg adf::AbstractDataFrame An AbstractDataFrame.
#' @@arg skip::Int The index of a column to skip.
#'
#' @@return missing::Vector{Int} The number of missing values in each column.
#' @@return success::Bool Did the operation complete successfully?
#'
#' @@examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "y", "z"])
#' m, s = count(df, 0)

function count(adf::AbstractDataFrame, skip::Int)
	return Int[], false
end

#' @@name sum
#'
#' @@description
#'
#' Sum the elements of each column of an AbstractDataFrame.
#'
#' @@arg adf::AbstractDataFrame An AbstractDataFrame.
#'
#' @@return sums::Vector{Float64} The sums of each column's entries.
#'
#' @@examples
#'
#' df = DataFrame(A = 1:3, B = ["x", "y", "z"])
#' ss = sum(df)

function sum(adf::AbstractDataFrame)
	return Float64[]
end


text:
module MethodOfParticularSolutions

using Nemo
using ArbTools
using Random
using StaticArrays
using CoordinateTransformations
using LinearAlgebra
using GenericSVD
using Optim
using OrderedCollections
import SpecialFunctions

using TimerOutputs
using Printf
using RecipesBase

export AbstractDomain,
    boundaries,
    area,
    boundary_parameterization,
    boundary_points,
    interior_points,
    SphericalTriangle,
    angles,
    vertex,
    center,
    greatcircleplane,
    greatcircle,
    LShape,
    Triangle,
    Polygon,
    TransformedDomain,
    IntersectedDomain,
    AbstractEigenfunction,
    coefficients,
    set_eigenfunction!,
    SphericalVertexEigenfunction,
    mu,
    SphericalInteriorEigenfunction,
    SphericalCombinedEigenfunction,
    basis_function,
    KrewerasEigenfunction,
    LShapeEigenfunction,
    StandaloneVertexEigenfunction,
    StandaloneInteriorEigenfunction,
    StandaloneLightningEigenfunction,
    LinkedEigenfunction,
    CombinedEigenfunction,
    sigma,
    sigma_coefficients,
    cartesian,
    spherical,
    mps!,
    iteratemps,
    certifyindices

include("arb.jl")
include("arb_series.jl")
include("utilities.jl")
include("special-functions.jl")

# Include types for domains and eigenfunctions
include("domain/domain_types.jl")
include("eigenfunction/eigenfunction_types.jl")

# Include methods for domains
include("domain/AbstractDomain.jl")
include("domain/SphericalTriangle.jl")
include("domain/LShape.jl")
include("domain/Triangle.jl")
include("domain/Polygon.jl")
include("domain/TransformedDomain.jl")
include("domain/IntersectedDomain.jl")
include("domain/plotting.jl")

# Include methods for eigenfunctions
include("eigenfunction/AbstractEigenfunction.jl")
include("eigenfunction/AbstractPlanarEigenfunction.jl")
include("eigenfunction/AbstractSphericalEigenfunction.jl")
include("eigenfunction/SphericalVertexEigenfunction.jl")
include("eigenfunction/SphericalInteriorEigenfunction.jl")
include("eigenfunction/SphericalCombinedEigenfunction.jl")
include("eigenfunction/KrewerasEigenfunction.jl")
include("eigenfunction/LShapeEigenfunction.jl")
include("eigenfunction/StandaloneVertexEigenfunction.jl")
include("eigenfunction/StandaloneInteriorEigenfunction.jl")
include("eigenfunction/StandaloneLightningEigenfunction.jl")
include("eigenfunction/LinkedEigenfunction.jl")
include("eigenfunction/CombinedEigenfunction.jl")
include("eigenfunction/plotting.jl")

include("sigma.jl")
include("enclose.jl")
include("mps.jl")
include("index.jl")
include("trace.jl")

include("examples/triangles.jl")
include("examples/denominators.jl")
include("examples/domains.jl")

end # module


text:
using ArrayInterfaceBlockBandedMatrices
using ArrayInterfaceCore
using BlockBandedMatrices
using Test

BB=BlockBandedMatrix(Ones(10,10),[1,2,3,4],[4,3,2,1],(1,0))
BB[Block(1,1)].=[1 2 3 4]
BB[Block(2,1)].=[5 6 7 8;9 10 11 12]
rowind,colind=ArrayInterfaceCore.findstructralnz(BB)
@test [BB[rowind[i],colind[i]] for i in 1:length(rowind)]==
    [1,5,9,2,6,10,3,7,11,4,8,12,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    1,1,1,1,1,1,1,1,1,1,1,1,1,1,
    1,1,1,1,1]
@test ArrayInterfaceCore.isstructured(typeof(BB))
@test ArrayInterfaceCore.has_sparsestruct(typeof(BB))
@test ArrayInterfaceCore.fast_matrix_colors(typeof(BB))

dense=collect(Ones(8,8))
for i in 1:8
    dense[:,i].=[1,2,3,4,5,6,7,8]
end
BBB=BandedBlockBandedMatrix(dense, [4, 4] ,[4, 4], (1, 1), (1, 1))
rowind,colind=ArrayInterfaceCore.findstructralnz(BBB)
@test [BBB[rowind[i],colind[i]] for i in 1:length(rowind)]==
    [1,2,3,1,2,3,4,2,3,4,5,6,7,5,6,7,8,6,7,8,
     1,2,3,1,2,3,4,2,3,4,5,6,7,5,6,7,8,6,7,8]
@test ArrayInterfaceCore.isstructured(typeof(BBB))
@test ArrayInterfaceCore.has_sparsestruct(typeof(BBB))
@test ArrayInterfaceCore.fast_matrix_colors(typeof(BBB))


text:
include("LatticeUtils.jl")
include("Crystal.jl")
include("Bond.jl")
include("SymmetryAnalysis.jl")
include("AllowedCouplings.jl")
include("Parsing.jl")
include("Printing.jl")


text:
import LinearAlgebra: lmul!, rmul!

function rmul!(A::AbstractVector, D::Diagonal)
    @assert !LinearAlgebra.has_offset_axes(A)
    A .= A .* transpose(D.diag)
    return A
end

function lmul!(D::Diagonal, B::AbstractVector)
    @assert !LinearAlgebra.has_offset_axes(B)
    B .= D.diag .* B
    return B
end


text:
using PowerMarketAdeq
using Test

@testset "Dummy test" begin
    @test PowerMarketAdeq.dummy_func(2,1) == 5
    @test PowerMarketAdeq.dummy_func(3,1) == 7
end


text:
#####################################################################################
# MPCCMCounters type ; based on NLPModels Counters:
# https://github.com/JuliaSmoothOptimizers/NLPModels.jl/blob/master/src/nlp_types.jl
#
# Problem both NLPModels and MPCC export neval_hess,...
#####################################################################################

"""
MPCCCounters

Initialization: `MPCCCounters()`
"""
mutable struct MPCCCounters

  neval_obj     :: Int  # Number of objective evaluations.
  neval_grad    :: Int  # Number of objective gradient evaluations.
  neval_cons    :: Int  # Number of constraint vector evaluations.
  neval_consG   :: Int  # Number of constraint vector evaluations of G(x).
  neval_consH   :: Int  # Number of constraint vector evaluations of H(x).
  neval_jcon    :: Int  # Number of individual constraint evaluations.
  neval_jgrad   :: Int  # Number of individual constraint gradient evaluations.
  neval_jac     :: Int  # Number of constraint Jacobian evaluations.
  neval_jacG    :: Int  # Number of constraint Jacobian evaluations of G(x).
  neval_jacH    :: Int  # Number of constraint Jacobian evaluations of H(x).
  neval_jprod   :: Int  # Number of Jacobian-vector products.
  neval_jGprod  :: Int  # Number of Jacobian-vector products.
  neval_jHprod  :: Int  # Number of Jacobian-vector products.
  neval_jtprod  :: Int  # Number of transposed Jacobian-vector products.
  neval_jGtprod :: Int  # Number of transposed Jacobian-vector products.
  neval_jHtprod :: Int  # Number of transposed Jacobian-vector products.
  neval_hess    :: Int  # Number of Lagrangian/objective Hessian evaluations.
  neval_hessG   :: Int  # Number of Lagrangian/objective Hessian evaluations.
  neval_hessH   :: Int  # Number of Lagrangian/objective Hessian evaluations.
  neval_hprod   :: Int  # Number of Lagrangian/objective Hessian-vector products.
  neval_hGprod  :: Int  # Number of Lagrangian/objective Hessian-vector products.
  neval_hHprod  :: Int  # Number of Lagrangian/objective Hessian-vector products.
  neval_jhprod  :: Int  # Number of individual constraint Hessian-vector products.

  function MPCCCounters()
    return new(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
  end
end

# simple default API for retrieving counters
for counter in fieldnames(MPCCCounters)
  @eval begin
    """
        $($counter)(nlp)
    Get the number of `$(split("$($counter)", "_")[2])` evaluations.
    """
    $counter(nlp :: AbstractMPCCModel) = nlp.counters.$counter
    export $counter
  end
end

"""
    increment!(nlp, s)
Increment counter `s` of problem `nlp`.
"""
function increment!(nlp :: AbstractMPCCModel, s :: Symbol)
  setfield!(nlp.counters, s, getfield(nlp.counters, s) + 1)
end

"""
    decrement!(nlp, s)
Increment counter `s` of problem `nlp`.
"""
function decrement!(nlp :: AbstractMPCCModel, s :: Symbol)
  if s in fieldnames(MPCCCounters)
    setfield!(nlp.counters, s, getfield(nlp.counters, s) - 1)
  end
end

"""
    sum_counters(counters)
Sum all counters of `counters`.
"""
sum_counters(c :: MPCCCounters) = sum(getfield(c, x) for x in fieldnames(MPCCCounters))

"""
    sum_counters(nlp)
Sum all counters of problem `nlp`.
"""
sum_counters(nlp :: AbstractMPCCModel) = sum_counters(nlp.counters)

"""
    reset!(counters)
Reset evaluation counters
"""
function reset!(counters :: MPCCCounters)
  for f in fieldnames(MPCCCounters)
    setfield!(counters, f, 0)
  end
  return counters
end

"""
    reset!(nlp)
Reset evaluation count in `nlp`
"""
function reset!(nlp :: AbstractMPCCModel)
  reset!(nlp.counters)
  return nlp
end


text:
module SCoBADomains

using Parameters
using StaticArrays
using Random
using Hungarian
using Distributions
using DataStructures
using ArgParse
using POMDPs
using Distances
using NearestNeighbors
using TOML

using Plots
using GeometryTypes

using SCoBA.SCoBASolver

# Conveyor belt types and simulator
export
    Arm,
    ArmProperties,
    ConveyorBelt,
    ArmAllocationServer,
    count_terminal_grabbed_boxes!,
    propagate_belt!,
    update_time_windows!,
    setup_mirror_gen!,
    mirror_gen_step!,
    add_box_to_belt!,
    propagate_arms!,
    setup_conveyor_scenario,
    parse_conveyor_commandline

# Conveyor belt algorithms - SCoBA and baseline stuff
export
    cumul_grasp_prob,
    scoba_conveyor!,
    expected_hungarian!,
    earliest_due_date!,
    ConvBeltSystemState,
    ConvBeltMDP,
    create_convbelt_mdp,
    get_current_belt_mdp_state,
    EDDPolicy,
    ConvBeltSingleArmMDP,
    create_convbelt_singlearm_mdp,
    get_singlearm_mdp_state,
    update_server_with_mdp_action!

# Routing types and simulator
export
    parse_city_params,
    Drone,
    DroneProperties,
    Package,
    RoutingAllocation,
    RoutingSimulator,
    NearestPkgPolicy,
    get_travel_time_estimate,
    sample_true_delivery_return_time,
    setup_routing_sim,
    update_routing_sim!,
    parse_routing_commandline,
    EuclideanLatLongMetric,
    LatLonCoords

# Routing algorithms - SCoBA and baselines
export
    scoba_routing!,
    RoutingMCTSMDP,
    update_routing_mcts_fullstate!,
    get_current_routing_mcts_state

include("conveyor_belt/conveyor_belt_types.jl")
include("conveyor_belt/conveyor_belt_simulator.jl")
include("conveyor_belt/conveyor_belt_scoba.jl")
include("conveyor_belt/conveyor_belt_baselines.jl")

include("routing/routing_types.jl")
include("routing/routing_simulator.jl")
include("routing/routing_scoba.jl")
include("routing/routing_baselines.jl")

end


text:
using Test, ePPR, FileIO, BenchmarkTools

## Simulated Data
# x: Natural Images
# y: Simulated Neuron with MFR=0.56, MIF=4.26
simdata=load(joinpath(@__DIR__,"simdata.jld2"));x=simdata["x"];y=simdata["y"];imagesize=simdata["imagesize"]
logdir = joinpath(@__DIR__,"log")

# ePPR with Linear Time Interaction
hp = ePPRHyperParams(imagesize...)
log = ePPRLog(debug=true,plot=true)
model,models = epprcv(x,y,hp,log)
plotmodel(model,hp)

# ePPR with Non-Linear Time Interaction
hp = ePPRHyperParams(imagesize...,ndelay=3,lambda=15,nft=[6])
log = ePPRLog(debug=true,plot=true)
model,models = epprcv(x,y,hp,log)
plotmodel(model,hp)

## ePPR with Hyper Parameter Search
hp = ePPRHyperParams(imagesize...,lambda=64,lambdascale=0.5)
log = ePPRLog(debug=true,plot=true,dir=logdir)
model,models = epprhypercv(x,y,hp,log)
plotmodel(model,hp)

## Data Recorded from an Anesthetized Cat
# x: Natural Images each present 40ms
# y1: Simple Cell Spike Sum
# y2: Complex Cell Spike Sum
data=load(joinpath(@__DIR__,"data.jld2"));x=data["x"];y1=data["y1"];y2=data["y2"];imagesize=data["imagesize"]

# ePPR with Linear Time Interaction
hp = ePPRHyperParams(imagesize...,lambda=15,blankcolor=mean(x))
log = ePPRLog(debug=true,plot=true)
model,models = epprcv(x,y1,hp,log)
plotmodel(model,hp)

# ePPR with Non-Linear Time Interaction
hp = ePPRHyperParams(imagesize...,ndelay=3,lambda=15,nft=[6],blankcolor=mean(x))
log = ePPRLog(debug=true,plot=true)
model,models = epprcv(x,y2,hp,log)
plotmodel(model,hp)

## Benchmark
hp = ePPRHyperParams(imagesize...)
@btime model,models = epprcv($x,$y1,$hp) samples=1 # ~19s


text:
using SoapySDR
using Test
using Unitful
using Unitful.DefaultSymbols
const dB = u"dB"
using Intervals

const sd = SoapySDR

const hardware = "loopback"

# Load dummy test harness or hardware
if hardware == "loopback"
    # build SoapyLoopback and dlopen it
    include("setup_loopback.jl")
elseif hardware == "rtlsdr"
    using SoapyRTLSDR_jll
else
    error("unknown test hardware")
end

@testset "SoapySDR.jl" begin
@testset "Ranges/Display" begin
    intervalrange = sd.SoapySDRRange(0, 1, 0)
    steprange = sd.SoapySDRRange(0, 1, 0.1)

    intervalrangedb = sd._gainrange(intervalrange)
    steprangedb = sd._gainrange(steprange) #TODO

    intervalrangehz = sd._hzrange(intervalrange)
    steprangehz = sd._hzrange(steprange)

    hztype = typeof(1.0*Hz)

    @test typeof(intervalrangedb) == Interval{Gain{Unitful.LogInfo{:Decibel, 10, 10}, :?, Float64}, Closed, Closed}
    @test typeof(steprangedb) == Interval{Gain{Unitful.LogInfo{:Decibel, 10, 10}, :?, Float64}, Closed, Closed}
    @test typeof(intervalrangehz) == Interval{hztype, Closed, Closed}
    @test typeof(steprangehz) == StepRangeLen{hztype, Base.TwicePrecision{hztype}, Base.TwicePrecision{hztype}}

    io = IOBuffer(read=true, write=true)

    sd.print_hz_range(io, intervalrangehz)
    @test String(take!(io)) == "00..0.001 kHz"
    sd.print_hz_range(io, steprangehz)
    @test String(take!(io)) == "00 Hz:0.0001 kHz:0.001 kHz"
end
@testset "High Level API" begin

    io = IOBuffer(read=true, write=true)

    # Device constructor, show, iterator
    @test length(Devices()) == 1
    show(io, Devices())
    deva = Devices()[1]
    deva["refclk"] = "internal"
    dev = open(deva)
    show(io, dev)
    for dev in Devices()
        show(io, dev)
    end
    @test typeof(dev) == sd.Device
    @test typeof(dev.info) == sd.OwnedKWArgs
    @test dev.driver == :LoopbackDriver
    @test dev.hardware == :LoopbackHardware
    @test dev.time_sources == SoapySDR.TimeSource[:sw_ticks,:hw_ticks]
    @test dev.time_source == SoapySDR.TimeSource(:sw_ticks)
    dev.time_source = "hw_ticks"
    @test dev.time_source == SoapySDR.TimeSource(:hw_ticks)
    dev.time_source = dev.time_sources[1]
    @test dev.time_source == SoapySDR.TimeSource(:sw_ticks)

    # Channels
    rx_chan_list = dev.rx
    tx_chan_list = dev.tx
    @test typeof(rx_chan_list) == sd.ChannelList
    @test typeof(tx_chan_list) == sd.ChannelList
    show(io, rx_chan_list)
    show(io, tx_chan_list)
    rx_chan = dev.rx[1]
    tx_chan = dev.tx[1]
    @test typeof(rx_chan) == sd.Channel
    @test typeof(tx_chan) == sd.Channel
    show(io, rx_chan)
    show(io, tx_chan)
    show(io, MIME"text/plain"(), rx_chan)
    show(io, MIME"text/plain"(), tx_chan)

    #channel set/get properties
    @test rx_chan.native_stream_format == SoapySDR.ComplexInt{12} #, fullscale
    @test rx_chan.stream_formats == [Complex{Int8}, SoapySDR.ComplexInt{12}, Complex{Int16}, ComplexF32]
    @test tx_chan.native_stream_format == SoapySDR.ComplexInt{12} #, fullscale
    @test tx_chan.stream_formats == [Complex{Int8}, SoapySDR.ComplexInt{12}, Complex{Int16}, ComplexF32]

    # channel gain tests
    @test rx_chan.gain_mode == false
    rx_chan.gain_mode = true
    @test rx_chan.gain_mode == true

    @test rx_chan.gain_elements == SoapySDR.GainElement[:IF1, :IF2, :IF3, :IF4, :IF5, :IF6, :TUNER]
    if1 = rx_chan.gain_elements[1]
    @show rx_chan[if1]
    rx_chan[if1] = 0.5u"dB"
    @test_broken rx_chan[if1] == 0.5u"dB"

    #@show rx_chan.gain_profile
    @show rx_chan.frequency_correction

    #@test tx_chan.bandwidth == 2.048e6u"Hz"
    #@test tx_chan.frequency == 1.0e8u"Hz"
    #@test tx_chan.gain == -53u"dB"
    #@test tx_chan.sample_rate == 2.048e6u"Hz"

    # setter/getter tests
    rx_chan.sample_rate = 1e5u"Hz"
    #@test rx_chan.sample_rate == 1e5u"Hz"

    rx_stream = sd.Stream(ComplexF32, [rx_chan])
    @test typeof(rx_stream) == sd.Stream{ComplexF32}
    tx_stream = sd.Stream(ComplexF32, [tx_chan])
    @test typeof(tx_stream) == sd.Stream{ComplexF32}

    rx_stream = sd.Stream([rx_chan])
    @test typeof(rx_stream) == sd.Stream{sd.ComplexInt{12}}
    tx_stream = sd.Stream([tx_chan])
    @test typeof(tx_stream) == sd.Stream{sd.ComplexInt{12}}

    #sd.activate!(rx_stream)
    #sd.activate!(tx_stream)
    #sd.deactivate!(rx_stream)
    #sd.deactivate!(tx_stream)

    # Stream Close
    close(rx_stream)
    close(tx_stream)
    #close(dev)
end
end


text:
mutable struct FunctionalMap
    mesh1
    mesh2
    descr1 # (n1,p) descriptors on the first mesh
    descr2 # (n2,p) descriptors on the second mesh
    D_a    # (k1,k1) area-based shape differnence operator
    D_c    # (k1,k1) conformal-based shape differnence operator
    FM_type # 'classic' | 'icp' | 'zoomout' which FM is currently used
    k1      # dimension of the first eigenspace (varies depending on the type of FM)
    k2      # dimension of the seconde eigenspace (varies depending on the type of FM)
    FM      # (k2,k1) current FM
    p2p     # (n2,) point to point map associated to the current functional map
end

function FunctionalMap(mesh1, mesh2)
    return FunctionalMap(mesh1, mesh2, undef, undef, undef, undef, "classic", undef, undef, undef, undef)
end

function preprocess!(fm::FunctionalMap; verbose=true)
end


text:
### All methods which PeaceVote users would use to interact with communities and communitcate.

function register(community::Community{T},certificate::Certificate) where T<:AbstractString
    CONFIG_DIR = community.config
    if iskey(community,:member)
        mname = fullname(community.m)[1]
        fname = "$CONFIG_DIR/keys/$mname/member"
        keyfile = open(fname)
        signer = deserialize(keyfile)

        @assert pubkeyid(certificate) == pubkeyid(signer,community)
    else
        error("The idkey had not been found.")
    end

    res = community.m.register(certificate)

    if res==true
        @info "Registration of the key had been succesfull"
    else
        error("Unsucessfull registration. $res")
    end
end

### Now I need to make an interface for registrating voters key

# using PeaceVote.Traits

function braidkey(community::Community,sign::Function)
    G = community.m.G
    signer = Signer(G)
    pid = pubkeyid(signer,community)
    res = community.m.braidkey(pid,sign) ### Depending on the signature pubkeyid server would say which gates to take.
    if res==true
        storekey(community,:voter,signer)
    else
        error("Unsuccesfull braiding. $res")
    end
end

function braidkey(community::Community)
    if iskey(community,:voter) ### voter key is only added if
        braidkey(community,community.voter)
    elseif iskey(community,:member)
        braidkey(community,community.member)
    else
        error("No keys found for the community.")
    end
end

### Now the last thing is just a function for vote

function vote(msg,community::Community)
    signature = community.voter(msg)
    community.m.vote(msg,signature)
end


text:
using Test
using OMEinsum
using OMEinsum: subindex

@testset "EinCode" begin
    code = EinCode(((1,2), (2,3)), (1,3))
    @test code isa EinCode
    @test OMEinsum.getixs(code) == ((1,2), (2,3))
    @test OMEinsum.getiy(code) == (1,3)
end

@testset "indexer" begin
    si = EinIndexer{()}(())
    @test subindex(si, (1,2,3)) == 1
    si = EinIndexer{(3,2)}((7,6))
    @test OMEinsum.getlocs(si) == (3,2)
    a = randn(7,6)
    @test a[subindex(si, (4,5,2))] == a[2,5]
end

@testset "EinArray" begin
    locs_xs = (EinIndexer{(1,2)}((8,8)), EinIndexer{(2,3)}((8,8)))
    ixs = ((1,2), (2,3))
    iy = (1,3)
    x1 = randn(8, 8)
    x2 = randn(8, 8)
    arr = einarray(EinCode(ixs, iy), (x1, x2), OMEinsum.get_size_dict(ixs, (x1, x2)))
    @test size(arr) == (8,8,8)
    # inner first, then outer
    @test arr[CartesianIndex((3,4,2))] == x1[4,3]*x2[3,2]
    @test arr[3,4,2] == x1[4,3]*x2[3,2]
end


text:
module LogRecords
  using HelloDB.Pages: Page, getint
  using HelloDB.Data: CHAR_SIZE, INT_SIZE, string_size

  abstract AbstractLogRecord

  type BasicLogRecord <: AbstractLogRecord
    page::Page
    pos::Int
  end

  function nextint(record::AbstractLogRecord)
    result = getint(record.page, record.pos)
    record.pos += INT_SIZE
    result
  end

  function next_string(record::AbstractLogRecord)
    result = getstring(record.page, record.pos)
    record.pos += string_size(length(result))
  end
end


text:
# test plotSLAM2D and plotPose features

using RoME, RoMEPlotting
using Test

##
@testset "test plotSLAM2D features" begin
##

fg = generateGraph_Hexagonal()

plotSLAM2D(fg, drawhist=false, drawPoints=false, drawContour=false, drawEllipse=false)
plotSLAM2D(fg, drawhist=false, drawPoints=false, drawContour=false, drawEllipse=true)
plotSLAM2D(fg, drawhist=false, drawPoints=false, drawContour=true, drawEllipse=false)
plotSLAM2D(fg, drawhist=false, drawPoints=false, drawContour=true, drawEllipse=true)
plotSLAM2D(fg, drawhist=false, drawPoints=true, drawContour=false, drawEllipse=false)
plotSLAM2D(fg, drawhist=false, drawPoints=true, drawContour=false, drawEllipse=true)
plotSLAM2D(fg, drawhist=false, drawPoints=true, drawContour=true, drawEllipse=false)
plotSLAM2D(fg, drawhist=false, drawPoints=true, drawContour=true, drawEllipse=true)

plotSLAM2D(fg, drawhist=true, drawPoints=false, drawContour=false, drawEllipse=false)
plotSLAM2D(fg, drawhist=true, drawPoints=false, drawContour=false, drawEllipse=true)
plotSLAM2D(fg, drawhist=true, drawPoints=false, drawContour=true, drawEllipse=false)
plotSLAM2D(fg, drawhist=true, drawPoints=false, drawContour=true, drawEllipse=true)
plotSLAM2D(fg, drawhist=true, drawPoints=true, drawContour=false, drawEllipse=false)
plotSLAM2D(fg, drawhist=true, drawPoints=true, drawContour=false, drawEllipse=true)
plotSLAM2D(fg, drawhist=true, drawPoints=true, drawContour=true, drawEllipse=false)
plotSLAM2D(fg, drawhist=true, drawPoints=true, drawContour=true, drawEllipse=true)

plotSLAM2D(fg; variableList=[:x1;:x5;:l1])

# plot runs but doesnt draw the levels at time of wrting the test.  Good to test the API regardless.
plotSLAM2D(fg, drawhist=false, drawPoints=true, drawContour=true, drawEllipse=false, levels=5)

## test with save and load of factor graph

saveDFG("/tmp/rp_test", fg)
fg_ = loadDFG("/tmp/rp_test.tar.gz")
# Base.rm("/tmp/rp_test.tar.gz")

plotSLAM2D(fg_, drawhist=true, drawPoints=true, drawContour=true, drawEllipse=true)

##
end

@testset "test CJL Docs example with plotSLAM2DPoses" begin
##

fg = generateGraph_Hexagonal()

# Draw the (x,y) marginal estimated belief contour for :x0, :x2, and Lx4
pl = plotKDE(fg, [:x0; :x2; :x4], c=["red";"green";"blue"], levels=2, dims=[1;2])

# add a few fun layers
pl3 = plotSLAM2DPoses(fg, regexPoses=r"x\d", from=3, to=3, drawContour=false, drawEllipse=true)
pl5 = plotSLAM2DPoses(fg, regexPoses=r"x\d", from=5, to=5, drawContour=false, drawEllipse=true, drawPoints=false)
pl_ = plotSLAM2DPoses(fg, drawContour=false, drawPoints=false, dyadScale=0.001, to=5)
union!(pl.layers, pl3.layers)
union!(pl.layers, pl5.layers)
union!(pl.layers, pl_.layers)

# change the plotting coordinates
pl.coord = Coord.Cartesian(xmin=-10,xmax=20, ymin=-1, ymax=25)

# save the plot to SVG and giving dedicated (although optional) sizing
pl |> SVG("/tmp/test.svg", 25cm, 15cm)

Base.rm("/tmp/test.svg")

##
end

##
@testset "test plotPose features" begin
##

fg = generateGraph_Hexagonal()

X1 = getBelief(fg, :x1)

plotPose(Pose2, X1)

plotPose(fg, :x1)
plotPose(fg, :x1, levels=1, scale=0.3)
plotPose(fg, [:x1;:x2], levels=2, c=["green"; "magenta"])

@warn "fix plotPose legend keyword"
# plotPose(fg, [:x1;:x2], levels=2, c=["green"; "magenta"], legend=["hello"; "world"])

##

fg_ = loadDFG("/tmp/rp_test.tar.gz")
Base.rm("/tmp/rp_test.tar.gz")

plotPose(fg, :x1)

##
end


text:
using ITensors
using ITensorMPOCompression
using Test
using Printf

#
# Single index tests
#
function test_qx(::Type{ElT},qx::Function,M::Int64,N::Int64,pos::Bool) where {ElT<:Number}
    ir=Index(M,"row")
    ic=Index(N,"col")
    A=randomITensor(ElT,ir,ic)
    Q,LR=qx(A,ir;positive=pos)
    norm(A-Q*LR)
end

function test_xq(::Type{ElT},xq::Function,M::Int64,N::Int64,pos::Bool) where {ElT<:Number}
    ir=Index(M,"row")
    ic=Index(N,"col")
    A=randomITensor(ElT,ir,ic)
    LR,Q=xq(A,ir;positive=pos)
    norm(A-LR*Q)
end

@testset "QR decomposition, single index" begin
    eps=1e-14
    M,N=5,5
    M1,N1=10,10
    @test test_qx(Float64   ,qr,M ,N,false)<eps
    @test test_qx(Float64   ,qr,M ,N1,false)<eps
    @test test_qx(Float64   ,qr,M1,N1,false)<eps
    @test test_qx(Float64   ,qr,M ,N ,true )<eps
    @test test_qx(Float64   ,qr,M ,N1,true )<eps
    @test test_qx(Float64   ,qr,M1,N1,true )<eps
    @test test_qx(ComplexF64,qr,M ,N ,false)<eps
    @test test_qx(ComplexF64,qr,M ,N1,false)<eps
    @test test_qx(ComplexF64,qr,M1,N1,false)<eps
    @test test_qx(ComplexF64,qr,M ,N ,true )<eps
    @test test_qx(ComplexF64,qr,M ,N1,true )<eps
    @test test_qx(ComplexF64,qr,M1,N1,true )<eps
end

@testset "QL decomposition, single index" begin
    eps=1e-14
    M,N=5,5
    M1,N1=10,10
    @test test_qx(Float64   ,ql,M ,N,false)<eps
    @test test_qx(Float64   ,ql,M ,N1,false)<eps
    @test test_qx(Float64   ,ql,M1,N1,false)<eps
    @test test_qx(Float64   ,ql,M ,N ,true )<eps
    @test test_qx(Float64   ,ql,M ,N1,true )<eps
    @test test_qx(Float64   ,ql,M1,N1,true )<eps
    @test test_qx(ComplexF64,ql,M ,N ,false)<eps
    @test test_qx(ComplexF64,ql,M ,N1,false)<eps
    @test test_qx(ComplexF64,ql,M1,N1,false)<eps
    @test test_qx(ComplexF64,ql,M ,N ,true )<eps
    @test test_qx(ComplexF64,ql,M ,N1,true )<eps
    @test test_qx(ComplexF64,ql,M1,N1,true )<eps
end

@testset "LQ decomposition, single index" begin
    eps=1e-14
    M,N=5,5
    M1,N1=10,10
    @test test_xq(Float64   ,lq,M ,N,false)<eps
    @test test_xq(Float64   ,lq,M ,N1,false)<eps
    @test test_xq(Float64   ,lq,M1,N1,false)<eps
    @test test_xq(Float64   ,lq,M ,N ,true )<eps
    @test test_xq(Float64   ,lq,M ,N1,true )<eps
    @test test_xq(Float64   ,lq,M1,N1,true )<eps
    @test test_xq(ComplexF64,lq,M ,N ,false)<eps
    @test test_xq(ComplexF64,lq,M ,N1,false)<eps
    @test test_xq(ComplexF64,lq,M1,N1,false)<eps
    @test test_xq(ComplexF64,lq,M ,N ,true )<eps
    @test test_xq(ComplexF64,lq,M ,N1,true )<eps
    @test test_xq(ComplexF64,lq,M1,N1,true )<eps
end

@testset "RQ decomposition, single index" begin
    eps=1e-14
    M,N=5,5
    M1,N1=10,10
    @test test_xq(Float64   ,rq,M ,N,false)<eps
    @test test_xq(Float64   ,rq,M ,N1,false)<eps
    @test test_xq(Float64   ,rq,M1,N1,false)<eps
    @test test_xq(Float64   ,rq,M ,N ,true )<eps
    @test test_xq(Float64   ,rq,M ,N1,true )<eps
    @test test_xq(Float64   ,rq,M1,N1,true )<eps
    @test test_xq(ComplexF64,rq,M ,N ,false)<eps
    @test test_xq(ComplexF64,rq,M ,N1,false)<eps
    @test test_xq(ComplexF64,rq,M1,N1,false)<eps
    @test test_xq(ComplexF64,rq,M ,N ,true )<eps
    @test test_xq(ComplexF64,rq,M ,N1,true )<eps
    @test test_xq(ComplexF64,rq,M1,N1,true )<eps
end

#
#  Test multiple indicies using MPO matrices
#
@testset "QR,QL,LQ,RQ decomposition fo MPO matrices" begin
    N=5
    NNN=4
    hx=0.5
    eps=2e-15
    sites = siteinds("SpinHalf", N)
    H=make_transIsing_MPO(sites,NNN,hx)
    W=H[2]
    d,n,r,c=parse_links(W)

    Lind=noncommoninds(inds(W),c)
    Rind=noncommoninds(inds(W),r)

    Q,R=qr(W,Lind;positive=true)
    iqr=commonindex(Q,R)
    Id=Q*prime(Q,iqr)
    Idq=delta(iqr,iqr')
    @test norm(Id-Idq)<eps
    @test norm(W-Q*R)<eps

    Q,L=ql(W,Lind;positive=true)
    iql=commonindex(Q,L)
    Id=Q*prime(Q,iql)
    Idq=delta(iql,iql')
    @test norm(Id-Idq)<eps
    @test norm(W-Q*L)<eps

    L,Q=lq(W,Rind;positive=true)
    iql=commonindex(Q,L)
    Id=Q*prime(Q,iql)
    Idq=delta(iql,iql')
    @test norm(Id-Idq)<eps
    @test norm(W-L*Q)<eps

    R,Q=rq(W,Rind;positive=true)
    iql=commonindex(Q,R)
    Id=Q*prime(Q,iql)
    Idq=delta(iql,iql')
    @test norm(Id-Idq)<eps
    @test norm(W-R*Q)<eps
end

#using Printf
#Base.show(io::IO, f::Float64) = @printf(io, "%1.3f", f)
#println("-----------Start--------------")

@testset "QR,QL,LQ,RQ decomposition with rank revealing" begin
    N=10
    NNN=6
    hx=0.5
    eps=2e-15
    sites = siteinds("SpinHalf", N)
    #
    #  use lower tri MPO to get some zero pivots for QL and RQ.
    #
    H=make_transIsing_MPO(sites,NNN,hx,lower)
    W=H[2]
    d,n,r,c=parse_links(W)

    Lind=noncommoninds(inds(W),c)
    Rind=noncommoninds(inds(W),r)
    @assert dim(c)==dim(r)
    #
    #  RQ decomp
    #
    R,Q=rq(W,Rind;positive=true,rank=true)
    iq=commonindex(Q,R)
    @printf "RQ decomposition %4i rows were removed from R\n" dim(c)-dim(iq)
    Id=Q*prime(Q,iq)
    Idq=delta(iq,iq')
    @test norm(Id-Idq)<eps
    @test norm(W-R*Q)<eps
    #
    #  QL decomp
    #
    Q,L=ql(W,Rind;positive=true,rank=true)
    iq=commonindex(Q,L)
    @printf "QL decomposition %4i rows were removed from L\n" dim(c)-dim(iq)
    Id=Q*prime(Q,iq)
    Idq=delta(iq,iq')
    @test norm(Id-Idq)<eps
    @test norm(W-L*Q)<eps

    #
    #  use upper tri MPO to get some zero pivots for LQ and QR.
    #
    H=make_transIsing_MPO(sites,NNN,hx,upper)
    W=H[2]
    d,n,r,c=parse_links(W)

    Lind=noncommoninds(inds(W),c)
    Rind=noncommoninds(inds(W),r)
    @assert dim(c)==dim(r)

    #
    #  QR decomp
    #
    Q,R=ITensorMPOCompression.qr(W,Rind;positive=true,rank=true)
    iq=commonindex(Q,R)
    @printf "QR decomposition %4i rows were removed from R\n" dim(c)-dim(iq)
    Id=Q*prime(Q,iq)
    Idq=delta(iq,iq')
    @test norm(Id-Idq)<eps
    @test norm(W-R*Q)<eps
    #
    #  LQ decomp
    #
    L,Q=lq(W,Rind;positive=true,rank=true)
    iq=commonindex(Q,L)
    @printf "LQ decomposition %4i rows were removed from L\n" dim(c)-dim(iq)
    Id=Q*prime(Q,iq)
    Idq=delta(iq,iq')
    @test norm(Id-Idq)<eps
    @test norm(W-L*Q)<eps

end
 #=
@testset "QR,QL,LQ,RQ decomposition for block sparse matrices" begin
    N=5
    NNN=2
    Dw=3
    hx=0.0 #can't make and sx op with QNs in play
    ul=lower
    lr=left
    ms=matrix_state(ul,lr)

    rowq=Index(QN("Sz",0)=>Dw;dir=ITensors.In,tags="Link,l=1") #link index with qn's
    colq=Index(QN("Sz",0)=>Dw;dir=ITensors.Out,tags="Link,l=2")
    sites = siteinds("S=1/2",N;conserve_qns=true) #sites with qn's
    H=make_transIsing_MPO(sites,NNN,hx,ul)
    W=H[2]
    pprint(W,1e-14)
    # Q,R,lq=block_qx(W,ul;orth=lr)

    off=V_offsets(ms)
    Vblock=getV(W,off)
    Rind=filterinds(Vblock,tags="l=2")
    Linds=noncommoninds(Vblock,Rind)
    Vd=dense(Vblock)
    Qd,Ld,iq=ql(Vd,Linds;positive=true,epsrr=0.0)
    set_scale!(Ld,Qd,off) #rescale so the L(n,n)==1.0
    pprint(Qd,1e-14)
    # @show Qd Ld

    # Vt=NDTensors.tensor(Vblock)
    # for b in eachnzblock(Vt)
    #     bl=NDTensors.blockview(Vt, b)
    #     @show Vt
    # end
    Q,L,iq=ql(Vblock,Linds;positive=true,epsrr=0.0)
    # U,s,V,iu,iv=svd(Vblock,Linds)
    #@show U V
end =#


text:
using IterTools

encode(str::String) = collect((length(g), first(g)) for g in groupby(first, str))
decode(cod::Vector) = join(repeat("$l", n) for (n, l) in cod)

for original in ["aaaaahhhhhhmmmmmmmuiiiiiiiaaaaaa", "WWWWWWWWWWWWBWWWWWWWWWWWWBBBWWWWWWWWWWWWWWWWWWWWWWWWBWWWWWWWWWWWWWW"]
    encoded = encode(original)
    decoded = decode(encoded)
    println("Original: $original\n -> encoded: $encoded\n -> decoded: $decoded")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
collater:
- Julia version: 1.1.0
- Author: bramb
- Date: 2019-03-15
=#

using MetaGraphs,Combinatorics

@enum(CollationState,
     needs_witness, ready_to_collate, is_collated
)

# import Base.isequal, Base.hash

mutable struct CollationGraph
    sigils::Vector{String}
    graph::MetaDiGraph

    CollationGraph() = new([], MetaDiGraph(SimpleDiGraph()))
end

function rank(cg::CollationGraph)
    ranking = Ranking()
    nodes_to_rank = []
    return ranking
end

function to_dot(cg::CollationGraph)
    digraph = _collationgraph_as_dot(cg)
    dot = """
    digraph CollationGraph {
        rankdir=LR
        labelloc=b
        color=white
        edge [arrowsize=0.5]
        $digraph
    }
    """
    return dot
end

function _collationgraph_as_dot(cg::CollationGraph)
    mg = cg.graph
    nodes_buf = IOBuffer()
    for n in 1:nv(mg)
        type = get_prop(mg,n,:type)
        if (type == TEXTNODE)
            text = get_prop(mg,n,:text)
            node_def = """v$n[shape=box;label="$text"]"""
        else
            node_def = """v$n[shape=circle;width=0.05;label=""]"""
        end
        println(nodes_buf,node_def)
    end
    nodes = String(take!(nodes_buf))

    edges_buf = IOBuffer()
    for e in edges(mg)
        edge_def = """v$(e.src) -> v$(e.dst)"""
        println(edges_buf,edge_def)
    end
    edgesstring = String(take!(edges_buf))

    dot = """
    $nodes
    $edgesstring
    """
    return dot
end

function to_html(cg::CollationGraph)
    stringbuf = IOBuffer()
    for s in reverse(cg.sigils)
        columns = "some text	"
        row = """

$s	$columns

        """
        println(stringbuf,row)
    end
    tablerows = String(take!(stringbuf))
    html =
    """

    $tablerows

    """
    return html
end

mutable struct Collation
    variantwitness_graphs::Dict{String,MetaDiGraph}
    state::CollationState
    graph::CollationGraph

    Collation() = new(Dict{String,MetaDiGraph}(),needs_witness,CollationGraph())
end

function add_witness!(collation::Collation, sigil::String, xml::String)
    wvg = to_graph(xml)
    set_props!(wvg,Dict(:sigil => sigil))
    collation.variantwitness_graphs[sigil] = wvg
    collation.state = length(collation.variantwitness_graphs) < 2 ? needs_witness : ready_to_collate
end

mutable struct Match
    witness_vertex_map::Dict{String,Integer}
    ranking_map::Dict{String,Integer}

    function Match(s1::String, v1::Integer, s2::String, v2::Integer)
        m = new(Dict{String,Integer}(),Dict{String,Integer}())
        m.witness_vertex_map[s1] = v1
        m.witness_vertex_map[s2] = v2
        return m
    end
end

set_rank!(match::Match, sigil::String, rank::Integer) = match.ranking_map[sigil] = rank

has_witness(m::Match, sigil::String) = haskey(m.witness_vertex_map,sigil)

# sigils(m::Match) = keys(m.witness_vertex_map)
function sigils(m::Match)
    return collect(keys(m.witness_vertex_map))
end

lowest_rank_for_witnesses_other_than(m::Match, s::String) = minimum([p[2] for p in m.ranking_map if p[1] != s])

mutable struct Ranking
    by_vertex::Dict{Integer,Integer}
    by_rank::Dict{Integer,Set{Integer}}

    Ranking() = new(Dict{Integer,Integer}(),Dict{Integer,Set{Integer}}())
end

apply(ranking::Ranking, vertex::Integer) = ranking.by_vertex[vertex]

function collate!(collation::Collation)
    sigils = sort(collect(keys(collation.variantwitness_graphs)))
    witnesses = [collation.variantwitness_graphs[s] for s in sigils]
    rankings = [ranking(wg) for wg in witnesses]

    matches = potential_matches(witnesses,rankings)
    @debug for m in matches
        @show(m)
    end
    collated_vertex_map = Dict{Integer,Integer}()
    first = popfirst!(witnesses)
    initialize(collation.graph,collated_vertex_map,first)
    matches_sorted_by_rank_per_witness = sort_and_filter_matches_by_witness(matches, sigils)
    @debug(matches_sorted_by_rank_per_witness)
    for s in sigils
        sorted_matches = matches_sorted_by_rank_per_witness[s]
        witness = collation.variantwitness_graphs[s]
        collate(collation.graph,witness,sorted_matches,collated_vertex_map)
    end
    collation.state = is_collated
    return collation
end

function potential_matches(witnesses, rankings)
    all_potential_matches = Set{Match}()
    vertex_to_match = Dict{Integer,Match}()

    for tuple in [t for t in permutations(1:length(witnesses), 2) if t[1] < t[2]]
        witness1 = witnesses[tuple[1]]
        ranking1 = rankings[tuple[1]]
        witness2 = witnesses[tuple[2]]
        ranking2 = rankings[tuple[2]]
        match!(all_potential_matches, vertex_to_match,witness1,witness2,ranking1,ranking2)
#         end_match = get_end_match(witness1,ranking1,witness2,ranking2)
#         push!(all_potential_matches, end_match)
    end
    return all_potential_matches
end

function match!(all_potential_matches, vertex_to_match, witness1, witness2, ranking1, ranking2)
    traversal1 = traversal(witness1)
    traversal2 = traversal(witness2)
    sigil1 = get_sigil(witness1)
    sigil2 = get_sigil(witness2)
    for v1 in traversal1
        text1 = get_prop(witness1,v1,:text)
        for v2 in traversal2
            text2 = get_prop(witness2,v2,:text)
#             println("$sigil1:$v1:$text1 == $sigil2:$v2:$text2 ?")
            if text1 == text2
                match = Match(sigil1, v1, sigil2, v2)
                set_rank!(match, sigil1, apply(ranking1,v1))
                set_rank!(match, sigil2, apply(ranking2,v2))
                push!(all_potential_matches, match)
                vertex_to_match[v1] = match
                vertex_to_match[v2] = match
            end
        end
    end
end

traversal(witness) = topological_sort_by_dfs(witness)

# function get_end_match(witness1, ranking1, witness2, ranking2)
#     endvertex1 = nv(witness1)
#     endvertex2 = nv(witness2)
#     endmatch = Match(endvertex1,endvertex2)
# end

function ranking(variantwitness_graph::MetaDiGraph)
    ranking = Ranking()
    for v in vertices(variantwitness_graph)
        rank = -1
        for inv in inneighbors(variantwitness_graph,v)
            if !haskey(ranking.by_vertex,inv)
                ranking.by_vertex[inv] = -1
            end
            rank = max(rank,ranking.by_vertex[inv])
        end
        rank += 1
        ranking.by_vertex[v] = rank;
        if !haskey(ranking.by_rank,rank)
            ranking.by_rank[rank] = Set()
        end
        push!(ranking.by_rank[rank],v);
    end
    return ranking
end

function initialize(collation_graph, collated_vertex_map, witnessgraph)
    sigil = get_sigil(witnessgraph)
    push!(collation_graph.sigils,sigil)
#     addMarkupNodes(collationGraph, markupNodeIndex, witnessGraph);
    collated_vertex_map[1] = 1
#     collatedTokenVertexMap.put(witnessGraph.getStartTokenVertex(), collationGraph.getTextStartNode());
    for v in [x for x in vertices(witnessgraph) if get_prop(witnessgraph,x,:type) == TEXTNODE]
        add_collation_node!(collation_graph,collated_vertex_map,v,witnessgraph)
    end
#     collated_vertex_map[witnessgraph]
#     collatedTokenVertexMap.put(witnessGraph.getEndTokenVertex(), collationGraph.getTextEndNode());
#     addEdges(collationGraph, collatedTokenVertexMap);
end

function get_sigil(mg::MetaDiGraph)
    get_prop(mg,:sigil)
end

function sort_and_filter_matches_by_witness(matches, sigils)
    return Dict{String,Vector{Match}}(s => filter_and_sort_matches_for_witness(matches,s) for s in sigils)
end

function filter_and_sort_matches_for_witness(matches,sigil)
    function my_isless(match1::Match,match2::Match)
        rank1 = match1.ranking_map[sigil]
        rank2 = match2.ranking_map[sigil]
        if rank1 == rank2
            rank1 = lowest_rank_for_witnesses_other_than(match1,sigil)
            rank2 = lowest_rank_for_witnesses_other_than(match2,sigil)
        end
        return rank1 < rank2
    end
    @debug(matches)
    filtered_matches = [m for m in matches if has_witness(m,sigil)]
    @debug(filtered_matches)
    return collect(sort(filtered_matches, lt=my_isless))
end

function add_collation_node!(collation_graph,collated_vertex_map,v,witnessgraph)
    if !haskey(collated_vertex_map,v)

#       TextNode collationNode = collationGraph.addTextNodeWithTokens(tokenVertex.getToken());
#       collationNode.addBranchPath(tokenVertex.getSigil(), tokenVertex.getBranchPath());
#         collated_vertex_map[v] = collationnode
#       collatedTokenVertexMap.put(tokenVertex, collationNode);
#       addMarkupHyperEdges(collationGraph, witnessGraph, markupNodeIndex, tokenVertex, collationNode);
    end
end

function collate(collation_graph::CollationGraph,witness,sorted_matches,collated_vertex_map)
    base_ranking = rank(collation_graph)

#     sigils = collation_graph.sigils
#     filtered_matches = filter(m -> !isempty(intersect(sigils(m),sigils)),sorted_matches)
    p(m::Match) = !isempty(sigils(m))
    filtered_matches = filter(p ,sorted_matches)
    witnesssigil = get_sigil(witness)
end