text:
# using Revise
using ClusterManagers
using Distributed

worker_ids = if "SLURM_JOBID" in keys(ENV)
    wids = ClusterManagers.addprocs_slurm(parse(Int, ENV["SLURM_NTASKS"]))
    wids
else
    addprocs(2)
end

Distributed.@everywhere begin
    using Pkg
    Pkg.activate(".")
    Pkg.instantiate()
end

Distributed.@everywhere begin
    using POMDPs
    using BasicPOMCP
    using POMCPOW
    using VOOTreeSearch
    using ParticleFilters
    using CPUTime
    using DataFrames
    using ARDESPOT
    using Distributions
    using VDPTag2
    using ProgressMeter
    using Dates
    using CSV
    using POMDPSimulators
    using Random
    using Statistics
    using StatsBase
    using LinearAlgebra
    using VDPTag2Experiments
    using LinearAlgebra: diagm, PosDefException, det, eigvals, eigvecs
    using Random: MersenneTwister, seed!
    using Distances: euclidean

    using BOMCP
    using Statistics
    using StatsBase
    include("./deps/BOMCP.jl/scripts/lander_po.jl")

    function BOMCP.vectorize!(v, dims, x::MvNormal)
        v = copy(mean(x))
        return v
    end
end

@info "Started $(Distributed.nworkers()) workers..."
@info "Precompiling simulation code..."

max_time = 1.0
max_depth = 250

start_mean = [20.0, 7.5, 4.0, 2.0, 6.0, 15.0, 0.5]
start_cov = diagm([10.0^2, 5.0^2, 5.0^2, 5.0^2, 5.0^2, 10.0^2, 0.5^2])
d = MvNormal(start_mean, start_cov)
on_cluster = true
save_data = true
parallelize = true

if on_cluster
  max_query = 10_000
  K = 150 # number of parameter samples
  n = 20 # number of evaluation simulations
  m = 30 # number of elite samples
  max_iters = 50
else
  max_query = 100
  K = 15 # number of parameter samples
  n = 1 # number of evaluation simulations
  m = 10 # number of elite samples
  max_iters = 3
end

# Lander
pomdp = LunarLander()

println(stderr, "Max time: ", max_time)
println(stderr, "Max depth: ", max_depth)
println(stderr, "Max query: ", max_query)
println(stderr, "Start mean: ", start_mean)
println(stderr, "Start cov: ", start_cov)
println(stderr, "POMDP: ", pomdp)

Distributed.@everywhere begin
    struct SimWrap
        x::Vector{Float64}
        k::Int64
        i::Int64
        max_time::Float64
        max_depth::Int64
        max_query::Int64
        pomdp::POMDP
    end

    function gen_sim(sim_wrap::SimWrap)
        x = sim_wrap.x
        k = sim_wrap.k
        i = sim_wrap.i
        max_time = sim_wrap.max_time
        max_depth = sim_wrap.max_depth
        max_query = sim_wrap.max_query
        pomdp = sim_wrap.pomdp

        c = x[1]
        @assert c >= 0.0
        k_act = x[2]
        @assert k_act >= 1.0
        inv_alpha_act = x[3]
        @assert inv_alpha_act >= 0.1
        k_obs = x[4]
        @assert k_obs >= 1.0
        inv_alpha_obs = x[5]
        @assert inv_alpha_obs >= 0.1
        log_kernel = x[6]
        @assert log_kernel >= 1.0
        lambda = x[7]
        @assert lambda >= 0.1

        rng = MersenneTwister(7)
        belief_updater = EKFUpdater(pomdp, pomdp.Q.^2, pomdp.R.^2)
        rollout_policy = LanderPolicy(pomdp)

        action_selector = BOActionSelector(3, # action dims
                                    6, #belief dims
                                    false, #discrete actions
                                    kernel_params=[log(log_kernel), 0.0],
                                    k_neighbors = 5,
                                    belief_λ = lambda,
                                    lower_bounds = [-5.0, 0.0, -1.0],
                                    upper_bounds = [5.0, 15.0, 1.0],
                                    buffer_size=100,
                                    initial_action=rollout_policy
                                    )

        solver = BOMCPSolver(action_selector, belief_updater,
                    depth=max_depth, n_iterations=max_query,
                    max_time = max_time,
                    exploration_constant=c,
                    k_belief = k_obs,
                    alpha_belief = 1/inv_alpha_obs,
                    k_action = k_act,
                    alpha_action = 1/inv_alpha_act,
                    estimate_value=BOMCP.RolloutEstimator(rollout_policy),
                    default_action=BOMCP.ReportWhenUsed([0.0, 0.0, 0.0]),
                    rng=rng
                    )

        planner = solve(solver, pomdp)

        seed!(planner, i+40000*k)
        sim = Sim(pomdp,
                  planner,
                  belief_updater,
                  rng=MersenneTwister(i+50_000*k),
                  max_steps=200,
                  metadata=Dict(:i=>i, :k=>k)
                 )

        result = POMDPs.simulate(sim)
        output = (reward=discounted_reward(result),)
        return merge(sim.metadata, output)
    end
end

datestring = Dates.format(now(), "e_d_u_Y_HH_MM")
results_csv = DataFrame(c=Float64[], k_act=Float64[], inv_alpha_act=Float64[], k_obs=Float64[], inv_alpha_obs=Float64[], log_kernel=Float64[], lambda=Float64[])
try
    for i in 1:max_iters
        sims = SimWrap[]
        params = Vector{Float64}[]
        println(stderr,"creating $K simulation sets")
        for k in 1:K
            p = rand(d)
            p[1] = max(0.0, p[1])
            p[2] = max(1.0, p[2])
            p[3] = max(0.1, p[3])
            p[4] = max(1.0, p[4])
            p[5] = max(0.1, p[5])
            p[6] = max(1.0, p[6])
            p[7] = max(0.1, p[7])
            push!(params, p)
            for j in 1:n
                sim_wrap = SimWrap(p, k, j, max_time, max_depth, max_query, deepcopy(pomdp))
                push!(sims, sim_wrap)
            end
        end
        @assert length(params) == K
        if parallelize
            results = @showprogress pmap(gen_sim, sims)
        else
            results = @showprogress map(gen_sim, sims)
        end
        results = POMDPSimulators.create_dataframe(results)
        combined = by(results, :k) do df
            DataFrame(mean_reward=mean(df[:reward]))
        end
        order = sortperm(combined[:mean_reward])
        elite = params[combined[:k][order[K-m:end]]]
        elite_matrix = Matrix{Float64}(undef, length(start_mean), m)
        for k in 1:m
            elite_matrix[:,k] = elite[k]
        end
        try
            global d = fit(typeof(d), elite_matrix)
        catch ex
            if ex isa PosDefException
                println(stderr,"pos def exception")
                global d = fit(typeof(d), elite_matrix += 0.01*randn(size(elite_matrix)))
            else
                rethrow(ex)
            end
        end
        println(stderr, "Iteration $i")
        println(stderr, "Reward (mean): ", mean(combined[:mean_reward]))
        println(stderr, "Mean: ", mean(d))
        println(stderr, "Cov: ", det(cov(d)))
        ev = eigvals(cov(d))
        println(stderr, "Cov (eig): ", ev)
        for j in 1:length(ev)
            println(stderr, "Eigvecs: ", eigvecs(cov(d))[:,j])
        end
        push!(results_csv, mean(d))
        if save_data
            fname_csv = joinpath(@__DIR__, "../data", "Lander_bomcp_opt"*"_"*datestring*".csv")
            CSV.write(fname_csv, results_csv)
            println(stderr,"Saving parameters to ", fname_csv)
        end
    end
finally
    Distributed.rmprocs(worker_ids)
end


text:
using Documenter
import Pkg

# Get root directory
root_dir = dirname(@__DIR__)

# Install dependencies and import modules...
Pkg.activate(root_dir)
Pkg.instantiate()
using SEAL

# Define module-wide setups such that the respective modules are available in doctests
DocMeta.setdocmeta!(SEAL,
                    :DocTestSetup,
                    :(push!(LOAD_PATH, ".."); using SEAL);
                    recursive=true)

# Make documentation
makedocs(
    # Specify modules for which docstrings should be shown
    modules = [SEAL],
    # Set sitename to SEAL
    sitename = "SEAL.jl",
    # Set authors
    authors = "Michael Schlottke-Lakemper",
    # Provide additional formatting options
    format = Documenter.HTML(
        # Disable pretty URLs during manual testing
        prettyurls = get(ENV, "CI", nothing) == "true",
        # Explicitly add favicon as asset
        # assets = ["assets/favicon.ico"],
        # Set canonical URL to GitLab pages URL
        # canonical = canonical
    ),
    # Explicitly specify documentation structure
    pages = [
        "Home" => "index.md",
        "Reference" => "reference.md",
        "Contributing" => "contributing.md",
        "License" => "license.md"
    ],
)

deploydocs(
    repo = "github.com/JuliaCrypto/SEAL.jl.git",
    devbranch = "main"
)

ext: jl


text:
"""
    runsim(config, seed, prerun)

Performs a simulation run using configuration file `config`, random seed `seed`
and other settings provided via commandline, configuration file or the defaults.
"""
function runsim(config::String = "", seed::Integer = 0)
    initlogsettings(defaultSettings()) #needed for log calls during `getsettings()`
    settings = getsettings(config, seed)
    Random.seed!(settings["seed"])
    initlogsettings(settings)
    setupdatadir(settings)
    world = Patch[]
    timesteps = 0
    timeoffset = 0
    correctmode!(settings)
    for i in 1:length(settings["maps"])
        timeoffset += timesteps
        timesteps, maptable = readmapfile(settings["maps"][i], settings)
        i == 1 && (world = createworld(maptable, settings))
        i > 1 && updateworld!(world, maptable, settings)
        i == 1 && writedata(world, settings, timeoffset)
        simulate!(world, settings, timesteps, timeoffset)
    end
    world
end

"""
    rungemm(config, seed)

Wrapper for `runsim()`
Runs a simulation using configuration file `config`, random seed `seed`
and other settings provided via commandline, configuration file or the defaults.
"""
function rungemm(config::String = "", seed::Integer = 0)
    @time world = runsim(config, seed)
end

"""
    correctmode!(settings)

A heuristic function to keep backward compatibility. (The `mode` setting was
only introduced after several different experiments had already been carried
out with GeMM, so old config files don't include it.)
"""
function correctmode!(settings::Dict{String, Any})
    if settings["mode"] == "default"
        # detect invasion mode
        if settings["global-species-pool"] > 0 || settings["propagule-pressure"] > 0
            simlog("Detected implicit invasion mode, updated settings.", 'w')
            settings["mode"] = "invasion"
        end
        # XXX further modes?
    end
end

max_stars_repo_name: lleiding/gemm

max_issues_repo_name: lleiding/gemm

max_forks_repo_name: lleiding/gemm


text:
using testPkg
using Test

@testset "testPkg.jl" begin
    @test 1 == 1
end

max_stars_repo_name: tseaward/testPkg

max_issues_repo_name: tseaward/testPkg

max_forks_repo_name: tseaward/testPkg


text:
# This file is a "Hello, world!" in Julia language for wandbox.

println("Hello, Wandbox!")

# Julia language references:
#   https://docs.julialang.org/en/v1/


text:
using Serialization, DataFrames, JLD2, Statistics, Plots
using GEAS:QTL, top_QTL

function sum_2021_07(rpt)
    @load "dat/run/base.jld2" base
    qtl = deserialize("dat/2021-06-30/qtl-$rpt.ser")
    tq = top_QTL(base[:hap], qtl[2], n=10) # top QTL loci
    prd = deserialize("dat/2021-06-30/prd-c-$rpt.ser")
    snp = deserialize("dat/2021-06-30/snp-c-$rpt.ser")
    st = snp[tq, :]             # only the top QTL
    frq = zeros(11, 10)
    i = start = 1
    for stop in 4000:4000:44000
        frq[i, :] = mean(st[:, start:stop], dims=2)'
        start = stop + 1
        i += 1
    end
    plot(frq, leg=false)
end

function improvement()
    m2 = zeros(11, 3)
    tt = ['a', 'b', 'c']

    for i in 1:3
        for j in 1:10
            w = lpad(string(j), 2, '0')
            prd = deserialize("dat/2021-06-30/prd-$(tt[i])-$w.ser")
            bt = combine(groupby(prd, :g8n), :t2c => mean => :m)
            m2[:, i] += bt.m
        end
    end
    m2 ./= 10
    plot(m2, dpi=300)
end


text:
using Distributed
@everywhere using Test
@everywhere using FileTrees
@everywhere using FileTrees: Thunk, NoValue
@everywhere using Dates
@everywhere import FilePathsBase: /

@testset "maketree" begin
    global t = maketree(["a" => ["b" => ["a"],
                                 "c" => ["b", "a", "c"=>[]]]])

    @test strip(repr(t)) == """
                            ./
                            └─ a/
                               ├─ b/
                               │  └─ a
                               └─ c/
                                  ├─ b
                                  ├─ a
                                  └─ c/"""
end

# Helper function to recursively check that all children of ft have ft as their parent
function isconsistent(ft)
    ok = true
    for c in children(ft)
        ok = parent(c) === ft && isconsistent(c)
    end
    return ok
end

@testset "indexing" begin
    @test name(t) == "."
    @test t["a"]["b"]["a"] isa File
    @test t["a"]["c"]["c"] isa FileTree

    @test Path(t["a"]["b"]["a"])  == p"." /"a"/"b"/"a"

    t1 = FileTrees.rename(t, "foo")
    @test Path(t1["a"]["b"]["a"])  == p"foo" /"a"/"b"/"a"

    @test isequal(t[r"a|b"], t)

    @test isempty(t[r"^c"])

    @test isequal(t["a"]["c"][["b", r"a|d"]],
                  maketree("c" => ["b","a"]))

    @test_throws ErrorException t["c"]

    @test isconsistent(t[r"a|c"])
    @test isconsistent(t[glob"a/*/*"])
end

@testset "filter" begin
    @test isequal(filter(f->f isa FileTree || name(f) == "c", t),
                  maketree(["a"=>["b"=>[],
                                  "c"=>["c"=>[]]]]))
end

@testset "merge" begin
    @test isequal(merge(t,t), t)

    # Check that we don't rearrange the tree more than necessary
    t2 =  maketree(["b" => ["x"], "a" => ["y"]])
    t3 =  maketree(["c" => ["x"], "a" => ["z"]])
    @test name.(children(merge(t2, t3))) == ["b", "a", "c"]
    @test name.(children(merge(t3, t2))) == ["c", "a", "b"]
    @test isequal(merge(t2, t3), merge(t3,t2))
end

@testset "diff" begin
    @test isempty(FileTrees.diff(t,t))
    @test isequal(FileTrees.diff(t,maketree([])), t)
end

import FileTrees: attach

@testset "touch mv and cp" begin
    global t1 = maketree([])
    for yr = 1:9
        global t1 = mkpath(t1, joinpath(string("0", yr)), value=yr)
        for mo = 1:6
            global t1 = attach(t1, joinpath(string("0", yr), string("0", mo)),
                        File(nothing, "data.csv", (yr, mo)) )
        end
    end

    @test isconsistent(t1)

    data1 = reducevalues(vcat, t1) |> exec

    t4 = mv(t1, r"^(.*)/(.*)/data.csv$", s"\1/\2.csv")

    @test isconsistent(t4)

    t5 = maketree([])
    for yr = 1:9
        for mo = 1:6
            t5 = attach(t5, joinpath(string("0", yr)),
                        File(nothing, string("0", mo) * ".csv", (yr, mo)))
        end
    end

    @test isconsistent(t5)

    @test isequal(t4, t5)

    data2 = reducevalues(vcat, t4) |> exec

    @test data1 == data2

    t7 = mapsubtrees(t1, glob"*") do subtree
        reducevalues(vcat, subtree)
    end

    @test isconsistent(t7)

    @test reducevalues(hcat, t7; dirs=true) == [(j, i) for i=1:6, j=1:9]

    t6 = cp(t1, r"^(.*)/(.*)/data.csv$", s"\1/\2.csv")

    @test isconsistent(t6)

    @test isequal(t6, merge(t1, t5))

    @testset "combine" begin
        tcombine = maketree(["a" => [(name="y", value="ay"), (name ="x", value="ax")], "b" => [(name ="x", value="bx"), (name="y", value="by")]])
        tcombined = mv(tcombine, r"^([^/]*)/([x|y])", s"\2"; combine=(v1,v2) -> v1 * "_" * v2)

        @test tcombined["x"][] == "ax_bx"
        @test tcombined["y"][] == "ay_by"
        # Also test that we maintained the same order as the first encountered node
        @test name.(children(tcombined)) == ["y", "x"]
    end
end

@testset "values" begin
    t1 = FileTrees.load(path, t)
    if isdir("test_dir")
        rm("test_dir", recursive=true)
    end

    @test get(t1["a/b/a"]) == string(p"." / "a" / "b" / "a")

    @test reducevalues(*, mapvalues(lowercase, t1)) == lowercase(reducevalues(*, t1))

    FileTrees.save(maketree("test_dir" => [t1])) do f
        @test f isa File
        open(path(f), "w") do io
            print(io, get(f))
        end
    end

    t2 = FileTree("test_dir")
    t3 = FileTrees.load(t2) do f
        open(path(f), "r") do io
            String(read(io))
        end
    end

    t4 = filter(!isempty, t1)

    @test isequal(t3, FileTrees.rename(t4, "test_dir"))
    if isdir("test_dir")
        rm("test_dir", recursive=true)
    end

    x1 = maketree("a"=>[(name="b", value=1)])
    x2 = mapvalues(x->NoValue(), x1, lazy=true)
    @test !isempty(values(x2))
    @test isempty(values(exec(x2)))
    x3 = mapvalues(x->rand(), x2)
    @test !isempty(values(x3))
    @test isempty(values(exec(x3)))

    # issue 16
    @test_throws ArgumentError reducevalues(+, maketree("." => []))
    @test reducevalues(+, maketree("." => []), init=0) === 0

    @test_throws ArgumentError reducevalues(+, maketree("." => []), associative=false)
    @test reducevalues(+, maketree("." => []), init=0, associative=false) === 0

    # issue 23
    @test FileTrees.save(identity, maketree([])) == nothing
end

@testset "lazy-exec" begin

    if isdir("test_dir_lazy")
        rm("test_dir_lazy", recursive=true)
    end

    t1 = FileTrees.load(uppercase∘path, t, lazy=true)

    @test get(t1["a/b/a"]) isa Thunk
    @test get(exec(t1)["a/b/a"]) == string(p"."/"A"/"B"/"A")

    @test exec(reducevalues(*, mapvalues(lowercase, t1))) == lowercase(exec(reducevalues(*, t1)))

    s = FileTrees.save(maketree("test_dir_lazy" => [t1])) do f
        open(path(f), "w") do io
            print(io, get(f))
        end
    end

    @test isdir("test_dir_lazy")
    @test isfile("test_dir_lazy/a/b/a")

    t2 = FileTree("test_dir_lazy")
    t3 = FileTrees.load(t2; lazy=true) do f
        open(path(f), "r") do io
            (String(read(io)), now())
        end
    end
    toc = now()
    sleep(0.01)
    tic = exec(reducevalues((x,y)->x, mapvalues(last, t3)))

    @test tic > toc

    t4 = filter(!isempty, t1) |> exec

    t5 = mapvalues(first, t3) |> exec
    @test isequal(t5, FileTrees.rename(t4, "test_dir_lazy"))

    if isdir("test_dir_lazy")
        rm("test_dir_lazy", recursive=true)
    end
end

@testset "exec with context" begin
    import Dagger

    struct SpecialContext end

    computespecial = Ref(false)
    function Dagger.compute(::SpecialContext, t::Dagger.Thunk, kws...)
        computespecial[] = true
        compute(Dagger.Context(procs()), t; kws...)
    end
    ncollectspecial = Ref(0)
    function Dagger.collect(::SpecialContext, c::Dagger.Chunk; kws...)
        ncollectspecial[] += 1
        collect(Dagger.Context(procs()), c; kws...)
    end

    t = mapvalues(identity, maketree("a" => ["b" => [(name="c", value=1)], "d" => ["e" => [(name="f", value=2), (name="g", value=3)]]]), lazy=true)

    @test exec(SpecialContext(), t) |> values == [1,2,3]
    @test computespecial[] == true # Dummy compare to make failed test outprint a little less confusing
    @test ncollectspecial[] == 3 # All values are collected with SpecialContext
end

@testset "iterators" begin
    @test values(t1) == map(get, filter(hasvalue, nodes(t1)))
    @test values(t1, dirs=false) != values(t1, dirs=true)
    @test values(t1, dirs=false) == reducevalues(vcat, t1)
    @test values(t1, dirs=true) == map(get, filter(hasvalue, nodes(t1)))
    @test nodes(t1) == collect(FileTrees.PostOrderDFS(t1))
    @test nodes(t1, dirs=false) == collect(Iterators.filter(x->x isa File, FileTrees.PostOrderDFS(t1)))
    @test files(t1) == collect(nodes(t1, dirs=false))
    @test dirs(t1) == collect(Iterators.filter(x->x isa FileTree, nodes(t1, dirs=true)))
end


text:
# FunctionPolicy
# A policy represented by a function
# maintained by @zsunberg

"""
FunctionPolicy

Policy `p=FunctionPolicy(f)` returns `f(x)` when `action(p, x)` is called.
"""
struct FunctionPolicy{F<:Function} <: Policy
    f::F
end

"""
FunctionSolver

Solver for a FunctionPolicy.
"""
mutable struct FunctionSolver{F<:Function} <: Solver
    f::F
end

solve(s::FunctionSolver, mdp::Union{MDP,POMDP}) = FunctionPolicy(s.f)

action(p::FunctionPolicy, x) = p.f(x)

updater(p::FunctionPolicy) = PreviousObservationUpdater()

max_stars_repo_name: JuliaPOMDP/POMDPPolicies

max_issues_repo_name: JuliaPOMDP/POMDPPolicies

max_forks_repo_name: JuliaPOMDP/POMDPPolicies


text:
stripunits(x) = uconvert(NoUnits, x)

@testset "spring models" begin

    springmodel = SpringModel([1.0], (0,4π), [1.0, 0.0])
    solslow = solve(odeproblem(springmodel))
    @test abs(solslow(π/2)[1]) < 1e-6
    @test abs(solslow(3π/2)[1]) < 1e-4

    accelerate!(springmodel, 4.0)
    solfast = solve(odeproblem(springmodel), atol=1e-8)
    @test abs(solfast(π/2)[1]+1) < 1e-5
    @test_skip abs(solfast(3π/2)[1]+1) < 1e-4
    @test abs(solfast(π/4)[1]) < 1e-5
    @test abs(solfast(3π/4)[1]) < 1e-4

end

@testset "SIR" begin
    function test1()
        sir = SemanticModels.BasicSIR()

        β = sir.equations[1].parameters[1]
        @test (4*β.unit + 6β.unit) * 3u"s" == u"30person"
        @test u"24person" + u"1kperson" == u"1024person"
        @test u"24person" + u"1kperson/s"*u"3s" == u"3024person"
        @test_throws DimensionError 4*β.unit + (6β.unit) * 3u"s"
    end
    test1()
end

function RealSIR()
    # β = NumParameter(:β, u"person/s", TransitionRate)
    # γ = NumParameter(:γ, u"person/s", TransitionRate)

    # S = NumVariable(:S, u"person", Amount)
    # I = NumVariable(:I, u"person", Amount)
    # R = NumVariable(:R, u"person", Amount)

    # sir = SIR([ODE(:(dS/dt -> -βSI/N),
    #             [β], [S]),
    #             ODE(:(dI/dt -> βSI/N -γI),
    #                 [β, γ], [S, I]),
    #         ODE(:(dR/dt -> γI), [γ], [I])],
    #         NumVariable.([:N,:S,:I,:R], u"person", Amount),
    #         NumVariable(:t, u"s", Amount))

    return odeproblem(SIRSimulation([100,1,0.0],
                      (0,100.0),
                       SIRParams(0.50,5.0)))
end
# test case where nothing happens because of *0.0
    # return odeproblem(SIRSimulation([100,0,0.0],
    #                   (0,10.0),
    #                    SIRParams(1.0,0.2)))
@testset "real sir" begin
    sir = RealSIR()
    solsir = solve(sir)
    @test solsir[end][1] < 1
    @test solsir[end][end] > 99
end

@testset "spring_units" begin
    springmodel = SpringModel([u"1.0s^-2"], (u"0s",4π*u"s"), [u"1.0m", u"0m/s"])
    prob = odeproblem(springmodel)
    sol = solve(springmodel)
    t = u"π/2*1s"
    v = (1e-6)*u"m"
    @test sol(u"4π*1s")[1] - u"1m" < u"1e-4*1m"
end

@testset "sir_meters" begin
    initialpop = [99, 1, 0.0].*u"m"
    prob = SIRSimulation(initialpop, (u"0.0s", u"200s"), SIRParams(u"40.0/s", u"20.0m/s"))
    sol = solve(odeproblem(prob), alg=Vern9(),  dt = u"0.1s")
    @test sol(sol.t[end])[1] < u"1e-4*1m"
    @test sol(sol.t[end])[end] > u"99m"
end

@testset "sir_units" begin
    initialpop = [99, 1, 0.0].*u"person"
    prob = SIRSimulation(initialpop, (u"0.0minute", u"75minute"), SIRParams(u"40.0/minute", u"20.0person/minute"))
    sol = solve(odeproblem(prob), alg=Vern9(),  dt = u"0.1minute")
    @test sol(sol.t[end])[1] < u"1e-4*1person"
    @test sol(sol.t[end])[end] > u"99*1person"
end

@testset "combined_meters" begin
    springmodel = SpringModel([u"1.0s^-2"], (u"0s",4π*u"s"), [u"1.0m", u"0m/s"])
    prob = odeproblem(springmodel)
    sol = solve(springmodel)
    @test sol(u"4π*1s")[1] - u"1m" < u"1e-4*1m"

    initialS = sol.u[end][1] * 100
    initialI = abs(sol.u[end][2] *u"200s" + u"1m")
    initialpop = [initialS, initialI, u"0.0m"]
    sirprob = SIRSimulation(initialpop, (u"0.0s", u"20s"), SIRParams(u"40.0/s", u"20.0m/s"))
    sirsol = solve(odeproblem(sirprob), alg=Vern9(),  dt = u"0.1s")

    # test that everything is converted from S to R
    @test sirsol(sirsol.t[end])[1] < u"1e-4*1m"
    @test sirsol(sirsol.t[end])[end] > 0.99 * sirsol(sirsol.t[1])[1]

end

@testset "CombinedModel" begin
    springmodel = SpringModel([u"1.0s^-2"], (u"0s",4π*u"s"), [u"1.0m", u"0m/s"])
    function create_sir(m, solns)
        sol = solns[1]
        initialS = sol.u[end][1] * 100
        initialI = abs(sol.u[end][2] *u"200s" + u"1m")
        initialpop = [initialS, initialI, u"0.0m"]
        sirprob = SIRSimulation(initialpop, (u"0.0s", u"20s"), SIRParams(u"40.0/s", u"20.0m/s"))
        return sirprob
    end
    cm = CombinedModel([springmodel], create_sir)

    sirsol = solve(cm)

    # test that everything is converted from S to R
    @test sirsol(sirsol.t[end])[1] < u"1e-4*1m"
    @test sirsol(sirsol.t[end])[end] > 0.99 * sirsol(sirsol.t[1])[1]
end

@testset "Regression Example" begin
    START_YR = 2000
    STOP_YR = 2005

    df = generate_synthetic_data(100, START_YR, STOP_YR)
    yhat = ols(df, rand(Uniform(50,100)), STOP_YR+1)
end

@testset "FluModel" begin

    tfinal = 4π*u"d"
    springmodel = SpringModel([u"1.0d^-2"], (u"0d",tfinal), [u"25.0C", u"0C/d"])
    function create_sir(m, solns)
        sol = solns[1]
        initialS = u"10000person" #sol.u[end][1] * 100
        initialI = abs(sol.u[end][1] *u"2person/C" + u"1person")
        initialpop = [initialS, initialI, u"0.0person"]
        γ = u"3.0person/d" / u"C" * sol(sol.t[end])[1]
        @show γ
        sirprob = SIRSimulation(initialpop, (u"0.0d", u"20d"), SIRParams(u"40.0/d", γ))
        return sirprob
    end

    function create_flu(cm, solns)
        sol = solns[1]
        finalI = stripunits(sol(sol.t[end])[2])
        population = stripunits(sol(sol.t[end])[2])
        # population = stripunits(sum(sol.u[end]))
        df = generate_synthetic_data(population, 0,100)
        f = @formula(vaccines_produced ~ flu_patients)
        model =  lm(f, df[2:length(df.year), [:year, :flu_patients, :vaccines_produced]])
        println("GLM Model:")
        println(model)

        year_to_predict = 50
        num_flu_patients_from_sim = finalI
        vaccines_produced = missing
        targetDF = DataFrame(year=year_to_predict, flu_patients=num_flu_patients_from_sim, vaccines_produced=missing)
        @show targetDF

        # predicted_num_vaccines = predict(model, targetDF)
        # println("Predicted number of vaccines based on simulated number of flu patients for year ", year_to_predict, " = ", ceil(predicted_num_vaccines[1]))
        # β = solve(RegressionProblem(f, model, targetDF, missing))
        return RegressionProblem(f, model, targetDF, missing)
    end
    cm = CombinedModel([springmodel], create_sir)
    flumodel = CombinedModel([cm], create_flu)
    sol = solve(flumodel)
    @test ceil(Int, sol[1]) >= 4393

    spsol = solve(springmodel)
    # plot(map(x->spsol(x)[1]/Unitful.C, collect(spsol.t[1]:0.1s:spsol.t[end])))
    print(spsol)
    print(solve(cm).u)
    # plot(solve(cm))
    print(solve(flumodel))
end

@testset "knowledge" begin
isa = "isa"
unit = "Unit"
knowledge = [
    :(u"m", isa, unit),
    (u"s", isa, unit),
    (u"m", "instantiates", "Distance"),
    (u"s", isa, "Time"),
    (u"person", isa, unit),
    (u"person", "instantiates", "Quantity"),
    ("sir","solvedby","ODEs"),
    :(springmodel,solvedby,ODEs),
    :(
        (x, represented, m),
        (springmodel, isa, ODEModel),
        (springmodel, var(1), x),
        (springmodel, var(2), dx/dt),
        (sirmodel, var(1), S),
        (sirmodel, var(2), I),
        (sirmodel, var(2), R)
    )
]
path(x, R) = :(
    (springmodel, var(1), x),
    (x, measured_in, m),
    (S, measured_in, m),
    (sirmodel, var(1), S),
    (sirmodel, var(3), R)
)
end


text:
struct Tabular{T}
    data::Observable{T}
end

render_row_value(x::AbstractFloat) = round(x, sigdigits=6)
render_row_value(x::Integer) = x
render_row_value(x::Missing) = "n/a"
render_row_value(x::String) = x
render_row_value(x) = string(x)

function jsrender(session::Session, tabular::Tabular)
    metadata = map(session, tabular.data) do data
        cols = Tables.columns(data)
        names = collect(Tables.columnnames(cols))
        columns = [map(render_row_value, Tables.getcolumn(cols, name)) for name in names]
        (; names, columns)
    end
    JSServe.register_resource!(session, metadata)

    table_div = DOM.div()
    create_table = js"""
        function create_table(div, m) {
            const names = m.names;
            const columns = m.columns;
            const indices = Object.keys(columns[0])

            const columnDefs = names.map(name => {
                return {field: name, headerName: name};
            });
            const rowData = indices.map(i => {
                const row = {}
                for (let colIdx in columns) {
                    row[names[colIdx]] = columns[colIdx][i];
                }
                return row
            });

            const rem2pix = parseFloat(getComputedStyle(document.documentElement).fontSize);
            const rowHeight = 2.5 * rem2pix;
            const headerHeight = 3.75 * rem2pix;

            // let the grid know which columns and what data to use
            const gridOptions = {
                defaultColDef: {
                    sortable: true,
                    resizable: true,
                    headerClass: "text-blue-800 text-xl font-semibold px-4 hover:bg-gray-200",
                    cellClass: "px-4 text-base",
                },
                columnDefs: columnDefs,
                rowData: rowData,
                headerHeight: headerHeight,
                getRowHeight: params => {
                    return rowHeight;
                },
                getRowClass: params => {
                    return "hover:bg-gray-200";
                },
                suppressFieldDotNotation: true,
            };

            new $(agGrid).Grid(div, gridOptions);
        }
    """
    onload(session, table_div, js"div => ($create_table)(div, JSServe.get_observable($metadata))")
    onjs(session, metadata, js"""
        function (m) {
            const table_div = $(table_div);
            while (table_div.firstChild) {
                table_div.removeChild(table_div.lastChild);
            }
            ($create_table)(table_div, m);
        }
    """)
    return table_div
end


text:
module Refactoring

    using MacroTools, MLStyle

    export extract_method, unassigned_variables, @search_and_replace

    include("typed_expressions.jl")
    include("extract_method.jl")
    include("search_and_replace.jl")

end # module


text:
"""
    module Faiss

An interface to the Faiss library for similarity-search of vectors (i.e. nearest-neighbour searching).

For basic usage, see [`Index`](@ref), [`add!`](@ref) and [`search`](@ref).
"""
module Faiss

using PythonCall

const faiss = PythonCall.pynew()

function __init__()
    PythonCall.pycopy!(faiss, pyimport("faiss"))
end

"""
    Index(dim, spec="Flat")

Create a Faiss index of the given dimension.

The `spec` is an index factory string describing the type of index to construct.
"""
struct Index
    py::Py
end

Base.size(x::Index) = (size(x, 1), size(x, 2))

Base.size(x::Index, i::Integer) = i == 1 ? pyconvert(Int, x.py.d) : i == 2 ? pyconvert(Int, x.py.ntotal) : error()

function Base.show(io::IO, ::MIME"text/plain", x::Index)
    print(io, typeof(x), " of ", size(x, 2), " vectors of dimension ", size(x, 1))
end

Index(d::Integer, str::AbstractString="Flat") = Index(faiss.index_factory(convert(Int, d), convert(String, str)))

"""
    add!(idx::Index, vs::AbstractMatrix)

Add the columns of `vs` to the index.
"""
function add!(x::Index, vs::AbstractMatrix)
    size(vs, 1) == size(x, 1) || error("expecting $(size(x, 1)) rows")
    vs_ = Py(convert(AbstractMatrix{Float32}, vs)').__array__()
    x.py.add(vs_)
    return x
end

"""
    search(idx::Index, vs::AbstractMatrix, k::Integer)

Search the index for the `k` nearest neighbours of each column of `vs`.

Return `(D, I)` where `I` is a matrix where each column gives the ids of the `k` nearest
neighbours of the corresponding column of `vs` and `D` is the corresponding matrix of
distances.
"""
function search(x::Index, vs::AbstractMatrix, k::Integer)
    size(vs, 1) == size(x, 1) || error("expecting $(size(x, 1)) rows")
    vs_ = Py(convert(AbstractMatrix{Float32}, vs)').__array__()
    k_ = convert(Int, k)
    D_, I_ = x.py.search(vs_, k_)
    D = PyArray{Float32,2,true,true,Float32}(D_.T)
    I = PyArray{Int64,2,true,true,Int64}(I_.T)
    return (D, I)
end

"""
    downcast(idx::Index)

Return the same index downcasted to its most specific type.
"""
downcast(x::Index) = Index(faiss.downcast_index(x.py))

end # module


text:
@testset "prearith $T" for T in (Double16, Double32, Double64)

    val = T(17)/3
    negval = -val

    @test -val == negval
    @test abs(val) == val
    @test abs(negval) == val
    @test flipsign(val, -1) == negval
    @test flipsign(negval, 1) == negval
    @test copysign(val, -1) == negval
    @test copysign(negval, 1) == val

    @test ldexp(frexp(val)...,) == val
   # @test ldexp(significand(val), exponent(val)) == val

    @test signs(val) == (sign(HI(val)), sign(LO(val)))

end

@testset "trunc $T" for T in (Double16, Double32, Double64)
   for I in (:Int16, :Int32, :Int64, :Int128)
       @test 3 == @eval trunc($I, $T(pi))
   end
end

@testset "floor $T" for T in (Double16, Double32, Double64)
   for I in (:Int16, :Int32, :Int64, :Int128)
       @test 3 == @eval floor($I, $T(pi))
   end
end

@testset "ceil $T" for T in (Double16, Double32, Double64)
   for I in (:Int16, :Int32, :Int64, :Int128)
       @test 4 == @eval ceil($I, $T(pi))
   end
end

@testset "round $T" for T in (Double16, Double32, Double64)
   for I in (:Int16, :Int32, :Int64, :Int128)
       @test 3 == @eval round($I, $T(pi), RoundNearest)
       @test 3 == @eval round($I, $T(pi), RoundNearestTiesAway)
       @test 4 == @eval round($I, $T(pi), RoundUp)
       @test 3 == @eval round($I, $T(pi), RoundDown)
       @test 3 == @eval round($I, $T(pi), RoundToZero)
   end
end

@testset "spread $T" for T in (Double16, Double32, Double64)
    @test 4 == @eval spread($T(pi))
end

@testset "sld $T" for T in (Double16, Double32, Double64)
    @test 7.0 == @eval sld($T(pi), 0.5)
end

@testset "tld $T" for T in (Double16, Double32, Double64)
    @test 6.0 == @eval tld($T(pi), 0.5)
end


text:
struct Hᶜ{T}
  ∑h::InfiniteSum{T}
  Hᴸ::InfiniteMPS
  Hᴿ::InfiniteMPS
  ψ::InfiniteCanonicalMPS
  n::Int
end

struct Hᴬᶜ{T}
  ∑h::InfiniteSum{T}
  Hᴸ::InfiniteMPS
  Hᴿ::InfiniteMPS
  ψ::InfiniteCanonicalMPS
  n::Int
end

function (H::Hᶜ)(v)
  Hᶜv = H * v
  ## return Hᶜv * δˡ * δʳ
  return noprime(Hᶜv)
end

function (H::Hᴬᶜ)(v)
  Hᴬᶜv = H * v
  ## return Hᶜv * δˡ⁻¹ * δˢ * δʳ
  return noprime(Hᴬᶜv)
end

# Struct for use in linear system solver
struct Aᴸ
  ψ::InfiniteCanonicalMPS
  n::Int
end

function (A::Aᴸ)(x)
  ψ = A.ψ
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  ψ̃ = prime(linkinds, ψᴴ)
  n = A.n

  N = length(ψ)
  #@assert n == N

  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)

  xT = translatecell(translator(ψ), x, -1)
  for k in (n - N + 1):n
    xT = xT * ψ.AL[k] * ψ̃.AL[k]
  end
  δˡ = δ(l[n], l′[n])
  δʳ = δ(r[n], r′[n])
  xR = x * ψ.C[n] * ψ′.C[n] * dag(δʳ) * denseblocks(δˡ)
  return xT - xR
end

function left_environment(hᴸ, 𝕙ᴸ, ψ; tol=1e-15)
  ψ̃ = prime(linkinds, dag(ψ))
  N = nsites(ψ)

  Aᴺ = Aᴸ(ψ, N)
  Hᴸᴺ¹, info = linsolve(Aᴺ, 𝕙ᴸ[N], 1, -1; tol=tol)
  # Get the rest of the environments in the unit cell
  Hᴸ = InfiniteMPS(Vector{ITensor}(undef, N))
  Hᴸ[N] = Hᴸᴺ¹
  for n in 1:(N - 1)
    Hᴸ[n] = Hᴸ[n - 1] * ψ.AL[n] * ψ̃.AL[n] + hᴸ[n]
  end
  return Hᴸ
end

# Struct for use in linear system solver
struct Aᴿ
  hᴿ::InfiniteMPS
  ψ::InfiniteCanonicalMPS
  n::Int
end

function (A::Aᴿ)(x)
  hᴿ = A.hᴿ
  ψ = A.ψ
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  ψ̃ = prime(linkinds, ψᴴ)
  n = A.n

  N = length(ψ)
  @assert n == N

  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)

  xT = x
  for k in reverse(1:N)
    xT = xT * ψ.AR[k] * ψ̃.AR[k]
  end
  xT = translatecell(translator(ψ), xT, 1)
  δˡ = δ(l[n], l′[n])
  δʳ = δ(r[n], r′[n])
  xR = x * ψ.C[n] * ψ′.C[n] * δˡ * denseblocks(dag(δʳ))
  return xT - xR
end

function right_environment(hᴿ, ψ; tol=1e-15)
  ψ̃ = prime(linkinds, dag(ψ))
  # XXX: replace with `nsites`
  #N = nsites(ψ)
  N = length(ψ)

  A = Aᴿ(hᴿ, ψ, N)
  Hᴿᴺ¹, info = linsolve(A, hᴿ[N], 1, -1; tol=tol)
  # Get the rest of the environments in the unit cell
  Hᴿ = InfiniteMPS(Vector{ITensor}(undef, N))
  Hᴿ[N] = Hᴿᴺ¹
  for n in reverse(1:(N - 1))
    Hᴿ[n] = Hᴿ[n + 1] * ψ.AR[n + 1] * ψ̃.AR[n + 1] + hᴿ[n]
  end
  return Hᴿ
end

# TODO Generate all environments, why? Only one is needed in the sequential version
function right_environment(hᴿ, 𝕙ᴿ, ψ; tol=1e-15)
  ψ̃ = prime(linkinds, dag(ψ))
  N = nsites(ψ)

  A = Aᴿ(hᴿ, ψ, N)
  Hᴿᴺ¹, info = linsolve(A, 𝕙ᴿ[N], 1, -1; tol=tol)
  # Get the rest of the environments in the unit cell
  Hᴿ = InfiniteMPS(Vector{ITensor}(undef, N))
  Hᴿ[N] = Hᴿᴺ¹
  for n in reverse(1:(N - 1))
    Hᴿ[n] = Hᴿ[n + 1] * ψ.AR[n + 1] * ψ̃.AR[n + 1] + hᴿ[n]
  end
  return Hᴿ
end

function tdvp_iteration(args...; multisite_update_alg="sequential", kwargs...)
  if multisite_update_alg == "sequential"
    return tdvp_iteration_sequential(args...; kwargs...)
  elseif multisite_update_alg == "parallel"
    return tdvp_iteration_parallel(args...; kwargs...)
  else
    error(
      "Multisite update algorithm multisite_update_alg = $multisite_update_alg not supported, use \"parallel\" or \"sequential\"",
    )
  end
end

function tdvp(
  solver::Function,
  ∑h,
  ψ;
  maxiter=10,
  tol=1e-8,
  outputlevel=1,
  multisite_update_alg="sequential",
  time_step,
  solver_tol=(x -> x / 100),
)
  N = nsites(ψ)
  (ϵᴸ!) = fill(tol, nsites(ψ))
  (ϵᴿ!) = fill(tol, nsites(ψ))
  if outputlevel > 0
    println("Running VUMPS with multisite_update_alg = $multisite_update_alg")
    flush(stdout)
    flush(stderr)
  end
  for iter in 1:maxiter
    ψ, (eᴸ, eᴿ) = tdvp_iteration(
      solver,
      ∑h,
      ψ;
      (ϵᴸ!)=(ϵᴸ!),
      (ϵᴿ!)=(ϵᴿ!),
      multisite_update_alg=multisite_update_alg,
      solver_tol=solver_tol,
      time_step=time_step,
    )
    ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
    maxdimψ = maxlinkdim(ψ[0:(N + 1)])
    if outputlevel > 0
      println(
        "VUMPS iteration $iter (out of maximum $maxiter). Bond dimension = $maxdimψ, energy = $((eᴸ, eᴿ)), ϵᵖʳᵉˢ = $ϵᵖʳᵉˢ, tol = $tol",
      )
      flush(stdout)
      flush(stderr)
    end
    if ϵᵖʳᵉˢ < tol
      println(
        "Precision error $ϵᵖʳᵉˢ reached tolerance $tol, stopping VUMPS after $iter iterations (of a maximum $maxiter).",
      )
      flush(stdout)
      flush(stderr)
      break
    end
  end
  return ψ
end

function vumps_solver(M, time_step, v₀, solver_tol)
  λ⃗, v⃗, info = eigsolve(M, v₀, 1, :SR; ishermitian=true, tol=solver_tol)
  return λ⃗[1], v⃗[1], info
end

return function tdvp_solver(M, time_step, v₀, solver_tol)
  v, info = exponentiate(M, time_step, v₀; ishermitian=true, tol=solver_tol)
  v = v / norm(v)
  return nothing, v, info
end

function vumps(
  args...; time_step=-Inf, eigsolve_tol=(x -> x / 100), solver_tol=eigsolve_tol, kwargs...
)
  @assert isinf(time_step) && time_step < 0
  println("Using VUMPS solver with time step $time_step")
  flush(stdout)
  flush(stderr)
  return tdvp(vumps_solver, args...; time_step=time_step, solver_tol=solver_tol, kwargs...)
end

function tdvp(args...; time_step, solver_tol=(x -> x / 100), kwargs...)
  solver = if !isinf(time_step)
    println("Using TDVP solver with time step $time_step")
    flush(stdout)
    flush(stderr)
    tdvp_solver
  elseif time_step < 0
    # Call VUMPS instead
    println("Using VUMPS solver with time step $time_step")
    flush(stdout)
    flush(stderr)
    vumps_solver
  else
    error("Time step $time_step not supported.")
  end
  return tdvp(solver, args...; time_step=time_step, solver_tol=solver_tol, kwargs...)
end


text:
#   Unit "ajd.jl" of the Diagonalization.jl Package for Julia language
#
#   MIT License
#   Copyright (c) 2019, 2020
#   Marco Congedo, CNRS, Grenoble, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements Approximate Joint Diagonalization.

"""
```
(1)
function ajd(𝐂::ℍVector;
             trace1    :: Bool   = false,
             w         :: Union{Tw, Function} = ○,
          algorithm :: Symbol = :NoJoB,
          preWhite  :: Bool   = false,
          sort      :: Bool   = true,
          init      :: Mato   = ○,
          tol       :: Real   = 1e-6,
          maxiter   :: Int    = _maxiter(algorithm, eltype(𝐂[1])),
          verbose   :: Bool   = false,
          threaded  :: Bool   = true,
        eVar     :: TeVaro    = _minDim(𝐂),
        eVarC    :: TeVaro    = ○,
        eVarMeth :: Function  = searchsortedfirst,
        simple   :: Bool      = false)

(2)
function ajd(𝐗::VecMat;
             covEst     :: StatsBase.CovarianceEstimator = SCM,
             dims       :: Into = ○,
             meanX      :: Into = 0,
          trace1     :: Bool = false,
          w          :: Twf  = ○,
       algorithm :: Symbol = :NoJoB,
       preWhite  :: Bool = false,
       sort      :: Bool = true,
       init      :: Mato = ○,
       tol       :: Real = 1e-6,
       maxiter   :: Int  = _maxiter(algorithm, eltype(𝐗[1])),
       verbose   :: Bool = false,
       threaded  :: Bool = true,
     eVar     :: TeVaro    = _minDim(𝐗),
     eVarC    :: TeVaro    = ○,
     eVarMeth :: Function  = searchsortedfirst,
     simple   :: Bool      = false)

```

Return a [LinearFilter](@ref) object:

**(1) Approximate joint diagonalization** of the set of ``k``
symmetric or Hermitian matrices `𝐂`, of type
[ℍVector](https://marco-congedo.github.io/PosDefManifold.jl/dev/MainModule/#%E2%84%8DVector-type-1)
using the given solving `algorithm` (*NoJoB* by default).

If `trace1` is true, all matrices in the set
`𝐂` are normalized so as to have trace equal to 1.
It is false by default. This option applies only for solving algorithms
that are not invariante by scaling, that is, those based on the least-squares
(Frebenius) criterion. See [Algorithms](@ref).

if `w` is a `StatsBase.AbstractWeights`, the weights are applied to the set
`𝐂`. If `w` is a `Function`, the weights are found passing each matrix
in the set to such function. An appropriate choice for AJD algorithms
minimizing a least-squares criterion, like *OJoB* and *NoJoB*, is the
[`nonDiagonality`](@ref) function (Congedo et al.(2008)[🎓](@ref)).
By default, no weights are applied.

If `preWhite` is true the solution is found by the two-step procedure
described here above in section [pre-whitening for AJD](@ref).
By default, it is false.
Dimensionality reduction can be obtained at this stage using arguments
`eVarC` and `eVarMeth`, in the same way they are used to find
the [subspace dimension](@ref) ``p``, but using the accumulated
regularized eigenvalues of

``\\frac{1}{k}\\sum_{l=1}^kC_k``.

The default values are:
- `eVarC` is set to 0.999
- `eVarMeth=searchsortedfirst`.

If `sort` is true (default), the vectors in `.F` are permuted as explained
here above in [permutation for AJD](@ref), otherwise they will be in arbitrary order.

Regarding arguments `init`, `tol` and `maxiter`, see [Algorithms](@ref).

If `verbose` is true (false by default), the convergence attained
at each iteration will be printed in the REPL.

`eVar` and `eVarMeth` are used to define a
[subspace dimension](@ref) ``p`` using the accumulated regularized
eigenvalues in Eq. [ajd.7].

The default values are:
- `eVar` is set to the dimension of the matrices in `𝐂`
- `eVarMeth=searchsortedfirst`.

Note that passing `nothing` or a real nummber as `eVar` (see
[subspace dimension](@ref)) is meningful only if `sort` is set to true
(default) and if the input matrices ``{C_1,...,C_k}`` are
positive definite.

If `simple` is set to `true`, ``p`` is set equal to the dimension
of the matrices ``{C_1,...,C_k}`` and only the fields `.F` and `.iF`
are written in the constructed object.
This corresponds to the typical output of AJD algorithms.

if `threaded`=true (default) and the number of threads Julia is instructed
to use (the output of Threads.nthreads()), is higher than 1,
AJD algorithms supporting multi-threading run in multi-threaded mode.
See [Algorithms](@ref) and
[these notes](https://marco-congedo.github.io/PosDefManifold.jl/dev/MainModule/#Threads-1)
on multi-threading.

**(2) Approximate joint diagonalization**
with a set of ``k`` data matrices `𝐗` as input; the
covariance matrices of the set are estimated using arguments
`covEst`, `dims` and `meanX`
(see [covariance matrix estimations](@ref))
and passed to method (1) with the
remaining arguments of method (2).

**See also:** [PCA](@ref), [CSP](@ref), [mAJD](@ref).

**Examples:**

```
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

const err=1e-6

# method (1) real
t, n, k=50, 10, 10
A=randn(n, n) # mixing matrix in model x=As
Xset = [genDataMatrix(t, n) for i = 1:k]
Xfixed=randn(t, n)./1
for i=1:length(Xset) Xset[i]+=Xfixed end
Cset = ℍVector([ℍ((Xset[s]'*Xset[s])/t) for s=1:k])
aC=ajd(Cset; algorithm=:OJoB, simple=true)
aC2=ajd(Cset; algorithm=:NoJoB, simple=true)
aC3=ajd(Cset; algorithm=:LogLike, simple=true)
aC4=ajd(Cset; algorithm=:LogLikeR, simple=true)
aC5=ajd(Cset; algorithm=:JADE, simple=true)
aC6=ajd(Cset; algorithm=:JADEmax, simple=true)
aC7=ajd(Cset; algorithm=:GAJD, simple=true)
aC8=ajd(Cset; algorithm=:QNLogLike, simple=true)

# a=ajd(Cset; algorithm=:GAJD2, simple=true, verbose=true)

# method (2) real
aX=ajd(Xset; algorithm=:OJoB, simple=true)
aX2=ajd(Xset; algorithm=:NoJoB, simple=true)
aX3=ajd(Xset; algorithm=:LogLike, simple=true)
aX4=ajd(Xset; algorithm=:LogLikeR, simple=true)
aX5=ajd(Xset; algorithm=:JADE, simple=true)
aX6=ajd(Xset; algorithm=:JADEmax, simple=true)
aX7=ajd(Xset; algorithm=:GAJD, simple=true)
aX8=ajd(Xset; algorithm=:QNLogLike, simple=true)

@test aX≈aC
@test aX2≈aC2
@test aX3≈aC3
@test aX4≈aC4
@test aX5≈aC5
@test aX6≈aC6
@test aX7≈aC7
@test aX8≈aC8

# method (1) complex
t, n, k=50, 10, 10
Ac=randn(ComplexF64, n, n) # mixing matrix in model x=As
Xcset = [genDataMatrix(ComplexF64, t, n) for i = 1:k]
Xcfixed=randn(ComplexF64, t, n)./1
for i=1:length(Xcset) Xcset[i]+=Xcfixed end
Ccset = ℍVector([ℍ((Xcset[s]'*Xcset[s])/t) for s=1:k])
aCc=ajd(Ccset; algorithm=:OJoB, simple=true)
aCc2=ajd(Ccset; algorithm=:NoJoB, simple=true)
aCc3=ajd(Ccset; algorithm=:LogLike, simple=true)
aCc4=ajd(Ccset; algorithm=:JADE, simple=true)
aCc5=ajd(Ccset; algorithm=:JADEmax, simple=true)

# method (2) complex
aXc=ajd(Xcset; algorithm=:OJoB, simple=true)
aXc2=ajd(Xcset; algorithm=:NoJoB, simple=true)
aXc3=ajd(Xcset; algorithm=:LogLike, simple=true)
aXc4=ajd(Xcset; algorithm=:JADE, simple=true)
aXc5=ajd(Xcset; algorithm=:JADEmax, simple=true)

@test aXc≈aCc
@test aXc2≈aCc2
@test aXc3≈aCc3
@test aXc4≈aCc4
@test aXc5≈aCc5

# create 20 REAL random commuting matrices
# they all have the same eigenvectors
Cset2=PosDefManifold.randP(3, 20; eigvalsSNR=Inf, commuting=true)
# estimate the approximate joint diagonalizer (AJD)
a=ajd(Cset2; algorithm=:OJoB)
# the orthogonal AJD must be equivalent to the eigenvector matrix
# of any of the matrices in Cset
@test norm([spForm(a.F'*eigvecs(C)) for C ∈ Cset2])/20 < err
# do the same for JADE algorithm
a=ajd(Cset2; algorithm=:JADE)
@test norm([spForm(a.F'*eigvecs(C)) for C ∈ Cset2])/20 < err
# do the same for JADEmax algorithm
a=ajd(Cset2; algorithm=:JADEmax)
@test norm([spForm(a.F'*eigvecs(C)) for C ∈ Cset2])/20 < err

# generate positive definite matrices with model A*D_κ*D, where
# A is the mixing matrix and D_κ, for all κ=1:k, are diagonal matrices.
# The estimated AJD matrix must be the inverse of A
# and all transformed matrices bust be diagonal
n, k=3, 10
Dest=PosDefManifold.randΛ(eigvalsSNR=10, n, k)
# make the problem identifiable
for i=1:k Dest[k][1, 1]*=i/(k/2) end
for i=1:k Dest[k][3, 3]/=i/(k/2) end
A=randn(n, n) # non-singular mixing matrix
Cset3=Vector{Hermitian}([Hermitian(A*D*A') for D ∈ Dest])
a=ajd(Cset3; algorithm=:NoJoB, eVarC=n)
@test spForm(a.F'*A)<√err
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err
a=ajd(Cset3; algorithm=:LogLike, eVarC=n)
@test spForm(a.F'*A)<√err
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err
a=ajd(Cset3; algorithm=:LogLikeR, eVarC=n)
@test spForm(a.F'*A)<√err
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err
a=ajd(Cset3; algorithm=:GAJD, eVarC=n)
@test spForm(a.F'*A)<√err
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err
a=ajd(Cset3; algorithm=:QNLogLike, eVarC=n)
@test spForm(a.F'*A)<√err
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<err

# Do the same thing for orthogonal diagonalizers:
# now A will be orthogonal
O=randU(n) # orthogonal mixing matrix
Cset4=Vector{Hermitian}([Hermitian(O*D*O') for D ∈ Dest])
a=ajd(Cset4; algorithm=:OJoB, eVarC=n)
@test spForm(a.F'*O)<√err
@test mean(nonD(a.F'*Cset4[i]*a.F) for i=1:k)<err
a=ajd(Cset4; algorithm=:JADE, eVarC=n)
@test spForm(a.F'*O)<√err
@test mean(nonD(a.F'*Cset4[i]*a.F) for i=1:k)<√err
a=ajd(Cset4; algorithm=:JADEmax, eVarC=n)
@test spForm(a.F'*O)<√err
@test mean(nonD(a.F'*Cset4[i]*a.F) for i=1:k)<√err

# repeat the test adding noise; now the model is no more exactly identifiable
for k=1:length(Cset3) Cset3[k]+=randP(n)/1000 end
a=ajd(Cset3; algorithm=:NoJoB, eVarC=n)
@test spForm(a.F'*A)<err^(1/6)
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<√err
a=ajd(Cset3; algorithm=:LogLike, eVarC=n)
@test spForm(a.F'*A)<err^(1/6)
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<√err
a=ajd(Cset3; algorithm=:LogLikeR, eVarC=n)
@test spForm(a.F'*A)<err^(1/6)
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<√err
a=ajd(Cset3; algorithm=:GAJD, eVarC=n)
@test spForm(a.F'*A)<err^(1/6)
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<√err
a=ajd(Cset3; algorithm=:QNLogLike, eVarC=n)
@test spForm(a.F'*A)<err^(1/6)
@test mean(nonD(a.F'*Cset3[i]*a.F) for i=1:k)<√err

# the same thing for orthogonal diagonalizers
for k=1:length(Cset4) Cset4[k]+=randP(n)/1000 end
a=ajd(Cset4; algorithm=:OJoB, eVarC=n)
@test spForm(a.F'*O)<err^(1/6)
@test mean(nonD(a.F'*Cset4[i]*a.F) for i=1:k)<√err
a=ajd(Cset4; algorithm=:JADE, eVarC=n)
@test spForm(a.F'*O)<err^(1/6)
@test mean(nonD(a.F'*Cset4[i]*a.F) for i=1:k)<√err
a=ajd(Cset4; algorithm=:JADEmax, eVarC=n)
@test spForm(a.F'*O)<err^(1/6)
@test mean(nonD(a.F'*Cset4[i]*a.F) for i=1:k)<√err

# create 20 COMPLEX random commuting matrices
# they all have the same eigenvectors
Ccset2=PosDefManifold.randP(ComplexF64, 3, 20; eigvalsSNR=Inf, commuting=true)
# estimate the approximate joint diagonalizer (AJD)
ac=ajd(Ccset2; algorithm=:OJoB)
# he AJD must be equivalent to the eigenvector matrix of any of the matrices in Cset
# just a sanity check as rounding errors appears for complex data
@test norm([spForm(ac.F'*eigvecs(C)) for C ∈ Ccset2])/20<√err
# do the same for JADE algorithm
ac=ajd(Ccset2; algorithm=:JADE)
@test norm([spForm(ac.F'*eigvecs(C)) for C ∈ Ccset2])/20<√err
# do the same for JADEmax algorithm
ac=ajd(Ccset2; algorithm=:JADEmax)
@test norm([spForm(ac.F'*eigvecs(C)) for C ∈ Ccset2])/3<√err

# the same thing using the NoJoB and LogLike algorithms. Require less precision
# as the NoJoB solution is not constrained in the orthogonal group
ac=ajd(Ccset2; algorithm=:NoJoB)
@test norm([spForm(ac.F'*eigvecs(C)) for C ∈ Ccset2])/20<√err
ac=ajd(Ccset2; algorithm=:LogLike)
@test norm([spForm(ac.F'*eigvecs(C)) for C ∈ Ccset2])/20<√err

# REAL data:
# normalize the trace of input matrices,
# give them weights according to the `nonDiagonality` function
# apply pre-whitening and limit the explained variance both
# at the pre-whitening level and at the level of final vector selection
Cset=PosDefManifold.randP(20, 80; eigvalsSNR=10, SNR=10, commuting=false)

a=ajd(Cset; algorithm=:OJoB, trace1=true, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:NoJoB, trace1=true, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:LogLike, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:LogLikeR, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:JADE, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:JADEmax, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:GAJD, w=nonD, preWhite=true, eVarC=4, eVar=0.99)
a=ajd(Cset; algorithm=:QNLogLike, w=nonD, preWhite=true, eVarC=4, eVar=0.99)

# AJD for plots below
a=ajd(Cset; algorithm=:QNLogLike, verbose=true, preWhite=true)

using Plots
# plot the original covariance matrices
# and their transformed counterpart
CMax=maximum(maximum(abs.(C)) for C ∈ Cset);
 h1 = heatmap(Cset[1], clim=(-CMax, CMax), title="C1", yflip=true, c=:bluesreds);
 h2 = heatmap(Cset[2], clim=(-CMax, CMax), title="C2", yflip=true, c=:bluesreds);
 h3 = heatmap(Cset[3], clim=(-CMax, CMax), title="C3", yflip=true, c=:bluesreds);
 h4 = heatmap(Cset[4], clim=(-CMax, CMax), title="C4", yflip=true, c=:bluesreds);
 📈=plot(h1, h2, h3, h4, size=(700,400))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigAJD1.png")

Dset=[a.F'*C*a.F for C ∈ Cset];
 DMax=maximum(maximum(abs.(D)) for D ∈ Dset);
 h5 = heatmap(Dset[1], clim=(-DMax, DMax), title="F'*C1*F", yflip=true, c=:bluesreds);
 h6 = heatmap(Dset[2], clim=(-DMax, DMax), title="F'*C2*F", yflip=true, c=:bluesreds);
 h7 = heatmap(Dset[3], clim=(-DMax, DMax), title="F'*C3*F", yflip=true, c=:bluesreds);
 h8 = heatmap(Dset[4], clim=(-DMax, DMax), title="F'*C4*F", yflip=true, c=:bluesreds);
 📉=plot(h5, h6, h7, h8, size=(700,400))
# savefig(📉, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigAJD2.png")

```

 ![Figure AJD1](assets/FigAJD1.png)

 ![Figure AJD2](assets/FigAJD2.png)

```

# COMPLEX data:
# normalize the trace of input matrices,
# give them weights according to the `nonDiagonality` function
# apply pre-whitening and limit the explained variance both
# at the pre-whitening level and at the level of final vector selection
Ccset=PosDefManifold.randP(3, 20; eigvalsSNR=10, SNR=2, commuting=false)

ac=ajd(Ccset; trace1=true, w=nonD, preWhite=true,
       algorithm=:OJoB, eVarC=8, eVar=0.99)
ac=ajd(Ccset; eVarC=8, eVar=0.99)
ac=ajd(Ccset; algorithm=:LogLike, eVarC=8, eVar=0.99)
ac=ajd(Ccset; algorithm=:JADE, eVarC=8, eVar=0.99)
ac=ajd(Ccset; algorithm=:JADEmax, eVarC=8, eVar=0.99)

```

"""
function ajd(𝐂::ℍVector;
             trace1    :: Bool   = false,
             w         :: Union{Tw, Function} = ○,
          algorithm :: Symbol = :QNLogLike,
          preWhite  :: Bool   = false,
          sort      :: Bool   = true,
          init      :: Mato   = ○,
          tol       :: Real   = 1e-06,
          maxiter   :: Int    = _maxiter(algorithm, eltype(𝐂[1])),
          verbose   :: Bool   = false,
          threaded  :: Bool   = true,
        eVar     :: TeVaro    = _minDim(𝐂),
        eVarC    :: TeVaro    = ○,
        eVarMeth :: Function  = searchsortedfirst,
        simple   :: Bool      = false)

   eVar isa Int && eVarC isa Int && eVar>eVarC && throw(ArgumentError(📌*", ajd function: `eVar` cannot be larger than `eVarC`"))
   args=("Approximate Joint Diagonalization", false)
   k, n=length(𝐂), size(𝐂[1], 1)
   simple && (whitening=false)

   if    algorithm == :QNLogLike
         U, V, λ, iter, conv=qnLogLike(𝐂;
            w=w, preWhite=preWhite, sort=sort,
            init=init, tol=tol, maxiter=maxiter, verbose=verbose,
            eVar=eVarC, eVarMeth=eVarMeth,
            threaded=threaded)
   elseif algorithm==:JADEmax
          U, V, λ, iter, conv=jademax(𝐂;
            trace1=trace1,
            w=w, preWhite=preWhite, sort=sort, init=init, tol=tol,
            maxiter=maxiter, verbose=verbose, eVar=eVarC, eVarMeth=eVarMeth)
   elseif algorithm==:JADE
          U, V, λ, iter, conv=jade(𝐂;
            trace1=trace1,
            w=w, preWhite=preWhite, sort=sort, init=init, tol=tol,
            maxiter=maxiter, verbose=verbose, eVar=eVarC, eVarMeth=eVarMeth)
   elseif algorithm ∈(:OJoB, :NoJoB)
          U, V, λ, iter, conv=JoB(reshape(𝕄Vector(𝐂), (k, 1, 1)), 1, k, :c, algorithm, eltype(𝐂[1]);
            trace1=trace1, threaded=threaded,
            w=w, preWhite=preWhite, sort=sort, init=init, tol=tol,
            maxiter=maxiter, verbose=verbose, eVar=eVarC, eVarMeth=eVarMeth)
   elseif algorithm == :LogLike
          U, V, λ, iter, conv=logLike(𝐂;
            w=w, preWhite=preWhite, sort=sort, init=init, tol=tol,
            maxiter=maxiter, verbose=verbose, eVar=eVarC, eVarMeth=eVarMeth)
   elseif algorithm == :LogLikeR
          U, V, λ, iter, conv=logLikeR(𝐂;
            w=w, preWhite=preWhite, sort=sort, init=init, tol=tol,
            maxiter=maxiter, verbose=verbose, eVar=eVarC, eVarMeth=eVarMeth)
   elseif algorithm==:GAJD
          U, V, λ, iter, conv=gajd(𝐂;
            trace1=trace1,
            w=w, preWhite=preWhite, sort=sort, init=init, tol=tol,
            maxiter=maxiter, verbose=verbose, eVar=eVarC, eVarMeth=eVarMeth)
   elseif algorithm==:GLogLike
          U, V, λ, iter, conv=gLogLike(𝐂;
            w=w, preWhite=preWhite, sort=sort, init=init, tol=tol,
            maxiter=maxiter, verbose=verbose, eVar=eVarC, eVarMeth=eVarMeth)
   else
      throw(ArgumentError(📌*", ajd constructor: invalid `algorithm` argument: $algorithm"))
   end

   # alert user if eVar passed as an integer exceeds the minimum whitening dim
   whiteDim=size(U, 2)
   if eVar isa Int && eVar>whiteDim
     @warn(📌*", ajd function: the whitening step reduced the rank to $(whiteDim); `eVar` has been lowered to this value.")
     eVar=whiteDim
   end

   λ = _checkλ(λ) # make sure no imaginary noise is present (for complex data)

   simple ? LF(U, V, Diagonal(λ), ○, ○, ○, args...) :
   begin
      # println(λ)
      p, arev = _getssd!(eVar, λ, n, eVarMeth) # find subspace
      LF(U[:, 1:p], V[1:p, :], Diagonal(λ[1:p]), arev[p], λ, arev, args...)
   end
end

function ajd(𝐗::VecMat;
             covEst     :: StatsBase.CovarianceEstimator = SCM,
             dims       :: Into = ○,
             meanX      :: Into = 0,
          trace1     :: Bool = false,
          w          :: Twf  = ○,
       algorithm :: Symbol = :QNLogLike,
       preWhite  :: Bool   = false,
       sort      :: Bool   = true,
       init      :: Mato   = ○,
       tol       :: Real   = 1e-06,
       maxiter   :: Int    = _maxiter(algorithm, eltype(𝐗[1])),
       verbose   :: Bool   = false,
       threaded  :: Bool   = true,
     eVar     :: TeVaro   = _minDim(𝐗),
     eVarC    :: TeVaro   = ○,
     eVarMeth :: Function = searchsortedfirst,
     simple   :: Bool     = false)

   dims===○ && (dims=_set_dims(𝐗))
   _check_data(𝐗, dims, covEst, meanX, w)===○ && return
   # check algo
   ajd(_cov(𝐗; covEst=covEst, dims=dims, meanX=meanX);
       trace1=trace1, w=w, algorithm=algorithm, preWhite=preWhite, sort=sort,
       init=init, tol=tol, maxiter=maxiter, verbose=verbose, threaded=threaded,
       eVar=eVar, eVarC=eVarC, eVarMeth=eVarMeth, simple=simple)
end

"""
```
function majd(𝑿::VecVecMat;
              covEst     :: StatsBase.CovarianceEstimator = SCM,
              dims       :: Into    = ○,
              meanX      :: Into    = 0,
          algorithm :: Symbol    = :NoJoB,
          fullModel :: Bool      = false,
          preWhite  :: Bool      = false,
          sort      :: Bool      = true,
          init      :: VecMato   = ○,
          tol       :: Real      = 0.,
          maxiter   :: Int       = _maxiter(algorithm, eltype(𝑿[1][1])),
          verbose   :: Bool      = false,
          threaded  :: Bool      = true,
        eVar     :: TeVaro   = _minDim(𝑿),
        eVarC    :: TeVaro   = ○,
        eVarMeth :: Function = searchsortedfirst,
        simple   :: Bool     = false)

```

Return a [LinearFilter](@ref) object.

**Multiple Approximate Joint Diagonalization** of the ``k`` sets
of ``m`` data matrices `𝐗`
using the given solving `algorithm` (*NoJoB* by default).

If `fullModel` is true, the [gmca.3] problem here above is solved,
otherwise (default), the [gmca.2] problem here above is solved.

If `preWhite` the two-step procedure explained here above in the
section [pre-whitening for MAJD](@ref) is used.
Dimensionality reduction can be obtained at this stage using arguments
`eVarC` and `eVarMeth`.

The default values are:
- `eVarC` is set to 0.999
- `eVarMeth=searchsortedfirst`.

If `sort` is true (default), the column vectors of the matrices ``F_1,...,F_m``
are signed and permuted
as explained here above in [permutation for MAJD](@ref),
otherwise they will have arbitrary sign and will be in arbitrary order.

If `verbose` is true (false by default), the convergence attained
at each iteration will be printed in the REPL.

`eVar` and `eVarMeth` are used to define a
[subspace dimension](@ref) ``p`` using the accumulated regularized
eigenvalues in Eq. [gmca.7]

The default values are:
- `eVar` is set to the minimum dimension of the matrices in `𝐗`
- `eVarMeth=searchsortedfirst`.

If `simple` is set to `true`, ``p`` is set equal to the dimension
of the covariance matrices that are computed on the matrices in `𝐗`,
which depends on the choice of `dims`,
and only the fields `.F` and `.iF`
are written in the constructed object.
This corresponds to the typical output of approximate diagonalization
algorithms.

if `threaded`=true (default) and the number of threads Julia is instructed
to use (the output of Threads.nthreads()), is higher than 1,
solving algorithms supporting multi-threading run in multi-threaded mode.
See [Algorithms](@ref) and
[these notes](https://marco-congedo.github.io/PosDefManifold.jl/dev/MainModule/#Threads-1)
on multi-threading.

**See also:** [gMCA](@ref), [gCCA](@ref), [AJD](@ref).

**Examples:**

```
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

##  Create data for testing the case k>1, m>1 ##
# `t` is the number of samples,
# `m` is the number of datasets,
# `k` is the number of observations,
# `n` is the number of variables,
# `noise` must be smaller than 1.0. The smaller the noise, the more data are correlated
# Output k vectors of m data data matrices
function getData(t, m, k, n, noise)
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(n) for i=1:m] # random orthogonal matrices
    # variables common to all subjects with unique variance profile across k
    X=[(abs2.(randn(n))).*randn(n, t) for s=1:k]
    # each subject has this common part plus a random part
    𝐗=[[𝐕[i]*((1-noise)*X[s] + noise*randn(n, t)) for i=1:m] for s=1:k]
    return 𝐗, 𝐕
end

function getData(::Type{Complex{T}}, t, m, k, n, noise) where {T<:AbstractFloat}
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(ComplexF64, n) for i=1:m] # random orthogonal matrices
    # variables common to all subjects with unique variance profile across k
    X=[(abs2.(randn(n))).*randn(ComplexF64, n, t) for s=1:k]
    # each subject has this common part plus a random part
    𝐗=[[𝐕[i]*((1-noise)*X[s] + noise*randn(ComplexF64, n, t)) for i=1:m] for s=1:k]
    return 𝐗, 𝐕
end

# REAL data
# do joint blind source separation of non-stationary data
t, m, n, k, noise = 200, 5, 4, 6, 0.1
Xset, Vset=getData(t, m, k, n, noise)
𝒞=Array{Matrix}(undef, k, m, m)
for s=1:k, i=1:m, j=1:m 𝒞[s, i, j]=(Xset[s][i]*Xset[s][j]')/t end

aX=majd(Xset; fullModel=true, algorithm=:OJoB)
# the spForm index of the estimated demixing matrices times the true
# mixing matrix must be low
@test mean(spForm(aX.F[i]'*Vset[i]) for i=1:m)<0.1

# test the same using NoJoB algorithm
aX=majd(Xset; fullModel=true, algorithm=:NoJoB)
@test mean(spForm(aX.F[i]'*Vset[i]) for i=1:m)<0.1

# plot the original cross-covariance matrices and the rotated
# cross-covariance matrices

# Get all products 𝐔[i]' * 𝒞[l, i, j] * 𝐔[j]
function _rotate_crossCov(𝐔, 𝒞, m, k)
    𝒮=Array{Matrix}(undef, k, m, m)
    @inbounds for l=1:k, i=1:m, j=1:m 𝒮[l, i, j]=𝐔[i]'*𝒞[l, i, j]*𝐔[j] end
    return 𝒮
end

# Put all `k` cross-covariances in a single matrix
# of dimension m*n x m*n for visualization
function 𝒞2Mat(𝒞::AbstractArray, m, k)
    n=size(𝒞[1, 1, 1], 1)
    C=Matrix{Float64}(undef, m*n, m*n)
    for i=1:m, j=1:m, x=1:n, y=1:n C[i*n-n+x, j*n-n+y]=𝒞[k, i, j][x, y] end
    return C
end

using Plots

Cset=[𝒞2Mat(𝒞, m, s) for s=1:k]
 Cmax=maximum(maximum(abs.(C)) for C ∈ Cset)
 h1 = heatmap(Cset[1], clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-cov, k=1")
 h2 = heatmap(Cset[2], clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-cov, k=2")
 h3 = heatmap(Cset[2], clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-cov, k=3")
 h4 = heatmap(Cset[2], clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-cov, k=4")
 h5 = heatmap(Cset[2], clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-cov, k=5")
 h6 = heatmap(Cset[2], clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-cov, k=6")
 📈=plot(h1, h2, h3, h4, h5, h6, size=(1200,550))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigmAJD1.png")

𝒮=_rotate_crossCov(aX.F, 𝒞, m, k)
 Sset=[𝒞2Mat(𝒮, m, s) for s=1:k]
 Smax=maximum(maximum(abs.(S)) for S ∈ Sset)
 h11 = heatmap(Sset[1], clim=(-Smax, Smax), yflip=true, c=:bluesreds, title="rotated cross-cov, k=1")
 h12 = heatmap(Sset[2], clim=(-Smax, Smax), yflip=true, c=:bluesreds, title="rotated cross-cov, k=2")
 h13 = heatmap(Sset[2], clim=(-Smax, Smax), yflip=true, c=:bluesreds, title="rotated cross-cov, k=3")
 h14 = heatmap(Sset[2], clim=(-Smax, Smax), yflip=true, c=:bluesreds, title="rotated cross-cov, k=4")
 h15 = heatmap(Sset[2], clim=(-Smax, Smax), yflip=true, c=:bluesreds, title="rotated cross-cov, k=5")
 h16 = heatmap(Sset[2], clim=(-Smax, Smax), yflip=true, c=:bluesreds, title="rotated cross-cov, k=6")
 📉=plot(h11, h12, h13, h14, h15, h16, size=(1200,550))
# savefig(📉, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FigmAJD2.png")

```
![Figure mAJD1](assets/FigmAJD1.png)

![Figure mAJD2](assets/FigmAJD2.png)

In the bottom figures here above, the rotated cross-covariance matrices have the
expected *strip-diagonal* form, that is, each block
``F_i^T\\frac{1}{T}(X_{li}X_{lj}^T)F_j``,
for ``l∈[1,...,k]``, ``i,j∈[1,...,m]``, is approximately diagonal.

```
# COMPLEX data
# do joint blind source separation of non-stationary data
t, m, n, k, noise = 200, 5, 4, 6, 0.1
Xcset, Vcset=getData(ComplexF64, t, m, k, n, noise)
𝒞=Array{Matrix}(undef, k, m, m)
for s=1:k, i=1:m, j=1:m 𝒞[s, i, j]=(Xcset[s][i]*Xcset[s][j]')/t end

aXc=majd(Xcset; fullModel=true, algorithm=:OJoB)
# the spForm index of the estimated demixing matrices times the true
# mixing matrix must be low
@test mean(spForm(aXc.F[i]'*Vcset[i]) for i=1:m)<0.1

# test the same using NoJoB algorithm
aXc=majd(Xcset; fullModel=true, algorithm=:NoJoB)
@test mean(spForm(aXc.F[i]'*Vcset[i]) for i=1:m)<0.1
```

"""
function majd(𝑿::VecVecMat;
              covEst     :: StatsBase.CovarianceEstimator = SCM,
              dims       :: Into    = ○,
              meanX      :: Into    = 0,
          algorithm :: Symbol    = :NoJoB,
          fullModel :: Bool      = false,
          preWhite  :: Bool      = false,
          sort      :: Bool      = true,
          init      :: VecMato   = ○,
          tol       :: Real      = 0.,
          maxiter   :: Int       = _maxiter(algorithm, eltype(𝑿[1][1])),
          verbose   :: Bool      = false,
          threaded  :: Bool      = true,
        eVar     :: TeVaro   = _minDim(𝑿),
        eVarC    :: TeVaro   = ○,
        eVarMeth :: Function = searchsortedfirst,
        simple   :: Bool     = false)

   if dims===○ dims=_set_dims(𝑿) end
   #TODO _check_data(𝑿, dims, covEst, meanX, ○)===○ && return
   eVar isa Int && eVarC isa Int && eVar>eVarC && throw(ArgumentError(📌*", majd function: `eVar` cannot be larger than `eVarC`"))
   (n, t)=dims==1 ? reverse(size(𝑿[1][1])) : size(𝑿[1][1])
   k, m=length(𝑿), length(𝑿[1])
   args=("Multiple Approximate Joint Diagonalization", false)

   if algorithm ∈(:OJoB, :NoJoB)
      𝐔, 𝐕, λ, iter, conv=JoB(𝑿, m, k, :d, algorithm, eltype(𝑿[1][1]);
               covEst=covEst, dims=dims, meanX=meanX,
               fullModel=fullModel, preWhite=preWhite, sort=sort, init=init,
               tol=tol, maxiter=maxiter, verbose=verbose, threaded=threaded,
               eVar=eVarC, eVarMeth=eVarMeth)
   # elseif...
   else
      throw(ArgumentError(📌*", majd constructor: invalid `algorithm` argument"))
   end

   # alert user if eVar passed as an integer exceeds the minimum whitening dim
   whiteDim=size(𝐔[1], 2)
   if eVar isa Int && eVar>whiteDim
     @warn(📌*", majd function: the whitening step reduced the rank to $(whiteDim); `eVar` has been lowered to this value.")
     eVar=whiteDim
   end

   λ = _checkλ(λ) # make sure no imaginary noise is present (for complex data)

   simple ? LF(𝐔, 𝐕, Diagonal(λ), ○, ○, ○, args...) :
   begin
      p, arev = _getssd!(eVar, λ, n, eVarMeth) # find subspace
      LF([𝐔[i][:, 1:p] for i=1:m], [𝐕[i][1:p, :] for i=1:m], Diagonal(λ[1:p]), arev[p], λ, arev, args...)
   end
end


text:
using Turing

@testset "MCMCChains and Turing utilities" begin
    using Distributions, Random
    Random.seed!(112)
    # simulated samples for μ
    samples = randn(100, 1, 1)
    data = randn(50)
    chain = Chains(samples)

    compute_loglike(μ, data) = logpdf(Normal(μ, 1), data)
    compute_loglike(μ, σ, data) = logpdf(Normal(μ, σ), data)

    pll1 = pointwise_log_likelihoods(compute_loglike, chain, data)
    pll2 = pointwise_log_likelihoods(compute_loglike, samples, data)
    # the pointwise log likehoods should be the same for both methods
    @test pll1 ≈ pll2 atol = 1e-6
    # test the dimensions: data points, samples, chains
    @test size(pll1) == (50, 100, 1)
    # test that sum of pointwise log likelihoods equals sum of log likelihoods
    @test sum(sum(map(s -> logpdf.(Normal(s, 1), data), samples))) ≈ sum(pll1) atol =
        1e-6
    # test that psis_loo works with MCMCChains and yields correct type
    psis_loo_output = psis_loo(compute_loglike, chain, data)
    @test isa(psis_loo_output, PsisLoo)
    # test that loo works with MCMCChains and yields correct type
    psis_output = loo(compute_loglike, chain, data)
    @test isa(psis_output, PsisLoo)
    # test that psis works with MCMCChains and yields correct type
    psis_output = psis(compute_loglike, chain, data)
    @test isa(psis_output, Psis)

    # ensure that methods work with r_eff argument
    r_eff = similar(pll2, 0)
    # test that psis_loo works with MCMCChains and yields correct type
    psis_loo_output = psis_loo(compute_loglike, chain, data; r_eff=r_eff)
    @test isa(psis_loo_output, PsisLoo)
    # test that loo works with MCMCChains and yields correct type
    psis_output = loo(compute_loglike, chain, data; r_eff=r_eff)
    @test isa(psis_output, PsisLoo)
    # test that psis works with MCMCChains and yields correct type
    psis_output = psis(compute_loglike, chain, data; r_eff=r_eff)
    @test isa(psis_output, Psis)

    @model function model(data)
        μ ~ Normal(0, 1)
        σ ~ truncated(Cauchy(0, 1), 0, Inf)
        for i in eachindex(data)
            data[i] ~ Normal(μ, σ)
        end
    end

    chain = sample(model(data), NUTS(1000, 0.9), MCMCThreads(), 1000, 4)
    pw_lls_turing = pointwise_log_likelihoods(model(data), chain)
    pw_lls_loglike = pointwise_log_likelihoods(compute_loglike, chain, data)

    # test the dimensions: data points, samples, chains
    @test size(pw_lls_turing) == (50, 1000, 4)
    # test that sum of pointwise log likelihoods equals sum of log likelihoods
    turing_samples = Array(Chains(chain, :parameters).value)
    # make this more terse with @tullio or other method later
    LL = 0.0
    n_samples, n_parms, n_chains = size(turing_samples)
    for s in 1:n_samples
        for c in 1:n_chains
            LL += sum(logpdf.(Normal(turing_samples[s, :, c]...), data))
        end
    end
    @test LL ≈ sum(pw_lls_turing) atol = 1e-6
    # Turing should work the same as compute_loglike
    @test pw_lls_loglike ≈ pw_lls_turing atol = 1e-6

    # test that psis_loo works with Turing model and MCMCChains and yields correct type
    psis_loo_output = psis_loo(model(data), chain)
    @test isa(psis_loo_output, PsisLoo)
    # test that loo works with Turing model and MCMCChains and yields correct type
    psis_output = loo(model(data), chain)
    @test isa(psis_output, PsisLoo)
    # test that psis works with Turing model and MCMCChains and yields correct type
    psis_output = psis(model(data), chain)
    @test isa(psis_output, Psis)

    # ensure that methods work with r_eff argument
    r_eff = similar(pw_lls_turing, 0)
    # test that psis_loo works with Turin and gives correct type
    psis_loo_output = psis_loo(model(data), chain; r_eff=r_eff)
    @test isa(psis_loo_output, PsisLoo)
    # test that loo works with Turing model and MCMCChains and gives correct type
    psis_output = loo(model(data), chain; r_eff=r_eff)
    @test isa(psis_output, PsisLoo)
    # test that psis works with Turing model and MCMCChains and gives correct type
    psis_output = psis(model(data), chain; r_eff=r_eff)
    @test isa(psis_output, Psis)

end


text:
using ProgressMeter

# directories
examples_dir = joinpath(Pkg.dir("Graft"), "examples/")
dataset_dir = joinpath(examples_dir, "gplus")

vertex_file = joinpath(examples_dir, "vertex_data.txt")
edge_file = joinpath(examples_dir, "gplus_combined.txt")

# Change dir if required
if pwd() != examples_dir
   cd(examples_dir)
end

# Download and uncompress the full dataset
if !isdir(dataset_dir)
   download("http://snap.stanford.edu/data/gplus.tar.gz", joinpath(examples_dir, "gplus.tar.gz"))
   run(`tar xzvf gplus.tar.gz`)
end

# Download and uncompress the edge list
if !isfile(edge_file)
   download("http://snap.stanford.edu/data/gplus_combined.txt.gz", joinpath(examples_dir, "gplus_combined.txt.gz"))
   run(`gzip -d gplus_combined.txt.gz`)
end

# List of all ego nodes in dataset
nodes = unique(map(x->split(x, '.')[1], readdir(dataset_dir)))

# Outfile
out = open("vertex_data.txt", "w")

# Function to extract linenums and linevals from *.featnames
function extract_featnames(node)
   N = countlines("$(dataset_dir)/$(node).featnames")
   open("$(dataset_dir)/$(node).featnames", "r") do featnames

      prop_map = Dict(
         "gender" => 1,
         "institution" => 2,
         "job_title" => 3,
         "last_name" => 4,
         "place" => 5,
         "university" => 6
      )

      linenums = sizehint!(Int[], N)
      linevals = sizehint!(String[], N)

      for line in eachline(featnames)
         prop, value = match(r"\d+\s([\w_]+):(.*)", rstrip(line)).captures
         push!(linenums, prop_map[prop])
         push!(linevals, join(value))
      end

      linenums, linevals
   end
end

# Iterate over nodes and build the vertex dataset
p = Progress(length(nodes), 1)
for (count,node) in enumerate(nodes)
   # Fetch the data compression scheme
   linenums, linevals = extract_featnames(node)

   # Compute the size of the data array stored in the file
   M = countlines("$(dataset_dir)/$(node).feat")
   N = length(linenums)

   # Initialize data arrays
   labels = sizehint!(Int128[], N+1)
   genders = sizehint!(Int[], N+1)
   institutions = sizehint!(Vector{String}[], N+1)
   job_titles = sizehint!(Vector{String}[], N+1)
   last_names = sizehint!(String[], N+1)
   places = sizehint!(Vector{String}[], N+1)
   universities = sizehint!(Vector{String}[], N+1)

   # Extract alter-ego data
   feats = readdlm("$(dataset_dir)/$(node).feat", Int128, dims=(M,N+1))

   for i in 1 : size(feats, 1)
      push!(labels, feats[i,1])

      arr = feats[i,2:end]
      props = linenums .* arr

      gs = filter(x->length(x) > 0, linevals[find(x-> x == 1, props)])
      push!(genders, length(gs) == 0 ? 0 : parse(gs[1]))
      push!(institutions, filter(x->length(x) > 0, linevals[find(x-> x == 2, props)]))
      push!(job_titles, filter(x->length(x) > 0, linevals[find(x-> x == 3, props)]))
      push!(last_names, join(linevals[find(x-> x == 4, props)], ' '))
      push!(places, filter(x->length(x) > 0, linevals[find(x-> x == 5, props)]))
      push!(universities, filter(x->length(x) > 0, linevals[find(x-> x == 6, props)]))
   end

   # Extract ego-node data
   feats = map(parse, split(readline(open("$(dataset_dir)/$(node).egofeat"))))
   props = linenums .* feats

   push!(labels, parse(Int128, node))

   gs = filter(x->length(x) > 0, linevals[find(x-> x == 1, props)])
   push!(genders, length(gs) == 0 ? 0 : parse(gs[1]))
   push!(institutions, filter(x->length(x) > 0, linevals[find(x-> x == 2, props)]))
   push!(job_titles, filter(x->length(x) > 0, linevals[find(x-> x == 3, props)]))
   push!(last_names, join(linevals[find(x-> x == 4, props)], ' '))
   push!(places, filter(x->length(x) > 0, linevals[find(x-> x == 5, props)]))
   push!(universities, filter(x->length(x) > 0, linevals[find(x-> x == 6, props)]))

   # Join data arrays and write to file+
   vdata = hcat(labels, genders, institutions, job_titles, last_names, places, universities)
   writedlm(out, vdata)

   update!(p, count)
end

# flush buffer
close(out)

# Write the header
open("Graph.txt", "w") do gfile
   println(gfile, "107614\t13673453\tInt128")
   println(gfile, join(["gender", "institution", "job_title", "last_name", "place", "university"], '\t'))
   println(gfile, join(["Int", "Vector{String}", "Vector{String}", "String", "Vector{String}", "Vector{String}"], '\t'))
   println(gfile)
   println(gfile)
end


text:
########################### Common/Utils ###########################

#Convenience struct to pass `conn` each time instead of the parts
mutable struct OmniSciConnection
    session::TSessionId
    c::MapDClient
end

#convert methods implicitly needed for squashbitmask
convert(::Type{DateTime}, x::Int64) =  Dates.unix2datetime(x)
convert(::Type{Date}, x::Int64) =  Date(Dates.unix2datetime(x))
convert(::Type{Time}, x::Int64) = Time(x/3600, x % 3600)

#conversions from WKT to GeoInterface types, to get typed df from sql_execute
convert(::Type{GeoInterface.Point}, x::String) = GeoInterface.Point(LibGEOS.readgeom(x))
convert(::Type{GeoInterface.LineString}, x::String) = GeoInterface.LineString(LibGEOS.readgeom(x))
convert(::Type{GeoInterface.Polygon}, x::String) = GeoInterface.Polygon(LibGEOS.readgeom(x))
convert(::Type{GeoInterface.MultiPolygon}, x::String) = GeoInterface.MultiPolygon(LibGEOS.readgeom(x))

#WKT from LibGEOS types
wkt(x::GeoInterface.AbstractPoint) =  writegeom(LibGEOS.Point(x))
wkt(x::GeoInterface.AbstractLineString) = writegeom(LibGEOS.LineString(x))
wkt(x::GeoInterface.AbstractPolygon) = writegeom(LibGEOS.Polygon(x))
wkt(x::GeoInterface.AbstractMultiPolygon) = writegeom(LibGEOS.MultiPolygon(x))
wkt(x::T) where T <: Union{String, Missing} = ""

#Define these methods to avoid type piracy
myDateTime(x::Missing) = missing
myDateTime(x) = DateTime(x)

mydatetime2unix(x::Missing) = missing
mydatetime2unix(x) = datetime2unix(floor(x, Dates.Second))

epochdays(x::Missing) = missing
epochdays(x) = (x - Dates.Date("1970-01-01")).value

myInt64(x::Missing) = missing
myInt64(x) = Int64(x)

#Find which field in the struct the data actually is
function findvalues(x::OmniSci.TColumn)
    for f in propertynames(x.data)
        n = length(getfield(x.data, f))
        if n > 0
            return (f, n)
        end
    end
end

#Take two vectors, values and nulls, make into a single vector
#TODO: Figure out if its possible to further minimize allocations
function squashbitmask(x::TColumn, typeinfo::Tuple{DataType, Bool})

    #Get location of data from struct, eltype of vector and its length
    valuescol, n = findvalues(x)

    #unpack type info
    ltype, nullable = typeinfo

    #Build/fill new vector based on missingness
    if nullable
        A = Vector{Union{ltype, Missing}}(undef, n)
        @simd for i = 1:n
            @inbounds A[i] = ifelse(x.nulls[i], missing, getfield(x.data, valuescol)[i])
        end
    else
        #Assumption here is that for columns without nulls, can cast directly
        #Since whatever is in the x.data column should be valid values
        #And appropriate convert methods defined
        A = convert(Vector{ltype}, getfield(x.data, valuescol))
    end

    return A
end

seconds_since_midnight(x::Time) = (hour(x) * 3600) + (minute(x) * 60) + second(x)
seconds_since_midnight(x::Missing) = missing

#TODO: ensure valid col names
#1. check if first character a number, prepend some valid string
#2. check if in reserved word list, append _ if reserved
function sanitizecolnames(x::Symbol)
    str = string(x)
end

#convert julia types to OmniSci types
function getsqlcoltype(x, precision::Tuple{Int, Int})

    dec = "DECIMAL$precision"

    lookup = Dict(
        #ints
        Union{Missing, Int8} => "TINYINT",
        Int8 => "TINYINT",
        Union{Missing, Int16} => "SMALLINT",
        Int16 => "SMALLINT",
        Union{Missing, Int32} => "INTEGER",
        Int32 => "INTEGER",
        Union{Missing, Int64} => "BIGINT",
        Int64 => "BIGINT",
        #floats
        Union{Missing, Float32} => "FLOAT",
        Float32 => "FLOAT",
        Union{Missing, Float64} => "DOUBLE",
        Float64 => "DOUBLE",
        #rational: convert to float, rational not avail in OmniSci
        Union{Missing, Rational} => "FLOAT",
        Rational => "FLOAT",
        #strings
        Union{String, Missing} => "TEXT ENCODING DICT",
        String => "TEXT ENCODING DICT",
        #bool
        Union{Missing, Bool} => "BOOLEAN",
        Bool => "BOOLEAN",
        #dates and time
        Union{Missing, Date} => "DATE",
        Date => "DATE",
        Union{Missing, Time} => "TIME",
        Time => "TIME",
        Union{Missing, DateTime} => "TIMESTAMP",
        DateTime => "TIMESTAMP",
        #geospatial
        Union{Missing, GeoInterface.Point} => "POINT",
        GeoInterface.Point => "POINT",
        Union{Missing, GeoInterface.LineString} => "LINESTRING",
        GeoInterface.LineString => "LINESTRING",
        Union{Missing, GeoInterface.Polygon} => "POLYGON",
        GeoInterface.Polygon => "POLYGON",
        Union{Missing, GeoInterface.MultiPolygon} => "MULTIPOLYGON",
        GeoInterface.MultiPolygon => "MULTIPOLYGON",
        #arrays
        Array{Union{Missing, Int8},1} => "TINYINT[]",
        Array{Int8,1} => "TINYINT[]",
        Array{Union{Missing, Int16},1} => "SMALLINT[]",
        Array{Int16,1} => "SMALLINT[]",
        Array{Union{Missing, Int32},1} => "INTEGER[]",
        Array{Int32,1} => "INTEGER[]",
        Array{Union{Missing, Int64},1} => "BIGINT[]",
        Array{Int64,1} => "BIGINT[]",
        Array{String,1} => "TEXT[]",
        Array{Union{Missing, Float32},1} => "FLOAT[]",
        Array{Float32,1} => "FLOAT[]",
        Array{Union{Missing, Float64},1} => "DOUBLE[]",
        Array{Float64,1} => "DOUBLE[]",
        Array{Union{Missing, Bool},1} => "BOOLEAN[]",
        Array{Bool,1} => "BOOLEAN[]",
        Array{Union{Missing, Date},1} => "DATE[]",
        Array{Date,1} => "DATE[]",
        Array{Union{Missing, Time},1} => "TIME[]",
        Array{Time,1} => "TIME[]",
        Array{Union{Missing, DateTime},1} => "TIMESTAMP[]",
        Array{DateTime,1} => "TIMESTAMP[]",
        #DecFP
        Union{DecFP.Dec32, Missing} => dec,
        DecFP.Dec32 => dec,
        Union{DecFP.Dec64, Missing} => dec,
        DecFP.Dec64 => dec,
        Union{DecFP.Dec128, Missing} => dec,
        DecFP.Dec128 => dec,
        #Decimals.jl
        Union{Decimals.Decimal, Missing} => dec,
        Decimals.Decimal => dec,
        #missing: default to TEXT ENCODING DICT as design decision
        #https://github.com/omnisci/OmniSci.jl/issues/76
        Missing => "TEXT ENCODING DICT"
    )

    get(lookup, x, "Unknown")

end

########################### Typedefs for load_table method ###########################

#For functions below, value for is_null should be known based on the dispatched type
#Left as keyword just in case my assumption incorrect

# convert vectors to string representations for atomic types only
function TStringValue(str_val::Vector{<:Union{Real, String, Char, TimeType, Missing}}, is_null::Bool = false)
  val = OmniSci.TStringValue()

  #Write values into buffer to avoid any weird display issues
  io = IOBuffer()
  write(io, "{")

  for val in str_val[1:end-1]
    ismissing(val) ? val_ = "NA" : val_ = val
    write(io, string(val_))
    write(io, ",")
  end

  #for last value in array, don't add trailing comma
  ismissing(str_val[end]) ? lastval = "NA" : lastval = str_val[end]
  write(io, string(lastval))
  write(io, "}")

  p = String(take!(io))

  Thrift.set_field!(val, :str_val, p)
  Thrift.set_field!(val, :is_null, is_null)
  return val
end

function TStringValue(str_val::T, is_null::Bool = false) where T <: Union{GeoInterface.AbstractLineString, GeoInterface.AbstractPoint, GeoInterface.AbstractPolygon, GeoInterface.AbstractMultiPolygon}
  val = OmniSci.TStringValue()
  p = wkt(str_val)
  Thrift.set_field!(val, :str_val, p)
  Thrift.set_field!(val, :is_null, is_null)
  return val
end

function TStringValue(str_val::Rational, is_null::Bool = false)
  val = OmniSci.TStringValue()
  Thrift.set_field!(val, :str_val, string(convert(Float64, str_val)))
  Thrift.set_field!(val, :is_null, is_null)
  return val
end

function TStringValue(str_val::T, is_null::Bool = true) where T <: Union{Missing, Nothing}
  val = OmniSci.TStringValue()
  Thrift.set_field!(val, :str_val, string(str_val))
  Thrift.set_field!(val, :is_null, is_null)
  return val
end

#generic fallback for any type serializable with string()
function TStringValue(str_val, is_null::Bool = false)
  val = OmniSci.TStringValue()
  Thrift.set_field!(val, :str_val, string(str_val))
  Thrift.set_field!(val, :is_null, is_null)
  return val
end

function TStringRow(cols::Vector{TStringValue})
    tsr = OmniSci.TStringRow()
    Thrift.set_field!(tsr, :cols, cols)
    return tsr
end

TStringRow(cols::AbstractVector) = TStringRow(TStringValue.(cols))
TStringRow(cols::DataFrameRow{DataFrame}) = TStringRow([TStringValue(x) for x in cols])

########################### Typedefs for load_table_binary_columnar method ###########################

function TColumn(x::AbstractVector{<:Union{Missing, T}}) where T <: Union{Int8, Int16, Int32, Int64}

    #Create TColumn, fill nulls column by checking for missingness
    tc = TColumn()
    Thrift.set_field!(tc, :nulls, convert(Vector{Bool}, ismissing.(x)))

    #Replace missing values with typed sentinel and convert to Vector{Int64} per API requirement
    tcd = TColumnData()
    Thrift.set_field!(tcd, :int_col, convert(Vector{Int64}, coalesce.(x, -1)))

    #Complete TColumn
    Thrift.set_field!(tc, :data, tcd)

    return tc
end

function TColumn(x::AbstractVector{<:Union{Missing, T}}) where T <: Union{AbstractFloat, Rational}

    #Create TColumn, fill nulls column by checking for missingness
    tc = TColumn()
    Thrift.set_field!(tc, :nulls, convert(Vector{Bool}, ismissing.(x)))

    #Replace missing values with typed sentinel and convert to Vector{Float64} per API requirement
    tcd = TColumnData()
    Thrift.set_field!(tcd, :real_col, convert(Vector{Float64}, coalesce.(x, -1.0)))

    #Complete TColumn
    Thrift.set_field!(tc, :data, tcd)

    return tc
end

function TColumn(x::AbstractVector{<:Union{Missing, AbstractString}})

    #Create TColumn, fill nulls column by checking for missingness
    tc = TColumn()
    Thrift.set_field!(tc, :nulls, convert(Vector{Bool}, ismissing.(x)))

    #Replace missing values with typed sentinel and convert to Vector{String} per API requirement
    tcd = TColumnData()
    Thrift.set_field!(tcd, :str_col, convert(Vector{String}, coalesce.(x, "")))

    #Complete TColumn
    Thrift.set_field!(tc, :data, tcd)

    return tc
end

function TColumn(x::AbstractVector{<:Union{Missing, Bool}})

    #Create TColumn, fill nulls column by checking for missingness
    tc = TColumn()
    Thrift.set_field!(tc, :nulls, convert(Vector{Bool}, ismissing.(x)))

    #Replace missing values with typed sentinel and convert to Vector{Int64} per API requirement
    tcd = TColumnData()
    Thrift.set_field!(tcd, :int_col, convert(Vector{Int64}, coalesce.(x, -1)))

    #Complete TColumn
    Thrift.set_field!(tc, :data, tcd)

    return tc
end

function TColumn(x::AbstractVector{<:Union{Missing, T}}) where T <: Union{GeoInterface.AbstractPoint, GeoInterface.AbstractPolygon, GeoInterface.AbstractLineString, GeoInterface.AbstractMultiPolygon}

    #Create TColumn, fill nulls column by checking for missingness
    tc = TColumn()
    Thrift.set_field!(tc, :nulls, convert(Vector{Bool}, ismissing.(x)))

    #Convert geo types to string
    tcd = TColumnData()
    Thrift.set_field!(tcd, :str_col, wkt.(x))

    #Complete TColumn
    Thrift.set_field!(tc, :data, tcd)

    return tc
end

#for a missing column, default to empty string column per https://github.com/omnisci/OmniSci.jl/pull/77
function TColumn(x::AbstractVector{Missing})

    #Create TColumn, fill nulls column by checking for missingness
    tc = TColumn()
    Thrift.set_field!(tc, :nulls, convert(Vector{Bool}, ismissing.(x)))

    #Replace missing values with typed sentinel and convert to Vector{String} per API requirement
    tcd = TColumnData()
    Thrift.set_field!(tcd, :str_col, convert(Vector{String}, coalesce.(x, "")))

    #Complete TColumn
    Thrift.set_field!(tc, :data, tcd)

    return tc
end

# Dispatches to Int after conversion function applied
TColumn(x::AbstractVector{<:Union{Missing, DateTime}}) = TColumn(myInt64.(mydatetime2unix.(x)))
TColumn(x::AbstractVector{<:Union{Missing, Time}}) = TColumn(seconds_since_midnight.(x))

#Before OmniSci 4.4, Dates were specified in epoch seconds
#TODO: consider making a warning about data loading?
TColumn(x::AbstractVector{<:Union{Missing, Date}}) = TColumn(myDateTime.(x))