text:
@testset "Thinning" begin
    # for □ like figure
    img = falses(20, 20)
    for i in 7:(7 + 4), j in 7:(7 + 4)
        img[i, j] = img[j, i] = true
    end
    thin = thinning(img)
    @test count(thin .== 1) == 1
    @test thin[9, 9] == 1

    # for + like figure
    img = falses(20, 20)
    for i in 8:(8 + 5), j in 4:(13 + 4)
        img[i, j] = true
        img[j, i] = true
    end
    ans = falses(size(img))
    ans[[111, 131, 151, 171, 191, 192, 193, 194, 195, 211, 231, 251, 271]] .= 1
    @test thinning(img) == ans

    img = Bool[
        1 1 1 1 1 1
        0 0 0 0 0 0
        1 1 1 1 1 1
        1 1 1 1 1 1
        0 0 1 1 0 0
        0 0 0 0 0 0
    ]
    ans = Bool[
        1 1 1 1 1 1
        0 0 0 0 0 0
        0 0 0 0 0 0
        1 1 1 1 1 0
        0 0 0 0 0 0
        0 0 0 0 0 0
    ]
    @test thinning(img) == ans

    # already thinned
    img = Bool[
        0 0 0 0 0 0
        0 0 0 1 1 0
        0 0 1 0 0 0
        0 0 1 0 0 0
        0 0 0 1 0 0
        0 0 0 0 0 0
        0 0 0 0 0 0
    ]
    @test thinning(img) == img
end


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb02kd(datfile, io=stdout)
    ZERO = 0.0e0
    ONE = 1.0e0
    NIN = 5
    NOUT = 6
    KMAX = 20
    LMAX = 20
    MMAX = 20
    NMAX = 20
    RMAX = 20
    LDB  = LMAX*NMAX
    LDC  = KMAX*MMAX
    LDTC = MMAX*KMAX
    LDTR = KMAX
    LDWORK = 2*( KMAX*LMAX + KMAX*RMAX + LMAX*RMAX + 1 )*( MMAX + NMAX )
    TC = Array{Float64,2}(undef, LDTC,LMAX)
    TR = Array{Float64,2}(undef, LDTR,NMAX*LMAX)
    B = Array{Float64,2}(undef, LDB,RMAX)
    C = Array{Float64,2}(undef, LDC,RMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    K = parse(BlasInt, vs[1])
    L = parse(BlasInt, vs[2])
    M = parse(BlasInt, vs[3])
    N = parse(BlasInt, vs[4])
    R = parse(BlasInt, vs[5])
    LDBLK =  vs[6][1]
    TRANS =  vs[7][1]
    if ( K<=0 || K>KMAX )
        @error "Illegal K=$K"
    elseif ( L<=0 || L>LMAX )
        @error "Illegal L=$L"
    elseif ( M<=0 || M>MMAX )
        @error "Illegal M=$M"
    elseif ( N<=0 || N>NMAX )
        @error "Illegal N=$N"
    elseif ( R<=0 || R>RMAX )
        @error "Illegal R=*R"
    end
    if ( LSAME( LDBLK, 'R' ) )

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

    end # if
    if ( LSAME( TRANS, 'N' ) )

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

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

    INFO = SLICOT.mb02kd!(LDBLK, TRANS, K, L, M, N, R, ALPHA, BETA, TC, TR, B, C)
    @test INFO == 0
    INFO == 0 || return
    if ( LSAME( TRANS, 'N' ) )
        # interp output 1
    println(io, "C:")
    _nc = R
    _nr = M*K
    show(io, "text/plain", C[1:_nr,1:_nc])
    println(io)

    else
# interp output 2
    println(io, "C:")
    _nc = R
    _nr = N*L
    show(io, "text/plain", C[1:_nr,1:_nc])
    println(io)

    end # if
end # run_mb02kd()

ext: jl


text:
#=
see description in native.jl
=#

using SparseArrays

struct MatrixCompletionJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    nuclear_obj::Bool # use nuclear norm objective, else spectral norm
    symmetric::Bool # use symmetric matrix, else rectangular matrix
    use_EF::Bool # in symmetric case, construct EF for spectral/nuclear cone
    nrow::Int # number of rows in matrix
    ncol::Int # not used in symmetric case
    sparsity::Real # fraction of values that are known
end

function build(inst::MatrixCompletionJuMP{T}) where {T <: Float64}
    @assert 0 < inst.sparsity < 1
    nrow = inst.nrow
    @assert nrow >= 1
    symmetric = inst.symmetric
    nuclear_obj = inst.nuclear_obj

    if symmetric
        ncol = nrow
        len = Cones.svec_length(nrow)
        U = triu!(sprandn(T, nrow, nrow, inst.sparsity))
        (rows, cols, Avals) = findnz(U)
    else
        ncol = inst.ncol
        @assert ncol >= nrow
        len = nrow * ncol
        (rows, cols, Avals) = findnz(sprandn(T, nrow, ncol, inst.sparsity))
    end
    num_unknown = len - length(Avals)

    model = JuMP.Model()
    JuMP.@variable(model, t)
    JuMP.@objective(model, Min, t)
    JuMP.@variable(model, x[1:num_unknown])
    JuMP.@constraint(model, vcat(1, x) in MOI.GeometricMeanCone(1 + length(x)))

    X = Matrix{JuMP.AffExpr}(undef, nrow, ncol)
    for (r, c, a) in zip(rows, cols, Avals)
        X[r, c] = a
    end
    k = 1
    for c in 1:ncol, r in 1:(symmetric ? c : nrow)
        if !isassigned(X, r, c)
            X[r, c] = x[k]
            k += 1
        end
    end
    @assert k == length(x) + 1

    if symmetric
        LinearAlgebra.copytri!(X, 'U')
        if inst.use_EF
            if nuclear_obj
                # EF for symmetric nuclear norm (L1 of eigvals)
                JuMP.@variable(model, X1[1:nrow, 1:nrow], PSD)
                JuMP.@variable(model, X2[1:nrow, 1:nrow], PSD)
                eq = [X[i, j] - X1[i, j] + X2[i, j] for j in 1:nrow for i in 1:j]
                JuMP.@constraint(model, eq .== 0)
                JuMP.@constraint(model, t >= tr(X1) + tr(X2))
            else
                # EF for symmetric spectral norm (Linf of eigvals)
                tI = t * Matrix(I, nrow, nrow)
                JuMP.@constraint(model, Symmetric(tI - X) in JuMP.PSDCone())
                JuMP.@constraint(model, Symmetric(tI + X) in JuMP.PSDCone())
            end
        else
            K = Hypatia.EpiNormSpectralTriCone{T, T}(1 + len, nuclear_obj)
            Xvec = Vector{JuMP.AffExpr}(undef, len)
            Cones.smat_to_svec!(Xvec, X, sqrt(T(2)))
            JuMP.@constraint(model, vcat(t, Xvec) in K)
        end
    else
        K = (nuclear_obj ? MOI.NormNuclearCone : MOI.NormSpectralCone)(
            nrow, ncol)
        JuMP.@constraint(model, vcat(t, vec(X)) in K)
    end

    # save for use in tests
    model.ext[:x_var] = x
    model.ext[:X_var] = X

    return model
end

function test_extra(inst::MatrixCompletionJuMP{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == MOI.OPTIMAL
    (stat == MOI.OPTIMAL) || return

    # check objective and feasibility
    tol = eps(T)^0.2
    x_opt = JuMP.value.(model.ext[:x_var])
    x_geom = exp(sum(log, x_opt) / length(x_opt))
    @test x_geom <= 1 + tol
    X_opt = JuMP.value.(model.ext[:X_var])
    s = (inst.symmetric ? abs.(eigvals(Symmetric(X_opt, :U))) : svdvals(X_opt))
    snorm = (inst.nuclear_obj ? sum(s) : maximum(s))
    @test JuMP.objective_value(model) ≈ snorm atol=tol rtol=tol
    return
end


text:
using UtilitiesForMRI, LinearAlgebra, Test

nt = 256
t = Float64.(1:nt)
D = derivative1d_motionpars_linop(t, 2; pars=(true,false,true,false,true,true))
θ = zeros(Float64, nt, 6); θ[129-30:129+30,:] .= 1; θ[1:10,1] .= 1; θ[end-9:end,1] .= 1
Dθ = reshape(D*vec(θ), :, 6)

nt = 256+1
t = Float64.(0:nt-1)
ti = Float64.(0:2:nt-1); nti = length(ti)
Itp = interpolation1d_motionpars_linop((ti,ti,t,ti,ti,t), t)
θ = randn(Float64, 4*nti+2*nt)
Iθ = reshape(Itp*θ, nt, 6)

max_stars_repo_name: grizzuti/UtilitiesForMRI

max_issues_repo_name: grizzuti/UtilitiesForMRI

max_forks_repo_name: grizzuti/UtilitiesForMRI


text:
using Documenter, Gomah

makedocs(;
    modules=[Gomah],
    format=Documenter.HTML(assets=String[]),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/terasakisatoshi/Gomah.jl/blob/{commit}{path}#L{line}",
    sitename="Gomah.jl",
    authors="SatoshiTerasaki <terasakisatoshi.math@gmail.com>",
)

deploydocs(;
    repo="github.com/terasakisatoshi/Gomah.jl",
)


text:
"""
    euclsq(X, Y)
Squared Euclidean distances
between the rows of `X` and `Y`.
* `X` : Data.
* `Y` : Data.

When `X` (n, p) and `Y` (m, p), it returns an object (n, m) with:
* i, j = distance between row i of `X` and row j of `Y`.
"""
function euclsq(X, Y)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    Distances.pairwise(SqEuclidean(), X', Y', dims = 2)
end

"""
    mahsq(X, Y)
    mahsq(X, Y, Sinv)
Squared Mahalanobis distances
between the rows of `X` and `Y`.
* `X` : Data.
* `Y` : Data.
* `Sinv` : Inverse of a covariance matrix S.
    If not given, this is the uncorrected covariance matrix of `X`.

For `X` (n, p) and `Y` (m, p), it returns an object (n, m) with:
* i, j = distance between row i of `X` and row j of `Y`.
"""
function mahsq(X, Y)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    S = Statistics.cov(X, corrected = false)
    LinearAlgebra.inv!(cholesky!(Hermitian(S)))
    Distances.pairwise(SqMahalanobis(S), X', Y', dims = 2)
end

function mahsq(X, Y, Sinv)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    Sinv = ensure_mat(Sinv)
    Distances.pairwise(SqMahalanobis(Sinv), X', Y', dims = 2)
end

"""
    mahsqchol(X, Y)
    mahsqchol(X, Y, Uinv)
Compute the squared Mahalanobis distances (with a Cholesky factorization)
between the observations (rows) of `X` and `Y`.
* `X` : Data.
* `Y` : Data.
* `Uinv` : Inverse of the upper matrix of a Cholesky factorization
    of a covariance matrix S.
    If not given, the factorization is done on S, the uncorrected covariance matrix of `X`.

When `X` (n, p) and `Y` (m, p), it returns an object (n, m) with:
* i, j = distance between row i of `X` and row j of `Y`.
"""
function mahsqchol(X, Y)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    p = size(X, 2)
    S = Statistics.cov(X, corrected = false)
    if p == 1
        Uinv = inv(sqrt(S))
    else
        Uinv = LinearAlgebra.inv!(cholesky!(Hermitian(S)).U)
    end
    zX = X * Uinv
    zY = Y * Uinv
    euclsq(zX, zY)
end

function mahsqchol(X, Y, Uinv)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    Uinv = ensure_mat(Uinv)
    zX = X * Uinv
    zY = Y * Uinv
    euclsq(zX, zY)
end


text:
# Parts of this code were taken / derived from NetworkX. See LICENSE for
# licensing details.

"""Calculates the [PageRank](https://en.wikipedia.org/wiki/PageRank) of the graph
`g`. Can optionally specify a different damping factor (`α`), number of
iterations (`n`), and convergence threshold (`ϵ`). If convergence is not
reached within `n` iterations, an error will be returned.
"""
function pagerank(g::DiGraph, α=0.85, n=100, ϵ = 1.0e-6)
    A = adjacency_matrix(g,:out,Float64)
    S = vec(sum(A,1))
    S = 1./S
    S[find(S .== Inf)]=0.0
    M = A' # need a separate line due to bug #17456 in julia
    M = (Diagonal(S) * M)'
    N = nv(g)
    x = repmat([1.0/N], N)
    p = repmat([1.0/N], N)
    dangling_weights = p
    is_dangling = find(S .== 0)

    for _ in 1:n
        xlast = x
        x = α * (M * x + sum(x[is_dangling]) * dangling_weights) + (1 - α) * p
        err = sum(abs(x - xlast))
        if (err < N * ϵ)
            return x
        end
    end
    error("Pagerank did not converge after $n iterations.")
end


text:
"""
    mutation!(model, de, group)

Mutates each particle by adding random noise y ~ normal(0, σ), where σ is
a small number.

# Arguments

- `model`: model containing a likelihood function with data and priors
- `de`: differential evolution object
- `group`: a group of particles
"""
function mutation!(model, de, group)
    @unpack σ,bounds = de
    noise = Normal(0.0, σ)
    for pt in group
        # add noise to each parameter
        proposal = pt + noise
        # add Loglikelihood of prior to particle weight
        de.evaluate_fitness!(de, model, proposal)
        # update whether proposal was accepted and add proposal if accepted
        de.update_particle!(de, pt, proposal)
    end
    return nothing
end


text:
export
    get_debug_file_id,
    reset_debug_file_id!,
    validate

DEBUG_ID_COUNTER = 0
get_debug_file_id() = Int(global DEBUG_ID_COUNTER += 1)
function reset_debug_file_id!()
    global DEBUG_ID_COUNTER = 0
end

function validate end

export
    ProblemSpec

"""
    ProblemSpec{G}

Encodes the the relaxed PC-TAPF problem that ignores collision-avoidance
constraints.

Elements:
- D::T - a distance matrix (or environment) that implements get_distance(D,x1,x2,args...)
- cost_function::F - the optimization objective (default is SumOfMakeSpans)
- Δt_collect::Dict{ObjectID,Int} # duration of COLLECT operations
- Δt_deposit::Dict{ObjectID,Int} # duration of DEPOSIT operations
"""
@with_kw struct ProblemSpec{T,F}
    D::T                            = zeros(0,0) # environment or distance matrix
    cost_function::F                = SumOfMakeSpans()
    Δt_collect::Dict{ObjectID,Int}  = Dict{ObjectID,Int}() # duration of COLLECT operations
    Δt_deposit::Dict{ObjectID,Int}  = Dict{ObjectID,Int}() # duration of DELIVER operations
end
GraphUtils.get_distance(spec::ProblemSpec,args...) = get_distance(spec.D,args...)
GraphUtils.get_distance(spec::ProblemSpec,a::LocationID,b::LocationID,args...) = get_distance(spec.D,get_id(a),get_id(b),args...)

export
    ProjectSpec,
    add_operation!,
    set_initial_condition!,
    set_final_condition!,
    get_initial_condition,
    get_final_condition,
    get_duration_vector,
    construct_operation

"""
    ProjectSpec <: AbstractCustomNDiGraph{Union{OBJECT_AT,Operation},AbstractID}

Encodes a set of operations and the prescribed initial and final locations
of the objects that form the inputs and outputs of those operations.

Elements:
- all standard fields of a concrete `GraphUtils.CustomNDiGraph` type
- initial_conditions::Dict{ObjectID,OBJECT_AT}
- final_conditions::Dict{ObjectID,OBJECT_AT}
"""
@with_kw struct ProjectSpec <: AbstractCustomNDiGraph{Union{OBJECT_AT,Operation},AbstractID}
    graph               ::DiGraph               = DiGraph()
    nodes               ::Vector{Union{OBJECT_AT,Operation}} = Vector{Union{OBJECT_AT,Operation}}()
    vtx_map             ::Dict{AbstractID,Int}  = Dict{AbstractID,Int}()
    vtx_ids             ::Vector{AbstractID}    = Vector{AbstractID}()
    initial_conditions  ::Dict{ObjectID,OBJECT_AT} = Dict{ObjectID,OBJECT_AT}()
    final_conditions    ::Dict{ObjectID,OBJECT_AT} = Dict{ObjectID,OBJECT_AT}()
end
get_initial_conditions(spec::ProjectSpec)       = spec.initial_conditions
get_final_conditions(spec::ProjectSpec)         = spec.final_conditions
get_initial_condition(spec::ProjectSpec,id::ObjectID)   = get_initial_conditions(spec)[id]
get_final_condition(spec::ProjectSpec,id::ObjectID)     = get_final_conditions(spec)[id]
initial_conditions_vector(spec::ProjectSpec)    = sort(collect(values(spec.initial_conditions));by=get_object_id)
final_conditions_vector(spec::ProjectSpec)      = sort(collect(values(spec.final_conditions));by=get_object_id)
function set_initial_condition!(spec,object_id,pred)
    spec.initial_conditions[object_id] = pred
end
function set_final_condition!(spec,object_id,pred)
    spec.final_conditions[object_id] = pred
end
for op in [:set_initial_condition!,:set_final_condition!]
    @eval $op(spec,pred) = $op(spec,get_object_id(pred),pred)
end
function ProjectSpec(ics::V,fcs::V) where {V<:Vector{OBJECT_AT}}
    spec = ProjectSpec(
        initial_conditions=Dict{ObjectID,OBJECT_AT}(get_object_id(p)=>p for p in ics),
        final_conditions=Dict{ObjectID,OBJECT_AT}(get_object_id(p)=>p for p in fcs),
        )
    for (k,o) in get_initial_conditions(spec)
        add_node!(spec,o,get_object_id(o))
        set_initial_condition!(spec,get_object_id(o),o)
    end
    return spec
end
get_pre_deps(spec::ProjectSpec, op_id::OperationID) = Set{ObjectID}(map(v->get_vtx_id(spec,v),inneighbors(spec,op_id))) #get(spec.pre_deps, op_id, Set{Int}())
get_post_deps(spec::ProjectSpec, op_id::OperationID) = Set{ObjectID}(map(v->get_vtx_id(spec,v),outneighbors(spec,op_id))) #get(spec.pre_deps, op_id, Set{Int}())
get_operations(spec::ProjectSpec) = Vector{Operation}(filter(n->isa(n,Operation), get_nodes(spec)))
get_num_delivery_tasks(spec::ProjectSpec) = length(get_nodes_of_type(spec,ObjectID))
function get_terminal_operation_ids(spec::ProjectSpec)
    ops = Set{OperationID}()
    for v in get_all_terminal_nodes(spec)
        id = get_vtx_id(spec,v)
        if isa(id,OperationID)
            push!(ops,id)
        end
    end
    return ops
end
function get_pickup_wait_time(spec::ProjectSpec,o::ObjectID)
    if !is_root_node(spec,o)
        parent = get_node(spec,inneighbors(spec,o)[1])
        @assert isa(parent,Operation)
        return duration(parent)
    end
    return 0.0
end
"""
    get_duration_vector(spec::ProjectSpec)

Return a vector `Δt` such that `Δt[i]` is the amount of time that must elapse
before object `i` can be picked up after its parent operation is performed.
"""
function get_duration_vector(spec::ProjectSpec)
    # Δt = Dict{ObjectID,Float64}()
    # for (v,n) in enumerate(get_nodes(spec))
    #     if isa(n,OBJECT_AT)
    #         Δt[get_object_id(n)] = get_pickup_wait_time(spec,get_object_id(n))
    #     end
    # end
    # Δt
    map(n->get_pickup_wait_time(spec,get_object_id(n)), initial_conditions_vector(spec))
end

"""
    construct_operation(spec::ProjectSpec, station_id, input_ids, output_ids, Δt, id=get_unique_operation_id())
"""
function construct_operation(
        spec::ProjectSpec,
        station_id::LocationID,
        input_ids::Vector{ObjectID},
        output_ids::Vector{ObjectID},
        Δt,
        id=get_unique_operation_id()
        )
    op = Operation(
        pre = Dict{ObjectID,OBJECT_AT}(map(id->id=>get_final_condition(spec,id), input_ids)),
        post = Dict{ObjectID,OBJECT_AT}(map(id->id=>get_initial_condition(spec,id), output_ids)),
        Δt = Δt,
        station_id = LocationID(station_id),
        id = id
    )
end
function construct_operation(
        spec::ProjectSpec,
        station_id::Int,
        input_ids::Vector,
        output_ids::Vector,args...)
    construct_operation(spec,
        LocationID(station_id),
        Vector{ObjectID}(convert.(ObjectID,input_ids)),
        Vector{ObjectID}(convert.(ObjectID,output_ids)),
        args...
        )
end
function add_operation!(spec::ProjectSpec, op::Operation)
    add_node!(spec,op,get_operation_id(op))
    op_id = get_operation_id(op)
    for (object_id,pred) in preconditions(op)
        if !has_vertex(spec,object_id)
            add_node!(spec,pred,object_id)
        end
        add_edge!(spec, object_id, op_id)
        set_final_condition!(spec,object_id,pred)
    end
    for (object_id,pred) in postconditions(op)
        if !has_vertex(spec,object_id)
            add_node!(spec,pred,object_id)
        else
            # want the graph nodes to represent the initial conditions of the objects
            replace_node!(spec,pred,object_id)
        end
        add_edge!(spec, op_id, object_id)
        set_initial_condition!(spec,object_id,pred)
    end
    # empty!(spec.terminal_vtxs)
    # union!(spec.terminal_vtxs, get_all_terminal_nodes(spec.graph))
    validate(spec,op)
    spec
end

export
    read_operation,
    read_project_spec

parse_object(pred::OBJECT_AT,dict=Dict{String,Any}()) = merge!(dict,Dict("type"=>"OBJECT_AT","data"=>[get_id(get_object_id(pred)),get_id(get_location_id(pred))]))
TOML.parse(pred::OBJECT_AT) = parse_object(pred)
read_object(dict) = OBJECT_AT(dict["data"]...)
read_object(arr::Vector) = OBJECT_AT(arr...)
function TOML.parse(op::Operation)
    toml_dict = Dict()
    toml_dict["pre"] = map(pred->TOML.parse(pred), collect(values(preconditions(op))))
    toml_dict["post"] = map(pred->TOML.parse(pred), collect(values(postconditions(op))))
    toml_dict["dt"] = duration(op)
    toml_dict["station_id"] = get_id(op.station_id)
    toml_dict["id"] = get_id(op.id)
    return toml_dict
end
function read_op_dict_from_array(toml_dict)
    dict = Dict{ObjectID,OBJECT_AT}()
    for arr in toml_dict
        o = read_object(arr)
        dict[get_object_id(o)] = o
    end
    return dict
end
function read_operation(toml_dict::Dict,keep_id=false)
    op_id = OperationID(get(toml_dict,"id",-1))
    if get_id(op_id) == -1 || keep_id == false
        op_id = get_unique_operation_id()
    end
    op = Operation(
        pre     = read_op_dict_from_array(toml_dict["pre"]),
        post    = read_op_dict_from_array(toml_dict["post"]),
        Δt      = get(toml_dict,"dt",get(toml_dict,"Δt",0.0)),
        station_id = LocationID(get(toml_dict,"station_id",-1)),
        id = op_id
        )
end
function TOML.parse(project_spec::ProjectSpec)
    toml_dict = Dict()
    toml_dict["operations"] = map(op->TOML.parse(op),get_operations(project_spec))
    toml_dict["initial_conditions"] = map(pred->TOML.parse(pred), initial_conditions_vector(project_spec))
    toml_dict["final_conditions"] = map(pred->TOML.parse(pred), final_conditions_vector(project_spec))
    toml_dict
end
function read_project_spec(toml_dict::Dict)
    project_spec = ProjectSpec()
    ics = toml_dict["initial_conditions"]
    fcs = toml_dict["final_conditions"]
    ics = isa(ics,Dict) ? collect(values(ics)) : ics
    fcs = isa(fcs,Dict) ? collect(values(fcs)) : fcs
    for arr in ics
        o = read_object(arr)
        set_initial_condition!(project_spec,get_object_id(o),o)
    end
    for arr in fcs
        o = read_object(arr)
        set_final_condition!(project_spec,get_object_id(o),o)
    end
    for op_dict in toml_dict["operations"]
        op = read_operation(op_dict)
        add_operation!(project_spec, op)
    end
    return project_spec
end
function read_project_spec(io)
    read_project_spec(TOML.parsefile(io))
end

export
    SimpleProblemDef,
    read_problem_def

@with_kw struct SimpleProblemDef
    project_spec::ProjectSpec       = ProjectSpec()
    r0::Vector{Int}                 = Int[]
    s0::Vector{Int}                 = Int[]
    sF::Vector{Int}                 = Int[]
    shapes::Vector{Tuple{Int,Int}}  = Vector{Tuple{Int,Int}}([(1,1) for o in s0])
end
SimpleProblemDef(project_spec,r0,s0,sF) = SimpleProblemDef(project_spec=project_spec,r0=r0,s0=s0,sF=sF)

function TOML.parse(def::SimpleProblemDef)
    toml_dict = TOML.parse(def.project_spec)
    toml_dict["r0"] = def.r0
    toml_dict["s0"] = def.s0
    toml_dict["sF"] = def.sF
    toml_dict["shapes"] = map(s->[s...], def.shapes)
    toml_dict
end
function read_problem_def(toml_dict::Dict)
    SimpleProblemDef(
        read_project_spec(toml_dict),
        toml_dict["r0"],
        toml_dict["s0"],
        toml_dict["sF"],
        # map(s->tuple(s...), toml_dict["shapes"])
        map(s->tuple(s...), get(toml_dict,"shapes",[[1,1] for o in toml_dict["s0"]]) )
    )
end
function read_problem_def(io)
    read_problem_def(TOML.parsefile(io))
end

# export
#     DeliveryTask,
#     DeliveryGraph

# """
#     DeliveryTask
# """
# struct DeliveryTask
#     o::ObjectID
#     s1::LocationID
#     s2::LocationID
# end

# """
#     DeliveryGraph{G}
# """
# struct DeliveryGraph{G}
#     tasks::Vector{DeliveryTask}
#     graph::G
# end

# export
#     construct_delivery_graph

# """
#     construct_delivery_graph(project_spec::ProjectSpec,M::Int)

# Assumes that initial station ids correspond to object ids.
# """
# function construct_delivery_graph(project_spec::ProjectSpec,M::Int)
#     delivery_graph = DeliveryGraph(Vector{DeliveryTask}(),MetaDiGraph(M))
#     object_id_map = Dict{ObjectID,Int}(
#         get_object_id(o)=>i for (i,o) in enumerate(initial_conditions_vector(project_spec)))
#     # for op in project_spec.operations
#     for op in get_operations(project_spec)
#         for (id,pred) in preconditions(op)
#             idx = object_id_map[id]
#             pre = get_initial_condition(project_spec,id)
#             post = get_final_condition(project_spec,id)
#             push!(delivery_graph.tasks,DeliveryTask(id,get_location_id(pre),get_location_id(post)))
#             for j in get_output_ids(op)
#                 idx2 = object_id_map[j]
#                 add_edge!(delivery_graph.graph, idx, idx2)
#             end
#         end
#     end
#     delivery_graph
# end

export
    PathSpec,
    generate_path_spec,
    get_t0,
    get_tF,
    set_t0!,
    set_tF!,
    get_slack,
    set_slack!,
    get_local_slack,
    set_local_slack!,
    get_min_duration,
    set_min_duration!

"""
    PathSpec

Encodes information about the path that must be planned for a particular
schedule node.

Fields:
* `node_type::Symbol = :EMPTY`
* `start_vtx::Int = -1`
* `final_vtx::Int = -1`
* `min_duration::Int = 0`
* `agent_id::Int = -1`
* `object_id::Int = -1`
* `plan_path::Bool = true` - flag indicating whether a path must be planned.
    For example, `Operation` nodes do not require any path planning.
* `tight::Bool = false` - if true, the path may not terminate prior to the
    beginning of successors. If `tight == true`, local slack == 0. For example,
    `GO` must not end before `COLLECT` can begin, because this would produce
    empty time between planning phases.
* `static::Bool = false` - if true, the robot must remain in place for this
    planning phase (e.g., COLLECT, DEPOSIT).
* `free::Bool = false` - if true, and if the node is a terminal node, the
    planning must go on until all non-free nodes are completed.
* `fixed::Bool = false` - if true, do not plan path because it is already fixed.
    Instead, retrieve the portion of the path directly from the pre-existing
    solution.
"""
@with_kw mutable struct PathSpec
    # temporal
    t0              ::Int               = 0
    min_duration    ::Int               = 0
    tF              ::Int               = t0 + min_duration
    slack           ::Vector{Float64}   = Float64[]
    local_slack     ::Vector{Float64}   = Float64[]
    # instructions
    plan_path       ::Bool              = true
    tight           ::Bool              = false
    static          ::Bool              = false
    free            ::Bool              = false
    fixed           ::Bool              = false
end
get_t0(spec::PathSpec) = spec.t0
set_t0!(spec::PathSpec,val) = begin spec.t0 = val end
get_tF(spec::PathSpec) = spec.tF
set_tF!(spec::PathSpec,val) = begin spec.tF = val end
get_min_duration(spec::PathSpec) = spec.min_duration
set_min_duration!(spec::PathSpec,val) = begin spec.min_duration = val end
get_slack(spec::PathSpec) = spec.slack
set_slack!(spec::PathSpec,val) = begin spec.slack = val end
get_local_slack(spec::PathSpec) = spec.local_slack
set_local_slack!(spec::PathSpec,val) = begin spec.local_slack = val end

const path_spec_accessor_interface = [
    :get_t0,
    :get_tF,
    :get_slack,
    :get_local_slack,
    :get_min_duration,
    ]
const path_spec_mutator_interface = [
    :set_min_duration!,
    :set_t0!,
    :set_slack!,
    :set_local_slack!,
    :set_tF!,
    ]

export
    ScheduleNode,
    get_path_spec,
    set_path_spec!

"""
    ScheduleNode{I<:AbstractID,V<:AbstractPlanningPredicate}

The node type of the `OperatingSchedule` graph.
"""
mutable struct ScheduleNode{I<:AbstractID,V<:AbstractPlanningPredicate}
    id::I
    node::V
    spec::PathSpec
end
ScheduleNode(id,node) = ScheduleNode(id,node,PathSpec())
get_path_spec(node::ScheduleNode) = node.spec
function set_path_spec!(node::ScheduleNode,spec)
    node.spec = spec
end
for op in path_spec_accessor_interface
    @eval $op(node::ScheduleNode) = $op(get_path_spec(node))
end
for op in path_spec_mutator_interface
    @eval $op(node::ScheduleNode,val) = $op(get_path_spec(node),val)
end
for op in predicate_accessor_interface
    @eval $op(node::ScheduleNode) = $op(node.node)
end
for op in [:(GraphUtils.matches_template)]
    @eval $op(template::Type{T},node::ScheduleNode) where {T} = $op(template,node.node)
end

const schedule_node_accessor_interface = [
    path_spec_accessor_interface...,
    :get_path_spec]
const schedule_node_mutator_interface = [
    path_spec_mutator_interface...,
    :set_path_spec!]

matches_node_type(n::ScheduleNode,b::Type{B}) where {B} = matches_node_type(n.node,b)

export
    OperatingSchedule,
    get_object_ICs,
    get_object_FCs,
    get_robot_ICs,
    get_actions,
    get_operations,
    get_vtx_ids,
    get_node_from_id,
    get_node_from_vtx,
    get_completion_time,
    get_duration,
    get_vtx,
    get_vtx_id,
    replace_in_schedule!,
    add_to_schedule!,
    process_schedule

"""
    OperatingSchedule

Encodes discrete events/activities that need to take place, and the precedence
constraints between them. Each `ScheduleNode` has a corresponding vertex index
and an `AbstractID`. An edge from node1 to node2 indicates a precedence
constraint between them.
"""
@with_kw struct OperatingSchedule <: AbstractCustomNDiGraph{ScheduleNode,AbstractID}
    graph               ::DiGraph               = DiGraph()
    nodes               ::Vector{ScheduleNode}  = Vector{ScheduleNode}()
    vtx_map             ::Dict{AbstractID,Int}  = Dict{AbstractID,Int}()
    vtx_ids             ::Vector{AbstractID}    = Vector{AbstractID}() # maps vertex uid to actual graph node
    terminal_vtxs       ::Vector{Int}           = Vector{Int}() # list of "project heads"
    weights             ::Dict{Int,Float64}     = Dict{Int,Float64}() # weights corresponding to project heads
end
get_terminal_vtxs(sched::P) where {P<:OperatingSchedule}     = sched.terminal_vtxs
get_root_node_weights(sched::P) where {P<:OperatingSchedule} = sched.weights

GraphUtils.get_vtx(sched::OperatingSchedule,node::ScheduleNode) = get_vtx(sched,node.id)
get_node_from_id(sched::OperatingSchedule,id)                   = get_node(sched,id).node
get_node_from_vtx(sched::OperatingSchedule,v)                   = get_node(sched,v).node

get_object_ICs(sched::P) where {P<:OperatingSchedule}  = Dict(k=>v.node for (k,v) in get_nodes_of_type(sched,ObjectID))
get_robot_ICs(sched::P) where {P<:OperatingSchedule}   = Dict(k=>v.node for (k,v) in get_nodes_of_type(sched,BotID))
get_actions(sched::P) where {P<:OperatingSchedule}     = Dict(k=>v.node for (k,v) in get_nodes_of_type(sched,ActionID))
get_operations(sched::P) where {P<:OperatingSchedule}  = Dict(k=>v.node for (k,v) in get_nodes_of_type(sched,OperationID))

for op in schedule_node_accessor_interface
    @eval $op(sched::OperatingSchedule,v) = $op(get_node(sched,v))
    @eval $op(sched::OperatingSchedule) = map(v->$op(get_node(sched,v)), vertices(sched))
end
for op in schedule_node_mutator_interface
    @eval $op(sched::OperatingSchedule,v,val) = $op(get_node(sched,v),val)
    @eval $op(sched::OperatingSchedule,val) = begin
        for v in vertices(sched)
            $op(get_node(sched,v),val)
        end
    end
end

# Defining default node parameters
is_tight(p) = false
is_tight(::BOT_GO) = true
is_free(p) = false
is_free(p::BOT_AT) = true
is_free(p::BOT_GO) = !CRCBS.is_valid(get_destination_location_id(p))
is_static(p) = false
is_static(p::Union{BOT_COLLECT,BOT_DEPOSIT}) = true
for op in [:is_free,:is_static,:is_tight]
    @eval $op(p::TEAM_ACTION) = $op(sub_nodes(p)[1])
end
needs_path(p) = false
needs_path(p::Union{BOT_AT,AbstractRobotAction}) = true
duration_lower_bound(args...) = 0
duration_lower_bound(op::Operation,spec) = duration(op)
duration_lower_bound(p::BOT_COLLECT,spec) = get(spec.Δt_collect,get_object_id(p),0)
duration_lower_bound(p::BOT_DEPOSIT,spec) = get(spec.Δt_deposit,get_object_id(p),0)
function duration_lower_bound(p::Union{BOT_AT,AbstractRobotAction},spec)
    x1 = get_initial_location_id(p)
    x2 = get_destination_location_id(p)
    return get_distance(spec,x1,x2)
end
function duration_lower_bound(p::TEAM_ACTION,spec)
    x1 = get_initial_location_id(p.instructions[1])
    x2 = get_destination_location_id(p.instructions[1])
    return get_distance(spec,x1,x2,team_configuration(p))
end

"""
    generate_path_spec(spec,node)

Generates a `PathSpec` struct that encodes information about the path to be
planned for `node`.

Arguments:
* spec::ProblemSpec
* node::T <: AbstractPlanningPredicate
"""
function generate_path_spec(spec::ProblemSpec,a)
    PathSpec(
        min_duration=duration_lower_bound(a,spec),
        tight=is_tight(a),
        static=is_static(a),
        free=is_free(a),
        plan_path=needs_path(a)
        )
end
generate_path_spec(sched::OperatingSchedule,spec::ProblemSpec,pred) = generate_path_spec(spec,pred)
generate_path_spec(sched::OperatingSchedule,pred) = generate_path_spec(sched,ProblemSpec(),pred)

"""
    replace_in_schedule!(schedule::OperatingSchedule,path_spec::T,pred,id::ID) where {T<:PathSpec,ID<:AbstractID}

Replace the `ScheduleNode` associated with `id` with the new node `pred`, and
the accompanying `PathSpec` `path_spec`.
"""
replace_in_schedule!(sched::OperatingSchedule,node::ScheduleNode,id::AbstractID=node.id) = replace_node!(sched,node,id)
function replace_in_schedule!(sched::P,path_spec::T,pred,id::ID) where {P<:OperatingSchedule,T<:PathSpec,ID<:AbstractID}
    replace_in_schedule!(sched,ScheduleNode(id,pred,path_spec))
end
function replace_in_schedule!(sched::P,spec::T,pred,id::ID) where {P<:OperatingSchedule,T<:ProblemSpec,ID<:AbstractID}
    replace_in_schedule!(sched,generate_path_spec(sched,spec,pred),pred,id)
end
function replace_in_schedule!(sched::P,pred,id::ID) where {P<:OperatingSchedule,ID<:AbstractID}
    replace_in_schedule!(sched,ProblemSpec(),pred,id)
end

GraphUtils.add_node!(sched::OperatingSchedule,node::ScheduleNode) = add_node!(sched,node,node.id)

function GraphUtils.make_node(g::OperatingSchedule,pred::AbstractPlanningPredicate,args...)
   make_node(g,PathSpec(),pred,args...)
end
function GraphUtils.make_node(g::OperatingSchedule,spec::ProblemSpec,pred,args...)
   make_node(g,generate_path_spec(g,spec,pred),pred,args...)
end
function GraphUtils.make_node(g::OperatingSchedule,spec::PathSpec,pred::BOT_AT,id=get_robot_id(pred))
    ScheduleNode(id,pred,spec)
end
function GraphUtils.make_node(g::OperatingSchedule,spec::PathSpec,pred::OBJECT_AT,id=get_object_id(pred))
    ScheduleNode(id,pred,spec)
end
function GraphUtils.make_node(g::OperatingSchedule,spec::PathSpec,pred::AbstractRobotAction,id=get_unique_action_id())
    ScheduleNode(id,pred,spec)
end
function GraphUtils.make_node(g::OperatingSchedule,spec::PathSpec,pred::Operation,id=get_operation_id(pred))
    ScheduleNode(id,pred,spec)
end

export
    get_leaf_operation_vtxs,
    set_leaf_operation_vtxs!

function get_leaf_operation_vtxs(sched::OperatingSchedule)
    terminal_vtxs = Int[]
    for v in get_all_terminal_nodes(sched)
        if matches_node_type(get_node(sched,v),Operation)
            push!(terminal_vtxs,v)
        end
    end
    return terminal_vtxs
end
function set_leaf_operation_vtxs!(sched::OperatingSchedule)
    empty!(get_terminal_vtxs(sched))
    empty!(get_root_node_weights(sched))
    for vtx in get_leaf_operation_vtxs(sched)
        push!(get_terminal_vtxs(sched),vtx)
        get_root_node_weights(sched)[vtx] = 1.0
    end
    sched
end

export backtrack_node

"""
    `backtrack_node(sched::OperatingSchedule,v::Int)`

Find the closest ancestor of `v` with overlapping `RobotID`s.
"""
function backtrack_node(sched::OperatingSchedule,v::Int)
    robot_ids = get_valid_robot_ids(sched,get_vtx_id(sched,v),v)
    vtxs = Int[]
    if isempty(robot_ids)
        return vtxs
    end
    for vp in inneighbors(sched,v)
        if !isempty(intersect(get_valid_robot_ids(sched,vp),robot_ids))
            push!(vtxs,vp)
        end
    end
    return vtxs
end

export
    validate,
    sanity_check

function validate(path::Path, msg::String)
    @assert( !any(convert_to_vertex_lists(path) .== -1), msg )
    return true
end
function validate(path::Path, v::Int)
    validate(path, string("invalid path with -1 vtxs: v = ",v,", path = ",convert_to_vertex_lists(path)))
end
function validate(path::Path, v::Int, cbs_env)
    validate(path, string("v = ",v,", path = ",convert_to_vertex_lists(path),", goal: ",cbs_env.goal))
end
function sanity_check(sched::OperatingSchedule,append_string="")
    try
        @assert !is_cyclic(sched) "is_cyclic(sched)"
        for v in vertices(sched)
            node = get_node_from_vtx(sched, v)
            id = get_vtx_id(sched,v)
            if matches_node_type(node,COLLECT)
                @assert(get_location_id(node) != -1, string("get_location_id(node) != -1 for node id ", id))
            end
            if matches_node_type(node,Operation)
                # input_ids = Set(map(o->get_id(get_object_id(o)),collect(preconditions(node)))))
                input_ids = Set(collect(keys(preconditions(node))))
                for v2 in inneighbors(sched,v)
                    node2 = get_node_from_vtx(sched, v2)
                    @assert(get_object_id(node2) in input_ids, string(string(node2), " should not be an inneighbor of ",string(node), " whose inputs should be ",input_ids))
                end
                @assert(
                    indegree(sched,v) == length(preconditions(node)),
                    string("Operation ",string(node),
                        " needs edges from objects ",collect(input_ids),
                        " but only has edges from objects ",map(v2->get_id(get_object_id(get_node_from_vtx(sched, v2))),inneighbors(sched,v))
                    )
                    )
            end
        end
    catch e
        if typeof(e) <: AssertionError
            bt = catch_backtrace()
            showerror(stdout,e,bt)
            print(string(e.msg, append_string))
        else
            rethrow(e)
        end
        return false
    end
    return true
end
function validate(spec::ProjectSpec,op)
    for (id,pred) in preconditions(op)
        @assert get_final_condition(spec,id) == pred
    end
    for (id,pred) in postconditions(op)
        @assert get_initial_condition(spec,id) == pred
    end
    return true
end
function validate(sched::OperatingSchedule)
    try
        @assert !is_cyclic(sched) "is_cyclic(G)"
        for e in edges(sched)
            node1 = get_node_from_id(sched, get_vtx_id(sched, e.src))
            node2 = get_node_from_id(sched, get_vtx_id(sched, e.dst))
            @assert(validate_edge(node1,node2), string(" INVALID EDGE: ", string(node1), " --> ",string(node2)))
        end
        for v in vertices(sched)
            id = get_vtx_id(sched, v)
            node = get_node_from_id(sched, id)
            if matches_node_type(node,COLLECT)
                @assert(get_location_id(node) != -1, string("get_location_id(node) != -1 for node id ", id))
            end
            @assert( outdegree(sched,v) >= sum([0, values(required_successors(node))...]) , string("node = ", string(node), " outdegree = ",outdegree(G,v), " "))
            @assert( indegree(sched,v) >= sum([0, values(required_predecessors(node))...]), string("node = ", string(node), " indegree = ",indegree(G,v), " ") )
            if matches_node_type(node, AbstractSingleRobotAction)
                for v2 in outneighbors(sched,v)
                    node2 = get_node_from_vtx(sched, v2)
                    if matches_node_type(node2, AbstractSingleRobotAction)
                        if length(intersect(resources_reserved(node),resources_reserved(node2))) == 0 # job shop constraint
                            @assert( get_robot_id(node) == get_robot_id(node2), string("robot IDs do not match: ",string(node), " --> ", string(node2)))
                        end
                    end
                end
            end
        end
    catch e
        if typeof(e) <: AssertionError
            bt = catch_backtrace()
            showerror(stdout,e,bt)
            print(e.msg)
        else
            rethrow(e)
        end
        return false
    end
    return true
end
function validate(node::ScheduleNode,paths::Vector{Vector{Int}})
    if matches_template(TEAM_ACTION,node)
        for n in sub_nodes(node)
            sub_node = ScheduleNode(node.id,n,node.spec)
            if !validate(sub_node,paths)
                return false
            end
        end
    else
        agent_id = get_id(get_default_robot_id(node))
        if agent_id != -1
            path = paths[agent_id]
            start_vtx = get_id(get_initial_location_id(node))
            final_vtx = get_id(get_destination_location_id(node))
            try
                @assert(length(path) > get_t0(node), string("length(path) == $(length(path)), should be greater than get_t0(node) == $(get_t0(node)) in node ",string(node)))
                @assert(length(path) > get_tF(node), string("length(path) == $(length(path)), should be greater than get_t0(node) == $(get_t0(node)) in node ",string(node)))
                if start_vtx != -1
                    if length(path) > get_t0(node)
                        @assert(path[get_t0(node) + 1] == start_vtx, string("node: ",string(node), ", start_vtx: ",start_vtx, ", t0+1: ",get_t0(node)+1,", path[get_t0(node) + 1] = ",path[get_t0(node) + 1],", path: ", path))
                    end
                end
                if final_vtx != -1
                    if length(path) > get_tF(node)
                        @assert(path[get_tF(node) + 1] == final_vtx, string("node: ",string(node), ", final vtx: ",final_vtx, ", tF+1: ",get_tF(node)+1,", path[get_tF(node) + 1] = ",path[get_tF(node) + 1],", path: ", path))
                    end
                end
            catch e
                if typeof(e) <: AssertionError
                    print(e.msg)
                else
                    throw(e)
                end
                return false
            end
        end
    end
    return true
end
function validate(sched::OperatingSchedule,paths::Vector{Vector{Int}})
    for v in vertices(sched)
        node = get_node(sched,v)
        if !validate(node,paths)
            return false
        end
    end
    return true
end

export
    add_single_robot_delivery_task!

function add_headless_delivery_task!(
        sched::OperatingSchedule,
        problem_spec::ProblemSpec,
        object_id::ObjectID,
        op_id::OperationID,
        robot_type::Type{R}=DeliveryBot,
        pickup_station_id::LocationID=get_location_id(get_node_from_id(sched,object_id)),
        dropoff_station_id::LocationID=get_dropoff(get_node_from_id(sched,op_id),object_id),
        ;
        t0::Int=0,
    ) where {R<:AbstractRobotType}

    robot_id = BotID{R}(-1)
    prev_id = object_id
    for p in [
        BOT_COLLECT(robot_id, object_id, pickup_station_id),
        BOT_CARRY(robot_id, object_id, pickup_station_id, dropoff_station_id),
        BOT_DEPOSIT(robot_id, object_id, dropoff_station_id),
        BOT_GO(robot_id, dropoff_station_id,LocationID(-1)),
        ]
        node = add_node!(sched, make_node(sched, problem_spec, p))
        set_t0!(sched,node.id,t0)
        add_edge!(sched, prev_id, node.id)
        if matches_node_type(node,BOT_DEPOSIT)
            add_edge!(sched, node.id, op_id)
        end
        prev_id = node.id
    end

    return
end
function add_headless_team_delivery_task!(
        sched::OperatingSchedule,
        problem_spec::ProblemSpec,
        object_id::ObjectID,
        op_id::OperationID,
        robot_type::Type{R}=DeliveryBot,
        pickup_station_ids::Vector{LocationID}=get_location_ids(get_node_from_id(sched,object_id)),
        dropoff_station_ids::Vector{LocationID}=get_dropoffs(get_node_from_id(sched,op_id),object_id)
        ) where {R<:AbstractRobotType}

    robot_id = BotID{R}(-1)
    @assert length(pickup_station_ids) == length(dropoff_station_ids)
    n = length(pickup_station_ids)

    object_node = get_node_from_id(sched,object_id)
    shape = object_node.shape
    prev_id = object_id
    for team_action in [
            TEAM_COLLECT(
            instructions = map(x->BOT_COLLECT(robot_id, object_id, x), pickup_station_ids),
            shape=shape
            ),
            TEAM_CARRY(
            instructions = [BOT_CARRY(robot_id, object_id, x1, x2) for (x1,x2) in zip(pickup_station_ids,dropoff_station_ids)],
            shape = shape
            ),
            TEAM_DEPOSIT(
            instructions = map(x->BOT_DEPOSIT(robot_id, object_id, x), dropoff_station_ids),
            shape=shape
            )
        ]
        node = add_node!(sched, make_node(sched,problem_spec,team_action))
        add_edge!(sched,prev_id,node.id)
        prev_id = node.id
        if matches_node_type(node,TEAM_COLLECT)
            for x in pickup_station_ids
                subnode = add_node!(sched,make_node(sched,problem_spec,BOT_GO(robot_id,x,x)))
                add_edge!(sched,subnode.id,node.id)
            end
        elseif matches_node_type(node,TEAM_DEPOSIT)
            add_edge!(sched,node.id,op_id)
            for x in dropoff_station_ids
                subnode = add_node!(sched,make_node(sched,problem_spec,BOT_GO(robot_id,x,-1)))
                add_edge!(sched, node.id, subnode.id)
            end
        end
    end
    return
end
function add_single_robot_delivery_task!(
        sched::OperatingSchedule,
        spec::ProblemSpec,
        robot_id::BotID{R},
        object_id::ObjectID,
        op_id::OperationID,
        go_id::ActionID,
        x1::LocationID=get_location_id(get_node_from_id(sched,object_id)),
        x2::LocationID=get_dropoff(get_node_from_id(sched,op_id),object_id),
        args...
    ) where {R<:AbstractRobotType}
    add_headless_delivery_task!(sched,spec,object_id,op_id,R,x1,x2)
    collect_id = get_vtx_id(sched,outneighbors(sched,get_vtx(sched,object_id))[1])
    add_edge!(sched,go_id,collect_id)
    propagate_valid_ids!(sched,spec)
end

export
    construct_partial_project_schedule

function add_new_robot_to_schedule!(sched,pred::BOT_AT,problem_spec=ProblemSpec())
    n1 = add_node!(sched, make_node(sched,problem_spec,pred)) #, pred, robot_id)
    action = BOT_GO(r=get_robot_id(pred),x1=get_location_id(pred))
    n2 = add_node!(sched, make_node(sched,problem_spec,action)) #, action, action_id)
    add_edge!(sched, n1.id, n2.id)
    return sched
end

"""
    construct_partial_project_schedule

Constructs a partial project graph
"""
function construct_partial_project_schedule(spec::ProjectSpec,problem_spec::ProblemSpec,robot_ICs=Vector{BOT_AT}())
    object_ICs = initial_conditions_vector(spec)
    object_FCs = final_conditions_vector(spec)
    operations = get_operations(spec)
    terminal_ops = sort(collect(get_terminal_operation_ids(spec)))

    # Construct Partial Project Schedule
    sched = OperatingSchedule()
    for op in operations
        add_node!(sched, make_node(sched,problem_spec,op))
    end
    for pred in object_ICs
        add_node!(sched, make_node(sched,problem_spec,pred))
    end
    for pred in robot_ICs
        add_new_robot_to_schedule!(sched,pred,problem_spec)
    end
    # add root nodes
    for op_id in terminal_ops
        v = get_vtx(sched, op_id)
        push!(sched.terminal_vtxs, v)
        sched.weights[v] = 1.0
    end
    # Fill in gaps in project schedule (except for GO assignments)
    for op in operations
        op_id = get_operation_id(op)
        for object_id in get_input_ids(op)
            # add action sequence
            object_ic           = get_node_from_id(sched, object_id)
            if length(get_location_ids(object_ic)) > 1 # COLLABORATIVE TRANSPORT
                add_headless_team_delivery_task!(sched,problem_spec,
                    ObjectID(object_id),op_id)
            else # SINGLE AGENT TRANSPORT
                add_headless_delivery_task!(sched,problem_spec,
                    ObjectID(object_id),op_id)
            end
        end
        for object_id in get_output_ids(op)
            add_edge!(sched, op_id, object_id)
        end
    end
    # NOTE: A hack to get speed up on SparseAdjacencyMILP. Not sure if it will work
    for (id, pred) in get_object_ICs(sched)
        v = get_vtx(sched, ObjectID(id))
        if indegree(get_graph(sched),v) == 0
            op_id = get_unique_operation_id()
            op = Operation(id=op_id)
            set_postcondition!(op,pred)
            add_node!(sched,make_node(sched,problem_spec,op))
            add_edge!(sched,op_id,ObjectID(id))
        end
    end
    set_leaf_operation_vtxs!(sched)
    process_schedule!(sched)
    sched
end
function construct_partial_project_schedule(
        robot_ICs::Vector{R},
        prob_spec::ProblemSpec=ProblemSpec(),
        ) where {R<:BOT_AT}
    construct_partial_project_schedule(ProjectSpec(),prob_spec,robot_ICs)
end

"""
    process_schedule(schedule::P) where {P<:OperatingSchedule}

Compute the optimistic start and end times, along with the slack associated
with each vertex in the `schedule`. Slack for each vertex is represented as
a vector in order to handle multi-headed projects.

Args:
* schedule::OperatingSchedule
* [OPTIONAL] t0::Vector{Int}: default = zeros(Int,nv(schedule))
* [OPTIONAL] tF::Vector{Int}: default = zeros(Int,nv(schedule))
"""
function process_schedule(sched::P,t0=zeros(Int,nv(sched)),
        tF=zeros(Int,nv(sched))
    ) where {P<:OperatingSchedule}

    G = get_graph(sched)
    traversal = topological_sort_by_dfs(G)
    n_roots = max(length(sched.terminal_vtxs),1)
    slack = map(i->Inf*ones(n_roots), vertices(G))
    local_slack = map(i->Inf*ones(n_roots), vertices(G))
    # max_deadlines = map(i->typemax(Int), vertices(G))
    # True terminal nodes
    for (i,v) in enumerate(sched.terminal_vtxs)
        slack[v] = Inf*ones(n_roots)
        slack[v][i] = 0 # only slack for corresponding head is set to 0
    end
    ########## Compute Lower Bounds Via Forward Dynamic Programming pass
    for v in traversal
        path_spec = get_path_spec(sched,v)
        Δt_min = path_spec.min_duration
        for v2 in inneighbors(G,v)
            t0[v] = max(t0[v], tF[v2])
        end
        tF[v] = max(tF[v], t0[v] + Δt_min)
    end
    ########### Compute Slack Via Backward Dynamic Programming pass
    for v in reverse(traversal)
        for v2 in outneighbors(G,v)
            local_slack[v] = min.(local_slack[v], (t0[v2] - tF[v]) )
            slack[v] = min.(slack[v], slack[v2] .+ (t0[v2] - tF[v]) )
        end
    end
    t0,tF,slack,local_slack
end
function update_schedule_times!(sched::OperatingSchedule)
    G = get_graph(sched)
    for v in topological_sort_by_dfs(G)
        t0 = get_t0(sched,v)
        for v2 in inneighbors(G,v)
            t0 = max(t0,get_tF(sched,v2))
        end
        set_t0!(sched,v,t0)
        set_tF!(sched,v,max(get_tF(sched,v),t0+get_min_duration(sched,v)))
    end
    return sched
end
function update_slack!(sched::OperatingSchedule)
    G = get_graph(sched)
    n_roots = max(length(sched.terminal_vtxs),1)
    slack = map(i->Inf*ones(n_roots), vertices(G))
    local_slack = map(i->Inf*ones(n_roots), vertices(G))
    for (i,v) in enumerate(sched.terminal_vtxs)
        slack[v][i] = 0 # only slack for corresponding head is set to 0
    end
    for v in reverse(topological_sort_by_dfs(G))
        for v2 in outneighbors(G,v)
            local_slack[v] = min.(local_slack[v], (get_t0(sched,v2) - get_tF(sched,v)))
            slack[v] = min.(slack[v], slack[v2] .+ (get_t0(sched,v2) - get_tF(sched,v)))
        end
    end
    for v in vertices(sched)
        set_slack!(sched,v,slack[v])
        set_local_slack!(sched,v,local_slack[v])
    end
    return sched
end
function process_schedule!(sched)
    update_schedule_times!(sched)
    update_slack!(sched)
end
function reset_schedule_times!(sched::OperatingSchedule,ref)
    for v in vertices(sched)
        set_t0!(sched,v,get_t0(ref,get_vtx_id(sched,v)))
        set_tF!(sched,v,get_tF(ref,get_vtx_id(sched,v)))
    end
    process_schedule!(sched)
end
makespan(sched::OperatingSchedule) = maximum(get_tF(sched))

export
    get_collect_node,
    get_deposit_node

function forward_search(sched,v,f)
    id = v
    node = get_node(sched,v)
    while !f(node)
        id = get_vtx_id(sched, outneighbors(sched,id)[1])
        node = get_node(sched, id)
    end
    return id, node
end
get_collect_node(sched::OperatingSchedule,id::ObjectID) = forward_search(sched,id,n->matches_node_type(n,BOT_COLLECT))
get_deposit_node(sched::OperatingSchedule,id::ObjectID) = forward_search(sched,id,n->matches_node_type(n,BOT_DEPOSIT))

export get_assignment_dict

"""
    get_assignment_dict(assignment_matrix,N,M)

Returns dictionary that maps each robot id to a sequence of tasks.
"""
function get_assignment_dict(assignment_matrix,N,M)
    assignment_dict = Dict{Int,Vector{Int}}()
    for i in 1:N
        vec = get!(assignment_dict,i,Int[])
        k = i
        j = 1
        while j <= M
            if assignment_matrix[k,j] == 1
                push!(vec,j)
                k = j + N
                j = 0
            end
            j += 1
        end
    end
    assignment_dict
end

export apply_assignment_dict

function apply_assignment_dict!(sched::OperatingSchedule,assignment_dict,prob_spec)
    tips = robot_tip_map(sched)
    collect_nodes = Dict(get_object_id(n.node)=>n for n in get_nodes(sched) if matches_template(BOT_COLLECT,n))
    for (robot_idx,itinerary) in assignment_dict
        r = RobotID(robot_idx)
        for object_idx in itinerary
            tip = get_node(sched,tips[r])
            o = ObjectID(object_idx)
            c = collect_nodes[o]
            add_edge!(sched,tip,c)
            propagate_valid_ids!(sched,prob_spec)
            # update tips
            tips = robot_tip_map(sched)
        end
    end
    process_schedule!(sched)
    sched
end


text:
@testset "457.circular-array-loop.jl" begin
    @test circular_array_loop([2,-1,1,2,2]) == true
    @test circular_array_loop([-1,2]) == false
    @test circular_array_loop([-2,1,-1,-2,-2]) == false
end


text:
using KiteConnect
using Test

@test_throws ArgumentError KiteConnect.ltp("INFY")


text:
module REPL

using Revise, SHA, Dates, Pkg
using Genie, Genie.Loggers, Genie.Configuration, Genie.Generator, Genie.Tester, Genie.Util, Genie.FileTemplates

const JULIA_PATH = joinpath(Sys.BINDIR, "julia")

"""
    secret_token() :: String

Generates a random secret token to be used for configuring the SECRET_TOKEN const.
"""
function secret_token() :: String
  sha256("$(randn()) $(Dates.now())") |> bytes2hex
end

"""
"""
function copy_fullstack_app(app_path::String = ".") :: Nothing
  cp(joinpath(@__DIR__, "..", "files", "new_app"), app_path)

  nothing
end

"""
"""
function copy_microstack_app(app_path::String = ".") :: Nothing
  mkdir(app_path)

  for f in ["bin", "config", "public", "src",
            ".gitattributes", ".gitignore",
            "bootstrap.jl", "env.jl", "genie.jl", "routes.jl"]
    cp(joinpath(@__DIR__, "..", "files", "new_app", f), joinpath(app_path, f))
  end

  remove_fingerprint_initializer(app_path)

  nothing
end

"""
"""
function copy_db_support(app_path::String = ".") :: Nothing
  cp(joinpath(@__DIR__, "..", "files", "new_app", "db"), joinpath(app_path, "db"))
  cp(joinpath(@__DIR__, "..", "files", "new_app", "config", "initializers", "searchlight.jl"), joinpath(app_path, "config", "initializers", "searchlight.jl"))

  nothing
end

"""
"""
function copy_mvc_support(app_path::String = ".") :: Nothing
  cp(joinpath(@__DIR__, "..", "files", "new_app", "app"), joinpath(app_path, "app"))

  nothing
end

"""
Generates a valid secrets.jl file with a random SECRET_TOKEN.
"""
function write_secrets_file(app_path::String = ".") :: Nothing
  open(joinpath(app_path, Genie.CONFIG_PATH, "secrets.jl"), "w") do f
    write(f, """const SECRET_TOKEN = "$(secret_token())" """)
  end

  nothing
end

"""
"""
function write_app_custom_files(path::String, app_path::String) :: Nothing
  moduleinfo = FileTemplates.appmodule(path)

  open(joinpath(app_path, "src", moduleinfo[1] * ".jl"), "w") do f
    write(f, moduleinfo[2])
  end

  chmod(joinpath(app_path, "bootstrap.jl"), 0o644)

  open(joinpath(app_path, "bootstrap.jl"), "w") do f
    write(f,
    """
      cd(@__DIR__)
      using Pkg
      pkg"activate ."

      function main()
        include(joinpath("src", "$(moduleinfo[1]).jl"))
      end; main()
    """)
  end

  nothing
end

"""
"""
function install_app_dependencies(app_path::String = ".") :: Nothing
  log("Installing app dependencies")
  pkg"activate ."

  pkg"add Genie#master"
  pkg"add JSON"
  pkg"add Millboard"
  pkg"add Revise"

  nothing
end

"""
"""
function autostart_app(autostart::Bool = true) :: Nothing
  if autostart
    log("Starting your brand new Genie app - hang tight!", :info)
    loadapp(".", autostart = autostart)
  else
    log("Your new Genie app is ready!
        Run \njulia> Genie.loadapp() \nto load the app's environment
        and then \njulia> Genie.startup() \nto start the web server on port 8000.")
  end

  nothing
end

"""
"""
function remove_fingerprint_initializer(app_path::String = ".") :: Nothing
  rm(joinpath(app_path, "config", "initializers", "fingerprint.jl"), force = true)

  nothing
end

"""
"""
function remove_searchlight_initializer(app_path::String = ".") :: Nothing
  rm(joinpath(app_path, "config", "initializers", "searchlight.jl"), force = true)

  nothing
end

"""
    newapp(path::String = "."; autostart::Bool = true, fullstack::Bool = false, dbsupport::Bool = false, mvcsupport::Bool = false) :: Nothing

Scaffolds a new Genie app, setting up the file structure indicated by the various arguments.

# Arguments
- `path::String`: the name of the app and the path where to bootstrap it
- `autostart::Bool`: automatically start the app once the file structure is created
- `fullstack::Bool`: the type of app to be bootstrapped. The fullstack app includes MVC structure, DB connection code, and asset pipeline files.
- `dbsupport::Bool`: bootstrap the files needed for DB connection setup via the SearchLight ORM
- `mvcsupport::Bool`: adds the files used for Flax view templates rendering and working with resources

# Examples
```julia-repl
julia> Genie.newapp("MyGenieApp")
2019-08-06 16:54:15:INFO:Main: Done! New app created at MyGenieApp
2019-08-06 16:54:15:DEBUG:Main: Changing active directory to MyGenieApp
2019-08-06 16:54:15:DEBUG:Main: Installing app dependencies
 Resolving package versions...
  Updating `~/Dropbox/Projects/GenieTests/MyGenieApp/Project.toml`
  [c43c736e] + Genie v0.10.1
  Updating `~/Dropbox/Projects/GenieTests/MyGenieApp/Manifest.toml`

2019-08-06 16:54:27:INFO:Main: Starting your brand new Genie app - hang tight!
 _____         _
|   __|___ ___|_|___
|  |  | -_|   | | -_|
|_____|___|_|_|_|___|

┌ Info:
│ Starting Genie in >> DEV << mode
└
[ Info: Logging to file at MyGenieApp/log/dev.log
[ Info: Ready!
2019-08-06 16:54:32:DEBUG:Main: Web Server starting at http://127.0.0.1:8000
2019-08-06 16:54:32:DEBUG:Main: Web Server running at http://127.0.0.1:8000
```
"""
function newapp(path::String = "."; autostart::Bool = true, fullstack::Bool = false, dbsupport::Bool = false, mvcsupport::Bool = false) :: Nothing
  app_path = abspath(path)

  fullstack ? copy_fullstack_app(app_path) : copy_microstack_app(app_path)

  dbsupport ? (fullstack || copy_db_support(app_path)) : remove_searchlight_initializer(app_path)

  mvcsupport && (fullstack || copy_mvc_support(app_path))

  write_secrets_file(app_path)

  write_app_custom_files(path, app_path)

  Sys.iswindows() ? setup_windows_bin_files(app_path) : setup_nix_bin_files(app_path)

  log("Done! New app created at $(abspath(path))", :info)

  log("Changing active directory to $app_path")
  cd(path)

  install_app_dependencies(app_path)

  autostart_app(autostart)

  nothing
end

"""
    loadapp(path::String = "."; autostart::Bool = false) :: Nothing

Loads an existing Genie app from the file system, within the current Julia REPL session.

# Arguments
- `path::String`: the path to the Genie app on the file system.
- `autostart::Bool`: automatically start the app upon loading it.

# Examples
```julia-repl
shell> tree -L 1
.
├── Manifest.toml
├── Project.toml
├── bin
├── bootstrap.jl
├── config
├── env.jl
├── genie.jl
├── log
├── public
├── routes.jl
└── src

5 directories, 6 files

julia> using Genie

julia> Genie.loadapp(".")
 _____         _
|   __|___ ___|_|___
|  |  | -_|   | | -_|
|_____|___|_|_|_|___|

┌ Info:
│ Starting Genie in >> DEV << mode
└
[ Info: Logging to file at MyGenieApp/log/dev.log
```
"""
function loadapp(path::String = "."; autostart::Bool = false) :: Nothing
  Core.eval(Main, Meta.parse("using Revise"))
  Core.eval(Main, Meta.parse("""include(joinpath("$path", "bootstrap.jl"))"""))
  Core.eval(Main, Meta.parse("Revise.revise()"))
  Core.eval(Main, Meta.parse("using Genie"))

  Core.eval(Main.UserApp, Meta.parse("Revise.revise()"))
  Core.eval(Main.UserApp, Meta.parse("$autostart && Genie.startup()"))

  nothing
end

"""
    setup_windows_bin_files(path::String = ".") :: Nothing

Creates the bin/server and bin/repl binaries for Windows
"""
function setup_windows_bin_files(path::String = ".") :: Nothing
  open(joinpath(path, "bin", "repl.bat"), "w") do f
    write(f, "$JULIA_PATH --color=yes --depwarn=no -q -i -- ../bootstrap.jl %*")
  end

  open(joinpath(path, "bin", "server.bat"), "w") do f
    write(f, "$JULIA_PATH --color=yes --depwarn=no -q -i -- ../bootstrap.jl s %*")
  end

  nothing
end

"""
    setup_nix_bin_files(app_path::String = ".") :: Nothing

Creates the bin/server and bin/repl binaries for *nix systems
"""
function setup_nix_bin_files(app_path::String = ".") :: Nothing
  chmod(joinpath(app_path, "bin", "server"), 0o700)
  chmod(joinpath(app_path, "bin", "repl"), 0o700)

  nothing
end

end


text:
# Set to get interior entry rates
function values_for_test(rng, J, nc, nl; typicalValue = 1.0)
    # vWork_jV = (0.8 * nl) .+ 1.0 .* rand(rng, Float64, J);
    vCollege_jcM = typicalValue .* (1.0 .+ 0.1 .* rand(rng, Float64, J, nc));
    vWork_jV = vec(sum(vCollege_jcM, dims = 2)) ./ nc .+ (0.6 * nl);
    return vWork_jV, vCollege_jcM
end

function make_test_symbol_table()
    st = SymbolTable();
    add_symbol!(st, SymbolInfo(:rankWt, "\\alpha", "Ranking weight", "Admissions"));
    add_symbol!(st, SymbolInfo(:collPrefScale, "\\pi_{c}", "College choice pref scale", "Admissions"));
    add_symbol!(st, SymbolInfo(:pScaleEntry, "\\pi", "Entry pref scale", "Admissions"));
    add_symbol!(st, SymbolInfo(:uLocal, "U_{local}", "Value local", "Admissions"));
    return st
end

# ------------


text:
module GnssWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const BeidouAlmanac = @jimport org.orekit.gnss.BeidouAlmanac
const BeidouOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.BeidouOrbitalElements
const Class = @jimport java.lang.Class
const CombinationType = @jimport org.orekit.estimation.measurements.gnss.CombinationType
const CombinedObservationData = @jimport org.orekit.gnss.CombinedObservationData
const CombinedObservationDataSet = @jimport org.orekit.gnss.CombinedObservationDataSet
const DOP = @jimport org.orekit.gnss.DOP
const DOPComputer = @jimport org.orekit.gnss.DOPComputer
const DataProvidersManager = @jimport org.orekit.data.DataProvidersManager
const ElevationMask = @jimport org.orekit.utils.ElevationMask
const Enum = @jimport java.lang.Enum
const Frequency = @jimport org.orekit.gnss.Frequency
const GLONASSAlmanac = @jimport org.orekit.gnss.GLONASSAlmanac
const GLONASSEphemeris = @jimport org.orekit.gnss.GLONASSEphemeris
const GLONASSOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.GLONASSOrbitalElements
const GNSSOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.GNSSOrbitalElements
const GPSAlmanac = @jimport org.orekit.gnss.GPSAlmanac
const GPSOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.GPSOrbitalElements
const GalileoAlmanac = @jimport org.orekit.gnss.GalileoAlmanac
const GalileoOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.GalileoOrbitalElements
const GeodeticPoint = @jimport org.orekit.bodies.GeodeticPoint
const HatanakaCompressFilter = @jimport org.orekit.gnss.HatanakaCompressFilter
const IRNSSAlmanac = @jimport org.orekit.gnss.IRNSSAlmanac
const IRNSSOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.IRNSSOrbitalElements
const InputStream = @jimport java.io.InputStream
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const MeasurementType = @jimport org.orekit.gnss.MeasurementType
const NamedData = @jimport org.orekit.data.NamedData
const Object = @jimport java.lang.Object
const ObservationData = @jimport org.orekit.gnss.ObservationData
const ObservationDataSet = @jimport org.orekit.gnss.ObservationDataSet
const ObservationType = @jimport org.orekit.gnss.ObservationType
const OneAxisEllipsoid = @jimport org.orekit.bodies.OneAxisEllipsoid
const Optional = @jimport java.util.Optional
const QZSSAlmanac = @jimport org.orekit.gnss.QZSSAlmanac
const QZSSOrbitalElements = @jimport org.orekit.propagation.analytical.gnss.QZSSOrbitalElements
const RinexHeader = @jimport org.orekit.gnss.RinexHeader
const RinexLoader = @jimport org.orekit.gnss.RinexLoader
const SEMParser = @jimport org.orekit.gnss.SEMParser
const SatelliteSystem = @jimport org.orekit.gnss.SatelliteSystem
const SignalCode = @jimport org.orekit.gnss.SignalCode
const TimeScale = @jimport org.orekit.time.TimeScale
const TimeScales = @jimport org.orekit.time.TimeScales
const Vector2D = @jimport org.hipparchus.geometry.euclidean.twod.Vector2D
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D
const YUMAParser = @jimport org.orekit.gnss.YUMAParser

export BeidouAlmanac
export CombinedObservationData
export CombinedObservationDataSet
export DOP
export DOPComputer
export Frequency
export GLONASSAlmanac
export GLONASSEphemeris
export GPSAlmanac
export GalileoAlmanac
export HatanakaCompressFilter
export IRNSSAlmanac
export MeasurementType
export ObservationData
export ObservationDataSet
export ObservationType
export QZSSAlmanac
export RinexHeader
export RinexLoader
export SEMParser
export SatelliteSystem
export SignalCode
export YUMAParser
export compare_to
export compute
export create
export describe_constable
export equals
export filter
export get_af0
export get_af1
export get_af2
export get_agency_name
export get_almanacs
export get_antenna_azimuth
export get_antenna_b_sight
export get_antenna_height
export get_antenna_number
export get_antenna_phase_center
export get_antenna_reference_point
export get_antenna_type
export get_antenna_zero_direction
export get_aodc
export get_aode
export get_approx_pos
export get_bgd
export get_bgde1e5a
export get_bgde5b_e1
export get_center_mass
export get_cic
export get_cis
export get_class
export get_clk_offset
export get_combination_type
export get_combined_m_hz_frequency
export get_crc
export get_crs
export get_cuc
export get_cus
export get_date
export get_declaring_class
export get_delta_i
export get_delta_t
export get_delta_t_dot
export get_e
export get_eccentricities
export get_elevation_mask
export get_frequency
export get_frequency_channel
export get_gamma_n
export get_gdop
export get_glo2utc
export get_glo_offset
export get_gnss_nb
export get_gps2_glo
export get_hdop
export get_header
export get_health
export get_health_e1
export get_health_e5a
export get_health_e5b
export get_i0
export get_i_dot
export get_interval
export get_iod
export get_iod_nav
export get_iodc
export get_iode
export get_iodec
export get_key
export get_lambda
export get_leap_seconds
export get_leap_seconds_day_num
export get_leap_seconds_future
export get_leap_seconds_week_num
export get_list_applied_dcbs
export get_list_applied_pcvs
export get_location
export get_loss_of_lock_indicator
export get_m0
export get_m_hz_frequency
export get_marker_name
export get_marker_number
export get_marker_type
export get_mean_motion
export get_measurement_type
export get_min_elevation
export get_n4
export get_na
export get_name
export get_observation_code
export get_observation_data
export get_observation_data_sets
export get_observation_type
export get_observer_name
export get_omega0
export get_omega_dot
export get_pa
export get_pdop
export get_phase_shift_corrections
export get_prn
export get_prn_number
export get_prn_numbers
export get_ratio
export get_rcvr_clk_offset
export get_receiver_number
export get_receiver_type
export get_receiver_version
export get_rinex_version
export get_sat_configuration
export get_satellite_system
export get_signal_code
export get_signal_strength
export get_signal_strength_unit
export get_sma
export get_source
export get_supported_names
export get_svn
export get_t_first_obs
export get_t_last_obs
export get_tdop
export get_tgd
export get_tgd1
export get_tgd2
export get_time
export get_tn
export get_toc
export get_ura
export get_used_observation_data
export get_value
export get_vdop
export get_wavelength
export get_week
export get_x
export get_x_dot
export get_x_dot_dot
export get_y
export get_y_dot
export get_y_dot_dot
export get_z
export get_z_dot
export get_z_dot_dot
export hash_code
export load_data
export name
export notify
export notify_all
export ordinal
export parse_satellite_system
export still_accepts_data
export to_string
export value_of
export values
export wait
export with_elevation_mask
export with_min_elevation

include("AntennaWrapper/AntennaWrapper.jl")
include("AttitudeWrapper/AttitudeWrapper.jl")
include("beidou_almanac.jl")
include("combined_observation_data.jl")
include("combined_observation_data_set.jl")
include("dop")
include("dop_computer.jl")
include("frequency.jl")
include("galileo_almanac.jl")
include("glonass_almanac.jl")
include("glonass_ephemeris.jl")
include("gps_almanac.jl")
include("hatanaka_compress_filter.jl")
include("irnss_almanac.jl")
include("measurement_type.jl")
include("observation_data.jl")
include("observation_data_set.jl")
include("observation_type.jl")
include("qzss_almanac.jl")
include("rinex_header.jl")
include("rinex_loader.jl")
include("satellite_system.jl")
include("sem_parser.jl")
include("signal_code.jl")
include("yuma_parser.jl")

end


text:
#========================================================================================#
"""
	NBodies

A system of N interacting bodies moving in dim dimensions over time T in res timesteps.

Author: Niall Palfreyman, 24/3/2022.
"""
module NBodies

# Externally callable methods of NBodies
export NBody, add!, simulate

#-----------------------------------------------------------------------------------------
# Module types:

"""
	NBody

An NBody system capable of containing multiple (N) bodies that gravitationally interact
with each other.
"""
mutable struct NBody
	N							# Number of bodies
	T							# Duration of simulation
	res							# Timestep resolution
	xi							# Initial positions of bodies
	pi							# Initial momenta of bodies
	m							# Masses of bodies

	"Construct a new NBody"
	function NBody( T=40, resolution=20000)
		# Initialise all fields of the decoding apparatus:
		new(
			0,					# Initially no bodies in the system
			T,					# Duration
			resolution,			# Timestep resolution
			[],					# Initial positions
			[],					# Initial momenta
			[],					# Masses of bodies
		)
	end
end

#-----------------------------------------------------------------------------------------
# Module methods:

"""
	addbody!( nbody::NBody, xi, pi, m=1)

Add to the system a new body with initial position and momentum xi, pi, and with mass m.
"""
function addbody!( nbody::NBody, xi, pi, m=1)
	push!( nbody.xi, xi)
	push!( nbody.pi, pi)
	push!( nbody.m, m)
	nbody.N  += 1
end

"""
	relpos( locations)

A useful helper function that calculates relative positions of a vector x of locations.
"""
function relpos( locations)
	locPerBody = repeat(locations,1,length(locations))
	locPerBody - permutedims(locPerBody)
end

"""
	simulate( nbody::NBody)

Run a simulation of the given NBody system over duration divided into res timesteps.
"""
function simulate( nbody)
	dt = nbody.T / nbody.res
	t = 0:dt:nbody.T
	x = Vector{typeof(nbody.xi)}(undef,nbody.res+1)
	p = Vector{typeof(nbody.pi)}(undef,nbody.res+1)

	# Pre-calculate products of masses:
	G = 1
	Gmm = G * nbody.m * nbody.m'

	# Initialisation:
	x[1] = nbody.xi
	p[1] = nbody.pi

	# Simulation:
	for n = 1:nbody.res
		# Simulate timestep n:
		x[n+1] = x[n] + dt*p[n]./nbody.m

		diff = relpos(x[n])
		invSq = abs.(diff'.*diff) .^ (3/2)
		for i in 1:length(x[n]) invSq[i,i] += 1 end
		forceFROMEachSource = -Gmm .* diff ./ invSq
		forceONEachSource = sum( forceFROMEachSource, dims=2)
		p[n+1] = p[n] + dt * forceONEachSource[:]
	end

	(t,x,p)
end

#-----------------------------------------------------------------------------------------
using GLMakie

"""
	demo()

Demonstrate simulation of a simple 2-body problem.
"""
function demo()
	# Initialise the system:
	nb = NBody( 2.25, 10000)
	addbody!( nb, -1, 0)
	addbody!( nb,  1, 0)

	# Run the simulation:
	(t,x,p) = simulate(nb)

	# Display some fancy graphics:
	fig = Figure(resolution=(600, 600))
	ax = Axis(fig[1, 1], xlabel = "Time (t)", title = "N-body Trajectories")

	lines!(ax, t, map(el->el[1],x), color=:blue)		# Body 1
	lines!(ax, t, map(el->el[2],x), color=:red)			# Body 2
	display(fig)
end

end		# of NBodies


text:
#=
Clustering methods. Each method needs
1) Type with fields necessary for algorithm
2) function "cluster" to operate on sort with type field defined above
3) any other necessary functions for clustering algorithm

=#

export ClusterOSort, ClusterNone, ClusterWindow, ClusterTemplate

function clusterprepare(c::Cluster, sort::Sorting,p)
end

#=
OSort

Rutishauser 2006
=#

mutable struct ClusterOSort <: Cluster
    clusters::Array{Float64,2}
    clusterWeight::Array{Int64,1}
    numClusters::Int64
    Tsm::Float64
    m_k::Float64
    k::Int64
    m_l::Float64
    s_k::Float64
    s_l::Float64
end

ClusterOSort()=ClusterOSort(50)

ClusterOSort(n::Int64)=ClusterOSort(10,n)

ClusterOSort(m::Int64,n::Int64)=ClusterOSort(zeros(Float64,n,m),zeros(Int64,m),0,1.0,0.0,1,0.0,0.0,0.0)

ClusterOSort(n::Int64,c::ClusterOSort)=ClusterOSort(zeros(Float64,n,size(c.clusters,2)),c.clusterWeight,c.numClusters,c.Tsm,c.m_k,c.k,c.m_l,c.s_k,c.s_l)

function cluster(c::ClusterOSort,sort::Sorting)

    x=getdist(sort)

    id=0

    #add new cluster or assign to old
    if x==0
        sort.c.numClusters+=1
        sort.c.clusters[:,sort.c.numClusters]=sort.features[:]
        sort.c.clusterWeight[sort.c.numClusters]=1
        id=sort.c.numClusters
    else
        if sort.c.clusterWeight[x]<20
            sort.c.clusterWeight[x]+=1
            for i=1:size(sort.c.clusters,1)
                sort.c.clusters[i,x]=(sort.c.clusterWeight[x]-1)/sort.c.clusterWeight[x]*sort.c.clusters[i,x]+1/sort.c.clusterWeight[x] .* sort.features[i]
            end
        else
            for i=1:size(sort.c.clusters,1)
                sort.c.clusters[i,x]=.95 .* sort.c.clusters[i,x]+.05 .* sort.features[i]
            end
        end
        id=x
    end

    if x>0
        findmerge!(sort,id)
    end

    id
end

function getdist(sort::Sorting)

    dist=zeros(Float64,sort.c.numClusters)
    for i=1:sort.c.numClusters
        dist[i]=norm(sort.features-sort.c.clusters[:,i])^2
    end
    if sort.c.numClusters<1
        return 0
    else
        ind=indmin(dist)
        if dist[ind]<sort.c.Tsm
            return ind
        else
            return 0
        end
    end
end

function findmerge!(sort::Sorting,id::Int64)

    if sort.c.numClusters>1

        dist=zeros(Float64,sort.c.numClusters)

        for i=1:sort.c.numClusters
            if i==id
                dist[i]=1.0e20
            else
                dist[i]=norm(sort.c.clusters[:,id]-sort.c.clusters[:,i])^2
            end
        end

        ind=indmin(dist)
        if dist[ind]<sort.c.Tsm #if less than threshold, combine clusters

            newid=min(ind,id)

            sort.c.clusters[:,newid]=(sort.c.clusters[:,ind]+sort.c.clusters[:,id])/2
            sort.c.clusterWeight[newid]=1

            delid=max(ind,id)

            for i=delid:sort.c.numClusters
                sort.c.clusters[:,i]=sort.c.clusters[:,i+1]
                sort.c.clusterWeight[i]=sort.c.clusterWeight[i+1]
            end

            sort.c.numClusters-=1

            findmerge!(sort,newid) #compare newly created cluster to existing clusters
        end
    end

    nothing
end

function clusterprepare(c::ClusterOSort,sort::Sorting,p)

    sort.c.k+=1

    sort.c.m_k=sort.c.m_l+(p-sort.c.m_l)/sort.c.k
    sort.c.s_k=sort.c.s_l+(p-sort.c.m_l)*(p-sort.c.m_k)

    myvar=sort.c.s_k/(sort.c.k-1)

    sort.c.Tsm=(myvar)*100

    sort.c.m_l=sort.c.m_k
    sort.c.s_l=sort.c.s_k

    nothing
end

#=
No clustering
=#

mutable struct ClusterNone <: Cluster
end

ClusterNone(n::Int64)=ClusterNone()

ClusterNone(n::Int64, c::ClusterNone)=ClusterNone()

cluster(c::ClusterNone,sort::Sorting)=1

function clusterprepare(c::ClusterNone, sort::Sorting,p)
end

#=
Window Discrimination
=#

struct mywin
    x1::Int64
    x2::Int64
    y1::Float64
    y2::Float64
end

mutable struct ClusterWindow <: Cluster
    win::Array{Array{mywin,1},1}
    hits::Array{UInt8,1}
end

ClusterWindow()=ClusterWindow(Array(Array{mywin,1},0),zeros(UInt8,0))

ClusterWindow(n::Int64)=ClusterWindow(Array(Array{mywin,1},0),zeros(UInt8,0))

ClusterWindow(n::Int64,c::ClusterWindow)=ClusterWindow(Array(Array{mywin,1},0),zeros(UInt8,0))

function cluster(c::ClusterWindow,sort::Sorting)

    @inbounds for i=1:length(c.win) #Loop over all clusters
        for j=1:length(c.win[i]) #Loop over all windows for cluster
            a1=c.win[i][j].x1
            a2=c.win[i][j].x2
            b1=c.win[i][j].y1
            b2=c.win[i][j].y2
            for k=(c.win[i][j].x1-1):(c.win[i][j].x2+1)
                if intersect(a1,a2,k,k+1,b1,b2,sort.features[k],sort.features[k+1])
                    c.hits[i]+=1
                    break
                end
            end
        end
    end

    if length(c.hits)==0
        return 1
    else
        myind=indmax(c.hits)
        if c.hits[myind]==0
            for i=1:length(c.hits)
                c.hits[i]=0
            end
            return 1
        else
            for i=1:length(c.hits)
                c.hits[i]=0
            end
            return myind+1
        end
    end
end

ccw(x1,x2,x3,y1,y2,y3)=(y3-y1) * (x2-x1) > (y2-y1) * (x3-x1)

function intersect(x1,x2,x3,x4,y1,y2,y3,y4)
    (ccw(x1,x3,x4,y1,y3,y4) != ccw(x2,x3,x4,y2,y3,y4))&&(ccw(x1,x2,x3,y1,y2,y3) != ccw(x1,x2,x4,y1,y2,y4))
end

function clusterprepare(c::ClusterWindow, sort::Sorting,p)
end

#=
Template matching
=#

mutable struct ClusterTemplate <: Cluster
    templates::Array{Float64,2}
    sig_min::Array{Float64,2}
    sig_max::Array{Float64,2}
    misses::Int64
    num::Int64
    tol::Array{Float64,1}
end

ClusterTemplate()=ClusterTemplate(48)

ClusterTemplate(n::Int64)=ClusterTemplate(zeros(Float64,n,10),zeros(Float64,n,10),zeros(Float64,n,10),3,0,ones(Float64,10))

ClusterTemplate(n::Int64,c::ClusterTemplate)=ClusterTemplate(n)

function cluster(c::ClusterTemplate,sort::Sorting)

    id=0

    for i=1:c.num
        mymisses=0
        for j=1:size(c.templates,1)
            if (sort.features[j]<(c.templates[j,i]-(c.sig_min[j,i]*c.tol[i])))|(sort.features[j]>(c.templates[j,i]+(c.sig_max[j,i]*c.tol[i])))
                mymisses+=1
                if mymisses>c.misses
                    break
                end
            end
        end
        if mymisses<c.misses
            id=i
            break
        end
    end

    id+1
end

function clusterprepare(c::ClusterTemplate,sort::Sorting,p)
end

#=
CLASSIT
=#

#=

type attribute
    S::Float64
    SS::Float64
    std::Float64
end

type node
    n::Int64 #count
    a::Array{attribute,1} #attributes
    instdsum::Float64
    c::Array{node,1} #children
    proc::Bool
end

function node(x::Array{Float64,1})

    a=[attribute(x[i],x[i]^2,acuity) for i=1:length(x)]
    instdsum=0.0
    for i=1:length(x)
        instdsum+=1/a[i].std
    end

    node(1,a,instdsum,Array(node,0),false)
end

function node(N1::node,N2::node)
    a=Array(attribute,length(N1.a))
    instdsum=0.0
    prior=N1.n+N2.n
    a=[attribute(N1.a[i].S+N2.a[i].S,N1.a[i].SS+N2.a[i].SS,0.0) for i=1:length(N1.a)]
    for j=1:length(N1.a)
        a[j].std=sqrt((1/(prior))*((a[j].SS+x[j]^2)-((a[j].S+x[j])^2/(prior+1))))
        if a[j].std<acuity
            a[j].std=acuity
        end
        instdsum+=1/a[j].std
    end
    node(prior,a,instdsum,[N1;N2],false)
end

type ClusterClassit <: Cluster
    clustree::node
end

function ClusterClassit()
    ClusterClassit(node(zeros(Float64,3)))
end

function ClusterClassit(N::Int64)
    ClusterClassit(node(zeros(Float64,N)))
end

function cluster(c::ClusterClassit,sort::Sorting)
    x=sort.features[:]
    cobweb(sort.c.clustree,x)
end

function updatestd(N::node,x::Array{Float64,1},ind::Int64)

    stdmat=0.0

    for j=1:length(N.a)

        std=sqrt((1/(N.c[ind].n))*((N.c[ind].a[j].SS+x[j]^2)-((N.c[ind].a[j].S+x[j])^2/(N.c[ind].n+1))))
        if std<acuity
            std=acuity
        end

        stdmat+=1/std
    end

    stdmat

end

function parentstd(N::node,x::Array{Float64,1})

    N.instdsum=0.0

    for j=1:length(N.a)
        N.a[j].S=N.a[j].S+x[j]
        N.a[j].SS=N.a[j].SS+x[j]^2
        N.a[j].std=sqrt((1/(N.n-1))*((N.a[j].SS)-((N.a[j].S)^2/(N.n))))
        if N.a[j].std<acuity
            N.a[j].std=acuity
        end

        N.instdsum+=1/N.a[j].std
    end

    nothing
end

function cob_incorp(N::node,ind::Int64,stdmat::Float64)

    CU=0.0

    for i=1:length(N.c)
        if i==ind
            CU += (N.c[i].n+1)/N.n*stdmat
        else
            CU += (N.c[i].n)/N.n*N.c[i].instdsum
        end
    end

    CU=CU/length(N.c)
end

function cob_create(N::node)
    CU=0.0

    for i=1:length(N.c)
        CU += (N.c[i].n/N.n)*N.c[i].instdsum
    end

    CU += 1/N.n*sum([1/acuity for i=1:length(N.a)])

    CU=CU/(length(N.c)+1)
end

function cob_merge(N::node,best1ind::Int64,best2ind::Int64,x::Array{Float64,1})

    CU=0.0
    count=0

    for i=1:length(N.c)
           if (i==best1ind | i==best2ind) & count==0
               stdmat=0.0
               prior=N.c[best1ind].n+N.c[best2ind].n
               for j=1:length(N.a)
                   SS=N.c[best1ind].a[j].SS+N.c[best2ind].a[j].SS
                   S=N.c[best1ind].a[j].S+N.c[best2ind].a[j].S
                   std=sqrt((1/(prior))*((SS+x[j]^2)-((S+x[j])^2/(prior+1))))
                   if std<acuity
                       std=acuity
                   end

                   stdmat+=1/std
               end
               CU += (prior+1)/N.n*stdmat
               count=1
           else
            CU += (N.c[i].n)/N.n*N.c[i].instdsum
           end
    end

    CU=CU/(length(N.c)-1)

end

function cob_split(N::node,best1ind::Int64)
    CU=0.0

    for i=1:length(N.c)
        if i==best1ind
            for j=1:length(N.c[i].c)
                CU += (N.c[i].c[j].n)/N.n*N.c[i].c[j].instdsum
            end
        else
            CU += (N.c[i].n)/N.n*N.c[i].instdsum
        end
    end

    CU=CU/(length(N.c)+length(N.c[best1ind].c))
end

function cobweb(N::node, x::Array{Float64,1})

    if length(N.c)==0 #If leaf

        CU=0.0
        instdsum=0.0
        for i=1:length(N.a)
            std=sqrt((1/(N.n))*((N.a[i].SS+x[i]^2)-((N.a[i].S+x[i])^2/(N.n+1))))
            if std<acuity
                std=acuity
            end
            instdsum+=1/std
        end

        CU=(1/N.n*1/acuity+N.n/(N.n+1)*N.instdsum - instdsum)/2

        if CU<cob_cutoff

        else
            N.proc=false
            newnode1=node(x)
            newnode2=deepcopy(N)
            push!(N.c,newnode1)
            push!(N.c,newnode2)
        end
        N.n+=1
        parentstd(N,x)
    else

        if N.proc==false
            N.n+=1
            parentstd(N,x)
        end

        bestoverallind=0
        bestoverall=0.0

        #add I to each child and get CUs
        best1=0.0
        best1ind=0
        best2=0.0
        best2ind=0
        for i=1:length(N.c)
            stdmat=updatestd(N,x,i)
            S1=cob_incorp(N,i,stdmat)
            if S1>best2
                if S1>best1
                    best1=S1
                    best1ind=i
                else
                    best2=S1
                    best2ind=i
                end
            end
        end
        bestoverallind=1
        bestoverall=best1

        #add I as new singleton child and get CU
        S2=cob_create(N)
        if S2>bestoverall
            bestoverallind=2
            bestoverall=S2
        end

        #merge best and second best CUs from 1 and add I to merged result
        if length(N.c)>2
            S3=cob_merge(N,best1ind,best2ind,x)
            if S3>bestoverall
                bestoverallind=3
                bestoverall=S3
            end
        end

        #promote children of best child to be children of P, and add instance
        if length(N.c[best1ind].c)>0
            S4=cob_split(N,best1ind)
            if S4>bestoverall
                bestoverallind=4
                bestoverall=S4
            end
        end

        if bestoverallind==1
            N.proc=false
            cobweb(N.c[best1ind],x) #go down tree at best index
        elseif bestoverallind==2
            N.proc=false
            newnode=node(x)
            push!(N.c,newnode)
        elseif bestoverallind==3
            newnode=node(N.c[best1ind],N.c[best2ind])
            deleteat!(N.c,best1ind)
            if best1ind>best2ind
                deleteat!(N.c,best2ind)
            else
                deleteat!(N.c,best2ind-1)
            end
            push!(N.c,newnode)
            N.c[end].proc=true
            cobweb(N.c[end],x)
        elseif bestoverallind==4
            append!(N.c,N.c[best1ind].c)
            deleteat!(N.c,best1ind)
            N.proc=true
            cobweb(N,x)
        end

    end

    return true
end

=#


text:
export AbstractLogNumber, PrimitiveLogNumber, WrappedLogNumber,
    LogFloat64, LogFloat32, LogFloat16,
    Log, LogNumber, logvalue,
    floattype, logtype,
    infty,
    LogZero64, LogInf64, LogNaN64,
    LogZero32, LogInf32, LogNaN32,
    LogZero16, LogInf16, LogNaN16,
    LogZero, LogNaN, LogInf

abstract type AbstractLogNumber{F<:AbstractFloat} <: AbstractFloat end

abstract type PrimitiveLogNumber{F<:AbstractFloat} <: AbstractLogNumber{F} end

primitive type LogFloat64 <: PrimitiveLogNumber{Float64} 64 end
primitive type LogFloat32 <: PrimitiveLogNumber{Float32} 32 end
primitive type LogFloat16 <: PrimitiveLogNumber{Float16} 16 end

logvalue(x::PrimitiveLogNumber{F}) where F = reinterpret(F, x)

mutable struct WrappedLogNumber{F<:AbstractFloat} <: AbstractLogNumber{F}
    log::F
    WrappedLogNumber{T}(x::T) where {T<:AbstractFloat} = new{T}(x)
end

logvalue(x::WrappedLogNumber) = x.log

# Pseudo-constructors and conversions
const PrimitiveFloat = Union{Float16, Float32, Float64}

floattype(::Type{<:AbstractLogNumber{F}}) where F = F
floattype(x) = floattype(typeof(x))

logtype(::Type{Float64}) = LogFloat64
logtype(::Type{Float32}) = LogFloat32
logtype(::Type{Float16}) = LogFloat16
logtype(::Type{F}) where {F<:AbstractFloat} = WrappedLogNumber{F}
logtype(::Type{N}) where N = logtype(float(N))
logtype(x) = logtype(typeof(x))

# reinterpret and LogNumber perform conversion without transforming to log space
Base.reinterpret(::Type{WrappedLogNumber{F}}, x::F) where {F} = WrappedLogNumber{F}(x)

LogNumber(x::F) where {F<:AbstractFloat} = reinterpret(logtype(F), x)
LogNumber(x) = LogNumber(float(x))
LogNumber(::Type{F}, x) where {F<:AbstractFloat} = LogNumber(convert(F, x))

# convert and log transform with implicit log transformation
Base.convert(::Type{L}, x::Number) where {L<:AbstractLogNumber} =
    reinterpret(L, log(convert(floattype(L), x)))

Log(x::F) where {F<:AbstractFloat} = convert(logtype(F), x)
Log(x) = Log(float(x))
Log(::Type{F}, x) where {F<:AbstractFloat} = Log(convert(F, x))

# converting back to normal space with exp
Base.convert(::Type{L}, x::AbstractLogNumber) where {L<:AbstractLogNumber} =
    LogNumber(convert(floattype(L), logvalue(x)))
Base.convert(::Type{T}, x::AbstractLogNumber) where {T<:Number} = convert(T, exp(logvalue(x)))

Base.promote_rule(::Type{S}, ::Type{T}) where {S<:AbstractLogNumber, T<:AbstractLogNumber} =
    logtype(promote_type(floattype(S), floattype(T)))
Base.promote_rule(::Type{L}, ::Type{R}) where {L<:AbstractLogNumber, R<:Real} =
    logtype(promote_type(floattype(L), R))

# Constants
const LogZero64 = LogNumber(-Inf64)
const LogInf64 = LogNumber(Inf64)
const LogNaN64 = LogNumber(NaN64)

const LogZero32 = LogNumber(-Inf32)
const LogNaN32 = LogNumber(NaN32)
const LogInf32 = LogNumber(Inf32)

const LogZero16 = LogNumber(-Inf16)
const LogNaN16 = LogNumber(NaN16)
const LogInf16 = LogNumber(Inf16)

const LogZero = LogZero64
const LogNaN = LogNaN64
const LogInf = LogInf64

infty(::Type{T}) where {T} = one(T) / zero(T)
infty(::T) where {T} = infty(T)

Base.zero(::Type{LogFloat64}) = LogZero64
Base.zero(::Type{LogFloat32}) = LogZero32
Base.zero(::Type{LogFloat16}) = LogZero16
Base.zero(::Type{WrappedLogNumber{F}}) where F = -infty(F)
Base.zero(x::AbstractLogNumber) = zero(typeof(x))

Base.one(::Type{LogFloat64}) = LogNumber(0e0)
Base.one(::Type{LogFloat32}) = LogNumber(0f0)
Base.one(::Type{LogFloat16}) = LogNumber(Float16(0f0))
Base.one(::Type{WrappedLogNumber{F}}) where F = zero(F)
Base.one(x::AbstractLogNumber) = one(typeof(x))


text:
using VectorizationBase: Vec, MM, mask, pick_vector_width_shift

using Base.Cartesian: @nexprs, @ntuple
using Base: IEEEFloat

import LinearAlgebra: lmul!, rmul!

@generated function rmul!(A::StructArray{MultiFloat{T,N}, 2}, G::Rotation2{MultiFloat{T,N},MultiFloat{T,N}}) where {T<:IEEEFloat,N}

    W, Wshift = pick_vector_width_shift(T)

    quote
        @nexprs $N k -> p_k = stridedpointer(reinterpret(T, getfield(A, 1)[k]))

        m = size(A, 1)

        i = 1

        for step = Base.OneTo(m >> $Wshift)
            @nexprs $N k -> v_k_1 = vload(p_k, (MM{$W}(i), G.i))
            @nexprs $N k -> v_k_2 = vload(p_k, (MM{$W}(i), G.i+1))

            a₁ = MultiFloat(@ntuple $N k -> v_k_1)
            a₂ = MultiFloat(@ntuple $N k -> v_k_2)

            a₁′ = a₁ * Vec{$W,$T}(G.c) + a₂ * Vec{$W,$T}(G.s)
            a₂′ = a₁ * -Vec{$W,$T}(G.s) + a₂ * Vec{$W,$T}(G.c)

            @nexprs $N k -> vstore!(p_k, a₁′._limbs[k], (MM{$W}(i), G.i))
            @nexprs $N k -> vstore!(p_k, a₂′._limbs[k], (MM{$W}(i), G.i+1))

            i += $W
        end

        remainder = m & ($W - 1)
        remainder == 0 && return A
        remainder_mask = mask($T, remainder)

        # Finally do the masked load, multiply and store.

        @nexprs $N k -> v_k_1 = vload(p_k, (MM{$W}(i), G.i), remainder_mask)
        @nexprs $N k -> v_k_2 = vload(p_k, (MM{$W}(i), G.i+1), remainder_mask)

        a₁ = MultiFloat(@ntuple $N k -> v_k_1)
        a₂ = MultiFloat(@ntuple $N k -> v_k_2)

        a₁′ = a₁ * Vec{$W,$T}(G.c) + a₂ * Vec{$W,$T}(G.s)
        a₂′ = a₁ * -Vec{$W,$T}(G.s) + a₂ * Vec{$W,$T}(G.c)

        @nexprs $N k -> vstore!(p_k, a₁′._limbs[k], (MM{$W}(i), G.i), remainder_mask)
        @nexprs $N k -> vstore!(p_k, a₂′._limbs[k], (MM{$W}(i), G.i+1), remainder_mask)

        A
    end
end

@generated function lmul!(G::Rotation2{MultiFloat{T,N},MultiFloat{T,N}}, A::StructArray{MultiFloat{T,N}, 2}) where {T<:IEEEFloat,N}

    W, Wshift = pick_vector_width_shift(T)

    quote
        @nexprs $N k -> p_k = stridedpointer(reinterpret(T, getfield(A, 1)[k]))

        m = size(A, 2)

        i = 1

        for step = Base.OneTo(m >> $Wshift)
            @nexprs $N k -> v_k_1 = vload(p_k, (G.i, MM{$W}(i)))
            @nexprs $N k -> v_k_2 = vload(p_k, (G.i+1, MM{$W}(i)))

            a₁ = MultiFloat(@ntuple $N k -> v_k_1)
            a₂ = MultiFloat(@ntuple $N k -> v_k_2)

            a₁′ = Vec{$W,$T}(G.c) * a₁ + Vec{$W,$T}(G.s) * a₂
            a₂′ = -Vec{$W,$T}(G.s) * a₁ + Vec{$W,$T}(G.c) * a₂

            @nexprs $N k -> vstore!(p_k, a₁′._limbs[k], (G.i, MM{$W}(i)))
            @nexprs $N k -> vstore!(p_k, a₂′._limbs[k], (G.i+1, MM{$W}(i)))

            i += $W
        end

        remainder = m & ($W - 1)
        remainder == 0 && return A
        remainder_mask = mask($T, remainder)

        # Finally do the masked load, multiply and store.
        @nexprs $N k -> v_k_1 = vload(p_k, (G.i, MM{$W}(i)), remainder_mask)
        @nexprs $N k -> v_k_2 = vload(p_k, (G.i+1, MM{$W}(i)), remainder_mask)

        a₁ = MultiFloat(@ntuple $N k -> v_k_1)
        a₂ = MultiFloat(@ntuple $N k -> v_k_2)

        a₁′ = Vec{$W,$T}(G.c) * a₁ + Vec{$W,$T}(G.s) * a₂
        a₂′ = -Vec{$W,$T}(G.s) * a₁ + Vec{$W,$T}(G.c) * a₂

        @nexprs $N k -> vstore!(p_k, a₁′._limbs[k], (G.i, MM{$W}(i)), remainder_mask)
        @nexprs $N k -> vstore!(p_k, a₂′._limbs[k], (G.i+1, MM{$W}(i)), remainder_mask)

        A
    end
end


text:
@testset "html tokens" begin
    t = """
        {{ A }} <script abc></script><script> <!-- -->
        """ |> FP.default_html_tokenizer
    deleteat!(t, 1)
    @test t[1].name == :DBB_OPEN
    @test t[2].name == :DBB_CLOSE
    @test t[3].name == :SCRIPT_OPEN
    @test t[4].name == :SCRIPT_CLOSE
    @test t[5].name == :SCRIPT_OPEN
    @test t[6].name == :COMMENT_OPEN
    @test t[7].name == :COMMENT_CLOSE
end


text:
"""
A handle for an individual database in the DB environment.
"""
mutable struct DBI
    handle::MDB_dbi
    name::String
end

"Check if database is open"
isopen(dbi::DBI) = dbi.handle != zero(Cuint)

"Open a database in the environment"
function open(txn::Transaction, dbname::String = ""; flags::Cuint = zero(Cuint))
    cdbname = length(dbname) > 0 ? dbname : Ptr{Cchar}(C_NULL)
    handle = Ref{MDB_dbi}()
    mdb_dbi_open(txn.handle, cdbname, flags, handle)
    return DBI(handle[], dbname)
end

"Wrapper of DBI `open` for `do` construct"
function open(f::Function, txn::Transaction, dbname::String = ""; flags::Cuint = zero(Cuint))
    dbi = open(txn, dbname, flags=flags)
    tenv = env(txn)
    try
        f(dbi)
    finally
        close(tenv, dbi)
    end
end

"Close a database handle"
function close(env::Environment, dbi::DBI)
    if !isopen(env)
        warn("Environment is closed")
    end
    _mdb_dbi_close(env.handle, dbi.handle)
    dbi.handle = zero(Cuint)
    return
end

"Retrieve the DB flags for a database handle"
function flags(txn::Transaction, dbi::DBI)
    flags = Ref{Cuint}(0)
    mdb_dbi_flags(txn.handle, dbi.handle, flags)
    return flags[]
end

"""Empty or delete+close a database.

If parameter `delete` is `false` DB will be emptied, otherwise
DB will be deleted from the environment and DB handle will be closed
"""
function drop(txn::Transaction, dbi::DBI; delete = false)
    del = Cint(delete)
    mdb_drop(txn.handle, dbi.handle, del)
end

toref(v) = isbitstype(typeof(v)) ? [v] : v
toref(v::Ptr{Nothing}) = v

"Store items into a database"
function put!(txn::Transaction, dbi::DBI, key, val; flags::Cuint = zero(Cuint))
    mdb_key_ref = Ref(MDBValue(toref(key)))
    mdb_val_ref = Ref(MDBValue(toref(val)))
    r = mdb_put(txn.handle, dbi.handle, mdb_key_ref, mdb_val_ref, flags)
    r
end

"Delete items from a database"
function delete!(txn::Transaction, dbi::DBI, key, val=C_NULL)
    mdb_key_ref = Ref(MDBValue(toref(key)))
    mdb_val_ref = val === C_NULL ? Ref(MDBValue()) : Ref(MDBValue(toref(val)))

    mdb_del(txn.handle, dbi.handle, mdb_key_ref, mdb_val_ref)
end

"Get items from a database"
function get(txn::Transaction, dbi::DBI, key, ::Type{T}) where T
    mdb_key_ref = Ref(MDBValue(toref(key)))
    mdb_val_ref = Ref(MDBValue())

    # Get value
    mdb_get(txn.handle, dbi.handle, mdb_key_ref, mdb_val_ref)

    # Convert to proper type
    return convert(T, mdb_val_ref)
end