text:
# OptiGas: Optimize Gas Network using Benders Algorithm
# Author: Edward J. Xu, edxu96@outlook.com
# Date: Aug 12, 2019

using LinearAlgebra
using MatrixOptim

function getData()
    vec_nameNodes = ["A", "B", "C", "D", "E", "F", "G", "H", "I", "J"]
    numNodes = length(vec_nameNodes)
    vec_xNode = [1    3    8    2    4    5    6    7    9    9]
    vec_yNode = [1    2    1    7    9    5    3    7    9    4]
    mat_arcTwoNodes = [
        0    1    0    0    0    0    0    0    0    0;
        0    0    0    0    0    1    0    0    0    0;
        0    0    0    0    0    0    0    0    0    1;
        1    0    0    0    1    0    0    0    0    0;
        0    0    0    0    0    1    0    0    0    0;
        1    0    0    1    0    0    0    1    0    0;
        0    1    1    0    0    0    0    0    0    0;
        0    0    0    0    0    0    0    0    1    0;
        0    0    0    0    1    0    0    0    0    0;
        0    0    0    0    0    1    1    0    1    0
    ]
    vec_netEject = [0    0   0    0    0   90    0    0   0   80 ]
    vec_netInject = [10   30  20   10   40   0   25   20  15    0 ]
    mat_distance = zeros(Float64, numNodes, numNodes)
    for m = 1: numNodes
        for n = 1: numNodes
            mat_distance[m, n] = floor(sqrt( (vec_xNode[m] - vec_xNode[n]) * (vec_xNode[m] - vec_xNode[n]) +
                                       (vec_yNode[m] - vec_yNode[n]) * (vec_yNode[m] - vec_yNode[n]) ) )
        end
    end
    mat_fixedCost = zeros(Float64, numNodes, numNodes)
    for m = 1: numNodes
        for n = 1: numNodes
            mat_fixedCost[m, n] = 10 * mat_distance[m, n]
        end
    end
    return numNodes, vec_netEject, vec_netInject, mat_distance, mat_fixedCost
end

function getModelQ4(numNodes, vec_netEject, vec_netInject, mat_distance, mat_fixedCost)
    # Fomulation of mat_aCap
    mat_aCap_1 = zeros(numNodes^2, numNodes^2)
    for mn = 1: numNodes^2
        mat_aCap_1[mn, mn] = -1
    end
    # mat_a2_1 = zeros(numNodes^2, numNodes^2)
    mat_aCap_3 = zeros(numNodes, numNodes^2)
    seq_zeroTenNinety = collect(0:10:90)
    for m = 1: numNodes
        mat_aCap_3[m, (10 * m - 9): (10 * m)] = repeat([-1], 10)  # sent out
        mat_aCap_3[m, (seq_zeroTenNinety + repeat([m],10))] = repeat([1], 10)  # sent in
        mat_aCap_3[:, (11 * m - 10)] = repeat([0], 10)  # self-sending doesn't count
    end
    mat_aCap = vcat(mat_aCap_1, mat_aCap_3)

    # Fomulation of Matrix B
    mat_bCap_1 = zeros(numNodes^2, numNodes^2)
    for mn = 1: numNodes^2
        mat_bCap_1[mn, mn] = 170
    end
    # mat_b2_1 = Diagnal(repeat([1], numNodes^2))
    mat_bCap_3 = zeros(numNodes, numNodes^2)
    mat_bCap = vcat(mat_bCap_1, mat_bCap_3)

    # Fomulation of Vector b
    vec_bCap_1 = zeros(numNodes^2)
    # vec_b2_1 = zeros(numNodes^2)
    # for m = 1: numNodes
    #     for n = 1: numNodes
    #         vec_b2_1[10 * (m - 1) + n] = mat_arcTwoNodes[m, n]
    #     end
    # end
    vec_bCap_3 = zeros(numNodes)
    for n = 1: numNodes
        vec_bCap_3[n] = vec_netInject[n] - vec_netEject[n]
    end
    vec_b = vcat(vec_bCap_1, vec_bCap_3)
    vec_b = hcat(vec_b)

    vec_y_max = repeat([1], numNodes^2)
    vec_y_max = hcat(vec_y_max)
    vec_y_min = zeros(numNodes^2)
    for m = 1: numNodes
        for n = 1: numNodes
            vec_y_min[10 * (m - 1) + n] = mat_arcTwoNodes[m, n]
        end
    end

    vec_c = zeros(numNodes^2)
    for m = 1: numNodes
        for n = 1: numNodes
            vec_c[10 * (m - 1) + n] = mat_distance[m, n]
        end
    end
    vec_c = hcat(vec_c)

    vec_f = zeros(numNodes^2)
    for m = 1: numNodes
        for n = 1: numNodes
            vec_f[10 * (m - 1) + n] = mat_fixedCost[m, n]
        end
    end
    vec_f = hcat(vec_f)

    return vec_c, mat_aCap, vec_b, vec_f, vec_bCap, vec_y_min, vec_y_max
end

function main()
    numNodes, vec_netEject, vec_netInject, mat_distance, mat_fixedCost = getData()
    vec_c, mat_aCap, vec_b, vec_f, vec_bCap, vec_y_min, vec_y_max = getModelQ4(
        numNodes, vec_netEject, vec_netInject, mat_distance, mat_fixedCost
    )
    solveBendersMilp(
        n_x = numNodes^2,
        n_y = numNodes^2,
        vec_min_y = vec_y_min,
        vec_max_y = vec_y_max,
        vec_c = vec_c,
        vec_f = vec_f,
        vec_b = vec_b,
        mat_a = mat_aCap,
        mat_b = mat_bCap,
        epsilon = 0.0001,
        timesIterationMax = 1000
    )
    # 1280 - sum(mat_arcTwoNodes .* mat_fixedCost)
end

main()


text:
#############################################################################
# trace_mpower(A, t, C) returns tr(C*A^t) where A and C are positive definite
# matrices and C is constant and t ∈ [-1, 2].
#
# When t ∈ [0,1], trace_mpower(A, t, C) is concave in A (for fixed
# positive semidefinite matrix C) and convex for t ∈ [-1,0] or [1,2].
#
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#
#REFERENCE
#   Ported from CVXQUAD which is based on the paper: "Lieb's concavity
#   theorem, matrix geometric means and semidefinite optimization" by Hamza
#   Fawzi and James Saunderson (arXiv:1512.03401)
#############################################################################

struct TraceMpower <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr,}
    size::Tuple{Int, Int}
    C::AbstractMatrix
    t::Rational

    function TraceMpower(A::AbstractExpr, t::Rational, C::AbstractMatrix)
        children = (A,)
        if size(A) != size(C)
            throw(DimensionMismatch("A and C must be the same size"))
        end
        n = size(A)[1]
        if size(A) != (n, n)
            throw(DimensionMismatch("A and C must be square"))
        end
        if norm(C - C') > 1e-6
            throw(DomainError(C, "C must be Hermitian"))
        end
        if any(eigvals(Hermitian(C)) .< -1e-6)
            throw(DomainError(C, "C must be positive semidefinite"))
        end
        if t < -1 || t > 2
            throw(DomainError(t, "t must be in the range [-1, 2]"))
        end
        return new(:trace_mpower, hash(children), children, (1, 1), C, t)
    end
end

function sign(atom::TraceMpower)
    return NoSign()
end

function monotonicity(atom::TraceMpower)
    return (NoMonotonicity(),)
end

function curvature(atom::TraceMpower)
    t = atom.t
    if t >= 0 && t <= 1
        return ConcaveVexity()
    else
        return ConvexVexity()
    end
end

function evaluate(atom::TraceMpower)
    A = evaluate(atom.children[1])
    return trace_mpower(A, atom.t, atom.C)
end

const MatrixOrConstant = Union{AbstractMatrix, Constant}

function trace_mpower(A::AbstractExpr, t::Rational, C::MatrixOrConstant)
    #println("trace_mpower general case")
    return TraceMpower(A, t, evaluate(C))
end

function trace_mpower(A::MatrixOrConstant, t::Rational, C::MatrixOrConstant)
    #println("trace_mpower constant A")
    return tr(C * A^t)
end

trace_mpower(A::AbstractExprOrValue, t::Integer, C::MatrixOrConstant) = trace_mpower(A, t//1, C)

function conic_form!(atom::TraceMpower, unique_conic_forms)
    if !has_conic_form(unique_conic_forms, atom)
        A = atom.children[1]
        C = atom.C
        t = atom.t
        eye = Matrix(1.0*I, size(A))

        is_complex = sign(A) == ComplexSign()
        if is_complex
            make_temporary = () -> HermitianSemidefinite(size(A)[1])
        else
            make_temporary = () -> Semidefinite(size(A)[1])
        end

        T = make_temporary()

        if t >= 0 && t <= 1
            conic_form!(T in GeomMeanHypoCone(eye, A, t, false), unique_conic_forms)
            # It's already a real mathematically, but need to make it a real type.
            u = real(tr(C * T))
            cache_conic_form!(unique_conic_forms, atom, maximize(u))
        else
            conic_form!(T in GeomMeanEpiCone(eye, A, t, false), unique_conic_forms)
            # It's already a real mathematically, but need to make it a real type.
            u = real(tr(C * T))
            cache_conic_form!(unique_conic_forms, atom, minimize(u))
        end
    end
    return get_conic_form(unique_conic_forms, atom)
end

ext: jl


text:
module ApplyTests

using Test
using InplaceArrays.Arrays
using FillArrays

a = rand(3,2,4)
test_array(a,a)

a = rand(3,2,4)
a = CartesianIndices(a)
test_array(a,a)

a = rand(3,2)
a = CartesianIndices(a)
c = apply(-,a)
test_array(c,-a)

a = rand(12)
c = apply(-,a)
test_array(c,-a)

a = rand(12)
b = rand(12)
c = apply(-,a,b)
test_array(c,a.-b)

c = apply(Float64,-,a,b)
test_array(c,a.-b)

a = rand(0)
b = rand(0)
c = apply(-,a,b)
test_array(c,a.-b)

a = fill(rand(2,3),12)
b = rand(12)
c = apply(bcast(-),a,b)
test_array(c,[ai.-bi for (ai,bi) in zip(a,b)])

a = fill(rand(2,3),0)
b = rand(0)
c = apply(bcast(-),a,b)
test_array(c,[ai.-bi for (ai,bi) in zip(a,b)])

a = fill(rand(2,3),12)
b = rand(12)
c = apply(bcast(-),a,b)
d = apply(bcast(+),a,c)
e = apply(bcast(*),d,c)
test_array(e,[((ai.-bi).+ai).*(ai.-bi) for (ai,bi) in zip(a,b)])

a = fill(rand(Int,2,3),12)
b = fill(rand(Int,1,3),12)
c = array_caches(a,b)
i = 1
v = getitems!(c,(a,b),i)
@test c == (nothing,nothing)
@test v == (a[i],b[i])

a = fill(rand(Int,2,3),12)
b = fill(rand(Int,1,3),12)
ai = testitem(a)
@test ai == a[1]
ai, bi = testitems(a,b)
@test ai == a[1]
@test bi == b[1]

a = fill(rand(Int,2,3),0)
b = fill(1,0)
ai = testitem(a)
@test ai == Array{Int,2}(undef,0,0)
ai, bi = testitems(a,b)
@test ai == Array{Int,2}(undef,0,0)
@test bi == zero(Int)

a = fill(+,10)
x = rand(10)
y = rand(10)
v = apply(a,x,y)
r = [(xi+yi) for (xi,yi) in zip(x,y)]
test_array(v,r)
v = apply(Float64,a,x,y)
test_array(v,r)

a = Fill(bcast(+),10)
x = [rand(2,3) for i in 1:10]
y = [rand(1,3) for i in 1:10]
v = apply(a,x,y)
r = [(xi.+yi) for (xi,yi) in zip(x,y)]
test_array(v,r)

a = Fill(bcast(+),10)
x = [rand(mod(i-1,3)+1,3) for i in 1:10]
y = [rand(1,3) for i in 1:10]
v = apply(a,x,y)
r = [(xi.+yi) for (xi,yi) in zip(x,y)]
test_array(v,r)

# Test the intermediate results caching mechanism

struct ArrayWithCounter{T,N,A,C} <: AbstractArray{T,N}
  array::A
  counter::C
  function ArrayWithCounter(a::AbstractArray{T,N}) where {T,N}
    c = zeros(Int,size(a))
    c[:] .= 0
    new{T,N,typeof(a),typeof(c)}(a,c)
  end
end

Base.size(a::ArrayWithCounter) = size(a.array)

function Base.getindex(a::ArrayWithCounter,i::Integer...)
  a.counter[i...] += 1
  a.array[i...]
end

Base.IndexStyle(::Type{<:ArrayWithCounter{T,N,A}}) where {T,A,N} = IndexStyle(A)

function reset_counter!(a::ArrayWithCounter)
  a.counter[:] .= 0
end

a = ArrayWithCounter(fill(rand(2,3),12))
b = ArrayWithCounter(rand(12))
c = apply(bcast(-),a,b)
d = apply(bcast(+),a,c)
e = apply(bcast(*),d,c)
r = [ (ai.-bi).*(ai.+(ai.-bi)) for (ai,bi) in zip(a,b)]
cache = array_cache(e)
reset_counter!(a)
reset_counter!(b)
for i in 1:length(e)
  ei = getindex!(cache,e,i)
  ei = getindex!(cache,e,i)
  ei = getindex!(cache,e,i)
end

@test all(a.counter .== 2)
@test all(b.counter .== 1)

l = 10
ai = 1.0
bi = 2.0
a = Fill(ai,l)
b = Fill(bi,l)
c = apply(+,a,b)
r = map(+,a,b)
test_array(c,r)
@test isa(c,Fill)

#f = rand(10)
#a = rand(10)
#@test apply(f,a) === f
#
#f = fill(rand(4),10)
#a = rand(10)
#@test apply(f,a) === f
#
#f = fill(rand(4),10)
#g = fill(rand(4),10)
#a = rand(10)
#h = apply_all((f,g),a)
#@test h[1] === f
#@test h[2] === g

#l = 10
#ai = 1.0
#bi = 2.0
#a = fill(ai,l)
#b = fill(bi,l)
#c = apply(+,a,b)
#d = apply(c,a,b)
#@test c == d
#@test typeof(d) == typeof(c)

end # module


text:
@symbol_func function P_T(cur_reactor::AbstractReactor)
  cur_P = P_F(cur_reactor)

  cur_P *= ( 1 + E_Li / E_F )

  cur_P
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Gumbo_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Gumbo")
JLLWrappers.@generate_main_file("Gumbo", UUID("528830af-5a63-567c-a44a-034ed33b8444"))
end  # module Gumbo_jll


text:
module Examples

using CombinatorialSpaces
using CombinatorialSpaces.ExteriorCalculus
using Catlab
using Catlab.Present
using Catlab.Programs
using Catlab.WiringDiagrams
using Catlab.Graphics
using Catlab.CategoricalAlgebra
using Decapodes.Simulations
using Decapodes.Diagrams
using Decapodes.Schedules
using CombinatorialSpaces: ∧

using LinearAlgebra

import Catlab.Graphics: to_graphviz

export dual, sym2func, to_graphviz, gen_dec_rules

function dual(s::EmbeddedDeltaSet2D{O, P}) where {O, P}
  sd = EmbeddedDeltaDualComplex2D{O, eltype(P), P}(s)
  subdivide_duals!(sd, Barycenter())
  sd
end

sym2func(sd) = begin
  Dict(:d₀=>Dict(:operator => d(Val{0}, sd), :type => MatrixFunc()),
       :d₁=>Dict(:operator => d(Val{1}, sd), :type => MatrixFunc()),
       :dual_d₀=>Dict(:operator => dual_derivative(Val{0}, sd),
                      :type => MatrixFunc()),
       :dual_d₁=>Dict(:operator => dual_derivative(Val{1}, sd),
                      :type => MatrixFunc()),
       :⋆₀=>Dict(:operator => hodge_star(Val{0}, sd), :type => MatrixFunc()),
       :⋆₁=>Dict(:operator => hodge_star(Val{1}, sd), :type => MatrixFunc()),
       :⋆₂=>Dict(:operator => hodge_star(Val{2}, sd), :type => MatrixFunc()),
       :⋆₀⁻¹=>Dict(:operator => inv_hodge_star(Val{0}, sd), :type => MatrixFunc()),
       :⋆₁⁻¹=>Dict(:operator => inv_hodge_star(Val{1}, sd; hodge=DiagonalHodge()),
                   :type => MatrixFunc()),
       :⋆₂⁻¹=>Dict(:operator => inv_hodge_star(Val{2}, sd), :type => MatrixFunc()),
       :∧₁₀=>Dict(:operator => (γ, α,β)->(γ .= ∧(Tuple{1,0}, sd, α, β)),
                  :type=>InPlaceFunc()),
       :∧₁₁=>Dict(:operator => (γ, α,β)->(γ .= ∧(Tuple{1,1}, sd, α, β)),
                  :type=>InPlaceFunc()),
       :plus => Dict(:operator => (x,y)->(x+y), :type => ElementwiseFunc()),
       :L₀ => Dict(:operator => (v,u)->(lie_derivative_flat(Val{2}, sd, v, u)),
                   :type => ArbitraryFunc()))
#       :Δ₀=>Dict(:operator => Δ(Val{0}, sd), :type => MatrixFunc()),
#       :Δ₁=>Dict(:operator => Δ(Val{1}, sd), :type => MatrixFunc()),
end

function expand_dwd(dwd_orig, patterns)
  dwd = deepcopy(dwd_orig)
  has_updated = true
  while(has_updated)
    updates = map(1:nboxes(dwd)) do b
      v = box(dwd, b).value
      if v ∈ keys(patterns)
        dwd.diagram[b, :value] = patterns[v]
        true
      else
        false
      end
    end
    dwd = substitute(dwd)
    has_updated = any(updates)
  end
  dwd
end

function contract_matrices!(dwd, s2f)
  has_updated = true
  is_matrix(b) = begin
    b >= 1 || return false
    v = Symbol(first(split("$(box(dwd, b).value)", "⋅")))
    v != :remove && (s2f[v][:type] isa MatrixFunc) &&
      length(input_ports(dwd, b)) == 1 &&
      length(output_ports(dwd, b)) == 1
  end
  while(has_updated)
    has_updated = false
    # Join parallel matrix boxes and mark for removal
    for b in 1:nboxes(dwd)
      if is_matrix(b)
        ows = out_wires(dwd, b)
        if length(ows) == 1
          ow = ows[1]
          if is_matrix(ow.target.box)
            has_updated = true
            tb = ow.target.box
            dwd.diagram[b, :value] = Symbol(dwd.diagram[tb, :value], :⋅, dwd.diagram[b, :value])
            dwd.diagram[tb, :value] = :remove
            add_wires!(dwd, map(w -> Wire(w.value, Port(b, OutputPort, 1), w.target),
                                out_wires(dwd, tb)))
            dwd.diagram[incident(dwd.diagram, b, :out_port_box)[1], :out_port_type] = output_ports(dwd, tb)[1]
          end
        end
      end
    end
    rem_boxes!(dwd, findall(b -> dwd.diagram[b, :value] == :remove, 1:nboxes(dwd)))
  end
end

function boundary_inds(::Type{Val{1}}, s)
  collect(findall(x -> x != 0, boundary(Val{2},s) * fill(1,ntriangles(s))))
end

function boundary_inds(::Type{Val{0}}, s)
  ∂1_inds = boundary_inds(Val{1}, s)
  unique(vcat(s[∂1_inds,:src],s[∂1_inds,:tgt]))
end

function boundary_inds(::Type{Val{2}}, s)
  ∂1_inds = boundary_inds(Val{1}, s)
  inds = map([:∂e0, :∂e1, :∂e2]) do esym
    vcat(incident(s, ∂1_inds, esym)...)
  end
  unique(vcat(inds...))
end

function bound_edges(s, ∂₀)
  te = vcat(incident(s, ∂₀, :tgt)...)
  se = vcat(incident(s, ∂₀, :src)...)
  intersect(te, se)
end

function adj_edges(s, ∂₀)
  te = vcat(incident(s, ∂₀, :tgt)...)
  se = vcat(incident(s, ∂₀, :src)...)
  unique(vcat(te, se))
end

function gen_dec_rules()
  @present ExtendedOperators(FreeExtCalc2D) begin
    X::Space
    neg₁̃::Hom(DualForm1(X), DualForm1(X)) # negative
    half₁̃::Hom(DualForm1(X), DualForm1(X)) # half
    proj₁_¹⁰₁::Hom(Form1(X) ⊗ Form0(X), Form1(X))
    proj₂_¹⁰₀::Hom(Form1(X) ⊗ Form0(X), Form0(X))
    proj₁_¹²₁::Hom(Form1(X) ⊗ Form2(X), Form1(X))
    proj₂_¹²₂::Hom(Form1(X) ⊗ Form2(X), Form2(X))
    proj₁_¹²̃₁::Hom(Form1(X) ⊗ DualForm2(X), Form1(X))
    proj₂_¹²̃₂̃::Hom(Form1(X) ⊗ DualForm2(X), DualForm2(X))
    proj₁_¹¹̃₁::Hom(Form1(X) ⊗ DualForm1(X), Form1(X))
    proj₂_¹¹̃₁̃::Hom(Form1(X) ⊗ DualForm1(X), DualForm1(X))
    proj₁_¹̃¹̃₁̃::Hom(DualForm1(X) ⊗ DualForm1(X), DualForm1(X))
    proj₂_¹̃¹̃₁̃::Hom(DualForm1(X) ⊗ DualForm1(X), DualForm1(X))
    proj₁_¹¹₁::Hom(Form1(X)⊗Form1(X), Form1(X))
    proj₂_¹¹₁::Hom(Form1(X)⊗Form1(X), Form1(X))
    sum₁̃::Hom(DualForm1(X)⊗DualForm1(X), DualForm1(X))
    sum₁::Hom(Form1(X)⊗Form1(X), Form1(X))
    L₀::Hom(Form1(X)⊗DualForm2(X), DualForm2(X))
    L₁::Hom(Form1(X)⊗DualForm1(X), DualForm1(X))
    i₀::Hom(Form1(X)⊗DualForm2(X), DualForm1(X))
    i₁::Hom(Form1(X)⊗DualForm1(X), DualForm0(X))
  end

  Lie0Imp = @decapode ExtendedOperators begin
    (dḞ2, dF2)::DualForm2{X}
    F1::Form1{X}

    dḞ2 ==  dual_d₁{X}(i₀(F1, dF2))
  end
  lie0_imp = diag2dwd(Lie0Imp, in_vars = [:F1, :dF2], out_vars = [:dḞ2])

  Lie1Imp = @decapode ExtendedOperators begin
    (dḞ1, dF1)::DualForm1{X}
    F1::Form1{X}
    dḞ1 == i₀(F1, dual_d₁{X}(dF1)) + dual_d₀{X}(i₁(F1, dF1))
  end
  lie1_imp = diag2dwd(Lie1Imp, in_vars = [:F1, :dF1], out_vars = [:dḞ1])

  I0Imp = @decapode ExtendedOperators begin
    F1::Form1{X}
    dF1::DualForm1{X}
    dF2::DualForm2{X}
    dF1 == ⋆₁{X}(∧₁₀{X}(F1, ⋆₀⁻¹{X}(dF2)))
  end
  i0_imp = diag2dwd(I0Imp, in_vars = [:F1, :dF2], out_vars = [:dF1])

  I1Imp = @decapode ExtendedOperators begin
    F1::Form1{X}
    dF1::DualForm1{X}
    dF0::DualForm0{X}
    dF0 == ⋆₂{X}(∧₁₁{X}(F1, ⋆₁⁻¹{X}(dF1)))
  end
  i1_imp = diag2dwd(I1Imp, in_vars = [:F1, :dF1], out_vars = [:dF0])

  δ₁Imp = @decapode ExtendedOperators begin
    F1::Form1{X}
    F0::Form0{X}
    F0 == ⋆₀⁻¹{X}(dual_d₁{X}(⋆₁{X}(F1)))
  end
  δ₁_imp = diag2dwd(δ₁Imp, in_vars = [:F1], out_vars = [:F0])

  δ₂Imp = @decapode ExtendedOperators begin
    F1::Form1{X}
    F2::Form2{X}
    F1 == ⋆₁⁻¹{X}(dual_d₀{X}(⋆₂{X}(F2)))
  end
  δ₂_imp = diag2dwd(δ₂Imp, in_vars = [:F2], out_vars = [:F1])

  Δ0Imp = @decapode ExtendedOperators begin
    (F0, Ḟ0)::Form0{X}
    Ḟ0 == δ₁{X}(d₀{X}(F0))
  end
  Δ0_imp = diag2dwd(Δ0Imp, in_vars = [:F0], out_vars = [:Ḟ0])

  Δ1Imp = @decapode ExtendedOperators begin
    (F1, Ḟ1)::Form1{X}
    Ḟ1 == δ₂{X}(d₁{X}(F1)) + d₀{X}(δ₁{X}(F1))
  end
  Δ1_imp = diag2dwd(Δ1Imp, in_vars = [:F1], out_vars = [:Ḟ1])

  Dict(:L₀ => lie0_imp, :i₀ => i0_imp, :L₁ => lie1_imp, :i₁ => i1_imp,
       :δ₁ => δ₁_imp, :δ₂ => δ₂_imp, :Δ₀ => Δ0_imp, :Δ₁ => Δ1_imp)
end
end


text:
using ITensors,
      Test

import ITensors.SmallString, ITensors.IntChar, ITensors.push
@testset "ctors" begin
    s = SmallString()
    @test length(s) == 0
    @test ITensors.isNull(s)
end

@testset "setindex" begin
    s = SmallString()
    @test ITensors.isNull(s)
    t = Base.setindex(s, IntChar(1), 1)
    @test length(t) == 1
    @test !ITensors.isNull(t)
end

@testset "push" begin
    s = SmallString()
    @test ITensors.isNull(s)
    t = push(s, IntChar(1))
    @test !ITensors.isNull(t)
end

@testset "comparison" begin
    s = SmallString()
    u = push(s, IntChar(1))
    t = push(s, IntChar(1))
    @test u == t
    t = push(s, IntChar(2))
    @test u < t
end


text:
using JuMP
using Ipopt
using CPLEX
using Gurobi
using PowerModels
include("prob.jl")
include("cut.jl")
include("misc.jl")

function master_problem(; file = "../data/nesta_case24_ieee_rts.m", k = 4, solver = IpoptSolver(), model_constructor = SOCWRPowerModel, cut_constructor = "AC", m_type="OPF")

    data = PowerModels.parse_file(file)
    buses = [ Int(bus["index"]) => bus for bus in data["bus"] ]
    num_buses = length(buses)
    buses = filter((i, bus) -> bus["bus_type"] != 4, buses)
    branches = [ Int(branch["index"]) => branch for branch in data["branch"] ]
    branches = filter((i, branch) -> branch["br_status"] == 1 && branch["f_bus"] in keys(buses) && branch["t_bus"] in keys(buses), branches)
    branch_indexes = collect(keys(branches))

    initial_feasible_solution = Dict{AbstractString, Any}()
    json_data_string = 0
    if file == "../data/nesta_case73_ieee_rts.m"
        json_data_string = readall(open("det_isolation_73.json"))
    elseif file == "../data/nesta_case24_ieee_rts.m"
        json_data_string = readall(open("det_isolation_24.json"))
    else
        println("no start solution")
    end

    if file == "../data/nesta_case73_ieee_rts__api.m"
        json_data_string = readall(open("det_isolation_73_api.json"))
    elseif file == "../data/nesta_case24_ieee_rts__api.m"
        json_data_string = readall(open("det_isolation_24_api.json"))
    else
        println("no start solution")
    end

    if file == "../data/nesta_case73_ieee_rts__sad.m"
        json_data_string = readall(open("det_isolation_73_sad.json"))
    elseif file == "../data/nesta_case24_ieee_rts__sad.m"
        json_data_string = readall(open("det_isolation_24_sad.json"))
    else
        println("no start solution")
    end

    if json_data_string != 0
        initial_feasible_solution = JSON.parse(json_data_string, dicttype = Dict{AbstractString,Any})
    end

    totalload = 0.0
    for (i,bus) in buses
        totalload += bus["pd"]
    end
    # println("total load: $totalload ...")

    if m_type == "OTS"
        @assert model_constructor == DCPPowerModel
    end
    if model_constructor == DCPPowerModel
        cut_constructor = "DC"
    end

    # master_problem = Model(solver=GurobiSolver(LogToConsole=0))
    master_problem = Model(solver=CplexSolver())
    @variable(master_problem, 0 <= eta <= 1e6)
    @variable(master_problem, x[i in branch_indexes], Bin)

    @constraint(master_problem, sum(x) == k)

    @objective(master_problem, Max, eta)

    # anonymous function
    x_val = i -> getvalue(x[i])

    zlb = 0
    zub = 1e10
    eps = 1e-6

    iteration = 0
    final_xvals = Dict{Any,Any}
    final_lines = Any[]
    current_lines = Any[]

    if length(initial_feasible_solution) != 0
        for value in initial_feasible_solution["attack"]
            if value["k"] == k
                current_lines = value["branch_ids"]
                break
            end
        end
        data = PowerModels.parse_file(file)
        for branch in data["branch"]
            for index in current_lines
                if Int(branch["index"]) == index
                    branch["br_status"] = 0
                end
            end
        end
        pm = model_constructor(data; solver = solver)
        if m_type != "OTS"
            post_pfls(pm)
        else
            post_ots(pm)
        end
        status, solve_time = solve(pm)
        result = PowerModels.build_solution(pm, status, solve_time)
        alpha = get_cut_coefficients(pm, cut_constructor = cut_constructor)
        zlb = result["objective"]
        @constraint(master_problem, eta <= result["objective"] + sum{alpha[i]*x[i], i in collect(keys(alpha))})
    end

    final_lines = current_lines
    tic()
    solve(master_problem)

    zub = getobjectivevalue(master_problem)

    println("ub: $zub, lb: $zlb")

    while (zub - zlb > 0.01*zlb)
        iteration += 1
        current_xvals = [ i => x_val(i) for i in branch_indexes]
        current_lines = collect(keys(filter( (i, xval) -> xval > 1e-4, current_xvals)))
        data = PowerModels.parse_file(file)
        for branch in data["branch"]
            for index in current_lines
                if Int(branch["index"]) == index
                    branch["br_status"] = 0
                end
            end
        end
        pm = model_constructor(data; solver = solver)
        if m_type != "OTS"
            post_pfls(pm)
        else
            post_ots(pm)
        end
        status, solve_time = solve(pm)
        result = PowerModels.build_solution(pm, status, solve_time)
        # println(result)
        alpha = get_cut_coefficients(pm, cut_constructor = cut_constructor)
        # println(getvalue(getvariable(pm.model, :line_z)))
        sub_objective = result["objective"]
        sub_lb = sub_objective
        if (sub_lb - zlb > eps)
            zlb = sub_lb
            final_xvals = current_xvals
            final_lines = current_lines
        end
        expr = zero(AffExpr)
        for i in branch_indexes
            if (current_xvals[i] <= eps)
                append!(expr, 1-x[i])
            else
                append!(expr, x[i])
            end
        end
        @constraint(master_problem, expr <= length(branch_indexes)-1)
        @constraint(master_problem, eta <= sub_objective + sum{alpha[i]*x[i], i in collect(keys(alpha))})
        solve(master_problem)
        zub = getobjectivevalue(master_problem)
        println("ub: $zub, lb: $zlb")
    end

    # independent_load_shed = get_independent_load_shed(file = file, solver = solver, model_constructor = model_constructor, lines = final_lines)
    println("###################################")
    println("k ... $k")
    println("solution = x* ... $(collect(final_lines))")
    println("model: $(model_constructor), cut_type: $(cut_constructor), model_type: $(m_type)")
    println("total number of iterations ... $iteration")
    println("solution cost ... $zlb")
    time_taken = toq()
    println("computation time ... $time_taken")
    #println("independent load shed values ... $independent_load_shed")
    #println("sum of n-1 load sheds ... $(sum(values(independent_load_shed)))")

    # solving the AC load shed model on the final set of lines
    data = PowerModels.parse_file(file)
    for branch in data["branch"]
        for index in final_lines
            if Int(branch["index"]) == index
                branch["br_status"] = 0
            end
        end
    end
    pm = ACPPowerModel(data; solver = IpoptSolver(print_level=0))
    post_pfls(pm)
    status, solve_time = solve(pm)
    result = PowerModels.build_solution(pm, status, solve_time)
    #println(result)
    println("full AC load shed ... $(result["objective"])")
    println("###################################")

    solution = Dict{AbstractString,Any}("k" => k, "lines" => collect(final_lines),
                                        "model" => model_constructor, "time" => time_taken,
                                        "cut" => cut_constructor, "load_shed" => zlb,
                                        "ac_load_shed" => result["objective"], "iterations" => iteration)
    push!(solution_vector, solution)

    return num_buses
end

solution_dict = Dict{AbstractString,Any}()
solution_vector = Any[]
args = Dict{Any,Any}()
println(ARGS[1])
args["file"] = string(ARGS[1])
args["k"] = parse(Int, ARGS[2])

# master_problem(file = args["file"], k = args["k"], solver = IpoptSolver(print_level=5), model_constructor = args["model"], cut_constructor = args["cut"])
num_buses = master_problem(file = args["file"], k = args["k"], solver = CplexSolver(), model_constructor = DCPPowerModel, cut_constructor = "DC")
num_buses = master_problem(file = args["file"], k = args["k"], solver = GurobiSolver(), model_constructor = SOCWRPowerModel, cut_constructor = "DC")
# num_buses = master_problem(file = args["file"], k = args["k"], solver = CplexSolver(), model_constructor = SOCWRPowerModel, cut_constructor = "AC")

# model_constructor = DCPPowerModel/DCPLLPowerModel/SOCWRPowerModel; cut_constructor = "DC"/"AC"

solution_dict["solution"] = solution_vector
json_string = JSON.json(solution_dict)

if endswith(args["file"], "api.m")
    write("./json_files/det_$(args["k"])_$(num_buses)_api.json", json_string)
elseif endswith(args["file"], "sad.m")
    write("./json_files/det_$(args["k"])_$(num_buses)_sad.json", json_string)
else
    write("./json_files/det_$(args["k"])_$(num_buses).json", json_string)
end

max_stars_repo_name: kaarthiksundar/acopf-n-k

max_issues_repo_name: kaarthiksundar/acopf-n-k

max_forks_repo_name: kaarthiksundar/acopf-n-k


text:
using Rotate
using Test

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


text:
using JuMP
using Cbc

function reopt_model(seconds::Float64)

    REopt = Model(with_optimizer(Cbc.Optimizer, logLevel=0, seconds=seconds))
	return REopt

end


text:
module LightXML
using XML2_jll

export
    # common
    name, free,

    # nodes
    AbstractXMLNode,
    XMLAttr, XMLAttrIter, XMLNode, XMLNodeIter, XMLElement, XMLElementIter,
    nodetype, value, content, attribute, has_attribute,
    is_elementnode, is_textnode, is_commentnode, is_cdatanode, is_blanknode, is_pinode,
    child_nodes, has_children, attributes, has_attributes, attributes_dict,
    child_elements, find_element, get_elements_by_tagname,

    new_element, add_child, new_child, new_textnode, add_text, add_cdata, add_pi,
    set_attribute, set_attributes, unlink, set_content,

    # document
    XMLDocument, version, encoding, compression, standalone, root,
    parse_file, parse_string, save_file, set_root, create_root

const Xchar = UInt8
const Xstr = Ptr{Xchar}

# opaque pointer type (do not dereference!) corresponding to xmlBufferPtr in C
struct xmlBuffer end
const Xptr = Ptr{xmlBuffer}

# pre-condition: p is not null
# (After tests, it seems that free in libc instead of xmlFree
#  should be used here.)
function _xcopystr(p::Xstr)
    r = unsafe_string(p)
    Libc.free(p)
    return r
end

include("errors.jl")

include("utils.jl")
include("nodes.jl")
include("document.jl")
include("cdata.jl")

end


text:
################################################################################

@testset "Backtrace" begin

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

  # declare model parameters
  params = HMMParams(
    penalty = 200,
    distance = euclDist,
    minimumFrequency = 20,
    verbosity = false,
  )

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

  # load small array into hidden Markov model
  include("data/smallArray.jl")

  # reset
  HiddenMarkovModelReaders.reset!(hmm)

  # feed frame
  for ι ∈ axes(x, 1)
    HiddenMarkovModelReaders.feed!(hmm, ι, x, params = params)
  end

  # backtrace
  HiddenMarkovModelReaders.backTrace(hmm)

  # testing
  @test hmm.data == [zeros(size(x, 2))]
  @test hmm.model == [[0.0, 11008.8486526594, 20778.203057470826, 30297.40010890268, 41255.56384672195, 53064.52909174562, 60668.05291217826, 70273.3913728092, 81813.54092914864, 90970.60619348944, 100636.00818678577]]
  @test hmm.traceback == [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

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

  # read file into hidden Markov model
  include("data/readfile.jl")

  # reset
  HiddenMarkovModelReaders.reset!(hmm)

  # feed frame
  for ι ∈ axes(x, 1)
    HiddenMarkovModelReaders.feed!(hmm, ι, x, params = params)
  end

  # backtrace
  HiddenMarkovModelReaders.backTrace(hmm)

  # testing
  @test hmm.data == [zeros(size(v, 2))]
  @test hmm.model == [[0.0, 11008.8486526594, 20778.203057470826, 30297.40010890268, 41255.56384672195, 53064.52909174562, 60668.05291217826, 70273.3913728092, 81813.54092914864, 90970.60619348944, 100636.00818678577, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0]]
  @test hmm.traceback == [-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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 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, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

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

end

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


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb04bd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 50
    LDA = NMAX÷2
    LDB = NMAX÷2
    LDC1 = NMAX÷2
    LDC2 = NMAX÷2
    LDDE = NMAX÷2
    LDF = NMAX÷2
    LDQ1 = NMAX
    LDQ2 = NMAX
    LDVW = NMAX÷2
    LDWORK = 2*NMAX*NMAX + max( 4*NMAX, 36 )
    LIWORK = max( NMAX + 12, 2*NMAX + 3 )
    A = Array{Float64,2}(undef, LDA,NMAX÷2)
    DE = Array{Float64,2}(undef, LDDE,NMAX÷2+1)
    C1 = Array{Float64,2}(undef, LDC1,NMAX÷2)
    VW = Array{Float64,2}(undef, LDVW,NMAX÷2+1)
    Q1 = Array{Float64,2}(undef, LDQ1,NMAX)
    Q2 = Array{Float64,2}(undef, LDQ2,NMAX)
    B = Array{Float64,2}(undef, LDB,NMAX÷2)
    F = Array{Float64,2}(undef, LDF,NMAX÷2)
    C2 = Array{Float64,2}(undef, LDC2,NMAX÷2)
    ALPHAR = Array{Float64,1}(undef, NMAX÷2)
    ALPHAI = Array{Float64,1}(undef, NMAX÷2)
    BETA = Array{Float64,1}(undef, NMAX÷2)
    IWORK = Array{BlasInt,1}(undef, LIWORK)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    JOB =  vs[1][1]
    COMPQ1 =  vs[2][1]
    COMPQ2 =  vs[3][1]
    N = parse(BlasInt, vs[4])
    if ( N<0 || N>NMAX )
        @error "Illegal N=$N"
    end

    M = N÷2

    vs = String[]
    _isz,_jsz = (M,M)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (M,M+1)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       DE[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (M,M)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       C1[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (M,M+1)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       VW[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    close(f)
# interp call 1

    INFO = SLICOT.mb04bd!(JOB, COMPQ1, COMPQ2, N, A, DE, C1, VW, Q1, Q2, B, F, C2, ALPHAR, ALPHAI, BETA, LIWORK, LDWORK)
    @test INFO == 0
    INFO == 0 || return

# interp output 1
    println(io, "A:")
    _nc = M
    _nr = M
    show(io, "text/plain", A[1:_nr,1:_nc])
    println(io)

# interp output 2
    println(io, "DE:")
    show(io, "text/plain", DE[1:M,2:M+1])
    println(io)

# interp output 3
    println(io, "B:")
    _nc = M
    _nr = M
    show(io, "text/plain", B[1:_nr,1:_nc])
    println(io)

# interp output 4
    println(io, "F:")
    _nc = M
    _nr = M
    show(io, "text/plain", UpperTriangular(F[1:_nr,1:_nc]))
    println(io)

# interp output 5
    println(io, "C1:")
    _nc = M
    _nr = M
    show(io, "text/plain", C1[1:_nr,1:_nc])
    println(io)

# interp output 6
    println(io, "C2:")
    _nc = M
    _nr = M
    show(io, "text/plain", C2[1:_nr,1:_nc])
    println(io)

# interp output 7
    println(io, "VW:")
    show(io, "text/plain", VW[1:M,2:M+1])
    println(io)

# interp output 8
    println(io, "ALPHAR:")
    _nr = M
    show(io, "text/plain", ALPHAR[1:_nr])
    println(io)

# interp output 9
    println(io, "ALPHAI:")
    _nr = M
    show(io, "text/plain", ALPHAI[1:_nr])
    println(io)

# interp output 10
    println(io, "BETA:")
    _nr = M
    show(io, "text/plain", BETA[1:_nr])
    println(io)

if ( !LSAME( COMPQ1, 'N' ) )
# interp output 11
    println(io, "Q1:")
    _nc = N
    _nr = N
    show(io, "text/plain", Q1[1:_nr,1:_nc])
    println(io)

end # if
if ( !LSAME( COMPQ2, 'N' ) )
# interp output 12
    println(io, "Q2:")
    _nc = N
    _nr = N
    show(io, "text/plain", Q2[1:_nr,1:_nc])
    println(io)
end
end # run_mb04bd()


text:
"""
    bellman_ford(graph::Graph, source::Int)::Tuple{Vector{AbstractFloat}, Vector{Int}}

The shortest path distances using `Bellman Ford` algorithm
"""
function bellman_ford(graph::Graph, source::Int)::Tuple{Vector{AbstractFloat}, Vector{Int}}
    node_size = size(graph, 1)
    dist = ones(Float64, node_size) * Inf
    dist[source] = 0.0
    pred = ones(Int64, node_size) * -1
    for k ∈ 1:node_size-1
        for node_ind ∈ 1:node_size
            curr_node = graph.nodes[node_ind]
            for (neigh, edge) ∈ adj(graph, curr_node)
                t_ind = neigh.index
                weight = edge.cost
                if dist[t_ind] > dist[node_ind] + weight
                    dist[t_ind] = dist[node_ind] + weight
                    pred[t_ind] = node_ind
                end
            end
        end
    end
    return dist, pred
end


text:
## Headers Handlers

function recv_stream_headers(connection::HTTPConnection, frame::HeadersFrame)
    stream_identifier = frame.stream_identifier
    stream = get_stream(connection, stream_identifier)

    if frame.is_priority
        handle_priority!(connection, stream_identifier, frame.exclusive.value,
                         frame.dependent_stream_identifier.value,
                         frame.weight.value)
    end

    block = copy(frame.fragment)
    headers = HPack.decode(connection.dynamic_table, IOBuffer(block))
    put!(connection.channel_evt, EvtRecvHeaders(stream_identifier, headers, frame.is_end_stream))
end

function send_stream_headers(connection::HTTPConnection, act::ActSendHeaders)
    stream_identifier = act.stream_identifier
    stream = get_stream(connection, stream_identifier)

    block = HPack.encode(connection.dynamic_table, act.headers; huffman=false)
    is_end_stream = act.is_end_stream

    # We don't use padding in this implementation
    if connection.settings.max_frame_size < (length(block) + 6)
        splitLength = connection.settings.max_frame_size - 6
        header = HeadersFrame(is_end_stream, false, false, stream_identifier, Nullable{Bool}(),
                              Nullable{UInt32}(), Nullable{UInt8}(), getindex(block, 1:splitLength))
        put!(connection.channel_act_raw, header)

        curPos = splitLength + 1
        while curPos <= length(block)
            endPos = max(curPos + splitLength - 1, length(block))
            continuation = ContinuationFrame(endPos >= length(block), stream_identifier, getindex(block, curPos:endPos))
            put!(connection.channel_act_raw, continuation)
            curPos = endPos + 1
        end
    else
        frame = HeadersFrame(is_end_stream, true, false, stream_identifier, Nullable{Bool}(),
                             Nullable{UInt32}(), Nullable{UInt8}(), block)
        put!(connection.channel_act_raw, frame)
    end

    return frame
end

## Data Handlers

function recv_stream_data(connection::HTTPConnection, frame::DataFrame)
    stream = get_stream(connection, frame.stream_identifier)

    data = copy(frame.data)
    put!(connection.channel_evt, EvtRecvData(frame.stream_identifier, data, frame.is_end_stream))
end

function send_stream_data(connection::HTTPConnection, act::ActSendData)
    stream = get_stream(connection, act.stream_identifier)

    is_end_stream = act.is_end_stream

    if connection.settings.max_frame_size < (length(act.data) + 1)
        splitLength = connection.settings.max_frame_size - 1
        curPos = 1

        while curPos <= length(act.data)
            endPos = max(curPos + splitLength - 1, length(act.data))
            frame = DataFrame(act.stream_identifier, is_end_stream, getindex(act.data, curPos:endPos))
            put!(connection.channel_act_raw, frame)
            curPos = endPos + 1
        end
    else
        frame = DataFrame(act.stream_identifier, is_end_stream, act.data)
        put!(connection.channel_act_raw, frame)
    end

    return frame
end

## Priority Handlers

function recv_stream_priority(connection::HTTPConnection, priority::PriorityFrame)
    handle_priority!(connection, priority.stream_identifier, priority.exclusive,
                     priority.dependent_stream_identifier,
                     priority.weight)
end

## Rst Stream Handlers

function recv_stream_rst_stream(connection::HTTPConnection, rstStream::RstStreamFrame)
    ## Nothing to do
end

## Push Promise Handlers

function recv_stream_push_promise(connection::HTTPConnection, frame::PushPromiseFrame)
    stream = get_stream(connection, frame.stream_identifier)

    block = copy(frame.fragment)
    headers = HPack.decode(connection.dynamic_table, IOBuffer(block))

    put!(connection.channel_evt, EvtPromise(frame.stream_identifier, frame.promised_stream_identifier,
                                            headers))
end

function send_stream_push_promise(connection::HTTPConnection, act::ActPromise)
    stream_identifier = act.stream_identifier
    stream = get_stream(connection, act.stream_identifier)

    block = HPack.encode(connection.dynamic_table, act.headers; huffman=false)

    # We don't use padding in this implementation
    if connection.settings.max_frame_size < (length(block) + 5)
        splitLength = connection.settings.max_frame_size - 5
        header = PushPromiseFrame(false, stream_identifier, act.promised_stream_identifier, getindex(block, 1:splitLength))
        put!(connection.channel_act_raw, header)

        curPos = splitLength + 1
        while curPos <= length(block)
            endPos = max(curPos + splitLength - 1, length(block))
            continuation = ContinuationFrame(endPos >= length(block), stream_identifier, getindex(block, curPos:endPos))
            put!(connection.channel_act_raw, continuation)
            curPos = endPos + 1
        end
    else
        frame = PushPromiseFrame(true, act.stream_identifier, act.promised_stream_identifier, block)
        put!(connection.channel_act_raw, frame)
    end

    return frame
end

## Window Update Handlers

function recv_stream_window_update(connection::HTTPConnection, frame::WindowUpdateFrame)
    stream = get_stream(connection, frame.stream_identifier)

    stream.window_size += frame.window_size_increment
end


text:
module TestSession

using Matte

title = "Test Session"

function ui()
    sidebar_layout(
        main_panel(
        p(text_output("counter"))
    ),
    side_panel(
        h1("hi there"),
        button("button", "Push me")
    ))
end

function register_session_vars()
    Dict(
        :counter => 0
    )
end

module Server

import DataFrames

function counter(button, session)
    if button
        session.counter += 1
        session.counter
    else
        session.counter
    end
end

end

end

run_app(TestSession, async = true)

HTTP.WebSockets.open("ws://127.0.0.1:8000") do ws
    payload = Dict(
        :channel => "matte",
        :message => "api",
        :payload => Dict(
            :id => "counter",
            :session_id => "session1",
            :input => Dict(
                :button => true
            )
        )
    )
    expected = Dict(
        "id" => "counter",
        "value" => 1
    )
    write(ws, JSON.json(payload))
    @test JSON.parse(String(readavailable(ws))) == expected
end

HTTP.WebSockets.open("ws://127.0.0.1:8000") do ws
    payload = Dict(
        :channel => "matte",
        :message => "api",
        :payload => Dict(
            :id => "counter",
            :session_id => "session1",
            :input => Dict(
                :button => true
            )
        )
    )
    expected = Dict(
        "id" => "counter",
        "value" => 2
    )
    write(ws, JSON.json(payload))
    @test JSON.parse(String(readavailable(ws))) == expected
end

HTTP.WebSockets.open("ws://127.0.0.1:8000") do ws
    payload = Dict(
        :channel => "matte",
        :message => "api",
        :payload => Dict(
            :id => "counter",
            :session_id => "session2",
            :input => Dict(
                :button => true
            )
        )
    )
    expected = Dict(
        "id" => "counter",
        "value" => 1
    )
    write(ws, JSON.json(payload))
    @test JSON.parse(String(readavailable(ws))) == expected
end

stop_app()


text:
using Plots
using MCMCChain

n_iter = 500
n_name = 3
n_chain = 2

val = randn(n_iter, n_name, n_chain) .+ [1, 2, 3]'
val = hcat(val, rand(1:2, n_iter, 1, n_chain))

chn = Chains(val)

# plotting singe plotting types
ps_trace = plot(chn, :trace)
ps_mean = plot(chn, :mean)
ps_density = plot(chn, :density)
ps_autocor = plot(chn, :autocor)
#ps_contour = plot(chn, :contour)
ps_hist = plot(chn, :histogram)
ps_mixed = plot(chn, :mixeddensity)

# plotting combinations
ps_trace_mean = plot(chn, [:trace, :mean])
ps_mixed_auto = plot(chn, [:mixeddensity, :autocor])


text:
printstyled("  demean, detrend, unscale\n", color=:light_green)

S = randSeisData(2, s=1.0)[2:2]
fs = 100.0
nx = 14400
t = floor(Int64, time()-60.0)*sμ

S.fs[1] = 100.0
S.t[1] = [1 t; nx 0]
x = (rand(nx) .- 0.5)
broadcast!(-, x, x, mean(x))
S.x[1] = x
S₀ = deepcopy(S)

# Test 1: de-mean
μ = 6.0
S.x[1] .+= μ
T = demean(S)

@test maximum(S.x[1] .- T.x[1]) ≤ 6.01
@test minimum(S.x[1] .- T.x[1]) ≥ 5.99
@test maximum(abs.(T.x[1]-S₀.x[1])) < 0.01
@test abs(mean(T.x[1])) < 0.01

# Test 2: de-trend
m = 0.01
demean!(S)
S.x[1] += (m.*(1.0:1.0:Float64(nx)))
U = detrend(S)
@test maximum(abs.(U.x[1]-S₀.x[1])) < 0.1
@test maximum(abs.(U.x[1]-T.x[1])) < 0.1
@test abs(mean(U.x[1])) < 1.0e-8

# de-mean and de-trend with NaNs
nx = length(S.x[1])
for i = 1:3
  j = randperm(nx)[1:rand(1:div(nx,2))]
  V = deepcopy(S)
  V.x[1][j] .= NaN
  demean!(V)
  @test length(j) == length(findall(isnan.(V.x[1])))

  k = randperm(nx)[1:rand(1:div(nx,2))]
  W = deepcopy(S)
  W.x[1][k] .= NaN
  detrend!(W)
  @test length(k) == length(findall(isnan.(W.x[1])))

  C = randSeisChannel()
  j = randperm(length(C.x))[1:rand(1:div(length(C.x),2))]
  C.x[j] .= NaN
  D = deepcopy(C)
  demean!(C)
  detrend!(D)
  @test length(j) == length(findall(isnan.(C.x))) == length(findall(isnan.(D.x)))
end

# Safe demean, detrend
C = randSeisChannel()
D = detrend(C)
S = randSeisData()
T = detrend(S)
# V = randSeisEvent()
# W = detrend(V)

C = randSeisChannel()
D = demean(C)
S = randSeisData()
T = demean(S)
# V = randSeisEvent()
# W = demean(V)

demean!(C)
detrend!(C)
# demean!(V)
# detrend!(V)

# Test in-place unscaling
C = randSeisChannel(s=true)
S = randSeisData()
# V = randSeisEvent()
unscale!(C)
unscale!(S, irr=true)
# unscale!(V, irr=true)

# Test for out-of-place unscaling
C = randSeisChannel()
S = randSeisData()
# V = randSeisEvent()
D = unscale(C)
T = unscale(S)
# W = unscale(V)

# tests for channel ranges
S = randSeisData(10, s=1.0)
U = deepcopy(S)
demean!(S, chans=1:5)
for i = 1:S.n
  @test (S[i] == U[i]) == (i < 6 ? false : true)
end

S = randSeisData(10, s=1.0)
U = deepcopy(S)
detrend!(S, chans=1:5)
for i = 1:S.n
  @test (S[i] == U[i]) == (i < 6 ? false : true)
end


text:
# low-level interface, meaning works with vectors and matrices, not data frames and estimation objects

"Structure to store test results"
struct BootTestResult{T}
  stat::T; stattype::String
  p::T; padj::T
  reps::Int64; repsfeas::Int64
  nbootclust::Int64
  dof::Int64; dof_r::T
  plot::Union{Nothing, NamedTuple{(:X, :p), Tuple{Tuple{Vararg{Vector{T}, N} where N},Vector{T}}}}
  peak::Union{Nothing, NamedTuple{(:X, :p), Tuple{Vector{T}, T}}}
  CI::Union{Nothing, Matrix{T}}
  dist::Matrix{T}
  b::Vector{T}
  V::Matrix{T}
  auxweights::Union{Nothing,Matrix{T}}
  # M::StrBootTest
end

"Return test statistic"
teststat(o::BootTestResult) = o.stat

"Return numerator of test statistic"
statnumer(o::BootTestResult) = o.b

"Return denominator of test statistic"
statvar(o::BootTestResult) = o.V

"""Return type of test statistic subject: "t", "z", "F", or "χ²" """
stattype(o::BootTestResult) = o.stattype

"Return p value"
p(o::BootTestResult) = o.p

"Returnp p value after multiple-hypothesis adjustment, if any"
padj(o::BootTestResult) = o.padj

"Return requested number of replications"
reps(o::BootTestResult) = o.reps

"Return actual number of replications, subject to enumeration of Rademacher draws"
repsfeas(o::BootTestResult) = o.repsfeas

"Return number of bootstrapping clusters in test"
nbootclust(o::BootTestResult) = o.nbootclust

"Return degrees of freedom test"
dof(o::BootTestResult) = o.dof

"Return residual degrees of freedom test"
dof_r(o::BootTestResult) = o.dof_r

"""
Return data for confidence plot of test.
Return value is a 2-tuple with named entries `X` and `p` holding
the confidence sampling locations and p values respectively. `X` is in turn
a 1- or 2-tuple of vectors of sampling coordinates for each
dimension of the tested hypothesis.
"""
plotpoints(o::BootTestResult) = o.plot

"Return parameter value with peak p value in test"
peak(o::BootTestResult) = o.peak

"Return confidence interval matrix from test, one row per disjoint piece"
CI(o::BootTestResult) = o.CI

"Return bootstrap distribution of statistic or statistic numerator in bootstrap test"
dist(o::BootTestResult) = o.dist

"Return auxilliary weight matrix for wild bootstrap"
auxweights(o::BootTestResult) = o.auxweights

using Printf
function Base.show(io::IO, o::BootTestResult{T}) where T
	print(io, "WildBootTests.BootTestResult{$T}\n\n")
	Printf.@printf(io, "%s = %5.3f\n", stattype(o)*repeat(' ',2-length(stattype(o))), teststat(o))
	Printf.@printf(io, "p  = %5.3f\n", p(o))
	isdefined(o, :CI) && !isnothing(o.CI) && length(o.CI)>0 && print(io, "CI = $(CI(o))\n")
end

# single entry point with arguments already converted to standardized types, to allow a small set of precompile() calls
function __wildboottest(
	R::Matrix{T},
	r::Vector{T};
	resp::Vector{T},
	predexog::Matrix{T},
	predendog::Matrix{T},
	inst::Matrix{T},
	R1::Matrix{T},
	r1::Vector{T},
	clustid::Matrix{Int64},
	nbootclustvar::Int8,
	nerrclustvar::Int8,
	issorted::Bool,
	hetrobust::Bool,
	feid::Vector{Int64},
	fedfadj::Bool,
	obswt::Vector{T},
	fweights::Bool,
	maxmatsize::Float16,
	ptype::PType,
	bootstrapc::Bool,
	LIML::Bool,
	Fuller::T,
	kappa::T,
	ARubin::Bool,
	small::Bool,
	scorebs::Bool,
	reps::Int64,
	imposenull::Bool,
	auxwttype::AuxWtType,
	rng::AbstractRNG,
	level::T,
	rtol::T,
	madjtype::MAdjType,
	NH0::Int16,
	ML::Bool,
	scores::Matrix{T},
	beta::Vector{T},
	A::Symmetric{T,Matrix{T}},
	gridmin::Vector{T},
	gridmax::Vector{T},
	gridpoints::Vector{Float32},
	diststat::DistStatType,
	getCI::Bool,
	getplot::Bool,
	getauxweights::Bool,
	turbo::Bool) where T

	M = StrBootTest{T}(R, r, R1, r1, resp, predexog, predendog, inst, obswt, fweights, LIML, Fuller, kappa, ARubin,
	                   reps, auxwttype, rng, maxmatsize, ptype, imposenull, scorebs, !bootstrapc, clustid, nbootclustvar, nerrclustvar, issorted, hetrobust, small,
	                   feid, fedfadj, level, rtol, madjtype, NH0, ML, beta, A, scores, getplot,
	                   gridmin, gridmax, gridpoints, turbo)

	if getplot || (level<1 && getCI)
		plot!(M)
		plot = getplot & isdefined(M, :plotX) ? (X=Tuple(M.plotX), p=M.plotY) : nothing
		peak = M.peak
		CI = level<1 & getCI ? M.CI : nothing
	else
		CI = plot = peak = nothing
	end

	padj = getp(M)  # trigger central (re)computation

	BootTestResult{T}(getstat(M),
	                  isone(nrows(R)) ? (small ? "t" : "z") : (small ? "F" : "χ²"),
	                  M.p, padj, M.B, M.BFeas, M.N✻, M.dof, M.dof_r, plot, peak, CI,
	                  getdist(M, diststat),
	                  getb(M), getV(M),
	                  getauxweights && reps>0 ? getv(M) : nothing #=, M=#)
end

vecconvert(T::DataType, X) = Vector(isa(X, AbstractArray) ? vec(    eltype(X)==T ? X : T.(X)) : X)
matconvert(T::DataType, X) = Matrix(isa(X, AbstractArray) ? reshape(eltype(X)==T ? X : T.(X), size(X,1), size(X,2)) : X)

function _wildboottest(T::DataType,
					  R::AbstractVecOrMat,
						r::AbstractVecOrMat;
					  resp::AbstractVector{<:Real},
					  predexog::AbstractVecOrMat{<:Real}=zeros(T,0,0),
					  predendog::AbstractVecOrMat{<:Real}=zeros(T,0,0),
					  inst::AbstractVecOrMat{<:Real}=zeros(T,0,0),
					  R1::AbstractVecOrMat=zeros(T,0,0),
						r1::AbstractVector=zeros(T,0),
					  clustid::AbstractVecOrMat{<:Integer}=zeros(Int,0,0),  # bootstrap-only clust vars, then boot&err clust vars, then err-only clust vars
					  nbootclustvar::Integer=1,
					  nerrclustvar::Integer=nbootclustvar,
						issorted::Bool=false,
					  hetrobust::Bool=true,
					  feid::AbstractVector{<:Integer}=Int8[],
					  fedfadj::Bool=true,
					  obswt::AbstractVector{<:Real}=T[],
					  fweights::Bool=false,
					  maxmatsize::Number=0,
					  ptype::PType=symmetric,
					  bootstrapc::Bool=false,
					  LIML::Bool=false,
					  Fuller::Number=0,
					  kappa::Number=NaN,
					  ARubin::Bool=false,
					  small::Bool=true,
					  scorebs::Bool=false,
					  reps::Integer=999,
					  imposenull::Bool=true,
					  auxwttype::AuxWtType=rademacher,
					  rng::AbstractRNG=MersenneTwister(),
					  level::Number=.95,
					  rtol::Number=1e-6,
					  madjtype::MAdjType=nomadj,
					  NH0::Integer=1,
					  ML::Bool=false,
					  scores::AbstractVecOrMat=Matrix{Float32}(undef,0,0),
					  beta::AbstractVector=T[],
					  A::AbstractMatrix=zeros(T,0,0),
					  gridmin::Union{Vector{S},Vector{Union{S,Missing}}} where S<:Number = T[],
					  gridmax::Union{Vector{S},Vector{Union{S,Missing}}} where S<:Number = T[],
					  gridpoints::Union{Vector{S},Vector{Union{S,Missing}}} where S<:Integer = Int32[],
					  diststat::DistStatType=nodist,
					  getCI::Bool=true,
					  getplot::Bool=getCI,
					  getauxweights::Bool=false,
						turbo::Bool=false)

	nrows(R)>2 && (getplot = getCI = false)

  @assert length(predexog)==0 || nrows(predexog)==nrows(resp) "All data vectors/matrices must have same height"
  @assert length(predendog)==0 || nrows(predendog)==nrows(resp) "All data vectors/matrices must have same height"
  @assert length(inst)==0 || nrows(inst)==nrows(resp) "All data vectors/matrices must have same height"
  @assert length(feid)==0 || nrows(feid)==nrows(resp) "feid vector must have same height as data matrices"
  @assert length(clustid)==0 || nrows(clustid)==nrows(resp) "clustid must have same height as data matrices"
  @assert nrows(obswt)==0 || nrows(obswt)==nrows(resp) "obswt must have same height as data matrices"
  @assert nrows(R)==nrows(r) "R and r must have same height"
  @assert ncols(R)==ncols(predexog)+ncols(predendog) && isone(ncols(r)) "Wrong number of columns in null specification"
  @assert nrows(R1)==nrows(r1) "R₁ and r₁ must have same height"
  @assert length(R1)==0 || ncols(R1)==ncols(predexog)+ncols(predendog) "Wrong number of columns in model constraint specification"
	@assert ncols(r)==1 "r should have one column"
	@assert length(R1)==0 || ncols(r1)==1 "r1 should have one column"
  @assert nbootclustvar ≤ ncols(clustid) "nbootclustvar > width of clustid"
  @assert nerrclustvar ≤ ncols(clustid) "nerrclustvar > width of clustid"
  @assert reps ≥ 0 "reps < 0"
  @assert level ≥ 0. && level≤1. "level must be in the range [0,1]"
  @assert rtol > 0. "rtol ≤ 0"
  @assert NH0 > 0 "NH0 ≤ 0"
	if getplot || getCI
		@assert iszero(length(gridmin   )) || length(gridmin   )==nrows(R) "Length of gridmin doesn't match number of hypotheses being jointly tested"
		@assert iszero(length(gridmax   )) || length(gridmax   )==nrows(R) "Length of gridmax doesn't match number of hypotheses being jointly tested"
		@assert iszero(length(gridpoints)) || length(gridpoints)==nrows(R) "Length of gridpoints doesn't match number of hypotheses being jointly tested"
	end

	_gridmin = Vector{T}(undef, length(gridmin))
	_gridmax = Vector{T}(undef, length(gridmax))
	_gridpoints = Vector{Float32}(undef, length(gridpoints))
	for i ∈ 1:length(gridmin)  # cumbersome loops because map() and list comprehensions mess up type inference(?!)
		_gridmin[i] = T(ismissing(gridmin[i]) ? NaN : gridmin[i])
	end
	for i ∈ 1:length(gridmax)
		_gridmax[i] = T(ismissing(gridmax[i]) ? NaN : gridmax[i])
	end
	for i ∈ 1:length(gridpoints)
		_gridpoints[i] = T(ismissing(gridpoints[i]) ? NaN : gridpoints[i])
	end

	__wildboottest(
		matconvert(T,R),
		vecconvert(T,r);
		resp=vecconvert(T,resp),
		predexog=matconvert(T,predexog),
		predendog=matconvert(T,predendog),
		inst=matconvert(T,inst),
		R1=matconvert(T,R1),
		r1=vecconvert(T,r1),
		clustid=matconvert(Int64,clustid),
		nbootclustvar=Int8(nbootclustvar),
		nerrclustvar=Int8(nerrclustvar),
		issorted,
		hetrobust,
		feid=vecconvert(Int64,feid),
		fedfadj,
		obswt=vecconvert(T,obswt),
		fweights,
		maxmatsize=Float16(maxmatsize),
		ptype,
		bootstrapc,
		LIML,
		Fuller=T(Fuller),
		kappa=T(kappa),
		ARubin,
		small,
		scorebs,
		reps=Int64(reps),
		imposenull,
		auxwttype,
		rng,
		level=T(level),
		rtol=T(rtol),
		madjtype,
		NH0=Int16(NH0),
		ML,
		scores=matconvert(T,scores),
		beta=vecconvert(T,beta),
		A=Symmetric(matconvert(T,A)),
		gridmin=_gridmin,
		gridmax=_gridmax,
		gridpoints=_gridpoints,
		diststat,
		getCI,
		getplot,
		getauxweights,
		turbo)
end

_wildboottest(T::DataType, R, r::Number; kwargs...) = _wildboottest(T, R, [r]; kwargs...)
_wildboottest(T::DataType, R::UniformScaling{Bool}, r; kwargs...) = _wildboottest(T, diagm(fill(T(R.λ),nrows(r))), r; kwargs...)
_wildboottest(R::Union{UniformScaling{Bool},AbstractVecOrMat}, r::Union{Number,AbstractVecOrMat}; kwargs...) = _wildboottest(Float32, R, r; kwargs...)

"""
wildboottest([T::DataType=Float32,] R::AbstractMatrix, r::AbstractVector;
             resp, <optional keyword arguments>) -> WildBootTests.BootTestResult

Function to perform wild-bootstrap-based hypothesis test

# Positional arguments
* `T::DataType`: data type for inputs, results, and computations: Float32 (default) or Float64
* `R::AbstractMatrix` and `r::AbstractVector`: required matrix and vector expressing the null Rβ=r; see notes below

# Required keyword argument
* `resp::AbstractVector`: response/dependent variable (y or y₁ in Roodman et al. (2019))

# Optional keyword arguments
* `predexog::AbstractVecOrMat`: exogenous predictors, including constant term, if any (X/X₁)
* `predendog::AbstractVecOrMat`: endogenous predictors (Y₂)
* `inst::AbstractVecOrMat`: instruments (X₂)
* `R1::AbstractMatrix` and `r1::AbstractVector`: model constraints; same format as for `R` and `r`
* `clustid::AbstractVecOrMat{<:Integer}`: data vector/matrix of error and bootstrapping cluster identifiers; see Notes
* `nbootclustvar::Integer=1`: number of bootstrap-clustering variables
* `nerrclustvar::Integer=nbootclustvar`: number of error-clustering variables
* `hetrobust::Bool=true`: true unless errors are treated as iid
* `feid::AbstractVector{<:Integer}`: data vector for fixed effect group identifier
* `fedfadj::Bool=true`: true if small-sample adjustment should reflect number of fixed effects (if any)
* `obswt::AbstractVector=[]`: observation weight vector; default is equal weighting
* `fweights::Bool=false`: true for frequency weights
* `maxmatsize::Number`: maximum size of auxilliary weight matrix (v), in gigabytes
* `ptype::PType=symmetric`: p value type (`symmetric`, `equaltail`, `lower`, `upper`)
* `bootstrapc::Bool=false`: true to request bootstrap-c instead of bootstrap-t
* `LIML::Bool=false`: true for LIML or Fuller LIML
* `Fuller::Number`: Fuller LIML factor
* `kappa::Number`: fixed κ for _k_-class estimation
* `ARubin::Bool=false`: true for Anderson-Rubin test
* `small::Bool=true`: true for small-sample corrections
* `scorebs::Bool=false`: true for score bootstrap instead of wild bootstrap
* `reps::Integer=999`: number of bootstrap replications; `reps` = 0 requests classical Rao (or Wald) test if `imposenull` = `true` (or `false`)
* `imposenull::Bool=true`: true to impose null
* `auxwttype::AuxWtType=rademacher`: auxilliary weight type (`rademacher`, `mammen`, `webb`, `normal`, `gamma`)
* `rng::AbstractRNG=MersenneTwister()`: randon number generator
* `level::Number=.95`: significance level (0-1)
* `rtol::Number=1e-6`: tolerance for CI bound determination
* `madjtype::MAdjType=nomadj`: multiple hypothesis adjustment (`nomadj`, `bonferroni`, `sidak`)
* `NH0::Integer=1`: number of hypotheses tested, including one being tested now
* `ML::Bool=false`: true for (nonlinear) ML estimation
* `scores::AbstractVecOrMat`: for ML, pre-computed scores
* `beta::AbstractVector`: for ML, parameter estimates
* `A::AbstractMatrix`: for ML, covariance estimates
* `gridmin`: vector of graph lower bounds; max length 2, `missing`/`NaN` entries ask wildboottest() to choose
* `gridmax`: vector of graph upper bounds; `missing`/`NaN` entries ask wildboottest() to choose
* `gridpoints`: vector of number of sampling points; `missing`/`NaN` entries ask wildboottest() to choose
* `diststat::DistStatType=nodiststat`: t to save bootstrap distribution of t/z/F/χ² statistics; numer to save numerators thereof
* `getCI::Bool=true`: whether to return CI
* `getplot::Bool=getCI`: whether to generate plot data
* `getauxweights::Bool=false`: whether to save auxilliary weight matrix (v)
* `turbo::Bool=false`: whether to exploit acceleration of the LoopVectorization package: slower on first use in a session, faster after
* `issorted::Bool=false`: time-saving flag that all data matrices are already sorted by variables 2, then 3, then 1, as numbered below

# Notes
`T`, `ptype`, `auxwttype`, `madjtype`, and `diststat` may also be strings. Examples: `"Float32"` and `"webb"`.

The columns of `R` in the statement of the null should correspond to those of the matrix [`predexog` `predendog`],
where `predendog` is non-empty only in regressions with instruments.

Order the columns of `clustid` this way:
1. Variables only used to define bootstrapping clusters, as in the subcluster bootstrap.
2. Variables used to define both bootstrapping and error clusters.
3. Variables only used to define error clusters.
`nbootclustvar` is then the number of columns of type 1 or 2; `nerrclustvar` is the number of columns of type 2 or 3. Typically `clustid` is a single column of type 2.

`wildboottest()` does not handle missing data values: all data and identifier matrices must be restricted to the estimation sample.

"""
wildboottest(   R, r; args...) = _wildboottest(                                                  R, r; Dict(a.first => (isa(a.second, AbstractString) ? eval(Meta.parse(a.second)) : a.second) for a ∈ args)...)
wildboottest(T, R, r; args...) = _wildboottest(isa(T, AbstractString) ? eval(Meta.parse(T)) : T, R, r; Dict(a.first => (isa(a.second, AbstractString) ? eval(Meta.parse(a.second)) : a.second) for a ∈ args)...)


text:
## --- Geolcont

    continentcolors = parse.(Color, ["#333399","#0066CC","#06A9C1","#66CC66","#FFCC33","#FFFF00","#FFFFFF"])
    export continentcolors

    continents = ["Africa","Eurasia","North America","South America","Australia","Antarctica","NA"]
    export continents

    # Find which continent a sample originates from
    function find_geolcont(lat,lon)
        # Interpret user input
        if length(lat) != length(lon)
            error("lat and lon must be equal length\n")
        end

        # Construct file path
        filedir = joinpath(resourcepath,"maps")
        filepath = joinpath(filedir,"geolcontwshelf.png")

        # Download png map from Google Cloud if necessary
        if ~isfile(filepath)
            print("Downloading map to $filedir\n")
            run(`mkdir -p $filedir`)
            download("https://storage.googleapis.com/statgeochem/geolcontwshelf.png", filepath)
        end

        img = load(filepath)

        ind = fill(7,size(img))
        for i=1:6
         ind[img .== continentcolors[i]] .= i
        end

        # Create and fill output vector
        contindex = Array{Int}(undef,size(lat))
        for i=1:length(lat)
         if isnan(lat[i]) || isnan(lon[i]) || lat[i]>90 || lat[i]<-90 || lon[i]>180 || lon[i]<-180
             # Result is unknown if either input is NaN or out of bounds
             contindex[i] = 7
         else
             # Convert latitude and longitude into indicies of the elevation map array
             # Note that STRTM15 plus has N+1 columns where N = 360*sf
             row = 1 + trunc(Int,(90-lat[i])*512/180)
             col = 1 + trunc(Int,(180+lon[i])*512/180)
             # Find result by indexing
             contindex[i] = ind[row,col]
         end
        end

        return contindex
    end
    export find_geolcont

## --- ETOPO1 (1 arc minute topography)

    # Read etopo file from HDF5 storage, downloading from cloud if necessary
    function get_etopo(varname="")
        # Available variable names: "elevation", "y_lat_cntr", "x_lon_cntr",
        # "cellsize", "scalefactor", and "reference". Units are meters of
        # elevation and decimal degrees of latitude and longitude

        # Construct file path
        filedir = joinpath(resourcepath,"etopo")
        filepath = joinpath(filedir,"etopo1.h5")

        # Download HDF5 file from Google Cloud if necessary
        if ~isfile(filepath)
            print("Downloading etopo1.h5 from google cloud storage to $filedir\n")
            run(`mkdir -p $filedir`)
            download("https://storage.googleapis.com/statgeochem/etopo1.references.txt", joinpath(filedir,"etopo1.references.txt"))
            download("https://storage.googleapis.com/statgeochem/etopo1.h5", filepath)
        end

        # Read and return the file
        return h5read(filepath, "vars/"*varname)
    end
    export get_etopo

    # Find the elevation of points at position (lat,lon) on the surface of the
    # Earth, using the ETOPO elevation model.
    find_etopoelev(lat,lon) = find_etopoelev(get_etopo(),lat,lon)
    find_etopoelev(etopo::Dict, lat, lon) = find_etopoelev(etopo["elevation"], lat, lon)
    function find_etopoelev(etopo::AbstractArray, lat, lon)
        # Interpret user input
        length(lat) != length(lon) && error("lat and lon must be of equal length\n")

        # Scale factor (cells per degree) = 60 = arc minutes in an arc degree
        sf = 60
        maxrow = 180 * sf
        maxcol = 360 * sf

        # Create and fill output vector
        result = Array{Float64}(undef,size(lat))
        for i=1:length(lat)
            if isnan(lat[i]) || isnan(lon[i]) || lat[i]>90 || lat[i]<-90 || lon[i]>180 || lon[i]<-180
                # Result is NaN if either input is NaN or out of bounds
                result[i] = NaN
            else
                # Convert latitude and longitude into indicies of the elevation map array
                row = 1 + trunc(Int,(90+lat[i])*sf)
                if row == (maxrow+1) # Edge case
                    row = maxrow
                end

                col = 1 + trunc(Int,(180+lon[i])*sf)
                if col == (maxcol+1) # Edge case
                    col = maxcol
                end

                # Find result by indexing
                result[i] = etopo[row,col]
            end
        end

        return result
    end
    export find_etopoelev

## --- SRTM15_PLUS (15 arc second topography)

    # Read srtm15plus file from HDF5 storage, downloading from cloud if necessary
    function get_srtm15plus(varname="")
        # Available variable names: "elevation", "y_lat_cntr", "x_lon_cntr",
        # "nanval", "cellsize", "scalefactor", and "reference". Units are
        # meters of elevation and decimal degrees of latitude and longitude

        # Construct file path
        filedir = joinpath(resourcepath,"srtm15plus")
        filepath = joinpath(filedir,"srtm15plus.h5")

        # Download HDF5 file from Google Cloud if necessary
        if ~isfile(filepath)
            print("Downloading srtm15plus.h5 from google cloud storage to $filedir\n")
            run(`mkdir -p $filedir`)
            download("https://storage.googleapis.com/statgeochem/srtm15plus.references.txt", joinpath(filedir,"srtm15plus.references.txt"))
            download("https://storage.googleapis.com/statgeochem/srtm15plus.h5", filepath)
        end

        # Read and return the file
        return h5read(filepath,"vars/"*varname)
    end
    export get_srtm15plus

    # Find the elevation of points at position (lat,lon) on the surface of the
    # Earth, using the SRTM15plus 15-arc-second elevation model.
    function find_srtm15plus(srtm15plus,lat,lon)

        # Interpret user input
        if length(lat) != length(lon)
            error("lat and lon must be equal length\n")
        elseif isa(srtm15plus,Dict)
            data = srtm15plus["elevation"]
        elseif isa(srtm15plus, Array)
            data = srtm15plus
        else
            error("wrong srtm15plus variable")
        end

        # Scale factor (cells per degree) = 60 * 4 = 240
        # (15 arc seconds goes into 1 arc degree 240 times)
        sf = 240

        # Create and fill output vector
        out = Array{Float64}(undef,size(lat))
        for i=1:length(lat)
            if isnan(lat[i]) || isnan(lon[i]) || lat[i]>90 || lat[i]<-90 || lon[i]>180 || lon[i]<-180
                # Result is NaN if either input is NaN or out of bounds
                out[i] = NaN
            else
                # Convert latitude and longitude into indicies of the elevation map array
                # Note that STRTM15 plus has N+1 columns where N = 360*sf
                row = 1 + round(Int,(90+lat[i])*sf)
                col = 1 + round(Int,(180+lon[i])*sf)
                # Find result by indexing
                res = data[row,col]
                if res == -32768
                    out[i] = NaN
                else
                    out[i] = res
                end
            end
        end
        return out
    end
    export find_srtm15plus

## --- End of File