text:
struct nSunsCAP3 <: AbstractProgramme end
typeProgramme = nSunsCAP3()
nameProgramme = "nSunsCAP3"
exerProg = Dict()
### Main lifts
BenchT1 = Exercise(
    name = bench["flat"],
    equipment = equipment["bb"],
    muscles = [muscles["pecs"], muscles["triceps"], muscles["front delts"]],
    trainingMax = 92.5,
)
BenchT2 = Exercise(
    name = bench["cg"],
    equipment = equipment["bb"],
    muscles = [muscles["pecs"], muscles["triceps"], muscles["front delts"]],
    trainingMax = 75,
)
push!(exerProg, "BenchT1" => (exercise = BenchT1, progression = CAP3_Bench_T1))
push!(exerProg, "BenchT2" => (exercise = BenchT2, progression = CAP3_Bench_T2))

DeadliftT1 = Exercise(
    name = deadlift["dl"],
    equipment = equipment["bb"],
    muscles = [
        muscles["hamstrings"],
        muscles["glutes"],
        muscles["spinal erectors"],
        muscles["quads"],
        muscles["traps"],
        muscles["lats"],
        muscles["finger flexors"],
    ],
    trainingMax = 162.5,
)
DeadliftT2 = Exercise(
    name = deadlift["sumo"],
    equipment = equipment["bb"],
    muscles = [
        muscles["quads"],
        muscles["glutes"],
        muscles["hamstrings"],
        muscles["spinal erectors"],
        muscles["traps"],
        muscles["lats"],
        muscles["finger flexors"],
    ],
    trainingMax = 140,
)
push!(exerProg, "DeadliftT1" => (exercise = DeadliftT1, progression = CAP3_Deadlift_T1))
push!(exerProg, "DeadliftT2" => (exercise = DeadliftT2, progression = CAP3_Deadlift_T2))
SquatT1 = Exercise(
    name = squat["low"],
    equipment = equipment["bb"],
    muscles = [
        muscles["quads"],
        muscles["glutes"],
        muscles["spinal erectors"],
        muscles["transverse abdominis"],
    ],
    trainingMax = 115,
)
SquatT2 = Exercise(
    name = squat["front"],
    equipment = equipment["bb"],
    muscles = [muscles["quads"], muscles["glutes"], muscles["core"]],
    trainingMax = 80,
)
push!(exerProg, "SquatT1" => (exercise = SquatT1, progression = CAP3_Squat_T1))
push!(exerProg, "SquatT2" => (exercise = SquatT2, progression = CAP3_Squat_T2))

RowT1 = Exercise(
    name = row["pendlay"],
    modality = modality["cheat"],
    equipment = equipment["bb"],
    muscles = [
        muscles["lats"],
        muscles["biceps"],
        muscles["rear delts"],
        muscles["rhomboids"],
        muscles["finger flexors"],
    ],
    trainingMax = 132.5,
)
RowT2 = Exercise(
    name = row["kroc"],
    equipment = equipment["db"],
    muscles = [
        muscles["lats"],
        muscles["biceps"],
        muscles["rear delts"],
        muscles["rhomboids"],
        muscles["finger flexors"],
    ],
    trainingMax = 42.5,
)
push!(exerProg, "RowT1" => (exercise = RowT1, progression = CAP3_Row_T1))
push!(exerProg, "RowT2" => (exercise = RowT2, progression = CAP3_Row_T2))

PressT1 = Exercise(
    name = press["ohp"],
    equipment = equipment["bb"],
    muscles = [muscles["front delts"], muscles["triceps"], muscles["traps"]],
    trainingMax = 60,
)
PressT2 = Exercise(
    name = press["ohp"],
    equipment = equipment["bb"],
    muscles = [muscles["front delts"], muscles["triceps"], muscles["traps"]],
    trainingMax = 60,
)
push!(exerProg, "PressT1" => (exercise = PressT1, progression = CAP3_OHP_T1))
push!(exerProg, "PressT2" => (exercise = PressT2, progression = CAP3_OHP_T2))

## Weeks
import Lifting: makeDays
function makeDays(::nSunsCAP3, dict::Dict)
    week1 = [[], [], [], [], [], [], []]
    week2 = [[], [], [], [], [], [], []]
    week3 = [[], [], [], [], [], [], []]

    push!(week1[1], dict["BenchT1"][1], dict["BenchT1"][2], 1)
    push!(week1[1], dict["BenchT2"][1], dict["BenchT2"][2], 1)

    push!(week1[2], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 1)
    push!(week1[2], dict["RowT2"][1], dict["RowT2"][2], 1)

    push!(week1[3], dict["SquatT1"][1], dict["SquatT1"][2], 1)
    push!(week1[3], dict["PressT2"][1], dict["PressT2"][2], 1)

    push!(week1[4], dict["BenchT1"][1], dict["BenchT1"][2], 2)
    push!(week1[4], dict["BenchT2"][1], dict["BenchT2"][2], 2)

    push!(week1[5], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 1)
    push!(week1[5], dict["RowT1"][1], dict["RowT1"][2], 1)

    push!(week1[6], dict["SquatT2"][1], dict["SquatT2"][2], 1)
    push!(week1[6], dict["PressT1"][1], dict["PressT1"][2], 1)

    # Week 2
    push!(week2[1], dict["BenchT1"][1], dict["BenchT1"][2], 3)
    push!(week2[1], dict["BenchT2"][1], dict["BenchT2"][2], 3)

    push!(week2[2], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 2)
    push!(week2[2], dict["RowT2"][1], dict["RowT2"][2], 2)

    push!(week2[3], dict["SquatT1"][1], dict["SquatT1"][2], 2)
    push!(week2[3], dict["PressT2"][1], dict["PressT2"][2], 2)

    push!(week2[4], dict["BenchT1"][1], dict["BenchT1"][2], 4)
    push!(week2[4], dict["BenchT2"][1], dict["BenchT2"][2], 4)

    push!(week2[5], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 2)
    push!(week2[5], dict["RowT1"][1], dict["RowT1"][2], 2)

    push!(week2[6], dict["SquatT2"][1], dict["SquatT2"][2], 2)
    push!(week2[6], dict["PressT1"][1], dict["PressT1"][2], 2)

    # Week3
    push!(week3[1], dict["BenchT1"][1], dict["BenchT1"][2], 5)
    push!(week3[1], dict["BenchT2"][1], dict["BenchT2"][2], 5)

    push!(week3[2], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 3)
    push!(week3[2], dict["RowT2"][1], dict["RowT2"][2], 3)

    push!(week3[3], dict["SquatT1"][1], dict["SquatT1"][2], 3)
    push!(week3[3], dict["PressT2"][1], dict["PressT2"][2], 3)

    push!(week3[4], dict["BenchT1"][1], dict["BenchT1"][2], 6)
    push!(week3[4], dict["BenchT2"][1], dict["BenchT2"][2], 6)

    push!(week3[5], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 3)
    push!(week3[5], dict["RowT1"][1], dict["RowT1"][2], 3)

    push!(week3[6], dict["SquatT2"][1], dict["SquatT2"][2], 3)
    push!(week3[6], dict["PressT1"][1], dict["PressT1"][2], 3)

    return [
        week1[1],
        week1[2],
        week1[3],
        week1[4],
        "Rest",
        week1[5],
        week1[6],
        week2[1],
        week2[2],
        "Rest",
        week2[3],
        week2[4],
        week2[5],
        week2[6],
        "Rest",
        week3[1],
        week3[2],
        week3[3],
        week3[4],
        "Rest",
        week3[5],
        week3[6],
        "Rest",
    ]
end

days = makeDays(nSunsCAP3(), exerProg)
push!(
    Lifting_Programmes,
    nameProgramme => Programme(typeProgramme, nameProgramme, exerProg, days),
)

ext: jl

max_stars_repo_name: dcelisgarza/Lifting


text:
@static if Base.VERSION <= v"1.0.5"
    # https://github.com/JuliaLang/julia/pull/29442
    _oneunit(::CartesianIndex{N}) where {N} = _oneunit(CartesianIndex{N})
    _oneunit(::Type{CartesianIndex{N}}) where {N} = CartesianIndex(ntuple(x -> 1, Val(N)))
else
    const _oneunit = Base.oneunit
end

# Once Base has colon defined here we can replace this
_colon(I::CartesianIndex{N}, J::CartesianIndex{N}) where N =
    CartesianIndices(map((i,j) -> i:j, Tuple(I), Tuple(J)))


text:
# Tests on Multivariate Normal distributions

import NumericExtensions
import PDMats: ScalMat, PDiagMat, PDMat

using Distributions
using Base.Test

##### construction, basic properties, and evaluation

mu = [1., 2., 3.]
va = [1.2, 3.4, 2.6]
C = [4. -2. -1.; -2. 5. -1.; -1. -1. 6.]

x1 = [3.2, 1.8, 2.4]
x = rand(3, 100)

# SGauss

gs = gmvnormal(mu, sqrt(2.0))
@test isa(gs, IsoNormal)
@test dim(gs) == 3
@test mean(gs) == mode(gs) == mu
@test_approx_eq cov(gs) diagm(fill(2.0, 3))
@test var(gs) == diag(cov(gs))
@test_approx_eq entropy(gs) 0.5 * logdet(2π * e * cov(gs))

gsz = gmvnormal(3, sqrt(2.0))
@test isa(gsz, IsoNormal)
@test dim(gsz) == 3
@test mean(gsz) == zeros(3)
@test gsz.zeromean

# DGauss

gd = gmvnormal(mu, sqrt(va))
@test isa(gd, DiagNormal)
@test dim(gd) == 3
@test mean(gd) == mode(gd) == mu
@test_approx_eq cov(gd) diagm(va)
@test var(gd) == diag(cov(gd))
@test_approx_eq entropy(gd) 0.5 * logdet(2π * e * cov(gd))

gdz = gmvnormal(PDiagMat(va))
@test isa(gdz, DiagNormal)
@test dim(gdz) == 3
@test mean(gdz) == zeros(3)
@test gdz.zeromean

# Gauss

gf = MvNormal(mu, C)
@test isa(gf, MvNormal)
@test dim(gf) == 3
@test mean(gf) == mode(gf) == mu
@test cov(gf) == C
@test var(gf) == diag(cov(gf))
@test_approx_eq entropy(gf) 0.5 * logdet(2π * e * cov(gf))

gfz = MvNormal(C)
@test isa(gfz, MvNormal)
@test dim(gfz) == 3
@test mean(gfz) == zeros(3)
@test gfz.zeromean

##### LogPDF/PDF evaluation

@test_approx_eq_eps logpdf(gs, x1) -5.106536370454 1.0e-8
@test_approx_eq_eps logpdf(gd, x1) -6.029399605174 1.0e-8
@test_approx_eq_eps logpdf(gf, x1) -5.680452770982 1.0e-8

@test_approx_eq_eps logpdf(gsz, x1) -8.606536370454 1.0e-8
@test_approx_eq_eps logpdf(gdz, x1) -9.788449378930 1.0e-8
@test_approx_eq_eps logpdf(gfz, x1) -9.621416626404 1.0e-8

n = size(x, 2)
r = zeros(n)

for i = 1:n; r[i] = logpdf(gs, x[:,i]); end
@test_approx_eq logpdf(gs, x) r
@test_approx_eq pdf(gs, x) exp(r)

for i = 1:n; r[i] = logpdf(gsz, x[:,i]); end
@test_approx_eq logpdf(gsz, x) r
@test_approx_eq pdf(gsz, x) exp(r)

for i = 1:n; r[i] = logpdf(gd, x[:,i]); end
@test_approx_eq logpdf(gd, x) r
@test_approx_eq pdf(gd, x) exp(r)

for i = 1:n; r[i] = logpdf(gdz, x[:,i]); end
@test_approx_eq logpdf(gdz, x) r
@test_approx_eq pdf(gdz, x) exp(r)

for i = 1:n; r[i] = logpdf(gf, x[:,i]); end
@test_approx_eq logpdf(gf, x) r
@test_approx_eq pdf(gf, x) exp(r)

for i = 1:n; r[i] = logpdf(gfz, x[:,i]); end
@test_approx_eq logpdf(gfz, x) r
@test_approx_eq pdf(gfz, x) exp(r)

##### Sampling

x = rand(gs)
@test isa(x, Vector{Float64})
@test length(x) == dim(gs)

x = rand(gd)
@test isa(x, Vector{Float64})
@test length(x) == dim(gd)

x = rand(gf)
@test isa(x, Vector{Float64})
@test length(x) == dim(gf)

x = rand(gsz)
@test isa(x, Vector{Float64})
@test length(x) == dim(gsz)

x = rand(gdz)
@test isa(x, Vector{Float64})
@test length(x) == dim(gdz)

x = rand(gfz)
@test isa(x, Vector{Float64})
@test length(x) == dim(gfz)

n = 10
x = rand(gs, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gs), n)

x = rand(gd, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gd), n)

x = rand(gf, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gf), n)

x = rand(gsz, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gsz), n)

x = rand(gdz, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gdz), n)

x = rand(gfz, n)
@test isa(x, Matrix{Float64})
@test size(x) == (dim(gfz), n)

##### MLE

# a slow but safe way to implement MLE for verification

function _gauss_mle(x::Matrix{Float64})
	mu = vec(mean(x, 2))
	z = x .- mu
	C = (z * z') * (1/size(x,2))
	return mu, C
end

function _gauss_mle(x::Matrix{Float64}, w::Vector{Float64})
	sw = sum(w)
	mu = (x * w) * (1/sw)
	z = x .- mu
	C = (z * scale(w, z')) * (1/sw)
	Base.LinAlg.copytri!(C, 'U')
	return mu, C
end

x = randn(3, 200) .+ randn(3) * 2.
w = rand(200)

g = fit(MvNormal, x)
mu, C = _gauss_mle(x)
@test_approx_eq mean(g) mu
@test_approx_eq cov(g) C

g = fit_mle(MvNormal, x, w)
mu, C = _gauss_mle(x, w)
@test_approx_eq mean(g) mu
@test_approx_eq cov(g) C

g = fit(IsoNormal, x)
mu, C = _gauss_mle(x)
@test_approx_eq g.μ mu
@test_approx_eq g.Σ.value mean(diag(C))

g = fit_mle(IsoNormal, x, w)
mu, C = _gauss_mle(x, w)
@test_approx_eq g.μ mu
@test_approx_eq g.Σ.value mean(diag(C))

assumed_g = Distributions.IsoNormalKnownSigma(3, 1.)
g = fit_mle(assumed_g, x, w)
mu, C = _gauss_mle(x, w)
@test_approx_eq g.μ mu
@test_approx_eq diag(g.Σ) diag(assumed_g.Σ)

g = fit(DiagNormal, x)
mu, C = _gauss_mle(x)
@test_approx_eq g.μ mu
@test_approx_eq g.Σ.diag diag(C)

g = fit_mle(DiagNormal, x, w)
mu, C = _gauss_mle(x, w)
@test_approx_eq g.μ mu
@test_approx_eq g.Σ.diag diag(C)

##### Sufficient statistics type
n = 3
X = reshape(Float64[1:12], 4, n)
w = rand(n)
Xw = X * diagm(w)

ss = suffstats(MvNormal, X)
ssw = suffstats(MvNormal, X, w)

s_t = sum(X, 2)
ws_t = sum(Xw, 2)
tmp = X .- (s_t ./ n)
ss_t = tmp*tmp'

tmp = X .- sum(Xw, 2) ./ sum(w)
wss_t = (tmp*diagm(w))*tmp'
tw_t = length(w)
wtw_t = sum(w)

@test_approx_eq ss.s s_t
@test_approx_eq ss.m (s_t ./ n)
@test_approx_eq ss.s2 ss_t
@test_approx_eq ss.tw tw_t

@test_approx_eq ssw.s ws_t
@test_approx_eq ssw.m (ws_t ./ wtw_t)
@test_approx_eq ssw.s2 wss_t
@test_approx_eq ssw.tw wtw_t


text:
#
# Autogenerated by Thrift Compiler (0.11.0)
#
# DO NOT EDIT UNLESS YOU ARE SURE THAT YOU KNOW WHAT YOU ARE DOING

struct _enum_TCompletionHintType
  COLUMN::Int32
  TABLE::Int32
  VIEW::Int32
  SCHEMA::Int32
  CATALOG::Int32
  REPOSITORY::Int32
  FUNCTION::Int32
  KEYWORD::Int32
end
const TCompletionHintType = _enum_TCompletionHintType(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7))

mutable struct TCompletionHint <: Thrift.TMsg
  _type::Int32
  hints::Vector{String}
  replaced::String
  TCompletionHint() = (o=new(); fillunset(o); o)
end # mutable struct TCompletionHint


text:
@testset "fit" begin
    dir = @__DIR__
    file = joinpath(dir,  "..", "example", "data", "Malitson.csv")
    data = SellmeierFit.read(file)

    (; mdl, err) = fit_sellmeier(data.λ, data.ε)
    @test err ≤ 1e-13

    mdl1 = fit_sellmeier(data.λ, data.ε, Val(3)).mdl
    @test mdl1 == mdl

    mdl2 = fit_sellmeier(data.λ, data.ε, 2, 1).mdl
    @test mdl2 == mdl
end  # @testset "fit"


text:
using Base.Cartesian
using BitBasis
using Combinatorics
using SparseArrays
using Arpack
using LinearAlgebra
using Plots

function occupation_ed(specs,Nₙs,temp::Float64)

    N=length(Nₙs)
    occ=zeros(Float64,length(Nₙs[1]))
    z=0.0
    for i in 1:N
        Nₙ=Nₙs[i]
        bolzman=sum(exp.(-(1/temp).*specs[i]))
        occ.+=bolzman.*Nₙ
        z+=bolzman
    end
    return occ./z
end

function energy_ed(specs,Nₙs,temp::Float64)

    N=length(Nₙs)
    enr=0.0
    z=0.0
    for i in 1:N
        enr+=sum(exp.(-(1/temp).*specs[i]) .*specs[i] )
        z+=sum(exp.(-(1/temp).*specs[i])  )
    end
    return enr./z
end

function energy_obs(Dₙ::Vector{Float64} ,β::Float64)

    D = Diagonal(Dₙ)
    expD = Diagonal(exp.(-β*Dₙ))
    enr = tr(D*expD)
    return enr
end

function number_obs(Dₙ::Vector{Float64},Nₙ ,β::Float64)

    D = Diagonal(Dₙ)
    expD = Diagonal(exp.(-β*Dₙ))
    N=sum(Nₙ)
    enr = tr(N*expD)
    return enr
end

function partition_obs(Dₙ::Vector{Float64} ,β::Float64)

    expD = Diagonal(exp.(-β*Dₙ))
    par = tr(expD)
    return par

end

function SC_corr_obs(spacial_dims::NTuple{DIMS2,Int64},index_i::CartesianIndex{DIMS},index_j::CartesianIndex{DIMS} , Nₙ , Dₙ::Vector{Float64} , Uₙ::Matrix{Float64} , β::Float64) where {DIMS} where {DIMS2}
    P = SC_corr_sub(spacial_dims,Nₙ, index_i , index_j)
    Pd = Uₙ'*P*Uₙ
    expD = Diagonal(exp.(-β*Dₙ))
    par = tr(Pd*expD)
    return par
end

function SC_corr_sw_tau_obs(spacial_dims::NTuple{DIMS2,Int64},num_species::Int64,index_i::CartesianIndex{DIMS},index_j::CartesianIndex{DIMS}
    , Nₙ , Dₙ::Vector{Float64} , Uₙ::Matrix{Float64} , Nₘ , Dₘ::Vector{Float64} , Uₘ::Matrix{Float64}
    , β::Float64,taus::Vector{Float64} ) where {DIMS} where {DIMS2}

    corr = zeros(Float64,length(taus))
    for c in 1:num_species
        for g in 1:2
            if g==1
                index_i_sp=CartesianIndex(index_i,c) #(index_i...,c)#
                index_j_sp=CartesianIndex(index_j,c) #(index_j...,c)#
            else
                index_i_sp=CartesianIndex(index_j,c) #(index_j...,c)#
                index_j_sp=CartesianIndex(index_i,c)
            end

            Δᵢ⁺ = SC_create_sub(spacial_dims,Nₙ,Nₘ, index_i_sp)
            Δⱼ  = SC_annihilate_sub(spacial_dims,Nₘ,Nₙ, index_j_sp)
            Δᵢ⁺d = Uₘ' *  Δᵢ⁺ * Uₙ
            Δⱼd =  Uₙ' *  Δⱼ  * Uₘ

            for i in 1:length(taus)
                τ=taus[i]
                e⁻⁽ᵝ⁻ᵀ⁾ᴰⁿ = Diagonal(exp.(-(β-τ).*Dₙ))
                e⁻ᵀᴰᵐ =     Diagonal(exp.(-τ*Dₘ))
                corr[i] += tr(e⁻⁽ᵝ⁻ᵀ⁾ᴰⁿ * Δⱼd * e⁻ᵀᴰᵐ * Δᵢ⁺d)
            end

        end
    end
    return corr
end

function SC_corr_dw_tau_obs(spacial_dims::NTuple{DIMS2,Int64},num_species::Int64,rᵢ::CartesianIndex{DIMS},rⱼ::CartesianIndex{DIMS}
    , Nₙ , Dₙ::Vector{Float64} , Uₙ::Matrix{Float64} , Nₘ , Dₘ::Vector{Float64} , Uₘ::Matrix{Float64}
    , β::Float64,taus::Vector{Float64} ) where {DIMS} where {DIMS2}

    dims = (spacial_dims... ,num_species)
    corr = zeros(Float64,2,2,length(taus))
    for c in 1:num_species

        index_i_sp=CartesianIndex(rᵢ,c) #(index_i...,c)#
        index_j_sp=CartesianIndex(rⱼ,c) #(index_j...,c)#
        Δᵢₕ⁺ = []
        Δᵢᵥ⁺ = []
        Δⱼₕ  = []
        Δⱼᵥ = []

        index_i = index_i_sp
        index_iₕ = hop(index_i,1,1,dims) # hop right (x direction)
        index_iᵥ = hop(index_i,2,1,dims) # hop up    (y direction)
        index_j = index_j_sp
        index_jₕ = hop(index_j,1,1,dims) # hop right (x direction)
        index_jᵥ = hop(index_j,2,1,dims) # hop up    (y direction)

        Δᵢₕ⁺  = SC_create_sub(spacial_dims,Nₙ,Nₘ, index_i,index_iₕ)# c^†_{index_i,↑} c^†_{index_iₕ,↓}
        Δᵢₕ⁺ += SC_create_sub(spacial_dims,Nₙ,Nₘ, index_iₕ,index_i)# c^†_{index_iₕ,↑} c^†_{index_i,↓}
        Δⱼₕ   = SC_annihilate_sub(spacial_dims,Nₘ,Nₙ, index_j,index_jₕ) # c_{index_j,↓} c^_{index_jₕ,↑}
        Δⱼₕ  += SC_annihilate_sub(spacial_dims,Nₘ,Nₙ, index_jₕ,index_j) # c_{index_jₕ,↓} c^_{index_j,↑}
        Δᵢᵥ⁺  = SC_create_sub(spacial_dims,Nₙ,Nₘ, index_i,index_iᵥ)# c^†_{index_i,↑} c^†_{index_iᵥ,↓}
        Δᵢᵥ⁺ += SC_create_sub(spacial_dims,Nₙ,Nₘ, index_iᵥ,index_i)# c^†_{index_iᵥ,↑} c^†_{index_i,↓}
        Δⱼᵥ   = SC_annihilate_sub(spacial_dims,Nₘ,Nₙ, index_j,index_jᵥ) # c_{index_j,↓} c^_{index_jᵥ,↑}
        Δⱼᵥ += SC_annihilate_sub(spacial_dims,Nₘ,Nₙ, index_jᵥ,index_j) # c_{index_jᵥ,↓} c^_{index_j,↑}

        Δᵢₕ⁺d  = Uₘ' *  Δᵢₕ⁺ * Uₙ
        Δᵢᵥ⁺d = Uₘ' *  Δᵢᵥ⁺ * Uₙ
        Δⱼₕd  =  Uₙ' *  Δⱼₕ  * Uₘ
        Δⱼᵥd  =  Uₙ' *  Δⱼᵥ  * Uₘ

        for i in 1:length(taus)
            τ=taus[i]
            e⁻⁽ᵝ⁻ᵀ⁾ᴰⁿ = Diagonal(exp.(-(β-τ).*Dₙ))
            e⁻ᵀᴰᵐ =     Diagonal(exp.(-τ*Dₘ))
            corr[1,1,i] += tr(e⁻⁽ᵝ⁻ᵀ⁾ᴰⁿ * Δⱼₕd * e⁻ᵀᴰᵐ * Δᵢₕ⁺d)
            corr[1,2,i] += tr(e⁻⁽ᵝ⁻ᵀ⁾ᴰⁿ * Δⱼₕd * e⁻ᵀᴰᵐ * Δᵢᵥ⁺d)
            corr[2,1,i] += tr(e⁻⁽ᵝ⁻ᵀ⁾ᴰⁿ * Δⱼᵥd * e⁻ᵀᴰᵐ * Δᵢₕ⁺d)
            corr[2,2,i] += tr(e⁻⁽ᵝ⁻ᵀ⁾ᴰⁿ * Δⱼᵥd * e⁻ᵀᴰᵐ * Δᵢᵥ⁺d)
        end

    end
    return corr
end

function observables(spacial_dims::NTuple{DIMS,Int64},num_species::Int64, t::Float64, U::Float64,μ::Float64
    ,βs::Vector{Float64},taus::Vector{Float64}) where {DIMS}

    if length(spacial_dims)==2 && spacial_dims[2]==1
        spacial_dims=spacial_dims=(spacial_dims[1],) #(spacial_dims[1]) #
    end
    num_spin=2
    N_max=prod(spacial_dims)
    nt=0:N_max
    n=[]
    for i in 1:num_spin*num_species

        push!(n,nt)
    end

    Nₙs=collect(Iterators.product(n...))[:]
    # define for the case of num_species=1, update later to general case
    #if (num_species==1)
    #    for Nₙ in Nₙs
    #        if Nₙ[1]!=Nₙ[2]
    #            Nₙ_temp=(Nₙ[2],Nₙ[1])
    #            deleteat!(Nₙs,Nₙs.==Nₙ_temp)
    #        end
    #    end
    #end
    b=length(βs)
    tt=length(taus)
    partitions = zeros(Float64,b)
    energys = zeros(Float64,b)
    numbers = zeros(Float64,b)
    SC_corrs=zeros(Float64,spacial_dims...,b)
    sw_tau=zeros(Float64,spacial_dims...,tt,b)
    dw_tau=zeros(Float64,spacial_dims...,2,2,tt,b)

    for Nₙ in Nₙs
        for i in 1:b #runs over βs
            β=βs[i]
            binos = binomial.(N_max,Nₙ)
            weight  = prod(binos)*exp(-β*μ*(sum(Nₙ)-1))
            if weight<10^(-5)
                break
            end
            if !isempty(taus)
                Nₘ =  zeros(Int64,length(Nₙ))
                for i in 1:length(Nₘ)
                     mod(i,num_species)==0 ? Nₘ[i]=Nₙ[i]+1 : Nₘ[i]=Nₙ[i]
                end
                Nₘ = (Nₘ...,)
                #(Nₙ[1]+1,Nₙ[2],Nₙ[3]+1,Nₙ[4])
                Dₘ =[]
                Uₘ=[]
                if Nₘ in Nₙs
                    Hₘ=hamiltonian_sub(spacial_dims,Nₘ,t,U,μ)
                    Dₘ , Uₘ = eigen(Hₘ)
                end
            end
            Hₙ=hamiltonian_sub(spacial_dims,Nₙ,t,U,μ)
            Dₙ , Uₙ = eigen(Hₙ)
            partitions[i] += partition_obs(Dₙ ,β)
            energys[i]    += energy_obs(Dₙ ,β)
            numbers[i]    += number_obs(Dₙ , Nₙ, β)
            ds=findall(ones(spacial_dims).!=0)
          #  for   xt in 1:spacial_dims[1]
            #    for yt  in 1:spacial_dims[2]
            for r in ds
                r0=ds[1]
                SC_corrs[r,i]+=SC_corr_obs(spacial_dims,CartesianIndex(r0) , CartesianIndex(r) , Nₙ , Dₙ , Uₙ, β)
                if !isempty(taus) && !isempty(Dₘ)
                    sw_tau[r,:,i]+=SC_corr_sw_tau_obs(spacial_dims,num_species,CartesianIndex(r0),CartesianIndex(r)
                    , Nₙ , Dₙ  , Uₙ , Nₘ , Dₘ  , Uₘ , β,taus )
                    if length(spacial_dims)>1
                        dw_tau[r,:,:,:,i]+=SC_corr_dw_tau_obs(spacial_dims,num_species,CartesianIndex(r0),CartesianIndex(r)
                        , Nₙ , Dₙ  , Uₙ , Nₘ , Dₘ  , Uₘ , β,taus )
                    end
                end
            end
               # end
           # end
        end
    end

    for i in 1:b
        SC_corrs[:,:,i] = SC_corrs[:,:,i]./partitions[i]
        sw_tau[:,:,:,i] = sw_tau[:,:,:,i]./partitions[i]
        dw_tau[:,:,:,:,:,i] = dw_tau[:,:,:,:,:,i]./partitions[i]
    end
    return partitions, energys./partitions, numbers./partitions,num_species*SC_corrs, num_species*sw_tau, num_species*dw_tau
end

function spectrum(spacial_dims::NTuple{DIMS,Int64},num_species::Int64, t::Float64, U::Float64,μ::Float64) where {DIMS}
    num_spin=2
    N_max=prod(spacial_dims)
    nt=0:N_max
    n=[]
    for i in 1:num_spin*num_species
        push!(n,nt)
    end
    Dₙ=[]
    #Uₙ=[]
    Nₙs=collect(Iterators.product(n...))[:]
    for Nₙ in Nₙs
        ham=hamiltonian_sub(spacial_dims,Nₙ,t,U,μ)
        E=eigen(ham)
        push!(Dₙ,E.values)
        #push!(Uₙ,E.vectors)
    end
    return Dₙ  , Nₙs
end

function hamiltonian_sub(spacial_dims::NTuple{DIMS,Int64},Nₙ, t::Float64, U::Float64,μ::Float64) where {DIMS}

   states , state_num =  build_states(spacial_dims,Nₙ)
   N=length(states)
   H_sub=zeros(Float64,N,N)
   for ket in 1:N
        state_ket=states[ket]
        state_sp , co = hamiltonian_operator(state_ket,t,U,μ)
        for j in 1:length(state_sp)
            state_bra=state_sp[j]
            bra = findall(x->x==packbits(state_bra[:]),state_num)
            H_sub[ket,bra...]+=co[j]
        end
    end
    return H_sub
end

function SC_corr_sub(spacial_dims::NTuple{DIMS2,Int64},Nₙ, index_i::CartesianIndex{DIMS} , index_j::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}

    states ,state_num =  build_states(spacial_dims,Nₙ)
    N=length(states)

    SC_sub=zeros(Float64,N,N)
    for ket in 1:N
     state_ket=states[ket]
         state_sp , co = SC_corr_operator(state_ket,index_i,index_j)
         for j in 1:length(state_sp)
             state_bra=state_sp[j]
             bra = findall(x->x==packbits(state_bra[:]),state_num)
             SC_sub[ket,bra].+=co[j]
         end
     end
     return SC_sub
 end

function SC_corr_operator(state::Array{Int8,DIMS2} , index_i::CartesianIndex{DIMS}, index_j::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}
    # index_i contains only the spacial index
    dims = size(state)
    spacial_dims = length(index_i)
    num_species = dims[spacial_dims+1]
    state_sp = Array{Int8,length(dims)}[]#state#
    co=Float64[]#0.0#

    for c in 1:num_species
        for g in 1:2
            if g==1
                index_i_sp=CartesianIndex(index_i,c) #(index_i...,c)#
                index_j_sp=CartesianIndex(index_j,c) #(index_j...,c)#
            else
                index_i_sp=CartesianIndex(index_j,c) #(index_j...,c)#
                index_j_sp=CartesianIndex(index_i,c)
            end
            state_temp , ex1 =  SC_create_operatr(state,index_i_sp)
            ex2=0
            ex1 == 0 ? state_f = state_temp : (state_f , ex2) = SC_annihilate_operatr(state_temp, index_j_sp)
            if ex1*ex2 != 0.0
                push!(state_sp,state_f)
                push!(co,ex1*ex2)
            end
        end
    end
    return state_sp , co
end

function hamiltonian_operator(state::Array{Int8,DIMS2} ,  t::Float64,U::Float64,μ::Float64) where {DIMS2}

    dims = size(state)
    spacial_dims = length(dims)-2
    num_species = dims[spacial_dims+1]
    state_sp = Array{Int8,length(dims)}[]
    co=Float64[]

    diag_term = U*sum( sum(state[ntuple(k->:,spacial_dims+1)...,1].-state[ntuple(k->:,spacial_dims+1)...,2] ,dims=spacial_dims+1 ).^2 )/2/num_species #interaction term
    diag_term = -diag_term - μ*sum(state) #chemical potential term
    push!(state_sp,state)
    push!(co,diag_term)

    #hopping term
    occupations = findall(x->x==1, state) # indices which are  occupied
    for occupation in occupations
        for dir in 1:spacial_dims
            for lr in 1:2
                index_i = occupation
                index_j = hop(index_i,dir,lr,dims)
                state_f , co_temp = hopping_operatr(state,index_i,index_j)
                if co_temp != 0 && state_f ∉ state_sp
                    push!(state_sp,state_f)
                    push!(co,-t*co_temp)
                end
            end
        end
    end
    return state_sp , co
end

function create(state::Array{Int8,DIMS2} ,index::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}
    # creation operator at (x,y,species,spin)
    new_state = copy(state)
    ex = 0  # denotes if the state was annihilated by the creation operator.
    if state[index] == 0
        new_state[index] = 1
        ex = check_fermion_comm(state ,index)#1
    end
        return new_state , ex
end

function annihilate(state::Array{Int8,DIMS2} ,index::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}
    # annihilation operator at (x,y,spin,species)
    new_state = copy(state)
    ex = 0  # denotes if the state was annihilated by the annihilation operator.
    if state[index] == 1
        new_state[index] = 0
        ex = check_fermion_comm(state ,index)#1
    end
        return new_state , ex
end

function check_fermion_comm(state::Array{Int8,DIMS2} ,index::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}
    state_l=LinearIndices(state)
    index_l=state_l[index]
    sign = sum(state[:][index_l+1:end])
    ex = (-1)^(sign-1)#(-1)^(sign-1)
    return ex
end

function hopping_operatr(state::Array{Int8,DIMS2} ,index_i::CartesianIndex{DIMS},index_j::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}
    #hopping from  (x_i,y_i,species_i,spin_i) to  (x_j,y_j,species_j,spin_j)
    state_temp , ex1 = annihilate(state,index_i)
    ex1 == 0 ? state_f = state_temp : (state_f , ex2) = create(state_temp, index_j)
    return state_f , ex1*ex2

end

function SC_create_operatr(state::Array{Int8,DIMS2} ,index_i::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}
    #sc operator c^†_{index_i,↑} c^†_{index_i,↓} where   index_i=(x_i,y_i,species_i)
    up_s   = CartesianIndex(index_i,1)
    down_s = CartesianIndex(index_i,2)

     state_temp , ex1 = create(state,down_s)
    ex2=0
    ex1 == 0 ? state_f = state_temp : (state_f , ex2) = create(state_temp, up_s)
    return state_f , ex1*ex2
end

function SC_create_operatr(state::Array{Int8,DIMS2} ,index_i::CartesianIndex{DIMS},index_j::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}
    #sc operator c^†_{index_i,↑} c^†_{index_j,↓} where   index_i=(x_i,y_i,species_i)
    up_s   = CartesianIndex(index_i,1)
    down_s = CartesianIndex(index_j,2)

     state_temp , ex1 = create(state,down_s)
    ex2=0
    ex1 == 0 ? state_f = state_temp : (state_f , ex2) = create(state_temp, up_s)
    return state_f , ex1*ex2
end

function SC_create_sub(spacial_dims::NTuple{DIMS2,Int64},Nₙ, Nₘ, index_i::CartesianIndex{DIMS} ) where {DIMS} where {DIMS2}

    statesₙ ,state_numₙ =  build_states(spacial_dims,Nₙ)
    statesₘ ,state_numₘ =  build_states(spacial_dims,Nₘ)
    N=length(statesₙ)
    M=length(statesₘ)
    mat=zeros(Float64,M,N)#spzeros(Float64,N,N)
    for ket in 1:N
        state_ket=statesₙ[ket]
        state_bra , co = SC_create_operatr(state_ket,index_i)
        bra = findall(x->x==packbits(state_bra[:]),state_numₘ)
        mat[bra,ket].+=co
     end
     return mat
 end

 function SC_create_sub(spacial_dims::NTuple{DIMS2,Int64},Nₙ, Nₘ, index_i::CartesianIndex{DIMS}, index_j::CartesianIndex{DIMS} ) where {DIMS} where {DIMS2}
    # matrix representation of the operator c^†_{index_i,↑} c^†_{index_j,↓} where   index_i=(x_i,y_i,species_i)
    statesₙ ,state_numₙ =  build_states(spacial_dims,Nₙ)
    statesₘ ,state_numₘ =  build_states(spacial_dims,Nₘ)
    N=length(statesₙ)
    M=length(statesₘ)
    mat=zeros(Float64,M,N)
    for ket in 1:N
        state_ket=statesₙ[ket]
        state_bra , co = SC_create_operatr(state_ket,index_i,index_j)
        bra = findall(x->x==packbits(state_bra[:]),state_numₘ)
        mat[bra,ket].+=co
     end
     return mat
 end

function SC_annihilate_operatr(state::Array{Int8,DIMS2} ,index_i::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}
    #sc operator c^†_{index_i,↓} c^_{index_i,↑} where   index_i=(x_i,y_i,species_i)
    up_s   = CartesianIndex(index_i,1)
    down_s = CartesianIndex(index_i,2)
    ex2=0
    state_temp , ex1 = annihilate(state,up_s)
    ex1 == 0 ? state_f = state_temp : (state_f , ex2) = annihilate(state_temp, down_s)
    return state_f , ex1*ex2
end

function SC_annihilate_operatr(state::Array{Int8,DIMS2} ,index_i::CartesianIndex{DIMS},index_j::CartesianIndex{DIMS}) where {DIMS} where {DIMS2}
    #sc operator c_{index_i,↓} c^_{index_j,↑} where   index_i=(x_i,y_i,species_i)
    up_s   = CartesianIndex(index_j,1)
    down_s = CartesianIndex(index_i,2)
    ex2=0
    state_temp , ex1 = annihilate(state,up_s)
    ex1 == 0 ? state_f = state_temp : (state_f , ex2) = annihilate(state_temp, down_s)
    return state_f , ex1*ex2
end

function SC_annihilate_sub(spacial_dims::NTuple{DIMS2,Int64},Nₙ, Nₘ, index_i::CartesianIndex{DIMS} ) where {DIMS} where {DIMS2}

    statesₙ ,state_numₙ =  build_states(spacial_dims,Nₙ)
    statesₘ ,state_numₘ =  build_states(spacial_dims,Nₘ)
    N=length(statesₙ)
    M=length(statesₘ)
    mat=zeros(Float64,M,N)#spzeros(Float64,N,N)
    for ket in 1:N
        state_ket=statesₙ[ket]
        state_bra , co = SC_annihilate_operatr(state_ket,index_i)
        bra = findall(x->x==packbits(state_bra[:]),state_numₘ)
        mat[bra,ket].+=co
     end
     return mat
 end

 function SC_annihilate_sub(spacial_dims::NTuple{DIMS2,Int64},Nₙ, Nₘ, index_i::CartesianIndex{DIMS}, index_j::CartesianIndex{DIMS} ) where {DIMS} where {DIMS2}

    statesₙ ,state_numₙ =  build_states(spacial_dims,Nₙ)
    statesₘ ,state_numₘ =  build_states(spacial_dims,Nₘ)
    N=length(statesₙ)
    M=length(statesₘ)
    mat=zeros(Float64,M,N)#spzeros(Float64,N,N)
    for ket in 1:N
        state_ket=statesₙ[ket]
        state_bra , co = SC_annihilate_operatr(state_ket,index_i,index_j)
        bra = findall(x->x==packbits(state_bra[:]),state_numₘ)
        mat[bra,ket].+=co
     end
     return mat
 end

function hop(index::CartesianIndex{DIMS},dir::Int64,lr::Int64,dims::NTuple{DIMS,Int64}) where {DIMS}
    # update index
    if (lr==1)
      hop_index= index[dir]==dims[dir] ? 1 : index[dir]+1
    else
      hop_index= index[dir]==1 ? dims[dir] : index[dir]-1
    end
    # generate a new CartesianIndex with updated index
    CartesianIndex(Base.setindex(Tuple(index), hop_index, dir))

end

function build_states(dims::NTuple{DIMS,Int64}) where {DIMS}#(N,L_x,L_y,num_species,num_spin)
    # returen all states in the Fock space,  each state is reshaped to a an array of dims=(L_x,L_y,num_spin,num_species)
    Nₙ_n=prod(dims)
    N=2^Nₙ_n
    states=[];
    for c in 0:N-1
        state = dec2bin(c,Nₙ_n)
        push!(states,reshape(state , dims...))
    end
    return states
end

function build_states(dims::NTuple{DIMS,Int64},Nₙ::NTuple{DIMS2,Int64}) where {DIMS} where {DIMS2}
    #Nₙ is a vector of  Nₙber operators for the different types of electrons (classified by their specie and spin)
    #dims is a vector of the spacial dimensions.
    # returen all states in the sub-space of the Fock space
    Nₙ_sites=prod(dims)
    states=[] # TO DO -set type to integer
    for Nsp in Nₙ
        st = [ones(Int8,Nsp);zeros(Int8,Nₙ_sites-Nsp)]
        bit_combs=unique(permutations(st) )
        Nₙ_comb = Int64[]
        for bit_comb in bit_combs
            push!(Nₙ_comb,packbits(bit_comb)  )
        end
        push!(states,Nₙ_comb)
    end
    states_comb = collect(Iterators.product(states... ))[:]
    states_bin=[]
    for sc in states_comb
        push!(states_bin,bin_states(sc,dims,Int8(length(Nₙ)/2)))
    end
    state_num=Int64[]
    N=length(states_bin)
    for c in 1:N#state in states
     state_c = states_bin[c]
         push!(state_num,packbits(state_c[:]))
    end
       return states_bin , state_num
end

function bin_states(state::NTuple{DIMS2,Int64},dims::NTuple{DIMS,Int64},Nspecies::Int8) where {DIMS} where {DIMS2}
    Nspin=2
     bs=zeros(Int8,dims...,Nspecies,Nspin)
     for ns in 1:Nspecies
         for np in 1:Nspin
            bs[ntuple(k->:,length(dims))...,ns,np]=reshape(bitarray(state[ns+(np-1)*Nspecies],prod(dims)),dims)
         end
    end
    return bs
end

#=
function SC_corr_dw_tau_obs(spacial_dims::NTuple{DIMS2,Int64},num_species::Int64,rᵢ::CartesianIndex{DIMS},rⱼ::CartesianIndex{DIMS}
    , Nₙ , Dₙ::Vector{Float64} , Uₙ::Matrix{Float64} , Nₘ , Dₘ::Vector{Float64} , Uₘ::Matrix{Float64}
    , β::Float64,taus::Vector{Float64} ) where {DIMS} where {DIMS2}

    dims = (spacial_dims... ,num_species)
    corr = zeros(Float64,length(taus))
    for c in 1:num_species
        for g in 1:2
            if g==1
                index_i_sp=CartesianIndex(rᵢ,c) #(index_i...,c)#
                index_j_sp=CartesianIndex(rⱼ,c) #(index_j...,c)#
            else
                index_i_sp=CartesianIndex(rⱼ,c) #for the hemition conjugate of ΔᵢΔⱼ^†
                index_j_sp=CartesianIndex(rᵢ,c)
            end
            Δᵢ⁺=[]
            Δⱼ=[]
           # Δᵢ⁺d=[]
           # Δⱼd=[]
            η=0.0
            for dir in 1:length(spacial_dims)
                for lr in 1:2
                    index_i = index_i_sp
                    index_iₓ = hop(index_i,dir,lr,dims)
                    index_j = index_j_sp
                    index_jₓ = hop(index_j,dir,lr,dims)
                    lr==1 ? η=0.25 : η=-0.25
                    if  (dir==1 && lr==1) #!isempty(Δᵢ⁺)
                        Δᵢ⁺  = η.*SC_create_sub(spacial_dims,Nₙ,Nₘ, index_i,index_iₓ)# c^†_{index_i,↑} c^†_{index_iₓ,↓}
                        Δᵢ⁺ += η.* SC_create_sub(spacial_dims,Nₙ,Nₘ, index_iₓ,index_i)# c^†_{index_iₓ,↑} c^†_{index_i,↓}
                        Δⱼ   = η.*SC_annihilate_sub(spacial_dims,Nₘ,Nₙ, index_j,index_jₓ) # c_{index_j,↓} c^_{index_jₓ,↑}
                        Δⱼ  += η.*SC_annihilate_sub(spacial_dims,Nₘ,Nₙ, index_jₓ,index_j) # c_{index_jₓ,↓} c^_{index_j,↑}
                    else
                        Δᵢ⁺ += η.*SC_create_sub(spacial_dims,Nₙ,Nₘ, index_i,index_iₓ)# c^†_{index_i,↑} c^†_{index_iₓ,↓}
                        Δᵢ⁺ += η.* SC_create_sub(spacial_dims,Nₙ,Nₘ, index_iₓ,index_i)# c^†_{index_iₓ,↑} c^†_{index_i,↓}
                        Δⱼ  += η.*SC_annihilate_sub(spacial_dims,Nₘ,Nₙ, index_j,index_jₓ) # c_{index_j,↓} c^_{index_jₓ,↑}
                        Δⱼ  += η.*SC_annihilate_sub(spacial_dims,Nₘ,Nₙ, index_jₓ,index_j) # c_{index_jₓ,↓} c^_{index_j,↑}
                    end
                end
            end

            Δᵢ⁺d = Uₘ' *  Δᵢ⁺ * Uₙ
            Δⱼd =  Uₙ' *  Δⱼ  * Uₘ

            for i in 1:length(taus)
                τ=taus[i]
                e⁻⁽ᵝ⁻ᵀ⁾ᴰⁿ = Diagonal(exp.(-(β-τ).*Dₙ))
                e⁻ᵀᴰᵐ =     Diagonal(exp.(-τ*Dₘ))
                corr[i] += tr(e⁻⁽ᵝ⁻ᵀ⁾ᴰⁿ * Δⱼd * e⁻ᵀᴰᵐ * Δᵢ⁺d)
            end

        end
    end
    return corr
end
=#


text:
using Revise, ApproxOperator, LaTeXStrings, GLMakie

elements, nodes = importmsh("./msh/patchtest.msh")
nₚ = length(nodes[:x])

type = (Node,:Quadratic2D,:□,:CubicSpline)
elements["Ω"] = ReproducingKernel{type...,:Tri3}(elements["Ω"])
s = 2.5*0.25*ones(nₚ)
push!(nodes,:s₁=>s,:s₂=>s,:s₃=>s)

set𝓖!(elements["Ω"],:TriGI3,:∂1)

# n = 12
n = 21
inte = 100
# as = collect(0.0:1/inte:1.0)
as = LinRange(0.0,1.0,inte)

𝗠 = elements["Ω"][1].𝗠
𝝭 = elements["Ω"][1].𝝭
ap = ReproducingKernel{type...,:Node}([Node(i,nodes) for i in 1:nₚ],Node[],𝗠,𝝭)

# x = [x for x in as for y in as]
# y = [y for x in as for y in as]
z = [get𝝭(ap,(x,y,0.0),[n])[1] for x in as, y in as]

# plot(x,y,z,
#     st=:surface,
#     c=:jet1,
#     xlims=(0.0,1.0),
#     ylims=(0.0,1.0),
#     zlims=(-0.1,1.0),
#     xlabel=L"x",
#     ylabel=L"y",
#     zlabel=L"\Psi_I",
#     ztickfontrotation=10.0,
#     framestyle=:box,
#     legend=:none,
#     aspect_ratio=0.8,
#     dpi=300
#     )
f = Figure(fontsize=14)
ax1 = Axis3(f[1,1],
      xlabel="",
      ylabel="",
      zlabel="",
      limits=(0.0,1.0,0.0,1.0,-0.1,1.0),
      aspect=(1,1,0.7)
)
# ax2 = Axis3(f[1,1],limits=(0.0,1.0,0.0,1.0,-0.1,1.0))
surface!(as,as,z,colormap = :jet1)

meshscatter!([nodes[:x][n]], [nodes[:y][n]], [nodes[:z][n]],markersize = 0.05, color=:red)

save("shape1.png",f)
# ax = Axis3(figure[1,1],xlabel="x")

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
"""
    servelazy(project, dir)

Serve `project` lazily, meaning documents will only be built when they're
requested.
"""
function servelazy(
        project::Project,
        dir;
        dstdir = Path(mktempdir()),
        builder = FileBuilder(HTML(), dstdir))

    foreach(reset!, project.rewriters)
    foreach(k -> delete!(project.outputs, k), keys(project.outputs))

    # Do an empty build once for assets etc.
    build(project, builder, Dict())

    server = LazyServer(project, dir, builder, getwatcher(project, dir, builder))
    try
        LiveServer.serve(callback = server, dir = string(dstdir))
    catch e
        rethrow(e)
    finally
    end
end

mutable struct LazyServer
    project::Project
    dir::AbstractPath
    builder::FileBuilder
    watcher::LiveServer.SimpleWatcher
end

function (server::LazyServer)(req::HTTP.Request)
    p = Path(req.target[2:end])
    return server(p)
end

function (server::LazyServer)(path::AbstractPath)
    # Check if document is not built but is physical file
    buildpath = joinpath(server.builder.dir, path)
    if isfile(buildpath)
        return
    elseif extension(path) == "html"
        docpath = joinpath(parent(path), filename(path))
        srcpath = joinpath(server.dir, docpath)
        if isfile(srcpath)
            @info "Building $docpath for the first time..."
            # Add it to sources and rebuild the project
            dirtypaths = addfiles!(server.project, Dict(docpath => Pollen.parse(srcpath)))
            LiveServer.stop(server.watcher)
            server.watcher = getwatcher(server.project, server.dir, server.builder)
            build(server.builder, server.project, dirtypaths)
            @assert isfile(buildpath)
            return String(read(buildpath))
        end
    end
end


text:
function FiniteDifferences.to_vec(C::Cholesky)
    C_vec, factors_from_vec = to_vec(C.factors)
    function cholesky_from_vec(v)
        return Cholesky(factors_from_vec(v), C.uplo, C.info)
    end
    return C_vec, cholesky_from_vec
end

function FiniteDifferences.to_vec(x::Val)
    Val_from_vec(v) = x
    return Bool[], Val_from_vec
end

@testset "Factorizations" begin
    @testset "svd" begin
        for n in [4, 6, 10], m in [3, 5, 10]
            X = randn(n, m)
            F, dX_pullback = rrule(svd, X)
            for p in [:U, :S, :V]
                Y, dF_pullback = rrule(getproperty, F, p)
                Ȳ = randn(size(Y)...)

                dself1, dF, dp = dF_pullback(Ȳ)
                @test dself1 === NO_FIELDS
                @test dp === DoesNotExist()

                dself2, dX = dX_pullback(dF)
                @test dself2 === NO_FIELDS
                X̄_ad = unthunk(dX)
                X̄_fd = only(j′vp(central_fdm(5, 1), X->getproperty(svd(X), p), Ȳ, X))
                @test all(isapprox.(X̄_ad, X̄_fd; rtol=1e-6, atol=1e-6))
            end
            @testset "Vt" begin
                Y, dF_pullback = rrule(getproperty, F, :Vt)
                Ȳ = randn(size(Y)...)
                @test_throws ArgumentError dF_pullback(Ȳ)
            end
        end

        @testset "Thunked inputs" begin
            X = randn(4, 3)
            F, dX_pullback = rrule(svd, X)
            for p in [:U, :S, :V]
                Y, dF_pullback = rrule(getproperty, F, p)
                Ȳ = randn(size(Y)...)

                _, dF_unthunked, _ = dF_pullback(Ȳ)

                # helper to let us check how things are stored.
                backing_field(c, p) = getproperty(ChainRulesCore.backing(c), p)
                @assert !(backing_field(dF_unthunked, p) isa AbstractThunk)

                dF_thunked = map(f->Thunk(()->f), dF_unthunked)
                @assert backing_field(dF_thunked, p) isa AbstractThunk

                dself_thunked, dX_thunked = dX_pullback(dF_thunked)
                dself_unthunked, dX_unthunked = dX_pullback(dF_unthunked)
                @test dself_thunked == dself_unthunked
                @test dX_thunked == dX_unthunked
            end
        end

        @testset "+" begin
            X = [1.0 2.0; 3.0 4.0; 5.0 6.0]
            F, dX_pullback = rrule(svd, X)
            X̄ = Composite{typeof(F)}(U=zeros(3, 2), S=zeros(2), V=zeros(2, 2))
            for p in [:U, :S, :V]
                Y, dF_pullback = rrule(getproperty, F, p)
                Ȳ = ones(size(Y)...)
                dself, dF, dp = dF_pullback(Ȳ)
                @test dself === NO_FIELDS
                @test dp === DoesNotExist()
                X̄ += dF
            end
            @test X̄.U ≈ ones(3, 2) atol=1e-6
            @test X̄.S ≈ ones(2) atol=1e-6
            @test X̄.V ≈ ones(2, 2) atol=1e-6
        end

        @testset "Helper functions" begin
            X = randn(10, 10)
            Y = randn(10, 10)
            @test ChainRules._mulsubtrans!!(copy(X), Y) ≈ Y .* (X - X')
            @test ChainRules._eyesubx!(copy(X)) ≈ I - X
        end
    end

    @testset "eigendecomposition" begin
        @testset "eigen/eigen!" begin
            # NOTE: eigen!/eigen are not type-stable, so neither are their frule/rrule

            # avoid implementing to_vec(::Eigen)
            asnt(E::Eigen) = (values=E.values, vectors=E.vectors)

            # NOTE: for unstructured matrices, low enough n, and this specific seed, finite
            # differences of eigen seems to be stable enough for direct comparison.
            # This allows us to directly check differential of normalization/phase
            # convention
            n = 10

            @testset "eigen!(::Matrix{$T}) frule" for T in (Float64,ComplexF64)
                X = randn(T, n, n)
                Ẋ = rand_tangent(X)
                F = eigen!(copy(X))
                F_fwd, Ḟ_ad = frule((Zero(), copy(Ẋ)), eigen!, copy(X))
                @test F_fwd == F
                @test Ḟ_ad isa Composite{typeof(F)}
                Ḟ_fd = jvp(_fdm, asnt ∘ eigen! ∘ copy, (X, Ẋ))
                @test Ḟ_ad.values ≈ Ḟ_fd.values
                @test Ḟ_ad.vectors ≈ Ḟ_fd.vectors
                @test frule((Zero(), Zero()), eigen!, copy(X)) == (F, Zero())

                @testset "tangents are real when outputs are" begin
                    # hermitian matrices have real eigenvalues and, when real, real eigenvectors
                    X = Matrix(Hermitian(randn(T, n, n)))
                    Ẋ = Matrix(Hermitian(rand_tangent(X)))
                    _, Ḟ = frule((Zero(), Ẋ), eigen!, X)
                    @test eltype(Ḟ.values) <: Real
                    T <: Real && @test eltype(Ḟ.vectors) <: Real
                end
            end

            @testset "eigen(::Matrix{$T}) rrule" for T in (Float64,ComplexF64)
                # NOTE: eigen is not type-stable, so neither are is its rrule
                X = randn(T, n, n)
                F = eigen(X)
                V̄ = rand_tangent(F.vectors)
                λ̄ = rand_tangent(F.values)
                CT = Composite{typeof(F)}
                F_rev, back = rrule(eigen, X)
                @test F_rev == F
                _, X̄_values_ad = @inferred back(CT(values = λ̄))
                @test X̄_values_ad ≈ j′vp(_fdm, x -> eigen(x).values, λ̄, X)[1]
                _, X̄_vectors_ad = @inferred back(CT(vectors = V̄))
                @test X̄_vectors_ad ≈ j′vp(_fdm, x -> eigen(x).vectors, V̄, X)[1]
                F̄ = CT(values = λ̄, vectors = V̄)
                s̄elf, X̄_ad = @inferred back(F̄)
                @test s̄elf === NO_FIELDS
                X̄_fd = j′vp(_fdm, asnt ∘ eigen, F̄, X)[1]
                @test X̄_ad ≈ X̄_fd
                @test @inferred(back(Zero())) === (NO_FIELDS, Zero())
                F̄zero = CT(values = Zero(), vectors = Zero())
                @test @inferred(back(F̄zero)) === (NO_FIELDS, Zero())

                T <: Real && @testset "cotangent is real when input is" begin
                    X = randn(T, n, n)
                    Ẋ = rand_tangent(X)

                    F = eigen(X)
                    V̄ = rand_tangent(F.vectors)
                    λ̄ = rand_tangent(F.values)
                    F̄ = Composite{typeof(F)}(values = λ̄, vectors = V̄)
                    X̄ = rrule(eigen, X)[2](F̄)[2]
                    @test eltype(X̄) <: Real
                end
            end

            @testset "normalization/phase functions are idempotent" for T in (Float64,ComplexF64)
                # this is as much a math check as a code check. because normalization when
                # applied repeatedly is idempotent, repeated pushforward/pullback should
                # leave the (co)tangent unchanged
                X = randn(T, n, n)
                Ẋ = rand_tangent(X)
                F = eigen(X)

                V̇ = rand_tangent(F.vectors)
                V̇proj = ChainRules._eigen_norm_phase_fwd!(copy(V̇), X, F.vectors)
                @test !isapprox(V̇, V̇proj)
                V̇proj2 = ChainRules._eigen_norm_phase_fwd!(copy(V̇proj), X, F.vectors)
                @test V̇proj2 ≈ V̇proj

                V̄ = rand_tangent(F.vectors)
                V̄proj = ChainRules._eigen_norm_phase_rev!(copy(V̄), X, F.vectors)
                @test !isapprox(V̄, V̄proj)
                V̄proj2 = ChainRules._eigen_norm_phase_rev!(copy(V̄proj), X, F.vectors)
                @test V̄proj2 ≈ V̄proj
            end

            # below tests adapted from /test/rulesets/LinearAlgebra/symmetric.jl
            @testset "hermitian matrices" begin
                function _eigvecs_stabilize_mat(vectors)
                    Ui = @view(vectors[end, :])
                    return Diagonal(conj.(sign.(Ui)))
                end

                function _eigen_stable(A)
                    F = eigen(A)
                    rmul!(F.vectors, _eigvecs_stabilize_mat(F.vectors))
                    return F
                end

                n = 10
                @testset "eigen!(::Matrix{$T})" for T in (Float64, ComplexF64)
                    A, ΔA = Matrix(Hermitian(randn(T, n, n))), Matrix(Hermitian(randn(T, n, n)))

                    F = eigen!(copy(A))
                    @test frule((Zero(), Zero()), eigen!, copy(A)) == (F, Zero())
                    F_ad, ∂F_ad = frule((Zero(), copy(ΔA)), eigen!, copy(A))
                    @test F_ad == F
                    @test ∂F_ad isa Composite{typeof(F)}
                    @test ∂F_ad.values isa typeof(F.values)
                    @test ∂F_ad.vectors isa typeof(F.vectors)

                    f = x -> asnt(eigen(Matrix(Hermitian(x))))
                    ∂F_fd = jvp(_fdm, f, (A, ΔA))
                    @test ∂F_ad.values ≈ ∂F_fd.values

                    f_stable = x -> asnt(_eigen_stable(Matrix(Hermitian(x))))
                    F_stable = f_stable(A)
                    ∂F_stable_fd = jvp(_fdm, f_stable, (A, ΔA))
                    C = _eigvecs_stabilize_mat(F.vectors)
                    @test ∂F_ad.vectors * C ≈ ∂F_stable_fd.vectors
                end

                @testset "eigen(::Matrix{$T})" for T in (Float64, ComplexF64)
                    A, ΔU, Δλ = Matrix(Hermitian(randn(T, n, n))), randn(T, n, n), randn(n)

                    F = eigen(A)
                    ΔF = Composite{typeof(F)}(; values=Δλ, vectors=ΔU)
                    F_ad, back = rrule(eigen, A)
                    @test F_ad == F

                    C = _eigvecs_stabilize_mat(F.vectors)
                    CT = Composite{typeof(F)}

                    @testset for nzprops in ([:values], [:vectors], [:values, :vectors])
                        ∂F = CT(; [s => getproperty(ΔF, s) for s in nzprops]...)
                        ∂F_stable = (; [s => copy(getproperty(ΔF, s)) for s in nzprops]...)
                        :vectors in nzprops && rmul!(∂F_stable.vectors, C)

                        f_stable = function(x)
                            F_ = _eigen_stable(Matrix(Hermitian(x)))
                            return (; (s => getproperty(F_, s) for s in nzprops)...)
                        end

                        ∂self, ∂A = @inferred back(∂F)
                        @test ∂self === NO_FIELDS
                        @test ∂A isa typeof(A)
                        ∂A_fd = j′vp(_fdm, f_stable, ∂F_stable, A)[1]
                        @test ∂A ≈ ∂A_fd
                    end
                end
            end
        end

        @testset "eigvals/eigvals!" begin
            # NOTE: eigvals!/eigvals are not type-stable, so neither are their frule/rrule
            @testset "eigvals!(::Matrix{$T}) frule" for T in (Float64,ComplexF64)
                n = 10
                X = randn(T, n, n)
                λ = eigvals!(copy(X))
                Ẋ = rand_tangent(X)
                frule_test(eigvals!, (X, Ẋ))
                @test frule((Zero(), Zero()), eigvals!, copy(X)) == (λ, Zero())

                @testset "tangents are real when outputs are" begin
                    # hermitian matrices have real eigenvalues
                    X = Matrix(Hermitian(randn(T, n, n)))
                    Ẋ = Matrix(Hermitian(rand_tangent(X)))
                    _, λ̇ = frule((Zero(), Ẋ), eigvals!, X)
                    @test eltype(λ̇) <: Real
                end
            end

            @testset "eigvals(::Matrix{$T}) rrule" for T in (Float64,ComplexF64)
                n = 10
                X = randn(T, n, n)
                X̄ = rand_tangent(X)
                λ̄ = rand_tangent(eigvals(X))
                rrule_test(eigvals, λ̄, (X, X̄))
                back = rrule(eigvals, X)[2]
                @inferred back(λ̄)
                @test @inferred(back(Zero())) === (NO_FIELDS, Zero())

                T <: Real && @testset "cotangent is real when input is" begin
                    X = randn(T, n, n)
                    λ = eigvals(X)
                    λ̄ = rand_tangent(λ)
                    X̄ = rrule(eigvals, X)[2](λ̄)[2]
                    @test eltype(X̄) <: Real
                end
            end

            # below tests adapted from /test/rulesets/LinearAlgebra/symmetric.jl
            @testset "hermitian matrices" begin
                n = 10
                @testset "eigvals!(::Matrix{$T})" for T in (Float64, ComplexF64)
                    A, ΔA = Matrix(Hermitian(randn(T, n, n))), Matrix(Hermitian(randn(T, n, n)))
                    λ = eigvals!(copy(A))
                    λ_ad, ∂λ_ad = frule((Zero(), copy(ΔA)), eigvals!, copy(A))
                    @test λ_ad ≈ λ # inexact because frule uses eigen not eigvals
                    @test ∂λ_ad isa typeof(λ)
                    @test ∂λ_ad ≈ jvp(_fdm, A -> eigvals(Matrix(Hermitian(A))), (A, ΔA))
                end

                @testset "eigvals(::Matrix{$T})" for T in (Float64, ComplexF64)
                    A, Δλ = Matrix(Hermitian(randn(T, n, n))), randn(n)
                    λ = eigvals(A)
                    λ_ad, back = rrule(eigvals, A)
                    @test λ_ad ≈ λ # inexact because rrule uses eigen not eigvals
                    ∂self, ∂A = @inferred back(Δλ)
                    @test ∂self === NO_FIELDS
                    @test ∂A isa typeof(A)
                    @test ∂A ≈ j′vp(_fdm, A -> eigvals(Matrix(Hermitian(A))), Δλ, A)[1]
                    @test @inferred(back(Zero())) == (NO_FIELDS, Zero())
                end
            end
        end
    end

    # These tests are generally a bit tricky to write because FiniteDifferences doesn't
    # have fantastic support for this stuff at the minute.
    @testset "cholesky" begin
        @testset "Real" begin
            C = cholesky(rand() + 0.1)
            ΔC = Composite{typeof(C)}((factors=rand_tangent(C.factors)))
            rrule_test(cholesky, ΔC, (rand() + 0.1, randn()))
        end
        @testset "Diagonal{<:Real}" begin
            D = Diagonal(rand(5) .+ 0.1)
            C = cholesky(D)
            ΔC = Composite{typeof(C)}((factors=Diagonal(randn(5))))
            rrule_test(cholesky, ΔC, (D, Diagonal(randn(5))), (Val(false), nothing))
        end

        X = generate_well_conditioned_matrix(10)
        V = generate_well_conditioned_matrix(10)
        F, dX_pullback = rrule(cholesky, X, Val(false))
        @testset "uplo=$p" for p in [:U, :L]
            Y, dF_pullback = rrule(getproperty, F, p)
            Ȳ = (p === :U ? UpperTriangular : LowerTriangular)(randn(size(Y)))
            (dself, dF, dp) = dF_pullback(Ȳ)
            @test dself === NO_FIELDS
            @test dp === DoesNotExist()

            # NOTE: We're doing Nabla-style testing here and avoiding using the `j′vp`
            # machinery from FiniteDifferences because that isn't set up to respect
            # necessary special properties of the input. In the case of the Cholesky
            # factorization, we need the input to be Hermitian.
            ΔF = unthunk(dF)
            _, dX = dX_pullback(ΔF)
            X̄_ad = dot(unthunk(dX), V)
            X̄_fd = central_fdm(5, 1)(0.000_001) do ε
                dot(Ȳ, getproperty(cholesky(X .+ ε .* V), p))
            end
            @test X̄_ad ≈ X̄_fd rtol=1e-4
        end

        # Ensure that cotangents of cholesky(::StridedMatrix) and
        # (cholesky ∘ Symmetric)(::StridedMatrix) are equal.
        @testset "Symmetric" begin
            X_symmetric, sym_back = rrule(Symmetric, X, :U)
            C, chol_back_sym = rrule(cholesky, X_symmetric, Val(false))

            Δ = Composite{typeof(C)}((U=UpperTriangular(randn(size(X)))))
            ΔX_symmetric = chol_back_sym(Δ)[2]
            @test sym_back(ΔX_symmetric)[2] ≈ dX_pullback(Δ)[2]
        end
    end
end


text:
function initialize_device(device::PSY.DynamicGenerator)
    #Obtain States
    device_states = zeros(PSY.get_n_states(device))

    #Initialize Machine and Shaft: δ and ω
    initialize_mach_shaft!(device_states, device)
    #Initialize extra Shaft states
    initialize_shaft!(device_states, device)
    #Initialize AVR
    initialize_avr!(device_states, device)
    #Initialize TG
    initialize_tg!(device_states, device)
    #Initialize PSS
    initialize_pss!(device_states, device)

    return device_states
end

function initialize_device(device::PSY.DynamicInverter)
    #Obtain States
    device_states = zeros(PSY.get_n_states(device))

    #Initialize Machine and Shaft: V and I
    initialize_filter!(device_states, device)
    #Initialize freq estimator
    initialize_frequency_estimator!(device_states, device)
    #Initialize OuterLoop
    initialize_outer!(device_states, device)
    #Initialize DCside
    initialize_DCside!(device_states, device)
    #Initialize InnerLoop
    initialize_inner!(device_states, device)
    #Initialize Converter
    initialize_converter!(device_states, device)
    return device_states
end

function initialize_device(device::PSY.Source)
    #PowerFlow Data
    P0 = PSY.get_active_power(device)
    Q0 = PSY.get_reactive_power(device)
    Vm = PSY.get_magnitude(PSY.get_bus(device))
    θ = PSY.get_angle(PSY.get_bus(device))
    S0 = P0 + Q0 * 1im
    V_R = Vm * cos(θ)
    V_I = Vm * sin(θ)
    V = V_R + V_I * 1im
    I = conj(S0 / V)
    I_R = real(I)
    I_I = imag(I)
    R_th = PSY.get_R_th(device)
    X_th = PSY.get_X_th(device)
    Zmag = R_th^2 + X_th^2

    function f!(out, x)
        V_R_internal = x[1]
        V_I_internal = x[2]

        out[1] =
            R_th * (V_R_internal - V_R) / Zmag + X_th * (V_I_internal - V_I) / Zmag - I_R
        out[2] =
            R_th * (V_I_internal - V_I) / Zmag - X_th * (V_R_internal - V_R) / Zmag - I_I
    end
    x0 = [V_R, V_I]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in Source failed")
    else
        sol_x0 = sol.zero
        #Update terminal voltages
        V_internal = sqrt(sol_x0[1]^2 + sol_x0[2]^2)
        θ_internal = angle(sol_x0[1] + sol_x0[2] * 1im)
        PSY.set_internal_voltage!(device, V_internal)
        PSY.set_internal_angle!(device, θ_internal)
    end
end

function initialize_device(branch::PSY.DynamicBranch)
    device_states = zeros(PSY.get_n_states(branch))
    #PowerFlow Data
    arc = PSY.get_arc(branch)
    Vm_from = PSY.get_magnitude(PSY.get_from(arc))
    θ_from = PSY.get_angle(PSY.get_from(arc))
    Vm_to = PSY.get_magnitude(PSY.get_to(arc))
    θ_to = PSY.get_angle(PSY.get_to(arc))
    V_R_from = Vm_from * cos(θ_from)
    V_I_from = Vm_from * sin(θ_from)
    V_R_to = Vm_to * cos(θ_to)
    V_I_to = Vm_to * sin(θ_to)
    R = PSY.get_r(branch)
    X = PSY.get_x(branch)
    Zmag = R^2 + X^2
    #Compute Current
    I_R = R * (V_R_from - V_R_to) / Zmag + X * (V_I_from - V_I_to) / Zmag
    I_I = R * (V_I_from - V_I_to) / Zmag - X * (V_R_from - V_R_to) / Zmag
    #Update Current
    device_states[1] = I_R
    device_states[2] = I_I
    return device_states
end


text:
axisTypes=Union{Number,Date,AbstractString}
struct DataAttrs
    data::Vector{T} where T<:axisTypes
end

mutable struct Series
   plot_type::String
   name::String
   data::Vector{DataAttrs}
   options::Dict{String,Any}
end

mutable struct Plot
   series::Vector{Series}
   options::Dict{String,Any}
end


text:
#=
function SetNextPlotTicksX(values::Vector{<:Real}, n_ticks::Integer, args...)
    eltype(values) !== Float64 && (values = Float64.(values))
    return SetNextPlotTicksX(values, n_ticks, args...)
end

function SetNextPlotTicksY(values::Vector{<:Real}, n_ticks::Integer, args...)
    eltype(values) !== Float64 && (values = Float64.(values))
    return SetNextPlotTicksY(values, n_ticks, args...)
end
=#


text:
using RoME
using Test

##

@testset "test DynPose2 and velocity..." begin

##

N = 100
fg = initfg()

# add first pose locations
addVariable!(fg, :x0, DynPose2; nanosecondtime=0)

# Prior factor as boundary condition
pp0 = DynPose2VelocityPrior(MvNormal(zeros(3), Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([10.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x0;], pp0)

# initialize the first pose
IncrementalInference.doautoinit!(fg, [getVariable(fg,:x0);])

addVariable!(fg, :x1, DynPose2;  nanosecondtime=1000_000_000)

# conditional likelihood between Dynamic Point2
dp2dp2 = VelPose2VelPose2(MvNormal([10.0;0;0], Matrix(Diagonal([0.01;0.01;0.001].^2))),
                          MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x0;:x1], dp2dp2)

pts = approxConv(fg, :x0x1f1, :x1)

# ensureAllInitialized!(fg)

tree, smt, hist = solveTree!(fg);

X1 = getVal(fg, :x1)

@test 0.9*N <= sum(abs.(X1[1,:] .- 10.0) .< 0.75)
@test 0.9*N <= sum(abs.(X1[2,:] .- 0.0) .< 0.75)
# @show TU.wrapRad.(X1[3,:])
@test 0.8*N <= sum(abs.(TU.wrapRad.(X1[3,:]) .- 0.0) .< 0.25)
# @warn "wrapRad issue, accepting 80% as good enough until issue JuliaRobotics/RoME.jl#90 is fixed."
@test 0.9*N <= sum(abs.(X1[4,:] .- 10.0) .< 0.5)
@test 0.9*N <= sum(abs.(X1[5,:] .- 0.0) .< 0.5)

# using RoMEPlotting
# # plotLocalProduct(fg, :x10, dims=[1;2])
# # plotSLAM2DPoses(fg)
# xx1 = marginal(getKDE(fg, :x1),[1;2;3])
# plotPose(Pose2(), [xx1])
# plotPose(fg, [:x1], levels=1, show=false)
#
# plotKDE(marginal(getKDE(fg, :x1),[4;5]), levels=5)

##

end

@testset "test distribution compare functions..." begin

##

mu = randn(6)
mv1 = MvNormal(deepcopy(mu), Matrix{Float64}(LinearAlgebra.I, 6,6))
mv2 = MvNormal(deepcopy(mu), Matrix{Float64}(LinearAlgebra.I, 6,6))
mv3 = MvNormal(randn(6), Matrix{Float64}(LinearAlgebra.I, 6,6))
@test RoME.compareDensity(mv1, mv2)
@test !RoME.compareDensity(mv1, mv3)
@test !RoME.compareDensity(mv2, mv3)

##

end

@testset "test DynPose2 packing converters..." begin

##

pp0 = DynPose2VelocityPrior(MvNormal(zeros(3), Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([10.0;0], Matrix(Diagonal([0.1; 0.1].^2))))

pp = convert(PackedDynPose2VelocityPrior, pp0)
ppu = convert(DynPose2VelocityPrior, pp)

@test RoME.compare(pp0, ppu)

dp2dp2 = VelPose2VelPose2(MvNormal([10.0;0;0], Matrix(Diagonal([0.01;0.01;0.001].^2))),
                          MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))

pp = convert(PackedVelPose2VelPose2, dp2dp2)
ppu = convert(VelPose2VelPose2, pp)

@test RoME.compare(dp2dp2, ppu)

##

end

@testset "test many DynPose2 chain stationary and 'pulled'..." begin

##

N = 100
fg = initfg()

# add first pose locations
addVariable!(fg, :x0, DynPose2; nanosecondtime=0)

# Prior factor as boundary condition
pp0 = DynPose2VelocityPrior(MvNormal(zeros(3), Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x0;], pp0)

sym = :x0
k = 0
for sy in Symbol[Symbol("x$i") for i in 1:10]

k+=1
addVariable!(fg, sy, DynPose2; nanosecondtime=1000_000_000*k)

# conditional likelihood between Dynamic Point2
dp2dp2 = VelPose2VelPose2(MvNormal([0.0;0;0], Matrix(Diagonal([1.0;0.1;0.001].^2))),
                          MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [sym;sy], dp2dp2)
sym =sy

end # for

x5 = KDE.getKDEMean(getKDE(getVariable(fg, :x5)))

@test abs(x5[1]) < 1.25
@test abs(x5[2]) < 1.25
@test abs(TU.wrapRad(x5[3])) < 0.4
@test abs(x5[4]) < 0.5
@test abs(x5[5]) < 0.5

ensureAllInitialized!(fg)

x10 = KDE.getKDEMean(getKDE(getVariable(fg, :x10)))

@test abs(x10[1]) < 1.25
@test abs(x10[2]) < 1.25
@test abs(TU.wrapRad(x10[3])) < 0.4
@test abs(x10[4]) < 0.5
@test abs(x10[5]) < 0.5

# drawGraph(fg, show=true)
# tree = buildTreeReset!(fg)
# drawTree(tree, show=true)

# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm, 25cm)
# plotSLAM2DPoses(fg)
# plotPose(fg, [:x10])

##

@error ".useMsgLikelihoods = false required until IIF #1010 completed."
getSolverParams(fg).useMsgLikelihoods = false

# solve
smtasks = Task[]
tree, smt, hist = solveTree!(fg, smtasks=smtasks); #, recordcliqs=ls(fg));

##

# hists = IIF.fetchCliqHistoryAll!(smtasks)

# IIF.printCSMHistorySequential(hists)
# IIF.printCSMHistorySequential(hists, 3=>1:50)
# IIF.printCSMHistoryLogical(hists)

# ##

# getLogPath(fg)

# using Caesar, Images, Graphs

# tree = hists[1][1].csmc.tree

# csmAnimateSideBySide(tree, hists, encode=true, nvenc=true, show=true)

# ##

# sfg = hists[3][8].csmc.cliqSubFg
# drawGraph(sfg, show=true)
# getSolverParams(hists[3][8].csmc.cliqSubFg).dbg = true
# fnc_, csmc_ = IIF.repeatCSMStep!(hists, 3, 8);
# # sfg_= csmc_.cliqSubFg

# sfg = hists[3][9].csmc.cliqSubFg
# drawGraph(sfg, show=true)

# IIF.getMessageBuffer(hists[3][8].csmc.cliq).downRx.belief

# getSolverParams(hists[3][9].csmc.cliqSubFg).dbg = true
# IIF.repeatCSMStep!(hists, 3, 9);

# ##

# tree_ = hists[3][9].csmc.tree
# drawTree(tree_, show=true)

# ## whats up with down msgs

# sfg9 = loadDFG(joinpath(getLogPath(fg), "logs/cliq3", "fg_beforedownsolve.tar.gz"))
# drawGraph(sfg9, show=true)

# ls(sfg9, :x0)

# getFactor(sfg9, :x0f2)

# ##

# sfg = hists[7][6].csmc.cliqSubFg
# drawGraph(sfg, show=true)
# IIF.repeatCSMStep!(hists, 7, 6)

# getSolverParams(hists[4][4].csmc.cliqSubFg).dbg = true
# IIF.repeatCSMStep!(hists, 4, 4);

# IIF.getMessageBuffer(hists[4][4].csmc.cliq).upRx[7].belief
# IIF.getMessageBuffer(hists[4][4].csmc.cliq).upRx[7]
# IIF.getMessageBuffer(hists[4][4].csmc.cliq).upRx[7].jointmsg.priors

# IIF.repeatCSMStep!(hists, 4, 5)

##

# tree = buildTreeReset!(fg);
# drawTree(tree, show=true);

##

x5 = getPPE(getVariable(fg, :x5)).suggested
# x5 = KDE.getKDEMean(getKDE(getVariable(fg, :x5)))

@test abs(x5[1]) < 1.5
@test abs(x5[2]) < 1.5
@test abs(TU.wrapRad(x5[3])) < 0.4
@test abs(x5[4]) < 0.5
@test abs(x5[5]) < 0.5

x10 = KDE.getKDEMean(getKDE(getVariable(fg, :x10)))

@test abs(x10[1]) < 2.75
@test abs(x10[2]) < 2.75
@test abs(TU.wrapRad(x10[3])) < 0.5
@test abs(x10[4]) < 0.5
@test abs(x10[5]) < 0.5

# pull the tail end out with position
pp10 = DynPose2VelocityPrior(MvNormal([10.0;0;0], Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x10;], pp10)

fg2 = deepcopy(fg)

tree, mst, hist = solveTree!(fg) # N=N

x10 = KDE.getKDEMean(getKDE(getVariable(fg, :x10)))

@test 5.0 < x10[1]
@test abs(x10[2]) < 1.0
@test abs(TU.wrapRad(x10[3])) < 0.6
@test -0.1 < x10[4] < 1.0
@test abs(x10[5]) < 0.5

for sym in [Symbol("x$i") for i in 2:9]

XX = KDE.getKDEMean(getKDE(getVariable(fg, sym)))

@show sym, round.(XX,digits=5)
@test -2.0 < XX[1] < 10.0
@test abs(XX[2]) < 1.0
@test abs(TU.wrapRad(XX[3])) < 1.3
@test -0.5 < XX[4] < 2.0
@test abs(XX[5]) < 0.5

end

##

end

# using RoMEPlotting
# plotLocalProduct(fg, :x10, dims=[1;2])
# drawPoses(fg)
# plotPose(fg, [:x9],levels=1);
#
# plotPose(fg, [:x1;:x2;:x3;:x4;:x5;:x6;:x7;:x8;:x9;:x10],levels=1);

# savejld(fg) # tempfg.jld

@testset "test many DynPose2 sideways velocity..." begin

global N = 100
global fg = initfg()

# add first pose locations
addVariable!(fg, :x0, DynPose2; nanosecondtime=0)

# Prior factor as boundary condition
global pp0 = DynPose2VelocityPrior(MvNormal([0.0;0.0;pi/2], Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([0.0;0], Matrix(Diagonal([0.5; 0.5].^2))))
addFactor!(fg, [:x0;], pp0)

addVariable!(fg, :x1, DynPose2; nanosecondtime=1000_000_000)

global pp0 = DynPose2VelocityPrior(MvNormal([1.0;0.0;pi/2], Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([0.0;0], Matrix(Diagonal([0.5; 0.5].^2))))
addFactor!(fg, [:x1;], pp0)

# conditional likelihood between Dynamic Point2
global dp2dp2 = VelPose2VelPose2(MvNormal([0.0;-1.0;0], Matrix(Diagonal([0.01;0.01;0.001].^2))),
                          MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x0;:x1], dp2dp2)

getSolverParams(fg).N = N
solveTree!(fg)

# test for velocity in the body frame
global x0 = KDE.getKDEMean(getKDE(getVariable(fg, :x0)))

@test -0.4 < x0[1] < 2.0
@test abs(x0[2]) < 0.5
@test abs(x0[3] - pi/2) < 0.1
@test abs(x0[4]) < 0.4
@test -1.5 < x0[5] < -0.5

global x1 = KDE.getKDEMean(getKDE(getVariable(fg, :x1)))

@test -0.1 < x1[1] < 2.0
@test abs(x1[2]) < 0.5
@test abs(x1[3] - pi/2) < 0.1
@test abs(x1[4]) < 0.4
@test -1.5 < x1[5] < -0.5

end

## debugging 458============================

# using RoMEPlotting, Gadfly
# Gadfly.set_default_plot_size(35cm,25cm)
# plotPose(fg, :x10)
#
# drawGraph(fg2, show=true)
# drawTree(tree, show=true, imgs=true)
#
#
# sfg = buildCliqSubgraph(fg, tree, :x9)
# sfg = buildCliqSubgraph(fg, tree, :x10)
# drawGraph(sfg)
#
# # fg = deepcopy(fg2)
#
# getSolverParams(fg).dbg = true
# getSolverParams(fg).showtree = true
# getSolverParams(fg).drawtree = true
# getSolverParams(fg).multiproc = false
#
# ##============================================================================
#
# getLogPath(fg)
# tree, mst, hist = solveTree!(fg, recordcliqs=ls(fg))
#
#
#
# printCliqHistorySummary(tree,:x10)
#
# getClique(tree, :x10)
#
# csmc1 = hist[1][6][4]
# csfg = csmc1.cliqSubFg
# drawGraph(csfg)
#
# stuff = sandboxCliqResolveStep(tree, :x10, 6)
#
#
# getKDE(hist[1][6][4].cliqSubFg, :x10) |> getPoints
# getKDE(hist[1][7][4].cliqSubFg, :x10) |> getPoints
# getKDE(stuff[4].cliqSubFg, :x10) |> getPoints
#
# getKDE(fg, :x10) |> getPoints
#
# tree = wipeBuildNewTree!(fg)
#
# getData(getClique(tree, :x9 ))
# getData(getClique(tree, :x10))

#


text:
export OrcJIT, OrcModule, OrcTargetAddress
export dispose, errormsg, compile!, remove!, add!,
       mangle, address, addressin, create_stub!, set_stub!,
       register!, unregister!, callback!

@checked struct OrcJIT
    ref::API.LLVMOrcJITStackRef
end

Base.unsafe_convert(::Type{API.LLVMOrcJITStackRef}, orc::OrcJIT) = orc.ref
Base.unsafe_convert(::Type{Ptr{Cvoid}}, orc::OrcJIT) = Base.unsafe_convert(Ptr{Cvoid}, orc.ref)

OrcJIT(ref::Ptr{Cvoid}) = OrcJIT(Base.unsafe_convert(API.LLVMOrcJITStackRef, ref))

"""
    OrcJIT(::TargetMachine)

Creates a OrcJIT stack based on the provided target machine.

!!! warning
    Takes ownership of the provided target machine.
"""
function OrcJIT(tm::TargetMachine)
    OrcJIT(API.LLVMOrcCreateInstance(tm))
end

function dispose(orc::OrcJIT)
    API.LLVMOrcDisposeInstance(orc)
end

function OrcJIT(f::Core.Function, tm::TargetMachine)
    orc = OrcJIT(tm)
    try
        f(orc)
    finally
        dispose(orc)
    end
end

function errormsg(orc::OrcJIT)
    # The error message is owned by `orc`, and will
    # be disposed along-side the OrcJIT.
    unsafe_string(LLVM.API.LLVMOrcGetErrorMsg(orc))
end

struct OrcModule
    handle::API.LLVMOrcModuleHandle
end
Base.convert(::Type{API.LLVMOrcModuleHandle}, mod::OrcModule) = mod.handle

"""
   resolver(name, ctx)

Lookup the symbol `name`. Iff `ctx` is passed to this function it should be a
pointer to the OrcJIT we are compiling for.
"""
function resolver(name, ctx)
    name = unsafe_string(name)
    ## Step 0: Should have already resolved it iff it was in the
    ##         same module
    ## Step 1: See if it's something known to the execution engine
    ptr = C_NULL
    if ctx != C_NULL
        orc = OrcJIT(ctx)
        ptr = pointer(address(orc, name))
    end

    ## Step 2: Search the program symbols
    if ptr == C_NULL
        #
        # SearchForAddressOfSymbol expects an unmangled 'C' symbol name.
        # Iff we are on Darwin, strip the leading '_' off.
        @static if Sys.isapple()
            if name[1] == '_'
                name = name[2:end]
            end
        end
        ptr = LLVM.find_symbol(name)
    end

    ## Step 4: Lookup in libatomic
    # TODO: Do we need to do this?

    if ptr == C_NULL
        error("OrcJIT: Symbol `$name` lookup failed. Aborting!")
    end

    return UInt64(reinterpret(UInt, ptr))
end

function compile!(orc::OrcJIT, mod::Module, resolver = @cfunction(resolver, UInt64, (Cstring, Ptr{Cvoid})), resolver_ctx = orc; lazy=false)
    r_mod = Ref{API.LLVMOrcModuleHandle}()
    if lazy
        API.LLVMOrcAddLazilyCompiledIR(orc, r_mod, mod, resolver, resolver_ctx)
    else
        API.LLVMOrcAddEagerlyCompiledIR(orc, r_mod, mod, resolver, resolver_ctx)
    end
    OrcModule(r_mod[])
end

function Base.delete!(orc::OrcJIT, mod::OrcModule)
    LLVM.API.LLVMOrcRemoveModule(orc, mod)
end

function add!(orc::OrcJIT, obj::MemoryBuffer, resolver = @cfunction(resolver, UInt64, (Cstring, Ptr{Cvoid})), resolver_ctx = orc)
    r_mod = Ref{API.LLVMOrcModuleHandle}()
    API.LLVMOrcAddObjectFile(orc, r_mod, obj, resolver, resolver_ctx)
    return OrcModule(r_mod[])
end

function mangle(orc::OrcJIT, name)
    r_symbol = Ref{Cstring}()
    LLVM.API.LLVMOrcGetMangledSymbol(orc, r_symbol, name)
    symbol = unsafe_string(r_symbol[])
    LLVM.API.LLVMOrcDisposeMangledSymbol(r_symbol[])
    return symbol
end

struct OrcTargetAddress
    ptr::API.LLVMOrcTargetAddress
end
Base.convert(::Type{API.LLVMOrcTargetAddress}, addr::OrcTargetAddress) = addr.ptr

Base.pointer(addr::OrcTargetAddress) = reinterpret(Ptr{Cvoid}, addr.ptr % UInt) # LLVMOrcTargetAddress is UInt64 even on 32-bit

OrcTargetAddress(ptr::Ptr{Cvoid}) = OrcTargetAddress(reinterpret(UInt, ptr))

function create_stub!(orc::OrcJIT, name, initial)
    LLVM.API.LLVMOrcCreateIndirectStub(orc, name, initial)
end

function set_stub!(orc::OrcJIT, name, new)
    LLVM.API.LLVMOrcSetIndirectStubPointer(orc, name, new)
end

function address(orc::OrcJIT, name)
    r_address = Ref{API.LLVMOrcTargetAddress}()
    API.LLVMOrcGetSymbolAddress(orc, r_address, name)
    OrcTargetAddress(r_address[])
end

function addressin(orc::OrcJIT, mod::OrcModule, name)
    r_address = Ref{API.LLVMOrcTargetAddress}()
    API.LLVMOrcGetSymbolAddressIn(orc, r_address, mod, name)
    OrcTargetAddress(r_address[])
end

function callback!(orc::OrcJIT, callback, ctx)
    r_address = Ref{API.LLVMOrcTargetAddress}()
    API.LLVMOrcCreateLazyCompileCallback(orc, r_address, callback, ctx)
    return OrcTargetAddress(r_address[])
end

function register!(orc::OrcJIT, listener::JITEventListener)
    LLVM.API.LLVMOrcRegisterJITEventListener(orc, listener)
end

function unregister!(orc::OrcJIT, listener::JITEventListener)
    LLVM.API.LLVMOrcUnregisterJITEventListener(orc, listener)
end


text:
using StrongFieldApproximation
using ElectricFields
using Unitful
using UnitfulAtomic
using Test

# write your own tests here
@test 1 == 1
@testset "Telescope iterators" begin
    generate_reference_common(ref) =
        map(I -> [Tuple(I)...], findall(ref))

    generate_1d_reference(N, memory) = [[i] for i in 1:N]

    function generate_2d_reference(N, memory)
        ref = [0 < i - j ≤ memory
               for i in 1:N, j in 1:N]
        generate_reference_common(ref)
    end

    function generate_3d_reference(N, memory)
        ref = [i > j > k && 0 < i - j ≤ memory && 0 < j - k ≤ memory
               for i in 1:N, j in 1:N, k in 1:N]
        generate_reference_common(ref)
    end

    function generate_4d_reference(N, memory)
        ref = [i > j > k > l && 0 < i - j ≤ memory && 0 < j - k ≤ memory && 0 < k - l ≤ memory
               for i in 1:N, j in 1:N, k in 1:N, l in 1:N]
        generate_reference_common(ref)
    end

    function test_telescope_iterator(N, n, memory)
        v = 1:N
        ti = StrongFieldApproximation.TelescopeIterator(v, n, memory)
        # At the moment we don't know how to compute length(ti), so
        # collect(ti) does not work.
        data = Vector{Int}[]
        for e in ti
            push!(data, e)
        end
        ref = if n == 1
            generate_1d_reference(N, memory)
        elseif n == 2
            generate_2d_reference(N, memory)
        elseif n == 3
            generate_3d_reference(N, memory)
        elseif n == 4
            generate_4d_reference(N, memory)
        end

        @test data == ref
    end
    @testset "N = $N, n = $n, memory = $memory" for N = 1:7, n = 1:4, memory = [0,1,3,100]
        test_telescope_iterator(N, n, memory)
    end
end

@testset "Diagrams" begin
    function display_diagram(system, diagram, expected_momenta, expected_unique, expected_indeterminate)
        display(diagram)
        ions, unique_momenta, momenta, indeterminate_momenta, order = StrongFieldApproximation.analyze_diagram(system, diagram)

        expected_ions = first.(diagram.path)
        if first(diagram.path)[2] == 0 && length(diagram) > 1
            expected_ions = expected_ions[2:end]
        end

        @test ions == expected_ions
        @test momenta == expected_momenta
        @test unique_momenta == expected_unique
        @test indeterminate_momenta == expected_indeterminate
    end

    @field(F) do
        λ = 800.0u"nm"
        I₀ = 1e14u"W/cm^2"
        cycles = 4.0
        ϕ = π
        env = :cos²
    end

    ndt = 238

    Iₚ = 14u"eV" # "Krypton"

    # Elastic scattering off a Yukawa potential
    cc = StrongFieldApproximation.CoulombCoupling((𝐤,𝐩) -> yukawa_fourier(𝐩-𝐤, 1, 0, 1))
    dc = StrongFieldApproximation.DipoleCoupling(0.1, F)
    couplings=Matrix{StrongFieldApproximation.AbstractCoupling}[reshape([dc],1,1),reshape([cc],1,1)]

    system = StrongFieldApproximation.System(Iₚ, F, ndt, couplings=couplings)

    for (path, expected_momenta, expected_unique, expected_indeterminate) in [
        ([(1,0)], [1], [(1,2)], []),
        ([(1,0),(1,0)], [1], [(1,2)], 1:1),
        ([(2,1),(1,0)], [1,1], [(1,3)], []),
        ([(2,2),(1,0)], [1,2], [(1,2),(2,3)], 2:2),
        ([(1,2),(2,2),(1,0)], [1,2,3], [(1,2),(2,3),(3,4)], 2:3),
        ([(1,0),(1,2),(2,2),(1,0)], [1,2,3], [(1,2),(2,3),(3,4)], 1:3),
        ([(1,1),(2,2),(1,0)], [1,1,2], [(1,3),(3,4)], 2:2),
        ([(1,0),(1,1),(2,2),(1,0)], [1,1,2], [(1,3),(3,4)], 1:2),
        ([(3,1),(1,1),(2,2),(1,0)], [1,1,1,2], [(1,4),(4,5)], 2:2),
        ([(3,0),(3,1),(1,1),(2,2),(1,0)], [1,1,1,2], [(1,4),(4,5)], 1:2),
        ([(1,2),(3,1),(1,1),(2,2),(1,0)], [1,2,2,2,3], [(1,2),(2,5),(5,6)], 2:3),
        ([(1,2),(1,2),(3,1),(2,2),(1,0)], [1,2,3,3,4], [(1,2),(2,3),(3,5),(5,6)], 2:4),
        ([(1,2),(1,2),(3,1),(1,1),(2,2),(1,0)], [1,2,3,3,3,4], [(1,2),(2,3),(3,6),(6,7)], 2:4),
        ([(1,2),(10,2),(3,2),(1,2),(2,2),(1,0)], [1,2,3,4,5,6], [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7)], 2:6),
        ([(1,0),(1,2),(10,2),(3,2),(1,2),(2,2),(1,0)], [1,2,3,4,5,6], [(1,2),(2,3),(3,4),(4,5),(5,6),(6,7)], 1:6),
        ([(1,0),(1,1),(10,1),(3,1),(1,1),(2,1),(1,0)], [1,1,1,1,1,1], [(1,7)], 1:1)
    ]
        diagram = Diagram(path, system)
        display_diagram(system, diagram, expected_momenta, expected_unique, expected_indeterminate)
    end
end


text:
"""
    AbstractPolyOptions

Abstract base type for all generic options.
"""
abstract type AbstractPolyOptions end

"""
    AbstractPolyChoice

Abstract base type for dispatching extracted JSONs.
"""
abstract type AbstractPolyChoice end

"""
    PolyOpts <: AbstractPolyOptions

User options for the PolygonIO.

# Options
 * api_key: String representing the API key from a registered polygon.io account.
 * table: Whether to covnert extracted results as a table or not. Compatible Table.jl or nothing for raw JSONs.

# Examples
```julia-repl
julia> using PolygonIO, DataFrames
julia> opts = PolyOpts(API_KEY, DataFrame)  # if user wants tabularised output
julia> opts = PolyOpts(API_KEY, nothing)    # if user wants raw JSON output
```
"""
struct PolyOpts{T1<:AbstractString, T2<:Any} <: AbstractPolyOptions
    api_key::T1
    sink::T2
end

########## No Sink Structs  ##################
"""
Internal function to dispatch on no sink no results scenarios.
"""
struct NoSinkNoResults <: AbstractPolyChoice end

"""
Internal function to dispatch on no sink yes results scenarios.
"""
struct NoSinkYesResults <: AbstractPolyChoice end

"""
Internal function to dispatch on no sink yes tickers scenarios.
"""
struct NoSinkYesTickers <: AbstractPolyChoice end

"""
Internal function to dispatch on no sink yes ticker scenarios.
"""
struct NoSinkYesTicker <: AbstractPolyChoice end

"""
Internal function to dispatch on no sink yes last scenarios.
"""
struct NoSinkYesLast <: AbstractPolyChoice end

"""
Internal function to dispatch on no sink yes data scenarios.
"""
struct NoSinkYesData <: AbstractPolyChoice end

"""
Internal function to dispatch on no sink yes ticks scenarios.
"""
struct NoSinkYesTicks <: AbstractPolyChoice end

########### Yes Sink Structs ##################
"""
Internal function to dispatch on yes sink no last scenarios.
"""
struct YesSinkYesLast <: AbstractPolyChoice end

"""
Internal function to dispatch on yes sink yes tickers scenarios.
"""
struct YesSinkYesTickers <: AbstractPolyChoice end

"""
Internal function to dispatch on yes sink yes ticker scenarios.
"""
struct YesSinkYesTicker <: AbstractPolyChoice end

"""
Internal function to dispatch on yes sink no results scenarios.
"""
struct YesSinkNoResults <: AbstractPolyChoice end

"""
Internal function to dispatch on yes sink yes results scenarios.
"""
struct YesSinkYesResults <: AbstractPolyChoice end

"""
Internal function to dispatch on yes sink yes data scenarios.
"""
struct YesSinkYesData <: AbstractPolyChoice end

"""
Internal function to dispatch on yes sink yes ticks scenarios.
"""
struct YesSinkYesTicks <: AbstractPolyChoice end