text:
using ConventionalApp

proj = Project(@__DIR__)
setup_env(proj)
@include_startup(proj)

#show(activeproject)

Main.include("src/AudioStick.jl")
Main.AudioStick.run_app()

ext: jl

max_stars_repo_name: ma-laforge/AudioStick

max_issues_repo_name: ma-laforge/AudioStick

max_forks_repo_name: ma-laforge/AudioStick


text:
module Reload

# package code goes here

println("\x1b[34m\ ")
print("Reload.jl loaded")
println("\x1b[00m\ ")

function watch(filename)

    start = stat(filename).mtime

    while start === stat(filename).mtime

        sleep(1);

    end

    eval(parse("try include(\"$(filename)\");catch e; println(\"\x1b[31m\ ERROR\");print(e);println(\"\x1b[00m\ \"); end;start(\"$(filename)\")"))

end

function start(filename)

    if (isfile(filename) == false)
        error("No such file: "*filename)
        return
    end

    @async (sleep(5); watch(filename))

    println("\x1b[34m\ ")
    println("watching $(filename)")
    println("\x1b[00m\ ")

end

end # module


text:
@testset "Arithmethic" begin
    @testset "Mag" begin
        @test Mag(4) <= Mag(1) + Mag(3) <= Mag(5)
        @test Mag(2) <= Mag(3) - Mag(1) <= Mag(3)
        @test Mag(0) <= Mag(3) - Mag(5) <= Mag(1)
        @test Mag(6) <= Mag(2) * Mag(3) <= Mag(7)
        @test Mag(3) <= Mag(6) / Mag(2) <= Mag(4)

        @test Mag(4) <= Mag(1) + 3 == 1 + Mag(3) <= Mag(5)
        @test Mag(6) <= Mag(2) * 3 == 3 * Mag(2) <= Mag(7)
        @test Mag(2) <= Mag(4) / 2 <= Mag(3)

        @test Mag(1 / 4) <= inv(Mag(4)) <= Mag(1 / 3)
        @test Mag(8) <= Mag(2)^3 <= Mag(9)
        @test Mag(2) <= sqrt(Mag(4)) <= Mag(3)
        @test Mag(1 / 2) <= Arblib.rsqrt(Mag(4)) <= Mag(1)
        @test Mag(5) <= hypot(Mag(3), Mag(4)) <= Mag(6)
        @test Mag(2) <= Arblib.root(Mag(8), 3) <= Mag(3)
        @test Mag(log(8)) <= log(Mag(8)) <= Mag(log(9))
        @test log(Mag(2)) <= Arblib.neglog(inv(Mag(2.0))) <= log(Mag(3))
        @test log(Mag(2)) <= log1p(Mag(1)) <= log(Mag(3))
        @test Mag(exp(2)) <= exp(Mag(2)) <= Mag(exp(3))
        @test Mag(exp(-2)) <= Arblib.expinv(Mag(2)) <= Mag(exp(-1))
        @test Mag(expm1(2)) <= expm1(Mag(2)) <= Mag(expm1(3))
        @test Mag(atan(2)) <= atan(Mag(2)) <= Mag(atan(3))
        @test Mag(cosh(2)) <= cosh(Mag(2)) <= Mag(cosh(3))
        @test Mag(sinh(2)) <= sinh(Mag(2)) <= Mag(sinh(3))

        @test min(Mag(1), Mag(2)) == Mag(1)
        @test max(Mag(1), Mag(2)) == Mag(2)
        @test minmax(Mag(1), Mag(2)) == minmax(Mag(2), Mag(1)) == (Mag(1), Mag(2))
        @test minimum(Mag[10:20; 0:9]) == Mag(0)
        @test maximum(Mag[10:20; 0:9]) == Mag(20)
        @test extrema(Mag[10:20; 0:9]) == (Mag(0), Mag(20))
    end

    @testset "Arf" begin
        @test sign(Arf(2)) == Arf(1)
        @test sign(Arf(-2)) == Arf(-1)
        @test sign(Arf(0)) == Arf(0)

        @test abs(Arf(-2)) == abs(Arf(2)) == Arf(2)
        @test -Arf(2) == Arf(-2)

        @test Arf(1) + Arf(2) ==
              Arf(1) + 2 ==
              2 + Arf(1) ==
              Arf(1) + UInt(2) ==
              UInt(2) + Arf(1) ==
              Arf(1) + UInt8(2) ==
              Arf(1) + UInt8(2) ==
              UInt8(2) + Arf(1) ==
              Arf(3)
        @test Arf(1) - Arf(2) == Arf(1) - 2 == Arf(1) - UInt(2) == Arf(-1)
        @test Arf(2) * Arf(3) ==
              Arf(2) * 3 ==
              3 * Arf(2) ==
              Arf(2) * UInt(3) ==
              UInt(3) * Arf(2) ==
              Arf(6)
        @test Arf(6) / Arf(2) ==
              Arf(6) / 2 ==
              6 / Arf(2) ==
              Arf(6) / UInt(2) ==
              UInt(6) / Arf(2) ==
              Arf(6) / UInt8(2) ==
              UInt8(6) / Arf(2) ==
              Arf(3)

        @test sqrt(Arf(4)) == Arf(2)
        @test Arblib.rsqrt(Arf(4)) == Arf(1 // 2)
        @test Arblib.root(Arf(8), 3) == Arf(2)

        @test min(Arf(1), Arf(2)) == Arf(1)
        @test max(Arf(1), Arf(2)) == Arf(2)
        @test minmax(Arf(1), Arf(2)) == minmax(Arf(2), Arf(1)) == (Arf(1), Arf(2))
        @test minimum(Arf[10:20; 0:9]) == Arf(0)
        @test maximum(Arf[10:20; 0:9]) == Arf(20)
        @test extrema(Arf[10:20; 0:9]) == (Arf(0), Arf(20))
    end

    @testset "$T" for T in [Arb, Acb]
        @test T(1) + T(2) ==
              T(1) + 2 ==
              2 + T(1) ==
              T(1) + UInt(2) ==
              UInt(2) + T(1) ==
              T(1) + UInt8(2) ==
              UInt8(2) + T(1) ==
              T(3)
        @test T(1) - T(2) == T(1) - 2 == T(1) - UInt(2) == T(-1)
        @test T(2) * T(3) ==
              T(2) * 3 ==
              3 * T(2) ==
              T(2) * UInt(3) ==
              UInt(3) * T(2) ==
              T(6)
        @test T(6) / T(2) == T(6) / 2 == T(6) / UInt(2) == T(3)

        @test Base.literal_pow(^, T(2), Val(-2)) ==
              T(2)^-2 ==
              Arblib.sqr(inv(T(2))) ==
              T(1 // 4)
        @test Base.literal_pow(^, T(2), Val(0)) == T(2)^0 == one(T)
        @test Base.literal_pow(^, T(2), Val(1)) == T(2)^1 == T(2)
        @test Base.literal_pow(^, T(2), Val(2)) == T(2)^2 == Arblib.sqr(T(2)) == T(4)

        @test Arblib.root(T(8), 3) == T(2)

        for f in [
            inv,
            sqrt,
            log,
            log1p,
            exp,
            expm1,
            sin,
            cos,
            tan,
            cot,
            sec,
            csc,
            atan,
            asin,
            acos,
            sinh,
            cosh,
            tanh,
            coth,
            sech,
            csch,
            atanh,
            asinh,
        ]
            @test f(T(0.5)) ≈ f(0.5)
        end
        @test acosh(T(2)) ≈ acosh(2)

        @test Arblib.rsqrt(T(1 // 4)) == T(2)
        @test Arblib.sqr(T(3)) == T(9)

        @test sinpi(T(1)) == T(0)
        @test cospi(T(1)) == T(-1)
        @test Arblib.tanpi(T(1)) == T(0)
        @test Arblib.cotpi(T(0.5)) == 0
        @test Arblib.cscpi(T(0.5)) == 1
        @test sinc(T(0.5)) ≈ sinc(0.5)

        @test isequal(sincos(T(1)), (sin(T(1)), cos(T(1))))
        if VERSION >= v"1.6"
            @test isequal(sincospi(T(1)), (sinpi(T(1)), cospi(T(1))))
        else
            @test isequal(Arblib.sincospi(T(1)), (sinpi(T(1)), cospi(T(1))))
        end
        @test isequal(Arblib.sinhcosh(T(1)), (sinh(T(1)), cosh(T(1))))
    end

    @testset "Arb - specific" begin
        @test Arb(1) + Arf(2) == Arf(2) + Arb(1) == Arb(3)
        @test Arb(1) - Arf(2) == Arb(-1)
        @test Arb(2) * Arf(3) == Arf(3) * Arb(2) == Arb(6)
        @test Arb(6) / Arf(2) == UInt(6) / Arb(2) == UInt8(6) / Arb(2) == Arb(3)

        @test Arb(2)^Arb(3) ==
              Arb(2)^3 ==
              Arb(2)^Int8(3) ==
              Arb(2)^UInt(3) ==
              Arb(2)^UInt8(3) ==
              (Arb(2)^-3)^-1 ==
              Arb(8)
        @test hypot(Arb(3), Arb(4)) == Arb(5)

        @test Arblib.sqrtpos(Arb(4)) == Arb(2)
        @test Arblib.sqrtpos(Arb(-4)) == Arb(0)
        @test Arblib.sqrt1pm1(Arb(3)) == Arb(1)

        @test atan(Arb(2), Arb(3)) ≈ atan(2, 3)

        @test min(Arb(1), Arb(2)) == Arb(1)
        @test max(Arb(1), Arb(2)) == Arb(2)
        @test minmax(Arb(1), Arb(2)) == minmax(Arb(2), Arb(1)) == (Arb(1), Arb(2))
        @test Arblib.contains(min(Arb((0, 2)), Arb((-1, 3))), -1)
        @test Arblib.contains(min(Arb((0, 2)), Arb((-1, 3))), 2)
        @test !Arblib.contains(min(Arb((0, 2)), Arb((-1, 3))), 3)
        @test Arblib.contains(max(Arb((0, 2)), Arb((-1, 3))), 0)
        @test Arblib.contains(max(Arb((0, 2)), Arb((-1, 3))), 3)
        @test !Arblib.contains(max(Arb((0, 2)), Arb((-1, 3))), -1)
        @test all(Arblib.contains.(minmax(Arb((0, 2)), Arb((-1, 3))), (-1, 0)))
        @test all(Arblib.contains.(minmax(Arb((0, 2)), Arb((-1, 3))), (2, 3)))
        @test all(.!Arblib.contains.(minmax(Arb((0, 2)), Arb((-1, 3))), (3, -1)))
        @test minimum(Arb[10:20; 0:9]) == Arb(0)
        @test maximum(Arb[10:20; 0:9]) == Arb(20)
        @test extrema(Arb[10:20; 0:9]) == (Arb(0), Arb(20))
        A = [Arb((i, i + 1)) for i = 0:10]
        @test Arblib.contains(minimum(A), Arb((0, 1)))
        @test Arblib.contains(minimum(reverse(A)), Arb((0, 1)))
        @test Arblib.contains(maximum(A), Arb((10, 11)))
        @test Arblib.contains(maximum(reverse(A)), Arb((10, 11)))
        @test all(Arblib.contains.(extrema(A), (Arb((0, 1)), Arb((10, 11)))))
        # These fails with the default implementation of minimum and maximum
        @test Arblib.contains(
            minimum([Arb((-i, -i + 1)) for i = 0:1000]),
            Arb((-1000, -999)),
        )
        @test Arblib.contains(maximum([Arb((i, i + 1)) for i = 0:1000]), Arb((1000, 1001)))
    end

    @testset "Acb - specific" begin
        @test Acb(1) + Arb(2) == Arb(2) + Acb(1) == Acb(3)
        @test Acb(1) - Arb(2) == Acb(-1)
        @test Acb(2) * Arb(3) == Arb(3) * Acb(2) == Acb(6)
        @test Acb(6) / Arb(2) == Acb(3)

        @test Acb(2)^Acb(3) ==
              Acb(2)^Arb(3) ==
              Acb(2)^3 ==
              Acb(2)^Int8(3) ==
              Acb(2)^UInt(3) ==
              Acb(2)^UInt8(3) ==
              (Acb(2)^-3)^-1 ==
              Acb(8)

        @test Acb(2, 3) * Complex{Bool}(0, 0) ==
              Complex{Bool}(0, 0) * Acb(2, 3) ==
              Acb(2, 3) * (0 + 0im)
        @test Acb(2, 3) * Complex{Bool}(0, 1) ==
              Complex{Bool}(0, 1) * Acb(2, 3) ==
              Acb(2, 3) * (0 + 1im)
        @test Acb(2, 3) * Complex{Bool}(1, 0) ==
              Complex{Bool}(1, 0) * Acb(2, 3) ==
              Acb(2, 3) * (1 + 0im)
        @test Acb(2, 3) * Complex{Bool}(1, 1) ==
              Complex{Bool}(1, 1) * Acb(2, 3) ==
              Acb(2, 3) * (1 + 1im)

        @test real(Acb(1, 2)) isa Arb
        @test real(Acb(1, 2)) == Arb(1)
        @test imag(Acb(1, 2)) isa Arb
        @test imag(Acb(1, 2)) == Arb(2)

        @test conj(Acb(1, 2)) == Acb(1, -2)

        @test abs(Acb(3, 4)) isa Arb
        @test abs(Acb(3, 4)) == Arb(5)
    end
end


text:
"""
Cost model typically used in an Objective function.
"""
abstract type CostModel end

"""
Cost that does not rely on any particular discretization.
"""
abstract type ContinuousCost <: CostModel end

"""
Punctual cost, that is only evaluated at designated points.
"""
abstract type PointCost <: CostModel end

"""
Linear Quadratic Regulator model.
"""
struct LQR <: ContinuousCost
    Q
    R
    nominal_state
end

function LQR(Q, R, model::Spacecraft; target_location=zeros(3))
    x_nom = RobotDynamics.build_state(model, target_location, UnitQuaternion(I), zeros(3), zeros(3))
    LQR(Q, R, x_nom)
end

TrajectoryOptimization.cost(cm::LQR) = LQRCost(cm.Q, cm.R, cm.nominal_state)

"""
LQR cost associated with waypoints, only at node points designed by integer indices in `times`.
"""
struct WaypointsLQR <: PointCost
    times
    costs
    function WaypointsLQR(waypoints, Q, R, model)
        costs = map(waypoints) do wp
            @show wp
            cost(LQR(Q, R, model; target_location=wp.location))
        end
        new(getproperty.(waypoints, :time), costs)
    end
end

"""
Cost structure holding a function that can be evaluated with the number of discretization points as argument to create a vector of costs.
"""
struct DiscretizeableCost <: CostModel
    discretization_f
end

function combine(cc::ContinuousCost, pc::PointCost)
    function discretize_cost(N)
        cost_1 = cost(cc)
        costs_all = map(1:N) do k
            i = findfirst(pc.times .== k)
            !isnothing(i) ? pc.costs[i] : cost_1
        end
    end
    DiscretizeableCost(discretize_cost)
end

"""
Structure holding a vector of point costs.
"""
struct TrajectoryCost
    costs
end
TrajectoryCost(N, cost_model::ContinuousCost) = TrajectoryCost(repeat([cost(cost_model)], N))
TrajectoryCost(N, cost_model::DiscretizeableCost) = TrajectoryCost(cost_model.discretization_f(N))

TrajectoryOptimization.Objective(t::TrajectoryCost) = Objective(t.costs)


text:
struct Constraints{TA,TB} <: AbstractConstraints
    before_update!::TA
    after_update!::TB
    n_before::Int
    n_after::Int
end

"""
    trivial_constraint!

Trival constraint function
"""
function trivial_constraint!(state::TS) where {TS<:AbstractSamplerState}
    state
end

"""
    Constraints

Set constraints on the sampler, if desired
### Fields
* `before! = trivial_constraint!` - Constraint to apply before a step of the sampler
* `after! = trivial_constraint!`  - Constraint to apply after a step of the sampler
* `n_before = 1` - Perform `before!` at `0, n_before-1, 2 n_before-1,...`
* `n_after = 1` - Perform `after!` at `n_after, 2 n_after,...`
"""

function Constraints(; before!::TA = trivial_constraint!, after!::TB = trivial_constraint!, n_before::TI = 1, n_after::TI = 1) where {TA,TB,TI<:Integer}
    return Constraints(before!, after!, n_before, n_after)
end


text:
module CGPrism3D
using LinearAlgebra
#using StatsBase

numchop(num::Real)= round(num,digits=14)#(abs(num) >= 1000*eps() ? num : zero(Real))

function numchop(num::Complex)
    numchop(imag(num)) == 0 ? numchop(real(num)) : Complex(numchop(real(num)), numchop(imag(num)))
end

#Define all functions needed for G-symbol
# Define delta{ijk} -> coupling rules
function delta(i::Float64,j::Float64,k::Float64,K::Int64)::Int64
    sol = 0
    if i <=(j+k) && j<=(i+k) && k<=(i+j) && i+j+k <= K && 2*(i+j+k)%2 ==0
        sol = 1
    end
    return sol
end

#Define quantum numbers qn (this is real)
function qn(n::Float64,K::Int64)
    sol = (exp(pi*n*im/(K+2)) - exp(-pi*n*im/(K+2))) / (exp(pi*im/(K+2)) - exp(-pi*im/(K+2)))
    return real(sol)
end

#Define qn factorial
function qnfact(n::Float64,K::Int64)
    sol = 1.
    for i in 1:n
        sol *= qn(i,K)
    end
    return sol
end

#Define square root of quantum dimension
function visqrt(i::Float64,K::Int64)
    sol = ((-1+0im)^i )*sqrt(qn(2*i+1,K))
    return sol
end

#Define triangle equality
function trian(i::Float64,j::Float64,k::Float64,K::Int64)
    sol = 0.
    if delta(i,j,k,K) == 1
        sol = delta(i,j,k,K)*sqrt(qnfact(i+j-k,K)*qnfact(i-j+k,K)*qnfact(-i+j+k,K)/qnfact(i+j+k+1,K))
    end
    return sol
end

# Define Racah-Wigner six-j symbol
function RacWig6j(i::Float64,j::Float64,m::Float64,k::Float64,l::Float64,n::Float64,K::Int64)
    a = i+j+m; b = i+l+n; c = k+j+n; d = k+l+m;  e = i+j+k+l; f = i+k+m+n; g = j+l+m+n
    sol = 0.
    if delta(i,j,m,K) != 0 && delta(i,l,n,K) != 0 && delta(k,j,n,K) != 0 && delta(k,l,m,K) != 0
        sumz = 0.
        for z in max(a,b,c,d):min(e,f,g)
            sumz += (-1)^z *qnfact(z+1,K)/
                ((qnfact(e-z,K)*qnfact(f-z,K)*qnfact(g-z,K))* (qnfact(z-a,K)*qnfact(z-b,K)*qnfact(z-c,K)*qnfact(z-d,K)))
        end
        sol = trian(i,j,m,K)*trian(i,l,n,K)*trian(k,j,n,K)*trian(k,l,m,K)*sumz
    end
    return sol
end

#Define F-symbol
function Fsymb(i::Float64,j::Float64,m::Float64,k::Float64,l::Float64,n::Float64,K::Int64)
    sol = 0.
    if delta(i,j,m,K) != 0 && delta(i,l,n,K) != 0 && delta(k,j,n,K) != 0 && delta(k,l,m,K) != 0
        sol = (-1+0im)^(i+j+k+l)*sqrt(qn(2*m+1,K)*qn(2*n+1,K)) * RacWig6j(i,j,m,k,l,n,K)
    end
    return sol
end

#Define G-symbol
function Gsymb(i::Float64,j::Float64,m::Float64,k::Float64,l::Float64,n::Float64,K::Int64)
    sol = 0.
    if delta(i,j,m,K) != 0 && delta(i,l,n,K) != 0 && delta(k,j,n,K) != 0 && delta(k,l,m,K) != 0
        sol = Fsymb(i,j,m,k,l,n,K) /(visqrt(m,K)*visqrt(n,K))
    end
    return sol
end

export Tetra6j, delta, visqrt
#Define G-symbol times dimensions
function Tetra6j(i::Float64,j::Float64,m::Float64,k::Float64,l::Float64,n::Float64,K::Int64)
    sol = 0.
    if delta(i,j,m,K) != 0 && delta(i,l,n,K) != 0 && delta(k,j,n,K) != 0 && delta(k,l,m,K) != 0
        dims = (prod(visqrt.([i,j,m,k,l,n],K)))
        sol =  (dims*Gsymb(i,j,m,k,l,n,K))
    end
    return sol
end

#Define G-symbol 2
function Gsymb2(i::Float64,j::Float64,m::Float64,k::Float64,l::Float64,n::Float64,K::Int64)
    sol = 0.
    if delta(i,j,m,K) != 0 && delta(i,l,n,K) != 0 && delta(k,j,n,K) != 0 && delta(k,l,m,K) != 0
        sol = (-1+0im)^(i+j+k+l+m+n)*RacWig6j(i,j,m,k,l,n,K)
    end
    return sol
end

#Define square root of quantum dimension
function visqrtT(a::Float64,J::Int64,K::Int64,alpha::Float64)
    ((-1+0im)^a )*sqrt(abs(TetraJK(0.,0.,0.,a,a,a,J,K,alpha)))
end

function visqrtU(Utensor,a::Float64)
    indx = Utensor[1]
    amps = Utensor[2]
    fc = findfirst(x-> x == [0.,0.,0.,a,a,a], indx)
    ans = 0
    if typeof(fc) == Int64
        ans = ((-1+0im)^a )*sqrt(abs(amps[fc]))
    end
    ans
    1
end

function TetraJK(i::Float64,j::Float64,m::Float64,k::Float64,l::Float64,n::Float64,J::Int64,K::Int64,alpha::Float64)
    #@assert J <= K
    alpha*Tetra6j(i,j,m,k,l,n,J) +(1-alpha)*Tetra6j(i,j,m,k,l,n,K)
end

function prismA(j1::Float64,j2::Float64,j3::Float64,j4::Float64,j5::Float64,j6::Float64,j7::Float64,
                j8::Float64,j9::Float64,j10::Float64,j11::Float64,j12::Float64,K::Int64)
    sol = 0
    if delta(j1,j2,j3,K) != 0 && delta(j1,j5,j6,K) != 0 && delta(j2,j4,j6,K) != 0 && delta(j3,j4,j5,K) != 0 && delta(j1,j8,j9,K) != 0 && delta(j3,j7,j8,K) != 0 && delta(j2,j7,j9,K) != 0 && delta(j9,j10,j11,K) != 0 && delta(j7,j10,j12,K) != 0 && delta(j2,j11,j12,K) != 0
        dims = prod(visqrt.([j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12],K))
        sol =  dims*Gsymb(j1,j2,j3,j4,j5,j6,K) * Gsymb(j1,j2,j3,j7,j8,j9,K) * Gsymb(j2,j7,j9,j10,j11,j12,K)
    end
    return sol
end

export dataTet, dataFsymb, dataPrsmA, dataGsymb

function dataPrsmA(J::Int64,K::Int64,alpha::Float64)
    @assert J <= K
    @assert 0. <= alpha <= 1.
    x = K+1
    y = K/2
    ampsInfo = Array{Float64,1}[]
    amps = Float64[]
    for j1 in 0.:0.5:y, j2 in 0.:0.5:y, j3 in 0.:0.5:y
        if delta(j1,j2,j3,K) != 0
            for j4 in 0.:0.5:y, j5 in 0.:0.5:y
                if delta(j3,j4,j5,K) != 0
                    for j6 in 0.:0.5:y
                        if delta(j1,j5,j6,K) != 0 && delta(j2,j4,j6,K) != 0
                            for j7 in 0.:0.5:y, j8 in 0.:0.5:y
                                if delta(j3,j7,j8,K) != 0
                                    for j9 in 0.:0.5:y
                                        if delta(j1,j8,j9,K) != 0 && delta(j2,j7,j9,K) != 0
                                            for j10 in 0.:0.5:y, j11 in 0.:0.5:y
                                                if delta(j9,j10,j11,K) != 0
                                                    for j12 in 0.:0.5:y
                                                        if delta(j7,j10,j12,K) != 0 && delta(j2,j11,j12,K) != 0 #
                                                            sol = numchop(alpha*prismA(j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12,J)+(1-alpha)*prismA(j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12,K))
                                                            if sol != 0
                                                                indx = [j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12]
                                                                push!(ampsInfo,indx)
                                                                push!(amps,sol)
                                                            end
                                                        end
                                                    end
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return ampsInfo, amps
end

function dataTet(J::Int64,K::Int64,alpha::Float64)
    @assert J <= K
    @assert 0. <= alpha <= 1.
    x = K+1
    y = K/2
    t6j = Array{Float64,1}[]
    t6s = Float64[]
    for ja in 0.:0.5:y, jb in 0.:0.5:y, jc in 0.:0.5:y
        if delta(ja,jb,jc,K) != 0
            for jd in 0.:0.5:y, je in 0.:0.5:y
                if delta(jc,jd,je,K) != 0
                    for jf in 0.:0.5:y
                        if delta(ja,je,jf,K) != 0 && delta(jb,jd,jf,K) != 0
            #dims = prod(visqrt.([ja,jb,jc,jd,je,jf]))
                            sol = numchop(TetraJK(ja,jb,jc,jd,je,jf,J,K,alpha))
                            #if sol != 0
                            push!(t6s,sol)
                            push!(t6j,[ja,jb,jc,jd,je,jf])
                            #end
                        end
                    end
                end
            end
        end
    end
    return t6j,t6s
end

function dataFsymb(J::Int64,K::Int64,alpha::Float64)
    @assert J <= K
    @assert 0. <= alpha <= 1.
    x = K+1
    y = K/2
    t6j = Array{Float64,1}[]
    t6s = Float64[]
    for ja in 0.:0.5:y, jb in 0.:0.5:y, jc in 0.:0.5:y
        if delta(ja,jb,jc,K) != 0
            for jd in 0.:0.5:y, je in 0.:0.5:y
                if delta(jc,jd,je,K) != 0
                    for jf in 0.:0.5:y
                        if delta(ja,je,jf,K) != 0 && delta(jb,jd,jf,K) != 0
            #dims = prod(visqrt.([ja,jb,jc,jd,je,jf]))
                            sol = numchop(alpha*Fsymb(ja,jb,jc,jd,je,jf,J) +(1-alpha)*Fsymb(ja,jb,jc,jd,je,jf,K))
                            if sol != 0
                                push!(t6s,sol)
                                push!(t6j,[ja,jb,jc,jd,je,jf])
                            end
                        end
                    end
                end
            end
        end
    end
    return t6j,t6s
end

function dataGsymb(J::Int64,K::Int64,alpha::Float64)
    @assert J <= K
    @assert 0. <= alpha <= 1.
    x = K+1
    y = K/2
    t6j = Array{Float64,1}[]
    t6s = Float64[]
    for ja in 0.:0.5:y, jb in 0.:0.5:y, jc in 0.:0.5:y
        if delta(ja,jb,jc,K) != 0
            for jd in 0.:0.5:y, je in 0.:0.5:y
                if delta(jc,jd,je,K) != 0
                    for jf in 0.:0.5:y
                        if delta(ja,je,jf,K) != 0 && delta(jb,jd,jf,K) != 0
            #dims = prod(visqrt.([ja,jb,jc,jd,je,jf]))
                            sol = numchop(alpha*Gsymb(ja,jb,jc,jd,je,jf,J) +(1-alpha)*Gsymb(ja,jb,jc,jd,je,jf,K))
                            if sol != 0
                                push!(t6s,sol)
                                push!(t6j,[ja,jb,jc,jd,je,jf])
                            end
                        end
                    end
                end
            end
        end
    end
    return t6j,t6s
end

function tensorBlockOld(tensorM,posnA,posnB,posnC,spinC)
    indx = tensorM[1]
    ampvals = tensorM[2]
    fc = findall(x-> x[posnC] == spinC, indx) #find spinsC in all spins length(posnC) = length(spinC)
    ampInfo = indx[fc] # this contains all spin configuration containing ampsC
    amps = ampvals[fc] # this contains the corresponding spin amplitude values for ampsC
    Acol = sort(unique(getindex.(ampInfo, [posnA])))
    Brow = sort(unique(getindex.(ampInfo, [posnB])))
    indxU = vcat.(Acol,repeat([spinC],length(Acol)))
    indxV = vcat.(Brow,repeat([spinC],length(Brow)))
    #indxU = []
    #indxV = []
    matM = zeros(length(Acol),length(Brow))
    if length(matM) != 0
        for i in 1:length(ampInfo)
            qu = findfirst(x-> x == ampInfo[i][posnA], Acol)[1]# [1] since it returns an array [1]
            qv = findfirst(x-> x == ampInfo[i][posnB], Brow)[1]
            matM[qu,qv] = amps[i]
            #push!(indxU,vcat(Acol[qu],spinC))
            #push!(indxV,vcat(Brow[qv],spinC))
        end
        #U, s, V = svd(matM)
        println(matM)
        return matM,indxU,indxV
    else
        return 0,0,0
    end
end

export tensorBlock

function tensorBlock(tensorM,posnA,posnB,posnC,spinC)
    indx = tensorM[1]
    ampvals = tensorM[2]
    fc = findall(x-> x[posnC] == spinC, indx) #find spinsC in all spins length(posnC) = length(spinC)
    indxI = indx[fc] # this contains all spin configuration containing ampsC
    ampsI = ampvals[fc] # this contains the corresponding spin amplitude values for ampsC
    Acol = sort(unique(getindex.(indxI, [posnA])))
    Brow = sort(unique(getindex.(indxI, [posnB])))
    indxU = vcat.(Acol,repeat([spinC],length(Acol)))
    indxV = vcat.(Brow,repeat([spinC],length(Brow)))
    #indxU = []
    #indxV = []
    #println(ampInfo)
    la = length(Acol)
    lb = length(Brow)
    lpa = length(posnA)
    lpb = length(posnB)
    matM = zeros(la,lb)
    if la != 0 || lb !=0
        #if length(ampsI) == la*lb
        indxM = []
        #pos = collect(Iterators.flatten(([posnA,posnB,posnC])))
        pos = vcat(posnA,posnB,posnC)
        for i in indxI
            push!(indxM,i[pos])
        end
        sp = sortperm(indxM)
        ampInfo = indxM[sp]
        amps = ampsI[sp]
        for i in 1:la, j in 1:lb
            if ampInfo[1][1:lpa] == Acol[i] && ampInfo[1][lpa+1:lpb+lpa] == Brow[j]
                #println("sucess")
                matM[i,j] = amps[1]
                deleteat!(amps,1)
                deleteat!(ampInfo,1)
            end
        end
        return matM,indxU,indxV#,ampInfo
    else
        return 0,0,0
    end
    #return ampInfo, matM
end

export tensorGlueTet3DO, tensorGlueTet3DN, tensorGlue

function tensorGlueTet3DO(tensorB,tensorA,posnB,posnA,Utensor)
    indx, amps = tensorGlue(tensorB,tensorA,posnB,posnA)
    face = getindex.(indx, [posnB])
    ans = complex(ones(length(face)))
    #visqt = visqJK(a,J,K,alpha)
    for i in 1:length(face)
        #ans[i] = prod(visqrtT.(face[i],J,K,alpha))
        pp = 1
        for j in face[i]
            pp *= visqrtU(Utensor,j)
        end
        ans[i] = pp
    end
    ampsN = @. amps / ans
    ampsN = real(ampsN)+imag(ampsN)
    return indx, ampsN
end

function tensorGlue(tensorB,tensorA,posnB,posnA)# glue two tensors TA,TB along posnA from TA and posnB from TB
    indxA = tensorA[1] # spins of TA
    indxB = tensorB[1] # spins of TB
    ampsA = tensorA[2] # amplitudes of TA
    ampsB = tensorB[2] # amplitudes of TB
    shrdA = getindex.(indxA,[posnA])
    shrdB = getindex.(indxB,[posnB])
    lena = collect(1:length(indxA[1]))
    deleteat!(lena,sort(posnA)) # assert needs posnA to be a vector
    #if length(posnA) > 1 deleteat!(lena,sort(posnA,dims=2) else  deleteat!(lena,sort(posnA)) end
    shrdAB = Tuple{Array{Float64,1},Int64,Int64}[]
    qq = Array{Array{Float64,1},1}[]
    amps = Array{Float64,1}[]
    sa = sortperm(shrdA)
    sb = sortperm(shrdB)
    shrdA = shrdA[sa]
    shrdB = shrdB[sb]
    sindxA = indxA[sa]
    sindxB = indxB[sb]
    sampsA = ampsA[sa]
    sampsB = ampsB[sb]
    shrdAu = unique(shrdA)
    shrdBu = unique(shrdB)
    shrd = shrdAu
    #println(shrdAu == shrdBu)
    if shrdAu == shrdBu
        shrd = shrdAu
    else
        for i in shrdBu
            shrd = shrdAu
            if !(i in shrdAu)
                push!(shrd,i)
            end
        end
    end
    #cshrdA = countmap(shrdA)
    #cshrdB = countmap(shrdB)
    for i in shrd
        ca = count(x->x==i,shrdA)
        cb = count(x->x==i,shrdB)
        push!(shrdAB, (i, ca,cb ))
    end
    #println(shrdAB)
    cnt = 1
    cntb = 1
    for i in shrdAB
        if i[2] == 0 || i[3] == 0 # Put condition for when one of a or b is zero
            nothing
        else
            for j in 1:i[2]
                ans = @. sampsA[j]*sampsB[1:i[3]]
                push!(amps,ans)
                indxa = repeat([sindxA[j][lena]],i[3])
                #for j in indxA[i][lena]
                indxa1 = vcat.(sindxB[1:i[3]],indxa)
                push!(qq,indxa1)
            end
        end
        deleteat!(sindxA,1:i[2])
        deleteat!(sindxB,1:i[3])
        deleteat!(sampsA,1:i[2])
        deleteat!(sampsB,1:i[3])
        #cntb = cnt
    end
    indxsol = collect(Iterators.flatten(qq))
    sol = collect(Iterators.flatten(amps))
    return indxsol,sol
end

export tensorGlueTet3DN2

function tensorGlueTet3DN2(tensorB,tensorA,posnB,posnA,Utensor)# glue two tensors TA,TB along posnA from TA and posnB from TB
    indxA = tensorA[1] # spins of TA
    indxB = tensorB[1] # spins of TB
    ampsA = tensorA[2] # amplitudes of TA
    ampsB = tensorB[2] # amplitudes of TB
    shrdA = getindex.(indxA,[posnA])
    shrdB = getindex.(indxB,[posnB])
    lena = collect(1:length(indxA[1]))
    deleteat!(lena,sort(posnA)) # assert needs posnA to be a vector
    #if length(posnA) > 1 deleteat!(lena,sort(posnA,dims=2) else  deleteat!(lena,sort(posnA)) end
    shrdAB = []#Tuple{Array{Float64,1},Int64,Int64}[]
    qq = []#Array{Array{Float64,1},1}[]
    amps = [] #Array{Float64,1}[]
    sa = sortperm(shrdA)
    sb = sortperm(shrdB)
    shrdA = shrdA[sa]
    shrdB = shrdB[sb]
    sindxA = indxA[sa]
    sindxB = indxB[sb]
    sampsA = ampsA[sa]
    sampsB = ampsB[sb]
    shrdAu = unique(shrdA)
    shrdBu = unique(shrdB)
    shrd = shrdAu
    #println(shrdAu == shrdBu)
    if shrdAu == shrdBu
        nothing
    else
        for i in shrdBu
            #shrd = shrdAu
            if !(i in shrdAu)
                push!(shrd,i)
            end
        end
    end
    for i in shrd
        cna = count(x->x==i,shrdA)
        cnb = count(x->x==i,shrdB)
        push!(shrdAB, [i, cna,cnb])
        #println("Ca,Cna = ",ca,"\t",cna)
        #println("Cb,Cnb = ",cb,"\t",cnb)
    end
    for i in shrdAB
        if i[2] == 0 || i[3] == 0 # Put condition for when one of a or b is zero
            nothing
        else
            for j in 1:i[2]
                for k in 1:i[3]
                    indx1 = vcat(sindxB[k],sindxA[j][lena])
                    pp=1
                    for k in i[1]
                        pp *= visqrtU(Utensor,k)
                    end
                    amps1 = (sampsA[j]*sampsB[k]/pp)
                    push!(qq,indx1)
                    push!(amps,amps1)
                    #ans = @. sampsA[j]*sampsB[1:i[3]]
                end
                #pp = 1
                #
                #ans /= pp
                #ampsN = numchop.(real(ans)+imag(ans))
                #push!(amps,ampsN)
                #indxa = repeat([sindxA[j][lena]],i[3])
                #for j in indxA[i][lena]
                #indxa1 = vcat.(sindxB[1:i[3]],indxa)
                #push!(qq,indxa1)
            end
        end
        deleteat!(sindxA,1:i[2])
        deleteat!(sindxB,1:i[3])
        deleteat!(sampsA,1:i[2])
        deleteat!(sampsB,1:i[3])
        #cntb = cnt
    end
    #indxsol = collect(Iterators.flatten(qq))
    #sol = collect(Iterators.flatten(amps))
    return qq, amps #indxsol,sol
end

function tensorGlueTet3DN(tensorB,tensorA,posnB,posnA,Utensor)# glue two tensors TA,TB along posnA from TA and posnB from TB
    indxA = tensorA[1] # spins of TA
    indxB = tensorB[1] # spins of TB
    ampsA = tensorA[2] # amplitudes of TA
    ampsB = tensorB[2] # amplitudes of TB
    shrdA = getindex.(indxA,[posnA])
    shrdB = getindex.(indxB,[posnB])
    lena = collect(1:length(indxA[1]))
    deleteat!(lena,sort(posnA)) # assert needs posnA to be a vector
    #if length(posnA) > 1 deleteat!(lena,sort(posnA,dims=2) else  deleteat!(lena,sort(posnA)) end
    shrdAB = []#Tuple{Array{Float64,1},Int64,Int64}[]
    qq = Array{Array{Float64,1},1}[]
    amps = Array{Float64,1}[]
    sa = sortperm(shrdA)
    sb = sortperm(shrdB)
    shrdA = shrdA[sa]
    shrdB = shrdB[sb]
    sindxA = indxA[sa]
    sindxB = indxB[sb]
    sampsA = ampsA[sa]
    sampsB = ampsB[sb]
    shrdAu = unique(shrdA)
    shrdBu = unique(shrdB)
    shrd = shrdAu
    #println(shrdAu == shrdBu)
    if shrdAu == shrdBu
        nothing
    else
        for i in shrdBu
            #shrd = shrdAu
            if !(i in shrdAu)
                push!(shrd,i)
            end
        end
    end
    for i in shrd
        cna = count(x->x==i,shrdA)
        cnb = count(x->x==i,shrdB)
        push!(shrdAB, [i, cna,cnb])
        #println("Ca,Cna = ",ca,"\t",cna)
        #println("Cb,Cnb = ",cb,"\t",cnb)
    end
    for i in shrdAB
        if i[2] == 0 || i[3] == 0 # Put condition for when one of a or b is zero
            nothing
        else
            for j in 1:i[2]
                ans = @. sampsA[j]*sampsB[1:i[3]]
                pp = 1
                for k in i[1]
                    pp *= visqrtU(Utensor,k)
                end
                ans /= pp
                ampsN = (real(ans)+imag(ans))
                push!(amps,ampsN)
                indxa = repeat([sindxA[j][lena]],i[3])
                #for j in indxA[i][lena]
                indxa1 = vcat.(sindxB[1:i[3]],indxa)
                push!(qq,indxa1)
            end
        end
        deleteat!(sindxA,1:i[2])
        deleteat!(sindxB,1:i[3])
        deleteat!(sampsA,1:i[2])
        deleteat!(sampsB,1:i[3])
        #cntb = cnt
    end
    indxsol = collect(Iterators.flatten(qq))
    sol = collect(Iterators.flatten(amps))
    return indxsol,numchop.(sol)
end

# can also improve tensorGlueTet3D maybe?
export fullSplitTet3DN

function fullSplitTet3DN(dataM,posnA,posnB,posnC,Utensor,J::Int64,K::Int64,alpha::Float64) # posnA,B,C must be vectors @assert length(dataM) = vcat(A,B,C)
    indx = dataM[1]
    ampvals = dataM[2]
    #posAC = vcat(posnA,posnC) # assert posnA,posnC must be both vectors
    #posBC = vcat(posnB,posnC) # assert posnB,posnC must be both vectors
    #indxU = unique(getindex.(indx, [posAC]))
    #indxV = unique(getindex.(indx, [posBC]))
    ampsC = unique(getindex.(indx, [posnC])) # get all unique spins C (jd1,jd2,je1)
    ampsU = [] # this will store all amplitudes for U^C_{A}
    ampsV = [] # this will store all amplitudes for V^C_{B}
    indxsU = []
    indxsV = []
    #indxUV = [] # want this to store all spins in the right order
    for i in 1:length(ampsC) # loop over the unique spins
        jd1 = ampsC[i][1]; jd2 = ampsC[i][2]; je1 = ampsC[i][3]  #jd1,jd2,je1 in spin C form a triangle
        if delta(jd1,jd2,je1,K) != 0
            spinC = [jd1,jd2,je1]
            block = tensorBlock(dataM,posnA,posnB,posnC,spinC)
            mat = block[1]
            blkU = block[2]
            blkV = block[3]
        #if mat != 0
            U, s, V = svd(mat)
            #truncate and use only first singular value
            trunc = 1#length(s)
            U1 = numchop.(U[:,1:trunc] )
            V1 = numchop.( V[:,1:trunc] )
            s1 = numchop.((s[1:trunc]) )
            pp = 1
            #println(ampsC[i])
            for j in ampsC[i]
                #println(j)
                pp *= visqrtU(Utensor,j)
            end
            #@show trunc
            valU = U1*sqrt.(s1)*sqrt(pp) # multiply by leftover dimension factors
            valV = V1*sqrt.(s1)*sqrt(pp)
            # valU,valV are all either real or purely imaginary
            valU = real(valU) - imag(valU) # take -ve of imaginary part if it gives only imag
            valV = real(valV) + imag(valV)
            ja1 = blkU[1][1]; jb1 = blkU[1][2]; jg = blkU[1][3] # pick the first spins in U
            #Fixing sign problem -- keep the same sign for U and Tetra6j symbol
            #fc = findfirst(x->x==[jd1,jd2,je1,ja1,jb1,jg],Utensor[1])
            #println(sign(Utensor[2][fc]),sign(valU[1]))
            if sign(valU[1]) == sign(TetraJK(jd1,jd2,je1,ja1,jb1,jg,J,K,alpha))
                push!(ampsU, valU)#/valU[1])
                push!(ampsV, valV)#/valV[1])
                push!(indxsU,blkU)
                push!(indxsV,blkV)
            else
                push!(ampsU, -valU)#/valU[1])
                push!(ampsV, -valV)#/valV[1])
                push!(indxsU,blkU)
                push!(indxsV,blkV)
            end
        end
    end
    #println("done")
    solU = collect(Iterators.flatten(ampsU))
    solV = collect(Iterators.flatten(ampsV))
    #@show (solU[1])
    ansU = solU/solU[1] # normalization condition
    ansV = solV/solV[1]
    indxUs = collect(Iterators.flatten(indxsU))
    indxVs = collect(Iterators.flatten(indxsV))
    return (indxUs,numchop.(ansU)),(indxVs,numchop.(ansV)) #,indxUV
end

export tensorSumTetO, tensorSumTetN, tensorSum

function tensorSumTetN(ampJ,posnN,Utensor)
    indx = ampJ[1]
    amps = ampJ[2]
    lenS = collect(1:length(indx[1]))
    deleteat!(lenS,sort(posnN))
    #if length(posnN) > 1   deleteat!(lenS,sort(posnN,dims=2)) else  deleteat!(lenS,sort(posnN)) end
    indxef = getindex.(indx,[lenS]) # get all spins without N
    #indxefN = getindex.(indx,[posnN])
    #indxefu = unique(indxef)
    tt = sortperm(indxef) # sort them spins
    indxef = indxef[tt] # apply sort on spins
    amps = amps[tt]  # apply sort on amps
    #fc = findall(x-> x[lenS]==ans ,indx)
    qq = []
    indq = []
    fc = sortperm(indxef)
    indxN = indxef[fc]
    amps = amps[fc]
    #indxN = indxefN[fc]
    if length(posnN) == 1
        blkdims = getindex.(indx,posnN)[fc]
        sqdm = []
        for j in blkdims
            push!(sqdm,visqrtU(Utensor,j))
        end
        ampsn = amps .*sqdm
        ampsN = real(ampsn) - imag(ampsn)
        #data = indx, ampsN2
        #ans = tensorSum(data,posnN)
    else
        blkdims = getindex.(indx,[posnN])[fc]
        sqdm = []
        for i in blkdims
            pp = 1
            for j in i
                pp *= (visqrtU(Utensor,j))
            end
            push!(sqdm,pp)
        end
        #println(sqdm)
        ampsn = amps .*sqdm
        ampsN = real(ampsn) - imag(ampsn)
        #data = indx, ampsN2
        #ans = tensorSum(data,posnN)
    end

    for i in 2:length(indxN)
        if indxN[i] == indxN[i-1]
            ampsN[i] += ampsN[i-1]
        elseif numchop(ampsN[i-1]) != 0
            push!(indq,indxN[i-1])
            push!(qq,ampsN[i-1])
        end
    end
    push!(indq,indxN[end])
    push!(qq,ampsN[end])
    return indq, qq
end

# Summing a tensor  over indices labelled posnN
function tensorSum(ampJ,posnN)
    indx = ampJ[1]
    amps = ampJ[2]
    lenS = collect(1:length(indx[1]))
    deleteat!(lenS,sort(posnN))
    #if length(posnN) > 1   deleteat!(lenS,sort(posnN,dims=2)) else  deleteat!(lenS,sort(posnN)) end
    indxef = getindex.(indx,[lenS]) # get all spins without N
    #indxefu = unique(indxef)
    tt = sortperm(indxef) # sort them spins
    indxef = indxef[tt] # apply sort on spins
    amps = amps[tt]  # apply sort on amps
    #fc = findall(x-> x[lenS]==ans ,indx)
    qq = []
    indq = []
    fc = sortperm(indxef)
    indxN = indxef[fc]
    ampsN = amps[fc]

    for i in 2:length(indxN)
        if indxN[i] == indxN[i-1]
            ampsN[i] += ampsN[i-1]
        elseif numchop(ampsN[i-1]) != 0
            push!(indq,indxN[i-1])
            push!(qq,ampsN[i-1])
        end
    end
    push!(indq,indxN[end])
    push!(qq,ampsN[end])
    return indq, numchop.(qq)
end

#maybe can also improve tensorSumTet

function tensorSumTetO(tensor,posnN,Utensor)
    indx = tensor[1]
    amps = tensor[2]
    if length(posnN) == 1
        blkdims = getindex.(indx,posnN)
        sqdm = []
        for j in blkdims
            push!(sqdm,visqrtU(Utensor,j))
        end
        ampsN = amps .*sqdm
        ampsN2 = real(ampsN) - imag(ampsN)
        data = indx, ampsN2
        ans = tensorSum(data,posnN)
    else
        blkdims = getindex.(indx,[posnN])
        sqdm = []
        for i in blkdims
            pp = 1
            for j in i
                pp *= (visqrtU(Utensor,j))
            end
            push!(sqdm,pp)
        end
        #println(sqdm)
        ampsN = amps .*sqdm
        ampsN2 = real(ampsN) - imag(ampsN)
        data = indx, ampsN2
        ans = tensorSum(data,posnN)
    end
    return ans
end

export tensor22moveT, tensor22moveF

function swap2(vec,a,b)
    vec[a],vec[b] = vec[b], vec[a]
    return vec
end

#Same tensor for gluing and adding dimesion factors
function tensor22moveT(tensor,posF,Utensor) # J,K,alpha determines what Fsymbol to use
    glu = tensorGlueTet3DN(tensor,Utensor,posF,[1,5,2,4,6],Utensor)
    n = posF[end]
    m = length(glu[1][1])
    swp = @. swap2(glu[1],n,m)
    tensorN = swp, glu[2]
    ans = tensorSumTetO(tensorN,[m],Utensor)
    return ans
end

export tensor22moveT2, tensor22moveT3

function tensor22moveF(tensor,posF,J::Int64,K::Int64,alpha::Float64)
    glu = tensorGlue(tensor,dataFsymb(J,K,alpha),posF,[1,5,2,4,6])
    n = posF[end]
    m = length(glu[1][1])
    swp = @. swap2(glu[1],n,m)
    tensorN = swp, glu[2]
    ans = tensorSum(tensorN,[m])
    return ans
end

function tensorSum2(ampJ,posnN,tri,K)
    indx = ampJ[1]
    amps = ampJ[2]
    lenS = collect(1:length(indx[1]))
    deleteat!(lenS,sort(posnN))
    #if length(posnN) > 1   deleteat!(lenS,sort(posnN,dims=2)) else  deleteat!(lenS,sort(posnN)) end
    indxef = getindex.(indx,[lenS]) # get all spins without N
    #indxefu = unique(indxef)
    tt = sortperm(indxef) # sort them spins
    indxef = indxef[tt] # apply sort on spins
    amps = amps[tt]  # apply sort on amps
    #fc = findall(x-> x[lenS]==ans ,indx)
    qq = []
    indq = []
    fc = sortperm(indxef)
    indxN = indxef[fc]
    ampsN = amps[fc]

    for i in 2:length(indxN)
        if delta(indxN[i-1][tri[1]],indxN[i-1][tri[2]],indxN[i-1][tri[3]],K) != 0.
            if indxN[i] == indxN[i-1]
                #println([indxN[i][tri[1]],indxN[i][tri[2]],indxN[i][tri[3]],K])
                ampsN[i] += ampsN[i-1]
            elseif numchop(ampsN[i-1]) != 0.
                push!(indq,indxN[i-1])
                push!(qq,ampsN[i-1])
            end
        end
    end
    push!(indq,indxN[end])
    push!(qq,ampsN[end])
    return indq, qq
end

function tensorSumTetO2(tensor,posnN,Utensor,tri,K)
    indx = tensor[1]
    amps = tensor[2]
    if length(posnN) == 1
        blkdims = getindex.(indx,posnN)
        sqdm = []
        for j in blkdims
            push!(sqdm,visqrtU(Utensor,j))
        end
    else
        blkdims = getindex.(indx,[posnN])
        sqdm = []
        for i in blkdims
            pp = 1
            for j in i
                pp *= numchop(visqrtU(Utensor,j))
            end
            push!(sqdm,pp)
        end
        #println(sqdm)
    end
    ampsN = amps .*sqdm
    ampsN2 = real(ampsN) - imag(ampsN)
    data = indx, ampsN2
    ans = tensorSum2(data,posnN,tri,K)
    return ans
end

function tensor22moveT2(tensor,posF,Utensor,tri,K::Int64) # J,K,alpha determines what Fsymbol to use
    glu = tensorGlueTet3DN(tensor,Utensor,posF,[1,5,2,4,6],Utensor)
    n = posF[end]
    m = length(glu[1][1])
    swp = @. swap2(glu[1],n,m)
    tensorN = swp, glu[2]
    #if length(tri) == 3
    ans = tensorSumTetO2(tensorN,[m],Utensor,tri,K)
    #elseif length(tri) == 0
    #    ans = tensorSumTetO(tensorN,[m],Utensor)
    #end
    return ans
end

function tensor22moveT3(tensor,posF,Utensor,tri,K::Int64) # J,K,alpha determines what Fsymbol to use
    glu = tensorGlueTet3DN(tensor,Utensor,posF,[1,5,2,4,6],Utensor)
    n = posF[end]
    m = length(glu[1][1])
    swp = @. swap2(glu[1],n,m)
    tensorN = swp, glu[2]
    indx, amps = tensorSumTetO(tensorN,[m],Utensor)
    ls = []
    la = []
    #if length(tri) == 3
    for i in 1:length(indx)
        if delta(indx[i][tri[1]],indx[i][tri[2]],indx[i][tri[3]],K) != 0
            push!(ls,indx[i])
            push!(la,amps[i])
        end
    end
    return ls,la
    #elseif length(tri) == 0
    #    return indx, amps
    #end
end

export permuteInd, permuteIndN, tensorPermute, tensorPermuteN

function permuteInd(vec)
    perm = [5,1,11,2,6,3,10,12,7,8,4,9]
    nvec = vec[perm]
    return nvec
end

function tensorPermute(tensor)
    indx = tensor[1]
    indxN = permuteInd.(indx)
    return indxN , tensor[2]
end

function permuteIndN(vec,perm)
    #perm = [5,1,11,2,6,3,10,12,7,8,4,9]
    nvec = vec[perm]
    return nvec
end

function tensorPermuteN(tensor,perm)
    indx = tensor[1]
    indxN = []
    @assert length(indx[1]) == length(perm)
    for i in indx
        push!(indxN,i[perm])
    end
    return indxN , tensor[2]
end

export tensorProjectPrism

function tensorProjectPrism(tensor,K)
    indx = tensor[1]
    amps = tensor[2]
    indxP = []
    ampsP = []
    for i in 1:length(indx)
        if delta(indx[i][1],indx[i][2],indx[i][3],K) != 0
            if delta(indx[i][3],indx[i][4],indx[i][5],K) != 0
                if delta(indx[i][1],indx[i][5],indx[i][6],K) != 0 && delta(indx[i][2],indx[i][4],indx[i][6],K) != 0
                    if delta(indx[i][3],indx[i][7],indx[i][8],K) != 0
                        if delta(indx[i][1],indx[i][8],indx[i][9],K) != 0 && delta(indx[i][2],indx[i][7],indx[i][9],K) != 0
                            if delta(indx[i][9],indx[i][10],indx[i][11],K) != 0
                                if delta(indx[i][7],indx[i][10],indx[i][12],K) != 0 && delta(indx[i][2],indx[i][11],indx[i][12],K) != 0
                                    push!(indxP,indx[i])
                                    push!(ampsP,amps[i])
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return indxP , ampsP
end

end # module


text:
"""
$(SIGNATURES)

Generate a code replica for a signal from satellite system `S`. The
replica contains `num_samples` prompt samples as well as an additional
number of early and late samples specified by `correlator_sample_shifts`.
The codefrequency is specified by `code_frequency`, while the sampling
rate is given by `sampling_frequency`. The phase of the first prompt
sample is given by `start_code_phase`. The generated signal is returned
in the array `code_replica` with the first generated sample written to
index start_sample.
"""
function gen_code_replica!(
    code_replica,
    system::AbstractGNSS,
    code_frequency,
    sampling_frequency,
    start_code_phase::AbstractFloat,
    start_sample::Integer,
    num_samples::Integer,
    correlator_sample_shifts::AbstractVector,
    prn::Integer
)
    earliest_sample_shift = correlator_sample_shifts[end]
    latest_sample_shift  = correlator_sample_shifts[1]
    total_samples = num_samples + earliest_sample_shift - latest_sample_shift
    gen_code!(
        view(code_replica, start_sample:start_sample + total_samples - 1),
        system,
        prn,
        sampling_frequency,
        code_frequency,
        start_code_phase,
        latest_sample_shift
    )
    code_replica
end

"""
$(SIGNATURES)

Updates the code phase.
"""
function update_code_phase(
    system::AbstractGNSS,
    num_samples,
    code_frequency,
    sampling_frequency,
    start_code_phase,
    secondary_code_or_bit_found
)
    secondary_code_or_bit_length = get_data_frequency(system) == 0Hz ?
        get_secondary_code_length(system) :
        Int(get_code_frequency(system) / (get_data_frequency(system) * get_code_length(system)))

    code_length = get_code_length(system) *
        (secondary_code_or_bit_found ? secondary_code_or_bit_length : 1)
    mod(code_frequency * num_samples / sampling_frequency + start_code_phase, code_length)
#    fixed_point = sizeof(Int) * 8 - 1 - min_bits_for_code_length(S)
#    delta = floor(Int, code_frequency * 1 << fixed_point / sampling_frequency)
#    fixed_point_start_phase = floor(Int, start_code_phase * 1 << fixed_point)
#    phase_fixed_point = delta * num_samples + fixed_point_start_phase
#    mod(phase_fixed_point / 1 << fixed_point, code_length)
end

"""
$(SIGNATURES)

Calculates the current code frequency.
"""
function get_current_code_frequency(system::AbstractGNSS, code_doppler)
    code_doppler + get_code_frequency(system)
end


text:
using PyPlot
using DelimitedFiles
using PyCall

mpl = pyimport("tikzplotlib")
d = readdlm("timing.txt")
idx = sortperm(d[:,1])
d = d[idx,:]

close("all")
plot(d[:,1], 3.693 * ones(length(d[:,1])), "--", label="Fortran")
loglog(d[:,1], d[:,2], "o-", label="ADSeismic MPI")
legend()
xlabel("Number of Processors")
ylabel("Time (sec)")
xticks(d[:,1], Int.(d[:,1]))
grid("on", which="both", linestyle=":")
savefig("acoustic_time_forward.png")

close("all")
loglog(d[:,1], d[:,3], "o-", color ="orange", label="ADSeismic MPI")
legend()
xlabel("Number of Processors")
ylabel("Time (sec)")
xticks(d[:,1], Int.(d[:,1]))
grid("on", which="both", linestyle=":")
savefig("acoustic_time_backward.png")

close("all")
loglog(d[:,1], d[:,2], "o-", label="Forward")
loglog(d[:,1], d[:,3], "o-", color ="orange", label="Backward")
legend()
xlabel("Number of Processors")
ylabel("Time (sec)")
xticks(d[:,1], Int.(d[:,1]))
grid("on", which="both", linestyle=":")
savefig("acoustic_time_forward_and_backward.png")
mpl.save("../figures/acoustic_time_forward_and_backward.tex")
close("all")
figure(figsize=(10,4))
subplot(121)
title("Forward")
loglog(d[:,1], d[1,2]./d[:,2], "o-", label="Speedup")
loglog(d[:,1], d[1,2]./(d[:,2].*d[:,1]), "o-", label="Efficiency")
legend()
xlabel("Number of Processors")
ylabel("Time (sec)")
xticks(d[:,1], Int.(d[:,1]))
grid("on", which="both", linestyle=":")

subplot(122)
title("Backward")
loglog(d[:,1], d[1,3]./d[:,3], "o-", label="Speedup")
loglog(d[:,1], d[1,3]./(d[:,3].*d[:,1]), "o-", label="Efficiency")
legend()
xlabel("Number of Processors")
ylabel("Time (sec)")
xticks(d[:,1], Int.(d[:,1]))
grid("on", which="both", linestyle=":")

savefig("acoustic_speedup_and_efficiency.png")
mpl.save("../figures/acoustic_speedup_and_efficiency.tex")


text:
function build_robot!(vis::Visualizer, model::PushBot; name::Symbol=:PushBot,
		r=0.05, n_leg::Int=30, α=1.0)
	r = convert(Float32, r)
	r_contact = convert(Float32, r * 1.5)
	r_arm = convert(Float32, r * 0.3)
	wall_thickness = 0.15
	l_joint = 0.16
	body_mat = MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, α))
	contact_mat = MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, α))
	wall_mat = MeshPhongMaterial(color = RGBA(0.7, 0.7, 0.7, 1.0))

	default_background!(vis)

	link = Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, model.l-3r_arm), r)
	setobject!(vis[name][:robot]["link"], link, body_mat)
	setobject!(vis[name][:robot]["linkbase"], Sphere(Point3f0(0.0), r), body_mat)

	setobject!(vis[name][:robot]["contact1"], Sphere(Point3f0(0.0), r_contact), contact_mat)
	setobject!(vis[name][:robot]["joint1"],
		Rect(Vec(-l_joint/2, -r/3, -4r_arm),Vec(l_joint, 2r/3, 2r_arm)), body_mat)
	setobject!(vis[name][:robot]["joint2"],
		Rect(Vec(-l_joint/2, -r/3, 2r_arm),Vec(l_joint, 2r/3, 2r_arm)), body_mat)

	for i = 1:n_leg
		setobject!(vis[name][:robot]["arm$i"], Sphere(Point3f0(0), r_arm), body_mat)
	end

	setobject!(vis[name][:env]["wall1"],
		Rect(Vec(0, 0, 0),Vec(wall_thickness, wall_thickness, 1.5)), wall_mat)

	setobject!(vis[name][:env]["wall2"],
		Rect(Vec(0, 0, 0),Vec(wall_thickness, wall_thickness, 1.5)), wall_mat)

	settransform!(vis[name][:env]["wall1"], Translation([0.5 + 1.5 * r;                 -wall_thickness/2; 0.0]))
	settransform!(vis[name][:env]["wall2"], Translation([-wall_thickness-0.5 - 1.5 * r; -wall_thickness/2; 0.0]))

	# settransform!(vis["/Cameras/default"],
	# 	compose(Translation(0.0, 0.5, -1.0), LinearMap(RotZ(-pi / 2.0))))
	return nothing
end

function set_robot!(vis::Visualizer, model::PushBot, q::AbstractVector;
		name::Symbol=:PushBot, n_leg::Int=30)

	pee = [_kinematics(model, q, mode = :ee)[1], 0.0, _kinematics(model, q, mode = :ee)[2]]
    p1 = [_kinematics(model, q, mode = :d)[1], 0.0, _kinematics(model, q, mode = :d)[2]]

	settransform!(vis[name][:robot]["joint1"], compose(cable_transform(zeros(3), pee), Translation(0,0,0)))
	settransform!(vis[name][:robot]["joint2"], cable_transform(zeros(3), pee))
	settransform!(vis[name][:robot]["link"], cable_transform(pee, zeros(3)))
	settransform!(vis[name][:robot]["contact1"], Translation(p1))

	r_range = range(0, stop = q[2], length = n_leg)
	for i = 1:n_leg
		q_tmp = Array(copy(q))
		q_tmp[2] = r_range[i]
		p_leg = [_kinematics(model, q_tmp, mode = :d)[1], 0.0, _kinematics(model, q_tmp, mode = :d)[2]]
        settransform!(vis[name][:robot]["arm$i"], Translation(p_leg))
    end

	return nothing
end

function pusher_arc(center::AbstractVector, radius::Real, α::Real; sense::Symbol=:trig)
	p = center + radius*[cos(α), 0, sin(α)]
	if sense == :trig
		θ = α + π/2
	elseif sense == :antitrig
		θ = α - π/2
	else
		error("Invalid sense.")
	end
	pθ = [p; θ]
	return pθ
end

function pusher_arc_traj(center::AbstractVector, radius::Real, α_start::Real, α_end::Real,
	 	N::Int; sense::Symbol=:trig)
	pθ = []
	for α in range(α_start, stop=α_end, length=N)
		pθ_ = pusher_arc(center, radius, α, sense=sense)
		push!(pθ, pθ_)
	end
	return pθ
end

function generate_pusher_traj(d::ImpulseDisturbance, traj::ContactTraj; side::Symbol=:right,
		center=zeros(3), radius::Real=0.8, N_arc::Int=17)
	H = traj.H
	if side == :right
		sense = :trig
		θ_rest = 0.0
	elseif side == :left
		sense = :antitrig
		θ_rest = π
	else
		error("Invalid side.")
	end

	pθ = [pusher_arc(center, radius, θ_rest; sense=sense) for t=1:H_sim]

	for (i,idx) in enumerate(d.idx)
		w = d.w[i][1]
		θ_impact = traj.q[idx+1][1] + π/2
		θ_next = traj.q[idx+2][1] + π/2

		if ((w > 0.0) && (side == :right)) || ((w < 0.0) && (side == :left))
			forward_arc = pusher_arc_traj(center, radius, θ_rest, θ_impact, N_arc; sense=sense)
			backward_arc = pusher_arc_traj(center, radius, θ_next, θ_rest, N_arc; sense=sense)
			arc = [forward_arc; backward_arc] # 2N
			pθ[idx-N_arc-0:idx+N_arc-1] .= arc
		end
	end
	return pθ
end

function build_disturbance!(vis::Visualizer, model::PushBot; name::Symbol=:Pusher, r=0.025, α=1.0)
	r1 = convert(Float32, r)
	r2 = convert(Float32, 4r)
	pusher_mat = MeshPhongMaterial(color = RGBA(1.0, 0.0, 0.0, α))

	# build pusher
	p1 = Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, 0.250), r1)
	p2 = Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, 0.075), r2)
	setobject!(vis[name, "pusher1"], p1, pusher_mat)
	setobject!(vis[name, "pusher2"], p2, pusher_mat)
	return anim
end

function set_disturbance!(vis::Visualizer, model::PushBot, pθ::AbstractVector; name::Symbol=:Pusher, offset::Real=0.0)
	p = pθ[1:3]
	θ = pθ[4]
	settransform!(vis[name],
		compose(Translation(p...),
		compose(LinearMap(RotY(-θ-pi/2)),
				Translation(0,0,offset))))
	return nothing
end

function animate_disturbance!(vis::Visualizer, anim::MeshCat.Animation, model::PushBot,
		pθ::AbstractVector; name::Symbol=:Pusher, offset::Real=0.0)
	H = length(pθ)
	for t in 1:H
		MeshCat.atframe(anim, t) do
			set_disturbance!(vis, model, pθ[t], name=name, offset=offset)
		end
	end
	setanimation!(vis, anim)
	return anim
end

function visualize_disturbance!(vis::Visualizer, model::ContactModel, pθ::AbstractVector;
		h=0.01, α=1.0,
		sample=max(1, Int(floor(length(pθ) / 100))),
		anim::MeshCat.Animation=MeshCat.Animation(Int(floor(1/h))),
		name::Symbol=:Pusher, offset::Real=0.0)

	build_disturbance!(vis, model, name=name, α=α)
	animate_disturbance!(vis, anim, model, pθ[1:sample:end], name=name, offset=offset)
	return anim
end


text:
function joinvec(embed::WordEmbedding, words::Vector, fn::Function, wv::Array=[])
    for word in words
        if !haskey(embed.embedding, word)
            warn("'$word' not present in the vocabulary")
            return wv
        end
        if isempty(wv)
            wv = vec(embed.embedding[word])
        else
            wv = fn(wv, vec(embed.embedding[word]))
        end
    end
    wv
end

function find_nearest_words(embed::WordEmbedding, words::AbstractString; k=5)
    positive_words = AbstractString[]
    negative_words = AbstractString[]
    wordlist = positive_words
    for tok in split(words)
        tok = strip(tok)
        isempty(tok) && continue
        if tok == "+"
            wordlist = positive_words
        elseif tok == "-"
            wordlist = negative_words
        else
            push!(wordlist, tok)
        end
    end
    find_nearest_words(embed, positive_words, negative_words; k=k)
end

function find_nearest_words(embed::WordEmbedding, positive_words::Vector, negative_words::Vector; k=5)
    pq = PriorityQueue(Base.Order.Reverse)

    wv = joinvec(embed, positive_words, .+, Array[])
    wv = joinvec(embed, negative_words, .-, wv)

    if !isempty(wv)
        for (w, embed_w) in embed.embedding
            ((w in positive_words) || (w in negative_words)) && continue
            dist = cosine_dist(wv, vec(embed_w))
            enqueue!(pq, w, dist)
            (length(pq) > k) && dequeue!(pq)
        end
    end
    sort(collect(pq), by = t -> t[2])
end


text:
#####
##### Time-stepping functionality that is independent of the TimeStepper
#####

"""
    precomputations!(diffusivities, pressures, velocities, tracers, model)

Perform precomputations necessary for an explicit timestep or substep.
"""
function precomputations!(diffusivities, pressures, velocities, tracers, model)

    # Fill halos for velocities and tracers
    fill_halo_regions!(merge(model.velocities, model.tracers), model.architecture,
                       model.clock, state(model))

    # Calculate diffusivities
    calculate_diffusivities!(diffusivities, model.architecture, model.grid, model.closure,
                             model.buoyancy, velocities, tracers)

    fill_halo_regions!(model.diffusivities, model.architecture, model.clock, state(model))

    # Calculate hydrostatic pressure
    pressure_calculation = launch!(model.architecture, model.grid, :xy, update_hydrostatic_pressure!,
                                   pressures.pHY′, model.grid, model.buoyancy, tracers,
                                   dependencies=Event(device(model.architecture)))

    # Fill halo regions for pressure
    wait(device(model.architecture), pressure_calculation)

    fill_halo_regions!(model.pressures.pHY′, model.architecture)

    return nothing
end

"""
    calculate_tendencies!(diffusivities, pressures, velocities, tracers, model)

Calculate the interior and boundary contributions to tendency terms without the
contribution from non-hydrostatic pressure.
"""
function calculate_tendencies!(tendencies, velocities, tracers, pressures, diffusivities, model)

    # Note:
    #
    # "tendencies" is a NamedTuple of OffsetArrays corresponding to the tendency data for use
    # in GPU computations.
    #
    # "model.timestepper.Gⁿ" is a NamedTuple of Fields, whose data also corresponds to
    # tendency data.

    # Calculate contributions to momentum and tracer tendencies from fluxes and volume terms in the
    # interior of the domain
    calculate_interior_tendency_contributions!(tendencies, model.architecture, model.grid, model.advection,
                                               model.coriolis, model.buoyancy, model.surface_waves,
                                               model.closure, velocities, tracers, pressures.pHY′,
                                               diffusivities, model.forcing, model.clock)

    # Calculate contributions to momentum and tracer tendencies from user-prescribed fluxes across the
    # boundaries of the domain
    calculate_boundary_tendency_contributions!(model.timestepper.Gⁿ, model.architecture, model.velocities,
                                               model.tracers, model.clock, state(model))

    return nothing
end

"""
    calculate_pressure_correction!(nonhydrostatic_pressure, Δt, tendencies, velocities, model)

Calculate the (nonhydrostatic) pressure correction associated `tendencies`, `velocities`, and step size `Δt`.
"""
function calculate_pressure_correction!(nonhydrostatic_pressure, Δt, velocities, model)

    fill_halo_regions!(model.velocities, model.architecture, model.clock, state(model))

    solve_for_pressure!(nonhydrostatic_pressure, model.pressure_solver, model.architecture, model.grid, Δt, velocities)

    fill_halo_regions!(model.pressures.pNHS, model.architecture)

    return nothing
end

#####
##### Fractional and time stepping
#####

"""
Update the predictor velocities u, v, and w with the non-hydrostatic pressure via

    `u^{n+1} = u^n - δₓp_{NH} / Δx * Δt`
"""
@kernel function _pressure_correct_velocities!(U, grid, Δt, pNHS)
    i, j, k = @index(Global, NTuple)

    @inbounds U.u[i, j, k] -= ∂xᶠᵃᵃ(i, j, k, grid, pNHS) * Δt
    @inbounds U.v[i, j, k] -= ∂yᵃᶠᵃ(i, j, k, grid, pNHS) * Δt
    @inbounds U.w[i, j, k] -= ∂zᵃᵃᶠ(i, j, k, grid, pNHS) * Δt
end

"Update the solution variables (velocities and tracers)."
function pressure_correct_velocities!(U, arch, grid, Δt, pNHS)
    event = launch!(arch, grid, :xyz, _pressure_correct_velocities!, U, grid, Δt, pNHS,
                    dependencies=Event(device(arch)))
    wait(device(arch), event)
    return nothing
end


text:
# Sequence

for _Seq in [:HSeq, :VSeq]
  @eval begin
    # constructor
    @test length($_Seq()) == 0
    @test length($_Seq(Event())) == 1
    @test length($_Seq(HSeq())) == 1
    @test length($_Seq(VSeq())) == 1
    @test length($_Seq(Event(), HSeq(), VSeq())) == 3
    @test length($_Seq([Event(), HSeq(), VSeq()])) == 3

    # conversion
    @test convert($_Seq, Event(:foo => 1)) == $_Seq(Event(:foo => 1))
    @test convert($_Seq, [Event(:foo => 1)]) == $_Seq(Event(:foo => 1))

    # show
    @test length("$($_Seq())") > 0
    @test length("$($_Seq(Event()))") > 0

    # merge
    @test merge($_Seq(), $_Seq()) == $_Seq()
    @test merge($_Seq(Event(:foo => 1)), $_Seq(Event(:bar => 1))) ==
      $_Seq(Event(:foo => 1, :bar => 1))
    @test merge(Event(), $_Seq()) == $_Seq()
    @test merge(Event(), $_Seq(Event())) == $_Seq(Event())
    @test merge($_Seq(), Event()) == $_Seq()
    @test merge($_Seq(Event()), Event()) == $_Seq(Event())
    @test merge($_Seq(Event()), [Dict(:foo => 1)]) == $_Seq(Event(:foo => 1))
    @test merge($_Seq(Event()), Dict(:foo => 1)) == $_Seq(Event(:foo => 1))
    @test merge($_Seq(Event()), :foo => 1) == $_Seq(Event(:foo => 1))

    # repeat
    @test repeat($_Seq(Event(:foo => 1), Event(:bar => 1)), inner=2) ==
      $_Seq(Event(:foo => 1), Event(:foo => 1), Event(:bar => 1), Event(:bar => 1))
    @test repeat($_Seq(Event(:foo => 1), Event(:bar => 1)), outer=2) ==
      $_Seq(Event(:foo => 1), Event(:bar => 1), Event(:foo => 1), Event(:bar => 1))

    # collection
    @test collect($_Seq, []) == $_Seq()
    @test collect($_Seq, [Event(), VSeq(), HSeq()]) ==
      $_Seq(Event(), VSeq(), HSeq())

    @test circshift($_Seq(Event(:i => 1), Event(:i => 2), Event(:i => 3)), 1) ==
      $_Seq(Event(:i => 3), Event(:i => 1), Event(:i => 2))
    @test circshift($_Seq(Event(:i => 1), Event(:i => 2), Event(:i => 3)), -1) ==
      $_Seq(Event(:i => 2), Event(:i => 3), Event(:i => 1))

    @test reverse($_Seq(Event(:i => 1), Event(:i => 2), Event(:i => 3))) ==
      $_Seq(Event(:i => 3), Event(:i => 2), Event(:i => 1))

    @test rot180($_Seq(Event(:i => 1), Event(:i => 2), Event(:i => 3))) ==
      $_Seq(Event(:i => 3), Event(:i => 2), Event(:i => 1))
  end
end

# specific conversion
@test convert(HSeq, VSeq(Event())) == HSeq(Event())
@test convert(VSeq, HSeq(Event())) == VSeq(Event())

# matrix
@test HSeq([Event(:i => 0) Event(:i => 1); Event(:i => 2) Event(:i => 3)]) ==
  HSeq(VSeq(Event(:i => 0),Event(:i => 1)), VSeq(Event(:i => 2),Event(:i => 3)))
@test VSeq([Event(:i => 0) Event(:i => 1); Event(:i => 2) Event(:i => 3)]) ==
  VSeq(HSeq(Event(:i => 0),Event(:i => 1)), HSeq(Event(:i => 2),Event(:i => 3)))

m =
  VSeq([Event(:i => 0) Event(:i => 1) Event(:i => 2);
        Event(:i => 3) Event(:i => 4) Event(:i => 5)])

r1 = rotr90(m)
r2 = rotr90(r1)
r3 = rotr90(r2)
r4 = rotr90(r3)
l1 = rotl90(m)
l2 = rotl90(l1)
l3 = rotl90(l2)
l4 = rotl90(l3)

@test r1 == l3 ==
  HSeq([Event(:i => 3) Event(:i => 4) Event(:i => 5);
        Event(:i => 0) Event(:i => 1) Event(:i => 2)])

@test r2 == l2 ==
  VSeq([Event(:i => 5) Event(:i => 4) Event(:i => 3);
        Event(:i => 2) Event(:i => 1) Event(:i => 0)])

@test r3 == l1 ==
  HSeq([Event(:i => 2) Event(:i => 1) Event(:i => 0);
        Event(:i => 5) Event(:i => 4) Event(:i => 3)])

@test r4 == l4 == m


text:
using LinearAlgebra
using LightGraphs
using SimpleWeightedGraphs
using GraphHierarchy
using SparseArrays
using Test

# (1,1)-(2,2)

@testset "GraphHierarchy.jl" begin
    # 3-chain
    gw = SimpleWeightedDiGraph(3)
    add_edge!(gw,1,2,2)
    add_edge!(gw,2,3,2)
    # forward: not weighted, not big
    HS = forward_hierarchical_structure(gw, weighted = false)
    @test norm(HS[1] - [-1,0,1]) < 1.0e-12
    @test norm(HS[2] - [1,0,0]) < 1.0e-12
    @test norm(HS[3] - [1,1,0]) < 1.0e-12
    @test norm(HS[4] - [ 0 1 0 ; 0 0 1 ; 0 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1]) < 1.0e-12
    @test abs(HC[2]) < 1.0e-12
    # forward: weighted, not big
    HS = forward_hierarchical_structure(gw)
    @test norm(HS[1] - [-1,0,1]) < 1.0e-12
    @test norm(HS[2] - [1,0,0]) < 1.0e-12
    @test norm(HS[3] - [1,1,0]) < 1.0e-12
    @test norm(HS[4] - [ 0 1 0 ; 0 0 1 ; 0 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1]) < 1.0e-12
    @test abs(HC[2]) < 1.0e-12
    # forward: not weighted, big
    HS = forward_hierarchical_structure(gw, weighted = false, big = true)
    @test norm(HS[1] - [-1,0,1]) < 1.0e-70
    @test norm(HS[2] - [1,0,0]) < 1.0e-70
    @test norm(HS[3] - [1,1,0]) < 1.0e-70
    @test norm(HS[4] - [ 0 1 0 ; 0 0 1 ; 0 0 0 ]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1]) < 1.0e-70
    @test abs(HC[2]) < 1.0e-70
    # forward: weighted, big
    HS = forward_hierarchical_structure(gw, big = true)
    @test norm(HS[1] - [-1,0,1]) < 1.0e-70
    @test norm(HS[2] - [1,0,0]) < 1.0e-70
    @test norm(HS[3] - [1,1,0]) < 1.0e-70
    @test norm(HS[4] - [ 0 1 0 ; 0 0 1 ; 0 0 0 ]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1]) < 1.0e-70
    @test abs(HC[2]) < 1.0e-70
    # backward: not weighted, not big
    HS = backward_hierarchical_structure(gw, weighted = false)
    @test norm(HS[1] - [1,0,-1]) < 1.0e-12
    @test norm(HS[2] - [0,0,1]) < 1.0e-12
    @test norm(HS[3] - [0,1,1]) < 1.0e-12
    @test norm(HS[4] - [ 0 1 0 ; 0 0 1 ; 0 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1]) < 1.0e-12
    @test abs(HC[2]) < 1.0e-12
    # backward: weighted, not big
    HS = backward_hierarchical_structure(gw)
    @test norm(HS[1] - [1,0,-1]) < 1.0e-12
    @test norm(HS[2] - [0,0,1]) < 1.0e-12
    @test norm(HS[3] - [0,1,1]) < 1.0e-12
    @test norm(HS[4] - [ 0 1 0 ; 0 0 1 ; 0 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1]) < 1.0e-12
    @test abs(HC[2]) < 1.0e-12
    # backward: not weighted, big
    HS = backward_hierarchical_structure(gw, weighted = false, big = true)
    @test norm(HS[1] - [1,0,-1]) < 1.0e-70
    @test norm(HS[2] - [0,0,1]) < 1.0e-70
    @test norm(HS[3] - [0,1,1]) < 1.0e-70
    @test norm(HS[4] - [ 0 1 0 ; 0 0 1 ; 0 0 0 ]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1]) < 1.0e-70
    @test abs(HC[2]) < 1.0e-70
    # backward: weighted, big
    HS = backward_hierarchical_structure(gw, big = true)
    @test norm(HS[1] - [1,0,-1]) < 1.0e-70
    @test norm(HS[2] - [0,0,1]) < 1.0e-70
    @test norm(HS[3] - [0,1,1]) < 1.0e-70
    @test norm(HS[4] - [ 0 1 0 ; 0 0 1 ; 0 0 0 ]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1]) < 1.0e-70
    @test abs(HC[2]) < 1.0e-70

    # 3-cycle
    gw = SimpleWeightedDiGraph(3)
    add_edge!(gw,1,2,1)
    add_edge!(gw,2,3,1)
    add_edge!(gw,3,1,2)
    # forward: not weighted, not big
    HS = forward_hierarchical_structure(gw, weighted = false)
    @test norm(HS[1]) < 1.0e-12
    @test norm(HS[2] - [1,1,1]) < 1.0e-12
    @test norm(HS[3]) < 1.0e-12
    @test norm(HS[4]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1] - 1) < 1.0e-12
    @test abs(HC[2]) < 1.0e-12
    # forward: weighted, not big
    HS = forward_hierarchical_structure(gw)
    @test norm(HS[1] - [1//3,0,-1//3]) < 1.0e-12
    @test norm(HS[2] - [1//3, 4//3, 4//3]) < 1.0e-12
    @test norm(HS[3] - [-1//3,-1//3,2//3]) < 1.0e-12
    @test norm(HS[4] - [ 0 -1//3 0 ; 0 0 -1//3 ; 2//3 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1] - 5//6) < 1.0e-12
    @test abs(HC[2] - 1//2) < 1.0e-12
    # forward: not weighted, big
    HS = forward_hierarchical_structure(gw, weighted = false, big = true)
    @test norm(HS[1]) < 1.0e-70
    @test norm(HS[2] - [1,1,1]) < 1.0e-70
    @test norm(HS[3]) < 1.0e-70
    @test norm(HS[4]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1] - 1) < 1.0e-70
    @test abs(HC[2]) < 1.0e-70
    # forward: weighted, big
    HS = forward_hierarchical_structure(gw, big = true)
    @test norm(HS[1] - [1//3,0,-1//3]) < 1.0e-70
    @test norm(HS[2] - [1//3, 4//3, 4//3]) < 1.0e-70
    @test norm(HS[3] - [-1//3,-1//3,2//3]) < 1.0e-70
    @test norm(HS[4] - [ 0 -1//3 0 ; 0 0 -1//3 ; 2//3 0 0 ]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1] - 5//6) < 1.0e-70
    @test abs(HC[2] - 1//2) < 1.0e-70
    # backward: not weighted, not big
    HS = backward_hierarchical_structure(gw, weighted = false)
    @test norm(HS[1]) < 1.0e-12
    @test norm(HS[2] - [1,1,1]) < 1.0e-12
    @test norm(HS[3]) < 1.0e-12
    @test norm(HS[4]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1] - 1) < 1.0e-12
    @test abs(HC[2]) < 1.0e-12
    # # backward: weighted, not big
    HS = backward_hierarchical_structure(gw)
    @test norm(HS[1] - [-1//3,0,1//3]) < 1.0e-12
    @test norm(HS[2] - [4//3,4//3,1//3]) < 1.0e-12
    @test norm(HS[3] - [2//3,-1//3,-1//3]) < 1.0e-12
    @test norm(HS[4] - [ 0 -1//3 0 ; 0 0 -1//3 ; 2//3 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1] - 5//6) < 1.0e-12
    @test abs(HC[2] - 1//2) < 1.0e-12
    # backward: not weighted, big
    HS = backward_hierarchical_structure(gw, weighted = false, big = true)
    @test norm(HS[1]) < 1.0e-70
    @test norm(HS[2] - [1,1,1]) < 1.0e-70
    @test norm(HS[3]) < 1.0e-70
    @test norm(HS[4]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1] - 1) < 1.0e-70
    @test abs(HC[2]) < 1.0e-70
    # backward: weighted, big
    HS = backward_hierarchical_structure(gw, big = true)
    @test norm(HS[1] - [-1//3,0,1//3]) < 1.0e-70
    @test norm(HS[2] - [4//3,4//3,1//3]) < 1.0e-70
    @test norm(HS[3] - [2//3,-1//3,-1//3]) < 1.0e-70
    @test norm(HS[4] - [ 0 -1//3 0 ; 0 0 -1//3 ; 2//3 0 0 ]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1] - 5//6) < 1.0e-70
    @test abs(HC[2] - 1//2) < 1.0e-70

    # pair with stem
    gw = SimpleWeightedDiGraph(3)
    add_edge!(gw,1,2,2)
    add_edge!(gw,2,3,1)
    add_edge!(gw,2,1,1)
    # forward: not weighted, not big
    HS = forward_hierarchical_structure(gw, weighted = false)
    @test norm(HS[1] - [-1//3, -1//3, 2//3]) < 1.0e-12
    @test norm(HS[2] - [1,1,0]) < 1.0e-12
    @test norm(HS[3] - [0,1//2,0]) < 1.0e-12
    @test norm(HS[4] - [ 0 0 0 ; 0 0 1 ; 0 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1] - 2//3) < 1.0e-12
    @test abs(HC[2] - sqrt(2)/3) < 1.0e-12
    # forward: weighted, not big
    HS = forward_hierarchical_structure(gw)
    @test norm(HS[1] - [-11//15,-2//15,13//15]) < 1.0e-12
    @test norm(HS[2] - [8//5,2//5,0]) < 1.0e-12
    @test norm(HS[3] - [3//5,1//5,0]) < 1.0e-12
    @test norm(HS[4] - [ 0 3//5 0 ; -3//5 0 1 ; 0 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1] - 3//5) < 1.0e-12
    @test abs(HC[2] - 3//5) < 1.0e-12
    # forward: not weighted, big
    HS = forward_hierarchical_structure(gw, weighted = false, big = true)
    @test norm(HS[1] - [-1//3, -1//3, 2//3]) < 1.0e-70
    @test norm(HS[2] - [1,1,0]) < 1.0e-70
    @test norm(HS[3] - [0,1//2,0]) < 1.0e-70
    @test norm(HS[4] - [ 0 0 0 ; 0 0 1 ; 0 0 0 ]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1] - 2//3) < 1.0e-70
    @test abs(HC[2] - sqrt(big(2))/3) < 1.0e-70
    # forward: weighted, big
    HS = forward_hierarchical_structure(gw, big = true)
    @test norm(HS[1] - [-11//15,-2//15,13//15]) < 1.0e-70
    @test norm(HS[2] - [8//5,2//5,0]) < 1.0e-70
    @test norm(HS[3] - [3//5,1//5,0]) < 1.0e-70
    @test norm(HS[4] - [ 0 3//5 0 ; -3//5 0 1 ; 0 0 0 ]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1] - 3//5) < 1.0e-70
    @test abs(HC[2] - 3//5) < 1.0e-70
    # backward: not weighted, not big
    HS = backward_hierarchical_structure(gw, weighted = false)
    @test norm(HS[1] - [5//3,2//3,-7//3]) < 1.0e-12
    @test norm(HS[2] - [0,0,1]) < 1.0e-12
    @test norm(HS[3] - [-1,1,3]) < 1.0e-12
    @test norm(HS[4] - [ 0 1 0 ; -1 0 3 ; 0 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1]) < 1.0e-12
    @test abs(HC[2] - 2*sqrt(2/3)) < 1.0e-12
    # backward: weighted, not big
    HS = backward_hierarchical_structure(gw)
    @test norm(HS[1] - [5//3,2//3,-7//3]) < 1.0e-12
    @test norm(HS[2] - [0,0,1]) < 1.0e-12
    @test norm(HS[3] - [-1,1,3]) < 1.0e-12
    @test norm(HS[4] - [ 0 1 0 ; -1 0 3 ; 0 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1]) < 1.0e-12
    @test abs(HC[2] - sqrt(2)) < 1.0e-12
    # backward: not weighted, big
    HS = backward_hierarchical_structure(gw, weighted = false, big = true)
    @test norm(HS[1] - [5//3,2//3,-7//3]) < 1.0e-70
    @test norm(HS[2] - [0,0,1]) < 1.0e-70
    @test norm(HS[3] - [-1,1,3]) < 1.0e-70
    @test norm(HS[4] - [ 0 1 0 ; -1 0 3 ; 0 0 0 ]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1]) < 1.0e-70
    @test abs(HC[2] - 2*sqrt(big(2)/3)) < 1.0e-70
    # backward: weighted, big
    HS = backward_hierarchical_structure(gw, big = true)
    @test norm(HS[1] - [5//3,2//3,-7//3]) < 1.0e-70
    @test norm(HS[2] - [0,0,1]) < 1.0e-70
    @test norm(HS[3] - [-1,1,3]) < 1.0e-70
    @test norm(HS[4] - [ 0 1 0 ; -1 0 3 ; 0 0 0 ]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1]) < 1.0e-70
    @test abs(HC[2] - sqrt(big(2))) < 1.0e-70

    # unbalanced triangle
    gw = SimpleWeightedDiGraph(3)
    add_edge!(gw,1,2,0.1)
    add_edge!(gw,2,3,1)
    add_edge!(gw,3,1,1)
    add_edge!(gw,2,1,1)
    # forward: not weighted, not big
    HS = forward_hierarchical_structure(gw, weighted = false)
    @test norm(HS[1] - [7//18,-5//18,-1//9]) < 1.0e-12
    @test norm(HS[2] - [5//12,5//3,5//6]) < 1.0e-12
    @test norm(HS[3] - [-2//3,5//12,1//2]) < 1.0e-12
    @test norm(HS[4] - [ 0 -2//3 0 ; 2//3 0 1//6 ; 1//2 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1] - 5//6) < 1.0e-12
    @test abs(HC[2] - sqrt(19/2)/6) < 1.0e-12
    # forward: weighted, not big
    HS = forward_hierarchical_structure(gw)
    @test norm(HS[1] - [799//1206,-995//1206,98//603]) < 1.0e-12
    @test norm(HS[2] - [5//804,500//201,5//402]) < 1.0e-12
    @test norm(HS[3] - [-299//201,995//804,1//2]) < 1.0e-12
    @test norm(HS[4] - [ 0 -299//201 0 ; 299//201 0 397//402 ; 1//2 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1] - 550//6231) < 1.0e-12
    @test abs(HC[2] - sqrt(13557755)/6231) < 1.0e-12
    # forward: not weighted, big
    HS = forward_hierarchical_structure(gw, weighted = false, big = true)
    @test norm(HS[1] - [7//18,-5//18,-1//9]) < 1.0e-70
    @test norm(HS[2] - [5//12,5//3,5//6]) < 1.0e-70
    @test norm(HS[3] - [-2//3,5//12,1//2]) < 1.0e-70
    @test norm(HS[4] - [ 0 -2//3 0 ; 2//3 0 1//6 ; 1//2 0 0 ]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1] - 5//6) < 1.0e-70
    @test abs(HC[2] - sqrt(big(19)/2)/6) < 1.0e-70
    # forward: weighted, big  ---  Note that since the weights are given as Float64, accuracy is low even when using BigFloat
    HS = forward_hierarchical_structure(gw, big = true)
    @test norm(HS[1] - [799//1206,-995//1206,98//603]) < 1.0e-16
    @test norm(HS[2] - [5//804,500//201,5//402]) < 1.0e-16
    @test norm(HS[3] - [-299//201,995//804,1//2]) < 1.0e-16
    @test norm(HS[4] - [ 0 -299//201 0 ; 299//201 0 397//402 ; 1//2 0 0 ]) < 1.0e-16
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1] - 550//6231) < 1.0e-16
    @test abs(HC[2] - sqrt(13557755)/6231) < 1.0e-16
    # backward: not weighted, not big
    HS = backward_hierarchical_structure(gw, weighted = false)
    @test norm(HS[1] - [-5//18,7//18,-1//9]) < 1.0e-12
    @test norm(HS[2] - [5//3,5//12,5//6]) < 1.0e-12
    @test norm(HS[3] - [5//12,-2//3,1//2]) < 1.0e-12
    @test norm(HS[4] - [ 0 -2//3 0 ; 2//3 0 1//2 ; 1//6 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1] - 5//6) < 1.0e-12
    @test abs(HC[2] - sqrt(19/2)/6) < 1.0e-12
    # # backward: weighted, not big
    HS = backward_hierarchical_structure(gw)
    @test norm(HS[1] - [-995//1206,799//1206,98//603]) < 1.0e-12
    @test norm(HS[2] - [500//201,5//804,5//402]) < 1.0e-12
    @test norm(HS[3] - [995//804,-299//201,1//2]) < 1.0e-12
    @test norm(HS[4] - [ 0 -299//201 0 ; 299//201 0 1//2 ; 397//402 0 0 ]) < 1.0e-12
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1] - 550//6231) < 1.0e-12
    @test abs(HC[2] - sqrt(13557755)/6231) < 1.0e-12
    # backward: not weighted, big
    HS = backward_hierarchical_structure(gw, weighted = false, big = true)
    @test norm(HS[1] - [-5//18,7//18,-1//9]) < 1.0e-70
    @test norm(HS[2] - [5//3,5//12,5//6]) < 1.0e-70
    @test norm(HS[3] - [5//12,-2//3,1//2]) < 1.0e-70
    @test norm(HS[4] - [ 0 -2//3 0 ; 2//3 0 1//2 ; 1//6 0 0 ]) < 1.0e-70
    HC = hierarchical_coefficients(gw, HS, weighted = false)
    @test abs(HC[1] - 5//6) < 1.0e-70
    @test abs(HC[2] - sqrt(big(19)/2)/6) < 1.0e-70
    # backward: weighted, big
    HS = backward_hierarchical_structure(gw, big = true)
    @test norm(HS[1] - [-995//1206,799//1206,98//603]) < 1.0e-16
    @test norm(HS[2] - [500//201,5//804,5//402]) < 1.0e-16
    @test norm(HS[3] - [995//804,-299//201,1//2]) < 1.0e-16
    @test norm(HS[4] - [ 0 -299//201 0 ; 299//201 0 1//2 ; 397//402 0 0  ]) < 1.0e-16
    HC = hierarchical_coefficients(gw, HS)
    @test abs(HC[1] - 550//6231) < 1.0e-16
    @test abs(HC[2] - sqrt(big(13557755))/6231) < 1.0e-16
end


text:
function read_and_prep_osm_map(;map_path::String
								, routes_dict::OrderedDict{Symbol, OrderedDict{Int,Tuple{Float64,Float64}}}
								, secs_between_stops::Dict{Tuple{Int,Int},Int64}
								, TTC_car_freq::Int = sim_params.TTC_car_freq
								, agents_speed_in_sec_per_m::Float64 = sim_params.agents_speed_in_sec_per_m
								, parset_id::Int64
								, serialize_results::Bool = false)
	"""
	Description:
		Reads the specified .osm map (e.g. central_toronto) with all related objects as nodes, graph etc.
		Next it will add TTC information (routes and stops as egdes and vertices) to the map.
		It will also generate structures which will be used in the simulation for routing (e.g. routes_path, routes_distances).
		Map, routes_path and routes_distances will be serialized and saved to the data_output folder.

		NOTE! Create data_output folder in your working directory.
	Returns:
		Map, number of nodes in the map before TTC info adding, routes_path, routes_distances
	"""

	println("By default TTC_car_freq will be used from object `sim_params`")
	println("TTC car freq specified: ", TTC_car_freq)
    m = OpenStreetMapX.get_map_data(map_path, use_cache=false, trim_to_connected_graph=true);
    orig_map_nodes_num = nv(m.g)
    node_id_stop_id_dict, stop_id_to_node_id_dict = adding_TTC_points_to_map!(routes_dict, m);
    routes_path, routes_distances = create_routes_dicts(routes_dict, secs_between_stops, orig_map_nodes_num);
	add_edges_and_weights_to_map!(m = m
									, TTC_graph = gg
									, stop_id_to_node_id_dict = stop_id_to_node_id_dict
									, node_id_stop_id_dict = node_id_stop_id_dict
									, wait_time = Int64(TTC_car_freq/2)
									, orig_map_nodes_num = orig_map_nodes_num
									, agents_speed_in_sec_per_m = agents_speed_in_sec_per_m);
	if serialize_results
		path_map = "./data_output/map_prepared_$TTC_car_freq.....$agents_speed_in_sec_per_m.bin"
		path_routes_path = "./data_output/routes_path_prepared_$TTC_car_freq.....$agents_speed_in_sec_per_m.bin"
		path_routes_distances = "./data_output/routes_distances_prepared_$TTC_car_freq.....$agents_speed_in_sec_per_m.bin"
		path_orig_map_nodes_num = "./data_output/orig_map_nodes_num_prepared_$TTC_car_freq.....$agents_speed_in_sec_per_m.bin"
		serialize(path_map, m)
		serialize(path_routes_path, routes_path)
		serialize(path_routes_distances, routes_distances)
		serialize(path_orig_map_nodes_num, orig_map_nodes_num)
		println("Map object `m` was saved in: ", path_map)
		println("Object `routes_path` was saved in: ", path_routes_path)
		println("Object `routes_distances` was saved in: ", path_routes_distances)
		println("Object `orig_map_nodes_num` was saved in: ", path_orig_map_nodes_num)
		println("OSM map prepared")
	end
    return m, orig_map_nodes_num, routes_path, routes_distances
end
println("Read and prepare OSM map functions were read.")

max_stars_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_issues_repo_name: NykPol/EpidemicInUrbanNetworkToronto

max_forks_repo_name: NykPol/EpidemicInUrbanNetworkToronto


text:
function add_noise(p::InfParams, l = 0.10)
    inf = copy(p.infectiousness)
    p_test = copy(p.pos_test_probs)

    for (i,d) in enumerate(inf)
        k,θ = Distributions.params(d)
        k′ = k*(rand()*l+(1-l/2))
        θ′ = θ*(rand()*l+(1-l/2))
        inf[i] = Gamma(k′,θ′)

        λ = (k′*θ′) / (k*θ)
        p_test[i] *= λ
        p_test[i] = min(p_test[i], 1.0)
    end

    return InfParams(
        p_test,
        p.symptom_probs,
        p.asymptomatic_prob,
        p.symptomatic_isolation_prob,
        inf
    )
end

function param_noise(s::S, l) where S <: CovidState
    return S(statevars(s)..., add_noise(s.params, l), s.prev_action)
end

mutable struct NoisyCovidFilter{PM, P, PMEM, RNG<:AbstractRNG} <: Updater
    predict_model::PM
    particles::P
    weights::Vector{Float64}
    _particle_mem::PMEM
    _unique::Int
    rng::RNG
end

function NoisyCovidFilter(pomdp::POMDP{S}, n::Int) where S<:CovidState
    return NoisyCovidFilter(
        pomdp,
        ParticleCollection(Vector{S}(undef,n)),
        Vector{Float64}(undef, n),
        ParticleCollection(Vector{S}(undef,n)),
        n,
        Random.GLOBAL_RNG
    )
end

function noise_width(N::Int, Un::Int)
    return ((0.05 - 0.25)/(N-1))*(Un-N) + 0.05
end

function ParticleFilters.predict!(up::NoisyCovidFilter, b, a, o)
    pomdp = up.predict_model
    pm = up.particles.particles
    rng = up.rng

    N = length(pm)
    un = up._unique
    noise_param = noise_width(N, un)

    for (i,s) in enumerate(particles(b))
        sp, o, r = POMDPs.gen(pomdp, s, a, rng)
        pm[i] = param_noise(sp, noise_param)
    end
    return up.particles
end

function ParticleFilters.reweight!(up::NoisyCovidFilter, b, a, bp, o)
    pomdp = up.predict_model
    p1 = b.particles
    p2 = bp.particles
    wm = up.weights
    rng = up.rng

    w_sum = 0.0
    for i in eachindex(p1)
        s = p1[i]
        sp = p2[i]
        dist = POMDPs.observation(pomdp, s, a, sp)
        w = pdf(dist, o)
        w_sum += w
        wm[i] = w
    end

    return wm ./= w_sum # normalized
end

function ParticleFilters.resample(up::NoisyCovidFilter, p::Vector{S}, w::Vector{Float64}) where S
    N = length(p)
    rng = up.rng

    ps = Array{S}(undef, N)
    r = rand(rng)/N
    c = w[1]
    i = 1
    U = r
    for m in 1:N
        while U > c && i < N
            i += 1
            c += w[i]
        end
        U += inv(N)
        s = p[i]
        ps[m] = s
    end
    up._unique = length(unique(ps))
    return ParticleCollection(ps)
end

function ParticleFilters.update(up::NoisyCovidFilter, b::ParticleCollection, a::CovidAction, o::Int)
    b′ = predict!(up, b, a, o)
    w = reweight!(up, b, a, b′, o)
    return resample(up, b′.particles, w)
end

ParticleFilters.initialize_belief(pf::NoisyCovidFilter, b::ParticleCollection) = b

function ParticleFilters.initialize_belief(pf::NoisyCovidFilter, b)
    p = pf.particles.particles
    for i in eachindex(p)
        p[i] = rand(b)
    end
    return deepcopy(pf.particles)
end


text:
using DataFrames
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/AdaptationCosts.jl")

    adaptationcosts_noneconomic = addadaptationcosts_noneconomic(m)
    adaptationcosts_noneconomic[:y_year_0] = 2015.
    adaptationcosts_noneconomic[:y_year] = Mimi.dim_keys(m.md, :time)
    adaptationcosts_noneconomic[:gdp] = readpagedata(m, "test/validationdata/$valdir/gdp.csv")

    adaptationcosts_economic = addadaptationcosts_economic(m)
    adaptationcosts_economic[:y_year_0] = 2015.
    adaptationcosts_economic[:y_year] = Mimi.dim_keys(m.md, :time)
    adaptationcosts_economic[:gdp] = readpagedata(m, "test/validationdata/$valdir/gdp.csv")

    adaptationcosts_sealevel = addadaptationcosts_sealevel(m)
    adaptationcosts_sealevel[:y_year_0] = 2015.
    adaptationcosts_sealevel[:y_year] = Mimi.dim_keys(m.md, :time)
    adaptationcosts_sealevel[:gdp] = readpagedata(m, "test/validationdata/$valdir/gdp.csv")

    p = load_parameters(m)
    set_leftover_params!(m, p)

    run(m)

    autofac = m[:AdaptiveCostsEconomic, :autofac_autonomouschangefraction]
    atl_economic = m[:AdaptiveCostsEconomic, :atl_adjustedtolerablelevel]
    imp_economic = m[:AdaptiveCostsEconomic, :imp_adaptedimpacts]
    acp_economic = m[:AdaptiveCostsEconomic, :acp_adaptivecostplateau]
    aci_economic = m[:AdaptiveCostsEconomic, :aci_adaptivecostimpact]

    autofac_compare = readpagedata(m, "test/validationdata/autofac_autonomouschangefraction.csv")
    atl_economic_compare = readpagedata(m, "test/validationdata/atl_adjustedtolerableleveloftemprise_market.csv")
    imp_economic_compare = readpagedata(m, "test/validationdata/$valdir/imp_actualreduction_market.csv")
    acp_economic_compare = readpagedata(m, "test/validationdata/acp_adaptivecostplateau_economic.csv")
    aci_economic_compare = readpagedata(m, "test/validationdata/aci_adaptivecostimpact_economic.csv")

    @test autofac ≈ autofac_compare rtol=1e-8
    @test atl_economic ≈ atl_economic_compare rtol=1e-8
    @test imp_economic ≈ imp_economic_compare rtol=1e-8
    @test acp_economic ≈ acp_economic_compare rtol=1e-3
    @test aci_economic ≈ aci_economic_compare rtol=1e-3

    ac_noneconomic = m[:AdaptiveCostsNonEconomic, :ac_adaptivecosts]
    ac_economic = m[:AdaptiveCostsEconomic, :ac_adaptivecosts]
    ac_sealevel = m[:AdaptiveCostsSeaLevel, :ac_adaptivecosts]

    ac_noneconomic_compare = readpagedata(m, "test/validationdata/ac_adaptationcosts_noneconomic.csv")
    ac_economic_compare = readpagedata(m, "test/validationdata/ac_adaptationcosts_economic.csv")
    ac_sealevel_compare = readpagedata(m, "test/validationdata/$valdir/ac_adaptationcosts_sealevelrise.csv")

    @test ac_noneconomic ≈ ac_noneconomic_compare rtol=1e-2
    @test ac_economic ≈ ac_economic_compare rtol=1e-3
    @test ac_sealevel ≈ ac_sealevel_compare rtol=1e-2
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
export ArduinoTemperatureSensor, ArduinoTemperatureSensorParams
# TODO comment relevant Arduino code once added to project
Base.@kwdef struct ArduinoTemperatureSensorParams <: DeviceParams
  portAdress::String
  numSensors::Int
  maxTemps::Vector{Int}
  selectSensors::Vector{Int}
  nameSensors::Vector{String}

  commandStart::String = "!"
  commandEnd::String = "*"
  pause_ms::Int = 30
  timeout_ms::Int = 1000
  delim::String = "#"
  baudrate::Integer = 115200
  ndatabits::Integer = 8
  parity::SPParity = SP_PARITY_NONE
  nstopbits::Integer = 1
end
ArduinoTemperatureSensorParams(dict::Dict) = params_from_dict(ArduinoTemperatureSensorParams, dict)

Base.@kwdef mutable struct ArduinoTemperatureSensor <: TemperatureSensor
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::ArduinoTemperatureSensorParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String,Union{Device,Missing}}

  ard::Union{SimpleArduino, Nothing} = nothing # Use composition as multiple inheritance is not supported
end

neededDependencies(::ArduinoTemperatureSensor) = []
optionalDependencies(::ArduinoTemperatureSensor) = []

function init(sensor::ArduinoTemperatureSensor)
  @info "Initializing ArduinoTemperatureSensor with ID $(sensor.deviceID)"

  params = sensor.params
  spTU = SerialPort(params.portAdress)
  open(spTU)
  set_speed(spTU, params.baudrate)
  set_frame(spTU,ndatabits=params.ndatabits,parity=params.parity,nstopbits=params.nstopbits)
  #set_flow_control(spTU,rts=rts,cts=cts,dtr=dtr,dsr=dsr,xonxoff=xonxoff)
  sleep(2) # TODO why this sleep?
  flush(spTU)
  write(spTU, "!ConnectionEstablished*#")
  response=readuntil(spTU, Vector{Char}(params.delim), params.timeout_ms);
  @info response
  if(!(response == "ArduinoTemperatureUnitV2") )
      close(spTU)
      throw(ScannerConfigurationError(string("Connected to wrong Device", response)))
    else
      @info "Connection to ArduinoTU established."
  end

  sd = SerialDevice(spTU, params.pause_ms, params.timeout_ms, params.delim, params.delim)
  ard = SimpleArduino(;commandStart = params.commandStart, commandEnd = params.commandEnd, delim = params.delim, sd = sd)
  sensor.ard = ard
  setMaximumTemps(sensor, params.maxTemps)
  sensor.present = true
end

function numChannels(sensor::ArduinoTemperatureSensor)
  return length(sensor.params.selectSensors)
end

function getTemperatures(sensor::ArduinoTemperatureSensor; names::Bool=false)
  temp = retrieveTemps(sensor)
  if length(temp) == sensor.params.numSensors
      if names
          return [temp[sensor.params.selectSensors] sensor.params.nameSensors[sensor.params.selectSensors]]
      else
          return temp[sensor.params.selectSensors]
      end
  else
      return zeros(length(sensor.params.selectSensors))
  end
end

function getTemperature(sensor::ArduinoTemperatureSensor, channel::Int)
  temp = retrieveTemps(sensor)
  return temp[channel]
end

function retrieveTemps(sensor::ArduinoTemperatureSensor)
  TempDelim = ","

  Temps = sendCommand(sensor.ard, "GET:ALLTEMPS")
  Temps = Temps[7:end]  #filter out "TEMPS:" at beginning of answer

  result =  tryparse.(Float64,split(Temps,TempDelim))
  result = map(x -> isnothing(x) ? 0.0 : x, result)
  return result
end

function setMaximumTemps(sensor::ArduinoTemperatureSensor, maxTemps::Array)
    if length(maxTemps) == sensor.params.numSensors
        maxTempString= join(maxTemps, ",")
        ack = sendCommand(sensor.ard, "SET:MAXTEMPS:<"*maxTempString*">")
        # TODO check ack?
        @info "acknowledge of MaxTemps from TempUnit."
    else
        @warn "Please parse a maximum temperature for each sensor" sensor.params.numSensors
    end
end

function getMaximumTemps(sensor::ArduinoTemperatureSensor)
    println(sendCommand(sensor.ard, "GET:MAXTEMPS"))
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((im+cot(x))^-1*cos(x)^4, x) == :(1//32*(im+-1*cot(x))^-2+5//32*(im+cot(x))^-2+im*(-8*cot(x)+8im)^-1+3*im*(16im+16*cot(x))^-1+-1//16*im*x+-1//24*im*(im+cot(x))^-3)
@test integrate((im+cot(x))^-1*cos(x)^3, x) == :(-1//5*cos(x)^5+-1//3*im*sin(x)^3+1//5*im*sin(x)^5)
@test integrate((im+cot(x))^-1*cos(x)^2, x) == :(1//8*(im+cot(x))^-2+im*(-8*cot(x)+8im)^-1+im*(4im+4*cot(x))^-1+-1//8*im*x)
@test integrate((im+cot(x))^-1*cos(x), x) == :(-1//3*cos(x)^3+-1//3*im*sin(x)^3)
@test integrate((im+cot(x))^-1*sec(x), x) == :(-1*cos(x)+im*sin(x)+-1*im*arctanh(sin(x)))
@test integrate((im+cot(x))^-1*sec(x)^2, x) == :(-1*log(sin(x))+im*x+-1*im*tan(x)+log(tan(x)))
@test integrate((im+cot(x))^-1*sec(x)^3, x) == :((1/2)*im*arctanh(sin(x))+-1//2*im*sec(x)*tan(x)+sec(x))
@test integrate((im+cot(x))^-1*sec(x)^4, x) == :((1/2)*tan(x)^2+-1//3*im*tan(x)^3)
@test integrate((im+cot(x))^-1*sec(x)^5, x) == :(1//3*sec(x)^3+1//8*im*arctanh(sin(x))+-1//4*im*sec(x)^3*tan(x)+1//8*im*sec(x)*tan(x))
@test integrate((im+cot(x))^-1*sec(x)^6, x) == :((1/2)*tan(x)^2+1//4*tan(x)^4+-1//3*im*tan(x)^3+-1//5*im*tan(x)^5)
@test integrate((a+b*cot(x))^-1*cos(x)^4, x) == :(-1*(4*a^2+4*b^2)^-1*sin(x)^4*(b+a*cot(x))+1//8*(a^2+b^2)^-2*sin(x)^2*(4*b*(b^2+2*a^2)+a*(b^2+5*a^2)*cot(x))+-1*b*a^4*(a^2+b^2)^-3*log(a*sin(x)+b*cos(x))+1//8*a*x*(a^2+b^2)^-3*(-1*b^4+3*a^4+-6*a^2*b^2))
@test integrate((a+b*cot(x))^-1*cos(x)^3, x) == :(a*(a^2+b^2)^-1*sin(x)+-1*a*(3*a^2+3*b^2)^-1*sin(x)^3+-1*b*(3*a^2+3*b^2)^-1*cos(x)^3+b*a^3*(a^2+b^2)^-5//2*arctanh((a^2+b^2)^-1//2*(a*cos(x)+-1*b*sin(x)))+-1*a*b^2*(a^2+b^2)^-2*sin(x)+-1*b*a^2*(a^2+b^2)^-2*cos(x))
@test integrate((a+b*cot(x))^-1*cos(x)^2, x) == :((2*a^2+2*b^2)^-1*sin(x)^2*(b+a*cot(x))+(1/2)*a*x*(a^2+b^2)^-2*(a^2+-1*b^2)+-1*b*a^2*(a^2+b^2)^-2*log(a*sin(x)+b*cos(x)))
@test integrate((a+b*cot(x))^-1*cos(x), x) == :(a*(a^2+b^2)^-1*sin(x)+-1*b*(a^2+b^2)^-1*cos(x)+a*b*(a^2+b^2)^-3//2*arctanh((a^2+b^2)^-1//2*(a*cos(x)+-1*b*sin(x))))
@test integrate((a+b*cot(x))^-1*sec(x), x) == :(a^-1*arctanh(sin(x))+b*a^-1*(a^2+b^2)^-1//2*arctanh((a^2+b^2)^-1//2*(a*cos(x)+-1*b*sin(x))))
@test integrate((a+b*cot(x))^-1*sec(x)^2, x) == :(a^-1*tan(x)+-1*b*a^-2*log(a+b*cot(x))+-1*b*a^-2*log(tan(x)))
@test integrate((a+b*cot(x))^-1*sec(x)^3, x) == :((1/2)*a^-1*arctanh(sin(x))+a^-3*b^2*arctanh(sin(x))+(1/2)*a^-1*sec(x)*tan(x)+-1*b*a^-2*sec(x)+b*a^-3*(a^2+b^2)^(1/2)*arctanh((a^2+b^2)^-1//2*(a*cos(x)+-1*b*sin(x))))
@test integrate((a+b*cot(x))^-1*sec(x)^4, x) == :(1//3*a^-1*tan(x)^3+a^-3*(a^2+b^2)*tan(x)+-1//2*b*a^-2*tan(x)^2+-1*b*a^-4*(a^2+b^2)*log(a+b*cot(x))+-1*b*a^-4*(a^2+b^2)*log(tan(x)))
@test integrate((1+2*cot(x))^-1*sec(x), x) == :(2//5*5^(1/2)*arctanh(1//5*5^(1/2)*(-2*sin(x)+cos(x)))+arctanh(sin(x)))

max_stars_repo_path: test/4 Trig functions/4.4 Cotangent/4.4.1.3 (d cos)^m (a+b cot)^n.jl


text:
# implementaion of gjk algorithm
# https://en.wikipedia.org/wiki/Gilbert%E2%80%93Johnson%E2%80%93Keerthi_distance_algorithm

using Base.Cartesian

type_immutable{T, n}(::Type{FlexibleSimplex{T}}, ::Type{Val{n}}) = Simplex{n,T}
type_immutable(s, v=Val{length(s)}) = type_immutable(typeof(s), v)
make_immutable(s, v=Val{length(s)}) = convert(type_immutable(s, v), s)

"""
with_immutable{n}(f, s, max_depth::Type{Val{n}}=Val{10})

Apply function f to immutable variant of s without introducing a type instability.
"""
@generated function with_immutable{n}(f, s, max_depth::Type{Val{n}}=Val{10})
    quote
        l = length(s)
        @nif $n d->(d == l) d->(f(make_immutable(s, Val{d}))) d->error("length(s) = $s < max_depth = {n}")
    end
end

function proj_sqdist(pt, fs::FlexibleSimplex)
    with_immutable(fs) do s
        proj_sqdist(pt, s)
    end::Tuple{eltype(fs), numtype(fs)}
end

function weights(pt, fs::FlexibleSimplex)
    with_immutable(fs) do s
        weights(pt, s)
    end
end

"""
MinkowskiDifference(c1, c2)

Represents the Minkowski difference of c1, c2.
"""
immutable MinkowskiDifference{S,T}
    c1::S
    c2::T
end

⊖(c1,c2) = MinkowskiDifference(c1,c2)
convert(::FlexibleConvexHull, md::MinkowskiDifference) = FlexibleConvexHull(vertices(md))

function vertices(m::MinkowskiDifference)
    # There might be redundant vertices, should we try
    # to get rid of them?
    c1, c2 = m.c1, m.c2
    T = promote_type(eltype(c1), eltype(c2))
    verts = T[]
    for v1 in vertices(c1), v2 in vertices(c2)
        push!(verts, v1-v2)
    end
    verts
end

any_inside(c::AbstractConvexHull) = first(vertices(c))
any_inside(m::MinkowskiDifference) = any_inside(m.c1) - any_inside(m.c2)
any_inside(c::Union{HyperCube, HyperRectangle}) = origin(c)
any_inside(v::Vec) = v

support_vector_max(ch::AbstractConvexHull, v) = argmax(x-> x⋅v, vertices(ch))
support_vector_max(w::Vec, v) = w, w⋅v
function support_vector_max(c::HyperRectangle, v)
    s = widths(c)
    takes = map(x -> x > zero(eltype(v)), s .* v)
    best_pt::typeof(v) = origin(c) + takes.*s
    score = (best_pt ⋅ v)
    return best_pt, score
end
support_vector_max(c::HyperCube) = support_vector_max(HyperRectangle(c))

function support_vector_max(m::MinkowskiDifference, v)
    v1, score1 =  support_vector_max(m.c1, v)
    v2, score2 =  support_vector_max(m.c2, -v)
    v1 - v2, score1 + score2
end
support_vector(ch, v) = support_vector_max(ch,v)[1]

"""
s = shrink!(s::FlexibleSimplex, pt, atol=0.)

Drop as many vertices from s as possible, while keeping pt inside.
"""
function shrink!(s::FlexibleSimplex, pt, atol=0.)
    (nvertices(s) <= 1) && return s
    w = weights(pt, s)
    @assert all(w .>= -atol)
    for i in length(w):-1:1
        w[i] .<= atol && deleteat!(s, i)
    end
    s
end

"""
pt_best, dist = gjk0(c)

Compute the distance of a convex set to the origin using gjk algorithm. If c
is not convex, the algorithm may or may not yield an optimal solution.

# Arguments
* c needs to implement the following methods:

1. pt = any_inside(c): returns a point that is guaranteed to lie inside c.

2. best_pt, best_score = support_vector_max(c, v): returns a point inside c,
which has maximal scalar product with v.

# Return
* dist: distance from c to the origin.
* pt_best: point inside c realizing dist.
"""
function gjk0(c,
    atol=1e-6,
    max_iter=100,
    pt_best=any_inside(c)
    )
    @assert atol >= 0

    T = typeof(pt_best)
    zero_point = zero(T)
    sim = FlexibleSimplex(T[])
    for k in 1:max_iter
        direction = -pt_best

        wk, score = support_vector_max(c, direction)
        if (score <= (pt_best ⋅ direction) + atol) # pt_best is already most extreme
            return pt_best, norm(pt_best)
        else
            push!(sim, wk)
            pt_best::T, sqd = proj_sqdist(zero_point, sim)
            sqd == 0 && return pt_best, norm(pt_best)
            #@show weights(pt_best, sim)
            shrink!(sim, pt_best, atol)
        end
    end
    pt_best, norm(pt_best)
end

"""
dist = gjk(c1,c2, args...)

Compute the euclidean between two convex sets using gjk algorithm. See gjk0.
"""
function gjk(c1,c2, args...)
    # can tweak this to also reconstruct pt1 in ch1, pt2 in ch2,
    # which realize dist?
    m = c1 ⊖ c2
    pt_best, dist = gjk0(m, args...)
    dist
end


text:
module Korg
    export synthesize, read_linelist, read_model_atmosphere

    _data_dir = joinpath(@__DIR__, "../data")

    include("constants.jl")                #physical constants
    include("atomic_data.jl")              #symbols and atomic weights
    include("isotopic_abundances.jl")      #self-explanatory
    include("linelist.jl")                 #parse linelists, define Line type
    include("line_opacity.jl")             #opacity, line profile, voigt function
    include("read_statmech_quantities.jl") #approximate Us, Ks, chis
    include("statmech.jl")                 #statistical mechanics, molecular equilibrium
    include("atmosphere.jl")               #parse model atmospheres
    include("transfer.jl")                 #radiative transfer integrals

    #load data when the package is imported.
    const ionization_energies = setup_ionization_energies()
    const partition_funcs = setup_partition_funcs()
    const equilibrium_constants = setup_equilibrium_constants()
    const hline_stark_profiles = setup_hydrogen_stark_profiles()

    include("utils.jl")                                #functions to apply LSF, vac<->air wls, etc.
    include("continuum_absorption/ContinuumAbsorption.jl")  #Define continuum absorption functions.
    include("synthesize.jl")                                #solve radiative transfer equation

end # module