text:
module TestLazyIterators

using VecMathOptInterface
using Test

const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function _test_EmptyVector(T)
    v = MOI.Utilities.EmptyVector{T}()
    @test size(v) == (0,)
    @test length(v) == 0
    @test isempty(v)
    @test eltype(v) == T
    @test iterate(v) === nothing
    c = collect(v)
    @test c isa Vector{T}
    @test isempty(c)
    return
end

test_EmptyVector_Int() = _test_EmptyVector(Int)
test_EmptyVector_Float64() = _test_EmptyVector(Float64)

function _test_LazyMap(T)
    v = MOI.Utilities.LazyMap{T}(x -> x^2, [2, 3])
    @test size(v) == (2,)
    @test length(v) == 2
    @test !isempty(v)
    @test eltype(v) == T
    c = collect(v)
    @test c isa Vector{T}
    @test c == [4, 9]
    return
end

test_LazyMap_Int() = _test_LazyMap(Int)
test_LazyMap_Float64() = _test_LazyMap(Float64)

end  # module

TestLazyIterators.runtests()


text:
include("../src/EDKit.jl")
include("deprecate/ExactDiagonalization.jl")
using .EDKit
using LinearAlgebra
using Test

#-------------------------------------------------------------------------------------------------------------------------
# Tesnsor Basis
#-------------------------------------------------------------------------------------------------------------------------
println("--------------------------------------------------")
println("                   Tensor Basis                   ")
println("--------------------------------------------------")

const X = [0 1; 1 0]
const Y = [0 -1; 1 0] * 1im
const Z = [1 0; 0 -1]

directkron(mat, i, j, L, base=2) = kron(I(base^(i-1)), mat, I(base^(L-j)))
boundarykron(m1, m2, i, j, base=2) = kron(m1, I(base^(j-i+1)), m2)

@testset "Base-2 Random Kron" begin
    L = 10
    for i=1:10
        mats = [rand(2^i, 2^i) for j=1:L-i+1]
        inds = [j:j+i-1 for j=1:L-i+1]
        M_test = operator(mats, inds, L) |> Array
        M_targ = sum(directkron(mats[j], j, j+i-1, L) for j = 1:L-i+1)
        @test M_test ≈ M_targ
    end
end

@testset "Base-3 Random Kron" begin
    L = 6
    for i=1:6
        mats = [rand(3^i, 3^i) for j=1:L-i+1]
        inds = [j:j+i-1 for j=1:L-i+1]
        M_test = operator(mats, inds, L) |> Array
        M_targ = sum(directkron(mats[j], j, j+i-1, L, 3) for j = 1:L-i+1)
        @test M_test ≈ M_targ
    end
end

@testset "Base-3 Boundary Kron" begin
    L = 6
    for i = 1:3, j = 1:3
        m1, m2 = rand(3^i, 3^i), rand(3^j, 3^j)
        mats = kron(m1, m2)
        inds = vcat(Vector(1:i), Vector(L-j+1:L))
        M_test = operator(mats, inds, L) |> Array
        M_targ = boundarykron(m1, m2, i+1, L-j, 3)
        @test M_test ≈ M_targ
    end
end

@testset "Operator Multiplication" begin
    L = 10
    mat = kron(I(2), X, Y, Z)
    v = rand(2^L)
    opt = trans_inv_operator(mat, 4, L)
    v2 = opt * v
    optm = Array(opt)
    v3 = optm * v
    @test v2 ≈ v3
end

#-------------------------------------------------------------------------------------------------------------------------
# Projected Basis
#-------------------------------------------------------------------------------------------------------------------------
println("--------------------------------------------------")
println("                 Projected Basis                  ")
println("--------------------------------------------------")

@testset "Spin-1/2 XY" begin
    L = 10
    mat = spin((1, "xx"), (1, "yy"))
    E = zeros(2^L)
    P = 0
    for n = 0:L
        basis = ProjectedBasis(x->sum(x)==n, L)
        if (l = size(basis, 1)) > 0
            vals = trans_inv_operator(mat, 2, basis) |> Array |> Hermitian |> eigvals
            E[P+1:P+l] = vals
            P += l
        end
    end
    @test P == 2^L
    vals = trans_inv_operator(mat, 2, L) |> Array |> Hermitian |> eigvals
    @test vals ≈ sort!(E)
end

@testset "Spin-1/2 Random" begin
    L = 10
    mat = begin
        M = zeros(ComplexF64, 8, 8)
        M[1,1] = rand()
        M[8,8] = rand()
        M[[2,3,5], [2,3,5]] = rand(ComplexF64, 3, 3) |> Hermitian
        M[[4,6,7], [4,6,7]] = rand(ComplexF64, 3, 3) |> Hermitian
        M
    end
    E = zeros(2^L)
    P = 0
    for n = 0:L
        basis = ProjectedBasis(x->sum(x)==n, L)
        if (l = size(basis, 1)) > 0
            vals = trans_inv_operator(mat, 3, basis) |> Array |> Hermitian |> eigvals
            E[P+1:P+l] = vals
            P += l
        end
    end
    @test P == 2^L
    vals = trans_inv_operator(mat, 3, L) |> Array |> Hermitian |> eigvals
    @test vals ≈ sort!(E)
end

@testset "PXP_OBC" begin
    function fib(n::Integer)
        if n == 0
            return 1
        elseif n == 1
            return 2
        end
        a, b = 1, 2
        for m = 2:n
            a, b = b, a+b
        end
        b
    end
    mat = begin
        P = Diagonal([1, 1, 1, 0, 1, 1, 0, 0])
        P * kron(I(2), X, I(2)) * P
    end
    pxpf(v::Vector{<:Integer}) = all(v[i]==0 || v[i+1]==0 for i=1:length(v)-1)
    for L = 3:20
        basis = ProjectedBasis(pxpf, L)
        @test size(basis, 1) == fib(L)
    end
    for L = 3:10
        basis = ProjectedBasis(pxpf, L)
        mats = fill(mat, L-2)
        inds = [[i,i+1, i+2] for i=1:L-2]
        E = operator(mats, inds, basis) |> Array |> Hermitian |> eigvals
        H = operator(mats, inds, L) |> Array
        vals = H[basis.I, basis.I] |> Hermitian |> eigvals
        @test norm(vals - E) ≈ 0.0
    end
end

@testset "PXP_PBC Projected Basis" begin
    mat = begin
        P = Diagonal([1, 1, 1, 0, 1, 1, 0, 0])
        P * kron(I(2), X, I(2)) * P
    end
    pxpf(v::Vector{<:Integer}) = all(v[i]==0 || v[mod(i, length(v))+1]==0 for i=1:length(v))
    for L = 3:14
        basis = ProjectedBasis(pxpf, L)
        E = trans_inv_operator(mat, 3, basis) |> Array |> Hermitian |> eigvals
        H = trans_inv_operator(mat, 3, L) |> Array
        vals = H[basis.I, basis.I] |> Hermitian |> eigvals
        @test norm(vals - E) ≈ 0.0
    end
end

@testset "Spin-1 Random" begin
    L = 6
    R3 = begin
        m2 = [2, 4, 10]
        m1 = [3, 5, 7, 11, 13, 19]
        m0 = [6, 8, 12, 14, 16, 20, 22]
        p1 = [9, 15, 17, 21, 23, 25]
        p2 = [18, 24, 26]
        mat = zeros(ComplexF64, 27, 27)
        mat[1, 1]    = rand()
        mat[27, 27]  = rand()
        mat[m2, m2] .= rand(ComplexF64, 3, 3) |> Hermitian
        mat[m1, m1] .= rand(ComplexF64, 6, 6) |> Hermitian
        mat[m0, m0] .= rand(ComplexF64, 7, 7) |> Hermitian
        mat[p1, p1] .= rand(ComplexF64, 6, 6) |> Hermitian
        mat[p2, p2] .= rand(ComplexF64, 3, 3) |> Hermitian
        mat
    end
    E = zeros(3^L)
    P = 0
    for n = 0:2L
        basis = ProjectedBasis(x->sum(x)==n, L, base=3)
        if (l = size(basis, 1)) > 0
            vals = trans_inv_operator(mat, 3, basis) |> Array |> Hermitian |> eigvals
            E[P+1:P+l] = vals
            P += l
        end
    end
    @test P == 3^L
    vals = trans_inv_operator(mat, 3, L) |> Array |> Hermitian |> eigvals
    @test vals ≈ sort!(E)
end

#-------------------------------------------------------------------------------------------------------------------------
# Translation Basis
#-------------------------------------------------------------------------------------------------------------------------
println("--------------------------------------------------")
println("               Translational Basis                ")
println("--------------------------------------------------")

@testset "Spin-1/2 XY" begin
    L = 10
    θ = 0.34
    expθ = exp(-1im*θ)
    mat = spin((expθ, "+-"), (1/expθ, "-+"), (1, "z1"), (1, "1z"))
    E = zeros(2^L)
    P = 0
    for n = 0:L, k = 0:L-1
        basis = TranslationalBasis(x->sum(x)==n, k, L)
        #println(size(basis,1))
        if (l = size(basis,1)) > 0
            vals = trans_inv_operator(mat, 2, basis) |> Array |> Hermitian |> eigvals
            E[P+1:P+l] = vals
            P += l
        end
    end
    @test P == 2^L
    vals = trans_inv_operator(mat, 2, L) |> Array |> Hermitian |> eigvals
    @test vals ≈ sort!(E)
end

@testset "Spin-1/2 Random" begin
    L = 10
    mat = begin
        M = zeros(ComplexF64, 8, 8)
        M[1,1] = rand()
        M[8,8] = rand()
        M[[2,3,5], [2,3,5]] = rand(ComplexF64, 3, 3) |> Hermitian
        M[[4,6,7], [4,6,7]] = rand(ComplexF64, 3, 3) |> Hermitian
        M
    end
    E = zeros(2^L)
    P = 0
    for n = 0:L, k = 0:L-1
        basis = TranslationalBasis(x->sum(x)==n, k, L)
        if (l = size(basis, 1)) > 0
            vals = trans_inv_operator(mat, 3, basis) |> Array |> Hermitian |> eigvals
            E[P+1:P+l] = vals
            P += l
        end
    end
    @test P == 2^L
    vals = trans_inv_operator(mat, 3, L) |> Array |> Hermitian |> eigvals
    @test norm(vals - sort(E)) ≈ 0.0 atol = 1e-10
end

@testset "Spin-1 XY" begin
    L = 6
    h = 1
    mat = spin((1, "xx"), (1, "yy"), (h/2, "1z"), (h/2, "z1"), D=3)
    E = zeros(3^L)
    P = 0
    for n = 0:2L, k = 0:L-1
        basis = TranslationalBasis(x->sum(x)==n, k, L, base=3)
        if (l = size(basis, 1)) > 0
            vals = trans_inv_operator(mat, 2, basis) |> Array |> Hermitian |> eigvals
            E[P+1:P+l] = vals
            P += l
        end
    end
    @test P == 3^L
    vals = trans_inv_operator(mat, 2, L) |> Array |> Hermitian |> eigvals
    @test vals ≈ sort!(E)
end

@testset "AKLT" begin
    L = 6
    mat = begin
        ss = spin((1, "xx"), (1, "yy"), (1, "zz"), D=3)
        1/2 * ss + 1/6 * ss^2 + 1/3 * I
    end
    E = zeros(3^L)
    P = 0
    for n = 0:2L, k = 0:L-1
        basis = TranslationalBasis(x->sum(x)==n, k, L, base=3)
        if (l = size(basis, 1)) > 0
            vals = trans_inv_operator(mat, 2, basis) |> Array |> Hermitian |> eigvals
            E[P+1:P+l] = vals
            P += l
        end
    end
    vals = trans_inv_operator(mat, 2, L) |> Array |> Hermitian |> eigvals
    @test vals ≈ sort!(E)
end

@testset "Spin-1 Random" begin
    L = 6
    mat = begin
        m2 = [2, 4, 10]
        m1 = [3, 5, 7, 11, 13, 19]
        m0 = [6, 8, 12, 14, 16, 20, 22]
        p1 = [9, 15, 17, 21, 23, 25]
        p2 = [18, 24, 26]
        M = zeros(ComplexF64, 27, 27)
        M[1, 1]    = rand()
        M[27, 27]  = rand()
        M[m2, m2] .= rand(ComplexF64, 3, 3) |> Hermitian
        M[m1, m1] .= rand(ComplexF64, 6, 6) |> Hermitian
        M[m0, m0] .= rand(ComplexF64, 7, 7) |> Hermitian
        M[p1, p1] .= rand(ComplexF64, 6, 6) |> Hermitian
        M[p2, p2] .= rand(ComplexF64, 3, 3) |> Hermitian
        M
    end
    E = zeros(3^L)
    P = 0
    for n = 0:2L, k = 0:L-1
        basis = TranslationalBasis(x->sum(x)==n, k, L, base=3)
        if (l = size(basis, 1)) > 0
            vals = trans_inv_operator(mat, 3, basis) |> Array |> Hermitian |> eigvals
            E[P+1:P+l] = vals
            P += l
        end
    end
    vals = trans_inv_operator(mat, 3, L) |> Array |> Hermitian |> eigvals
    @test vals ≈ sort!(E)
end

#-------------------------------------------------------------------------------------------------------------------------
# Translation + Parity
#-------------------------------------------------------------------------------------------------------------------------
println("--------------------------------------------------")
println("             Translation Parity Basis             ")
println("--------------------------------------------------")

@testset "Spin-1/2 XY" begin
    mat = spin((1, "+-"), (1, "-+"), (1, "z1"), (1, "1z"))
    for L = 2:2:10
        for n = 0:L, k in [0, L ÷ 2]
            be = TranslationParityBasis(x->sum(x)==n, k, +1, L)
            bo = TranslationParityBasis(x->sum(x)==n, k, -1, L)
            ba = TranslationalBasis(x->sum(x)==n, k, L)
            ve = trans_inv_operator(mat, 2, be) |> Array |> Hermitian |> eigvals
            vo = trans_inv_operator(mat, 2, bo) |> Array |> Hermitian |> eigvals
            va = trans_inv_operator(mat, 2, ba) |> Array |> Hermitian |> eigvals
            E = sort(vcat(ve, vo))
            @test norm(E-va) ≈ 0.0 atol = 1e-12
        end
    end
end

@testset "Spin-1 XY" begin
    mat = spin((1, "+-"), (1, "-+"), (1, "z1"), (1, "1z"), D=3)
    for L = 2:2:6
        for n = 0:2L, k in [0, L ÷ 2]
            be = TranslationParityBasis(x->sum(x)==n, k, +1, L, base=3)
            bo = TranslationParityBasis(x->sum(x)==n, k, -1, L, base=3)
            ba = TranslationalBasis(x->sum(x)==n, k, L, base=3)
            ve = trans_inv_operator(mat, 2, be) |> Array |> Hermitian |> eigvals
            vo = trans_inv_operator(mat, 2, bo) |> Array |> Hermitian |> eigvals
            va = trans_inv_operator(mat, 2, ba) |> Array |> Hermitian |> eigvals
            E = sort(vcat(ve, vo))
            @test norm(E-va) ≈ 0.0 atol = 1e-12
        end
    end
end

@testset "PXP" begin
    mat = begin
        P = Diagonal([1, 1, 1, 0, 1, 1, 0, 0])
        P * kron(I(2), X, I(2)) * P
    end
    pxpf(v::Vector{<:Integer}) = all(v[i]==0 || v[mod(i, length(v))+1]==0 for i=1:length(v))
    for L = 4:2:20
        for k in [0, L ÷ 2]
            be = TranslationParityBasis(pxpf, k, +1, L)
            bo = TranslationParityBasis(pxpf, k, -1, L)
            ba = TranslationalBasis(pxpf, k, L)
            ve = trans_inv_operator(mat, 2, be) |> Array |> Hermitian |> eigvals
            vo = trans_inv_operator(mat, 2, bo) |> Array |> Hermitian |> eigvals
            va = trans_inv_operator(mat, 2, ba) |> Array |> Hermitian |> eigvals
            E = sort(vcat(ve, vo))
            @test norm(E-va) ≈ 0.0 atol = 1e-12
        end
    end
end

@testset "Ising Spin-flip" begin
    zz = spin((1, "zz"))
    x = spin((0.3,"x"))
    for L = 4:2:12
        E = zeros(2^L)
        for k in 0:L-1
            ba = TranslationalBasis(k, L)
            be = TranslationFlipBasis(k, 1, L)
            bo = TranslationFlipBasis(k, -1, L)
            va = trans_inv_operator(zz, 2, ba) + trans_inv_operator(x, 1, ba) |> Array |> Hermitian |> eigvals
            ve = trans_inv_operator(zz, 2, be) + trans_inv_operator(x, 1, be) |> Array |> Hermitian |> eigvals
            vo = trans_inv_operator(zz, 2, bo) + trans_inv_operator(x, 1, bo) |> Array |> Hermitian |> eigvals
            veo = sort(vcat(ve, vo))
            @test norm(va-veo) ≈ 0.0 atol = 1e-12
        end
    end
end

ext: jl


text:
import Base: read, write

"""
    read(fil::String,x::MeshArray)

Read binary file to MeshArray. Other methods:

```
read(xx::Array,x::MeshArray) #from Array
```
"""
function read(fil::String,x::MeshArray)

  grTopo=x.grid.class
  facesSize=x.grid.fSize
  (n1,n2)=x.grid.ioSize
  (nFaces,n3)=nFacesEtc(x)

  fid = open(fil)
  xx = Array{eltype(x),2}(undef,(n1*n2,n3))
  read!(fid,xx)
  xx = hton.(xx)
  close(fid)

  return read(xx,x)

end

function read(xx::Array,x::MeshArray)

  grTopo=x.grid.class
  facesSize=x.grid.fSize
  (n1,n2)=x.grid.ioSize
  (nFaces,n3)=nFacesEtc(x)

  size(xx)!=(n1*n2,n3) ? xx=reshape(xx,(n1*n2,n3)) : nothing

  y=similar(x)
  i0=0; i1=0;
  for iFace=1:nFaces
    i0=i1+1;
    nn=facesSize[iFace][1]; mm=facesSize[iFace][2];
    i1=i1+nn*mm;

    if n3>1 && ndims(x.f)==1
      y.f[iFace]=reshape(xx[i0:i1,:],(nn,mm,n3))
    elseif n3>1 && ndims(x.f)==2
      for i3=1:n3
        y.f[iFace,i3]=reshape(xx[i0:i1,i3],(nn,mm))
      end
    else
      y.f[iFace]=reshape(xx[i0:i1,:],(nn,mm))
    end
  end

  return y

end

"""
    write(fil::String,x::MeshArray)

Write MeshArray to binary file. Other methods:

```
write(xx::Array,x::MeshArray) #to Array
```
"""
function write(fil::String,x::MeshArray)
  y=write(x)
  fid = open(fil,"w")
  write(fid,ntoh.(y))
  close(fid)
end

function write(x::MeshArray)

  grTopo=x.grid.class
  facesSize=x.grid.fSize
  (n1,n2)=x.grid.ioSize
  (nFaces,n3)=nFacesEtc(x)

  y = Array{eltype(x),2}(undef,(n1*n2,n3))
  i0=0; i1=0;
  for iFace=1:nFaces;
    i0=i1+1;
    nn=facesSize[iFace][1];
    mm=facesSize[iFace][2];
    i1=i1+nn*mm;
    if n3>1 && ndims(x.f)==2
      for i3=1:n3
        y[i0:i1,i3]=reshape(x.f[iFace,i3],(nn*mm,1))
      end
    else
      y[i0:i1,:]=reshape(x.f[iFace],(nn*mm,n3))
    end
  end

  y=reshape(y,(n1,n2,n3));
  n3==1 ? y=dropdims(y,dims=3) : nothing

  return y

end


text:
# using packages

# Pkg is designed around "environments"
# this significantly reduces dependencies

julia> Pkg.status()
ERROR: UndefVarError: Pkg not defined
Stacktrace:
 [1] top-level scope at none:0

# enter ']' on the julia repl to go to 'pkg' prompt

julia> ]

# adding the package 'Example' to the current project
(v1.0) pkg> add Example
 Resolving package versions...
 Installed Example ─ v0.5.1
  Updating `C:\Users\praveen\.julia\environments\v1.0\Project.toml`
  [7876af07] + Example v0.5.1
  Updating `C:\Users\praveen\.julia\environments\v1.0\Manifest.toml`
  [7876af07] + Example v0.5.1
  [2a0f44e3] + Base64
  [8ba89e20] + Distributed
  [b77e0a4c] + InteractiveUtils
  [8f399da3] + Libdl
  [37e2e46d] + LinearAlgebra
  [56ddb016] + Logging
  [d6f4376e] + Markdown
  [9a3f8284] + Random
  [9e88b42a] + Serialization
  [6462fe0b] + Sockets
  [8dfed614] + Test

# manifest: manifest status, includes the dependencies of explicitly added packages

# loading a package in julia
(v1.0) pkg> ^C

julia> using Example
[ Info: Precompiling Example [7876af07-990d-54b4-ab0e-23690620f79a]

julia>

julia> Example.hello("User")
"Hello, User"


text:
"""
Return the underlying flattened problem structure (includes generated code).
"""
function get_flatproblem end

"""
A flattened representation of a continuation problem.
"""
struct FlatProblem{G,O}
    var::Vector{Var}
    var_names::Dict{String,Int64}
    data::Vector{Data}
    data_names::Dict{String,Int64}
    func::Vector{Func}
    func_names::Dict{String,Int64}
    group::Vector{Vector{Func}}
    group_names::Dict{Symbol,Int64}
    problem::Vector{Problem}
    problem_names::Dict{String,Int64}
    mfuncs::MonitorFunctions
    call_group::G
    call_owner::O
end

function FlatProblem(mfuncs::MonitorFunctions)
    return FlatProblem(
        Var[],
        Dict{String,Int64}(),
        Data[],
        Dict{String,Int64}(),
        Func[],
        Dict{String,Int64}(),
        Vector{Func}[],
        Dict{Symbol,Int64}(),
        Problem[],
        Dict{String,Int64}(),
        mfuncs,
        nothing,
        nothing,
    )
end

function Base.show(io::IO, mime::MIME"text/plain", flat::FlatProblem)
    println(io, "$FlatProblem()")
    println(io, "    var → $(collect(keys(flat.var_names)))")
    println(io, "    data → $(collect(keys(flat.data_names)))")
    println(io, "    func → $(collect(keys(flat.func_names)))")
    println(io, "    problem → $(collect(keys(flat.problem_names)))")
    print(io, "    group → $(collect(keys(flat.group_names)))")
    return nothing
end

get_flatproblem(flat::FlatProblem) = flat

"""
Take a hierarchical `Problem` structure and return a flattened version that contains
generated code for each function group.
"""
function flatten(problem::Problem)
    flat = FlatProblem()
    _flatten!(flat, problem, "")
    call_group = Tuple(eval(_gen_call_group(flat, i)) for i in eachindex(flat.group))
    call_owner = eval(_gen_call_owner(flat))
    return FlatProblem(
        flat.var,
        flat.var_names,
        flat.data,
        flat.data_names,
        flat.func,
        flat.func_names,
        flat.group,
        flat.group_names,
        flat.problem,
        flat.problem_names,
        flat.mfuncs,
        call_group,
        call_owner,
    )
end

const NAME_SEP = "."

function _flatten!(flat::FlatProblem, problem::Problem, basename)
    basename_sep =
        (isempty(basename) || basename[end] == NAME_SEP) ? basename : basename * NAME_SEP
    push!(flat.problem, problem)
    flat.problem_names[basename] = lastindex(flat.problem)
    # Iterate over the sub-problems (depth first)
    for subproblem in problem.problem
        if isempty(subproblem.name)
            _flatten!(flat, subproblem, basename)
        else
            _flatten!(flat, subproblem, basename_sep * subproblem.name)
        end
    end
    # Iterate over functions in the problem
    for func in problem.func
        if !(func in flat.func)
            # Func has not been previously added, so add it
            push!(flat.func, func)
            fidx = lastindex(flat.func)
            if !isempty(func.name)
                fname = basename_sep * func.name
                if haskey(flat.func_names, fname)
                    @warn "Duplicate Func name" fname func
                else
                    flat.func_names[fname] = fidx
                end
            end
            # Iterate over the groups the Func belongs to
            for group in func.group
                if !haskey(flat.group_names, group)
                    push!(flat.group, Func[])
                    gidx = lastindex(flat.group)
                    flat.group_names[group] = gidx
                else
                    gidx = flat.group_names[group]
                end
                push!(flat.group[gidx], func)
            end
            # Iterate over the Vars belonging to the Func
            for var in func.var
                if !(var in flat.var)
                    # Var has not been previously added, so add it
                    push!(flat.var, var)
                    vidx = lastindex(flat.var)
                    if !isempty(var.name)
                        vname = var.top_level ? var.name : basename_sep * var.name
                        if haskey(flat.var_names, vname)
                            @warn "Duplicate Var name" vname var
                        else
                            flat.var_names[vname] = vidx
                        end
                    end
                end
            end
            # Iterate over the Data belonging to the Func
            for data in func.data
                if !(data in flat.data)
                    # Data has not been previously added, so add it
                    push!(flat.data, data)
                    didx = lastindex(flat.data)
                    if !isempty(data.name)
                        dname = basename_sep * data.name
                        if haskey(flat.data_names, dname)
                            @warn "Duplicate Data name" dname data
                        else
                            flat.data_names[dname] = didx
                        end
                    end
                end
            end
        else
            throw(ErrorException("Duplicate Func within the problem structure"))
        end
    end
    return flat
end

function _gen_call_group(flat::FlatProblem, group::Integer)
    func = :(function (res, u, data, problem) end)
    for (i, f) in enumerate(flat.group[group])
        func_f = :($(f.func)(res[$i]))
        u = :(())
        for var in f.var
            idx = findfirst(==(var), flat.var)
            push!(u.args, :(u[$idx]))
        end
        if length(u.args) == 0
            push!(func_f.args, :(u[:]))
        elseif length(u.args) == 1
            push!(func_f.args, u.args[1])
        else
            push!(func_f.args, u)
        end
        d = :(())
        for data in f.data
            idx = findfirst(==(data), flat.data)
            push!(d.args, :(data[$idx]))
        end
        if length(d.args) == 1
            push!(func_f.args, d.args[1])
        elseif length(d.args) > 1
            push!(func_f.args, d)
        end
        if f.pass_problem
            push!(func_f.args, :problem)
        end
        push!(func.args[2].args, func_f)
    end
    push!(func.args[2].args, :nothing)
    return func
end
function _gen_call_group(flat::FlatProblem, group::Symbol)
    return _gen_call_group(flat, flat.group_names[group])
end

function _gen_call_owner(flat::FlatProblem)
    func = :(function (signal::Signal, problem) end)
    for problem in flat.problem
        if problem.owner !== nothing
            indices = Int64[]
            for item in problem.owner_pass_indices
                if item isa Var
                    idx = findfirst(==(item), flat.var)
                elseif item isa Data
                    idx = findfirst(==(item), flat.data)
                elseif item isa Func
                    idx = findfirst(==(item), flat.func)
                elseif item isa Problem
                    idx = findfirst(==(item), flat.problem)
                else
                    idx = nothing
                end
                if idx === nothing
                    throw(ErrorException("Requested item does not exist in the problem structure: $item"))
                else
                    push!(indices, idx)
                end
            end
            # Using the call signature owner(signal, problem) means that we don't have to
            # worry about whether we are adding methods to the correct function. For
            # example, we could have two covering algorithms that provide an adapt signal -
            # which function do you overload?
            if isempty(indices)
                push!(func.args[2].args, :($(problem.owner)(signal, problem)))
            else
                push!(
                    func.args[2].args,
                    :($(problem.owner)(signal, problem, $((indices...,)))),
                )
            end
        end
    end
    push!(func.args[2].args, :nothing)
    return func
end


text:
function acyclic(network::AbstractNetwork, linkcosts::Vector{<:Real}, order=orderednodes(network); kind=:short)
    nnodes = length(order)
    r = order[1]

    nodecosts = fill((kind == :short) ? Inf : (kind == :long) ? -Inf : error("Only short and long kinds allowed!"), nnodes)
    parentvec = fill(-1, nnodes)
    nodecosts[r] = 0; parentvec[r] = r

    condn = ((kind == :short) ? (i,j) -> (nodecosts[j] > nodecosts[i] + linkcosts[idx(network, i, j)])
                              : (i,j) -> (nodecosts[j] < nodecosts[i] + linkcosts[idx(network, i, j)]))

    for j in order[2:end]
        for i in inneighbors(network, j)
            if condn(i,j)
                nodecosts[j] = nodecosts[i] + linkcosts[idx(network, i, j)]
                parentvec[j] = i
            end
        end
    end

    return nodecosts, parentvec
end


text:
#=

  Utilities submodule package

=#

module Utils

# package dependencies
using RCall: reval, rcopy

include("utils/utils.jl")
include("utils/tree_utils.jl")
include("utils/density_functions.jl")

export rtree, read_tree, make_ape_tree, maketriads, abs_time_branches,
  brts, tree_height, postorderedges, logdexp, logdunifU, logdunif, llrdexp_x,
  logdbeta, llrdbeta_x, logdnorm, logdnorm_tc, llrdnorm_ωxm, llrdnorm_σ²,
  llrdnorm_μ, llrdnorm_x, llrdnorm_xμ, logdtnorm, llrdtnorm_x, erf_custom,
  logdhcau, logdhcau1

end # module Utils

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
function getproba(proba::Vector{Float64}, npaths::Vector{Int}, Ktot)
    if Ktot == -1
        # We have thrown one walk per possible path so npaths does not represent the probability.
        # Therefore we use proba instead.
        proba
    else
        # npaths is representative since the paths were sampled using the probabilities.
        @assert sum(npaths) == Ktot
        npaths ./ Ktot
    end
end
function meanstdpaths(z::Vector{Float64}, proba::Vector{Float64}, npaths::Vector{Int}, Ktot)
    proba = getproba(proba, npaths, Ktot)
    μ = dot(proba, z)
    σ = sqrt(dot(proba, (z .- μ).^2))
    μ, σ
end

mutable struct SolutionPool{TT<:SOI.AbstractTransition, SolT<:SOI.AbstractSolution} <: SOI.AbstractSolutionPool
    # Parent solution
    parent_solution::SolT
    # Are all the child jobs solved feasible ?
    children_feasible::Bool
    # Are all the child jobs solved bounded ?
    children_bounded::Bool
    # List of children solutions
    children_solutions::Dict{TT, SolT}
end
SolutionPool{TT}(sol) where TT = SolutionPool{TT, typeof(sol)}(sol, true, true, Dict{TT, typeof(sol)}())
SOI.getsolution(pool::SolutionPool) = pool.parent_solution
SOI.allfeasible(pool::SolutionPool) = pool.children_feasible
SOI.allbounded(pool::SolutionPool) = pool.children_bounded
SOI.hassolution(pool::SolutionPool{TT}, tr::TT) where TT = haskey(pool.children_solutions, tr)
SOI.getsolution(pool::SolutionPool{TT}, tr::TT) where TT = pool.children_solutions[tr]

mutable struct SDDPPath{TT<:SOI.AbstractTransition, SolT<:SOI.AbstractSolution}
    # Solution pool for the last node of the path
    pool::SolutionPool{TT, SolT}

    # The three following arguments are vectors since when we merge paths,
    # we need to keep track of the difference of `z` that they have accumulated
    # while diverging for computing the variance `σ^2` in [`meanstdpaths`](@ref) at the end.
    # Sum of the objective solution (without θ) at each node of the path
    z::Vector{Float64}
    # Probability of the path
    proba::Vector{Float64}
    # Number of walks in the path
    K::Vector{Int}

    function SDDPPath{TT}(sol::SolT, z, proba, K) where {TT, SolT}
        new{TT, SolT}(SolutionPool{TT}(sol), z, proba, K)
    end
end

function meanstdpaths(paths::Vector{<:SDDPPath}, Ktot)
    z = Compat.reduce(append!, Vector{Float64}[x.z for x in paths], init=Float64[])
    proba = Compat.reduce(append!, Vector{Float64}[x.proba for x in paths], init=Float64[])
    npaths = Compat.reduce(append!, Vector{Int}[x.K for x in paths], init=Int[])
    meanstdpaths(z, proba, npaths, Ktot)
end

function Base.isapprox(p::SDDPPath, q::SDDPPath)
    Base.isapprox(SOI.getnodevalue(SOI.getsolution(p.pool)), SOI.getnodevalue(SOI.getsolution(q.pool)))
end

function canmerge(p::SDDPPath, q::SDDPPath, ztol)
    _isapprox(SOI.getnodevalue(SOI.getsolution(p.pool)), SOI.getnodevalue(SOI.getsolution(q.pool)), ztol)
end

function merge!(p::SDDPPath, q::SDDPPath)
    @assert p.pool.children_feasible == q.pool.children_feasible
    append!(p.z, q.z)
    append!(p.proba, q.proba)
    append!(p.K, q.K)
end

mutable struct Job{SolT<:SOI.AbstractSolution, NodeT, TT<:SOI.AbstractTransition}
    sol::Union{Nothing, SolT}
    proba::Vector{Float64}
    K::Vector{Int}
    parentnode::NodeT
    parent::SDDPPath{TT, SolT}
    tr::TT

    function Job(proba::Vector{Float64}, K::Vector{Int}, parentnode::NodeT, parent::SDDPPath{TT, SolT}, tr::TT) where {SolT, NodeT, TT}
        new{SolT, NodeT, TT}(nothing, proba, K, parentnode, parent, tr)
    end
end

function addjob!(jobsd::Dict{NodeT, Vector{Job{SolT, NodeT, TT}}}, node::NodeT, job::Job{SolT, NodeT, TT}) where {SolT, NodeT, TT}
    if node in keys(jobsd)
        push!(jobsd[node], job)
    else
        jobsd[node] = [job]
    end
end

"""
    mergesamepaths(pathsd::Vector{Tuple{NodeT, Vector{<:SDDPPath}}}, to::TimerOutput, ztol) where NodeT

Find paths that are at the same node with the same parent solution and merge them.
It could happend that two path diverges (go to different node) but then meet again at the same node, if their parent solution is the same then each path will do exactly the same computation so merging them should save time.
"""
function mergesamepaths(pathsd::Vector{Tuple{NodeT, Vector{SDDPPath{TT, SolT}}}}, to::TimerOutput, ztol) where {SolT, NodeT, TT}
    before = sum([sum([sum(path.K) for path in paths]) for (node, paths) in pathsd])
    newpathsd = Tuple{NodeT, Vector{SDDPPath{TT, SolT}}}[]
    @timeit to "merged" for (node, paths) in pathsd
        keep = ones(Bool, length(paths))
        merged = false
        for i in 1:length(paths)
            for j in 1:(i-1)
                if keep[j] && canmerge(paths[i], paths[j], ztol)
                    #println("Merging path since ||x_i - x_j||_∞ = $(norm(getnodevalue(paths[j].sol) - getnodevalue(paths[i].sol), Inf))")
                    merge!(paths[i], paths[j])
                    keep[j] = false
                    merged = true
                    break
                end
            end
        end
        push!(newpathsd, (node, paths[keep]))
    end
    after = sum([sum([sum(path.K) for path in paths]) for (node, paths) in newpathsd])
    @assert before == after
    newpathsd
end

"""
    childjobs(g::SOI.AbstractStochasticProgram, pathsd::Vector{Tuple{NodeT, Vector{<:SDDPPath}}}, pathsampler::AbstractPathSampler, t, num_stages, endedpaths) where SolT

Given paths in `pathsd`, put the paths that have no child in `endedpaths` and sample child jobs using `pathsample` for other paths.
"""
function childjobs(sp::SOI.AbstractStochasticProgram, pathsd::Vector{Tuple{NodeT, Vector{SDDPPath{TT, SolT}}}}, pathsampler::AbstractPathSampler, t, num_stages, endedpaths) where {SolT, NodeT, TT}
    jobsd = Dict{NodeT, Vector{Job{SolT, NodeT, SOI.get(sp, SOI.TransitionType())}}}()
    for (node, paths) in pathsd
        if !isempty(SOI.get(sp, SOI.OutTransitions(), node))
            for path in paths
                # Adding Jobs
                npaths = samplepaths(pathsampler, sp, node, path.K, t, num_stages)
                for (i, tr) in enumerate(SOI.get(sp, SOI.OutTransitions(), node))
                    if !iszero(sum(npaths[i])) || SOI.get(sp, SOI.NeedAllSolutions(), node) # || t == 2
                        addjob!(jobsd, SOI.get(sp, SOI.Target(), tr), Job(path.proba * SOI.get(sp, SOI.Probability(), tr), npaths[i], node, path, tr))
                    end
                end
            end
        else
            append!(endedpaths, paths)
        end
    end
    jobsd
end

"""
    jobstopath(jobsd::Dict{NodeT, Vector{<:Job}}, g::SOI.AbstractStochasticProgram)

Transforms the jobs `jobsd` created by [`childjobs`](@ref) to to paths.
"""
function jobstopaths(jobsd::Dict{NodeT, Vector{Job{SolT, NodeT, TT}}}, g::SOI.AbstractStochasticProgram) where {SolT, NodeT, TT}
    pathsd = Tuple{NodeT, Vector{SDDPPath{TT, SolT}}}[]
    for (node, jobs) in jobsd
        # We create a job even if there is no path going to the node in case
        # we want to create an AveragedCut (since in this case we need to solve all children).
        # However we do not want to create a path for these jobs so we filter them out.
        K = [findall(job.K .!= 0) for job in jobs]
        keep = findall(Bool[jobs[i].parent.pool.children_feasible && !isempty(K[i]) for i in 1:length(jobs)])
        if !isempty(keep)
            paths = SDDPPath{TT, SolT}[SDDPPath{TT}(jobs[i].sol, jobs[i].parent.z[K[i]] .+ jobs[i].sol.objvalx, jobs[i].proba[K[i]], jobs[i].K[K[i]]) for i in keep]
            push!(pathsd, (node, paths))
        end
    end
    pathsd
end

"""
    solvejob!(sp::SOI.AbstractStochasticProgram, job::Job, node, to::TimerOutput)

Solves the job `job` of node `node`.
"""
function solvejob!(sp::SOI.AbstractStochasticProgram, job::Job, node, to::TimerOutput)
    @timeit to "setx" SOI.set!(sp, SOI.SourceSolution(), job.tr, SOI.getsolution(job.parent.pool))
    @timeit to "solve" job.sol = SOI.get(sp, SOI.Solution(), node)
    job.parent.pool.children_solutions[job.tr] = job.sol
    if job.sol.status == :Infeasible
        job.parent.pool.children_feasible = false
    elseif job.sol.status == :Unbounded
        job.parent.pool.children_bounded = false
    end
end


text:
@testset "Opt common" begin
    function graph_test_setting(d = 0)
        coeffs = [4.0 - 0.5im, 0.1im, 2.0 - 0.2im]
        discr = collect(0:9.0) .+ d
        n = length(discr)
        A = hcat(ones(n), discr, discr .^ 2)
        b = coeffs[1] .+ coeffs[2] * discr + coeffs[3] * discr .^ 2
        return (A, b, coeffs, discr)
    end

    # Testing generic complex version
    for linlsqr in [:backslash, :svd, :nrmeq]
        (A, b, coeffs, discr) = graph_test_setting()
        x = solve_linlsqr!(A, b, linlsqr, 0.0)
        @test x ≈ coeffs
        @test !isreal(x)
    end
    (A, b, coeffs, discr) = graph_test_setting()
    #testing droptolerance, high = garbage result
    x = solve_linlsqr!(A, b, :svd, 1e-1)
    @test norm(x - [4.0, 0.0, 2.0]) > 1

    # Testing to force real solution
    for linlsqr in [:real_backslash, :real_svd, :real_nrmeq]
        (A, b, coeffs, discr) = graph_test_setting()
        A = complex(A)
        x = solve_linlsqr!(A, b, linlsqr, 0.0)
        @test x ≈ real(coeffs)
        @test isreal(x)

        (A, b, coeffs, discr) = graph_test_setting(0.1im)
        A = complex(A)
        x = solve_linlsqr!(A, b, linlsqr, 0.0)
        @test isreal(x)
        xx = real(A \ b)
        @test norm(A * x - b) < norm(A * xx - b)
    end
    (A, b, coeffs, discr) = graph_test_setting()
    A = complex(A)
    #testing droptolerance, high = garbage result
    x = solve_linlsqr!(A, b, :real_svd, 1e-1)
    @test norm(x - [4.0, 0.0, 2.0]) > 1

    # Testing BigFloat arithmetic
    (A, b, coeffs, discr) = graph_test_setting()
    A = big.(A)
    b = big.(b)
    x = solve_linlsqr!(A, b, :svd, 0.0)
    @test x ≈ coeffs

    (A, b, coeffs, discr) = graph_test_setting()
    A = big.(A)
    b = big.(b)
    x = solve_linlsqr!(A, b, :real_svd, 0.0)
    @test x ≈ real(coeffs)
    @test isreal(x)
end


text:
@safetestset "isethereumaddress" begin
    using Validator

    @test Validator.isethereumaddress("0x0000000000000000000000000000000000000001") === true
    @test Validator.isethereumaddress("0x683E07492fBDfDA84457C16546ac3f433BFaa128") === true
    @test Validator.isethereumaddress("0x88dA6B6a8D3590e88E0FcadD5CEC56A7C9478319") === true
    @test Validator.isethereumaddress("0x8a718a84ee7B1621E63E680371e0C03C417cCaF6") === true
    @test Validator.isethereumaddress("0xFCb5AFB808b5679b4911230Aa41FfCD0cd335b42") === true

    @test Validator.isethereumaddress("0xGHIJK05pwm37asdf5555QWERZCXV2345AoEuIdHt") === false
    @test Validator.isethereumaddress("0xFCb5AFB808b5679b4911230Aa41FfCD0cd335b422222") === false
    @test Validator.isethereumaddress("0xFCb5AFB808b5679b4911230Aa41FfCD0cd33") === false
    @test Validator.isethereumaddress("0b0110100001100101011011000110110001101111") === false
    @test Validator.isethereumaddress("683E07492fBDfDA84457C16546ac3f433BFaa128") === false
    @test Validator.isethereumaddress("1C6o5CDkLxjsVpnLSuqRs1UBFozXLEwYvU") === false
end


text: const GENIE_PLUGINS = []

max_stars_repo_name: anirudhmurali/hello-julia-genie

max_issues_repo_name: anirudhmurali/hello-julia-genie

max_forks_repo_name: anirudhmurali/hello-julia-genie


text:
#=

  Euler #35 in Julia.

  Problem 35
  """
  The number, 197, is called a circular prime because all rotations
  of the digits: 197, 971, and 719, are themselves prime.

  There are thirteen such primes below 100:
  2, 3, 5, 7, 11, 13, 17, 31, 37, 71, 73, 79, and 97.

  How many circular primes are there below one million?
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# Rotate an array/string
function rotate(n,i)
    s = string(n)
    ii = i+1
    return parse(Int,join(vcat(s[ii:end],s[1:ii-1])))
end

#
# Note: It's a little faster when using prime_set as
#       a parameter for is_circular_prime than
#       relying of the global prime_set in the function.
#       That's surprising...
#
function is_circular_prime(n, prime_set)
    # s = split(string(n),"") # num_to_list(n)
    # len = length(s)
    len = nlen(n)
    v = 0
    for i = 1:len
        v = rotate(n,i)
        if !(v in prime_set)
            return false;
        end
    end
    return v in prime_set
end

# slightly slower than is_circular_prime
function is_circular_prime2(n, prime_set)
    s = split(string(n),"") # num_to_list(n)
    len = length(s)
    v = n
    for i = 2:len
        v = parse(Int,join(vcat([s[j] for j in i:len], [s[j] for j in 1:i-1]),""))
        if !(v in prime_set)
            return false;
        end
    end
    return v in prime_set
end

# Faster if global
prime_set = Set(primes(1_000_000))

# 0.08751774s
function euler35a()
    # prime_set = Set(primes(1_000_000))
    numCircularPrimes = 0
    for n in prime_set
        if is_circular_prime(n,prime_set)
            numCircularPrimes += 1
        end
    end

    return numCircularPrimes
end

# 0.11947861s
function euler35b()
    # prime_set = Set(primes(1_000_000))
    return [1 for n in prime_set if is_circular_prime(n,prime_set)]|>length
end

# Using is_circular_prime2
# 0.24161784s
function euler35c()
    # prime_set = Set(primes(1_000_000))
    # return filter(n->is_circular_prime2(n,prime_set), prime_set)|>length
    return [1 for n in prime_set if is_circular_prime2(n,prime_set)]|>length
end

run_euler(euler35a)
# run_euler(euler35b)
# run_euler(euler35c)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
"""
utils.jl

Provides helper functions to `particle_filtering.jl` such as PGFPlots, rmse compute
Also providese helpers used by `scripts/helpers.jl` such as `truncate_vecs`
"""

# function: get frenet s
"""
    function get_frenet_s(scene;car_id=-1)

# Examples
```julia
true_next_pos = get_frenet_s(true_next_scene,car_id=1)
```
"""
function get_frenet_s(scene;car_id=-1)
    if car_id==-1 print("get_frenet_s says: Give valid car id") end
    veh = scene[findfirst(car_id,scene)]
    return veh.state.posF.s
end

"""
function get_veh_info(scene;car_id = -1)
- Get position and velocity of specific vehicle from scene
"""
function get_veh_info(scene;car_id = -1)
    @assert car_id>0
    pos=scene[findfirst(car_id,scene)].state.posF.s
    vel = scene[findfirst(car_id,scene)].state.v
    return pos,vel
end

# function: get lane id
"""
    function get_lane_id(scene,car_id)
# Examples
```julia
get_lane_id(scene,1)
```
"""
function get_lane_id(scene,car_id)
    veh = scene[findfirst(car_id,scene)]
    return veh.state.posF.roadind.tag.lane
end

# function: get lane change probability
"""
    function get_lane_change_prob(start_scene,particle;car_id=-1,num_samplings=10)
- Probability of lane changing start from `start_scene`
- hallucinating using `particle` for `car_id` using `num_samplings` hallucinations

# Examples
```julia
lp = get_lane_change_prob(scene,particle,car_id = 1,roadway=road_ext)
```
"""
function get_lane_change_prob(f::FilteringEnvironment,start_scene,particle;car_id=-1,
    num_samplings=10)
    if car_id==-1 @show "get_lane_change_prob says: Please give valid car_id" end
    start_lane = get_lane_id(start_scene,car_id)
    changed_count = 0; unchanged_count = 0
    for i in 1:num_samplings
        hpos,hlane = hallucinate_a_step(f,start_scene,particle,car_id=car_id)
        if hlane == start_lane
            unchanged_count += 1
	else
	    changed_count += 1
	end
    end
    return (changed_count+1)/(num_samplings+2)
end

# function: Generate uniform sampling to start the initial particle matrix
"""
    function initial_pmat(;limits,num_particles,seed)
- Generate initial particle matrix with `num_particles` particles with every col being a diff particle
- Range of values that parameters can take is specified in `limits`. Should be num_params rows x 2 cols

# Examples
```julia
limits = [10. 40.;0.1 10.;0.5 5.;1. 10.;0. 1.;-1. 1.;-20. 20.;0. 1.]
initial_pmat(limits=limits,num_particles=10,seed=4)
```
"""
function initial_pmat(;limits,num_particles,seed)
    Random.seed!(seed)
    num_params = size(limits,1)
    p_mat = fill(0.,num_params,num_particles)

    for i in 1:num_params
        p_mat[i,:] = rand(Uniform(limits[i,1],limits[i,2]),1,num_particles)
    end
    return p_mat
end

# function: pgfplots 2 gif using Reel
"""
    function pgfplots2gif(plots;filename="output.gif")
- Make a video using an array of plots. Uses the Reel library.

# Caveat
- Only works with .gif output type. If use .mp4, errors saying process exited ffmpeg
- Relies on the setting Reel.extension
Reel.extension(m::MIME"image/svg+xml") = "svg"
"""
function pgfplots2gif(plots;filename="output.gif")
    @assert typeof(filename) == String

    frames = Frames(MIME("image/svg+xml"), fps=1)

    for (i,plt) in enumerate(plots)
        push!(frames, plt)
        #PGFPlots.save("test/media/$i.pdf",plt) # Store the pic for debugging
    end
    write(filename, frames)
    print("pgfplots2gif: Making animation called $(filename)")
    return nothing
end # End of the reel gif writing function

# function: make pairwise particle variation videos, what pairs to make is the question
"""
- We have a list with corresponding particle matrix at every iteration of filtering
- We shall scatter two columns of said matrix
- We shall make a video with said scatter at every timestep

# Examples
```julia
seed = 2; num_p = 500
Random.seed!(seed)
final_p_mat,iterwise_p_mat = multistep_update(car_id=1,start_frame=2,last_frame=99,num_p=num_p);
plot_pairwise_particles(iterwise_p_mat,filename="test/media/seed2.gif")
```
"""
function plot_pairwise_particles(iterwise_p_mat;filename)
    num_params = size(iterwise_p_mat[1],1) # Number of rows tells us the number of parameters
    print("num_params = $num_params\n")
    pairs = collect(combinations(collect(1:num_params),2))
    num_pairs = length(pairs)
    print("num_pairs = $num_pairs\n")
    groupplots = PGFPlots.GroupPlot[]

        # Parameter names. These better match the top of particle_filtering.jl
    param_names = Dict(1=>"Desired velocity",2=>"Acceleration output noise",
        3=>"Min time headway",4=>"Min separation",5=>"Politeness",6=>"Advantage threshold",
        7=>"Headway sensor noise",8=>"Cooperation");

    for i in 1:length(iterwise_p_mat)
        print("plot_pairwise_particles: iternum = $i\n")
        g = PGFPlots.GroupPlot(7,4)
        p_mat = iterwise_p_mat[i]

        for j in 1:num_pairs
            kk = pairs[j]
            param_1 = kk[1];param_2 = kk[2]
            name_1 = param_names[param_1];name_2=param_names[param_2]
            p = PGFPlots.Axis([PGFPlots.Plots.Scatter(p_mat[param_1,:],p_mat[param_2,:])],
                    xlabel=name_1,ylabel=name_2)
            push!(g,p)
        end
        push!(groupplots,g)
    end

    pgfplots2gif(groupplots,filename=filename)
    return nothing
end

"""
Plot specific pair of parameters eg: v_des vs cooperation

# Examples
```julia
plotspecific(iterwise_p_mat;filename="test/media/vdes_coop.gif")
```
"""
function plotspecific(iterwise_p_mat;filename)
    plots=PGFPlots.Axis[]
    param_names = Dict(1=>"Desired velocity",2=>"Acceleration output noise",
        3=>"Min time headway",4=>"Min separation",5=>"Politeness",6=>"Advantage threshold",
        7=>"Headway sensor noise",8=>"Cooperation");

    name_1 = param_names[1] # 1 is for desired velocity
    name_2 = param_names[8] # 8 is for cooperation

    for i in 1:length(iterwise_p_mat)
        p_mat = iterwise_p_mat[i]

        p = PGFPlots.Axis([PGFPlots.Plots.Scatter(p_mat[1,:],p_mat[8,:])],
        xlabel=name_1,ylabel=name_2)
        push!(plots,p)
    end

    pgfplots2gif(plots,filename=filename)
    return nothing
end

"""
    function avg_dist_particles(p_mat,p_fin)
- Goal is to show that particle filering converges
- `p_mat` is a matrix with particles in every column, say at a certain iteration before convergence
- `p_fin` is the mean of the final particle set i.e. after convergence
- For every particle in `p_mat`, find its norm to `p_fin`. And return the mean of these distances

# Arguments
- `p_mat`: Matrix with particles in every column
- `p_fin`: num_paramx1 vector with the final particle

# Example
```julia
num_iter = length(iterwise_p_mat)
mean_dist_array = fill(0.,num_iter,1)
for i in 1:num_iter
    p_mat = iterwise_p_mat[i]
    mean_dist_array[i] = avg_dist_particles(p_mat,p_fin)
end
```
"""
function avg_dist_particles(p_mat,p_fin)
    val = 0
    num_particles = size(p_mat,2)
    for i in 1:num_particles # loop over the particles
        val+= norm(p_mat[:,i]-p_fin) # add norm distance to final particle
    end
    return val/num_particles # find mean of the norm distance to final particles
    #return sum(sqrt.(sum((p_mat .- p_fin).^2,dims=1)))*1/size(p_mat,2)
end

# function: generate imitation trajectory
"""
function gen_imitation_traj(f::FilteringEnvironment,models;
    id_list=[],start_frame=1,duration=10.)

- Use driver models obtained from filtering to generate imitation trajectory

# Example
```julia
veh_id_list = [6]
f = FilteringEnvironment()
new_models,final_particles,mean_dist = obtain_driver_models(f,veh_id_list,500,1,5)
imit_scene_list = gen_imitation_traj(f,new_models,id_list=veh_id_list)
```
"""
function gen_imitation_traj(f::FilteringEnvironment,models;
    id_list=[],start_frame=1,duration=10.)

    scene_real = f.traj[start_frame]
    if !isempty(id_list) keep_vehicle_subset!(scene_real,id_list) end

    nticks = Int(ceil(duration/f.timestep))
    scene_list = simulate(scene_real,f.roadway,models,nticks,f.timestep)

    return scene_list # Note: simulate appends scenes to the start scene
end

"""
- Idea is to have leader vehicles that are replayed from the data
- models has driver model associated to vehicles in id_list
- Not using `simulate` because we want to customize it. Thus use observe, propogate
- And then inserting the replay vehicles within the scene

# Example
```julia
f=FilteringEnvironment()
new_models,final_particles,mean_dist = obtain_driver_models(f,egoids,30,1,30)
imit_scene_list = imitation_with_replay(f,new_models,egoids=[28,29],replay_ids=[6,8,13])
```
"""
function imitation_with_replay(f::FilteringEnvironment,models;
    egoids=[],replay_ids=[],start_frame=1,duration=10.)

    print("egoids = $egoids\n")
    nticks = Int(ceil(duration/f.timestep))

    start_scene = deepcopy(f.traj[start_frame])

    ego_replay_ids = deepcopy(egoids)
    append!(ego_replay_ids,replay_ids)
    print("ego_replay_ids = $(ego_replay_ids)\n")
    print("start_scene = $start_scene\n")
    keep_vehicle_subset!(start_scene,ego_replay_ids)

    print("start_scene = $start_scene\n")

        # Populate list of ego vehicles from egoids
    ego_vehs = Vector{Entity}(undef,length(egoids))
    for (i, egoid) in enumerate(egoids)
        print("egoid = $egoid\n")
        vehidx = findfirst(egoid, start_scene)
        print("id = $vehidx\n")
        ego_vehs[i] = start_scene[vehidx]
    end

    scenes = [Scene(Entity{VehicleState,VehicleDef,Int64},
                    length(start_scene)) for i=1:nticks+1]

    copyto!(scenes[1], start_scene)

    for tick in 1:nticks
        print("tick = $tick")
        #empty!(scenes[tick + 1])
        traj_scene = deepcopy(f.traj[start_frame+tick])
        keep_vehicle_subset!(traj_scene,replay_ids)

        copyto!(scenes[tick+1],traj_scene)

        for (i, ego_veh) in enumerate(ego_vehs)
            print("scenes[tick] = $(scenes[tick])\n")
            observe!(models[ego_veh.id], scenes[tick], f.roadway, ego_veh.id)
            a = rand(models[ego_veh.id])

            veh_state_p  = propagate(ego_veh, a, f.roadway, f.timestep)

            push!(scenes[tick + 1], Entity(veh_state_p, ego_veh.def, ego_veh.id))
            ego_vehs[i] = Entity(ego_veh,veh_state_p)

        end
    end
    return scenes
end

# function: compute rmse of generated trajectory vs true trajectory
"""
    function compute_rmse(true_scene_list,halluc_scene_list;id_list)
- Compute rmse position and velocity between `halluc_scene_list` and `true_scene_list`
- `true_scene_list` is the ground truth demonstration trajectory
- `imit_scene_list` is generated by running vehicles using parameters that we want for the driver model

# Returns
- `rmse_pos::Dict{Int64,Vector{Float64}}`: Key is vehid. Value is array. Each elem is timewise rmse pos value for that veh_id
- `rmse_vel::Dict{Int64,Vector{Float64}}`: Key is vehid. Value is array. Each elem is timewise rmse vel value for that veh_id

# Examples
```julia
imit_scene_list = gen_imitation_traj(f,new_models,id_list=veh_id_list)
true_scene_list = f.traj[]
rmse_pos_dict,rmse_vel_dict = compute_rmse(true_scene_list,imit_scene_list,id_list=[1,2,3])
```
"""
function compute_rmse(true_scene_list,imit_scene_list;id_list=[])
    @assert length(true_scene_list) == length(imit_scene_list)
    rmse_pos = Dict{Int64,Vector{Float64}}()
    rmse_vel = Dict{Int64,Vector{Float64}}()

    for veh_id in id_list
        rmse_pos[veh_id] = []
        rmse_vel[veh_id] = []
    end

    for i in 1:length(true_scene_list)
        scene_halluc = imit_scene_list[i]
        demo_scene_target = true_scene_list[i]

        for veh_id in id_list
            demo_veh = demo_scene_target[findfirst(veh_id,demo_scene_target)]
            ego_veh = scene_halluc[findfirst(veh_id,scene_halluc)]

            push!(rmse_pos[veh_id],norm(demo_veh.state.posG[1:2]-ego_veh.state.posG[1:2]))
            push!(rmse_vel[veh_id],norm(demo_veh.state.v - ego_veh.state.v))
        end
    end
    return rmse_pos,rmse_vel
end

# function: Combine carwise rmse into one metric by averaging over cars
"""
    function rmse_dict2mean(rmse_dict)
- Take dict of carwise rmse value and return an array with mean of rmse taken over cars

# Returns
- `carmean_rmse`: Vector with length same as num timesteps.
Each element is the rmse value averaged over all cars at that particular timestep

# Examples
```julia
rmse_pos_dict,rmse_vel_dict = compute_rmse(scene_list_1,scene_list_2,id_list=egoids)
rmse_pos = rmse_dict2mean(rmse_pos_dict);rmse_vel = rmse_dict2mean(rmse_vel_dict)
rmse_pos = reshape(rmse_pos,length(rmse_pos),);rmse_vel = reshape(rmse_vel,length(rmse_vel),)
p_pos = PGFPlots.Plots.Linear(collect(1:length(rmse_pos)),rmse_pos,legendentry="rmse pos")
p_vel = PGFPlots.Plots.Linear(collect(1:length(rmse_vel)),rmse_vel,legendentry="rmse vel")
rmse_axis = PGFPlots.Axis([p_pos,p_vel],xlabel="timestep",ylabel="rmse",title="rmse pos and vel")
display(rmse_axis)
```
"""
function rmse_dict2mean(rmse_dict)
    num_veh = length(collect(keys(rmse_dict))) # Find length of the vector of keys
    num_iter = length(rmse_dict[collect(keys(rmse_dict))[1]]) # Find length of value contained in 1st key
    rmse_array = fill(0.,num_iter,num_veh)

    i = 0
    for (k,v) in rmse_dict
        i = i+1
        rmse_vals = reshape(v,length(v),1)
        rmse_array[:,i] = rmse_vals
    end

    carmean_rmse = mean(rmse_array,dims=2)
    carmean_rmse = reshape(carmean_rmse,length(carmean_rmse),)
    return carmean_rmse
end

# Function: Test collisons for list of vehicle ids
"""
- Adapt the get collision to work with vehicle ids and give us collision metrics within notebook
- Source: `AutomotiveDrivingModels/src/Collision-Checkers/Minkowski.jl`

# Example
```julia
f = FilteringEnvironment()
c = collision_check(f.traj[1],[6,8])
c.is_colliding
```
"""
function collision_check(scene, veh_id_list, mem::CPAMemory=CPAMemory())

    N = length(veh_id_list)
    for (a,A) in enumerate(veh_id_list)
        vehA = scene[findfirst(A,scene)]
        to_oriented_bounding_box!(mem.vehA, vehA)
        for b in a +1 : length(veh_id_list)
            B = veh_id_list[b]
            vehB = scene[findfirst(B,scene)]
            if is_potentially_colliding(vehA, vehB)
                to_oriented_bounding_box!(mem.vehB, vehB)
                if is_colliding(mem)
                    return CollisionCheckResult(true, A, B)
                end
            end
        end
    end

    CollisionCheckResult(false, 0, 0)
end

# function: Find the collisions instances to extract collision metrics
"""
    function test_collision
- Assess list of scenes for collision information

# Returns
- Binary array with 0 if no collision at that timestep, and 1 if any collision at that timestep

# Examples
```julia

```
"""
function test_collision(scenes_list,id_list)

    collisions_array = fill(0.,length(scenes_list))

    for (i,scene) in enumerate(scenes_list)
        if collision_check(scene,id_list).is_colliding
            collisions_array[i] = 1
        end
    end
    return collisions_array
end

"""
function frac_colliding_timesteps(coll_mat)
- `coll_mat` has different scenarios in different columns
- Each column has 1 in every timestep where there is at least one pair of cars colliding
- We take all the timesteps i.e. N scenarios x H timesteps in each scenario
- Find the fraction of timesteps that have a collision
"""
function frac_colliding_timesteps(coll_mat)
        n = size(coll_mat,2) # num columns i.e. scenarios
        h = size(coll_mat,1) # num timesteps per scenario
        return sum(coll_mat)/(n*h)
end

"""
function pgfplot_vector(vec;leg = "nolegend")

- Plot a vector using PGFPlots

# Example
```julia
plot_rmse_pf = pgfplot_vector(pos_pf,leg="pf")
```
"""
function pgfplot_vector(vec;leg = "nolegend")
    vec = reshape(vec,length(vec),)
    return PGFPlots.Plots.Linear(collect(1:length(vec)),vec,legendentry=leg)
end

# function to truncate vectors to shortest length one
"""
function truncate_vecs(list_of_vectors)

- Rmse vec lengths will not be same as simulation durations not the same
- Select the shortest length rmse and truncate other scenarios to same length

# Example
```julia
intset = collect(1:10) #integer set to draw random numbers from
l = [rand(intset,4,1),rand(intset,5,1),rand(intset,6,1),rand(intset,10,1),rand(intset,3,1)]
a = truncate_vecs(l)
```
"""
function truncate_vecs(list_of_vectors)
    shortest_length = 10000
    numvec = length(list_of_vectors)
    for i in 1:numvec
        curr_len = length(list_of_vectors[i])
        if curr_len < shortest_length
            shortest_length = curr_len
        end
    end
    #print("longest length = $(longest_length)\n")

    new_list = []

    for j in 1:numvec
        #print("j = $j\n")
        # pad with zeros
        curr_vec = list_of_vectors[j]
        truncated_vec = curr_vec[1:shortest_length]
        #print("padded_vec = $(padded_vec)\n")
        push!(new_list,truncated_vec)
    end

    return hcat(new_list...)
end

# Elongate vectors to longest in the list by padding zeros
"""
function pad_zeros(list_of_vectors)
- We need to append zeros to the shorter rmse vectors to make all lengths equal

# Example
```julia
intset = collect(1:10) #integer set to draw random numbers from
l = [rand(intset,4,1),rand(intset,5,1),rand(intset,6,1),rand(intset,10,1),rand(intset,3,1)]
a = pad_zeros(l)
```
"""
function pad_zeros(list_of_vectors)
    longest_length = -1
    numvec = length(list_of_vectors)
    for i in 1:numvec
        curr_len = length(list_of_vectors[i])
        if curr_len > longest_length
            longest_length = curr_len
        end
    end
    #print("longest length = $(longest_length)\n")

    new_list = []

    for j in 1:numvec
        #print("j = $j\n")
        # pad with zeros
        curr_vec = list_of_vectors[j]
        lengthdiff = longest_length-length(curr_vec)
        padded_vec = vcat(curr_vec,fill(0.,lengthdiff,1))
        #print("padded_vec = $(padded_vec)\n")
        push!(new_list,padded_vec)
    end

    return hcat(new_list...)
end

"""
- Get a list of scenes from replay trajectory
- This will be used to compute metrics of real driving behavior

# Example
```julia
f = FilteringEnvironment()
models,id_list,ts,te = JLD.load(filename,"m","veh_id_list","ts","te")
scene_list = replay_scenelist(f,id_list=id_list,ts=ts,te=te)
```
"""
function replay_scenelist(f::FilteringEnvironment;id_list=[],ts,te)
    scene_list = deepcopy(f.traj[ts:te])
    for scene in scene_list
        if !isempty(id_list) keep_vehicle_subset!(scene,id_list) end
    end
    return scene_list
end

"""
function reality_metrics(scene_list)
- Compute things about the driving trajectory to assess how `real` the traffic is

# Example
```julia
f = FilteringEnvironment()
filename = "media/upper_4.jld"
models,id_list,ts,te = JLD.load(filename,"m","veh_id_list","ts","te")

scene_list_true = replay_scenelist(f,id_list=id_list,ts=ts,te=te)
stopfrac_true,brakefrac_true = reality_metrics(f,scene_list_true,id_list=id_list)

scene_real = deepcopy(f.traj[ts])
if !isempty(id_list) keep_vehicle_subset!(scene_real,id_list) end
nticks = te-ts+1
scene_list_generated = simulate(start_scene,f.roadway,models,nticks,f.timestep)
stopfrac_gen,brakefrac_gen = reality_metrics(f,scene_list_generated,id_list=id_list)
```
"""
function reality_metrics(f::FilteringEnvironment,scene_list;id_list=[])
    # Calculate number of timesteps where more than one car almost stopped
    # This is based on visual in May 7 update slides turing test

    numscenes = length(scene_list)
    stopped_scenes = fill(0.,numscenes,)
    for (i,scene) in enumerate(scene_list)
        laggards = 0
        for veh in scene
            if veh.state.v < 2
                laggards += 1
            end
        end
        if laggards>1
            stopped_scenes[i] = 1
        end
    end
    stopfrac = sum(stopped_scenes)/numscenes

    # Calculate hard deceleration instances
    numvehs = length(id_list)
    acc_mat = fill(0.,numscenes-1,numvehs) #-1 because 1st elem of acc is `missing`
    #print("size(acc_mat) = $(size(acc_mat))\n")
    for (j,veh_id) in enumerate(id_list)
        print("veh_id = $veh_id\n")
        acc_trace = acc(f.roadway,scene_list,veh_id)

        trace = Float64.(acc_trace[2:end]) # Convert Union missing type to just Float
        #popfirst!(trace) # Remove the first element i.e. missing

        #print("size(trace) = $(size(trace))\n")
        trace = reshape(trace,numscenes-1,)
        #print("size(trace) = $(size(trace))\n")

        acc_mat[:,j] .= trace
    end

    num_harddecel = 0
    for jj in 1:numscenes-1
        if any(x->x<-2,acc_mat[jj,:]) # If any vehicle shows accel less than -2, its hard decel
            num_harddecel += 1
        end
    end
    #print("num_harddecel = $(num_harddecel)\n")
    harddecelfrac = num_harddecel/(numscenes-1)

    return stopfrac,harddecelfrac
end

max_stars_repo_name: sisl/AutomotiveInteraction

max_issues_repo_name: sisl/AutomotiveInteraction

max_forks_repo_name: sisl/AutomotiveInteraction


text:
# ---
# title: 721. Accounts Merge
# id: problem721
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Depth-first Search, Union Find
# link: <https://leetcode.com/problems/accounts-merge/description/>
# hidden: true
# ---
#
# Given a list `accounts`, each element `accounts[i]` is a list of strings,
# where the first element `accounts[i][0]` is a _name_ , and the rest of the
# elements are _emails_ representing emails of the account.
#
# Now, we would like to merge these accounts. Two accounts definitely belong to
# the same person if there is some email that is common to both accounts. Note
# that even if two accounts have the same name, they may belong to different
# people as people could have the same name. A person can have any number of
# accounts initially, but all of their accounts definitely have the same name.
#
# After merging the accounts, return the accounts in the following format: the
# first element of each account is the name, and the rest of the elements are
# emails **in sorted order**. The accounts themselves can be returned in any
# order.
#
# **Example 1:**
#
#
#
#     Input:
#     accounts = [["John", "johnsmith@mail.com", "john00@mail.com"], ["John", "johnnybravo@mail.com"], ["John", "johnsmith@mail.com", "john_newyork@mail.com"], ["Mary", "mary@mail.com"]]
#     Output: [["John", 'john00@mail.com', 'john_newyork@mail.com', 'johnsmith@mail.com'],  ["John", "johnnybravo@mail.com"], ["Mary", "mary@mail.com"]]
#     Explanation:
#     The first and third John's are the same person as they have the common email "johnsmith@mail.com".
#     The second John and Mary are different people as none of their email addresses are used by other accounts.
#     We could return these lists in any order, for example the answer [['Mary', 'mary@mail.com'], ['John', 'johnnybravo@mail.com'],
#     ['John', 'john00@mail.com', 'john_newyork@mail.com', 'johnsmith@mail.com']] would still be accepted.
#
#
# **Note:**
#
# * The length of `accounts` will be in the range `[1, 1000]`.
# * The length of `accounts[i]` will be in the range `[1, 10]`.
# * The length of `accounts[i][j]` will be in the range `[1, 30]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module SkewDist

using Distributions, StatsFuns, MultivariateStats, PDMats, Optim, Roots
import Base: length, mean, show, rand, var, cov, *
import Distributions: sampler, _rand!, pdf, _logpdf, _logpdf!, dof, cdf, quantile, minimum, maximum, mgf

export MvSkewNormal, MvSkewTDist, SkewNormal, SkewTDist, fit_MvSkewNormal, fit_MvSkewTDist, marginals, mgf
export pdf, dof, cdf, quantile

include("utils.jl")
include("mvskewnormal.jl")
include("mvskewtdist.jl")
include("skewnormal.jl")
include("skewtdist.jl")
include("operators.jl")
include("fit_mvskewtdist.jl")

end # module


text:
function test_filters()
    test_fill_missing_availability()
    test_downgrade_availability_noop()
    test_downgrade_shifts()
    test_downgrade_hours()
    test_downgrade_hours_and_shifts()
    test_only_available_when_business_closed_downgrade()
end

function test_fill_missing_availability()
    schedule = {
        "demand" => {
            "sunday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "monday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "tuesday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "wednesday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "thursday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "friday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "saturday" => [1, 2, 3, 3, 2, 1, 1, 1],
        },
    }

    organization = {
        "day_week_starts" => "sunday",
    }

    workers = Any[
        {
            "id" => 1,
            "name" => "Nikola Tesla",
            "max_shifts_per_week" => 2,
        },
        {
            "id" => 2,
            "name" => "Richard Feynman",
            "max_shifts_per_week" => 2,
        },
        {
            "id" => 3,
            "name" => "Dean Quatrano",
            "max_shifts_per_week" => 0,
        },
    ]

    availability = Any[
        {
            "user_id" => 2,
            "availability" => {
                "sunday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "monday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "tuesday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "wednesday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "thursday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "friday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "saturday" => [1, 0, 1, 1, 0, 1, 1, 1],
            },
        }
        {
            "user_id" => 3,
            "availability" => {
                "sunday" => [1, 0, 0, 0, 0, 0, 1, 1],
                "monday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "tuesday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "wednesday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "thursday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "friday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "saturday" => [1, 0, 1, 1, 0, 1, 1, 1],
            },
        }
    ]

    expected = Any[
        {
            "id" => 1,
            "name" => "Nikola Tesla",
            "max_shifts_per_week" => 2,
            "availability" => {
                "sunday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "monday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "tuesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
        {
            "id" => 2,
            "name" => "Richard Feynman",
            "max_shifts_per_week" => 2,
            "availability" => {
                "sunday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "monday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "tuesday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "wednesday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "thursday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "friday" => [1, 0, 1, 1, 0, 1, 1, 1],
                "saturday" => [1, 0, 1, 1, 0, 1, 1, 1],
            },
        },
    ]

    actual = fill_missing_availability(schedule, availability, workers, organization, Dict(), false)

    @test expected == actual
end

# Everybody meets constraints - don't touch them
function test_downgrade_availability_noop()
    schedule = {
        "demand" => {
            "sunday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "monday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "tuesday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "wednesday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "thursday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "friday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "saturday" => [1, 2, 3, 3, 2, 1, 1, 1],
        },
    }

    organization = {
        "day_week_starts" => "sunday",
    }

    workers = Any[
        {
            "id" => 1,
            "name" => "Nikola Tesla",
            "min_shift_length" => 1,
            "max_shift_length" => 9,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 2,
            "min_hours_per_week" => 1,
            "max_hours_per_week" => 20,
            "availability" => {
                "sunday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "monday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "tuesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
        {
            "id" => 2,
            "name" => "Richard Feynman",
            "min_shift_length" => 1,
            "max_shift_length" => 9,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 2,
            "min_hours_per_week" => 1,
            "max_hours_per_week" => 20,
            "availability" => {
                "sunday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "monday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "tuesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
    ]

    actual = downgrade_availability(schedule, workers, organization, Dict(), false)

    @test workers == actual
end

function test_downgrade_shifts()
    schedule = {
        "demand" => {
            "sunday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "monday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "tuesday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "wednesday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "thursday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "friday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "saturday" => [1, 2, 3, 3, 2, 1, 1, 1],
        },
    }

    organization = {
        "day_week_starts" => "sunday",
    }

    workers = Any[
        # Adjust this guy's shift max
        {
            "id" => 1,
            "name" => "Nikola Tesla",
            "min_shift_length" => 2,
            "max_shift_length" => 9,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 7,
            "min_hours_per_week" => 1,
            "max_hours_per_week" => 20,
            "availability" => {
                "sunday" => [0, 1, 0, 0, 0, 0, 0, 0],
                "monday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "tuesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
        {   # Adjust this guy's shift min and max
            "id" => 2,
            "name" => "Richard Feynman",
            "min_shift_length" => 1,
            "max_shift_length" => 9,
            "min_shifts_per_week" => 5,
            "max_shifts_per_week" => 7,
            "min_hours_per_week" => 1,
            "max_hours_per_week" => 20,
            "availability" => {
                "sunday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "monday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "tuesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
    ]

    expected = Any[
        # Adjust this guy's shift max
        {
            "id" => 1,
            "name" => "Nikola Tesla",
            "min_shift_length" => 2,
            "max_shift_length" => 9,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 6,
            "min_hours_per_week" => 2,
            "max_hours_per_week" => 20,
            "availability" => {
                # Annealed
                "sunday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "monday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "tuesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
        {   # Adjust this guy's shift min and max
            "id" => 2,
            "name" => "Richard Feynman",
            "min_shift_length" => 1,
            "max_shift_length" => 9,
            "min_shifts_per_week" => 4,
            "max_shifts_per_week" => 4,
            "min_hours_per_week" => 4,
            "max_hours_per_week" => 20,
            "availability" => {
                "sunday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "monday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "tuesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
    ]

    actual = downgrade_availability(schedule, workers, organization, Dict(), false)

    @test expected == actual
end

function test_downgrade_hours()
    schedule = {
        "demand" => {
            "sunday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "monday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "tuesday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "wednesday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "thursday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "friday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "saturday" => [1, 2, 3, 3, 2, 1, 1, 1],
        },
    }

    organization = {
        "day_week_starts" => "sunday",
    }

    workers = Any[
        # Adjust this guy's shift max
        {
            "id" => 1,
            "name" => "Nikola Tesla",
            "min_shift_length" => 2,
            "max_shift_length" => 3,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 7,
            "min_hours_per_week" => 1,
            "max_hours_per_week" => 80,
            "availability" => {
                "sunday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "monday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "tuesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
        {   # Adjust this guy's shift min and max
            "id" => 2,
            "name" => "Richard Feynman",
            "min_shift_length" => 5,
            "max_shift_length" => 5,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 1,
            "min_hours_per_week" => 6,
            "max_hours_per_week" => 20,
            "availability" => {
                "sunday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "monday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "tuesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
    ]

    expected = Any[
        # Adjust this guy's shift max
        {
            "id" => 1,
            "name" => "Nikola Tesla",
            "min_shift_length" => 2,
            "max_shift_length" => 3,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 7,
            "min_hours_per_week" => 2,
            "max_hours_per_week" => 21,
            "availability" => {
                # Annealed
                "sunday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "monday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "tuesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
        {   # Adjust this guy's shift min and max
            "id" => 2,
            "name" => "Richard Feynman",
            "min_shift_length" => 5,
            "max_shift_length" => 5,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 1,
            "min_hours_per_week" => 5,
            "max_hours_per_week" => 5,
            "availability" => {
                "sunday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "monday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "tuesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
    ]

    actual = downgrade_availability(schedule, workers, organization, Dict(), false)

    @test expected == actual
end

function test_downgrade_hours_and_shifts()
    schedule = {
        "demand" => {
            "sunday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "monday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "tuesday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "wednesday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "thursday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "friday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "saturday" => [1, 2, 3, 3, 2, 1, 1, 1],
        },
    }

    organization = {
        "day_week_starts" => "sunday",
    }

    workers = Any[
        {
            "id" => 1,
            "name" => "Nikola Tesla",
            "min_shift_length" => 2,
            "max_shift_length" => 3,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 7,
            "min_hours_per_week" => 1,
            "max_hours_per_week" => 80,
            "availability" => {
                "sunday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "monday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "tuesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
        {
            "id" => 2,
            "name" => "Richard Feynman",
            "min_shift_length" => 4,
            "max_shift_length" => 5,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 2,
            "min_hours_per_week" => 3,
            "max_hours_per_week" => 20,
            "availability" => {
                "sunday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "monday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "tuesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "wednesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "thursday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "friday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
        {
            "id" => 3,
            "name" => "That weird guy who smells like onions",
            "min_shift_length" => 5,
            "max_shift_length" => 8,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 4,
            "min_hours_per_week" => 16,
            "max_hours_per_week" => 32,
            "availability" => {
                "sunday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "monday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "tuesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
    ]

    expected = Any[
        # Adjust this guy's shift max
        {
            "id" => 1,
            "name" => "Nikola Tesla",
            "min_shift_length" => 2,
            "max_shift_length" => 3,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 4,
            "min_hours_per_week" => 2,
            "max_hours_per_week" => 12,
            "availability" => {
                "sunday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "monday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "tuesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
        {   # Adjust this guy's shift min and max
            "id" => 2,
            "name" => "Richard Feynman",
            "min_shift_length" => 4,
            "max_shift_length" => 5,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 1,
            "min_hours_per_week" => 4,
            "max_hours_per_week" => 5,
            "availability" => {
                "sunday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "monday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "tuesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "wednesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "thursday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "friday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
        {
            "id" => 3,
            "name" => "That weird guy who smells like onions",
            "min_shift_length" => 5,
            "max_shift_length" => 8,
            "min_shifts_per_week" => 2,
            "max_shifts_per_week" => 4,
            "min_hours_per_week" => 16,
            "max_hours_per_week" => 32,
            "availability" => {
                "sunday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "monday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "tuesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
    ]

    actual = downgrade_availability(schedule, workers, organization, Dict(), false)
    @test expected == actual
end

function test_only_available_when_business_closed_downgrade()
    # Basically same as last test, but some unavailability comes from
    # the business being closed and the worker not having enough hours
    # on the day to have it count.

    schedule = {
        "demand" => {
            "sunday" => [0, 0, 0, 0, 2, 1, 1, 1],
            "monday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "tuesday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "wednesday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "thursday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "friday" => [1, 2, 3, 3, 2, 1, 1, 1],
            "saturday" => [1, 2, 3, 3, 2, 1, 1, 1],
        },
    }

    organization = {
        "day_week_starts" => "sunday",
    }

    workers = Any[
        {
            "id" => 1,
            "name" => "Nikola Tesla",
            "min_shift_length" => 2,
            "max_shift_length" => 3,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 7,
            "min_hours_per_week" => 1,
            "max_hours_per_week" => 80,
            "availability" => {
                "sunday" => [1, 1, 1, 1, 0, 0, 0, 0],
                "monday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "tuesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
        {
            "id" => 2,
            "name" => "Richard Feynman",
            "min_shift_length" => 4,
            "max_shift_length" => 5,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 2,
            "min_hours_per_week" => 3,
            "max_hours_per_week" => 20,
            "availability" => {
                "sunday" => [1, 1, 1, 0, 0, 0, 0, 0],
                "monday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "tuesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "wednesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "thursday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "friday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
    ]

    expected = Any[
        # Adjust this guy's shift max
        {
            "id" => 1,
            "name" => "Nikola Tesla",
            "min_shift_length" => 2,
            "max_shift_length" => 3,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 4,
            "min_hours_per_week" => 2,
            "max_hours_per_week" => 12,
            "availability" => {
                "sunday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "monday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "tuesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "wednesday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "thursday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "friday" => [1, 1, 1, 1, 1, 1, 1, 1],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
        {   # Adjust this guy's shift min and max
            "id" => 2,
            "name" => "Richard Feynman",
            "min_shift_length" => 4,
            "max_shift_length" => 5,
            "min_shifts_per_week" => 1,
            "max_shifts_per_week" => 1,
            "min_hours_per_week" => 4,
            "max_hours_per_week" => 5,
            "availability" => {
                "sunday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "monday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "tuesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "wednesday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "thursday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "friday" => [0, 0, 0, 0, 0, 0, 0, 0],
                "saturday" => [1, 1, 1, 1, 1, 1, 1, 1],
            },
        },
    ]

    actual = downgrade_availability(schedule, workers, organization, Dict(), false)
    @test expected == actual
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
"""
    poweriter(A,numiter)

Perform `numiter` power iterations with the matrix `A`, starting
from a random vector. Returns a vector of eigenvalue estimates
and the final eigenvector approximation.
"""
function poweriter(A,numiter)
    n = size(A,1)
    x = normalize(randn(n),Inf)
    β = zeros(numiter)
    for k in 1:numiter
        y = A*x
        m = argmax(abs.(y))
        β[k] = y[m]/x[m]
        x = y/y[m]
    end
    return β,x
end

"""
    inviter(A,s,numiter)

Perform `numiter` inverse iterations with the matrix `A` and shift
`s`, starting from a random vector. Returns a vector of
eigenvalue estimates and the final eigenvector approximation.
"""
function inviter(A,s,numiter)
    n = size(A,1)
    x = normalize(randn(n),Inf)
    β = zeros(numiter)
    fact = lu(A - s*I)
    for k in 1:numiter
        y = fact\x
        normy,m = findmax(abs.(y))
        β[k] = x[m]/y[m] + s
        x = y/y[m]
    end
    return β,x
end

"""
    arnoldi(A,u,m)

Perform the Arnoldi iteration for `A` starting with vector `u`, out
to the Krylov subspace of degree `m`. Returns the orthonormal basis
(`m`+1 columns) and the upper Hessenberg `H` of size `m`+1 by `m`.
"""
function arnoldi(A,u,m)
    n = length(u)
    Q = zeros(n,m+1)
    H = zeros(m+1,m)
    Q[:,1] = u/norm(u)
    for j in 1:m
        # Find the new direction that extends the Krylov subspace.
        v = A*Q[:,j]
        # Remove the projections onto the previous vectors.
        for i in 1:j
            H[i,j] = dot(Q[:,i],v)
            v -= H[i,j]*Q[:,i]
        end
        # Normalize and store the new basis vector.
        H[j+1,j] = norm(v)
        Q[:,j+1] = v/H[j+1,j]
    end
    return Q,H
end

"""
    gmres(A,b,m)

Do `m` iterations of GMRES for the linear system `A`*x=`b`. Returns
the final solution estimate x and a vector with the history of
residual norms. (This function is for demo only, not practical use.)
"""
function gmres(A,b,m)
    n = length(b)
    Q = zeros(n,m+1)
    Q[:,1] = b/norm(b)
    H = zeros(m+1,m)

    # Initial solution is zero.
    x = 0
    residual = [norm(b);zeros(m)]

    for j in 1:m
        # Next step of Arnoldi iteration.
        v = A*Q[:,j]
        for i in 1:j
            H[i,j] = dot(Q[:,i],v)
            v -= H[i,j]*Q[:,i]
        end
        H[j+1,j] = norm(v)
        Q[:,j+1] = v/H[j+1,j]

        # Solve the minimum residual problem.
        r = [norm(b); zeros(j)]
        z = H[1:j+1,1:j] \ r
        x = Q[:,1:j]*z
        residual[j+1] = norm( A*x - b )
    end
    return x,residual
end

### These are not part of the text proper, but help to construct
### interesting examples.

"""
    sprandsym(n,density,λ)
    sprandsym(n,density,rcond)

Construct a randomized `n` by `n` symmetric sparse matrix of approximate
density `density`. For vector `λ`, the matrix has eigenvalues as
prescribed by λ. For scalar `rcond`, the matrix has condition number
equal to 1/`rcond`.
"""
function sprandsym(n,density,rcond::Number)
    λ = [ rcond^(i/(n-1)) for i in 0:n-1 ]
    sprandsym(n,density,λ)
end

function sprandsym(n,density,λ::Vector)
    # Apply a random Jacobi rotation similarity transformation.
    function randjr(A)
        theta = 2π*rand()
        c,s = cos(theta),sin(theta)
        i = j = 0
        while i==j
          i,j = rand(1:n,2)
        end
        A[[i,j],:] = [c s;-s c]*A[[i,j],:]
        A[:,[i,j]] = A[:,[i,j]]*[c -s;s c]
        return A
    end

    targetnz = ceil(min(0.98,density)*n^2)
    A = spdiagm(λ)

    while nnz(A) < targetnz
        A = randjr(A)
    end
    return A
end

"""
    poisson(n)

Construct the finite difference Laplacian matrix for an `n` by `n`
interior lattice.
"""
function poisson(n)
    D = spdiagm(-1=>fill(-1,n-1),0=>fill(2,n),1=>fill(-1,n-1)) * (n+1)^2/pi^2
    return kron(D,I(n)) + kron(I(n),D)
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ dd29cb1a-4770-11eb-0db2-272859d3b045
begin
	using FSM
	using CSV
	using DataFrames
	using Plots
	using PlutoUI
end

# ╔═╡ a2eb46ba-4770-11eb-1d86-4d2e18f25bc1
md"""

### Flexible snow model

Here are some experiments with the Julia implementation of FSM...

"""

# ╔═╡ f126375a-4772-11eb-11f8-a71c6c86ba36
md"""
Snow albedo model
"""

# ╔═╡ 125fc900-4772-11eb-3d6b-d121f32c31c5
@bind am Select(["0", "1"])

# ╔═╡ 2b7b1710-49fd-11eb-072a-e320f2a80a3c
md"""
Snow conductivity model
"""

# ╔═╡ 36b71ff0-49fd-11eb-0cc8-9f7f0fdac201
@bind cm Select(["0", "1"])

# ╔═╡ 96f8c1f2-49fd-11eb-19ff-8db150ad53a3
md"""
Snow density model
"""

# ╔═╡ 9fd8aa12-49fd-11eb-32cb-bdd6f962d075
@bind dm Select(["0", "1"])

# ╔═╡ bfe40edc-49fd-11eb-0268-570ff06b4020
md"""
Surface exchange model
"""

# ╔═╡ c9ef08a0-49fd-11eb-07ed-252c5d4da577
@bind em Select(["0", "1"])

# ╔═╡ 1207461a-4773-11eb-3b2e-a181f28c29ba
md"""
Snow hydraulics model
"""

# ╔═╡ 1e2df4f4-4773-11eb-3f1f-dd30364bc006
@bind hm Select(["0", "1"])

# ╔═╡ 3db3b4e0-5825-11eb-22e3-b1fe7436794a
md"""
## Result plots...
"""

# ╔═╡ 556f9b8e-4773-11eb-2b6b-abe07b0efdc1
md"""Setup simulations..."""

# ╔═╡ 980e8ac4-4771-11eb-0a52-c35a9a0a00ea
begin

	# Load forcing data

	data_force = CSV.File("../data/met_CdP_0506.csv") |> DataFrame

	input = Input{Float64}(
		data_force.year,
		data_force.month,
		data_force.day,
		data_force.hour,
		data_force.SW,
		data_force.LW,
		data_force.Sf,
		data_force.Rf,
		data_force.Ta,
		data_force.RH,
		data_force.Ua,
		data_force.Ps,
		)

	# Initilize model for experiments

	ebm = EBM{Float64}(
		am=parse(Int, am),
		cm=parse(Int, cm),
		dm=parse(Int, dm),
		em=parse(Int, em),
		hm=parse(Int, hm),
		zT=1.5,
		zvar=false,
		Tsoil=[282.98, 284.17, 284.70, 284.70]
	)

	cn = Constants{Float64}()

	hs_exp = similar(input.Ta)
	swe_exp = similar(input.Ta)
	tsurf_exp = similar(input.Ta)

	run!(ebm, cn, hs_exp, swe_exp, tsurf_exp, input)

	hs_base = similar(input.Ta)
	swe_base = similar(input.Ta)
	tsurf_base = similar(input.Ta)

	ebm = EBM{Float64}(
		am=0,
		cm=0,
		dm=0,
		em=0,
		hm=0,
		zT=1.5,
		zvar=false,
		Tsoil=[282.98, 284.17, 284.70, 284.70]
	)

	run!(ebm, cn, hs_base, swe_base, tsurf_base, input)

end

# ╔═╡ 7c6e503c-55aa-11eb-040e-69123898c27b
begin
	file_ref = "out_CdP_0506_$(am)$(cm)$(dm)$(em)$(hm).txt"

	data_ref = CSV.File(joinpath("../test/data", file_ref), header=["year", "month", "day", "hour", "alb", "Roff", "snowdepth", "SWE", "Tsurf", "Tsoil"], delim=" ", ignorerepeated=true) |> DataFrame

	hs_ref = data_ref.snowdepth
	swe_ref = data_ref.SWE
	tsurf_ref = data_ref.Tsurf
end

# ╔═╡ 1a0e63e2-4771-11eb-0e89-cdac7083e1ae
begin
	plot(hs_ref, ylabel="Snow depth (m)", linewidth=2, label="Ref")
	plot!(hs_exp, label="Exp")
	plot!(hs_base, label="Base")
end

# ╔═╡ 4a1e7488-49ff-11eb-22f1-997033d28ca0
begin
	plot(swe_ref, ylabel="Snow water equivalent (mm)", linewidth=2, label="Ref")
	plot!(swe_exp, label="Exp")
	plot!(swe_base, label="Base")
end

# ╔═╡ 062a89c0-5825-11eb-30c0-1d2e3d30e3e6
begin
	plot(tsurf_ref .+ 273.15, ylabel="Surface temperature (K)", linewidth=2, label="Ref")
	plot!(tsurf_exp, label="Exp")
	plot!(tsurf_base, label="Base")
end

# ╔═╡ 0c4d506a-4a00-11eb-00bf-e92fb70b348c

# ╔═╡ Cell order:
# ╟─a2eb46ba-4770-11eb-1d86-4d2e18f25bc1
# ╟─f126375a-4772-11eb-11f8-a71c6c86ba36
# ╟─125fc900-4772-11eb-3d6b-d121f32c31c5
# ╟─2b7b1710-49fd-11eb-072a-e320f2a80a3c
# ╟─36b71ff0-49fd-11eb-0cc8-9f7f0fdac201
# ╟─96f8c1f2-49fd-11eb-19ff-8db150ad53a3
# ╟─9fd8aa12-49fd-11eb-32cb-bdd6f962d075
# ╟─bfe40edc-49fd-11eb-0268-570ff06b4020
# ╟─c9ef08a0-49fd-11eb-07ed-252c5d4da577
# ╟─1207461a-4773-11eb-3b2e-a181f28c29ba
# ╟─1e2df4f4-4773-11eb-3f1f-dd30364bc006
# ╟─3db3b4e0-5825-11eb-22e3-b1fe7436794a
# ╟─1a0e63e2-4771-11eb-0e89-cdac7083e1ae
# ╟─4a1e7488-49ff-11eb-22f1-997033d28ca0
# ╟─062a89c0-5825-11eb-30c0-1d2e3d30e3e6
# ╟─556f9b8e-4773-11eb-2b6b-abe07b0efdc1
# ╟─980e8ac4-4771-11eb-0a52-c35a9a0a00ea
# ╟─7c6e503c-55aa-11eb-040e-69123898c27b
# ╟─dd29cb1a-4770-11eb-0db2-272859d3b045
# ╠═0c4d506a-4a00-11eb-00bf-e92fb70b348c


text:
using LinearAlgebraicRepresentation
using Plasm
using SparseArrays
Lar = LinearAlgebraicRepresentation
Lara = Lar.Arrangement

V = [
	0.0 0.4 0.5 0.6 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 0.4 0.5 0.6 1.0;
	0.0 0.0 0.0 0.0 0.0 0.4 0.4 0.5 0.5 0.6 0.6 1.0 1.0 1.0 1.0 1.0
]

EV = [
	[ 1,  2], [ 1,  6], [ 2,  6],
	[ 4,  5], [ 4,  7], [ 5,  7],
	[10, 12], [10, 13], [12, 13],
	[11, 15], [11, 16], [15, 16],
	[ 3,  8], [ 3,  9], [ 3, 14], [ 8,  9], [ 8, 14], [ 9, 14]
]

## 1 - Input Reading
Plasm.view(Plasm.numbering(0.5)((V,[[[k] for k=1:size(V,2)], EV])));

# Planar Arrangement 1
W = convert(Lar.Points, V'); # Infering type for W = V'
copEV = Lar.coboundary_0(EV::Lar.Cells);
W1, copEV1 = Lar.planar_arrangement_1(W::Lar.Points, copEV::Lar.ChainOp);

## 2 - 2-cells fragmentation
EV1 = Lar.cop2lar(copEV1);
V1 = convert(Lar.Points, W1');
Plasm.view(Plasm.numbering(0.5)((V1,[[[k] for k=1:size(V1,2)], EV1])));

# Biconnected COmponent Evaluation
bicon_comps = Lar.Arrangement.biconnected_components(copEV1);

## 3 - Biconnected Components
hpcs = [ Plasm.lar2hpc(V1,[EV1[e] for e in comp]) for comp in bicon_comps ]
Plasm.view([
	Plasm.color(Plasm.colorkey[(k%12)==0 ? 12 : k%12])(hpcs[k])
	for k = 1 : (length(hpcs))
])

# computation of 2-cells and 2-boundary

n = size(bicon_comps, 1)
shells = Array{Lar.Chain, 1}(undef, n)
boundaries = Array{Lar.ChainOp, 1}(undef, n)
EVs = Array{Lar.ChainOp, 1}(undef, n)
# for each component
for p=1:n
	ev = copEV1[sort(bicon_comps[p]), :]
	fe = Lara.minimal_2cycles(W1, ev)
	global shell_num = Lara.get_external_cycle(W1, ev, fe)
	EVs[p] = ev
	global tokeep = setdiff(1:fe.m, shell_num)
	boundaries[p] = fe[tokeep, :]
	shells[p] = fe[shell_num, :]
end

shell_bboxes = []
for i in 1:n
	vs_indexes = (abs.(EVs[i]')*abs.(shells[i])).nzind
	push!(shell_bboxes, Lar.bbox(W1[vs_indexes, :]))
end

containment_graph = Lara.pre_containment_test(shell_bboxes)
containment_graph = Lara.prune_containment_graph(
	n, W1, EVs, shells, containment_graph
)
Lara.transitive_reduction!(containment_graph)
W2 = W1;
copEV2, copFE2 = Lara.cell_merging(
	n, containment_graph, W2, EVs, boundaries, shells, shell_bboxes
)

## 4 - 3-cells identification & dangling 1-cells elimination
Plasm.view( Plasm.numbering1(0.5)((W2, copEV2, copFE2)) )

# 5 - Colorfull Representation
triangulated_faces = Lar.triangulate2D(W2, [copEV2, copFE2])
V2 = convert(Lar.Points, W2')
FVs2 = convert(Array{Lar.Cells}, triangulated_faces)
Plasm.viewcolor(V2::Lar.Points, FVs2::Array{Lar.Cells})

# polygonal face boundaries
EVs2 = Lar.FV2EVs(copEV2, copFE2)
EVs2 = convert(Array{Array{Array{Int64,1},1},1}, EVs2)
Plasm.viewcolor(V2::Lar.Points, EVs2::Array{Lar.Cells})

# 6 - Exploded Representation
model = V2,EVs2
Plasm.view(Plasm.lar_exploded(model)(1.2,1.2,1.2))


text:
using Test
x = [1, 2, 3, missing,4]
y =[missing,1, 2,3,4]
u = [missing,missing,1,2]
v = [3,4,missing,missing]

X = [1  2;
     missing  3;
     4 missing;
     missing missing;
     5 6]
Y = [1  2;
    missing  3;
    4 missing;
    missing missing;
    5 6]

@test KendallTau.skipmissingpairs(x,y) == ([2, 3, 4], [1, 2, 4])
@test KendallTau.skipmissingpairs(float.(x),y) == ([2.0, 3.0, 4.0], [1, 2,4])
@test KendallTau.skipmissingpairs(x,float.(y)) == ([2, 3, 4], [1.0, 2.0, 4.0])
@test KendallTau.skipmissingpairs(u,v) == (Int64[], Int64[])
@test KendallTau.skipmissingpairs(X) == [1 2;5 6]
@test KendallTau.skipmissingpairs(float.(X)) == [1.0 2.0;5.0 6.0]
@test KendallTau.skipmissingpairs(X,X) == ([1 2;5 6],[1 2;5 6])
@test KendallTau.skipmissingpairs(X,x) == ([1 2; 5 6], [1, 4])
@test KendallTau.skipmissingpairs(X,y) == ([5 6], [4])
@test KendallTau.skipmissingpairs(x,X) == ([1,4],[1 2; 5 6])
@test KendallTau.skipmissingpairs(y,X) == ([4],[5 6])