text:
using Test
using TightBindingLattice

using TightBindingLattice: simplify_name

using LinearAlgebra
using YAML

@testset "PointSymmetry" begin
    @testset "simplify_name" begin
        @test simplify_name("1") == "1"
        @test simplify_name("-1") == "-1"
        @test simplify_name("2011") == "2"
        @test simplify_name("3+011") == "3"
        @test simplify_name("m01") == "m"
    end

    @testset "constructor failure" begin
        group = FiniteGroup([1 2; 2 1])
        generators = [2]
        conjugacy_classes = [[1], [2]]
        character_table = [1 1; 1 -1]
        irreps = [
                [ones(ComplexF64, 1, 1), ones(ComplexF64, 1, 1)],
                [ones(ComplexF64, 1, 1), -ones(ComplexF64, 1, 1)],
        ]
        element_names = ["1", "-1"]
        matrix_representations = [[1 0; 0 1], [-1 0; 0 -1]]
        hermann_mauguin = "-1"
        schoenflies = "Ci"
        psym = PointSymmetry(group, generators,
                             conjugacy_classes, character_table, irreps,
                             element_names, matrix_representations,
                             hermann_mauguin, schoenflies)

        let TBL = TightBindingLattice
            @test eltype(psym) == PointOperation{Int}
            @test valtype(psym) == PointOperation{Int}
            @test all(TBL.element(psym, i) == PointOperation(matrix_representations[i]) for i in 1:2)
            @test TBL.elements(psym) == PointOperation.(matrix_representations)
            @test all(TBL.element_name(psym, i) == element_names[i] for i in 1:2)
            @test TBL.element_names(psym) == element_names
            @test TBL.group(psym) == group
            @test TBL.group_order(psym) == group_order(group)
            @test TBL.group_multiplication_table(psym) == group.multiplication_table
            @test TBL.character_table(psym) == character_table
            @test TBL.irreps(psym) == irreps
            @test all(TBL.irrep(psym, i) == irreps[i] for i in 1:2)
            @test num_irreps(psym) == 2
            @test all(TBL.irrep_dimension(psym, i) == 1 for i in 1:2)
            @test generator_indices(psym) == [2]
            @test generator_elements(psym) == [PointOperation([-1 0; 0 -1])]
        end

        let generators = [1]
            @test_throws ArgumentError PointSymmetry(group, generators,
                        conjugacy_classes, character_table, irreps,
                        element_names, matrix_representations, hermann_mauguin, schoenflies)
        end
        let conjugacy_classes = [[1], [1,2]]
            @test_throws ArgumentError PointSymmetry(group, generators,
                        conjugacy_classes, character_table, irreps,
                        element_names, matrix_representations, hermann_mauguin, schoenflies)
        end
        let character_table = [1 2 3; 4 5 6; 7 8 9]
            @test_throws ArgumentError PointSymmetry(group, generators,
                        conjugacy_classes, character_table, irreps,
                        element_names, matrix_representations, hermann_mauguin, schoenflies)
        end
        let irreps = [
                    [ones(ComplexF64, 1, 1), ones(ComplexF64, 1, 1)],
            ]
            @test_throws ArgumentError PointSymmetry(group, generators,
                        conjugacy_classes, character_table, irreps,
                        element_names, matrix_representations, hermann_mauguin, schoenflies)
        end
        let irreps = [
                    [ones(ComplexF64, 1, 1), ones(ComplexF64, 1, 1)],
                    [ones(ComplexF64, 1, 1)],
            ]
            @test_throws ArgumentError PointSymmetry(group, generators,
                        conjugacy_classes, character_table, irreps,
                        element_names, matrix_representations, hermann_mauguin, schoenflies)
        end
        let irreps = [
                    [ones(ComplexF64, 1, 1), ones(ComplexF64, 1, 1)],
                    [ones(ComplexF64, 1, 1), -ones(ComplexF64, 2, 2)],
            ]
            @test_throws ArgumentError PointSymmetry(group, generators,
                        conjugacy_classes, character_table, irreps,
                        element_names, matrix_representations, hermann_mauguin, schoenflies)
        end
        let irreps = [
                    [ones(ComplexF64, 1, 1), ones(ComplexF64, 1, 1)],
                    [2*ones(ComplexF64, 1, 1), -ones(ComplexF64, 1, 1)],
            ]
            @test_throws ArgumentError PointSymmetry(group, generators,
                        conjugacy_classes, character_table, irreps,
                        element_names, matrix_representations, hermann_mauguin, schoenflies)
        end
        let irreps = [
                    [ones(ComplexF64, 1, 1), ones(ComplexF64, 1, 1)],
                    [ones(ComplexF64, 1, 1), 2*ones(ComplexF64, 1, 1)],
            ]
            @test_throws ArgumentError PointSymmetry(group, generators,
                        conjugacy_classes, character_table, irreps,
                        element_names, matrix_representations, hermann_mauguin, schoenflies)
        end
        let element_names = ["1", "2", "3"]
            @test_throws ArgumentError PointSymmetry(group, generators,
                        conjugacy_classes, character_table, irreps,
                        element_names, matrix_representations, hermann_mauguin, schoenflies)
        end
        let matrix_representations = [[1 0; 0 1], [-1 0 0; 0 -1 0]]
            @test_throws ArgumentError PointSymmetry(group, generators,
                        conjugacy_classes, character_table, irreps,
                        element_names, matrix_representations, hermann_mauguin, schoenflies)
        end
        let matrix_representations = [[1 0 0; 1 0 0], [-1 0 0; -1 0 0]]
            @test_throws ArgumentError PointSymmetry(group, generators,
                        conjugacy_classes, character_table, irreps,
                        element_names, matrix_representations, hermann_mauguin, schoenflies)
        end
        let matrix_representations = [[1 0 0; 1 0 0]]
            @test_throws ArgumentError PointSymmetry(group, generators,
                        conjugacy_classes, character_table, irreps,
                        element_names, matrix_representations, hermann_mauguin, schoenflies)
        end
        let matrix_representations = [[1 0 ; 0 1], [1 0; 0 1]]
            @test_throws ArgumentError PointSymmetry(group, generators,
                        conjugacy_classes, character_table, irreps,
                        element_names, matrix_representations, hermann_mauguin, schoenflies)
        end
    end

    # D_4 (422)
    file_path = abspath(@__DIR__, "..", "data", "PointGroup3D", "PointGroup3D-12.yaml")
    data_yaml = YAML.load_file(file_path)

    let data_yaml_2 = deepcopy(data_yaml)
        data_yaml_2["Generators"] = [1]
        @test_throws ArgumentError read_point_symmetry(data_yaml_2)
    end

    psym = read_point_symmetry(data_yaml)
    ord_group = group_order(psym.group)

    @testset "elements" begin
        @test group_order(psym) == 8
        for i in 1:ord_group
            @test element_names(psym)[i] == element_name(psym, i)
        end
    end

    @testset "generators" begin
        # Generators completely generate the group
        @test generate_subgroup(psym.group, psym.generators) == BitSet(1:ord_group)
        @test prod(psym.group.period_lengths[x] for x in psym.generators) != ord_group # nonabelian
    end

    @testset "multiplication table" begin
        matrep_lookup = Dict(m=>i for (i, m) in enumerate(psym.matrix_representations))
        matrep_mtab = zeros(Int, (ord_group, ord_group))
        for i1 in 1:ord_group, i2 in 1:ord_group
            m1 = psym.matrix_representations[i1]
            m2 = psym.matrix_representations[i2]
            m3 = m1 * m2
            i3 = matrep_lookup[m3]
            matrep_mtab[i1, i2] = i3
        end
        @test matrep_mtab == psym.group.multiplication_table
        @test matrep_mtab == group_multiplication_table(psym)
        @test matrep_mtab == group_multiplication_table(psym.group)
    end

    @testset "symmetry_product" begin
        p = symmetry_product(psym)
        for x in elements(psym), y in elements(psym)
            z = x*y
            @test isa(z, PointOperation{Int})
            @test z.matrix == x.matrix * y.matrix
        end
    end

    @testset "iterate" begin
        @test [x for x in psym] == collect(elements(psym))
        @test all(x ∈ psym for x in psym)
        @test 100 ∉ psym
        @test PointOperation([1 2 0; -1 1 0; 0 0 1]) ∉ psym
        @test IdentityOperation(Int, 2) ∉ psym
        @test IdentityOperation(Int, 3) ∈ psym
        @test TranslationOperation([100, 100, 0]) ∉ psym
        @test TranslationOperation([0, 0]) ∉ psym
        @test TranslationOperation([0, 0, 0]) ∈ psym
        @test PointOperation([1 0 0; 0 1 0; 0 0 1]) ∈ psym
        @test PointOperation([0 1 0; 1 0 0; 0 0 -1]) ∈ psym
        @test SpaceOperation([0 1 0; 1 0 0; 0 0 -1], [0, 0, 0]) ∈ psym
        @test SpaceOperation([0 1 0; 1 0 0; 0 0 -1], [1, 0, 0]) ∉ psym
    end

    @testset "irreps" begin
        @test num_irreps(psym) == 5
        let χ = [1  1  1  1  1;
                 1  1  1 -1 -1;
                 1  1 -1 -1  1;
                 1  1 -1  1 -1;
                 2 -2  0  0  0]
            @test character_table(psym) ≈ χ
            @test size(character_table(psym)) == (5, 5)
        end

        @test length(irreps(psym)) == 5
        for i in 1:5
            @test irreps(psym)[i] == irrep(psym, i)
        end

        for (idx_irrep, irrep) in enumerate(irreps(psym))
            d = psym.character_table[idx_irrep, 1]
            @test irrep_dimension(psym, idx_irrep) == d
            for m in irrep
                @test size(m) == (d, d)
            end # for irrep
            for (idx_cc, cc) in enumerate(psym.conjugacy_classes)
                character = character_table(psym)[idx_irrep, idx_cc]
                for idx_elem in cc
                    @test isapprox(tr(irrep[idx_elem]), character; atol=Base.rtoldefault(Float64))
                end # for cc
            end # for enumerate(psym.conjugacy_classes)
        end # for enumerate(irreps(psym))
    end

    @testset "project" begin
        project(psym, [1 0 0; 0 1 0])
        @test_throws ArgumentError project(psym, [1 0; 0 1])
        @test_throws ArgumentError project(psym, [1 1 1; 1 1 1])
        @test_throws ArgumentError project(psym, [1 0 0;]) # projection not faithful
    end

    @testset "iscompatible" begin
        hc1 = OrthoCube([4 0; 0 4])
        hc2 = OrthoCube([4 0; 0 3])
        tsym1 = TranslationSymmetry(hc1)
        tsym2 = TranslationSymmetry(hc2)

        @test_throws DimensionMismatch iscompatible(hc1, psym)
        @test_throws DimensionMismatch iscompatible(tsym1, psym)

        psym_proj = project(psym, [1 0 0; 0 1 0])

        @test iscompatible(hc1, psym_proj)
        @test !iscompatible(hc2, psym_proj)
        @test iscompatible(tsym1, psym_proj)
        @test !iscompatible(tsym2, psym_proj)
        @test little_symmetry(tsym1, psym_proj).hermann_mauguin == "422"
        @test little_symmetry(tsym2, psym_proj).hermann_mauguin == "222"
        @test little_symmetry(TranslationSymmetry([4 1; 0 3]), psym_proj).hermann_mauguin == "2"
        @test_throws ArgumentError little_symmetry(tsym2, 1, psym_proj) # when specifying irrep, tsym and psym have to be compatible
    end

    @testset "symmetry_name" begin
        let n = lowercase(symmetry_name(psym))
            @test occursin("point", n)
            @test occursin("422", n)
        end
    end

    @testset "two-band model" begin
        unitcell = make_unitcell([1.0 0.0; 0.0 1.0]; SiteType=String)
        addsite!(unitcell, "Ox", FractCoord([0,0], [0.5, 0.0]))
        addsite!(unitcell, "Oy", FractCoord([0,0], [0.0, 0.5]))

        psym = project(PointSymmetryDatabase.get(13), [1 0 0; 0 1 0])

        @testset "lattice permutations" begin
            idx_C4 = 3
            @test element_name(psym, idx_C4) == "4+001"
            @test findsitemap(unitcell, element(psym,idx_C4)) == [(2, [0,0]), (1, [-1,0])]
            @test findsitemap(unitcell, psym)[idx_C4] == [(2, [0,0]), (1, [-1,0])]

            lattice = make_lattice(unitcell, [2 0; 0 2])
            psymbed = embed(lattice, psym)
            tsym = TranslationSymmetry(lattice)
            tsymbed = embed(lattice, tsym)

            perms = [embed(lattice, op) for op in elements(psym)]
            @test length(perms) == length(psym.element_names)
            @test perms[1] == SitePermutation(1:8) # identity

            perms2 = SitePermutation[]
            for pop in elements(psym)
                push!(perms2, embed(lattice, pop))
                #push!(perms2, SitePermutation(get_site_permutation(lattice, mat, map)))
            end
            @test perms == perms2

            # |       |                 |       |
            # 6       8                 3       7
            # |       |         C4      |       |
            # . - 5 - . - 7 -   =>      . - 8 - . - 4 -
            # |       |                 |       |
            # 2       4                 1       5
            # |       |                 |       |
            # o - 1 - . - 3 -           o - 6 - . - 2 -
            @test perms[idx_C4] == SitePermutation([2,3,6,7,4,1,8,5])

        end # testset "lattice permutations"

        @testset "little group" begin
            lattice = make_lattice(unitcell, [2 0; 0 2])
            tsym = TranslationSymmetry(lattice)
            @test little_group_elements(tsym, 1, psym) == collect(1:8)
            @test little_group_elements(tsym, 2, psym) == [1,2,5,6]
            @test little_group_elements(tsym, 3, psym) == [1,2,5,6]
            @test little_group_elements(tsym, 4, psym) == collect(1:8)

            @test little_group(tsym, 1, psym) == little_group(tsym, psym, 1:8)
            @test little_group(tsym, 2, psym) == little_group(tsym, psym, [1,2,5,6])
            @test little_group(tsym, 3, psym) == little_group(tsym, psym, [1,2,5,6])
            @test little_group(tsym, 4, psym) == little_group(tsym, psym, 1:8)
        end

        @testset "little_symmetry" begin
            for LSYM in [little_symmetry, little_symmetry_iso]
                lattice = make_lattice(unitcell, [4 0; 0 4])
                tsym = TranslationSymmetry(lattice)
                for tsym_irrep in 1:num_irreps(tsym)
                    psym_little = LSYM(tsym, tsym_irrep, psym)
                    k = tsym.coordinates[tsym_irrep]
                    @test iscompatible(tsym, tsym_irrep, psym) == (k in [[0,0], [2,2]])
                    @test iscompatible(tsym, tsym_irrep, psym_little)
                    lg_matrep = psym.matrix_representations[little_group_elements(tsym, tsym_irrep, psym)]
                    @test !isnothing(group_isomorphism(little_group(tsym, tsym_irrep, psym),
                                                       FiniteGroup(group_multiplication_table(lg_matrep))
                                                       ))
                end # for tsym_irrep
            end
        end # testset little_symmetry
    end
end # @testset "PointSymmetry"

@testset "PointSymmetryDatabase" begin
    psym1 = PointSymmetryDatabase.get(13)
    psym2 = PointSymmetryDatabase.find("4mm")
    @test length(psym1) == length(psym2) == 8
    @test psym1.hermann_mauguin == "4mm"
    @test psym2.hermann_mauguin == "4mm"
    @test_throws ArgumentError PointSymmetryDatabase.get(999)
    @test isnothing(PointSymmetryDatabase.find("blah blah"))
    @test dimension(psym1) == dimension(psym2) == 3

    psym1p = PointSymmetryDatabase.get3d(13)
    psym2p = PointSymmetryDatabase.find3d("4mm")
    @test length(psym1p) == length(psym2p) == 8
    @test psym1p.hermann_mauguin == "4mm"
    @test psym2p.hermann_mauguin == "4mm"
    @test_throws ArgumentError PointSymmetryDatabase.get3d(999)
    @test isnothing(PointSymmetryDatabase.find3d("blah blah"))
    @test dimension(psym1p) == dimension(psym2p) == 3

    psym1p = PointSymmetryDatabase.get2d(6)
    psym2p = PointSymmetryDatabase.find2d("4mm")
    @test length(psym1p) == length(psym2p) == 8
    @test psym1p.hermann_mauguin == "4mm"
    @test psym2p.hermann_mauguin == "4mm"
    @test_throws ArgumentError PointSymmetryDatabase.get2d(999)
    @test isnothing(PointSymmetryDatabase.find2d("blah blah"))
    @test dimension(psym1p) == dimension(psym2p) == 2
end


text:
using ArgParse
using DrWatson
@quickactivate
using BSON
using Random
using FileIO
using DataFrames
using Base.Threads: @threads
using GenerativeAD

# pkgs which come from deserialized BSONs
# have to be present in the Main module
using ValueHistories
using LinearAlgebra

s = ArgParseSettings()
@add_arg_table! s begin
    "source_prefix"
		arg_type = String
		default = "experiments/images"
		help = "Data prefix of experiment files."
	"target_prefix"
		arg_type = String
		default = "evaluation/images"
		help = "Data prefix of generated files."
	"-f", "--force"
    	action = :store_true
		help = "Overwrite all generated files."
end

"""
	generate_stats(source_prefix::String, target_prefix::String; force=true)

Collects all the results from experiments in datadir prefix `source_prefix`,
computes evaluation metrics and stores results in datadir prefix `target_prefix`
while retaining the folder structure. If `force=true` the function overwrites
already precomputed results.
"""
function generate_stats(source_prefix::String, target_prefix::String; force=true)
	(source_prefix == target_prefix) && error("Results have to be stored in different folder.")

	source = datadir(source_prefix)
	@info "Collecting files from $source folder."
	files = GenerativeAD.Evaluation.collect_files_th(source)
	# filter out model files
	filter!(x -> !startswith(basename(x), "model"), files)
	@info "Collected $(length(files)) files from $source folder."
	# it might happen that when appending results some of the cores
	# just go over already computed files
	files = files[randperm(length(files))]

	@threads for f in files
		target_dir = dirname(replace(f, source_prefix => target_prefix))
		target = joinpath(target_dir, "eval_$(basename(f))")
		try
			if (isfile(target) && force) || ~isfile(target)
				df = GenerativeAD.Evaluation.compute_quantile_stats(f)
				wsave(target, Dict(:df => df))
				@info "Saving evaluation results at $(target)"
			end
		catch e
			# remove old files in order to ensure consistency
			if (isfile(target) && force)
				rm(target)
			end
			@warn "Processing of $f failed due to \n$e"
		end
	end
end

function main(args)
	generate_stats(
		args["source_prefix"],
		args["target_prefix"];
		force=args["force"])
	@info "---------------- DONE -----------------"
end

main(parse_args(ARGS, s))

ext: jl


text:
# Generate list of LaTeX tab-completions supported by the Julia REPL
# as documented in doc/manual/unicode-input-table.rst
# The output will be rendered as a reStructuredText document to STDOUT
# Note: This script will download a file called UnicodeData.txt from the Unicode
# Consortium

include("../base/latex_symbols.jl")
include("../base/emoji_symbols.jl")

# Create list of different tab-completions for a given character
# Sometimes there is more than one way
vals = Dict()
for symbols in [latex_symbols, emoji_symbols], (k, v) in symbols
    vals[v] = push!(get!(vals, v, String[]), "\\"*k)
end

# Join with Unicode names to aid in lookup
isfile("UnicodeData.txt") || download(
    "http://www.unicode.org/Public/UCD/latest/ucd/UnicodeData.txt", "UnicodeData.txt")
unicodenames = Dict()
open("UnicodeData.txt") do unidata
    for line in readlines(unidata)
        tokens = split(line, ";")
        codepoint = parse(UInt32, "0x"tokens[1])
        name = tokens[ 2]=="" ? tokens[11] :
               tokens[11]=="" ? tokens[ 2] :
               tokens[ 2]*" / "*tokens[11]
        unicodenames[codepoint] = name
    end
end

# Render list
# Need to do this in two passes since ReST complains if the tables aren't exactly aligned
# Pass 1. Generate strings
entries = Any[("Code point(s)", "Character(s)", "Tab completion sequence(s)", "Unicode name(s)")]
maxlen = [map(length, entries[1])...]

# Prepend a non-breakable space to combining characters
# (ranges obtained from https://en.wikipedia.org/wiki/Combining_character)
function fix_combining_chars(uni)
    u = uni[1]
    isc = ('\u0300' ≤ u ≤ '\u036F') ||
          ('\u1AB0' ≤ u ≤ '\u1AFF') ||
          ('\u1DC0' ≤ u ≤ '\u1DFF') ||
          ('\u20D0' ≤ u ≤ '\u20FF') ||
          ('\uFE20' ≤ u ≤ '\uFE2F')
    return isc ? string('\u00A0', uni) : uni
end

for (chars, inputs) in sort!([x for x in vals], by=first)
    # Find all keys with this value
    entry = (
            join(map(c->"U+"*uppercase(hex(c, 5)), collect(chars)), " + "),
            fix_combining_chars(chars),
            join(inputs, ", "),
            join(map(c->get(unicodenames, UInt32(c), "(No Unicode name)"), collect(chars)), " + ")
        )

    currentlength = map(length, entry)
    for i=1:length(entry)
        maxlen[i] = max(maxlen[i], currentlength[i])
    end

    push!(entries, entry)
end

# Pass 2. Print table in ReST simple table format
function underline(str, maxlen)
    join(map(n->str^n, maxlen), " ")
end

isheader = true
println(underline("=", maxlen))
for entry in entries
    thisline = ""
    for (i, col) in enumerate(entry)
        thisline *= rpad(col, maxlen[i], " ") * " "

        # Hack round JuliaLang/julia#10825
        if i==2 && any(x->charwidth(x)==2, collect(col))
            thisline *=" "
        end
    end
    println(rstrip(thisline))

    if isheader
        println(underline("-", maxlen))
        isheader = false
    end
end
println(underline("=", maxlen))

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
using Dagger
Dagger.use_shared_array[] = false

f = i -> (i+1)
g = (i, j) -> (i + j)
h = (i, j, k) -> floor(Int, mean([i,j,k]))

const T = Any[Thunk(()->10), Thunk(()->10), Thunk(()->10)]

function gendag(U,V,W)
    fd = delayed(f)(U)
    fg = delayed(g)(fd, V)
    fh = delayed(h)(fd, fg, W)
    fh
end

function gendag(T)
    f1 = gendag(T[end-2:end]...)
    f2 = gendag(T[end-2:end]...)
    f3 = gendag(T[end-2:end]...)
    push!(T, f1)
    push!(T, f2)
    push!(T, f3)
    gendag(f1, f2, f3)
end

function pushthunk(t, a = Any[])
   isa(t, Thunk) && map(x->pushthunk(x, a), t.inputs)
   (t in a) || push!(a, t)
   a
end

function gen_cross_dag()
    x = gendag(T)
    for idx in 1:5
        x = gendag(T)
    end
    x
end

show_dag(dag) = println(show_plan(dag))

function gen_straight_dag(inp)
    nextinp = []
    if length(inp) == 6
        push!(nextinp, delayed(f)(inp[1]))
        push!(nextinp, delayed(g)(inp[2], inp[3]))
        push!(nextinp, delayed(h)(inp[4], inp[5], inp[6]))
        return delayed(h)(nextinp...)
    else
        for idx in 1:6:length(inp)
            f1 = delayed(f)(inp[idx])
            f2 = delayed(g)(inp[idx+1], inp[idx+2])
            f3 = delayed(h)(inp[idx+3], inp[idx+4], inp[idx+5])
            push!(nextinp, delayed(h)(f1, f2, f3))
        end
        return gen_straight_dag(nextinp)
    end
end

#const deepdag = gen_straight_dag(T)
#const ordereddag = pushthunk(deepdag)

import Dagger: dsort_chunks

function gen_sort_dag(totalsz=10^6, nchunks=40, batchsize=4, nresultchunks=1)
    cs = compute(rand(Blocks(ceil(Int,totalsz/nchunks)), totalsz)).chunks
    lt = Base.isless
    by = identity
    rev = false
    ord = Base.Sort.ord(lt,by,rev,Dagger.default_ord)

    dag1 = dsort_chunks(cs, nresultchunks, batchsize=batchsize, merge=(x,y)->Dagger.merge_sorted(ord, x, y))
    if nresultchunks == 1
        return dag1[1]
    else
        foreach(Dagger.persist!, dag1)
        return delayed((x...)->[x...], meta=true)(dag1...)
    end
end


text:
using Random

export JTable

mutable struct JTable
    data
end

JTable(x::AbstractArray) = JTable(DataFrame(x))

function Base.open(tbl::JTable)
    csv = joinpath(homedir(), ".cache", randstring() * ".csv")
    tbl.data.to_csv(csv, encoding = "gbk")
    finalizer(x -> rm(csv), tbl)
    Atom.msg("openFile", csv)
end


text:
"""
    focal_heuristic(env::MAPFEnvironment, solution::Vector{PlanResult})

Computes the (potentially inadmissible) focal heuristic value for the high-level node
in ECBS; see the original ECBS paper for details - https://www.aaai.org/ocs/index.php/SOCS/SOCS14/paper/view/8911
"""
function focal_heuristic end

# TODO: focal heuristics called on current solution?
@with_kw mutable struct ECBSHighLevelNode{S <: MAPFState, A <: MAPFAction, C <: Number, CNR <: MAPFConstraints}
    solution::Vector{PlanResult{S,A,C}} = Vector{PlanResult{S,A,C}}(undef, 0)
    constraints::Vector{CNR}            = Vector{CNR}(undef,0)
    cost::C                             = zero(C)
    lb::C                               = zero(C)
    focal_heuristic_value::C                   = zero(C)
    id::Int64                           = 0
end

Base.isless(hln1::ECBSHighLevelNode, hln2::ECBSHighLevelNode) = hln1.cost < hln2.cost

struct CompareFocalHeuristic
end

DataStructures.compare(comp::CompareFocalHeuristic, n1::ECBSHighLevelNode, n2::ECBSHighLevelNode) = (n1.focal_heuristic_value, n1.cost) < (n2.focal_heuristic_value, n2.cost)

@with_kw mutable struct ECBSSolver{S <: MAPFState, A <: MAPFAction, C <: Number, HC <: HighLevelCost, F <: MAPFConflict, CNR <: MAPFConstraints, E <: MAPFEnvironment}
    env::E
    hlcost::HC                                                                      = HC()
    weight::Float64                                                                 = 1.0
    heap::MutableBinaryMinHeap{ECBSHighLevelNode{S,A,C,CNR}}                            = MutableBinaryMinHeap{ECBSHighLevelNode{S,A,C,CNR}}()
    hmap::Dict{Int,Int}                                                             = Dict{Int,Int}()
    focal_heap::MutableBinaryHeap{ECBSHighLevelNode{S,A,C,CNR},CompareFocalHeuristic}   = MutableBinaryHeap{ECBSHighLevelNode{S,A,C,CNR},CompareFocalHeuristic}()
    focal_hmap::Dict{Int,Int}                                                       = Dict{Int,Int}()
    num_global_conflicts::Int64                                                     = 0
end

"""
    search!(solver::ECBSSolver{S,A,C,HC,F,CNR,E}, initial_states::Vector{S}) where {S <: MAPFState, A <: MAPFAction, C <: Number, HC <: HighLevelCost,
                                                                               F <: MAPFConflict, CNR <: MAPFConstraints, E <: MAPFEnvironment}

Calls the ECBS Solver on the given problem.
"""
function search!(solver::ECBSSolver{S,A,C,HC,F,CNR,E}, initial_states::Vector{S}) where {S <: MAPFState, A <: MAPFAction, C <: Number,  HC <: HighLevelCost,
                                                                                F <: MAPFConflict, CNR <: MAPFConstraints, E <: MAPFEnvironment}

    num_agents = length(initial_states)

    start = ECBSHighLevelNode{S,A,C,CNR}(solution = Vector{PlanResult{S,A,C}}(undef, num_agents),
                              constraints = Vector{CNR}(undef, num_agents))

    for idx = 1:num_agents

        start.constraints[idx] = get_empty_constraint(CNR)

        # Get constraints from low level context
        set_low_level_context!(solver.env, idx, start.constraints[idx])

        # Calls get_plan_result_from_astar within
        new_solution = low_level_search!(solver, idx, initial_states[idx], start.constraints[idx],
                                         Vector{PlanResult{S,A,C}}(undef, 0))

        # Return empty solution if cannot find
        if isempty(new_solution)
          return Vector{PlanResult{S,A,C}}(undef, 0)
        end

        start.solution[idx] = new_solution
        start.lb += start.solution[idx].fmin
    end
    start.cost = compute_cost(solver.hlcost, start.solution)
    start.focal_heuristic_value = focal_heuristic(solver.env, start.solution)

    # Push on to heaps
    solver.hmap[start.id] = push!(solver.heap, start)
    solver.focal_hmap[start.id] = push!(solver.focal_heap, start)

    best_cost = start.cost

    id = 1

    while ~(isempty(solver.heap))

        old_best_cost = best_cost
        best_cost = top(solver.heap).cost

        if best_cost > old_best_cost

            for node in sort(solver.heap.nodes, by = x->x.value.cost)
                val = node.value.cost

                if val > solver.weight * old_best_cost && val <= solver.weight * best_cost &&
                    ~(haskey(solver.focal_hmap, node.value.id))
                    solver.focal_hmap[node.value.id] = push!(solver.focal_heap, node.value)
                end

                if val > solver.weight * best_cost
                    break
                end
            end
        end

        # TODO : Check consistency?

        # Pop from focal list and do rest
        # @show solver.focal_heap
        focal_entry, focal_handle = top_with_handle(solver.focal_heap)
        heap_handle = solver.hmap[focal_entry.id]

        # @show focal_entry.id

        pop!(solver.focal_heap)
        delete!(solver.focal_hmap, focal_entry.id)
        delete!(solver.heap, heap_handle)
        delete!(solver.hmap, focal_entry.id)

        # Handle conflict
        conflict = get_first_conflict(solver.env, focal_entry.solution)

        # If no conflict, we are done
        if conflict == nothing
            @info "SOLVED! Cost: ",focal_entry.cost
            # IMP - One or more solutions might be empty - need to check at higher level
            return focal_entry.solution
        end

        solver.num_global_conflicts += 1

        # VECTOR OF CONSTRAINTS
        constraints = create_constraints_from_conflict(solver.env, conflict)

        for constraint in constraints
            for (i, c) in constraint

                # @show (i, c)

                new_node = deepcopy(focal_entry)

                add_constraint!(new_node.constraints[i], c)

                new_node.cost = deaccumulate_cost(solver.hlcost, new_node.cost, new_node.solution[i].cost)
                new_node.lb -= new_node.solution[i].fmin

                set_low_level_context!(solver.env, i, new_node.constraints[i])
                new_solution = low_level_search!(solver, i, initial_states[i], new_node.constraints[i], new_node.solution)

                if ~(isempty(new_solution))

                    # Can enter a new node
                    new_node.id = id
                    new_node.solution[i] = new_solution
                    new_node.cost = accumulate_cost(solver.hlcost, new_node.cost, new_solution.cost)
                    new_node.lb += new_solution.fmin
                    new_node.focal_heuristic_value = focal_heuristic(solver.env, new_node.solution)

                    solver.hmap[id] = push!(solver.heap, new_node)

                    if new_node.cost <= solver.weight * best_cost
                        solver.focal_hmap[id] = push!(solver.focal_heap, new_node)
                    end
                    id += 1
                end
            end
        end
    end

    # Return an empty solution
    return Vector{PlanResult{S,A,C}}(undef, 0)

end


text:
function orderTransitions!( mpSim::MPsim )

    transitions = vcat( getindex.( Ref( mpSim.transitionsByName ),
        collect( keys( mpSim.transitionsByName ) ) )... )
    reflextrans = filter( trans -> trans.sourceNode == trans.targetNode,
        transitions )
    filter!( trans -> trans ∉ reflextrans, transitions )
    recruitment = vcat( getindex.( Ref( mpSim.recruitmentByName ),
        collect( keys( mpSim.recruitmentByName ) ) )... )

    transGraph = MetaDiGraph( length( transitions ) + length( recruitment ) )

    # Don't take any action if there are no transitions.
    nv( transGraph ) == 0 && return
    jj = 1

    for trans in transitions
        set_prop!( transGraph, jj, :name, trans.name )
        set_prop!( transGraph, jj, :source, trans.sourceNode )
        set_prop!( transGraph, jj, :target,
            trans.isOutTransition ? "OUT" : trans.targetNode )
        set_prop!( transGraph, jj, :trans, trans )
        jj += 1
    end  # for trans in transitions

    for rec in recruitment
        set_prop!( transGraph, jj, :name, rec.name )
        set_prop!( transGraph, jj, :source, "IN" )
        set_prop!( transGraph, jj, :target, rec.targetNode )
        set_prop!( transGraph, jj, :trans, rec )
        jj += 1
    end  # for rec in recruitment

    transitions = vcat( transitions, recruitment )

    # Out transitions have priority over in transitions.
    for node in keys( mpSim.baseNodeList )
        outTransitions = findall( trans -> ( trans isa Transition ) &&
            ( trans.sourceNode == node ), transitions )
        inTransitions = findall( trans -> trans.targetNode == node,
            transitions )

        for inTrans in inTransitions, outTrans in outTransitions
            add_edge!( transGraph, outTrans, inTrans )
        end  # for inTrans in inTransitions
    end  # for node in mpSim.baseNodeList

    # Prioritise transitions by transition type.
    names = get_prop.( Ref( transGraph ), vertices( transGraph ), :name )
    uniqueNames = unique( names )
    sources = get_prop.( Ref( transGraph ), vertices( transGraph ), :source )
    uniqueSources = unique( sources )
    targets = get_prop.( Ref( transGraph ), vertices( transGraph ), :target )
    uniqueTargets = unique( targets )
    transOrder = mpSim.transitionTypeOrder
    nodeOrder = mpSim.baseNodeOrder

    for source in uniqueSources
        transFromSource = findall( sources .== source )
        length(transFromSource) < 2 && continue

        for trans1 in transFromSource, trans2 in transFromSource
            if haskey( transOrder, names[trans1] ) &&
                haskey( transOrder, names[trans2] ) &&
                ( transOrder[names[trans1]] <
                    transOrder[names[trans2]] )
                add_edge!( transGraph, trans1, trans2 )
            end  # if haskey( transOrder, names[trans1] ) && ...
        end  # for trans1 in transFromSource, trans2 in transFromSource

        transNames = names[transFromSource]

        for name in uniqueNames
            transNameSource = transFromSource[transNames .== name]
            length(transNameSource) < 2 && continue

            for trans1 in transNameSource, trans2 in transNameSource
                if haskey( nodeOrder, targets[trans1] ) &&
                    haskey( nodeOrder, targets[trans2] ) &&
                    ( nodeOrder[targets[trans1]] <
                        nodeOrder[targets[trans2]] )
                    add_edge!( transGraph, trans1, trans2 )
                end  # if haskey( nodeOrder, targets[trans1] ) && ...
            end  # for trans1 in transNameSource, trans2 in transNameSource
        end  # for name in uniqueNames
    end  # for source in uniqueSources

    for target in uniqueTargets
        transToTarget = findall( targets .== target )
        length(transToTarget) < 2 && continue

        for trans1 in transToTarget, trans2 in transToTarget
            if haskey( transOrder, names[trans1] ) &&
                haskey( transOrder, names[trans2] ) &&
                ( transOrder[names[trans1]] <
                    transOrder[names[trans2]] )
                add_edge!( transGraph, trans1, trans2 )
            end  # if haskey( transOrder, names[trans1] ) && ...
        end  # for trans1 in transToTarget, trans2 in transToTarget

        transNames = names[transToTarget]

        for name in uniqueNames
            transNameTarget = transToTarget[transNames .== name]

            length( transNameTarget ) < 2 && continue

            for trans1 in transNameTarget, trans2 in transNameTarget
                if haskey( nodeOrder, sources[trans1] ) &&
                    haskey( nodeOrder, sources[trans2] ) &&
                    ( nodeOrder[sources[trans1]] <
                        nodeOrder[sources[trans2]] )
                    add_edge!( transGraph, trans1, trans2 )
                end  # if haskey( nodeOrder, sources[trans1] ) && ...
            end  # for trans1 in transNameTarget, trans2 in transNameTarget
        end  # for name in uniqueNames
    end  # for target in uniqueTargets

    # println( "Is the graph connected? ", is_connected( transGraph ) )
    # println( "Is the graph cyclic? ", is_cyclic( transGraph ) )

    subgraphNodes = connected_components( transGraph )
    graphRanks = zeros( vertices( transGraph ) )

    for nodeList in subgraphNodes
        subgraph = induced_subgraph( transGraph, nodeList )[1]
        generateGraphNodeRanks!( subgraph )
        graphRanks[nodeList] = get_prop.( Ref( subgraph ),
            vertices( subgraph ), :rank )
    end  # for nodeList in subgraphNodes

    setfield!.( transitions, :priority, -Int.( graphRanks ) )
    mpSim.nPriorities = maximum( graphRanks )

end  # orderTransitions!( mpSim )

@resumable function transitionProcess( sim::Simulation, transition::Transition,
    mpSim::MPsim )

    processTime = Dates.Millisecond( 0 )
    tStart = now()

    # Initialize the schedule.
    priority = transition.priority
    priorityShift = ( transition.hasPriority ? 0 : 1 ) * mpSim.nPriorities
    maxAttempts = transition.maxAttempts == 0 ?
        length( transition.probabilityList ) : transition.maxAttempts
    nAttempts = Dict{String,Int}()

    timeToWait = transition.offset

    # If an initial population snapshot is uploaded, and it contains zero-time
    #   events, don't execute zero-time transition.
    if ( timeToWait == 0 ) && !mpSim.isVirgin
        timeToWait += transition.freq
    end  # if ( timeToWait == 0 ) && ...

    while now( sim ) + timeToWait <= mpSim.simLength
        processTime += now() - tStart
        @yield( timeout( sim, timeToWait, priority = priority -
            ( transition.minFlux == 0 ? priorityShift : 0 ) ) )
        tStart = now()
        timeToWait = transition.freq

        # Identify all persons who're in the start state long enough and are not
        #   already going to transition to another state.
        # eligibleIDs = getEligibleIDs( transition, nAttempts, maxAttempts,
        #     mpSim )
        # checkedIDs = checkExtraConditions( transition, eligibleIDs, mpSim )
        checkedIDs = getEligibleIDs( transition, nAttempts, maxAttempts, mpSim )
        updateAttemptsAndIDs!( nAttempts, maxAttempts, checkedIDs )

        # Assign a modified probability to each person.
        transitionLevels = determineTransitionLevels( transition, checkedIDs,
            nAttempts )

        # If the transition has a minimum flux, transition those IDs first, and
        #   set the priority of the rest of the transition to the proper level.
        if transition.minFlux > 0
            mandatoryIDs = determineMandatoryIDs( transition, checkedIDs,
                transitionLevels )
            performTransitions( mpSim, transition, mandatoryIDs )
            @yield( timeout( sim, 0, priority = priority - priorityShift ) )
            sourceNode = mpSim.baseNodeList[transition.sourceNode]
            filter!( id -> haskey( sourceNode.inNodeSince, id ), checkedIDs )
        end  # if transition.minFlux > 0

        transitionIDs = determineAdditionalIDs( transition, checkedIDs,
            transitionLevels, mpSim )
        performTransitions( mpSim, transition, transitionIDs )
    end  # while now( sim ) + timeToWait <= mpSim.simLength

    processTime += now() - tStart

    if mpSim.showInfo
        println( "Transition process for '", transition.name, "' took ",
            processTime.value / 1000, " seconds." )
    end  # if mpSim.showOutput

end  # transitionProcess( sim, transition, mpSim )

function getEligibleIDs( transition::Transition, nAttempts::Dict{String,Int},
    maxAttempts::Int, mpSim::MPsim )
    sourceNode = mpSim.baseNodeList[transition.sourceNode]
    eligibleIDs = collect( keys( sourceNode.inNodeSince ) )

    # Check if attempts have been exhausted (important if transition doesn't
    #   fire after exhaustion of attempts).
    # The first condition means 0 attempts so far.
    maxAttempts > 0 && filter!( id -> !haskey( nAttempts, id ) ||
        ( nAttempts[id] < maxAttempts ), eligibleIDs )
    isempty(eligibleIDs) && return eligibleIDs

    # Check extra conditions.
    queryCmd = string( "SELECT *, ageAtRecruitment + ", now(mpSim),
        " - timeEntered age, ", now(mpSim), " - timeEntered tenure, ",
        now(mpSim), " - inNodeSince `time in node` FROM `",
        mpSim.persDBname, "` WHERE ",
        "\n    `", mpSim.idKey, "` IN ('", join( eligibleIDs, "', '" ), "')",
        "\n        AND `time in node` > 0",
        join(string.( "\n        AND ",
            conditionToSQLite.( transition.extraConditions, Ref(mpSim) ) )) )
    Vector{String}( DataFrame(DBInterface.execute( mpSim.simDB,
        queryCmd ))[:, Symbol(mpSim.idKey)] )
end  # getEligibleIDs( transition, nAttempts, maxAttempts, mpSim )
#=
function getEligibleIDs( transition::Transition, nAttempts::Dict{String,Int},
    maxAttempts::Int, mpSim::MPsim )::Vector{String}

    sourceNode = mpSim.baseNodeList[transition.sourceNode]
    eligibleIDs = collect( keys( sourceNode.inNodeSince ) )

    # Check if attempts have been exhausted (important if transition doesn't
    #   fire after exhaustion of attempts).
    if maxAttempts > 0
        # The first condition means 0 attempts so far.
        filter!( id -> !haskey( nAttempts, id ) ||
            ( nAttempts[id] < maxAttempts ), eligibleIDs )
    end  # if maxAttempts > 0

    # filter!( id -> now( mpSim ) > sourceNode.inNodeSince[id], eligibleIDs )  # Needed???

    for cond in filter( cond -> cond.attribute == "time in node",
        transition.extraConditions )
        filter!( id -> cond.operator( now( mpSim ) -
            sourceNode.inNodeSince[id], cond.value ), eligibleIDs )
    end  # for cond in filter( ...

    return eligibleIDs

end  # getEligibleIDs( transition, mpSim )

function checkExtraConditions( transition::Transition,
    eligibleIDs::Vector{String}, mpSim::MPsim )

    if isempty( eligibleIDs ) || isempty( transition.extraConditions )
        return eligibleIDs
    end  # if isempty( eligibleIDs ) || ...

    # Check extra conditions.
    queryCmd = string( "SELECT *, ageAtRecruitment + ", now( mpSim ),
        " - timeEntered age, ", now( mpSim ), " - timeEntered tenure FROM `",
        mpSim.persDBname, "` WHERE ",
        "\n    `", mpSim.idKey, "` IN ('", join( eligibleIDs, "', '" ), "')" )
    eligibleIDsState = DataFrame( DBInterface.execute( mpSim.simDB, queryCmd ) )
    checkedIDs = Vector{String}( eligibleIDsState[:, Symbol( mpSim.idKey )] )
    isIDokay = trues( length( checkedIDs ) )

    for cond in filter( cond -> cond.attribute != "time in node",
        transition.extraConditions )
        isIDokay = isIDokay .& map( ii -> cond.operator(
            eligibleIDsState[:, Symbol( cond.attribute )][ii], cond.value ),
            eachindex( isIDokay ) )
    end  # for cond in filter( ... )

    return checkedIDs[isIDokay]

end  # checkExtraConditions( transition, eligibleIDs, mpSim )
=#

function updateAttemptsAndIDs!( nAttempts::Dict{String,Int}, maxAttempts::Int,
    eligibleIDs::Vector{String} )

    # Increase number of attempts for all eligible IDs. 0 attempts means
    #   the personnel member has already executed the transition
    #   succesfully while stil being eligible for it.
    for id in eligibleIDs
        if haskey( nAttempts, id ) && ( ( maxAttempts == -1 ) ||
            ( nAttempts[id] <= maxAttempts ) ) && ( nAttempts[id] > 0 )
            nAttempts[id] += 1
        else
            nAttempts[id] = 1
        end  # if haskey( nAttempts, id ) && ...
    end  # for id in eligibleIDs

end  # updateAttemptsAndIDs!( nAttempts, eligibleIDs )

function determineTransitionLevels( transition::Transition,
    eligibleIDs::Vector{String}, nAttempts::Dict{String,Int} )

    probDict = Dict{String,Float64}()

    if isempty( eligibleIDs )
        return probDict
    end  # if isempty( eligibleIDs )

    # The probabilities of undergoing the transition for each person.
    tresholds =  transition.probabilityList[
        min.( length( transition.probabilityList ),
        getindex.( Ref( nAttempts ), eligibleIDs ) )]
    modProbs = rand( transition.probRNG, length( eligibleIDs ) ) ./ tresholds
    setindex!.( Ref( probDict ), modProbs, eligibleIDs )
    return probDict

end  # determineTransitionLevels( trans, eligibleIDs, nAttempts )

function determineMandatoryIDs( transition::Transition,
    eligibleIDs::Vector{String}, transitionLevels::Dict{String,Float64} )

    if length( eligibleIDs ) <= transition.minFlux
        return eligibleIDs
    end  # if length( eligibleIDs ) <= transition.minFlux

    # Get the min flux IDs which have the smallest value for their transLevel.
    modProbs = getindex.( Ref( transitionLevels ), eligibleIDs )
    return eligibleIDs[sortperm( modProbs )[1:(transition.minFlux)]]

end  # determineMandatoryIDs( transition, eligibleIDs, transLevels )

function determineAdditionalIDs( transition::Transition,
    eligibleIDs::Vector{String}, transitionLevels::Dict{String,Float64},
    mpSim::MPsim )

    if transition.maxFlux == transition.minFlux
        return Vector{String}()
    end  # if transition.maxFlux == transition.minFlux

    if isempty( eligibleIDs )
        return eligibleIDs
    end  # if isempty( eligibleIDs )

    # Determine the max number of personnel members left to transfer.
    maxToTransfer = transition.maxFlux == -1 ? length( eligibleIDs ) :
        transition.maxFlux - transition.minFlux
    vacancies = -1

    targetNode = transition.isOutTransition ? nothing :
        mpSim.baseNodeList[transition.targetNode]

    if !transition.isOutTransition && !transition.hasPriority &&
        ( transition.sourceNode != transition.targetNode ) &&
        ( targetNode.target >= 0 )
        vacancies = max( 0, targetNode.target -
            length( targetNode.inNodeSince ) )
    end  # if !( transition.isOutTransition ) && ...

    if vacancies >= 0
        maxToTransfer = min( maxToTransfer, vacancies )
    end  # if vacancies >= 0

    if maxToTransfer == 0
        return Vector{String}()
    end  # if maxToTransfer == 0

    modProbs = getindex.( Ref( transitionLevels ), eligibleIDs )
    hasPassed = modProbs .< 1

    # If the number of people who have passed the test is smaller than the
    #   number of people who can transfer, return the list of people who passed.
    if count( hasPassed ) <= maxToTransfer
        return eligibleIDs[hasPassed]
    end  # if sum( hasPassed ) <= maxToTransfer

    # Otherwise, only select people with highest scores.
    return eligibleIDs[sortperm( modProbs )[1:maxToTransfer]]

end  # determineAdditionalIDs( transition, eligibleIDs, transitionLevels,
     #   mpSim )

function performTransitions( mpSim::MPsim, transition::Transition,
    idList::Vector{String} )

    if isempty( idList )
        return
    end  # if isempty( idList )

    if transition.isOutTransition
        removePersons( idList, transition.sourceNode, transition.name, mpSim )
    else
        # Execute the transition for the selected personnel members and
        #   perform necessary updates.
        executeTransitions( transition, idList, mpSim )

        if transition.sourceNode !== transition.targetNode
            updateNodes( transition, idList, mpSim )
        end  # if transition.sourceNode !== transition.targetNode
    end  # if transition.isOutTrans

end  # performTransitions( mpSim, transition, idList )

function executeTransitions( transition::Transition, idList::Vector{String},
    mpSim::MPsim )

    # Get current state of the personnel members undergoing the
    #   transition.
    queryCmd = string( "SELECT * FROM `", mpSim.persDBname, "` WHERE",
        "\n    `", mpSim.idKey, "` IN ('", join( idList, "', '" ), "')" )
    idListState = DataFrame( DBInterface.execute( mpSim.simDB, queryCmd ) )
    tmpIDs = idListState[:, Symbol( mpSim.idKey )]

    # Make changes in personnel database.
    persChangesCmd = [string( "currentNode = '", transition.targetNode, "'" ),
        string( "inNodeSince = ", now(mpSim) )]
    changedAttrVals = Dict{String,String}()

    # Target node attributes.
    if !transition.isOutTransition &&
        ( transition.sourceNode !== transition.targetNode )
        targetNode = mpSim.baseNodeList[transition.targetNode]
        changedAttrVals = deepcopy( targetNode.requirements )
    end  # if !transition.isOutTransition && ...

    # Extra changes because of the transition.
    changedAttrs = collect( keys( transition.extraChanges ) )
    setindex!.( Ref( changedAttrVals ),
        getindex.( Ref( transition.extraChanges ), changedAttrs ),
        changedAttrs )

    # Change attributes in personnel database.
    changedAttrs = collect( keys( changedAttrVals ) )
    append!( persChangesCmd, string.( "`", changedAttrs, "` = '",
        getindex.( Ref( changedAttrVals ), changedAttrs ), "'" ) )
    persChangesCmd = string( "UPDATE `", mpSim.persDBname, "`",
        "\n    SET ", join( persChangesCmd, ", " ),
        "\n    WHERE `", mpSim.idKey, "` IN ('", join( tmpIDs, "', '" ), "')" )
    DBInterface.execute( mpSim.simDB, persChangesCmd )

    # Record attribute changes in history database.
    histChangesCmd = Vector{String}()

    for attr in changedAttrs
        oldAttrValues = idListState[:, Symbol( attr )]
        changedIDs = tmpIDs[oldAttrValues .!= changedAttrVals[attr]]

        if !isempty( changedIDs )
            push!( histChangesCmd, join( string.( "\n    ('", changedIDs,
                "', '", attr, "', ", now( mpSim ), ", '",
                changedAttrVals[attr], "')" ), ", " ) )
        end  # if !isempty( changedIDs )
    end  # for attr in changedAttrs

    if !isempty( histChangesCmd )
        histChangesCmd = string( "INSERT INTO `", mpSim.histDBname, "` (`",
            mpSim.idKey, "`, attribute, timeIndex, value) VALUES",
            join( histChangesCmd, "," ) )
        DBInterface.execute( mpSim.simDB, histChangesCmd )
    end  # if !isempty( histChangesCmd )

end  # executeTransitions( trans, idList, mpSim )

function updateNodes( transition::Transition, transIDs::Vector{String},
    mpSim::MPsim )

    # Update the inNodeSince fields of the source and target nodes.
    sourceNode = mpSim.baseNodeList[transition.sourceNode]
    tName = transition.isOutTransition ? "NULL" :
        string( "'", transition.targetNode, "'" )
    delete!.( Ref( sourceNode.inNodeSince ), transIDs )

    if !transition.isOutTransition
        targetNode = mpSim.baseNodeList[transition.targetNode]
        setindex!.( Ref( targetNode.inNodeSince ), now( mpSim ), transIDs )
    end  # if !transition.isOutTransition

    # Update the database of node transitions.
    transCmd = string.( "\n    ('", transIDs, "', ", now( mpSim ), ", '",
        transition.name, "', '", transition.sourceNode, "', ", tName, ")" )
    transCmd = string( "INSERT INTO `", mpSim.transDBname, "` (`", mpSim.idKey,
        "`, timeIndex, transition, sourceNode, targetNode) VALUES",
        join( transCmd, "," ) )
        DBInterface.execute( mpSim.simDB, transCmd )

    # Update the expected time of attrition.
    updateTimeOfAttrition( transIDs, transition.targetNode, mpSim )

end  # function updateNodes( transition, transIDs, mpSim )

function updateTimeOfAttrition( transIDs::Vector{String}, targetName::String,
    mpSim::MPsim  )

    # Update the expected times to attrition.
    targetNode = mpSim.baseNodeList[targetName]
    attrition = mpSim.attritionSchemes[targetNode.attrition]
    timeOfAttrition = now( mpSim ) .+ generateTimeToAttrition( attrition,
        length( transIDs ) )
    timeOfAttrition = map( toa -> toa == +Inf ? "NULL" : toa, timeOfAttrition )

    updateCmds = string.( "UPDATE `", mpSim.persDBname, "`",
        "\n    SET expectedAttritionTime = ", timeOfAttrition,
        "\n    WHERE `", mpSim.idKey, "` IS '", transIDs, "'" )
    SQLite.execute.( Ref( mpSim.simDB ), updateCmds )

    # Check if any of the new attrition times occur before the next attrition
    #   process check.
    # newProcessCheckTime = ceil( now( mpSim ) / mpSim.attritionTimeSkip ) *
    #     mpSim.attritionTimeSkip

    # for ii in eachindex( idList )
    #     if newAttrTimeList[ii] <= newProcessCheckTime
    #         executeAttritionProcess( mpSim.sim, idList[ii],
    #             newAttrTimeList[ii], mpSim )
    #     end  # if newAttrTimeList[ii] <= newProcessCheckTime
    # end  # for ii in eachindex( idList )

end  # updateTimeOfAttrition( transIDs, targetName, mpSim )

include(joinpath( simPrivPath, "graphops.jl" ))


text:
include("generar-vecino.jl")

"""
PrimeroMejor(numero_maximo_vecinos_sin_mejora::Int,
            evaluaciones_maximas_funcion_objetivo::Int,
            numero_atributos::Int,
            umbral_atributo::Real,
            F)::Vector{<:Real}

Devuelve mejor vector de peso encontrando de acorde a función
de evaluación `F`
"""
function PrimeroMejor(numero_maximo_vecinos_sin_mejora::Int,
    evaluaciones_maximas_funcion_objetivo::Int,
    numero_atributos::Int,
    umbral_atributo::Real,
    F,
    w = 0)::Vector{<:Real}
    if (w == 0)
        w = rand(numero_atributos) # vector de tamaño atributo uniformement inicializado
    end
    # Ponemos a cero atributos menores que umbral atributo
    w = map(x-> (x<umbral_atributo) ? 0 : x , w)

    iteraciones = 0
    evaluaciones = 0
    F_w = F(w)
    indice::Int = 1
    while iteraciones < numero_maximo_vecinos_sin_mejora &&
        evaluaciones < evaluaciones_maximas_funcion_objetivo
        w_n = GenNeighbourhoodIesimo(w,0, 0.3, indice)
        F_w_vecino = F(w_n)

        if(F_w < F_w_vecino)
            w = w_n
            F_w = F_w_vecino
            iteraciones = 0
        else
            # Si no hay mejora se cambia la exploración del atributo
            indice = indice % (numero_atributos)+1
        end
        iteraciones += 1
        evaluaciones += 1
    end
    return w
end


text:
function perform_step!(integ::ODEFilterIntegrator)
    @unpack t, dt = integ
    @unpack E0, Precond, InvPrecond = integ.constants
    @unpack x_pred, u_pred, x_filt, u_filt, err_tmp = integ.cache

    P = Precond(dt)
    PI = InvPrecond(dt)
    integ.cache.x = P * integ.cache.x

    t = t + dt
    integ.t_new = t

    x_pred = predict!(integ)
    mul!(u_pred, E0, PI*x_pred.μ)

    measure!(integ, x_pred, t)

    if isdynamic(integ.sigma_estimator)
        # @info "Before dynamic sigma:" x_pred.Σ
        σ_sq = dynamic_sigma_estimation(integ.sigma_estimator, integ)

        # Adjust prediction and measurement accordingly
        x_pred.Σ .+= (σ_sq .- 1) .* integ.cache.Qh
        integ.cache.measurement.Σ .+= integ.cache.H * ((σ_sq .- 1) .* integ.cache.Qh) * integ.cache.H'

        integ.cache.σ_sq = σ_sq
    end

    x_filt = update!(integ, x_pred)
    mul!(u_filt, E0, PI*x_filt.μ)

    if isstatic(integ.sigma_estimator)
        # E.g. estimate the /current/ MLE sigma; Needed for error estimation
        σ_sq = static_sigma_estimation(integ.sigma_estimator, integ)
        integ.cache.σ_sq = σ_sq
    end

    err_est_unscaled = estimate_errors(integ)
    # Scale the error with old u-values and tolerances
    DiffEqBase.calculate_residuals!(
        err_tmp,
        dt * err_est_unscaled, integ.u, u_filt, integ.opts.abstol, integ.opts.reltol, integ.opts.internalnorm, t)
    err_est_combined = integ.opts.internalnorm(err_tmp, t)  # Norm over the dimensions
    integ.EEst = err_est_combined

    integ.cache.x = PI * integ.cache.x
    integ.cache.x_pred = PI * integ.cache.x_pred
    integ.cache.x_filt = PI * integ.cache.x_filt
end

function predict!(integ::ODEFilterIntegrator)

    @unpack dt = integ
    @unpack A!, Q!, InvPrecond = integ.constants
    @unpack x, Ah, Qh, x_pred = integ.cache
    PI = InvPrecond(dt)

    A!(Ah, dt)
    Q!(Qh, dt)

    mul!(x_pred.μ, Ah, x.μ)
    x_pred.Σ .= Symmetric(Ah * x.Σ * Ah' .+ Qh)

    return x_pred
end

function measure_h!(integ::ODEFilterIntegrator, x_pred, t)

    @unpack p, f, dt = integ
    @unpack E0, h!, InvPrecond = integ.constants
    @unpack du, h, u_pred = integ.cache
    PI = InvPrecond(dt)

    IIP = isinplace(integ)
    if IIP
        f(du, u_pred, p, t)
    else
        du .= f(u_pred, p, t)
    end
    integ.destats.nf += 1

    h!(h, du, PI*x_pred.μ)
end

function measure_H!(integ::ODEFilterIntegrator, x_pred, t)

    @unpack p, f, dt = integ
    @unpack jac, H!, InvPrecond = integ.constants
    @unpack u_pred, ddu, H = integ.cache
    PI = InvPrecond(dt)

    if !isnothing(jac)
        if isinplace(integ)
            jac(ddu, u_pred, p, t)
        else
            ddu .= jac(u_pred, p, t)
        end
        integ.destats.njacs += 1
    end
    H!(H, ddu)
    H .= H * PI
end

function measure!(integ, x_pred, t)
    measure_h!(integ, x_pred, t)
    measure_H!(integ, x_pred, t)

    @unpack dt = integ
    @unpack R, q, d = integ.constants
    @unpack measurement, h, H = integ.cache

    v, S = measurement.μ, measurement.Σ
    v .= 0 .- h
    # R .= Diagonal(eps.(v) .^ 2)
    R .= Diagonal(eps.(v))
    # R .= Diagonal(dt^(2q) * ones(d))
    S .= Symmetric(H * x_pred.Σ * H' .+ R)

    return nothing
end

function update!(integ::ODEFilterIntegrator, prediction)

    @unpack dt = integ
    @unpack R, q, d, Precond, InvPrecond, E1 = integ.constants
    @unpack measurement, h, H, K, x_filt = integ.cache
    P, PI = Precond(dt), InvPrecond(dt)

    v, S = measurement.μ, measurement.Σ

    m_p, P_p = prediction.μ, prediction.Σ

    S_inv = inv(S)
    K .= P_p * H' * S_inv

    x_filt.μ .= m_p .+ K*v
    # x_filt.μ .= m_p .+ P_p * H' * (S\v)

    # Vanilla
    # KSK = K*S*K'
    # KSK = P_p * H' * (S \ (H * P_p'))
    # x_filt.Σ .= P_p .- KSK
    # approx_diff!(x_filt.Σ, P_p, KSK)

    # Joseph Form
    x_filt.Σ .= Symmetric(X_A_Xt(PDMat(P_p), (I-K*H)))
    if !iszero(R)
        x_filt.Σ .+= Symmetric(X_A_Xt(PDMat(R), K))
    end

    assert_nonnegative_diagonal(x_filt.Σ)
    # cholesky(Symmetric(x_filt.Σ))

    return x_filt
end

function approx_diff!(A, B, C)
    @assert size(A) == size(B) == size(C)
    @assert eltype(A) == eltype(B) == eltype(C)
    # If B_ij ≈ C_ij, then A_ij = 0
    # But, only do this if the value in A is actually negative
    @simd for i in 1:length(A)
        @inbounds if B[i] ≈ C[i]
            A[i] = 0
        else
            A[i] = B[i] - C[i]
        end
    end
end

function fix_negative_variances(g::Gaussian, abstol::Real, reltol::Real)
    for i in 1:length(g.μ)
        if (g.Σ[i,i] < 0) && (sqrt(-g.Σ[i,i]) / (abstol + reltol * abs(g.μ[i]))) .< eps(eltype(g.Σ))
            g.Σ[i,i] = eps(eltype(g.Σ))*(abstol+reltol*abs(g.μ[i]))
            g.Σ[i,i] = 0
        end
    end
end

function estimate_errors(integ)
    @unpack dt = integ
    @unpack InvPrecond = integ.constants
    @unpack σ_sq, Qh, H = integ.cache

    if σ_sq isa Real && isinf(σ_sq)
        return Inf
    end

    error_estimate = sqrt.(diag(H * σ_sq*Qh * H'))

    return error_estimate
end

max_stars_repo_name: nathanaelbosch/capos

max_issues_repo_name: nathanaelbosch/capos

max_forks_repo_name: nathanaelbosch/capos


text:
# define custom types required for the test
type TestCustomPersonType
	firstname::String
	lastname::String
end

type TestKnowsType
	someproperty::String
end

function savetest()
	# create a graph
	g = Graph()
	g.tracker = ChangeTracker()

	# create a vertex and add it to the graph
	fred = add!(g, Vertex("Person", @compat Dict{String,Any}("name"=>"Fred","age"=>45,"other"=>50)))

	# make some more vertices
	sally = add!(g, Vertex("Person", @compat Dict{String,Any}("name"=>"Sally","age"=>25)))

	johanobject = TestCustomPersonType("Johan", "Smith")
	johan = vertexforobject("Person", johanobject,  @compat Dict{String,Any}("name"=>"Johan","age"=>28))
	add!(g, johan)

	# create some edges
	e = Edge(fred,sally)
	e2 = Edge("knows", fred, sally, @compat Dict{String,Any}("since"=>2004))

	#including one with an associated object
	fredknowsjohanobject = TestKnowsType("SomeValue")
	e3 = edgeforobject("knows", fredknowsjohanobject, fred,johan)

	# add the edges to the graph
	add!(g, e)
	add!(g, e2)
	add!(g, e3)

	s = MongoStore(g)

	result = persist(s)

	@test result.insertverticesresult.issuccessful == true
	@test result.insertedgesresult.issuccessful == true
	@test result.issuccessful == true

	setpropertyvalues!(fred, @compat Dict{String,Any}("other"=>101,"new"=>"new value"))
	setpropertyvalues!(e2, @compat Dict{String,Any}("since"=>2003))

	result = persist(s)

	@test result.updateverticesresult.issuccessful == true
	@test result.updateedgesresult.issuccessful == true
	@test result.issuccessful == true

	remove!(g, e)
	remove!(g, johan)

	result = persist(s)

	@test result.removeverticesresult.issuccessful == true
	@test result.removeedgesresult.issuccessful == true
	@test result.issuccessful == true

end

savetest()


text:
atol = 1e-4
rtol = 1e-2
rtol_insol = 0.1

## Test zero insolation at night
# sunrise at equator
date = Dates.DateTime(2020, 1, 1, 6, 0, 0)
lon, lat = [FT(0.0), FT(0.0)]
sza, azi, d = instantaneous_zenith_angle(date, lon, lat, param_set)
F = insolation(sza, d, param_set)
@test F ≈ 0.0 atol=atol

S, mu = solar_flux_and_cos_sza(date, lon, lat, param_set)
@test S ≈ tot_solar_irrad(param_set) rtol=rtol_insol
@test mu ≈ 0.0 atol=atol

# polar night NH 1
date = Dates.DateTime(2020, 12, 20, 11, 0, 0)
lon, lat = [FT(0.0), FT(80.0)]
sza, azi, d = instantaneous_zenith_angle(date, lon, lat, param_set)
F = insolation(sza, d, param_set)
@test F ≈ 0.0 atol=atol

S, mu = solar_flux_and_cos_sza(date, lon, lat, param_set)
@test S ≈ tot_solar_irrad(param_set) rtol=rtol_insol
@test mu ≈ 0.0 atol=atol

# polar night NH 2
date = Dates.DateTime(2020, 12, 20, 23, 0, 0)
sza, azi, d = instantaneous_zenith_angle(date, lon, lat, param_set)
F = insolation(sza, d, param_set)
@test F ≈ 0.0 atol=atol

S, mu = solar_flux_and_cos_sza(date, lon, lat, param_set)
@test S ≈ tot_solar_irrad(param_set) rtol=rtol_insol
@test mu ≈ 0.0 atol=atol

## Test symmetry of insolation at equinox
nlats = 181
date = Dates.DateTime(2021, 3, 20, 9, 37, 0) # vernal equinox 2021
l_arr = FT.(collect(range(-90, stop = 90, length = nlats)))
F_arr = zeros(nlats)
for (i, lat) in enumerate(l_arr)
    θ, dist = daily_zenith_angle(date, lat, param_set)
    F_arr[i] = insolation(θ, dist, param_set)
end
F_NH = sort(F_arr[l_arr .>= 0])
F_SH = sort(F_arr[l_arr .<= 0])
@test F_NH ≈ F_SH rtol=rtol

## Test globally averaged insolation ≈ TSI
ndays, nlats = [365, 361]
d_arr = Array{Int}(round.(collect(range(0, stop = 365, length = ndays))))
l_arr = FT.(collect(range(-90, stop = 90, length = nlats)))
F_arr = zeros(ndays, nlats)

for (i, d) in enumerate(d_arr)
    for (j, lat) in enumerate(l_arr)
        datei = Dates.DateTime(2000,1,1) + Dates.Day(d)
        θ, dist = daily_zenith_angle(datei, lat, param_set)
        F_arr[i, j] = insolation(θ, dist, param_set)
    end
end

zonal_mean_insol = mean(F_arr, dims=1)
area_fac = abs.(cosd.(l_arr))
global_mean_insol = sum(zonal_mean_insol * area_fac) / sum(area_fac)
@test global_mean_insol ≈ tot_solar_irrad(param_set) / 4 rtol=rtol

## Test invariance of zonal-mean insolation under rotation of ϖ
ϖ0 = lon_perihelion_epoch(param_set)
CLIMAParameters.Planet.lon_perihelion_epoch(::EarthParameterSet) = ϖ0 + π

for (i, d) in enumerate(d_arr)
    for (j, lat) in enumerate(l_arr)
        datei = Dates.DateTime(2000,1,1) + Dates.Day(d)
        θ, dist = daily_zenith_angle(datei, lat, param_set)
        F_arr[i, j] = insolation(θ, dist, param_set)
    end
end

zonal_mean_insol_rotate = mean(F_arr, dims=1)
@test zonal_mean_insol_rotate ≈ zonal_mean_insol rtol=rtol

CLIMAParameters.Planet.lon_perihelion_epoch(::EarthParameterSet) = ϖ0


text:
using Formatting
using PrettyTables

# --------------- Instrument ---------------

function Base.show(io::IO, inst::Instrument)
    data_str = isnothing(inst.data) ? "" : "  data=<object>"
    print(io, "[Instrument] symbol=$(inst.symbol)$data_str")
end

# --------------- BidAsk ---------------

function Base.show(io::IO, ba::BidAsk)
    print(io, "[BidAsk] dt=$(ba.dt)  bid=$(ba.bid)  ask=$(ba.ask)")
end

# --------------- Position ---------------

function Base.show(io::IO, pos::Position)
    size_str = fmt(".2f", pos.size)
    if sign(pos.size) != -1
        size_str = " " * size_str
    end
    pnl_str = fmt("+.2f", pos.pnl)
    data_str = isnothing(pos.data) ? "nothing" : "<object>"
    print(io, "[Position] $(pos.inst.symbol) $(pos.dir) $size_str  "*
        "open=($(Dates.format(pos.open_dt, "yyyy-mm-dd HH:MM:SS")), $(fmt(".2f", pos.open_price)))  "*
        "last=($(Dates.format(pos.last_dt, "yyyy-mm-dd HH:MM:SS")), $(fmt(".2f", pos.last_price)))  "*
        "pnl=$pnl_str  stop_loss=$(fmt(".2f", pos.stop_loss))  take_profit=$(fmt(".2f", pos.take_profit))  "*
        "close_reason=$(pos.close_reason)  data=$data_str")
end

function print_positions(positions::Vector{Position};
    max_print=25, volume_digits=1, price_digits=2,
    data_renderer::Union{Function, Nothing}=nothing, kwargs...
)
    print_positions(stdout, positions::Vector{Position}; max_print, volume_digits, price_digits, data_renderer, kwargs...)
end

function print_positions(io::IO, positions::Vector{Position};
    max_print=25, volume_digits=1, price_digits=2,
    data_renderer::Union{Function, Nothing}=nothing
)
    n = length(positions)
    if n == 0
        return
    end

    df = dateformat"yyyy-mm-dd HH:MM:SS"

    columns = [
        "Symbol"; "Volume"; "Open time"; "Open price";
        "Last quote"; "Last price"; "P&L";
        "Stop loss"; "Take profit"; "Close reason"
    ]
    if !isnothing(data_renderer)
        push!(columns, "Data")
    end

    idx_volume = findfirst(columns .== "Volume")
    idx_open_time = findfirst(columns .== "Open time")
    idx_open_price = findfirst(columns .== "Open price")
    idx_last_price = findfirst(columns .=="Last price")
    idx_last_quote = findfirst(columns .== "Last quote")
    idx_pnl = findfirst(columns .== "P&L")
    idx_stop_loss = findfirst(columns .== "Stop loss")
    idx_take_profit = findfirst(columns .== "Take profit")
    idx_data = findfirst(columns .== "Data")
    idx_close_reason = findfirst(columns .== "Close reason")

    formatter = (v, i, j) -> begin
        o = v
        if j == idx_volume
            o = sprintf1("%.$(volume_digits)f", v)
        elseif j == idx_open_price || j == idx_last_price || j == idx_pnl || j == idx_stop_loss || j == idx_take_profit
            o = isnan(v) ? "—" : sprintf1("%.$(price_digits)f", v)
        elseif j == idx_open_time || j == idx_last_quote
            o = Dates.format(v, df)
        elseif j == idx_data
            # position "data" field renderer
            o = data_renderer(positions[j], v)
        elseif j == idx_close_reason
            o = v == Unspecified::CloseReason ? "—" : v
        end
        o
    end

    if !isnothing(idx_data)
        data = map(pos -> [pos.inst.symbol pos.size pos.open_dt pos.open_price pos.last_dt pos.last_price pos.pnl pos.stop_loss pos.take_profit pos.close_reason pos.data], positions)
    else
        data = map(pos -> [pos.inst.symbol pos.size pos.open_dt pos.open_price pos.last_dt pos.last_price pos.pnl pos.stop_loss pos.take_profit pos.close_reason], positions)
    end
    data = reduce(vcat, data)

    h_pos_green = Highlighter((data, i, j) -> j == idx_pnl && data[i, j] > 0, bold=true, foreground=:green)
    h_neg_red = Highlighter((data, i, j) -> j == idx_pnl && data[i, j] < 0, bold=true, foreground=:red)

    if !isnan(max_print) && size(data, 1) > max_print
        pretty_table(io, data[1:max_print, :], columns;
            highlighters = (h_pos_green, h_neg_red),
            formatters = formatter,
            compact_printing = false)
        println(io, " [...] $(n - max_print) more positions")
    else
        pretty_table(io, data, columns;
            highlighters = (h_pos_green, h_neg_red),
            formatters = formatter,
            compact_printing = false)
    end
end

# --------------- OpenOrder ---------------

function Base.show(io::IO, order::OpenOrder)
    print(io, "[OpenOrder] $(order.inst.symbol) $(order.size) $(order.dir)  stop_loss=$(fmt(".2f", order.stop_loss))  take_profit=$(fmt(".2f", order.take_profit))")
end

# --------------- CloseOrder ---------------

function Base.show(io::IO, order::CloseOrder)
    print(io, "[CloseOrder] $(order.pos)  $(order.close_reason)")
end

# --------------- CloseAllOrder ---------------

function Base.show(io::IO, order::CloseAllOrder)
    print(io, "[CloseAllOrder]")
end

# --------------- Account ---------------

# TODO: Digits
function Base.show(io::IO, acc::Account; volume_digits=1, price_digits=2, kwargs...)
    x, y = displaysize(io)
    get_color(val) = val >= 0 ? (val == 0 ? :black : :green) : :red
    title = " ACCOUNT SUMMARY "
    border_char = '━'
    eq_width = y - length(title)
    eqs = border_char^(floor(Int64, eq_width/2))
    println(io, "")
    println(io, eqs * title * eqs)
    println(io, " ", "Open positions:     $(length(acc.open_positions))")
    print_positions(io, acc.open_positions; kwargs...)
    println(io, "")
    println(io, " ", "Closed positions:   $(length(acc.closed_positions))")
    print_positions(io, acc.closed_positions; kwargs...)
    println(io, "")
    println(io, " ", "Initial balance:    $(fmt(".2f", acc.initial_balance))")
    print(io,   " ", "Balance:            $(fmt(".2f", acc.balance))")
    print(io, " (")
    printstyled(io, "$(fmt("+.2f", balance_ret(acc)*100))%"; color=get_color(balance_ret(acc)))
    print(io, ")\n")
    print(io, " ", "Equity:             $(fmt(".2f", acc.equity))")
    print(io, " (")
    printstyled(io, "$(fmt("+.2f", equity_ret(acc)*100))%"; color=get_color(equity_ret(acc)))
    print(io, ")\n")
    println(io, border_char^y)
    println(io, "")
end


text:
function rawscore(tree::CSGNode, points, normals, params)
    @unpack ϵ_d, α = params
    # λ should be unpacked?
    λ = log(size(points, 1))
    score = 0.

    for i in eachindex(points)
        v, n = valueandnormal(tree, points[i])
        d_i = v/ϵ_d
        ddot = clamp(dot(n, normals[i]), -1, 1)
        θ_i = acos(ddot)/α
        score += exp(-d_i^2) + exp(-θ_i^2)
    end
    return score - λ*treesize(tree)
end

function rankpopulation(population, points, normals, params)
    sum = 0.
    score = Array{Float64,1}(undef, size(population))
    normed = similar(score)

    @threads for i in eachindex(population)
    #for i in eachindex(population)
        score[i] = rawscore(population[i], points, normals, params)
    end
    for i in eachindex(score)
        normed[i] = 1/(1+score[i])
        sum += normed[i]
    end
    for i in eachindex(score)
        score[i] = normed[i]/sum
    end
    # this is the rescaled score

    p = sortperm(score)
    return population[p], score[p]
end

function mutate(creature, surfaces, params)
    @unpack μ, μ_0, maxdepth = params

    if rand() < μ
        # return mutated
        if rand() < μ_0
            # return new tree
            return randomtree(surfaces, maxdepth)
        else
            # mutate tree
            creat_copy = deepcopy(creature)
            podfs = PreOrderDFS(creat_copy)
            nofnodes = numberofnodes(podfs)
            isempty(1:nofnodes) && return creature
            randnode = rand(1:nofnodes)
            selected_node = selectfirstchildnode(podfs, randnode)
            nd = depth(selected_node)
            newdepth = nd > 0 ? nd : 1
            newnode = randomtree(surfaces, newdepth)
            n = 0
            for j in podfs
                isempty(j.children) && continue
                n += 1
                if n == randnode
                    j.children[1] = newnode
                    break
                end
            end
            if depth(creat_copy) > maxdepth
                return creature
            end
            return creat_copy
        end
    else
        # don't mutate tree
        return creature
    end
end

function geneticbuildtree(surfaces, points, normals, params)
    @unpack itermax, maxdepth, populationsize = params
    @unpack keepbestn = params
    notifit = itermax > 10 ? div(itermax,10) : 1
    snodes = [CSGNode(s, []) for s in surfaces]
    population = [randomtree(snodes, maxdepth) for _ in 1:populationsize]
    npopulation = similar(population)
    start_time = time_ns()
    for i in 1:itermax
        if i%notifit == 0
            @info "$i-th iteration"
        end
        population, _ = rankpopulation(population, points, normals, params)
        @debug "ranked population"
        # save the best
        npopulation[1:keepbestn] = population[1:keepbestn]
        n = keepbestn+1
        while true
            # select 2 creatures by tournament
            t1 = deepcopy(population[tournament(params)])
            t2 = deepcopy(population[tournament(params)])
            # crossover them
            newc = crossover([t1, t2], params)
            # mutate first
            cv = mutate(newc[1], snodes, params)
            npopulation[n] = cv
            n += 1
            # finish if we have enough
            n > populationsize && break
            # mutate second
            cv = mutate(newc[2], snodes, params)
            npopulation[n] = cv
            n += 1
            n > populationsize && break
        end
        # npopulation is populated with a new set of creature
        population = npopulation
    end
    fint = trunc((time_ns() - start_time)/1_000_000_000, digits=2)
    @info "Finished: $itermax iteration in $fint seconds."
    return population
end

function rawscorecached(tree, cpoints, cnormals, normals, params)
    @unpack ϵ_d, α = params
    # λ should be unpacked?
    λ = log(size(cpoints, 1))
    score = 0.
    for i in eachindex(cpoints)
        v, n = valueandnormal(tree, cpoints, cnormals, i)
        d_i = v/ϵ_d
        ddot = clamp(dot(n, normals[i]), -1, 1)
        θ_i = acos(ddot)/α
        score += exp(-d_i^2) + exp(-θ_i^2)
    end
    return score - λ*treesize(tree)
end

function rankcachedpopulation(population, cpoints, cnormals, normals, params)
    sum = 0.
    score = Array{Float64,1}(undef, size(population))
    normed = similar(score)

    @threads for i in eachindex(population)
    #for i in eachindex(population)
        score[i] = rawscorecached(population[i], cpoints, cnormals, normals, params)
    end
    for i in eachindex(score)
        normed[i] = 1/(1+score[i])
        sum += normed[i]
    end
    for i in eachindex(score)
        score[i] = normed[i]/sum
    end
    # this is the rescaled score

    p = sortperm(score)
    return population[p], score[p]
end

function cachedgeneticbuildtree(surfaces, points, normals, params)
    @unpack itermax, maxdepth, populationsize = params
    @unpack keepbestn = params
    notifit = itermax > 10 ? div(itermax,10) : 1
    cnodes, cvalues, cnormals = cachenodes(surfaces, points)
    population = [randomcachedtree(cnodes, maxdepth) for _ in 1:populationsize]
    npopulation = similar(population)
    start_time = time_ns()
    @info "Iteration in progress..."
    for i in 1:itermax
        if i%notifit == 0
            @info "$i-th iteration"
        end
        population, _ = rankcachedpopulation(population, cvalues, cnormals, normals, params)
        @debug "ranked population"
        # save the best
        npopulation[1:keepbestn] = population[1:keepbestn]
        n = keepbestn+1
        while true
            # select 2 creatures by tournament
            t1 = deepcopy(population[tournament(params)])
            t2 = deepcopy(population[tournament(params)])
            # crossover them
            newc = crossover([t1, t2], params)
            # mutate first
            cv = mutatecached(newc[1], cnodes, params)
            npopulation[n] = cv
            n += 1
            # finish if we have enough
            n > populationsize && break
            # mutate second
            cv = mutatecached(newc[2], cnodes, params)
            npopulation[n] = cv
            n += 1
            n > populationsize && break
        end
        # npopulation is populated with a new set of creature
        population = npopulation
    end
    fint = trunc((time_ns() - start_time)/1_000_000_000, digits=2)
    @info "Finished: $itermax iteration in $fint seconds."
    return population, cached2normaltree(population[1], surfaces)
end


text:
using DiffEqJump, DelaySSAToolkit
using Random, Distributions

rn = @reaction_network begin
    kon, Goff --> Gon
    koff, Gon --> Goff
    ρ, Gon --> Gon + N
end kon koff ρ
jumpsys = convert(JumpSystem, rn, combinatoric_ratelaws = false)

u0 = [1,0,0]
de_chan0 = [[]]
tf = 2000.
tspan = (0,tf)
p = [0.0282,0.609,2.11]
dprob = DiscreteProblem(u0, tspan, p)

delay_trigger_affect! = function (integrator, rng)
    τ=rand(LogNormal(1,sqrt(2)))+120
    append!(integrator.de_chan[1], τ)
end
delay_trigger = Dict(3=>delay_trigger_affect!)
delay_complete = Dict(1=>[3=>-1])
delay_interrupt = Dict()
delayjumpset = DelayJumpSet(delay_trigger,delay_complete,delay_interrupt)

djprob = DelayJumpProblem(jumpsys, dprob, DelayRejection(), delayjumpset, de_chan0, save_positions=(false,false))

ensprob = EnsembleProblem(djprob)
@time ens = solve(ensprob, SSAStepper(), EnsembleThreads(), trajectories=10^5)

using DifferentialEquations.EnsembleAnalysis, Plots; theme(:vibrant)
last_slice = componentwise_vectors_timepoint(ens, tf)
histogram(last_slice[3],bins=0:1:60,normalize=:pdf, label = "LogNormal(1,sqrt(2))",fillalpha = 0.6, linecolor = :orange,fmt=:svg,xlabel = "# of products", ylabel = "Probability")
# savefig("docs/src/assets/stochastic_delay1.svg")

delay_trigger_affect2! = function (integrator, rng)
    τ=rand(LogNormal(0,2))+120
    append!(integrator.de_chan[1], τ)
end
delay_trigger2 = Dict(3=>delay_trigger_affect2!)
delayjumpset2 = DelayJumpSet(delay_trigger2,delay_complete,delay_interrupt)

djprob2 = DelayJumpProblem(dprob, DelayRejection(), jumpset, delayjumpset2, de_chan0, save_positions=(false,false))

@time ens2 = solve(EnsembleProblem(djprob2), SSAStepper(), EnsembleThreads(), trajectories=10^5)

last_slice2 = componentwise_vectors_timepoint(ens2, tf)
histogram(last_slice2[3],bins=0:1:60,normalize=:pdf, label = "LogNormal(0,2)",fillalpha = 0.6, linecolor = :orange,fmt=:svg,xlabel = "# of products", ylabel = "Probability")
# savefig("docs/src/assets/stochastic_delay2.svg")


text:
# Copyright 2022 RelationalAI, Inc.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License

# Fetch details for the given OAuth client.

using RAI: Context, HTTPError, load_config, get_oauth_client

include("parseargs.jl")

function run(id; profile)
    cfg = load_config(; profile = profile)
    ctx = Context(cfg)
    rsp = get_oauth_client(ctx, id)
    println(rsp)
end

function main()
    args = parseargs(
        "id", Dict(:help => "OAuth client id", :required => true),
        "--profile", Dict(:help => "config profile (default: default)"))
    try
        run(args.id; profile = args.profile)
    catch e
        e isa HTTPError ? show(e) : rethrow()
    end
end

main()

max_stars_repo_name: RelationalAI/rai-sdk-julia

max_issues_repo_name: RelationalAI/rai-sdk-julia

max_forks_repo_name: RelationalAI/rai-sdk-julia


text:
@testset "FFT" begin
    ndv = NamedDimsArray(zeros(ComplexF64, 16), :μ)
    ndv[μ=2] = sqrt(2)

    @testset "vector" begin
        @test dimnames(fft(ndv, :μ => :κ)) == (:κ,)
        @test dimnames(ifft(ndv, :μ => :κ)) == (:κ,)
        @test dimnames(bfft(ndv, :μ => :κ)) == (:κ,)
        @test_throws ArgumentError fft(ndv, :nope => :κ)
        @test_throws ArgumentError ifft(ndv, :μ => :κ, :μ => :again)

        @test dimnames(fft(ndv)) == (:μ∿,)
        @test dimnames(ifft(ndv)) == (:μ∿,)
        @test dimnames(bfft(ndv)) == (:μ∿,)
        @test dimnames(ifft(fft(ndv))) == (:μ,)

        @test ifft(fft(ndv, :μ => :κ), :κ => :μ) ≈ ndv
        @test ifft(fft(ndv)) ≈ ndv
        @test bfft(fft(ndv)) ≈ 16 * ndv
    end

    @testset "vector plan" begin
        pv1 = plan_fft(ndv, :μ)
        pv2 = plan_fft(randn(ComplexF64, 16))
        @test pv2 isa FFTW.FFTWPlan

        @test dimnames(pv1 * ndv) == (:μ∿,)
        @test dimnames(pv2 * ndv) == (:μ∿,)
        @test dimnames(pv1 * randn(16)) == (:_,)

        pv3 = plan_ifft(pv2 * ndv)
        @test dimnames(pv3 * (pv2 * ndv)) == (:μ,)
        @test dimnames(inv(pv3) * ndv) == (:μ∿,)
        @test dimnames(inv(pv1) * (pv1 * ndv)) == (:μ,)

        @test dimnames(pv3 \ ndv) == (:μ∿,)
    end

    nda = NamedDimsArray(zeros(Float32, 4,4,4), (:a, :b′, :c))
    nda[1,2,3] = nda[2,3,4] = 1

    @testset "three dims" begin
        @test dimnames(fft(nda, :a => :k)) == (:k, :b′, :c)
        @test dimnames(ifft(nda, :a => :k, :c => :k′)) == (:k, :b′, :k′)

        @test dimnames(fft(nda)) == (:a∿, :b′∿, :c∿)
        @test dimnames(ifft(nda, 1)) == (:a∿, :b′, :c)
        @test dimnames(bfft(nda, :b′)) == (:a, :b′∿, :c)
        @test dimnames(fft(nda, 1:2)) == (:a∿, :b′∿, :c)
        @test dimnames(fft(nda, (:a, :c))) == (:a∿, :b′, :c∿)
    end

    @testset "three plan" begin
        p1 = plan_fft(nda, :a)
        p2 = plan_fft(zeros(Float32, 4,4,4), 2:2)
        @test_throws ArgumentError plan_fft(nda, :c => :C)
        p3 = plan_fft(nda, :c)

        @test dimnames(p1 * nda) == (:a∿, :b′, :c)
        @test dimnames(p2 * nda) == (:a, :b′∿, :c)
        @test dimnames(p3 * nda) == (:a, :b′, :c∿)
        @test dimnames(p2 * randn(4,4,4)) == (:_, :_, :_)

        @test dimnames(p3 * (p1 * nda)) == (:a∿, :b′, :c∿)

        @test_throws ArgumentError plan_fft(nda, :z => :Z)
    end

    @testset "wave_name" begin
        @test wave_name(:k) == :k∿
        @test wave_name((:k1, :k2∿)) == (:k1∿, :k2)
        @test wave_name((:k1, :k2, :k3), 2) == (:k1, :k2∿, :k3)
        @test wave_name((:k1, :k2, :k3), 1:2) == (:k1∿, :k2∿, :k3)
        @test wave_name((:k1, :k2, :k3), (1,3)) == (:k1∿, :k2, :k3∿)
    end

    @testset "other functions" begin
        nda = NamedDimsArray(rand(16, 16), (:x, :y))
        ndv = NamedDimsArray(randn(32), :z)

        for f in [fft, ifft, bfft, rfft]
            @test dimnames(f(nda)) == (:x∿, :y∿)
            @test dimnames(f(nda, :x)) == (:x∿, :y)
            @test dimnames(f(nda, :x => :k)) == (:k, :y)

            @test dimnames(f(ndv, 1)) == (:z∿,)
        end

        ndc = rfft(nda, 1) # size 9x16
        @test dimnames(irfft(ndc, 16, 1)) == (:x, :y)
        @test dimnames(brfft(ndc, 16, 1)) == (:x, :y)

        for f in [plan_fft, plan_ifft, plan_bfft]
            plan = f(nda, :x)
            @test plan isa AbstractFFTs.Plan
            @test dimnames(plan * nda) == (:x∿, :y)
        end

        @test dimnames(fftshift(nda)) == (:x, :y)
        @test fftshift(nda)[1,1] == nda[9,9]
        @test fftshift(nda, :x)[1,1] == nda[9,1]
    end
end

@testset "allocations: FFT" begin
    @test 0 == @ballocated wave_name(:k)
    @test 0 == @ballocated wave_name((:k1, :k2∿))
    if VERSION >= v"1.1"
        @test 0 == @ballocated wave_name((:k1, :k2, :k3), 2)
        @test 0 == @ballocated wave_name((:k1, :k2, :k3), (1,3))
    end
end


text:
function compute_tss(y::Vector{Float64}, hasintercept::Bool, sqrtw::AbstractVector)
    if hasintercept
        tss = zero(Float64)
        m = (mean(y) / sum(sqrtw) * length(y))::Float64
        @inbounds @simd for i in eachindex(y)
            tss += abs2(y[i] - sqrtw[i] * m)
        end
    else
        tss = sum(abs2, y)
    end
    return tss
end


text:
module CUSOLVER

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cuComplex, cuDoubleComplex, libraryPropertyType, cudaDataType
using ..CUDA: libcusolver, @allowscalar, assertscalar, unsafe_free!, @retry_reclaim

using ..CUBLAS: cublasFillMode_t, cublasOperation_t, cublasSideMode_t, cublasDiagType_t
using ..CUSPARSE: cusparseMatDescr_t

using CEnum

# core library
include("libcusolver_common.jl")
include("error.jl")
include("libcusolver.jl")

# low-level wrappers
include("util.jl")
include("wrappers.jl")

# high-level integrations
include("linalg.jl")

# thread cache for task-local library handles
const thread_dense_handles = Vector{Union{Nothing,cusolverDnHandle_t}}()
const thread_sparse_handles = Vector{Union{Nothing,cusolverSpHandle_t}}()

function dense_handle()
    tid = Threads.threadid()
    if @inbounds thread_dense_handles[tid] === nothing
        ctx = context()
        thread_dense_handles[tid] = get!(task_local_storage(), (:CUSOLVER, :dense, ctx)) do
            handle = cusolverDnCreate()
            finalizer(current_task()) do task
                CUDA.isvalid(ctx) || return
                context!(ctx) do
                    cusolverDnDestroy(handle)
                end
            end

            handle
        end
    end
    @inbounds thread_dense_handles[tid]
end

function sparse_handle()
    tid = Threads.threadid()
    if @inbounds thread_sparse_handles[tid] === nothing
        ctx = context()
        thread_sparse_handles[tid] = get!(task_local_storage(), (:CUSOLVER, :sparse, ctx)) do
            handle = cusolverSpCreate()
            finalizer(current_task()) do task
                CUDA.isvalid(ctx) || return
                context!(ctx) do
                    cusolverSpDestroy(handle)
                end
            end

            handle
        end
    end
    @inbounds thread_sparse_handles[tid]
end

function __init__()
    resize!(thread_dense_handles, Threads.nthreads())
    fill!(thread_dense_handles, nothing)

    resize!(thread_sparse_handles, Threads.nthreads())
    fill!(thread_sparse_handles, nothing)

    CUDA.atcontextswitch() do tid, ctx
        thread_dense_handles[tid] = nothing
        thread_sparse_handles[tid] = nothing
    end

    CUDA.attaskswitch() do tid, task
        thread_dense_handles[tid] = nothing
        thread_sparse_handles[tid] = nothing
    end
end

end


text:
using Test, DiffEqOperators

# Set up coefficient functions to test with
vec_fcn = Vector{Function}(undef,0)
f1(x::Float64) = sin(x)
f2(x::Vector{Float64}) = sin.(x)
push!(vec_fcn, f1)
push!(vec_fcn, f2)

@testset "Test coefficient functions when current_coeffs exists" begin

    vec_fcn_ans = Vector{Vector{Float64}}(undef,0)
    for i = 1:2
        push!(vec_fcn_ans, [sin(1.1), sin(2.3), sin(4.5)])
    end

    current_coeffs = [1.1, 2.3, 4.5]

    # Test inplace versions
    current_coeffs1 = Vector{Float64}(undef,3)
    @test DiffEqOperators.compute_coeffs!(2.54, current_coeffs1) ≈ 2.54 .* ones(3)
    current_coeffs1 = Vector{Float64}(undef,3)
    @test DiffEqOperators.compute_coeffs!([1.1, 2.3, 4.5], current_coeffs1) ≈ current_coeffs
    @test DiffEqOperators.compute_coeffs!(2.54, current_coeffs1) ≈ [3.64, 4.84, 7.04]
    @test all(current_coeffs .!= current_coeffs1)
    current_coeffs1 = copy(current_coeffs)
    @test DiffEqOperators.compute_coeffs!([2.54, 1.03, .18], current_coeffs1) ≈ [3.64, 3.33, 4.68]
    @test all(current_coeffs .!= current_coeffs1)
    for (fcn, fcn_ans) in zip(vec_fcn, vec_fcn_ans)
        current_coeffs1 = copy(current_coeffs)
        @test DiffEqOperators.compute_coeffs!(fcn, current_coeffs1) ≈ fcn_ans
        @test all(current_coeffs1 .!= current_coeffs)
    end

    # Test not in place versions
    current_coeffs1 = Vector{Float64}(undef,3)
    @test DiffEqOperators.compute_coeffs(2.54, current_coeffs) ≈ [3.64, 4.84, 7.04]
    @test DiffEqOperators.compute_coeffs([2.54, 1.03, 0.18], current_coeffs) ≈ [3.64, 3.33, 4.68]
    for (fcn, fcn_ans) in zip(vec_fcn, vec_fcn_ans)
        @test DiffEqOperators.compute_coeffs(fcn, current_coeffs) ≈ fcn_ans
    end
end

@testset "Check coefficients of DerivativeOperators are properly computed" begin

    # Set up operators to test construction with different coeff_func
    vec_op = Vector{DerivativeOperator}(undef,0)
    dxv = [1.4, 1.1, 2.3, 3.6, 1.5]
    push!(vec_op, UpwindDifference(1, 1, 1., 3, 1.))
    push!(vec_op, UpwindDifference(1, 1, 1., 3, [1., 2.25, 3.5]))
    push!(vec_op, UpwindDifference(1, 1, 1., 3, cos))
    push!(vec_op, UpwindDifference(1, 1, dxv, 3, 1.))
    push!(vec_op, UpwindDifference(1, 1, dxv, 3, [1., 2.25, 3.5]))
    push!(vec_op, UpwindDifference(1, 1, dxv, 3, cos))

    vec_op_ans = Vector{Any}(undef,0)
    push!(vec_op_ans, ones(3))
    push!(vec_op_ans, [1., 2.25, 3.5])
    push!(vec_op_ans, ones(3))
    push!(vec_op_ans, ones(3))
    push!(vec_op_ans, [1., 2.25, 3.5])
    push!(vec_op_ans, ones(3))

    # Check constructors
    @test UpwindDifference(1, 1, 1., 3).coefficients + ones(3) - ones(3) ≈ zeros(3)
    @test CenteredDifference(2, 2, 1., 3, 0.).coefficients + ones(3) - ones(3) ≈ zeros(3)
    for (L, op_ans) in zip(vec_op, vec_op_ans)
        @test L.coefficients ≈ op_ans
    end
    push!(vec_op, CenteredDifference(2, 2, 1., 3)) == nothing

    # Compute answers to coeff_func * operator
    func_mul_op_ans1 = Vector{Any}(undef,0)
    for i = 1:2
        push!(func_mul_op_ans1, zeros(3))
    end
    func_mul_op_ans2 = Vector{Any}(undef,0)
    for i = 1:2
        push!(func_mul_op_ans2, sin(1) .* ones(3))
    end
    func_mul_op_ans3 = Vector{Any}(undef,0)
    for i = 1:2
        push!(func_mul_op_ans3, [0., sin(1.5), 0.])
    end

    # Check coeff_func * L
    L1 = CenteredDifference(2, 2, 1., 3)
    L2 = UpwindDifference(1, 1, 1., 3, 1.)
    L3 = UpwindDifference(1, 1, 1., 3, [0., 1.5, π])
    for (fcn, fcn_ans) in zip(vec_fcn, func_mul_op_ans1)
        tmp = fcn * L1
        @test tmp.coefficients + ones(3) - ones(3) ≈ fcn_ans
    end
    for (fcn, fcn_ans) in zip(vec_fcn, func_mul_op_ans2)
        tmp = fcn * L2
        @test tmp.coefficients ≈ fcn_ans
    end
    for (fcn, fcn_ans) in zip(vec_fcn, func_mul_op_ans3)
        tmp = fcn * L3
        @test tmp.coefficients ≈ fcn_ans
    end

end

nothing


text:
@testset "Only T axis" begin
    planes = []

    dims = (X=512, Y=512, Z=1, C=1, T=10, P=1)
    axes = [Axis{:x}(), Axis{:y}(), Axis{:z}(), Axis{:channel}(1), Axis{:time}(1:10), Axis{:position}(1)]

    for t=0:9
        push!(planes, """<Plane DeltaT="$(t*10.0)" DeltaTUnit="ms" TheC="0" TheT="$t" TheZ="0"/>""")
    end
    container = root(parsexml(wrap(join(planes, "\n\t"))))

    @test all(OMETIFF.get_elapsed_times([container], dims, axes; default_unit=Unitful.ms).data .== collect(0:10:90)*u"ms")
end

@testset "Missing DeltaT" begin
    planes = []

    dims = (X=512, Y=512, Z=1, C=1, T=10, P=1)
    axes = [Axis{:x}(), Axis{:y}(), Axis{:z}(), Axis{:channel}(1), Axis{:time}(1:10), Axis{:position}(1)]

    for t=0:9
        push!(planes, """<Plane TheC="0" TheT="$t" TheZ="0"/>""")
    end
    container = root(parsexml(wrap(join(planes, "\n\t"))))

    @test all(isnan.(OMETIFF.get_elapsed_times([container], dims, axes; default_unit=Unitful.ms)))
end

@testset "Multiple dimensions" begin
    planes = []

    dims = (X=512, Y=512, Z=1, C=2, T=10, P=3)
    axes = [Axis{:x}(), Axis{:y}(), Axis{:z}(), Axis{:channel}(1:2), Axis{:time}(1:10), Axis{:position}(1:3)]

    containers = EzXML.Node[]
    for p in 1:3
        for t in 1:10, c in 1:2
            push!(planes, """<Plane DeltaT="$(t*p*c)" TheC="$(c-1)" TheT="$(t-1)" TheZ="0"/>""")
        end
        push!(containers, root(parsexml(wrap(join(planes, "\n\t")))))
    end

    output = OMETIFF.get_elapsed_times(containers, dims, axes; default_unit=Unitful.ms)
    @test all(size(output) .== (2, 10, 3))
    @test all(output[Axis{:position}(2), Axis{:channel}(1)] .== collect(2:2:20.0)*u"ms")
end

@testset "Mixed units" begin
    planes = []

    dims = (X=512, Y=512, Z=1, C=1, T=10, P=1)
    axes = [Axis{:x}(), Axis{:y}(), Axis{:z}(), Axis{:channel}(1), Axis{:time}(1:10), Axis{:position}(1)]

    for t in 1:10
        unit = "ms"
        (iseven(t)) && (unit = "s")
        push!(planes, """<Plane DeltaT="$((t+1)÷2*10.0)" DeltaTUnit="$unit" TheC="0" TheT="$(t-1)" TheZ="0"/>""")
    end
    container = root(parsexml(wrap(join(planes, "\n\t"))))

    output = OMETIFF.get_elapsed_times([container], dims, axes; default_unit=Unitful.ms)
    @test all(output[1:2:end]*1000 .== output[2:2:end]) # the evens are in seconds so they should be 1000 times the odd values
end