text:
using DimensionalData, Test, Unitful
using DimensionalData: X, Y, Z, Time, Forward, Reverse,
      reversearray, reverseindex, slicebounds, slicegrid

@testset "sampling" begin

end

@testset "reverse" begin
    @test reverse(Reverse()) == Forward()
    @test reverse(Forward()) == Reverse()

    @test reversearray(Unordered(Forward())) ==
        Unordered(Reverse())
    @test reversearray(Ordered(Forward(), Reverse(), Forward())) ==
        Ordered(Forward(), Forward(), Reverse())

    @test reverseindex(Unordered(Forward())) ==
        Unordered(Reverse())
    @test reverseindex(Ordered(Forward(), Reverse(), Forward())) ==
        Ordered(Reverse(), Reverse(), Reverse())

    @test order(reverseindex(AlignedGrid(order=Ordered(Forward(), Reverse(), Forward())))) ==
        Ordered(Reverse(), Reverse(), Reverse())
    @test order(reversearray(AlignedGrid(order=Ordered(Forward(), Reverse(), Forward())))) ==
        Ordered(Forward(), Forward(), Reverse())
    @test order(reverseindex(RegularGrid(order=Ordered(Forward(), Reverse(), Reverse())))) ==
        Ordered(Reverse(), Reverse(), Forward())
    @test order(reverseindex(CategoricalGrid(order=Ordered(Forward(), Reverse(), Reverse())))) ==
        Ordered(Reverse(), Reverse(), Forward())
end

@testset "slice bounds" begin
    index = [10.0, 20.0, 30.0, 40.0, 50.0]
    bound = (10.0, 60.0)
    @test slicebounds(Start(), bounds, index, 2:3) == (20.0, 40.0)
    bound = (0.0, 50.0)
    @test slicebounds(End(), bounds, index, 2:3) == (10.0, 30.0)
    bound = (0.5, 55.0)
    @test slicebounds(Center(), bounds, index, 2:3) == (15.0, 35.0)
    grid = BoundedGrid(; locus=Start(), bounds=(10.0, 60.0))
    @test bounds(slicegrid(grid, index, 3)) == (30.0, 40.0)

    dim = X(index; grid=RegularGrid(; locus=Start(), span=10.0))
    @test bounds(dim) == (10.0, 60.0)
    dim = X(index; grid=RegularGrid(; locus=End(), span=10.0))
    @test bounds(dim) == (0.0, 50.0)
    dim = X(index; grid=RegularGrid(; locus=Center(), span=10.0))
    @test bounds(dim) == (5.0, 55.0)

    index = [:a, :b, :c, :d]
    dim = X(index; grid=CategoricalGrid(; order=Ordered()))
    @test bounds(dim) == (:a, :d)
    dim = X(index; grid=CategoricalGrid(; order=Ordered(;index=Reverse())))
    @test bounds(dim) == (:d, :a)
    dim = X(index; grid=CategoricalGrid(; order=Unordered()))
    @test_throws ErrorException bounds(dim)
end


text:
# Hello world example, similar to the Boost.Python hello world

using CxxWrap
using Base.Test
using Compat

# Wrap the functions defined in C++
wrap_modules(CxxWrap._l_types)

using CppTypes
using CppTypes.World

# Stress test
for i in 1:1000000
  d = CppTypes.DoubleData()
end

# Default constructor
@test World <: CxxWrap.CppAny
@test supertype(World) == CxxWrap.CppAny
w = World()
println("Dumping type w...")
dump(w)
@test CppTypes.greet(w) == "default hello"

@show fw = CppTypes.world_factory()
@test CppTypes.greet(fw) == "factory hello"

@show swf = CppTypes.shared_world_factory()
@test CppTypes.greet(swf) == "shared factory hello" # Uses the shared ptr overload
@test CppTypes.greet(CppTypes.get(swf)) == "shared factory hello" # Explicitly get the shared ptr

@show uwf = CppTypes.unique_world_factory()
@test CppTypes.greet(CppTypes.get(uwf)) == "unique factory hello"

CppTypes.set(w, "hello")
@show CppTypes.greet(w)
@test CppTypes.greet(w) == "hello"

w = World("constructed")
@test CppTypes.greet(w) == "constructed"

w_assigned = w
w_deep = deepcopy(w)

@test w_assigned == w
@test w_deep != w

# Destroy w: w and w_assigned should be dead, w_deep alive
finalize(w)
println("finalized w")
if !(is_windows() && Sys.WORD_SIZE == 32)
  @test_throws ErrorException CppTypes.greet(w)
  println("throw test 1 passed")
  @test_throws ErrorException CppTypes.greet(w_assigned)
  println("throw test 2 passed")
end
@test CppTypes.greet(w_deep) == "constructed"
println("completed deepcopy test")

noncopyable = CppTypes.NonCopyable()
println("noncopyable constructed")
if !(is_windows() && Sys.WORD_SIZE == 32)
  @test_throws ErrorException other_noncopyable = deepcopy(noncopyable)
end

println("completed noncopyable test")

import CppTypes.ImmutableDouble

@test sizeof(ImmutableDouble) == 8
@test isbits(ImmutableDouble)
@test length(fieldnames(ImmutableDouble)) == 1
println("creating bitsval1")
bitsval1 = ImmutableDouble(1)
println("created bitsval1")
@test bitsval1.value == 1.
@test bitsval1 == 1.
@test CppTypes.getvalue(bitsval1) == 1
bitsval2 = CppTypes.ImmutableDouble(2)
@test bitsval2 == 2
@test typeof(bitsval1 + bitsval2) == CppTypes.ImmutableDouble
@test (bitsval1 + bitsval2) == 3.

println("creating bc")
bc = make_bits(1, 2)
println("created bc")
@test sizeof(bc)==16
@test get_bits_a(bc)==1
@test get_bits_b(bc)==2

@test value(value(ReturnConstRef())) == 42

@test CppTypes.greet(ConstPtrConstruct(World())) == "default hello"

call_op = CallOperator()
@test call_op() == 43

type JuliaTestType
    a::Float64
    b::Float64
end

function julia_test_func(data)
  println("a: ", data.a, ", b: ", data.b)
  @test data.a == 2.
  @test data.b == 3.
end

CppTypes.call_testype_function()

ext: jl


text:
# test the "tangent projection" approach in the somewhat trival 3-site uniform mps
# roughly equivalent to the appraoch used in PhysRevB.94.035133

using TensorKit
using TensorOperations
using TensorKitAD
using Zygote
using KrylovKit
using LinearAlgebra
using Revise
using statmech_tm_solver
using Optim

function fovlp(A::TensorMap{ComplexSpace, 2, 1}, B::TensorMap{ComplexSpace, 2, 1})
    @tensor A[3, 7, 1] * A[1, 8, 2] * A[2, 9, 3] * B'[4, 6, 7] * B'[5, 4, 8] * B'[6, 5, 9]
end

function fidel(A::TensorMap{ComplexSpace, 2, 1}, B::TensorMap{ComplexSpace, 2, 1})
    fovlp(A, B) * fovlp(B, A) / fovlp(A, A) / fovlp(B, B) |> real
end

Afull = 0.1*TensorMap(rand, ComplexF64, ℂ^8*ℂ^2, ℂ^8)
A = TensorMap(rand, ComplexF64, ℂ^2*ℂ^2, ℂ^2)
Iph = id(ℂ^2)

A = A / fovlp(A, A)^(1/6)

for ix in 1:75
    @tensor G3[-1, -2, -3; -4, -5, -6] :=
        A[-4, 1, 2] * A'[3, -3, 1] * A[2, 4, -5] * A'[-1, 3, 4] * Iph[-2, -6] +
        A[2, 1, -5] * A'[3, -3, 1] * A'[-1, 3, -6] * A[-4, -2, 2] +
        A'[3, -3, -6] * A[-4, 1, 2] * A'[-1, 3, 1] * A[2, -2, -5] ;

    _, S, _ = tsvd(G3)

    G3_inv = pinv(G3, rtol=1e-12, atol=1e-12);
    G3_inv = permute(G3_inv, (6, 2, 3), (4, 5, 1))

    @tensor A_new[-1, -2; -3] := Afull[6, 7, 4] * Afull[4, 8, 5] * Afull[5, 9, 6] * A'[1, 3, 7] * A'[3, 2, 9] * G3_inv[2, -1, -2, 1, 8, -3]

    fidel0 = fidel(Afull, A)
    coeff = 1.0

    while fidel(Afull, A_new*coeff+A*(1-coeff)) < fidel0
        coeff = 0.5*coeff
    end
    #@show ix, coeff, fidel(Afull, A_new*coeff+A*(1-coeff))

    A = coeff*A_new + A*(1-coeff)
end
@show fidel(Afull, A)
@show A

A2_arr = rand(ComplexF64, (2,2,2))
function f2(X::Array{<:Number, 3})
    X_TensorMap = convert_to_tensormap(X, 2)
    return -fidel(Afull, X_TensorMap)
end
function g2!(G::Array{<:Number, 3}, X::Array{<:Number, 3})
    G .= gradient(f2, X)[1]
end

res = optimize(f2, g2!, A2_arr,  LBFGS(), Optim.Options(iterations=100) )


text:
module NuclearModelTests
using AtomicMiscellany.NuclearModels
using Test
using QuadGK

maxnonzero(::AbstractNuclearModel) = Inf
maxnonzero(m::UniformSphericalNucleus) = m.R

@testset "Nuclear models" begin
    # Test that any AbstractNuclearModel broadcasts like a scalar
    @test length(potential.(PointNucleus(), [1,2,3])) == 3
    @test length(density.(GaussianNucleus(1), [1,2,3,4])) == 4

    # Check that if we construct the nuclear models using the from_rms method, that the
    # resulting objects have the correct RMS value.
    function check_rms_constructor(T::Type{<:AbstractNuclearModel})
        @testset "check_rms_constructor: $T" begin
            for _rms in [3.1e-6, 1e-4, 0.4, 1.0, 10.0]
                m = from_rms(T, _rms)
                @test rms(m) ≈ _rms
            end
        end
    end

    @testset "from_rms" begin
        @test_throws MethodError from_rms(PointNucleus, 0.0)
        check_rms_constructor(UniformShellNucleus)
        check_rms_constructor(UniformSphericalNucleus)
        check_rms_constructor(GaussianNucleus)
        #check_rms_constructor(FermiNucleus)
    end

    # Use direct integration to calculate the RMS values and compare against the implemented
    # equations.
    function compare_integrals(m::AbstractNuclearModel)
        rmax, rscale = maxnonzero(m), rms(m)
        @testset "compare_integrals: $m" begin
            # We need to scale and split the integrals: see the comments for potential_quadk()
            # below.
            ρ₀, δ = let f(t) = 4π * rscale^3 * t^2 * density(m, rscale*t)
                rmax/rscale < 10.0 ? quadgk(f, 0, rmax/rscale) : quadgk(f, 0, 10.0, rmax/rscale)
            end
            @test ρ₀ ≈ 1.0

            r², δ = let f(t) = 4π * rscale^5 * t^4 * density(m, rscale*t)
                rmax/rscale < 10.0 ? quadgk(f, 0, rmax/rscale) : quadgk(f, 0, 10.0, rmax/rscale)
            end
            @debug "$r² ± $δ" sqrt(r²) rms(m) rms(m)^2
            @test sqrt(r²) ≈ rms(m)
        end
    end

    @testset "r^k integrals" begin
        for R in [1, 1e-4, 5.1e-6, 35, π, 5π^4]
            compare_integrals(UniformSphericalNucleus(R))
            compare_integrals(GaussianNucleus(R))
            #compare_integrals(FermiNucleus(R))
        end
    end

    # Test to make sure that the density and the potential correspond to each other, for the
    # nuclear models that we have the density implemented. We'll use QuadGK to numerically
    # integrate the potential values from the density directly and then check against the
    # implementation of the potential.
    function potential_quadk(ρ, r; rmax = Inf, rscale = 1.0)
        # Based on equation (22) in Andrae. The potential is split into two separate integrals
        # on the intervals [0, r] and [r, ∞).
        #
        # However, often the density can have its non-trivial parts at very small r values,
        # which means that the adaptive quadgk algorithm might not realize that it even has
        # non-zero parts. To work around that, we do 2 things:
        #
        # 1. We scale the r-axis with rscale (so r = rscale*t), and actually integrate in t.
        #    In practice, we'll use the RMS value for the scaling, since we'd expect that to
        #    correspond to the scale at which the density has a non-trivial shape.
        # 2. If r is much greater than rscale, we split the [0, r] integration into two
        #    sub-intervals, hoping that we capture all the non-trivial part in the first
        #    interval.
        #
        # The change of variable r = rscale*t leads to
        #
        #   ∫ f(r) dr = rscale ∫ f(rscale*t) dt
        #
        # So we have to also multiply the resulting integral with rscale.
        _I1(t) = rscale^3 * t^2 * ρ(rscale*t)
        I1, δ1 = min(r, rmax)/rscale > 10.0 ? quadgk(_I1, 0, 10.0, min(r, rmax)/rscale) :
            quadgk(_I1, 0, min(r, rmax)/rscale)
        I2, δ2 = (r < rmax) ? quadgk(t -> rscale^2 * t * ρ(rscale*t), r/rscale, rmax/rscale) : (0.0, 0.0)
        # add up the subintegral based on equation (22):
        -4π*(I1/r + I2)
    end

    test_potential(m::AbstractNuclearModel) = @testset "test_potential: $m" begin
        for r0 in 10 .^ range(-6, 2, length=21)
            V = potential_quadk(r -> density(m, r), r0, rmax = maxnonzero(m), rscale = rms(m))
            @test V ≈ potential(m, r0)
        end
    end

    @testset "V(r) from ρ(r)" begin
        for R in [1, 1e-4, 5.1e-6, 35, π, 5π^4]
            test_potential(UniformSphericalNucleus(R))
            test_potential(GaussianNucleus(R))
            #test_potential(FermiNucleus(R))
        end
    end
end

end


text:
@testset "CenteredAxis" begin
    ca = @inferred(CenteredAxis(1:10))
    @test @inferred(keys(ca)) == -5:4
    @test @inferred(parent(ca)) == 1:10
    @test typeof(ca)(parent(ca)) isa typeof(ca)
    ca = @inferred(CenteredAxis{Int32}(UnitRange(1, 10)))
    @test typeof(ca)(parent(ca)) isa typeof(ca)
    @test keytype(ca) <: Int32
    centered_axis = @inferred(CenteredAxis{Int32}(UnitRange(1, 10)))
    @test eltype(ca) <: Int32
    ca2 = ca[-1:1]
    @test @inferred(keys(ca2)) == -1:1
    @test @inferred(parent(ca2)) == 5:7
end


text:
using Plots, LaTeXStrings, Plots.Measures
using LinearAlgebra, MAT

(@isdefined TARGET_FOLDER) ? nothing : TARGET_FOLDER = ""

# load the model
include("Heat1D_Model.jl")

import IntervalArithmetic
const IA = IntervalArithmetic

file = joinpath(@__DIR__, "heat1d_99.mat")

δ = 1e-5
k = 1_000
tol = 0.1

prob, alg, NSTEPS = _solve_heat1d(NSTEPS=1, δ=δ , file=file, model=StepIntersect(Forward()))
n = statedim(prob);

LGG_solu_dirs = LGG_dirs_solucion(n)
LGG_grad_dirs = LGG_dirs_gradiente(n)

L = 1.0
Δx = L/(n+2)

# warm up
solve(prob, alg=LGG09(δ=δ, template=LGG_grad_dirs, approx_model=Forward(setops=:box)), NSTEPS=1)
varias_sol(file=file, δ=δ, tol=tol, NSTEPS=1, NCURVAS=1);

println("Flowpipe solution with gradient directions (LGG)")
@time grad_lgg = solve(prob, alg=LGG09(δ=δ, template=LGG_grad_dirs, approx_model=Forward(setops=:box)), NSTEPS=k);

NCURVAS = 1_000
println("Trajectories solution (Euler's method, $NCURVAS Trajectories)")
@time Temps, gradsT = varias_sol(file=file, δ=δ, tol=tol, NSTEPS=k, NCURVAS=NCURVAS);

# spatial gradient
function plot_gradiente_espacial(; paso = 1_000, δ = 1E-5, NCURVAS=length(gradsT), extremos_y = 6.0)

    fig = plot(xlab=L"\textrm{Position}",
               ylab=L"\textrm{Temperature gradient}",
               legendfontsize=18, legend=:top,
               tickfont=font(25, "Times"), guidefontsize=20,
               xguidefont=font(30, "Times"),
               yguidefont=font(30, "Times"),
               ylims=(-extremos_y-0.5, extremos_y),
               xlims=(0.0, 1.0),
               xtick=([0.0, 0.25, 0.50, 0.75, 1.0], [L"0.0", L"0.25", L"0.50", L"0.75", L"1.0"]),
               ytick=([-10.0, -5.0, 0.0, 5.0, 10.0], [L"-10", L"-5", L"0", L"5", L"10"]),
               bottom_margin=6mm, left_margin=6mm, right_margin=4mm, top_margin=6mm, size=(900, 600))

    fp_lgg_grad = gradiente_LGG_graddirs(grad_lgg = grad_lgg, paso_temporal = paso)
    plot!(fig, fp_lgg_grad, vars=(0,1), lw=.5, c=:lightblue)

    xsnumer = range(0, L, length=length(gradsT[1][1])+1)
    xsnumer = 0.5*(xsnumer[1:end-1]+xsnumer[2:end])
    for j = 1:NCURVAS
        plot!(fig, xsnumer, gradsT[j][paso], c=:magenta, lw=1.0, lab="")
    end
    fig

end

num_paso = 100
fig = plot_gradiente_espacial(paso = num_paso, extremos_y = 10.0)

savefig(fig, joinpath(TARGET_FOLDER, "grad_temp_barra_paso$num_paso.pdf"))

# Temporal gradient
function plot_gradiente_temporal(; nodo_grad = 66, δ = 1E-5, kmax = length(grad_lgg), NCURVAS=length(gradsT), extremos_y = [-20, 25])

    fig = plot(xlab=L"\textrm{Time } (\times 10^{-4})",
               ylab=L"\textrm{Temperature gradient}",
               legendfontsize=18, legend=:top,
               tickfont=font(25, "Times"), guidefontsize=20,
               xguidefont=font(30, "Times"),
               yguidefont=font(30, "Times"),
               ylims=(extremos_y[1], extremos_y[2]),
               xlims=(0.0, kmax*δ + 0.1e-4),
               xticks=([0.0, 0.00025, 0.0005, 0.00075, 0.001], [L"0", L"2.5", L"5.0", L"7.5", L"10.0"]),
               yticks=([-20, -10, 0, 10, 20], [L"-20", L"-10", L"0", L"10", L"20"]),
               bottom_margin=6mm, left_margin=6mm, right_margin=4mm, top_margin=3mm, size=(900, 600))

    kmax = length(grad_lgg)
    xsnumer = range(0, L, length=length(gradsT[1][1])+1)
    xsnumer = 0.5*(xsnumer[1:end-1]+xsnumer[2:end])
    x_nodo_grad = xsnumer[nodo_grad]

    fp_lgg_grad = gradiente_LGG_graddirs(grad_lgg = grad_lgg, paso_temporal = 1)

    fp_grad_nodo = [ReachSet(set(fp_lgg_grad(x_nodo_grad)), tspan(grad_lgg[1]))]
    for j in 2:kmax
        fp_lgg_grad = gradiente_LGG_graddirs(grad_lgg = grad_lgg, paso_temporal = j)
        push!(fp_grad_nodo, ReachSet(set(fp_lgg_grad(x_nodo_grad)), tspan(grad_lgg[j])))
    end
    fff = Flowpipe(fp_grad_nodo)

    plot!(fig, fff, vars=(0,1), lw=0.5, lc=:black, c=:lightblue)

    tiempo = range(0., step=δ, length=kmax)
    long_t = length(gradsT[1])
    for j in 1:NCURVAS
        y = [gradsT[j][i][nodo_grad] for i in 1:long_t]
        plot!(fig, tiempo, y[1:end-1], c=:magenta, lab="", lw=0.5)
    end

    fig
end

nnodo = 66
fig = plot_gradiente_temporal(nodo_grad=nnodo, kmax=100)
savefig(fig, joinpath(TARGET_FOLDER, "grad_temp_tiempo_paso$nnodo.pdf"))

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
"""
# Module RWG
This module defines the modified Rao, Wilton, and Glisson triangle
subdomain basis function derived data type.  The definition of a basis
function has been generalized slightly to allow the "plus" and "minus"
triangles to be nonadjacent, and to identify those half-basis functions
whose defining edge lies at the ξ=1 or η=1 boundaries of the unit cell.
These modifications allow one to handle the analysis of a structure in a
periodic unit cell, as required in phased array and frequency selective
surface formulations.  Reference: P. S. Simon, "Modified RWG basis functions
for analysis of periodic structures," *2002 IEEE MTT-S International Microwave
Symposium Digest* (Cat. No.02CH37278), Seattle, WA, USA, 2002, pp. 2029-2032
vol. 3, doi: 10.1109/MWSYM.2002.1012266.
"""
module RWG
export RWGData, setup_rwg, edge_current_unit_vector, rwgbfft!

using ..Sheets: RWGSheet, SV2
using StaticArrays: SVector, MVector, SArray
using LinearAlgebra: ⋅, norm
using NearestNeighbors: KDTree, inrange
using OffsetArrays
using Statistics: mean
using ..PSSFSSLen

mutable struct RWGData
    #  Basis function edge indices.  The value in bfe[1,i] is index of the
    #  edge of the "plus" triangle associated with basis function i.
    #  The value in bfe[2,i] is the index of the edge of the "minus" triangle
    #  associated with basis function i. In most cases, these
    #  two values (edges) will be identical.
    bfe::Array{Int,2}

    #  Basis function face indices.  The value in bff[1,i] is index of the
    #  "plus" triangle face associated with basis function i.
    #  The value in bff[2,i] is the index of the "minus" triangle face
    #  associated with basis function i.
    bff::Array{Int,2}

    #  Edge basis function map.
    #  ebf[i], if nonzero, is the index of the basis function associated with
    #  edge i.
    ebf::Array{Int,1}

    #  Edge cell index.
    #  eci[i] takes on the values 0, 1, 2, 3, or 4.   The values have the
    #  following meanings:
    #  0  The edge does not lie along a unit cell boundary.
    #  1  The edge lies along the ξ=0 unit cell boundary.
    #  2  The edge lies along the ξ=1 unit cell boundary.
    #  3  The edge lies along the η=0 unit cell boundary.
    #  4  The edge lies along the η=1 unit cell boundary.
    eci::Array{Int,1}

    #  Unique face pair matrix.
    #  ufpm[i,j] contains the unique face pair index for observation face i
    #  with respect to source face j.  Two sets of face pairs are considered
    #  to be equivalent if they can be made to superimpose via a rigid
    #  translation, and the nodes in the corresponding triangles of each pair
    #  are numbered in the same order. Since the periodic Green's functions are
    #  translationally invariant, the integrals involving equivalent face pairs
    #  will have identical values and thus need be computed only a single time.
    ufpm::Array{Int,2}

    #  ufp2fp[i] contains the vector of face pair indices for equivalence class i. The
    #  face pair index uses column major ordering to enumerate the elements
    #  of a matrix of order Nface × Nface.
    ufp2fp::Array{Array{Int,1},1}

    zorymat::Array{ComplexF64,2} # MoM matrix
    rhs::Vector{ComplexF64}  # MoM right-hand side
    bfftstore::Array{SArray{Tuple{2},ComplexF64,1,2}, 3} # Basis function Fourier Transforms

    nufp::Int  # Number of unique face pairs.
end # mutable struct

"""
    setup_rwg(sheet::RWGSheet)::RWGdata

This function accepts the sheet geometry data structure as created
by the function `get_sheet_data` and creates a instance of `RWGdata` as the
function return value.  When `sheet.fufp` is `true`, it directs this
function to search for unique face pairs.  This search can be time consuming
and thus the default action, when `sheet.fufp` is .false.,
is to skip the search.  The tradeoff is the greater time needed to fill the
interaction matrix when all face pairs are considered unique.
"""
function setup_rwg(sheet::RWGSheet, leafsize::Int=9)::RWGData
    tol = 1.e-5 # Comparison tolerance
    ieη0 = Int[] # List edges at η=0
    ieη1 = Int[] # List edges at η=1
    ieξ0 = Int[] # List edges at ξ=0
    ieξ1 = Int[] # List edges at ξ=1

    nface = size(sheet.fe, 2) # Number of faces in triangulated sheet.
    nedge = length(sheet.e1) # Number of edges in triangulated sheet.

    eci = zeros(Int, nedge)

    # Count the number of basis functions:
    nbf = 0
    # To begin with, we locate any edge that is adjacent to a pair of faces...
    for fp in 1:nface, fm in fp+1:nface  # Loop over "plus" and "minus" triangles
        for gep in @view sheet.fe[:,fp]  # Global index of "plus" edge.
            for gem in @view sheet.fe[:,fm] # Global index of "minus" edge.
                gep == gem && (nbf += 1)    # Faces share common edge
            end
        end
    end
    # We must now add to this count the number of edges that are located at
    # the ξ=0  or the η=0 boundaries of the unit cell.  In the process of
    # locating such edges, we will also set the correct values in the eci,
    # ieξ0, ieξ1, ieη0, and  ieη1 arrays.
    if sheet.ξη_check
        for gep in 1:nedge  # Loop over each edge of the structure.
            # Obtain the ξ and η values for each vertex of edge ie:
            (ξinit, ξterm, ηinit, ηterm) = get_ie_ξη(gep, sheet)
            if abs(ξinit) < tol  &&  abs(ξterm) < tol
                # Edge is located on the ξ=0 boundary of unit cell.
                nbf += 1         # Bump count of # basis functions.
                eci[gep] = 1     # Store edge code.
                push!(ieξ0, gep) # Update list of edges at ξ = 0.
            elseif abs(ξinit-1) < tol  && abs(ξterm-1) < tol
                # Edge is located on the ξ=1 boundary of unit cell.
                eci[gep] = 2
                push!(ieξ1, gep)  # Update list of edges at ξ = 1.
            elseif abs(ηinit) < tol  && abs(ηterm) < tol
                # Edge is located on the η=0 boundary of unit cell.
                nbf += 1
                eci[gep] = 3
                push!(ieη0, gep)  # Update list of edges at η = 0.
            elseif abs(ηinit-1) < tol && abs(ηterm-1) < tol
                # Edge is located on the η=1 boundary of unit cell.
                eci[gep] = 4          # Store edge code.
                push!(ieη1, gep)  # Update list of edges at η = 1.
            end
        end # for
    end # if
    length(ieξ0) == length(ieξ1) || error("Inconsistent # edges at ξ=0 and ξ=1")
    length(ieη0) == length(ieη1) || error("Inconsistent # edges at η=0 and η=1")

    bfe = zeros(Int, 2,nbf)
    bff = zeros(Int, 2,nbf)
    ebf = zeros(Int, nedge)

    # Loop over pairs of triangles.  Note that the order of the loops dictates
    # that the face with lower index will be the "plus" face.
    i = 0  # Basis function index
    for fp in 1:nface, fm = fp+1:nface            # Loop over "plus" and "minus" triangles.
        for gep in view(sheet.fe, :, fp), gem in view(sheet.fe, :, fm) # Global edge index for "+" and "-" triangles.
            if gep == gem   # Faces share common edge:
                i += 1        # Bump basis function index.
                bfe[1,i] = gep # "Plus" triangle edge index for i'th basis funct
                bfe[2,i] = gem # "Minus" triangle edge index for i'th basis funct
                bff[1,i] = fp  # "Plus" triangle face index for i'th basis funct
                bff[2,i] = fm  # "Minus" triangle face index for i'th basis funct
                ebf[gep] = i
            end
        end
    end

    #  Now treat the basis functions located at unit cell boundaries.
    for ge0 in ieη0  # Loop over edges at η=0 unit cell boundary.
        # Search for the triangle face containing edge ge0:
        face0 = 0  # Initialize index of face containing edge at η=0.
        for j0 in 1:nface
            if ge0 in @view sheet.fe[:,j0]
                face0 = j0  # Save face index.
                break
            end
        end
        face0 == 0 && error("Unable to locate face containing η=0 edge $ge0")

        # Obtain the initial and terminal ξ values of edge at η = 0:
        (ξinit, ξterm, _, _) = get_ie_ξη(ge0, sheet)
        ξmin0, ξmax0 = extrema((ξinit, ξterm))

        # Search for the corresponding edge at η = 1:
        ge1 = 0   # Initialize global index of corresponding edge at η = 1.
        for ge in ieη1  # Loop over each edge at the η=1 unit cell boundary.
            # Obtain the initial and terminal ξ values:
            (ξinit, ξterm, _, _) = get_ie_ξη(ge, sheet)
            ξmin1, ξmax1 = extrema((ξinit, ξterm))
            # Search for edges occupying the same ξ interval:
            if abs(ξmin1-ξmin0) < tol && abs(ξmax1-ξmax0) < tol
                ge1 = ge     # Save global edge index for edge at η = 1.
                break
            end
        end
        ge1 == 0 && error("Unable to locate edge at η=1 corresponding to edge $ge0")
        # Search for the triangle face containing edge ge1:
        face1 = 0  # Initialize index of face containing edge at η=1.
        for j1 in 1:nface
            if ge1 in @view sheet.fe[:,j1]
                face1 = j1  # Save face index.
                break        # Jump out of search loop.
            end
        end
        face1 == 0 && error("Unable to locate face containing η=1 edge $ge1")
        i += 1      # Bump basis function index.
        bfe[1,i] = ge0 # "Plus" triangle edge index for i'th basis function.
        bfe[2,i] = ge1 # "Minus" triangle edge index for i'th basis function.
        bff[1,i] = face0  # "Plus" triangle face index for i'th basis function.
        bff[2,i] = face1  # "Minus" triangle face index for i'th basis function.
        ebf[ge0] = i
        ebf[ge1] = i
    end # for edges at η=0
    #
    for ge0 in ieξ0  # Loop over edges at ξ=0 unit cell boundary.
        # Search for the triangle face containing edge ge0:
        face0 = 0  # Initialize index of face containing edge at η=0.
        for j0 = 1:nface
            if ge0 in @view sheet.fe[:,j0]
                face0 = j0  # Save face index.
                break
            end
        end
        face0 == 0 && error("Unable to locate face containing ξ=0 edge $ge0")
        # Obtain the initial and terminal η values of edge at ξ = 0:
        (_, _, ηinit, ηterm) = get_ie_ξη(ge0, sheet)
        ηmin0, ηmax0 = extrema((ηinit, ηterm))
        # Search for the corresponding edge at ξ = 1:
        ge1 = 0   # Initialize global index of corresponding edge at ξ = 1.
        for ge in ieξ1  # Loop over each edge at the ξ=1 unit cell boundary.
            # Obtain the initial and terminal η values:
            (_, _, ηinit, ηterm) = get_ie_ξη(ge, sheet)
            ηmin1, ηmax1 = extrema((ηinit, ηterm))
            # Search for edges occupying the same η interval:
            if abs(ηmin1-ηmin0) < tol && abs(ηmax1-ηmax0) < tol
                ge1 = ge     # Save global edge index for edge at ξ = 1.
                break
            end
        end
        ge1 == 0 && error("Unable to locate edge at ξ=1 corresponding to edge $ge0")
        # Search for the triangle face containing edge ge1:
        face1 = 0  # Initialize index of face containing edge at η=1.
        for j1 = 1:nface
            if ge1 in @view sheet.fe[:,j1]
                face1 = j1  # Save face index.
                break
            end
        end
        face1 == 0 && error("Unable to locate face containing ξ=1 edge $ge1")
        i += 1      # Bump basis function index.
        bfe[1,i] = ge0 # "Plus" triangle edge index for i'th basis function.
        bfe[2,i] = ge1 # "Minus" triangle edge index for i'th basis function.
        bff[1,i] = face0  # "Plus" triangle face index for i'th basis function.
        bff[2,i] = face1  # "Minus" triangle face index for i'th basis function.
        ebf[ge0] = i
        ebf[ge1] = i
    end # ξ0_loop

    i == nbf || error("Inconsistent number of basis functions")
    zorymat = zeros(ComplexF64, nbf, nbf)
    rhs = zeros(ComplexF64, nbf)
    bfftstore = zeros(SArray{Tuple{2},ComplexF64,1,2}, 0,0,0)
    nufp =  nface*nface
    ufpm = reshape(collect(1:nufp), (nface,nface))
    if !sheet.fufp  # Don't search for unique face pairs. Assume all are unique.
        ufp2fp = [ [i] for i in 1:nufp]
        return RWGData(bfe, bff, ebf, eci, ufpm, ufp2fp, zorymat, rhs, bfftstore, nufp)
    end

    #  The remaining code in this function sets up the
    #  two arrays ufpm and ufp2fp.  These are used to identify
    #  redundant face/pairs.  ufpm(i) contains the unique face/pair
    #  index (or the equivalence class (E.C.) index) of face/pair i.
    #  The face pairs are numbered sequentially in the same manner as
    #  the elements of a square matrix of dimension nface (column major
    #  ordering). The row index is used for the observation face number and
    #  the column index is the source face number.  The array ufp2fp is
    #  of length nufp (the number of unique face pairs, or equivalence
    #  classes). Each element of ufp2fp contains a pointer to an allocated
    #  array of face/pair indices.  Row i contains a list of the face/pairs
    #  that are members of equivalence class i.  Two face pairs belong to
    #  the same equivalence class if the source triangles can be overlaid
    #  using a rigid translation, the source triangle nodes are listed
    #  in the same order, and if the observation point (centroid
    #  of the match triangle) is in the same relative position wrt the
    #  source triangle.  This version uses the NearestNeighbors package:

      #  Allocate the unique face pairs matrix and some scratch arrays.
    ufpm = zeros(Int, (nface,nface))

    r1 = @view sheet.ρ[sheet.fv[1,:]]
    r2 = @view sheet.ρ[sheet.fv[2,:]]
    r3 = @view sheet.ρ[sheet.fv[3,:]]
    centroid = [(r1[n] + r2[n] + r3[n]) / 3 for n in eachindex(r1)]
    data = Array{SVector{6,Float64}}(undef, nface^2)
    #  Characterize each face pair:
    mn = 0
    for n in 1:nface, m in 1:nface
        mn += 1 # Bump single face/face index.
        # Calculate test vectors for face pair (m,n):
        rmn1 = centroid[m] - r1[n]
        rmn2 = centroid[m] - r2[n]
        rmn3 = centroid[m] - r3[n]
        data[mn] = SVector{6,Float64}(rmn1[1], rmn1[2], rmn2[1], rmn2[2], rmn3[1], rmn3[2])
    end

    kdtree = KDTree(data, leafsize=leafsize)
    ufp2fp = Vector{Int}[]
    mn = 0  # Initialize face/pair index.
    nufp = 0
    found = fill(false, nface^2)
    r = 1e-6 * norm(centroid[1] - centroid[2])
    for n in 1:nface, m in 1:nface
        mn += 1  # Bump global face/pair counter.
        found[mn] && continue
        nufp += 1
        idxs = inrange(kdtree, data[mn], r, true)
        found[idxs] .= true
        push!(ufp2fp, idxs)
        ufpm[idxs] .= nufp
    end
    nufp == length(ufp2fp) || error("Miscount of nufp")

    return RWGData(bfe, bff, ebf, eci, ufpm, ufp2fp, zorymat, rhs, bfftstore, nufp)

end # function setup_rwg

zhatcross(t) = SV2([-t[2], t[1]])
zhatcross(t::SV2) = SV2(-t[2], t[1])

"""
    edge_current_unit_vector(ie::Integer, rwgdat::RWGData, metal::RWGSheet)::SV2

Evaluate a unit vector u = SV2([ux,uy]) in the positive reference direction for
the basis function associated with edge ie of the triangulated sheet.
"""
function edge_current_unit_vector(ie::Integer, rwgdat::RWGData, metal::RWGSheet)::SV2
    bf = rwgdat.ebf[ie]    #  Basis function index associated with edge ie
    bf == 0 && error("No basis function for edge $ie")
    # Begin by assuming the normal is parallel to \zhat \cross (\r_2 - \r_1):
    n1 = metal.e1[ie]   # Initial node of edge ie.
    n2 = metal.e2[ie]   # Terminal node of edge ie.
    ρ21 = metal.ρ[n2] - metal.ρ[n1]
    d = norm(ρ21)
    u = zhatcross(ρ21/d)
    # Now check that dot product of assumed unit vector with \vecrho^+ evaluated
    # at one of the edge nodes is positive:
    f = rwgdat.bff[1,bf]# Obtain the "plus" face adjacent to edge ie.
    # Find the free vertex of this face:
    nfree = 0 # Establish scope
    for i in 1:3
      nfree = metal.fv[i,f]
      nfree != n1  &&  nfree != n2 && break
    end
    ρ = metal.ρ[n1] - metal.ρ[nfree]
    if ρ ⋅ u < 0
      return -u
    end
    return u
end

"""
    get_ie_ξη(ie::Int, sheet::RWGSheet) -> (ξinit, ξterm, ηinit, ηterm)

Evaluate the ξ and η coordinates of the initial and terminal vertices
of a given edge. The position vector of a point is represented as
ξ*s1 + η*s2, where s1 and s2 are the direct lattice vectors.

## Arguments:

- `ie`:  Global edge index.
- `sheet`: The RWGSheet object containing the triangulation info.

## Return value:

A 4-tuple containing

- `ξinit`: The ξ value of the initial vertex of edge `ie`.
- `ξterm`: The ξ value of the terminal vertex of edge `ie`.
- `η_init`: The η value of the initial vertex of edge ie.
- `η_term`: The η value of the terminal vertex of edge ie.

"""
@inline function get_ie_ξη(ie::Int, sheet::RWGSheet)
    node_init = sheet.e1[ie]
    node_term = sheet.e2[ie]
    ρinit = sheet.ρ[node_init]
    ρterm = sheet.ρ[node_term]
    ξinit = (sheet.β₁ ⋅ ρinit) / (2π)
    ηinit = (sheet.β₂ ⋅ ρinit) / (2π)
    ξterm = (sheet.β₁ ⋅ ρterm) / (2π)
    ηterm = (sheet.β₂ ⋅ ρterm) / (2π)
    return (ξinit, ξterm, ηinit, ηterm)
end

"""
    rwgbfft!(ft, rwgdat::RWGData, sheet::RWGSheet, k::AbstractVector, ψ₁::Real, ψ₂::Real) -> ft::Vector

Compute the 2D fourier transform (FT) of the set of modified Rao-Wilton-Glisson basis functions
defined in `rwgdat` and `sheet`, evaluated at the transform variable `k`. `ψ₁` and `ψ₂` are the
unit cell incremental phase shifts in radians.

## Arguments:

- `ft`: A vector of length `nbf` whose elements are 2-vectors, where `nbf` is the number of
    basis functions.  The contents of `ft` will be overwritten with the basis function FTs.
- `rwgdat`: Contains the basis function definitions.
- `sheet`: The contains the triangulation info.
- `k`: A 2-vector containing the transform variable value at which the fourier transforms are to be
    evaluated. Units are (1/meter).
- `ψ₁`, `ψ₂`:  Unit cell incremental phase shifts in units of radians.

## Return value:

A length `nbf` complex vector containing the Fourier transforms of the basis functions.  Units are meters^2.
`nbf` is the number of basis functions.

## Reference
Kim McInturff and Peter S. Simon, "The Fourier transform of linearly varying functions with polygonal
support," IEEE Trans. Antennas Propagat., Vol. 39, no. 9, Sept. 1991, pp. 1441-1443.

"""
function rwgbfft!(ft, rwgdat::RWGData, sheet::RWGSheet, k::AbstractVector, ψ₁::Real, ψ₂::Real)
    next = SVector{3, Int}(2,3,1)
    nbf = size(rwgdat.bfe,2) # Number of basis functions
    length(ft) == nbf || error("length(ft) is $(length(ft)) instead of $nbf")
    ft .*= 0
    floquet_factor = OffsetArray(SVector{5, ComplexF64}(1, 1, cis(-ψ₁), 1, cis(-ψ₂)),   0:4)
    # Meanings:
    # floquet_factor[0] Edges not on unit cell boundary
    # floquet_factor[1] Edges at ξ=0 boundary.
    # floquet_factor[2] Edges at ξ=1 boundary.
    # floquet_factor[3] Edges at η=0 boundary.
    # floquet_factor[4] Edges at η=1 boundary.

    nface = size(sheet.fv,2) # Number of triangular faces
    kmagsq = k ⋅ k
    kmag = sqrt(kmagsq)
    one_meter = ustrip(Float64, sheet.units, 1u"m")

    lvec = zeros(SV2, 3)
    rc = zeros(SV2, 3)
    j0kl2 = MVector{3,Float64}(0.0, 0.0, 0.0)
    rtrm2 = zeros(SV2, 3)
    cphasv = MVector{3,ComplexF64}(0.0im, 0.0im, 0.0im)
    csum = MVector{2,ComplexF64}(0.0+0.0im, 0.0+0.0im)
    centroid = SVector{2,Float64}(0.0, 0.0)
    ft0 = SVector{2,ComplexF64}(0.0+0.0im, 0.0+0.0im)
    for iface in 1:nface
        r =  vtxcrd_m(iface, sheet, one_meter)
        lvec .= (r[next[i]] - r[i] for i in 1:3) # Edge vectors
        rc .= (r[i] + 0.5 * lvec[i] for i in 1:3)  # Edge centers
        if kmag * maximum(norm.(lvec)) < 1e-4 # small k
            centroid = mean(r)   #  Compute centroid coordinates.
            cphase = cis(k ⋅ centroid) #  Phase factor at centroid.
            for i in 1:3 #  Loop over three edges of the triangle
                ie = sheet.fe[i,iface] # Global index for edge opposite r(i).
                ib = rwgdat.ebf[ie]    # Global basis function index.
                ib == 0 && continue
                ft0 = 0.5 * (centroid - r[i])
                if rwgdat.bff[1,ib] == iface  # Plus triangle
                    ft[ib] += ft0 * (cphase * floquet_factor[rwgdat.eci[ie]])
                elseif rwgdat.bff[2,ib] == iface # Minus triangle
                    ft[ib] -= ft0 * (cphase * floquet_factor[rwgdat.eci[ie]])
                else
                    error("Impossible situation!")
                end
            end
        else # Case where k is not small
            kfact = 2k / kmagsq
            darea = zdotaxb(lvec[1], lvec[2]) # Twice the directed area
            denom = darea * kmagsq
            zdotlxk = [zdotaxb(l, k) for l in lvec]
            ctrm1 = [im * rci - kfact for rci in rc]
            for i in 1:3  # Loop over three edges
                dotkl2 = 0.5 * (k ⋅ lvec[i])
                j0kl2[i] = j₀(dotkl2)
                rtrm2[i] = 0.5 * zdotlxk[i] * j₁(dotkl2) * lvec[i]
                cphasv[i] = cis(k ⋅ rc[i])
            end
            #
            #  Loop over edges of triangle associated with basis functions.
            #
            for i in 1:3
                ie = sheet.fe[i,iface] # Global index for edge opposite r[i].
                ib = rwgdat.ebf[ie]    # Global basis function index.
                ib == 0 && continue # Skip if no basis func for this edge.
                cjr = im * r[i]
                csum .= complex(0.0,0.0)
                for n in 1:3  # Perform sum over n as shown in Equation (2-11):
                    ctrm3 = (zhatcross(lvec[n]) + zdotlxk[n]*(ctrm1[n]-cjr)) * j0kl2[n]
                    csum += cphasv[n] * (ctrm3 - rtrm2[n])
                end
                #csum = csum * norm(lvec(next(i))) # Needed for orig. defn. of RWG
                #                                          # basis funct.
                #  Add to sum total fourier transform with proper sign:
                if rwgdat.bff[1,ib] == iface  # Plus triangle.
                    ft[ib] += csum * (floquet_factor[rwgdat.eci[ie]] / denom)
                else
                    ft[ib] -= csum * (floquet_factor[rwgdat.eci[ie]] / denom)
                end
            end
        end
    end

    return ft
end # function

"""
    j₀(x)

Spherical Bessel function of the first kind of order 0 and argument x.
"""
j₀(x::Real) =     abs(x) < 1e-3 ? begin x² = x*x; 1 - x²/6*(1-x²/20) end : sin(x)/x
j₀(x::Complex) = abs2(x) < 1e-6 ? begin x² = x*x; 1 - x²/6*(1-x²/20) end : sin(x)/x

"""
    j₁(x)

Spherical Bessel function of the first kind of order 1 and argument x.
"""
j₁(x::Real) =     abs(x) < 1e-3 ? x * (1/3 - x*x/30) : begin s,c = sincos(x); (s/x - c) / x end
j₁(x::Complex) = abs2(x) < 1e-6 ? x * (1/3 - x*x/30) : begin s,c = sincos(x); (s/x - c) / x end

"""
    vtxcrd_m(iface::Int, sheet::Sheet, one_meter::Real)

Return the coordinates (in meters) of the triangle vertices for face iface.
"""
@inline function vtxcrd_m(iface, sheet, one_meter)
    vi = @view sheet.fv[:,iface] # Vertex indices
    sheet.ρ[vi] ./ one_meter
end

zdotaxb(a,b) = a[1] * b[2] - a[2] * b[1]

end # module


text:
"""
`preview(df::DataFrame, n::Integer; show_row_number=true)`
gives preview of first, middle and last columns of total n rows.
"""
function preview(df::DataFrame, n::Integer; show_row_number=true)
	nrow, ncol = size(df);
	if show_row_number
		rows = DataFrame("Row"=>rownumber.(eachrow(df)));
		df = hcat(rows,df);
	end

	hsplit = similar(df, 1);
	for col in names(df)
		hsplit[!,col] = ["..."];
	end

	nsub = Int64(floor(n/3));
	df_prev = [first(df, nsub); hsplit; middle(df, nsub-1); hsplit];
	hsofar = size(df_prev)[1];
	df_prev = [df_prev; last(df,n - hsofar)];
	# For `pretty_table()`'s option, see DataFrames.jl\src\abstractdataframe\show.jl
	DataFrames.pretty_table(df_prev; nosubheader=true, show_row_number=false);
	println("($nrow rows × $ncol columns)");
	nothing
end

"""
`middle(df::DataFrame, n::Integer)`
Similar to `first(df, n)` or `last(df, n)`, but get the dataframe of middle n rows.
"""
function middle(df::DataFrame, n::Integer)
	height = size(df)[1];
	mid0 = Int64(floor(height/2-n/2));
	mid1 = mid0 + n - 1;
	return df[mid0:mid1,:]
end


text:
module TestSummary

using Test
using TMLEEpistasis
using CSV
using DataFrames

include("helper_fns.jl")

sieve_filename(prefix) = string(prefix, "_sieve_variance.hdf5")
summary_filename(prefix) = string(prefix, "_summary.csv")

function build_sieve_file(prefix)
    sieve_file = jldopen(sieve_filename(prefix), "w")
    sieve_file["STDERRORS"] = [0.3, nothing]
    sieve_file["SOURCEFILE_REPORTID_PAIRS"] =
        ["$(prefix)_batch_1_Real.hdf5" => "2_1", "$(prefix)_batch_1_Real.hdf5" => "2_2"]
    close(sieve_file)
end

@testset "Test sieve_result_" begin
    @test all(x===missing for x in TMLEEpistasis.sieve_results_(0.3, nothing))
    @test all(x isa Float64 for x in TMLEEpistasis.sieve_results_(0.3, 0.1))
end

@testset "Test build_summary" begin
    grm_ids = TMLEEpistasis.GRMIDs("data/grm/test.grm.id")
    prefix = "rs12_rs45"
    summaryfilename = summary_filename(prefix)
    build_results_files(grm_ids, prefix)
    # With sieve is false
    parsed_args = Dict(
        "prefix" => prefix,
        "out" => summary_filename(prefix),
        "sieve" => false
    )
    build_summary(parsed_args)

    summary_no_sieve = CSV.read(summaryfilename, DataFrame)

    rm(summaryfilename)

    # With sieve is true
    build_sieve_file(prefix)
    parsed_args = Dict(
        "prefix" => prefix,
        "out" => summary_filename(prefix),
        "sieve" => true
    )

    build_summary(parsed_args)

    summary = CSV.read(summaryfilename, DataFrame)

    @test summary.PHENOTYPE == ["cancer", "cancer", "height", "height", "bmi", "bmi"]
    @test summary.PHENOTYPE_TYPE == ["BINARY", "BINARY", "CONTINUOUS", "CONTINUOUS", "CONTINUOUS", "CONTINUOUS"]
    @test summary.QUERYNAME == ["QUERY_1", "QUERY_2", "QUERY_1", "QUERY_2", "QUERY_1", "QUERY_2"]
    @test summary.FILENAME_ORIGIN == ["rs12_rs45_batch_1_Bool.hdf5",
                                      "rs12_rs45_batch_1_Bool.hdf5",
                                      "rs12_rs45_batch_1_Real.hdf5",
                                      "rs12_rs45_batch_1_Real.hdf5",
                                      "rs12_rs45_batch_1_Real.hdf5",
                                      "rs12_rs45_batch_1_Real.hdf5"]
    @test summary.REPORT_KEY == ["1_1", "1_2", "1_1", "1_2", "2_1", "2_2"]
    @test eltype(summary.INITIAL_ESTIMATE) == Float64
    @test eltype(summary.ESTIMATE) == Float64
    @test eltype(summary.PVAL) == Float64
    @test eltype(summary.LWB) == Float64
    @test eltype(summary.UPB) == Float64
    @test summary.SIEVE_STDERR[5] == 0.3
    @test summary.SIEVE_PVAL[5] isa Float64
    @test summary.SIEVE_LWB[5] isa Float64
    @test summary.SIEVE_UPB[5] isa Float64

    clean_hdf5estimates_files(prefix)
    rm(sieve_filename(prefix))
    rm(summaryfilename)

    # Check the summary with and without sieve match on the same columns
    sievecolnames = names(TMLEEpistasis.init_sieve("toto"))
    @test summary_no_sieve == summary[!, Not(sievecolnames)]

end

end

true


text:
using Gadfly

function plot_matrix( tiles_matrix, path )
     light_blue = color("blue")
     dark_blue  = color("red")
     sc = Scale.ContinuousColorScale(Scale.lab_gradient( dark_blue,
                                                         light_blue ) )
     num_rows = size(tiles_matrix,1)
     Gadfly.draw(PNG( path, 10inch, 10inch), spy2( tiles_matrix,  sc, Scale.y_continuous(maxvalue=num_rows) ) )
     # or call the R code, because the time this takes can be excessive??
end

function plot_matrix2(tiles_matrix,path; colour_scale =Scale.discrete_color_manual("black","red","green") )
    #border1 = color("blue")
    #border2 = color("green")
    num_rows = size(tiles_matrix,1)
    #sc = Scale.ContinuousColorScale(Scale.lab_gradient( border1,
    #                                                    border2 ) )
    #dc =Scale.discrete_color_manual("black","red","green")
    is, js, values = findnz(tiles_matrix)
    df = DataFrame(i=is, j=js, value=values)
    T = Theme(panel_fill=color("black"))
    p = Gadfly.plot(df, x="j", y="i", color="value",
        T,
        Coord.cartesian(yflip=true),
        #Scale.continuous_color,
        colour_scale,
        Scale.x_continuous,
        Scale.y_continuous(maxvalue=num_rows),
        Geom.rectbin,
        Stat.identity)

    Gadfly.draw(PNG( path, 10inch, 10inch), p )
end

# So you want to get the number of genes in each cluster and make a bar plot of those
# But the bar plot should be coloured by the average methylation level
# So we have to get a value for the average number of tiles methylated

# A. Get the number in each cluster
# B. Get the average methylation for each cluster

function plot_cluster_bar_chart(tiles_matrix, cluster_assignments,path,id=""; gadfly_theme = Theme(major_label_color = color("black"), minor_label_color=color("black")) )

    # GET NUMBERS OF REGIONS IN EACH CLUSTER
    # Make dictionary with keys as cluster_id
    # and value how many regions in each cluster
    cluster_counts = Dict{Int64,Int64}()
    for c in cluster_assignments
        if haskey( cluster_counts,c )
            cluster_counts[c] +=1
        else
            cluster_counts[c]=1
        end
     end

    # GET SUM OF TILES VALUES FOR EACH REGION
    i=1
    cluster_value_sum =Dict{Int64,Float64}()
    tiles_matrix_transpose = tiles_matrix
    for c in cluster_assignments
         sum_of_tiles = sum(tiles_matrix_transpose[i,:])
         if haskey(cluster_value_sum,c)
             cluster_value_sum[c] +=sum_of_tiles
         else
             cluster_value_sum[c] = sum_of_tiles
         end
         i +=1
     end

     # CALCULATE AVERAGE VALUE TILES FROM EACH REGION
     ave_cluster_value = Float64[]
     cluster_num = Int64[]
     for k in sort(collect(keys(cluster_value_sum)))
        println(cluster_value_sum[k])
        println(cluster_counts[k])
        ave_value_tiles = cluster_value_sum[k] / cluster_counts[k]
        push!(ave_cluster_value, ave_value_tiles)
        push!(cluster_num,cluster_counts[k])
     end

     # PLOT
     #sort_by_clusterN = sortperm(cluster_num,rev=true)
     df = DataFrame( average_cluster_value = ave_cluster_value, cluster_num = cluster_num)
     sort!(df, cols = [:cluster_num],rev=true)
     nrows = nrow(df)
     if nrows > 15
         nrows = 15
     end
     df = df[1:nrows,:]
     gp1 = plot(df, x=1:nrows, y="cluster_num", color="average_cluster_value",
         gadfly_theme,Geom.bar,
         Guide.title("Genes in Cluster $id"),
         Guide.xlabel("Cluster Number"),
         Guide.ylabel("Num of Regions"),
         Guide.colorkey("Avg Tiles Value"),
         Scale.y_continuous(format=:plain))
     Gadfly.draw(PNG(path,20cm, 15cm), gp1)
     return gp1
end


text:
#
#=
Examples for developing the JAC tools
=#
println("A)  Tests & development of the electronic structure part.")
println("---------------------------------------------------------")
println("Aa) Test of the radial integration in CI computations (Randolf).")
println("Ab) Read in a full multiplet representation from files compatible to Ratip2012, Grasp2018, ....")
println("Ac) Test of several radial e-e potentials.")
println("Ad) Tests towards a SCF field: B-spline primitives and one-particle spectra in a local potential.")
println("Ae) Test of the CI part for an internally generated neon multiplet without Breit interaction.")
println("Af) Test of the frequency-independent Breit interaction for an internally generated neon multiplet.")
println("Ag) Test of the QED model corrections to the level structure of atoms and ions.")
println("Ah) Generate and normalize continuum orbitals in a local potential.")
println("")
println("B)  Tests & development of atomic amplitudes")
println("--------------------------------------------")
println("Ba) Tests of the dipole and em multipole amplitudes.")
println("Bb) Tests of the momentum-transfer amplitudes.")
println("Bc) Tests of the parity non-conservation, Schiff moment and anapole moment amplitudes.")
println("")
println("C)   Tests & development of atomic properties")
println("---------------------------------------------")
println("Ca)  Test of the Einstein module with ASF from an internally generated multiplet.")
println("Cb)  Test of the Hfs module for HFS A,B parameters and hyperfine representation with ASF from an internally generated multiplet.")
println("Cc)  Test of the IsotopeShift module with ASF from an internally generated multiplet.")
println("Cd)  Test of the FormFactor module with ASF from an internally generated multiplet.")
println("Ce)  Test of the PlasmaShift module with ASF from an internally generated multiplet.")
println("Cfg) Test of the LandeZeeman module with ASF from an internally generated multiplet.")
println("Ch)  Test of the AlphaVariation module with ASF from an internally generated multiplet.")
println("Ci)  Test of the DecayYields module with ASF from an internally generated multiplet.")
println("Cj)  Test of the GreenFunction module with ASF from an internally generated multiplet.")
println("Ck)  Test of the MultipolePolarizibility module with ASF from an internally generated multiplet.")
println("")
println("D)   Tests & development of atomic processes")
println("--------------------------------------------")
println("Da)  Test of the PhotoEmission module with ASF from an internally generated initial- and final-state multiplet.")
println("Db)  Test of the PhotoExcitation module with ASF from an internally generated initial- and final-state multiplet.")
println("Dc)  Test of the Photoionization module with ASF from an internally generated initial- and final-state multiplet.")
println("Dcp) Test of the Photoionization and PlasmaShift module with ASF from an internally generated initial- and final-state multiplet.")
println("Dd)  Test of the PhotoRecombination module with ASF from an internally generated initial- and final-state multiplet.")
println("De)  Test of the Auger module with ASF from an internally generated initial- and final-state multiplet.")
println("Dep) Test of the Auger and PlasmaShift modules with ASF from an internally generated initial- and final-state multiplet.")
println("Df)  Test of the Dielectronic module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
println("Dg)  Test of the PhotoExcitationFluores module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
println("Dh)  Test of the PhotoExcitationAutoIon module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
println("Di)  Test of the RayleighCompton module with ASF from an internally generated initial- and final-state multiplets.")
println("Dj)  Test of the MultiPhotonDeExcitation module with ASF from an internally generated initial- and final-state multiplet.")
println("Dk)  Test of the CoulombExcitation module with ASF from an internally generated initial- and final-state multiplet.")
println("Dl)  Test of the PhotoIonizationFluores module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
println("Dm)  Test of the PhotoIonizationAutoion module with ASF from an internally generated initial-, intermediate and final-state multiplets.")
println("Dn)  Test of the ImpactExcitation module with ASF from an internally generated initial- and final-state multiplet.")
println("Dp)  Test of the ImpactExcitationAutoion module with ASF from an internally generated initial-, intermediate and final-state multiplet.")
println("Dq)  Test of the RadiativeAuger module with ASF from an internally generated initial- and final-state multiplet.")
println("Dr)  Test of the MultiPhotonIonization module with ASF from an internally generated initial- and final-state multiplet.")
println("Ds)  Test of the MultiPhotonDoubleIon module with ASF from an internally generated initial- and final-state multiplet.")
println("Dt)  Test of the InternalConversion module with ASF from an internally generated initial- and final-state multiplet.")

println("The next tests need first to be adapted.")

println("")
println("Dx9) Test of the PairAnnihilation1Photon.jl module with ASF from an internally generated initial- and final-state multiplet.")
println("Dy5) Test of the CoulombIonization module with ASF from an internally generated initial- and final-state multiplet.")
println("")
println("E)  Tests & development of atomic cascades")
println("------------------------------------------")
println("E1) Three-step cascade after 1s-3p photo-excitation of Si^+: Configuration model.")
println("E2) Simulation of the four-step cascade after 1s-3p photo-excitation of Si^+: Configuration model.")
println("")
println("F)  Tests & development of semiempirical computations")
println("-----------------------------------------------------")
println("F1) Tests of the semi-empirical computation of ADK rates.")


text:
create(c::Client, ::Type{ApplicationCommand}, name::AbstractString, description::AbstractString; kwargs...) = create_application_command(c; name=name, description=description, kwargs...)
create(c::Client, ::Type{ApplicationCommand}, name::AbstractString, description::AbstractString, g::Snowflake; kwargs...) = create_application_command(c, g; name=name, description=description, kwargs...)
create(c::Client, ::Type{Message}, interaction::Interaction; kwargs...) = create_followup_message(c, interaction.id, interaction.token; compkwfix(; kwargs...)...)
create(c::Client, ::Type{Vector{ApplicationCommand}}, cmds::Vector{ApplicationCommand}) = bulk_overwrite_application_commands(c, cmds)
create(c::Client, ::Type{Vector{ApplicationCommand}}, g::Snowflake, cmds::Vector{ApplicationCommand}) = bulk_overwrite_application_commands(c, g, cmds)

retrieve(c::Client, ::Type{Vector{ApplicationCommand}}) = get_application_commands(c)
retrieve(c::Client, ::Type{Vector{ApplicationCommand}}, g::Snowflake) = get_application_commands(c, g)

update(c::Client, ctx::Context, m::Message; kwargs...) = edit_interaction(c, ctx.interaction.token, m.id; compkwfix(; kwargs...)...)
update(c::Client, ctx::Context; kwargs...) = hasproperty(ctx, :interaction) ?  update_message_int(c, ctx.interaction.id, ctx.interaction.token; compkwfix(; kwargs...)...) : update(c, ctx.message; compkwfix(; kwargs...)...)


text:
using DelimitedFiles
using Formatting

include("fibonacci.jl")
include("matrix_multiplication.jl")
include("quicksort.jl")

function chop_suffix_from_path(path, suffix)
    if endswith(path, suffix)
        return path[1:(end-length(suffix))]
    end
    return path
end

function run_test(suite, ntests)
    if suite == 1
        test_fibonacci(ntests)
    elseif suite == 2
        test_quicksort(ntests)
    elseif suite == 3
        test_matrix_multiplication(ntests)
    elseif suite == 4
        test_all(ntests)
    end
end

function test_fibonacci(ntests)
    fibarg = 20
    base_path = chop_suffix_from_path(pwd(), "/benchmarks/Julia")
    path = string(base_path, "/results/fibonacci/julia_fibonacci_benchmark.txt")

    @assert fib(20) == 6765

    open(path, "w") do fd
        for _ in 1:ntests
            time = @elapsed fib(fibarg)
            write(fd, string(fmt(".6f", time), "\n"))
        end
    end
end

function test_quicksort(ntests)
    arraysize = 1000000
    base_path = chop_suffix_from_path(pwd(), "/benchmarks/Julia")
    read_path = string(base_path, "/utils/data/integers.txt")

    open(read_path, "r") do fdread
        test_array = [parse(Int, i) for i in readlines(read_path)]
        quicksort!(test_array, 1, arraysize)
        @assert issorted(test_array)

        write_path = string(base_path, "/results/quicksort/julia_quicksort_benchmark.txt")
        open(write_path, "w") do fdwrite
            for _ in 1:ntests
                seekstart(fdread)
                array = [parse(Int, i) for i in readlines(read_path)]
                time = @elapsed quicksort!(array, 1, arraysize)
                write(fdwrite, string(fmt(".6f", time), "\n"))
            end
        end
    end
end

function test_matrix_multiplication(ntests)
    matrix_size = 200
    base_path = chop_suffix_from_path(pwd(), "/benchmarks/Julia")
    write_path = string(base_path, "/results/matrix_multiplication/julia_matrix_multiplication_benchmark.txt")
    read_path = string(base_path, "/utils/data/matmul_assert.txt")

    A = B = [col-1 for row in 1:matrix_size, col in 1:matrix_size]
    R = zeros(Int, matrix_size, matrix_size)
    matmul!(A, B, R)
    C = readdlm(read_path, ' ', BigInt, '\n')
    @assert R == C

    open(write_path, "w") do fd
        for _ in 1:ntests
            R = zeros(Int, matrix_size, matrix_size)
            time = @elapsed matmul!(A, B, R)

            write(fd, string(fmt(".6f", time), "\n"))
        end
    end
end

function test_all(ntests)
    test_fibonacci(ntests)
    test_quicksort(ntests)
    test_matrix_multiplication(ntests)
end

max_stars_repo_name: mjaglarz/BA-thesis

max_issues_repo_name: mjaglarz/BA-thesis

max_forks_repo_name: mjaglarz/BA-thesis


text:
module Meas

struct Accumulator
    count::Dict{String,UInt64}
    data::Dict{String,Any}
end

"""
Constructor
"""
function Accumulator()
    Accumulator(Dict{String,UInt64}(), Dict{String,Any}())
end

"""
Add a sample
"""
function add!(acc::Accumulator, name::String, data)
    if haskey(acc.count, name)
        acc.count[name] += 1
        acc.data[name] += data
    else
        acc.count[name] = 1
        acc.data[name] = copy(data)
    end
end

"""
Compuate mean
"""
function mean(acc::Accumulator, name::String)
    return acc.data[name]/acc.count[name]
end

end


text:
#=
  From BLOG examples/mixture-of-gaussian-full.blog

  Summary Statistics
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

           p    0.3237    0.3005     0.0030    0.0033   8144.0941    0.9999      464.8986
           a    0.1618    0.5440     0.0054    0.0063   8092.3273    0.9999      461.9436
           b    0.2979    0.4666     0.0047    0.0052   7836.0821    0.9999      447.3160
        z[1]    0.3249    0.4684     0.0047    0.0048   8429.8487    1.0000      481.2107
        z[2]    0.3224    0.4674     0.0047    0.0045   8398.4099    0.9999      479.4160
        z[3]    0.3204    0.4667     0.0047    0.0054   8177.2874    1.0000      466.7934
        z[4]    0.3247    0.4683     0.0047    0.0050   8229.2433    0.9999      469.7593
      min_ab   -0.0771    0.4745     0.0047    0.0058   7699.9809    0.9999      439.5468

  Cf ~/webppl/mixture_of_gaussian2.wppl

=#

using Turing, StatsPlots
include("jl_utils.jl")

@model function mixture_of_gaussian2(x=[0.2,1.0,0.5,0.6])
    n = length(x)
    p ~ Beta(0.5, 1)

    a ~ Uniform(-1, 1)
    b ~ Uniform(-1, 1)
    z ~ filldist(flip(p), n)
    for i in 1:n
        if z[i] == true
            x[i] ~ Normal(a, 1.0)
        else
            x[i] ~ Normal(b, 1.0)
        end
    end

   min_ab ~ Dirac(min(a,b))

end

x=[0.2,1.0,0.5,0.6]
model = mixture_of_gaussian2(x)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
mutable struct Comm
    ptr::Ptr{Cvoid}
end

Base.unsafe_convert(::Type{Ptr{Cvoid}}, comm::Comm) = comm.ptr

function Comm(nranks::Int, commid, rank::Int)
    ref = Ref{Ptr{Cvoid}}()
    @nccl :ncclCommInitRank (Ptr{Ptr{Cvoid}},Cint,Ptr{UInt8},Cint) ref nranks commid rank
    Comm(ref[])
end

function Comm(ndev::Int, devlist::Vector{Int})
    ref = Ref{Ptr{Cvoid}}()
    @nccl :ncclCommInitAll (Ptr{Ptr{Cvoid}},Cint,Ptr{Cint}) ref ndev devlist
    Comm(ref[])
end

function destroy(comm::Comm)
    @nccl :ncclCommDestroy (Ptr{Cvoid},) comm
end

function get_unique_id()
    ref = Array{UInt8}(NCCL_UNIQUE_ID_BYTES)
    @nccl :ncclGetUniqueId (Ptr{UInt8},) ref
    unsafe_string(ref[])
end

function cudevice(comm::Comm)
    ref = Ref{Cint}()
    @nccl :ncclCommCuDevice (Ptr{Cvoid},Ptr{Cint}) comm ref
    Int(ref[])
end

function userrank(comm::Comm)
    ref = Ref{Cint}()
    @nccl :ncclCommUserRank (Ptr{Cvoid},Ptr{Cint}) comm ref
    Int(ref[])
end


text:
export iauAtoiq
"""
Quick observed place to CIRS, given the star-independent astrometry
parameters.

Use of this function is appropriate when efficiency is important and
where many star positions are all to be transformed for one date.
The star-independent astrometry parameters can be obtained by
calling iauApio[13] or iauApco[13].

Status:  support function.

Given:
   type   char[]     type of coordinates: "R", "H" or "A" (Note 1)
   ob1    double     observed Az, HA or RA (radians; Az is N=0,E=90)
   ob2    double     observed ZD or Dec (radians)
   astrom iauASTROM* star-independent astrometry parameters:
      pmt    double       PM time interval (SSB, Julian years)
      eb     double[3]    SSB to observer (vector, au)
      eh     double[3]    Sun to observer (unit vector)
      em     double       distance from Sun to observer (au)
      v      double[3]    barycentric observer velocity (vector, c)
      bm1    double       sqrt(1-|v|^2): reciprocal of Lorenz factor
      bpn    double[3][3] bias-precession-nutation matrix
      along  double       longitude + s' (radians)
      xpl    double       polar motion xp wrt local meridian (radians)
      ypl    double       polar motion yp wrt local meridian (radians)
      sphi   double       sine of geodetic latitude
      cphi   double       cosine of geodetic latitude
      diurab double       magnitude of diurnal aberration vector
      eral   double       "local" Earth rotation angle (radians)
      refa   double       refraction constant A (radians)
      refb   double       refraction constant B (radians)

Returned:
   ri     double*    CIRS right ascension (CIO-based, radians)
   di     double*    CIRS declination (radians)

Notes:

   1. "Observed" Az,El means the position that would be seen by a
   perfect geodetically aligned theodolite.  This is related to
   the observed HA,Dec via the standard rotation, using the geodetic
   latitude (corrected for polar motion), while the observed HA and
   RA are related simply through the Earth rotation angle and the
   site longitude.  "Observed" RA,Dec or HA,Dec thus means the
   position that would be seen by a perfect equatorial with its
   polar axis aligned to the Earth's axis of rotation.  By removing
   from the observed place the effects of atmospheric refraction and
   diurnal aberration, the CIRS RA,Dec is obtained.

   2. Only the first character of the type argument is significant.
   "R" or "r" indicates that ob1 and ob2 are the observed right
   ascension and declination;  "H" or "h" indicates that they are
   hour angle (west +ve) and declination;  anything else ("A" or
   "a" is recommended) indicates that ob1 and ob2 are azimuth (north
   zero, east 90 deg) and zenith distance.  (Zenith distance is used
   rather than altitude in order to reflect the fact that no
   allowance is made for depression of the horizon.)

   3. The accuracy of the result is limited by the corrections for
   refraction, which use a simple A*tan(z) + B*tan^3(z) model.
   Providing the meteorological parameters are known accurately and
   there are no gross local effects, the predicted observed
   coordinates should be within 0.05 arcsec (optical) or 1 arcsec
   (radio) for a zenith distance of less than 70 degrees, better
   than 30 arcsec (optical or radio) at 85 degrees and better than
   20 arcmin (optical) or 30 arcmin (radio) at the horizon.

   Without refraction, the complementary functions iauAtioq and
   iauAtoiq are self-consistent to better than 1 microarcsecond all
   over the celestial sphere.  With refraction included, consistency
   falls off at high zenith distances, but is still better than
   0.05 arcsec at 85 degrees.

   4. It is advisable to take great care with units, as even unlikely
   values of the input parameters are accepted and processed in
   accordance with the models used.

Called:
   iauS2c       spherical coordinates to unit vector
   iauC2s       p-vector to spherical
   iauAnp       normalize angle into range 0 to 2pi

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# void iauAtoiq(const char *type,
#               double ob1, double ob2, iauASTROM *astrom,
#               double *ri, double *di)

function iauAtoiq(type_::Char, ob1::Real, ob2::Real, astrom::iauASTROM)

   # Allocate return values
   ref_type_  = Ref{UInt8}(convert(UInt8, type_))
   ref_astrom = Ref{iauASTROM}(astrom)
   ref_ri     = Ref{Float64}(0.0)
   ref_di     = Ref{Float64}(0.0)

   ccall((:iauAtoiq, libsofa_c), Cvoid,
       (Ref{UInt8}, Cdouble, Cdouble, Ref{iauASTROM},
       Ref{Cdouble}, Ref{Cdouble}),
       ref_type_, convert(Float64, ob1), convert(Float64, ob2),
       ref_astrom, ref_ri, ref_di)

   return ref_ri[], ref_di[]
end


text:
@testset "to_unsigned" begin

    @test BTCParser.to_unsigned((0x01, 0x02)) == reinterpret(UInt16, UInt8[0x01, 0x02])[1]

    @test BTCParser.to_unsigned((0x01, 0x02, 0x03, 0x04)) ==
        reinterpret(UInt32, UInt8[0x01, 0x02, 0x03, 0x04])[1]

    @test BTCParser.to_unsigned((0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08)) ==
        reinterpret(UInt64, UInt8[0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08])[1]

    @test BTCParser.to_unsigned((
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08
    )) == reinterpret(UInt128, UInt8[
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08
    ])[1]

    @test BTCParser.to_unsigned((
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
    )) == reinterpret(BTCParser.UInt256, UInt8[
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
    ])[1]
end

@testset "to_byte_tuple" begin

    @test (0x01, 0x02) == BTCParser.to_byte_tuple(reinterpret(UInt16, UInt8[0x01, 0x02])[1])

    @test (0x01, 0x02, 0x03, 0x04) ==
        BTCParser.to_byte_tuple(reinterpret(UInt32, UInt8[0x01, 0x02, 0x03, 0x04])[1])

    @test (0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08) ==
        BTCParser.to_byte_tuple(reinterpret(UInt64, UInt8[0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08])[1])

    @test (
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08
    ) == BTCParser.to_byte_tuple(reinterpret(UInt128, UInt8[
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08
    ])[1])

    @test (
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
    ) == BTCParser.to_byte_tuple(reinterpret(UInt256, UInt8[
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
        0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08,
    ])[1])
end


text:
function nordsieck_adjust!(integrator, cache::T) where T
  @unpack nextorder, order = cache
  if nextorder != order
    # TODO: optimize?
    nordsieck_adjust_order!(cache, nextorder-order)
    cache.order = cache.nextorder
    cache.L = cache.order + 1
    cache.n_wait = cache.L
  end
  nordsieck_rescale!(cache)
  return nothing
end

# TODO: SUNDIALS NLsolve handling

function nordsieck_finalize!(integrator, cache::T) where T
  isconst = T <: OrdinaryDiffEqConstantCache
  @unpack order, dts = cache
  update_nordsieck_vector!(cache)
  cache.n_wait -= 1
  if is_nordsieck_change_order(cache, 1) && cache.order != 12
    if isconst
      cache.z[end] = cache.Δ
    else
      @. cache.z[end] = cache.Δ
    end
    cache.prev_𝒟 = cache.c_𝒟
  end
end

function nordsieck_prepare_next!(integrator, cache::T) where T
  isconst = T <: OrdinaryDiffEqConstantCache
  @unpack maxη, order, L = cache
  # TODO: further clean up
  @unpack bias1, bias2, bias3, addon = integrator.alg
  if integrator.EEst > one(integrator.EEst)
    nordsieck_rewind!(cache)
    cache.n_wait = max(2, cache.n_wait)
    cache.nextorder = order
    cache.η = inv( (bias2*integrator.EEst)^inv(L) + addon )
    return nothing
  end
  cache.ηq = inv( (bias2*integrator.EEst)^inv(L) + addon )
  stepsize_η!(integrator, cache, cache.order)
  if !is_nordsieck_change_order(cache)
    cache.η = cache.ηq
    cache.nextorder = order
    setη!(integrator, cache)
    return nothing
  end
  # On an order change (cache.n_wait == 0), we are going to compute the η for
  # order q+1 and q-1, where η = dt_next/dt
  cache.n_wait = 2
  stepsize_η₊₁!(integrator, cache, order)
  stepsize_η₋₁!(integrator, cache, order)
  chooseη!(integrator, cache)
  setη!(integrator, cache)
  # TODO: Maybe not here
  if isconst
    cache.Δ = cache.c_LTE * cache.Δ
  else
    @. cache.Δ = cache.c_LTE * cache.Δ
  end
  return nothing
end

##############################################################
# Lower level functions
##############################################################

# This function computes the integral, from -1 to 0, of a polynomial
# `P(x)` from the coefficients of `P` with an offset `k`.
function ∫₋₁⁰dx(a, deg, k)
  @inbounds begin
    int = zero(eltype(a))
    sign = 1
    for i in 0:deg
      int += flipsign(a[i+1]/(i+k), sign)
      sign = -sign
    end
    return int
  end
end

# `l` is the coefficients of the polynomial `Λ` that satisfies conditions
# Λ(0) = 1, Λ(-1) = 0, and Λ̇(-ξᵢ) = 0, where ξᵢ = (tₙ-tₙ₋₁)/dt.
# It is described in the paper "A Polyalgorithm for the Numerical Solution
# of Ordinary Differential Equations" by G. D. Byrne and A. C. Hindmarsh in
# the page 86.
# https://dl.acm.org/citation.cfm?id=355636

# More implementation details are in the
# https://github.com/JuliaDiffEq/DiffEqDevMaterials repository
function calc_coeff!(cache::T) where T
  isvode = ( T <: JVODECache || T <: JVODEConstantCache )
  @inbounds begin
    isconst = T <: OrdinaryDiffEqConstantCache
    isvarorder = is_nordsieck_change_order(cache, 1)
    @unpack m, l, dts, order = cache
    dtsum = dt = dts[1]
    if order == 1
      l[1] = l[2] = cache.c_LTE₋₁ = cache.c_𝒟 = 1
      cache.c_LTE = 1//2
      cache.c_LTE₊₁ = 1//12
      cache.c_conv = 1//10 / cache.c_LTE
      return nothing
    end
    m[1] = 1
    for i in 2:order+1
      m[i] = 0
    end
    # initialize ξ_inv
    ξ_inv = dt / dtsum
    # compute coefficients from the Newton polynomial
    # check the `JuliaDiffEq/DiffEqDevMaterials` repository for more details
    for j in 1:order-1
      if isvarorder && j == order-1
        M₋₁ = ∫₋₁⁰dx(m, order-2, 2)
        # It is the same with `tq[1]` in SUNDIALS cvode.c
        cache.c_LTE₋₁ = order * M₋₁ / m[order-1]
      end
      ξ_inv = dt / dtsum
      for i in j:-1:1
        m[i+1] = muladd(m[i], ξ_inv, m[i+1])
      end
      dtsum += dts[j+1]
    end
    ξ_inv = dt / dtsum

    M0 = ∫₋₁⁰dx(m, order-1, 1)
    M1 = ∫₋₁⁰dx(m, order-1, 2)
    M0_inv = inv(M0)
    l[1] = 1
    for i in 1:order
      l[i+1] = M0_inv * m[i] / i
    end
    # TODO: simplify LTE calculation
    # This is the error estimation coefficient for the current order `q`
    # ||Δ||⋅c_LTE yields the difference between a `q` degree interpolating
    # polynomial and a `q+1` degree interpolating polynomial at time `t`.
    # It is the same with `tq[2]` in SUNDIALS cvode.c
    cache.c_LTE = M1 * M0_inv * ξ_inv
    # It is the same with `tq[5]` in SUNDIALS cvode.c
    isvode && (cache.c_𝒟 = inv(ξ_inv) / l[order+1])
    if isvarorder
      for i in order:-1:1
        m[i+1] = muladd(ξ_inv, m[i], m[i+1])
      end
      M2 = ∫₋₁⁰dx(m, order, 2)
      # It is the same with `tq[3]` in SUNDIALS cvode.c
      cache.c_LTE₊₁ = M2 * M0_inv / (order+1)
    end # endif isvarorder
    # It is the same with `tq[4]` in SUNDIALS cvode.c
    cache.c_conv = 1//10 / cache.c_LTE
    return nothing
  end # end @inbounds
end

# Apply the Pascal linear operator
function perform_predict!(cache::T, rewind=false) where T
  @inbounds begin
    isconst = T <: OrdinaryDiffEqConstantCache
    @unpack z, order = cache
    # This can be parallelized
    if !rewind
      if isconst
        for i in 1:order, j in order:-1:i
          z[j] = z[j] + z[j+1]
        end
      else
        for i in 1:order, j in order:-1:i
          @. z[j] = z[j] + z[j+1]
        end
      end # endif const cache
    else
      if isconst
        for i in 1:order, j in order:-1:i
          z[j] = z[j] - z[j+1]
        end
      else
        for i in 1:order, j in order:-1:i
          @. z[j] = z[j] - z[j+1]
        end
      end # endif const cache
    end # endif !rewind
  end # end @inbounds
end

# Apply corrections on the Nordsieck vector
function update_nordsieck_vector!(cache::T) where T
  isvode = ( T <: JVODECache || T <: JVODEConstantCache )
  @inbounds begin
    isconst = T <: OrdinaryDiffEqConstantCache
    @unpack z,Δ,l,order = cache
    if isconst
      for i in 1:order+1
        z[i] = muladd.(l[i], Δ, z[i])
      end
    else
      for i in 1:order+1
        @. z[i] = muladd(l[i], Δ, z[i])
      end
    end # endif not const cache
  end # end @inbounds
end

function nlsolve_functional!(integrator, cache::T) where T
  @unpack f, dt, t, p = integrator
  isconstcache = T <: OrdinaryDiffEqConstantCache
  @unpack z, l, c_conv, Δ = cache
  if isconstcache
    ratetmp = integrator.f(z[1], p, dt+t)
  else
    @unpack ratetmp = cache
    integrator.f(ratetmp, z[1], p, dt+t)
  end
  integrator.destats.nf += 1
  max_iter = 3
  div_rate = 2
  # Zero out the difference vector
  isconstcache ? ( cache.Δ = zero(cache.Δ) ) : ( Δ .= zero(eltype(Δ)) )
  # `k` is a counter for convergence test
  k = 0
  # `conv_rate` is used in convergence rate estimation
  conv_rate = 1.
  # initialize `δ_prev`
  δ_prev = 0
  # Start the functional iteration & store the difference into `Δ`
  for k in 1:max_iter
    if isconstcache
      ratetmp = inv(l[2])*muladd.(dt, ratetmp, -z[2])
      integrator.u = ratetmp + z[1]
      cache.Δ = ratetmp - cache.Δ
    else
      @. integrator.u = -z[2]
      @. ratetmp = inv(l[2])*muladd(dt, ratetmp, integrator.u)
      @. integrator.u = ratetmp + z[1]
      @. cache.Δ = ratetmp - cache.Δ
    end
    # @show norm(dt*ratetmp - ( z[2] + (integrator.u - z[1])*l[2] ))
    # @show norm(cache.Δ - (integrator.u - z[1]))
    # It only makes sense to calculate convergence rate in the second iteration
    δ = integrator.opts.internalnorm(cache.Δ,t)
    isconstcache ? ( cache.Δ = copy(ratetmp) ) : copyto!(cache.Δ, ratetmp)
    if k >= 1
      conv_rate = max(1//10*conv_rate, δ/δ_prev)
    end
    test_rate = δ * min(one(conv_rate), conv_rate) / c_conv
    if test_rate <= one(test_rate)
      return true
    end
    # Divergence criteria
    if ( (k == max_iter) || (k >= 2 && δ > div_rate * δ_prev) )
      return false
    end
    δ_prev = δ
    integrator.destats.nf += 1
    isconstcache ? (ratetmp = integrator.f(integrator.u, p, dt+t)) :
                    integrator.f(ratetmp, integrator.u, p, dt+t)
  end
end

function nordsieck_rescale!(cache::T, rewind=false) where T
  isconstcache = T <: OrdinaryDiffEqConstantCache
  @unpack z, dts, order = cache
  eta = rewind ? dts[2]/dts[1] : dts[1]/dts[2]
  factor = eta
  for i in 2:order+1
    if isconstcache
      z[i] = z[i]*factor
    else
      rmul!(z[i], factor)
    end
    factor *= eta
  end
  return nothing
end

function nordsieck_rewind!(cache)
  perform_predict!(cache, true)
  nordsieck_rescale!(cache, true)
end

function is_nordsieck_change_order(cache::T, n=0) where T
  isconstcache = T <: OrdinaryDiffEqConstantCache
  isvode = ( T <: JVODECache || T <: JVODEConstantCache )
  isvode || return false
  cache.n_wait == 0+n
end

function nordsieck_decrement_wait!(cache::T) where T
  isvode = ( T <: JVODECache || T <: JVODEConstantCache )
  isvode || return nothing
  isconstcache = T <: OrdinaryDiffEqConstantCache
  cache.n_wait = max(0, cache.n_wait-1)
  return nothing
end

function nordsieck_adjust_order!(cache::T, dorder) where T
  isconstcache = T <: OrdinaryDiffEqConstantCache
  @unpack order, dts = cache
  # WIP: uncomment when finished
  #@inbound begin
  begin
    # Adams order increase
    if dorder == 1
      if isconstcache
        cache.z[order+2] = zero(cache.z[order+2])
      else
        cache.z[order+2] .= 0
      end
    else
      # Adams order decrease
      # One needs to rescale the Nordsieck vector on an order decrease
      cache.l .= 0
      cache.l[2] = 1
      dt = dts[1]
      hsum = zero(eltype(cache.dts))
      for j in 1:order-2
        hsum += cache.dts[j+1]
        # TODO: `hscale`?
        ξ = hsum / dt
        for i in j+1:-1:1
          cache.l[i+1] = cache.l[i+1] * ξ + cache.l[i]
        end # for i
      end # for j

      for j in 2:order-1
        cache.l[j+1] = order * cache.l[j] / j
      end
      for j in 3:order
        if isconstcache
          cache.z[j] = muladd.(-cache.l[j], cache.z[order+1], cache.z[j])
        else
          @. cache.z[j] = muladd(-cache.l[j], cache.z[order+1], cache.z[j])
        end
      end # for j
    end # else
  end # @inbound
end

# `η` is `dtₙ₊₁/dtₙ`
function setη!(integrator, cache::T) where T
  if cache.η < integrator.opts.qsteady_max
    cache.η = 1
  else
    # TODO: Not the same with SUNDIALS
    ( integrator.iter == 1 || integrator.u_modified ) && ( cache.η = min(1e5, cache.η); return nothing )
    cache.η = min(integrator.opts.qmax, max(integrator.opts.qmin, cache.η))
  end
  return nothing
end

function chooseη!(integrator, cache::T) where T
  isconst = T <: OrdinaryDiffEqConstantCache
  @unpack ηq, η₋₁, η₊₁, order, z, Δ = cache
  η = max(ηq, η₋₁, η₊₁)
  if η < integrator.opts.qsteady_max
    cache.η = 1
    cache.nextorder = order
  end

  if η == ηq
    cache.η = cache.ηq
    cache.nextorder = order
  elseif η == η₋₁
    cache.η = cache.η₋₁
    cache.nextorder = order - 1
  else
    cache.η = cache.η₊₁
    cache.nextorder = order + 1
    # TODO: BDF
    if integrator.alg.algorithm == :BDF
      if isconst
        z[end] = Δ
      else
        @. z[end] = Δ
      end
    end #endif BDF
  end # endif η == ηq
  return nothing
end

function stepsize_η!(integrator, cache, order)
  bias2 = integrator.alg.bias2
  addon = integrator.alg.addon
  L = order+1
  cache.ηq = inv( (bias2*integrator.EEst)^inv(L) + addon )
  return cache.ηq
end

# TODO: Check them
function stepsize_η₊₁!(integrator, cache::T, order) where T
  isconstcache = T <: OrdinaryDiffEqConstantCache
  isconstcache || ( @unpack atmp, ratetmp = cache )
  @unpack uprev, t, u = integrator
  @unpack z, c_LTE₊₁, dts, c_𝒟  = cache
  bias3 = integrator.alg.bias3
  addon = integrator.alg.addon
  q = order
  cache.η₊₁ = 0
  qmax = length(z)-1
  L = q+1
  if q != qmax
    cache.prev_𝒟 == 0 && return cache.η₊₁
    cquot = (c_𝒟 / cache.prev_𝒟) * (dts[1]/dts[2])^L
    if isconstcache
      atmp = muladd.(-cquot, z[end], cache.Δ)
      atmp = calculate_residuals(atmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
    else
      @. ratetmp = muladd(-cquot, z[end], cache.Δ)
      calculate_residuals!(atmp, ratetmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
    end
    dup = integrator.opts.internalnorm(atmp,t) * c_LTE₊₁
    cache.η₊₁ = inv( (bias3*dup)^inv(L+1) + addon )
  end
  return cache.η₊₁
end

function stepsize_η₋₁!(integrator, cache::T, order) where T
  isconstcache = T <: OrdinaryDiffEqConstantCache
  isconstcache || ( atmp = cache.atmp )
  @unpack uprev, t, u = integrator
  @unpack z, c_LTE₋₁ = cache
  bias1 = integrator.alg.bias1
  addon = integrator.alg.addon
  cache.η₋₁ = 0
  if order > 1
    if isconstcache
      atmp = calculate_residuals(z[order+1], uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
    else
      calculate_residuals!(atmp, z[order+1], uprev, u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm, t)
    end
    approx = integrator.opts.internalnorm(atmp,t) * c_LTE₋₁
    cache.η₋₁ = inv( (bias1*approx)^inv(order) + addon )
  end
  return cache.η₋₁
end


text:
export expected_n,
       expected_pq

"""
    expected_n(basis::Vector{Int},
               phi::Vector{T},
               psi::Vector{T},
               p::Vector{Int}) where {T<:Number}

Compute the expectation value <phi|n_p[1]*n_p[2]*...*n_p[n]|psi>.
"""
function expected_n(basis::Vector{Int},
                    phi::Vector{T},
                    psi::Vector{T},
                    p::Vector{Int}) where {T<:Number}
    o = zero(Float64)
    for i in eachindex(basis)
        bs = basis[i]
        has_all_n = true
        for pi in p
            (bs>>(pi-1))&1 == 0 && (has_all_n = false)
        end
        if has_all_n
            o += conj(phi[i])*psi[i]
        end
    end
    return o
end

"""
    expected_n(basis::Vector{Int},
               psi::Vector{T},
               p::Vector{Int}) where {T<:Number}

Compute the expectation value <psi|n_p[1]*n_p[2]*...*n_p[n]|psi>.
"""
function expected_n(basis::Vector{Int},
                    psi::Vector{T},
                    p::Vector{Int}) where {T<:Number}
    return expected_n(basis, psi, psi, p)
end

expected_n(basis, psi, p::Int) = expected_n(basis, psi, Int[p])
expected_n(basis, phi, psi, p::Int) = expected_n(basis, phi, psi, Int[p])

"""
    expected_pq(basis::Vector{Int},
                phi::Vector{T},
                psi::Vector{T},
                p::Vector{Int},
                q::Vector{Int}) where {T<:Number}

Compute <phi|b^+_p[1]*...*b^+_p[n]*b_q[1]*...*b_q[m]|psi>.
"""
function expected_pq(basis::Vector{Int},
                     phi::Vector{T},
                     psi::Vector{T},
                     p::Vector{Int},
                     q::Vector{Int}) where {T<:Number}

    p = deepcopy(p)
    q = deepcopy(q)
    (allunique(p) && allunique(q)) || return zero(T)
    # Number operators.
    n_i = p ∩ q
    # Remove those indices where there are number ops from p and q.
    setdiff!(p, n_i)
    setdiff!(q, n_i)

    o = zero(T)
    for i in eachindex(basis)
        bs = basis[i]
        has_all_pq = true
        # See if b^+_p*b_q acting on bs is not zero.
        for pi in p
            (bs>>(pi-1))&1 != 0 && (has_all_pq = false)
        end
        for qi in q
            (bs>>(qi-1))&1 != 1 && (has_all_pq = false)
        end
        for ni in n_i
            (bs>>(ni-1))&1 != 1 && (has_all_pq = false)
        end
        # If b^+_q*b_q|bs> not zero, compute the resulting state.
        if has_all_pq
            rs = bs
            for pi in p
                rs += 1<<(pi-1)
            end
            for qi in q
                rs -= 1<<(qi-1)
            end
            j = searchsortedfirst(basis, rs)
            o += psi[i]*conj(phi[j])
        end
    end
    return o
end

"""
    expected_pq(basis::Vector{Int},
                psi::Vector{<:Number},
                phi::Vector{<:Number},
                p::Int,
                q::Int)

Compute <phi|b^+_p*b_q|psi>.
"""
function expected_pq(basis::Vector{Int},
                     psi::Vector{T},
                     phi::Vector{T},
                     p::Int,
                     q::Int) where {T<:Number}

    p == q && return expected_n(basis, psi, phi, Int[p])

    o = zero(T)
    for i in eachindex(basis)
        bs = basis[i]
        if (bs>>(p-1))&1 == 0 && (bs>>(q-1))&1 == 1
            rs = bs + 1<<(p-1) - 1<<(q-1)
            j = searchsortedfirst(basis, rs)
            o += psi[i]*conj(phi[j])
        end
    end
    return o
end

expected_pq(basis, psi, p::Int, q::Int) = expected_pq(basis, psi, psi, p, q)
expected_pq(basis, psi, p::Vector{Int}, q::Vector{Int}) = expected_pq(basis, psi, psi, p, q)


text:
using Distributed
addprocs(2)
@everywhere begin
    using Pkg
    Pkg.activate("..")
    using Test
    using DistributedFFT
end

include("FFTtests.jl")


text:
n = 10 # rand(1:10)
@test matrixdepot("chow", n) == matrixdepot("chow", Float64, n)
@test matrixdepot("chow",Int, 5, 5, 0)[5,1] == 3125
@test diag(matrixdepot("chow", n, 2, 4)) == ones(Float64, n) * 6
println("'chow' passed test...")


text:
function AbstractTrees.printnode(io::IO, node::Node)
    print(io, join(["Node -> ID:", node.id, "Name: ", node.name, ", Link: ", node.MTG.link, "Index: ", node.MTG.index]))
end

"""
Print a node to io using an UTF-8 formatted representation of the `tree`.
Most of the code from [DataTrees.jl](https://github.com/vh-d/DataTrees.jl/blob/master/src/printing.jl)

# Examples

```julia
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)
mtg
# / 1: \$
# └─ / 2: Individual
#    └─ / 3: Axis
#       └─ / 4: Internode
#          ├─ + 5: Leaf
#          └─ < 6: Internode
#             └─ + 7: Leaf
```
"""
function Base.print(node::Node; leading::AbstractString = "", io::IO = stdout, limit = true)
    node_vec = get_printing(node; leading = leading)

    for (i, j) in enumerate(node_vec)
        print(io, j * "\n")
        limit && i >= displaysize(io)[1] && (print(io, "…"); break)
    end
end

function Base.print(node::Node, vars; leading::AbstractString = "", io::IO = stdout)
    DataFrame(node, vars)
end

function Base.show(io::IO, node::Node)
    print(node; io = io)
end

"""
    get_printing(node::Node; leading::AbstractString = "")

Format the printing of the tree according to link: follow or branching
"""
function get_printing(node::Node; leading::AbstractString = "")
    node_vec = [node.MTG.link * " " * string(node.id) * ": " * node.MTG.symbol]
    print_below = get_printing_(node; leading = "")
    if print_below !== nothing
        append!(node_vec, print_below)
    else
        node_vec
    end
end

function get_printing_(node::Node; leading::AbstractString = "")
    child_vec = Array{String,1}()
    if !isleaf(node)
        last = length(node.children)
        i = 0
        for chnode in ordered_children(node)
            i += 1
            if i != last
                to_print = leading * "\u251C\u2500 " * chnode.MTG.link * " " * string(chnode.id) * ": " * chnode.MTG.symbol
                new_leading = leading * "\u2502  "
            else
                to_print = leading * "\u2514\u2500 " * chnode.MTG.link * " " * string(chnode.id) * ": " * chnode.MTG.symbol
                new_leading = leading * "   "
            end
            append!(child_vec, [to_print])
            if length(children(chnode)) != 0
                children_vec = get_printing_(chnode; leading = new_leading)
                append!(child_vec, children_vec)
            end
        end
        return child_vec
    end
end