text:
@testset "Stamp" begin
    # Test Stamo and VStamp
    using Test
    N = 10
    p = 3

    vstamp1 = VStamp(:opt, 1, :x, 1, 5)
    @test valid(vstamp1,N,p)

    vstamp1 = VStamp(:opt, 4, :x, 1, 5)
    @test !valid(vstamp1,N,p)

    vstamp1 = VStamp(:opt, 2, :u, 2, 10)
    @test !valid(vstamp1,N,p)

    vstamp1 = VStamp(:opt, 2, :u, 2, 9)
    @test valid(vstamp1,N,p)

    vstamp1 = VStamp(:opt, 2, :u, 3, 9)
    @test !valid(vstamp1,N,p)

    vstamp1 = VStamp(:opt, 2, :z, 1, 1)
    @test !valid(vstamp1,N,p)

    vstamp1 = VStamp(:dyn, 1, :u, 2, 1)
    @test !valid(vstamp1,N,p)

    vstamp1 = VStamp(:dyn, 1, :x, 1, 1)
    @test valid(vstamp1,N,p)

    vstamp1 = VStamp(:dyn, 2, :x, 1, 1)
    @test !valid(vstamp1,N,p)

    vstamp1 = VStamp(:dyn, 2, :x, 2, 1)
    @test !valid(vstamp1,N,p)

    vstamp1 = VStamp(:opt, 2, :u, 2, 4)
    @test valid(vstamp1,N,p)

    hstamp = HStamp()
    @test !valid(hstamp, N, p)

    hstamp1 = HStamp(:u, 2, 3)
    @test valid(hstamp1,N,p)

    hstamp1 = HStamp(:u, 0, 3)
    @test !valid(hstamp1,N,p)

    hstamp1 = HStamp(:u, 1, 100)
    @test !valid(hstamp1,N,p)

    hstamp1 = HStamp(:x, 2, 2)
    @test !valid(hstamp1,N,p)

    hstamp1 = HStamp(:x, 1, 1)
    @test !valid(hstamp1,N,p)

    hstamp1 = HStamp(:x, 1, 2)
    @test valid(hstamp1,N,p)

    hstamp1 = HStamp(:λ, 1, 2)
    @test valid(hstamp1,N,p)

    hstamp1 = HStamp(:λ, 1, 2)
    @test valid(hstamp1,N,p)

    hstamp1 = HStamp(:λ, 0, 2)
    @test !valid(hstamp1,N,p)

    hstamp1 = HStamp(:λ, 2, 4)
    @test valid(hstamp1,N,p)

    hstamp1 = HStamp(:λ, 2, 40)
    @test !valid(hstamp1,N,p)

    stamp1 = Stamp(:opt, 1, :x, 1, 5, :x, 1, 3)
    @test valid(stamp1,N,p)

    stamp1 = Stamp(:opt, 1, :x, 1, 5, :u, 0, 3)
    @test !valid(stamp1,N,p)

    stamp1 = Stamp(:opt, 1, :x, 1, 5, :u, 3, 3)
    @test valid(stamp1,N,p)

    stamp1 = Stamp(:opt, 1, :x, 1, 5, :u, 1, 3)
    @test valid(stamp1,N,p)

    stamp1 = Stamp(:opt, 1, :x, 1, 5, :λ, 3, 3)
    @test !valid(stamp1,N,p)

    stamp1 = Stamp(:opt, 1, :x, 1, 5, :λ, 1, 3)
    @test valid(stamp1,N,p)

    stamp1 = Stamp(:opt, 1, :x, 1, 5, :λ, 1, 3)
    @test valid(stamp1,N,p)

    stamp1 = Stamp(:dyn, 1, :x, 1, 5, :x, 1, 3)
    @test valid(stamp1,N,p)

    stamp1 = Stamp(:dyn, 1, :x, 1, 5, :u, 2, 3)
    @test valid(stamp1,N,p)

    # Test equality
    stamp1 = Stamp(:opt, 1, :x, 1, 5, :λ, 1, 3)
    stamp2 = Stamp(:opt, 1, :x, 1, 5, :λ, 1, 3)
    stamp3 = Stamp(:opt, 2, :x, 1, 5, :λ, 1, 3)
    @test stamp1 == stamp2
    @test !(stamp1 == stamp3)
    @test !(Stamp() == VStamp())
    vstamp1 = VStamp()
    vstamp2 = VStamp()
    vstamp3 = VStamp(:opt, 2, :x, 1, 3)
    @test vstamp1 == vstamp2
    @test !(vstamp1 == vstamp3)
    hstamp1 = HStamp()
    hstamp2 = HStamp()
    hstamp3 = HStamp(:x, 1, 3)
    @test hstamp1 == hstamp2
    @test !(hstamp1 == hstamp3)

    # Test Stampify
    stamp0 = stampify(:opt, 1, :x_1, 1, 6, :λ, 1, 3)
    stamp1 = stampify(:opt, 1, :x, 1, 5, :λ, 1, 3)
    stamp2 = stampify(:opt, 1, :x1, 1, 4, :λ, 1, 3)
    @test stamp0 == stamp1
    @test stamp1 == stamp2

    stampify!(stamp2, :opt, 2, :x1, 1, 4, :λ, 1, 3)
    @test !(stamp1 == stamp2)

    vstamp0 = VStamp()
    vstamp1 = stampify(:opt, 1, :x, 1, 5)
    stampify!(vstamp0, :opt, 1, :x, 1, 5)
    @test vstamp0 == vstamp1
    @test isequal(vstamp0, vstamp1)

    hstamp0 = HStamp()
    hstamp1 = stampify(:x, 1, 5)
    stampify!(hstamp0, :x, 1, 5)
    @test hstamp0 == hstamp1
    @test isequal(hstamp0, hstamp1)

    @test !isequal(vstamp0, hstamp1)

end


text:
module HTSLIB

const debug = true
if !debug
@linux_only const LIBHTS = Libdl.find_library(["libhts.so"],[joinpath(dirname(dirname(@__FILE__)),"deps/usr/lib")])
@linux_only const LIBKSTR = joinpath(dirname(@__FILE__), "c/libkstr.so")
@osx_only const LIBHTS = Libdl.find_library(["libhts.dylib"],[joinpath(dirname(dirname(@__FILE__)),"deps/usr/lib")])
@osx_only const LIBKSTR = joinpath(dirname(@__FILE__), "c/libkstr.dylab")
end

function __init__()
    if !isfile(LIBKSTR)
        if OS_NAME == ":Windows"
            error("Windows is not supported now")
        else
            c_path = joinpath(dirname(@__FILE__),"c")
            cd(c_path) do
                run(`make libkstr`)
            end
        end
    end
end

import Base:readline,eof,write,close

export sam_open,
       bam_open,
       readline,
       eof,
       write,
       close

include("HT.jl")

end

ext: jl


text:
using ClearSky
using Documenter

DocMeta.setdocmeta!(ClearSky, :DocTestSetup, :(using ClearSky); recursive=true)

makedocs(;
    modules=[ClearSky],
    authors="Mark Baum <markmbaum@protonmail.com>",
    repo="https://github.com/markmbaum/ClearSky.jl/blob/{commit}{path}#{line}",
    sitename="ClearSky.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://markmbaum.github.io/ClearSky.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Absorption Data" => "absorption_data.md",
        "Line Shapes" => "line_shapes.md",
        "Gas Objects" => "gas_objects.md",
        "Atmospheric Profiles" => "atmospheric_profiles.md",
        "Radiation" => "radiation.md",
        "Modeling" => "modeling.md",
        "Orbits & Insolation" => "orbits_insolation.md"
    ],
)

deploydocs(;
    repo="github.com/markmbaum/ClearSky.jl",
    devbranch="main",
    versions=nothing
)


text:
using Base.Test
using XPP

#Define custom handlers to get consone-output in the event of a successful test
custom_handler(r::Test.Success) = println("Success on $(r.expr)")
custom_handler(r::Test.Failure) = error("test failed: $(r.expr)")
custom_handler(r::Test.Error)   = rethrow(r)

#Run test suite
Test.with_handler(custom_handler) do

    #Does the Module load properly?
    M = fromOdeFile("CaChannel.ode")
    loading_model = typeof(M) == XPP.Model
    @test loading_model

    #Is the modelstructure correct?
    #..i.e. do keys of odes-Dict and init-Dict match up?
    keys_ode = sort([k for k in keys(M.odes)])
    keys_init = sort([k for k in keys(M.init)])
    match_odes_init = keys_ode == keys_init
    @test match_odes_init

    #Can we run a simulation?
    runSimulation!(M, "test1")
    running_simulation = typeof(M.sims["test1"]) == XPP.SimulationData
    @test running_simulation

    #Are datasets from multiple simulations actually different?
    # ... there used to be problems with the way julia stores arrays!
    pars = [p for p in keys(M.pars)]
    selected_parameter = pars[rand(1:length(pars))]
    M.pars[selected_parameter] = M.pars[selected_parameter] * rand()
    runSimulation!(M, "test2")
    correct_storage_overall = M.sims["test1"] != M.sims["test2"]
    @test correct_storage_overall
    correct_storage_pars = M.sims["test1"].P[selected_parameter] != M.sims["test2"].P[selected_parameter]
    @test correct_storage_pars

    #Can we produce a plot?
    myPlot = plotModel(M, "test1")
    plot_model = myPlot != Union{}
    @test plot_model

    #Can we save to a json file?
    saveModel(M, "MSave.json")
    #size checkd
    save_json = sizeof("MSave.json") != 0
    @test save_json

    #Can we load from a json file?
    M2 = loadModel("MSave.json")
    load_json = typeof(M2) == XPP.Model
    @test load_json

    #Cleanup
    run(`rm MSave.json`)

end

max_stars_repo_name: novakgroupoxford/XPPjl

max_issues_repo_name: novakgroupoxford/XPPjl

max_forks_repo_name: novakgroupoxford/XPPjl


text:
# Binary heap (non-mutable)

include("arrays_as_heaps.jl")

#################################################
#
#   heap type and constructors
#
#################################################

"""
    FasterForward()

FasterForward enables 2x faster float comparison versus Base.ForwardOrdering,
but ordering is undefined if the data contains NaN values.
Enable this higher-performance option by calling the BinaryHeap
constructor instead of the BinaryMinHeap helper constructor.
"""
struct FasterForward <: Base.Ordering end
Base.lt(o::FasterForward, a, b) = a < b

"""
    FasterReverse()

FasterReverse enables 2x faster float comparison versus Base.ReverseOrdering,
but ordering is undefined if the data contains NaN values.
Enable this higher-performance option by calling the BinaryHeap
constructor instead of the BinaryMaxHeap helper constructor.
"""
struct FasterReverse <: Base.Ordering end
Base.lt(o::FasterReverse, a, b) = a > b

mutable struct BinaryHeap{T, O <: Base.Ordering} <: AbstractHeap{T}
    ordering::O
    valtree::Vector{T}

    function BinaryHeap{T}(ordering::Base.Ordering) where T
        new{T, typeof(ordering)}(ordering, Vector{T}())
    end

    function BinaryHeap{T}(ordering::Base.Ordering, xs::AbstractVector) where T
        valtree = heapify(xs, ordering)
        new{T, typeof(ordering)}(ordering, valtree)
    end
end

BinaryHeap(ordering::Base.Ordering, xs::AbstractVector{T}) where T = BinaryHeap{T}(ordering, xs)

# Constructors using singleton order types as type parameters rather than arguments
BinaryHeap{T, O}() where {T, O<:Base.Ordering} = BinaryHeap{T}(O())
BinaryHeap{T, O}(xs::AbstractVector) where {T, O<:Base.Ordering} = BinaryHeap{T}(O(), xs)

# These constructors needed for BinaryMaxHeap, until we have https://github.com/JuliaLang/julia/pull/37822
BinaryHeap{T, DefaultReverseOrdering}() where {T} = BinaryHeap{T}(Base.Reverse)
BinaryHeap{T, DefaultReverseOrdering}(xs::AbstractVector) where {T} = BinaryHeap{T}(Base.Reverse, xs)

# Forward/reverse ordering type aliases
const BinaryMinHeap{T} = BinaryHeap{T, Base.ForwardOrdering}
const BinaryMaxHeap{T} = BinaryHeap{T, DefaultReverseOrdering}

BinaryMinHeap(xs::AbstractVector{T}) where T = BinaryMinHeap{T}(xs)
BinaryMaxHeap(xs::AbstractVector{T}) where T = BinaryMaxHeap{T}(xs)

#################################################
#
#   interfaces
#
#################################################

"""
    length(h::BinaryHeap)

Returns the number of elements in heap `h`.
"""
Base.length(h::BinaryHeap) = length(h.valtree)

"""
    isempty(h::BinaryHeap)

Returns whether the heap `h` is empty.
"""
Base.isempty(h::BinaryHeap) = isempty(h.valtree)

"""
    push!(h::BinaryHeap, value)

Adds the `value` element to the heap `h`.
"""
@inline function Base.push!(h::BinaryHeap, v)
    heappush!(h.valtree, v, h.ordering)
    return h
end

"""
    first(h::BinaryHeap)

Returns the element at the top of the heap `h`.
"""
@inline Base.first(h::BinaryHeap) = h.valtree[1]

"""
    pop!(h::BinaryHeap)

Removes and returns the element at the top of the heap `h`.
"""
Base.pop!(h::BinaryHeap) = heappop!(h.valtree, h.ordering)

"""
    sizehint!(h::BinaryHeap, n::Integer)

Suggest that heap `h` reserve capacity for at least `n` elements. This can improve performance.
"""
function Base.sizehint!(h::BinaryHeap, n::Integer)
    sizehint!(h.valtree, n)
    return h
end


text:
summary_group = [:range, :any, :all, :eltype, :unique, :minimum, :maximum, :extrema, :first, :last]

for op in summary_group
    @eval begin
        function ($op)(x::RLEVector)
            ($op)(x.runvalues)
        end
    end
end

struct RLEVectorStyle <: Broadcast.AbstractArrayStyle{1} end
Base.BroadcastStyle(::Type{<:RLEVector}) = RLEVectorStyle()
RLEVectorStyle(::Val{0}) = RLEVectorStyle()
RLEVectorStyle(::Val{1}) = RLEVectorStyle()
RLEVectorStyle(::Val{N}) where {N} = Broadcast.DefaultArrayStyle{N}()
function Base.similar(bc::Broadcast.Broadcasted{RLEVectorStyle}, ::Type{ElType}) where {N,ElType}
    RLEVector(Vector{ElType}(undef, 1), [size(bc)[1]])
end
Base.map(f, x::RLEVector) = RLEVector(map(f, x.runvalues), ends(x))

## Methods that take two arguments, delegate to rle.runvalues and return something other than an RLEVector
Base.in(y::T1, x::RLEVector{T1,T2}) where {T1,T2<:Integer} = in(y, x.runvalues)
Base.in(x::RLEVector) = in(x.runvalues)

# Defaulting to fun(itr) for some things
for op in [:findmin, :findmax]
    @eval begin
        function ($op)(x::RLEVector)
            m = ($op)(x.runvalues)
            (m[1], starts(x, m[2]))
        end
    end
end

function indexin(a, b::RLEVector)
    inds = starts(b)
    bdict = Dict{eltype(b),eltype(inds)}()
    for (val, ind) in zip(b.runvalues, inds)
        get!(bdict, val, ind)
    end
    return Union{eltype(inds),Nothing}[get(bdict, i, nothing) for i in a]
end

function median(x::RLEVector)
    len = length(x)
    len <= 2 && return (middle(x.runvalues))
    sorted = sort(x)
    mid = cld(len, 2)
    mid_run = ind2run(sorted, mid)
    if mod(len, 2) == 0 && mid == sorted.runends[mid_run] # even numbered and at end of run, avg with next value
        median = middle(sorted.runvalues[mid_run], sorted.runvalues[mid_run+1])
    else
        median = middle(sorted.runvalues[mid_run])
    end
    return (median)
end

function sum(x::RLEVector{T1,T2}) where {T1,T2}
    rval = zero(T1)
    @simd for i = 1:nrun(x)
        @inbounds rval = rval + (x.runvalues[i] * widths(x, i))
    end
    return (rval)
end

mean(x::RLEVector) = rval = sum(x) / length(x)


text:
using StorageGraphs
using Test

testdir = dirname(@__FILE__)

tests = [
    "interface",
    "walk",
    "creation",
    "query",
    "persistence",
    "readme",
    "chaos_mwe"
]

@testset "StorageGraphs.jl" begin
    for t in tests
        tp = joinpath(testdir, "$(t).jl")
        include(tp)
    end
end


text:
mutable struct Sun3DReader <: RgbdReader
    frame_num::Int64
    f_id::Int64
    ids::Vector{Int64} # actually using image ids
    timestamps::Vector{Int64}
    data_folder::String
    depth_folder::String
    image_folder::String
    image_id2fname::Dict{Int64, String}
    image_id2time::Dict{Int64, Int64}
    depth_id2fname::Dict{Int64, String}
    depth_id2time::Dict{Int64, Int64}
    img_id_to_depth_id::Dict{Int64, Int64}
    cam::CameraModel.RgbdCamParams
    depth_factor::Float64
end

function Sun3DReader(data_folder::String)
    depth_folder = joinpath(data_folder, "depth")
    image_folder = joinpath(data_folder, "image")
    # get all depth file names
    depth_fnames = readdir(depth_folder)
    filter!(x->splitext(x)[2]==".png", depth_fnames)
    image_fnames = readdir(image_folder)
    filter!(x->splitext(x)[2]==".jpg", image_fnames)

    image_id2fname = Dict{Int64, String}()
    image_id2time = Dict{Int64, Int64}()

    for fname in image_fnames
        fname_wo_ext = splitext(fname)[1]
        id, time = parse.(Int64, split(fname_wo_ext, "-"))
        # id_time = parse(Int64,replace(splitext(fname)[1], "-"=>""))
        image_id2fname[id] = fname
        image_id2time[id] = time
    end

    depth_id2fname = Dict{Int64, String}()
    depth_id2time = Dict{Int64, Int64}()
    depth_ids = Vector{Int64}(undef, length(depth_fnames))
    depth_times = Vector{Int64}(undef, length(depth_fnames))

    for (idx, fname) in enumerate(depth_fnames)
        fname_wo_ext = splitext(fname)[1]
        id, time = parse.(Int64, split(fname_wo_ext, "-"))
        # id_time = parse(Int64,replace(splitext(fname)[1], "-"=>""))
        depth_id2fname[id] = fname
        depth_id2time[id] = time
        depth_ids[idx] = id
        depth_times[idx] = time
    end

    # synchronize: find a depth for each image
    img_id_to_depth_id = Dict{Int64, Int64}()
    for (img_id, img_time) = image_id2time
        time_diff = abs.(img_time .- depth_times)
        matched_depth_idx = argmin(time_diff)
        # @assert time_diff[matched_depth_idx] < 3000 # 3000 is 3ms
        if time_diff[matched_depth_idx] > 3000
            @warn "sun3d data image depth time diff is $(time_diff[matched_depth_idx]/1000) ms"
        end
        img_id_to_depth_id[img_id] = depth_ids[matched_depth_idx]
    end

    # load camera intrinsic
    intrinsic_fname = joinpath(data_folder, "intrinsics.txt")
    k = DelimitedFiles.readdlm(intrinsic_fname, ' ', Float64)
    fx = k[1,1]
    fy = k[2,2]
    cx = k[1,3]
    cy = k[2,3]
    w, h = 640, 480

    depth_factor=1000.0

    cam = CameraModel.RgbdCamParams(
        fx=fx,
        fy=fy,
        cx=cx,
        cy=cy,
        width=w,
        height=h,
        depth_factor=depth_factor,
        baseline=50.0/fx
    )

    ids = sort!(collect(keys(img_id_to_depth_id)))
    f_id::Int64 = ids[1]

    # for sun2d, timestamp = id
    timestamps = copy(ids)

    res = Sun3DReader(
        length(ids),
        f_id,
        ids, # actually using image ids
        timestamps,
        data_folder,
        depth_folder,
        image_folder,
        image_id2fname,
        image_id2time,
        depth_id2fname,
        depth_id2time,
        img_id_to_depth_id,
        cam,
        depth_factor
    )

    return res
end

function Sun3DReader(
    data_folder::String, index_fpath::String; one_indexed::Bool=true
)
    depth_folder = joinpath(data_folder, "depth")
    image_folder = joinpath(data_folder, "image")
    # # get all depth file names
    # depth_fnames = readdir(depth_folder)
    # filter!(x->splitext(x)[2]==".png", depth_fnames)
    # image_fnames = readdir(image_folder)
    # filter!(x->splitext(x)[2]==".jpg", image_fnames)
    ids, id2rgb, id2depth = load_boundlefusion_index(
        index_fpath,
        seq_folder=""
    )

    image_id2fname = Dict{Int64, String}()
    image_id2time = Dict{Int64, Int64}()

    # for fname in image_fnames
    for (img_id, img_path) in id2rgb
        # fname_wo_ext = splitext(splitpath(img_path)[end])[1]
        fname_wo_ext = splitext(img_path)[1]
        id, time = parse.(Int64, split(fname_wo_ext, "-"))
        # id_time = parse(Int64,replace(splitext(fname)[1], "-"=>""))
        image_id2time[img_id+one_indexed] = time
        image_id2fname[img_id+one_indexed] = img_path
    end

    depth_id2fname = Dict{Int64, String}()
    depth_id2time = Dict{Int64, Int64}()
    # TODO what is the use of this ids
    depth_ids = Vector{Int64}(undef, length(id2depth))
    depth_times = Vector{Int64}(undef, length(id2depth))

    # for (idx, fname) in enumerate(depth_fnames)
    for (depth_id, depth_fpath) in id2depth
        # fname_wo_ext = splitext(splitpath(depth_fpath)[end])[1] #splitext(fname)[1]
        fname_wo_ext = splitext(depth_fpath)[1]
        id, time = parse.(Int64, split(fname_wo_ext, "-"))
        # id_time = parse(Int64,replace(splitext(fname)[1], "-"=>""))
        depth_id2fname[depth_id+one_indexed] = depth_fpath
        depth_id2time[depth_id+one_indexed] = time
        # depth_ids[depth_id+1] = id
        # depth_times[depth_id] = time
    end

    # synchronize: find a depth for each image
    ids .+= one_indexed
    img_id_to_depth_id = Dict{Int64, Int64}()
    for fid in ids
        img_id_to_depth_id[fid] = fid
    end
    # for (img_id, img_time) = image_id2time
    #     time_diff = abs.(img_time .- depth_times)
    #     matched_depth_idx = argmin(time_diff)
    #     # @assert time_diff[matched_depth_idx] < 3000 # 3000 is 3ms
    #     if time_diff[matched_depth_idx] > 3000
    #         @warn "sun3d data image depth time diff is $(time_diff[matched_depth_idx]/1000) ms"
    #     end
    #     img_id_to_depth_id[img_id] = depth_ids[matched_depth_idx]
    # end

    # load camera intrinsic
    intrinsic_fname = joinpath(data_folder, "intrinsics.txt")
    k = DelimitedFiles.readdlm(intrinsic_fname, ' ', Float64)
    fx = k[1,1]
    fy = k[2,2]
    cx = k[1,3]
    cy = k[2,3]
    w, h = 640, 480

    depth_factor=1000.0

    cam = CameraModel.RgbdCamParams(
        fx=fx,
        fy=fy,
        cx=cx,
        cy=cy,
        width=w,
        height=h,
        depth_factor=depth_factor,
        baseline=50.0/fx
    )

    # ids = sort!(collect(keys(img_id_to_depth_id)))
    f_id::Int64 = ids[1]

    # for sun3d, timestamp = id
    timestamps = copy(ids)

    res = SLAMData.Sun3DReader(
        length(ids),
        f_id,
        ids, # actually using image ids
        timestamps,
        data_folder,
        depth_folder,
        image_folder,
        image_id2fname,
        image_id2time,
        depth_id2fname,
        depth_id2time,
        img_id_to_depth_id,
        cam,
        depth_factor
    )

    return res
end

function readid(reader::Sun3DReader, id::Int64)
    # id = tum_data.ids[k]
    img = FileIO.load(joinpath(reader.image_folder, reader.image_id2fname[id]))

    depth_id = reader.img_id_to_depth_id[id]
    depth = FileIO.load(joinpath(reader.depth_folder, reader.depth_id2fname[depth_id]))
    depth_view = Images.channelview(depth)
    depth_uint16 = if ndims(depth_view) ==  3
        copy(reinterpret(UInt16, depth_view[1,:,:]))
    else
        copy(reinterpret(UInt16, depth_view))
    end
    # depth_uint16 = copy(reinterpret(UInt16, depth_view))
    depth_uint16 = (depth_uint16 .<< 13) .| (depth_uint16 .>>> 3)
    return img, depth_uint16
end

# function depth = depthRead(filename)
#     depth = imread(filename);
#     depth = bitor(bitshift(depth,-3), bitshift(depth,16-3));
#     depth = single(depth)/1000;
# end
"""
readRGBD(reader::Sun3DReader; id::Int64)
read rgbd images with id as input
"""
readRGBD(reader::Sun3DReader; id::Int64) = readid(reader, id)

"""
sun3d traj use Twc
"""
function load_sun3d_extrinics(traj_fpath::String; start_id::Int64=0)
    sun3d_raw_data = DelimitedFiles.readdlm(traj_fpath)

    Twcs = Dict{Int64,Matrix{Float64}}()

    @assert size(sun3d_raw_data, 1) % 3 == 0

    n = Int64(size(sun3d_raw_data, 1) / 3)

    for i1 = 1:n
        row_offset = (i1-1)*3
        T = Matrix{Float64}(LinearAlgebra.I, 4, 4)
        T[1:3, 1:4] = sun3d_raw_data[row_offset+1:row_offset+3, 1:4]
        Twcs[i1-1+start_id] = T
    end
    return Twcs
end


text:
" Returns the difference between two strings, in
  reverse order to the input. "
function reversed_non_matching(ionestr::String, itwostr::String)
  unmatchedone = ""
  unmatchedtwo = ""

  l1 = length(ionestr)
  l2 = length(itwostr)

  #Comparison part, build strings of non-matches
  lrcmplen = min(l1, l2)
  for i in 1:lrcmplen
    if ionestr[i] != itwostr[i]
      unmatchedone = string(ionestr[i]) * unmatchedone
      unmatchedtwo = string(itwostr[i]) * unmatchedtwo
    end
  end

  #Round off longer string with remaining characters
  for i in lrcmplen+1:max(l1,l2)
    if i <= l1
      unmatchedone = string(ionestr[i]) * unmatchedone
    elseif i <= l2
      unmatchedtwo = string(itwostr[i]) * unmatchedtwo
    end
  end
  return [unmatchedone,unmatchedtwo]
end

" Table lookup which builds an array, checking all the bins. "
function table_lookup_plural(chr::Char, table::Array{String, 1})
  return map(x -> '0'+x, filter(pos -> chr in table[pos], 1:length(table)))
end

"""
  `match_rating(onestr, twostr)`

  Return the similarity measure between two strings as an integer `0:6`.

  The higher output is better. Typically, it would be compared against a table
  of minimum match ratings, based on the original length of the strings. This
  is implemented in `meets_match_rating(onestr, twostr)`.

  If the length of the two strings differs by more than `2`, the result will be
  `-1`.

  - For the encoding used by the system, call `match_rating_encode(str)`

  - For an automatic binary response about the closeness of strings, call
    `meets_match_rating(onestr, twostr)`
"""
function match_rating(onestr::String, twostr::String)

  #encode input
  ionestr = match_rating_encode(onestr)
  itwostr = match_rating_encode(twostr)

  #If length difference is 3 or greater, no comparison.
  if abs(length(ionestr) - length(itwostr)) > 2
    return -1
  end

  #Remove identical characters in encoded string
  unmatchedone, unmatchedtwo = reversed_non_matching(ionestr, itwostr)

  #Remove identical characters in leftovers
  reunmatchedone, reunmatchedtwo = reversed_non_matching(unmatchedone, unmatchedtwo)

  #Pick longest string
  cmpstr = reunmatchedone
  if length(reunmatchedone) < length(reunmatchedtwo)
    cmpstr = reunmatchedtwo
  end

  #Return 6-remaining error count
  return 6-length(cmpstr)
end

"""
  `meets_match_rating(onestr, twostr)`

  Returns a Bool indicating whether two strings are sufficiently similar to be
  matched, according the the Match Rating Approach.

  A match rating value from `match_rating(onestr, twostr)` is compared to a minimum
  rating threshold based on the combined length of the input strings.

  - For the actual similarity measure, use `match_rating(onestr, twostr)`.
  - For the encoding used by the system, call `match_rating_encode(str)`
"""
function meets_match_rating(onestr::String, twostr::String)

  ionestr = match_rating_encode(onestr)
  itwostr = match_rating_encode(twostr)

  rating(x) = if x <= 4
      return 5
    elseif x <= 7
      return 4
    elseif x <= 12
      return 3
    else
      return 2
  end

  #calculate minimum rating based on length
  lensum = length(ionestr) + length(itwostr)
  minrating = rating(lensum)

  #get similarity of strings
  similarity = match_rating(onestr, twostr)

  #return comparison
  return similarity >= minrating
end

"""
  `editex(onestr, twostr)`

  Return the Editex distance between two strings.

  Editex is an approximate phonetic string comparison algorithm. It can be
  thought of as a kind of combination of Soundex and edit distance, though
  the table used differs from that of Soundex and the edit distance cost is
  calculated differently.

  In brief, a dynamic programming solution is used to minimise the edit cost,
  which is comprised at base of:

  - `0` where two characters are the same.
  - `1` where two characters are from the same phonetic group.
  - `2` otherwise.

  Complexity will be `O(mn)` for lengths `m` and `n` of input strings.

  See also: `code_similarity`.
"""
function editex(onestr::String, twostr::String)

  #Prepare soundexily
  ionestr = prep(onestr)
  itwostr = prep(twostr)

  l1 = length(ionestr)
  l2 = length(itwostr)

  editcost = editex_internal(l1, l2, ionestr, itwostr)
  return editcost
end

editex_table = ["aeiouwy", "bp", "ckq", "dt", "lr", "mn", "gj", "fpv", "sxz", "csz"]

" r-hand-side of editex internals "
function editex_r(charone, chartwo)
  if charone == chartwo
    return 0
  elseif length(intersect(table_lookup_plural(charone, editex_table), table_lookup_plural(chartwo, editex_table))) > 0
    return 1
  else
    return 2
  end
end

" d-hand side of editex internals. Mostly identical."
function editex_d(charone, chartwo)
  if charone == chartwo
    return 0
  elseif charone in "hw" || length(intersect(table_lookup_plural(charone, editex_table), table_lookup_plural(chartwo, editex_table))) > 0
    return 1
  else
    return 2
  end
end

" Dynamic programming implementation. "
function editex_internal(i, j, onestr, twostr)
  #terminate at start of strings
  if i == j && i == 1
    # Cost of both deletions (other move) would always be >= cost of comparison.
    return editex_r(onestr[i], twostr[j])
  #finished one string, eat other
  elseif j == 1
    return editex_internal(i - 1, 1, onestr, twostr) + editex_d(onestr[i-1], onestr[i])
  #finished the other
  elseif i == 1
    return editex_internal(1, j - 1, onestr, twostr) + editex_d(twostr[j-1], twostr[j])
  # minimum of deleting from either string or cost of cross-comparison
  else
    return min(editex_internal(i - 1, j, onestr, twostr) + editex_d(onestr[i-1], onestr[i]),
               editex_internal(i, j - 1, onestr, twostr) + editex_d(twostr[j-1], twostr[j]),
               editex_internal(i - 1, j - 1 , onestr, twostr) + editex_r(onestr[i], twostr[j]))
  end
end

export match_rating, meets_match_rating, editex


text:
using BrainSim
using Test

@testset "BrainSim.jl" begin
    # Write your tests here.
end


text:
# manually
@test compile_css("a") == Spiderman.IsTag("a")
@test compile_css(".a") == Spiderman.HasClass(r"\ba\b")
@test compile_css("#a") == Spiderman.HasId("a")
@test compile_css("a#b.c") == Spiderman.MatchGroup([Spiderman.IsTag("a"), Spiderman.HasId("b"), Spiderman.HasClass("c")])
@test compile_css("a b") == [Spiderman.IsTag("a"), Spiderman.IsTag("b")]

# macro'd
@test css"a" == Spiderman.IsTag("a")
@test css".a" == Spiderman.HasClass(r"\ba\b")
@test css"#a" == Spiderman.HasId("a")
@test css"a#b.c" == Spiderman.MatchGroup([Spiderman.IsTag("a"), Spiderman.HasId("b"), Spiderman.HasClass("c")])
@test css"a b" == [Spiderman.IsTag("a"), Spiderman.IsTag("b")]

# testing data
basic = parsehtml("""

    Some text here
    and more here

""")

# one-level matches
@test select(css"#one", basic) == [t("

")]
@test select(css"#two", basic) == [t("

")]
@test select(css".e", basic) == [t("Some text here")]
@test select(css".b", basic) == [t("

"), t("")]

@test select(css"a.b", basic) == [t("")]
@test select(css"#one.b", basic) == [t("

")]
@test select(css"p.b", basic) == [t("

")]

@test select(css"p", basic) == [t("

"), t("

")]
@test select(css"span", basic) == [t(""),
                                   t("Some text here"),
                                   t("and more here")]

# multi-level matches
@test select(css".d2 .d", basic) == [t("")]
@test select(css".d2 span.e", basic) == [t("Some text here")]


text:
"""
```
nearest_spd(A::Matrix)
```
Calculates the nearest (in Frobenius norm) symmetric positive semidefinite matrix to A.

From Higham: The nearest symmetric positive semidefinite matrix in the
Frobenius norm to an arbitrary real matrix A is shown to be (B + H)/2,
 where H is the symmetric polar factor of B=(A + A')/2.

http://www.sciencedirect.com/science/article/pii/0024379588902236

### Inputs

- `A::Matrix`: A square matrix

### Outputs

- `Ahat::Matrix`: The symmetric, square positive semidefinite matrix closest to A
"""
function nearest_spd(A::Matrix)
    @show "WHAT ARE YOU DOING? WHY ARE YOU CALLING THIS?"
    return "WHAT ARE YOU DOING? WHY ARE YOU CALLING THIS?"
end
#=function nearest_spd(A::Matrix)

    if size(A,1) != size(A,2)
        error("A must be a square matrix")
    end
    # symmetrize A into B
    B = (A + A')/2;

    # Compute the symmetric polar factor of B. Call it H.
    # Clearly H is itself SPD.
    U, Σ, V = svd(B);
    H = V*diagm(Σ)*V';

    Ahat = (B+H)/2;
    # Ensure symmetry
    Ahat = (Ahat + Ahat')/2;

    # Test that Ahat is in fact PD. if it is not so, then tweak it just a bit.
    k = 0
    chol_success = false

    eig_matrix = eigvals(Ahat)
    while !chol_success
        try
            chol(Ahat)
            chol_success = true
        catch
        end
        k = k + 1
        if !chol_success
            # mineig_factor is used in the event that the matrix has very small, slightly
            # negative eigenvalues, then the identity scaled by the absolute value of the eigenvalue must
            # be added so as to make the matrix "more" PD and not less.
            # If the minimum eigenvalue is 0 then mineig_factor will permanently remain at 0
            # and the Ahat matrix will cease to change in the iteration, so the ternary
            # condition sets mineig_factor to be the machine epsilon for Float64 numbers
            # so that the algorithm will continue to perturb Ahat in the right direction and not
            # get stuck.

            eig_list = eigvals(Ahat)
            mineig = eig_list[1]
            mineig_factor = mineig != 0 ? abs(mineig) : eps(Float64)
            Ahat = Ahat + (mineig_factor*k^2 + eps(Float64))*eye(size(A)...)
            eig_matrix = hcat(eig_matrix, eigvals(Ahat))
        end
    end

    return Ahat
end =#


text:
# test_old_interface.jl

@testset "DB Connection" begin
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    @test isa(con, DuckDB.Connection)
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

@testset "Test append DataFrame" begin
    # Open the database
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)

    # Create the table the data is appended to
    DuckDB.execute(
        con,
        "CREATE TABLE dtypes(bool BOOLEAN, tint TINYINT, sint SMALLINT, int INTEGER, bint BIGINT, utint UTINYINT, usint USMALLINT, uint UINTEGER, ubint UBIGINT, float FLOAT, double DOUBLE, date DATE, time TIME, vchar VARCHAR, nullval INTEGER)"
    )

    # Create test DataFrame
    input_df = DataFrame(
        bool = [true, false],
        tint = Int8.(1:2),
        sint = Int16.(1:2),
        int = Int32.(1:2),
        bint = Int64.(1:2),
        utint = UInt8.(1:2),
        usint = UInt16.(1:2),
        uint = UInt32.(1:2),
        ubint = UInt64.(1:2),
        float = Float32.(1:2),
        double = Float64.(1:2),
        date = [Dates.Date("1970-04-11"), Dates.Date("1970-04-12")],
        time = [Dates.Time(0, 0, 0, 100, 0), Dates.Time(0, 0, 0, 200, 0)],
        vchar = ["Foo", "Bar"],
        nullval = [missing, Int32(2)]
    )

    # append the DataFrame to the table
    DuckDB.appendDataFrame(input_df, con, "dtypes")

    # Output the data from the table
    output_df = DuckDB.toDataFrame(con, "select * from dtypes;")

    # Compare each column of the input and output dataframe with each other
    for (col_pos, input_col) in enumerate(eachcol(input_df))
        @test isequal(input_col, output_df[:, col_pos])
    end

    # Disconnect and close the database
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

@testset "Test README" begin
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    res = DuckDB.execute(con, "CREATE TABLE integers(date DATE, jcol INTEGER)")
    res = DuckDB.execute(con, "INSERT INTO integers VALUES ('2021-09-27', 4), ('2021-09-28', 6), ('2021-09-29', 8)")
    res = DuckDB.execute(con, "SELECT * FROM integers")
    df = DuckDB.toDataFrame(res)
    @test isa(df, DataFrame)
    df = DuckDB.toDataFrame(con, "SELECT * FROM integers")
    println(typeof(df))
    @test isa(df, DataFrame)
    DuckDB.appendDataFrame(df, con, "integers")
    DuckDB.disconnect(con)
    DuckDB.close(db)
end
#
@testset "HUGE Int test" begin
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    res = DuckDB.execute(con, "CREATE TABLE huge(id INTEGER,data HUGEINT);")
    res = DuckDB.execute(con, "INSERT INTO huge VALUES (1,NULL), (2, 1761718171), (3, 171661889178);")
    res = DuckDB.toDataFrame(con, "SELECT * FROM huge")
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

@testset "Interval type" begin
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    res = DuckDB.execute(con, "CREATE TABLE interval(interval INTERVAL);")
    res = DuckDB.execute(
        con,
        """
INSERT INTO interval VALUES
(INTERVAL 5 HOUR),
(INTERVAL 12 MONTH),
(INTERVAL 12 MICROSECOND),
(INTERVAL 1 YEAR);
"""
    )
    res = DuckDB.toDataFrame(con, "SELECT * FROM interval;")
    @test isa(res, DataFrame)
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

@testset "Timestamp" begin
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    res = DuckDB.execute(con, "CREATE TABLE timestamp(timestamp TIMESTAMP , data INTEGER);")
    res = DuckDB.execute(
        con,
        "INSERT INTO timestamp VALUES ('2021-09-27 11:30:00.000', 4), ('2021-09-28 12:30:00.000', 6), ('2021-09-29 13:30:00.000', 8);"
    )
    res = DuckDB.execute(con, "SELECT * FROM timestamp;")
    res = DuckDB.toDataFrame(res)
    @test isa(res, DataFrame)
    DuckDB.disconnect(con)
    DuckDB.close(db)
end

@testset "Items table" begin
    db = DuckDB.open(":memory:")
    con = DuckDB.connect(db)
    res = DuckDB.execute(con, "CREATE TABLE items(item VARCHAR, value DECIMAL(10,2), count INTEGER);")
    res = DuckDB.execute(con, "INSERT INTO items VALUES ('jeans', 20.0, 1), ('hammer', 42.2, 2);")
    res = DuckDB.toDataFrame(con, "SELECT * FROM items;")
    @test isa(res, DataFrame)
    DuckDB.disconnect(con)
end

@testset "Integers and dates table" begin
    db = DuckDB.DB()
    res = DBInterface.execute(db, "CREATE TABLE integers(date DATE, data INTEGER);")
    res =
        DBInterface.execute(db, "INSERT INTO integers VALUES ('2021-09-27', 4), ('2021-09-28', 6), ('2021-09-29', 8);")
    res = DBInterface.execute(db, "SELECT * FROM integers;")
    res = DuckDB.toDataFrame(res)
    @test res.date == [Date(2021, 9, 27), Date(2021, 9, 28), Date(2021, 9, 29)]
    @test isa(res, DataFrame)
    DBInterface.close!(db)
end

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


text:
function cluster2p(x::SVector{3,T}, cluster::Cluster{T}, particles::Particles, parindices::Vector{Int64}; p_kernel, n, eta, stretch=SVector(1.0,1.0,1.0)) where {T}
    efield = SVector(0.0,0.0,0.0)
    bfield = SVector(0.0,0.0,0.0)
    if cluster.children === nothing
        index_start = cluster.pindex_start
        index_end = index_start + cluster.npar - 1
        efield, bfield = cluster2p_brutal(x, particles, parindices, index_start, index_end)
    elseif admissible(x, cluster.bbox; eta=eta, stretch=stretch)
        xcoords = cluster.xcoords
        ycoords = cluster.ycoords
        zcoords = cluster.zcoords
        gamma_hat = cluster.gamma_hat
        mom_hat = cluster.mom_hat
        for k in 1:(n+1)
            sz = zcoords[k]
            for j in 1:(n+1)
                sy = ycoords[j]
                for i in 1:(n+1)
                    sx = xcoords[i]
                    p_hat = mom_hat[i,j,k]
                    K = kernel_relativity(x, SVector(sx,sy,sz), p_kernel)
                    efield += gamma_hat[i,j,k]*K
                    bfield += cross(p_hat, K)
                end
            end
        end
    else
        efield1, bfield1 = cluster2p(x, cluster.children[1], particles, parindices; p_kernel=p_kernel, n=n, eta=eta, stretch=stretch)
        efield2, bfield2 = cluster2p(x, cluster.children[2], particles, parindices; p_kernel=p_kernel, n=n, eta=eta, stretch=stretch)
        efield = efield1 + efield2
        bfield = bfield1 + bfield2
    end
    return efield, bfield
end

function admissible(x::SVector{3,T}, bbox::BBox{T}; stretch::SVector{3,T}=SVector(1.0,1.0,1.0), eta::T) where {T}
    bmin = bbox.bmin
    bmax = bbox.bmax
    xc = (bmin + bmax)/2.0
    r = norm(stretch .* (bmax - xc))
    R = norm(stretch .* (x - xc))
    return r/R < eta
end

function cluster2p_brutal(x::SVector{3,T}, particles::Particles, parindices::Vector{Int64}, index_start, index_end) where {T}
    efield = SVector(0.0,0.0,0.0)
    bfield = SVector(0.0,0.0,0.0)
    pos = particles.positions
    mom = particles.momentums
    for k in index_start:index_end
        xj = pos[parindices[k]]
        pj = mom[parindices[k]]
        K = kernel_relativity(x, xj, pj)
        efield += sqrt(1.0 + dot(pj, pj))*K
        bfield += cross(pj, K)
    end
    return efield, bfield
end


text:
export Integral

"""
A one-form is a map `V → K`, i.e., it is a generic function of one argument.

A one-form is applied to an argument using the `F(f)` syntax, which itself
invokes `form(F::OneForm, f)`. The `form` function can be specialized for
concrete types of `F` and `f`.

The main use cases in `DomainIntegrals.jl` are a lazy `Integral` object and a
quadratic form that derives from a bilinear form.
"""
abstract type OneForm end

(F::OneForm)(f) = form(F, f)

"`Integral` is a functional that represents integration of a function."
abstract type AbstractIntegral <: OneForm end

domain(F::AbstractIntegral) = F.domain
measure(F::AbstractIntegral) = F.measure
properties(F::AbstractIntegral) = F.properties

form(F::AbstractIntegral, f) = integral(f, domain(F), measure(F), properties(F)...)

struct Integral <: AbstractIntegral
    integrand
    domain
    measure
    properties
end

Integral(f, args...) = Integral(f, process_arguments(args...)...)

Integral(integrand, domain::Domain, measure::Measure, properties::Property...) =
    Integral(integrand, domain, measure, properties)

"""
A two-form is a map `V x W → K`, i.e., it is a generic function of two arguments.

The two-form is applied to two arguments using the `F(f,g)` syntax, which itself
invokes `form(F::TwoForm, f, g)`. The `form` function can be specialized for
concrete types of `F`, `f` and `g`.
"""
abstract type TwoForm end

(F::TwoForm)(f, g) = form(F, f, g)

"A bilinear form is linear in its two arguments."
abstract type BilinearForm <: TwoForm end

"A sesquilinear form is linear in the first and conjugate linear in the second argument."
abstract type SesquilinearForm <: TwoForm end

isbilinear(F::BilinearForm) = true
isbilinear(F::SesquilinearForm) = false

issesquilinear(F::BilinearForm) = false
issesquilinear(F::SesquilinearForm) = true

"""
Supertype of quadratic forms that derive from a two-form.

If `F` is a TwoForm, then we call `F(f,f)` the associated quadratic form.
"""
abstract type AbstractQuadraticForm <: OneForm end

# for convenience, assume the underlying form is stored in the field F
superform(F::AbstractQuadraticForm) = F.F

form(F::AbstractQuadraticForm, f) = quadraticform(superform(F), f)
quadraticform(F, f) = form(F, f, f)

"A quadratic form derived from a stored two-form."
struct QuadraticForm{FORM} <: AbstractQuadraticForm
    F   ::  FORM
end

"""
Apply a two-form to the given arguments.

`form(F::TwoForm, f, g)`

By default, `form` invokes `form1` which in turn invokes `form2`. The `form`
function itself can be specialized on the type of `F` without ambiguity.
Similarly, `form1` and `form2` can be specialized without ambiguity on the types
of `f` and `g` respectively.

Each of these functions can be specialized on all three types involved, as long
as the main type is at least as specific as any other specialization.
"""
form(F::TwoForm, f, g) = form1(F, f, g)

"Helper function for `form(F, f, g)` that can be specialized on the type of `f`."
form1(F, f, g) = form2(F, f, g)

"Helper function for `form(F, f, g)` that can be specialized on the type of `g`."
form2(F, f, g) = default_form(F, f, g)

functionproduct(f, g) = x -> f(x)*g(x)
functionconjugateproduct(f, g) = x -> f(x)*conj(g(x))

# By default we assume that the form is induced by a measure.
default_form(F::BilinearForm, f, g) =
    integral(functionproduct(f,g), measure(F))
default_form(F::SesquilinearForm, f, g) =
    integral(functionconjugateproduct(f,g), measure(F))

##################
## Inner products
##################

"An inner-product is a positive definite sesquilinear form."
abstract type InnerProduct <: SesquilinearForm end

# For convenience, assume the underlying form is stored in the field F.
superform(F::InnerProduct) = F.F

ishermitian(F::InnerProduct) = true
ispositivesemidefinite(F::InnerProduct) = true
ispositivedefinite(F::InnerProduct) = true

# an inner product is typically implemented in terms of another form
form(F::InnerProduct, f, g) = form(superform(F), f, g)

# An example concrete type follows. Other concrete types may be defined that
# return a form on-the-fly by defining `superform`, for example.

"A generic inner product that stores a positive definite Hermitian sesquilinear form."
struct GenericInnerProduct{FORM} <: InnerProduct
    F   ::  FORM
end

#########################
## Measure induced forms
#########################

"""
A bilinear form induced by a measure.

It is the case that `F(f,g)` equals `integral(x -> f(x)*g(x), measure(F))`.
"""
abstract type InducedBilinearForm <: BilinearForm end

issymmetric(F::InducedBilinearForm) = true
isconjugatesymmetric(F::InducedBilinearForm) = false
ishermitian(F::InducedBilinearForm) = false

"""
A sesquilinear form induced by a measure.

It is the case that `F(f,g)` equals `integral(x -> f(x)*conj(g(x)), measure(F))`.
"""
abstract type InducedSesquilinearForm <: SesquilinearForm end

issymmetric(F::InducedSesquilinearForm) = false
isconjugatesymmetric(F::InducedSesquilinearForm) = true
ishermitian(F::InducedSesquilinearForm) = isreal(measure(F))
ispositivesemidefinite(F::InducedSesquilinearForm) = ishermitian(F)

# Two concrete types follow. Other concrete types may be defined that do not
# store a measure, but merely return it on-the-fly by implementing `measure(F)`.

"A bilinear form induced by a generic measure."
struct MeasureBilinearForm{M} <: InducedBilinearForm
    measure     ::  M
end

"A sesquilinear form induced by a generic measure."
struct MeasureSesquilinearForm{M} <: InducedSesquilinearForm
    measure     ::  M
end


text:
# Заполняет словарь d индексами индивидуальных значений
#=
function indsdict!(d::Dict{T}, cdata::Tuple) where T
    @inbounds for (i, element) in enumerate(zip(cdata...))
        ind = ht_keyindex(d, element)
        if ind > 0
            push!(d.vals[ind], i)
        else
            v = Vector{Int}(undef, 1)
            v[1] = i
            d[element] = v
        end
    end
    d
end
=#
nonunique(v) = [k for (k, v) in StatsBase.countmap(v) if v > 1]

function floatparse(data)
    v = Vector{Float64}(undef, length(data))
    @inbounds for i = 1:length(data)
        if !isa(data[i], AbstractFloat) && !ismissing(data[i])
            if isa(data[i], AbstractString)
                try
                    v[i] = parse(Float64, data[i])
                catch
                    v[i] = NaN
                end
            else
                try
                    v[i] = float(data[i])
                catch
                    v[i] = NaN
                end
            end
        elseif ismissing(data[i])
            v[i] = NaN
        else
            v[i] = data[i]
        end
    end
    identity.(v)
end

function checkvalues(timevals_sp, concvals_sp)
    if !(eltype(timevals_sp) <: Number)
        timevals_sp = identity.(timevals_sp)
        eltype(timevals_sp) <: Number || error("Some time values not a number!")
    end
    if !(eltype(concvals_sp) <: Union{Number, Missing})
        concvals_sp = identity.(concvals_sp)
        if !(eltype(concvals_sp) <: Union{Number, Missing})
            @warn "Some concentration values not a number, try to fix."
            concvals_sp = floatparse(concvals_sp)
        end
    end
    timevals_sp, concvals_sp
end

"""
    pkimport(data, time, conc, sort; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime())

Import PK data from table `data`.

* `time` - time column;
* `conc` - concentration column;
* `sort` - subject sorting columns.

keywords:

* `kelauto` - if `true` auto range settings, if `false` used `kelstart`/`kelend` from `elimrange`;
* `elimrange` - set elimination range settings;
* `dosetime` - set dose and dose time, by default dosetime = 0, dose is `NaN`.

!!! note

    If time column have non-unique values - last pair time-concentration will be used.

"""
function pkimport(data, time, conc, sort; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime())
    if isa(sort, String) sort = [Symbol(sort)] end
    if isa(sort, Symbol) sort = [sort] end

    Tables.istable(data) || error("Data not a table!")

    cols   = Tables.columns(data)
    cdata  = Tuple(Tables.getcolumn(cols, y) for y in sort)
    d      = Dict{Tuple{eltype.(cdata)...}, Vector{Int}}()
    indsdict!(d, cdata)

    timec = Tables.getcolumn(data, time)
    concc = Tables.getcolumn(data, conc)

    any(isnanormissing, timec) && error("Some time values is NaN or Missing!")

    sdata = Vector{PKSubject}(undef, length(d))
    i = one(Int)
    @inbounds for (k, v) in d
        timevals = timec[v]
        concvals = concc[v]
        if !allunique(timevals)
            @warn "Not all time values is unique, last observation used! ($k)"
            nuv = nonunique(timevals)
            nuvinds = findall(x -> x == first(nuv), timevals)[1:end-1]
            if length(nuv) > 1
                for cnt = 2:length(nuv)
                    append!(nuvinds, findall(x -> x == nuv[cnt], timevals)[1:end-1])
                end
            end
            sort!(nuvinds)
            deleteat!(timevals, nuvinds)
            deleteat!(concvals, nuvinds)
        end
        sp = sortperm(timevals)
        timevals_sp = timevals[sp]
        concvals_sp = concvals[sp]
        timevals_sp, concvals_sp = checkvalues(timevals_sp, concvals_sp)
        sdata[i] = PKSubject(timevals_sp, concvals_sp, kelauto, elimrange,  dosetime, Dict(sort .=> k))
        i += one(Int)
    end
    return DataSet(identity.(sdata))
end
"""
    pkimport(data, time, conc; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime())

Import PK data from tabular data `data`, `time` - time column, `conc` - concentration column.
"""
function pkimport(data, time, conc; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime())
    timevals_sp, concvals_sp = checkvalues(copy(Tables.getcolumn(data, time)), copy(Tables.getcolumn(data, conc)))
    pkimport(timevals_sp, concvals_sp; kelauto = kelauto,  elimrange = elimrange, dosetime = dosetime)
end
"""
    pkimport(time, conc; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime(), id = Dict{Symbol, Any}())

Import PK data from time vector `time` and concentration vector `conc`.
"""
function pkimport(time, conc; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime(), id = Dict{Symbol, Any}())
    timevals_sp, concvals_sp = checkvalues(copy(time), copy(conc))
    PKSubject(timevals_sp, concvals_sp, kelauto, elimrange,  dosetime, id)
end

"""
    upkimport(data, stime, etime, conc, vol, sort; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime())

Import urine PK data from table `data`.

* `stime` - start time column;
* `etime` - end time column;
* `conc` - concentration column;
* `vol` - volume column;
* `sort` - subject sorting columns.
"""
function upkimport(data, stime, etime, conc, vol, sort; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime())
    if isa(sort, String) sort = [Symbol(sort)] end
    if isa(sort, Symbol) sort = [sort] end

    cols   = Tables.columns(data)
    cdata  = Tuple(Tables.getcolumn(cols, y) for y in sort)
    d      = Dict{Tuple{eltype.(cdata)...}, Vector{Int}}()
    indsdict!(d, cdata)

    tnames = Symbol.(names(data))
    isa(stime, Symbol) || stime in tnames || error("column Start Time ($stime) not found")
    isa(etime, Symbol) || etime in tnames || error("column End Time ($etime) not found")
    isa(conc, Symbol) || conc in tnames || error("column Concentration ($conc) not found")
    isa(vol, Symbol) || vol in tnames || error("column Volume ($vol) not found")

    stimec = Tables.getcolumn(data, stime)
    etimec = Tables.getcolumn(data, etime)
    concc = Tables.getcolumn(data, conc)
    volc = Tables.getcolumn(data, vol)

    any(isnanormissing, stimec) && error("Some Start Time values is NaN or Missing!")
    any(isnanormissing, etimec) && error("Some End Time values is NaN or Missing!")

    sdata = Vector{UPKSubject}(undef, length(d))
    i = one(Int)
    @inbounds for (k, v) in d
        stimevals = stimec[v]
        etimevals = etimec[v]
        concvals  = concc[v]
        volvals   = volc[v]
        #=
        timeranges = collect(zip(stimevals, etimevals))
        sp = sortperm(stimevals)

        timevals_sp = timeranges[sp]
        concvals_sp = concvals[sp]
        volvals_sp  = volvals[sp]

        if length(timevals_sp) > 1
            for c = 2:length(timevals_sp)
                timevals_sp[c][1] == timevals_sp[c-1][2] || error("Start time ($(timevals_sp[c][1])) for observation $c not equal End time ($(timevals_sp[c-1][2])) for observation $(c-1)!")
            end
        end
        sdata[i] = UPKSubject(timevals_sp, concvals_sp, volvals_sp, kelauto, elimrange,  dosetime, Dict(sort .=> k))
        =#
        sdata[i] = upkimport(stimevals, etimevals, concvals, volvals; kelauto = kelauto,  elimrange = elimrange, dosetime = dosetime, id = Dict(sort .=> k))
        i += one(Int)
    end
    return DataSet(identity.(sdata))
end
"""
    upkimport(data, stime, etime, conc, vol; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime())

Import single urine PK data from table `data`.

* `stime` - start time column;
* `etime` - end time column;
* `conc` - concentration column;
* `vol` - volume column.
"""
function upkimport(data, stime, etime, conc, vol; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime())
    upkimport(Tables.getcolumn(data, stime), Tables.getcolumn(data, etime), Tables.getcolumn(data, conc), Tables.getcolumn(data, vol); kelauto = kelauto,  elimrange = elimrange, dosetime = dosetime)
end
"""
    upkimport(stime, etime, conc, vol; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime())

Import urine PK data from time vectors:
* `stime` - start times;
* `etime` - end times;
* `conc` - concentrations;
* `vol` - volumes.
"""
function upkimport(stime, etime, conc, vol; kelauto = true,  elimrange = ElimRange(), dosetime = DoseTime(), id = Dict{Symbol, Any}())
    any(isnanormissing, stime) && error("Some Start Time values is NaN or Missing!")
    any(isnanormissing, etime) && error("Some End Time values is NaN or Missing!")
    timeranges = collect(zip(stime, etime))
    sp = sortperm(stime)
    timevals_sp = timeranges[sp]
    concvals_sp = conc[sp]
    volvals_sp  = vol[sp]
    if length(timevals_sp) > 1
        for c = 2:length(timevals_sp)
            timevals_sp[c][1] == timevals_sp[c-1][2] || error("Start time ($(timevals_sp[c][1])) for observation $c not equal End time ($(timevals_sp[c-1][2])) for observation $(c-1)!")
        end
    end
    UPKSubject(timevals_sp, concvals_sp, volvals_sp, kelauto, elimrange,  dosetime, id)
end

function pdimport(data, time, obs, sort; bl = 0, th = NaN)
    if isa(sort, String) sort = [Symbol(sort)] end
    if isa(sort, Symbol) sort = [sort] end

    Tables.istable(data) || error("Data not a table!")

    cols   = Tables.columns(data)
    cdata  = Tuple(Tables.getcolumn(cols, y) for y in sort)
    d      = Dict{Tuple{eltype.(cdata)...}, Vector{Int}}()
    indsdict!(d, cdata)

    timec = Tables.getcolumn(data, time)
    obsc  = Tables.getcolumn(data, obs)

    any(isnanormissing, timec) && error("Some time values is NaN or Missing!")

    sdata = Vector{PDSubject}(undef, length(d))
    i = one(Int)
    @inbounds for (k, v) in d
        timevals = timec[v]
        obsvals  = obsc[v]

        if !allunique(timevals)
            @warn "Not all time values is unique, last observation used! ($k)"
            nuv = nonunique(timevals)
            nuvinds = findall(x -> x == first(nuv), timevals)[1:end-1]
            if length(nuv) > 1
                for cnt = 2:length(nuv)
                    append!(nuvinds, findall(x -> x == nuv[cnt], timevals)[1:end-1])
                end
            end
            sort!(nuvinds)
            deleteat!(timevals, nuvinds)
            deleteat!(obsvals, nuvinds)
        end
        sp = sortperm(timevals)
        timevals_sp = timevals[sp]
        obsvals_sp = obsvals[sp]
        timevals_sp, obsvals_sp = checkvalues(timevals_sp, obsvals_sp)
        sdata[i] = PDSubject(timevals_sp, obsvals_sp, bl, th, Dict(sort .=> k))
        i += one(Int)
    end
    return DataSet(identity.(sdata))
end


text:
    refVals=[]
    refPrecs=[]

# SC ========== Test number 1 reference values and precision match template. =======
# SC ========== /home/jmc/cliMa/cliMa_new_jmc/test/Ocean/SplitExplicit/simple_box_2dt.jl test reference values ======================================
# BEGIN SCPRINT
# varr - reference values (from reference run)
# parr - digits match precision (hand edit as needed)
#
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]
varr = [
 [  "oce Q_3D",   "u[1]", -1.50962361753736729053e-01,  1.73787194589605287209e-01,  4.68098176666432563148e-03,  1.95918702087896703934e-02 ],
 [  "oce Q_3D",   "u[2]", -1.41693277014211693743e-01,  1.46256903927605874660e-01, -2.51527289903215496569e-03,  1.75577307352173875299e-02 ],
 [  "oce Q_3D",       :η, -7.42525448756723882582e-01,  3.69952349991580053956e-01, -2.12300579216294864728e-04,  2.72196784419864945548e-01 ],
 [  "oce Q_3D",       :θ,  4.64352816670868251414e-03,  9.93366845172767831684e+00,  2.49647184360694174288e+00,  2.18037711611593776340e+00 ],
 [   "oce aux",       :w, -1.93642118079277015088e-04,  1.64388494599151930133e-04,  4.30134987337706410355e-07,  1.44731816723112416059e-05 ],
 [   "oce aux",    :pkin, -9.02450734444008029200e-01,  0.00000000000000000000e+00, -3.32591606145584972598e-01,  2.55055561435138855586e-01 ],
 [   "oce aux",     :wz0, -2.13130590203792726805e-05,  3.03764295307291628048e-05, -1.29351489286388951209e-10,  7.79133385457225507325e-06 ],
 [   "oce aux", "u_d[1]", -1.47678995808530610923e-01,  1.16175707854701368293e-01, -3.40162514353149251730e-05,  1.18326667626959050605e-02 ],
 [   "oce aux", "u_d[2]", -1.41297097563442025647e-01,  1.40410041592192808002e-01, -8.66960554174275520279e-06,  1.16477781860160160832e-02 ],
 [   "oce aux", "ΔGu[1]", -1.82462180035313031949e-06,  1.75730579970564754281e-06, -3.58905562801251356298e-08,  2.63341240419211030793e-07 ],
 [   "oce aux", "ΔGu[2]", -1.23605953229188910373e-06,  2.21985915696832233044e-06,  1.31099407764425287042e-06,  6.56785741293597450970e-07 ],
 [   "oce aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.15573163901915703900e+06 ],
 [ "baro Q_2D",   "U[1]", -1.28633175515245561371e+01,  6.10318317172543416405e+01,  4.70740761170336963204e+00,  1.49615297140118084229e+01 ],
 [ "baro Q_2D",   "U[2]", -2.82520952983562132488e+01,  6.47619508637900196391e+01, -2.47935794730033354227e+00,  1.34335968079767127392e+01 ],
 [ "baro Q_2D",       :η, -7.42731313350159738640e-01,  3.69996965336395311486e-01, -2.12333995604797106261e-04,  2.72347141148035598590e-01 ],
 [  "baro aux",  "Gᵁ[1]", -1.75730579970564764446e-03,  1.82462180035313035337e-03,  3.58905562801251301771e-05,  2.63354276791606873361e-04 ],
 [  "baro aux",  "Gᵁ[2]", -2.21985915696832226268e-03,  1.23605953229188902411e-03, -1.31099407764425279418e-03,  6.56818254634440362900e-04 ],
 [  "baro aux", "U_c[1]", -1.28413474945461345555e+01,  6.09612852044860531464e+01,  4.71376723747856551938e+00,  1.49474605677143106419e+01 ],
 [  "baro aux", "U_c[2]", -2.82757649462843900778e+01,  6.47162674795010843809e+01, -2.50650840410464237351e+00,  1.34316239381179141077e+01 ],
 [  "baro aux",     :η_c, -7.42525448756723882582e-01,  3.69952349991580053956e-01, -2.12300579216321237946e-04,  2.72210259174677282612e-01 ],
 [  "baro aux", "U_s[1]", -1.28415899595582310155e+01,  6.09623267104856623178e+01,  4.71378782923472794408e+00,  1.49476602108240328448e+01 ],
 [  "baro aux", "U_s[2]", -2.82766214528259105521e+01,  6.47161195504508555132e+01, -2.50657419606653197874e+00,  1.34316536525727965312e+01 ],
 [  "baro aux",     :η_s, -7.42533259412902157948e-01,  3.69956621986456846152e-01, -2.12295179177931456769e-04,  2.72211064871600239012e-01 ],
 [  "baro aux",  "Δu[1]", -4.52584773336763403293e-04,  4.01399962579148425305e-04, -1.17814846783320024927e-05,  1.01681656567342138366e-04 ],
 [  "baro aux",  "Δu[2]", -2.43941531101083075800e-04,  3.86401032853338948295e-04,  5.11421508443407068253e-05,  7.52601602185662730366e-05 ],
 [  "baro aux",  :η_diag, -7.39952355750089774133e-01,  3.68437192305505589740e-01, -2.11882685629188936925e-04,  2.71995932052603417439e-01 ],
 [  "baro aux",      :Δη, -3.57682889898169875664e-03,  4.28352100974604965700e-03, -4.17893587097645250512e-07,  1.12668312540401493148e-03 ],
 [  "baro aux",       :y,  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.15578885204060329124e+06 ],
]
parr = [
 [  "oce Q_3D",   "u[1]",    12,    12,    12,    12 ],
 [  "oce Q_3D",   "u[2]",    12,    12,    12,    12 ],
 [  "oce Q_3D",       :η,    12,    12,     8,    12 ],
 [  "oce Q_3D",       :θ,    12,    12,    12,    12 ],
 [   "oce aux",       :w,    12,    12,     8,    12 ],
 [   "oce aux",    :pkin,    12,    12,    12,    12 ],
 [   "oce aux",     :wz0,    12,    12,     8,    12 ],
 [   "oce aux", "u_d[1]",    12,    12,    12,    12 ],
 [   "oce aux", "u_d[2]",    12,    12,    12,    12 ],
 [   "oce aux", "ΔGu[1]",    12,    12,    12,    12 ],
 [   "oce aux", "ΔGu[2]",    12,    12,    12,    12 ],
 [   "oce aux",       :y,    12,    12,    12,    12 ],
 [ "baro Q_2D",   "U[1]",    12,    12,    12,    12 ],
 [ "baro Q_2D",   "U[2]",    12,    12,    12,    12 ],
 [ "baro Q_2D",       :η,    12,    12,     8,    12 ],
 [  "baro aux",  "Gᵁ[1]",    12,    12,    12,    12 ],
 [  "baro aux",  "Gᵁ[2]",    12,    12,    12,    12 ],
 [  "baro aux", "U_c[1]",    12,    12,    12,    12 ],
 [  "baro aux", "U_c[2]",    12,    12,    12,    12 ],
 [  "baro aux",     :η_c,    12,    12,     8,    12 ],
 [  "baro aux", "U_s[1]",    12,    12,    12,    12 ],
 [  "baro aux", "U_s[2]",    12,    12,    12,    12 ],
 [  "baro aux",     :η_s,    12,    12,     8,    12 ],
 [  "baro aux",  "Δu[1]",    12,    12,    12,    12 ],
 [  "baro aux",  "Δu[2]",    12,    12,    12,    12 ],
 [  "baro aux",  :η_diag,    12,    12,     8,    12 ],
 [  "baro aux",      :Δη,    12,    12,     8,    12 ],
 [  "baro aux",       :y,    12,    12,    12,    12 ],
]
# END SCPRINT
# SC ====================================================================================

    append!(refVals ,[ varr ] )
    append!(refPrecs,[ parr ] )

#! format: on


text:
function showinfo(tree::Tree)
    i = 1
    for n in values(tree.lnodes)
        println("Node $i: $(n.label)")
        n.isroot ? println("Root") : println("Ancestor: $(n.anc.label)")
        if n.isleaf
            println("Leaf")
        else
            print("Children: ")
            for c in n.child
                print(" $(c.label), ")
            end
            println()
        end
        i+=1
        println()
    end
end

"""
    show(io::IO, tree::Tree, maxnodes=40; kwargs...)
    show(t::Tree, maxnodes=40; kwargs...)
"""
function Base.show(io::IO, tree::Tree, maxnodes=40; kwargs...)
    if length(tree.lnodes) < maxnodes
        print_tree(io, tree; kwargs...)
    end
end
Base.show(t::Tree, maxnodes=40; kwargs...) = show(stdout, t, maxnodes; kwargs...)

function Base.show(io::IO, n::TreeNode)
    if !get(io, :compact, false)
        nodeinfo(io, n)
    end
end
Base.show(n::TreeNode) = show(stdout, n)

"""
    nodeinfo(io, node)

Print information about `node`.
"""
function nodeinfo(io, node)
    println(io, "Node $(node.label): ")
    node.isroot ? println(io, "Ancestor : none (root)") : println(io, "Ancestor: $(node.anc.label), tau = $(node.tau)")
    println(io, "$(length(node.child)) children: $([x.label for x in node.child])")
end

function print_tree_(io, node, cdepth; vindent=2, hindent=5, hoffset=0, maxdepth=5)
    hspace = ""
    for i in 1:hindent
        hspace *= "-"
    end
    offset = ""
    for i in 1:hoffset
        offset *= " "
    end
    if cdepth < maxdepth
    	println(io, "$offset $hspace $(node.label):$(node.tau)")
    elseif cdepth == maxdepth
    	if node.isleaf
    		println(io, "$offset $hspace $(node.label):$(node.tau)")
    	else
    		println(io, "$offset $hspace $(node.label):$(node.tau) ...")
    	end
    end
        #
    if cdepth <= maxdepth
        if !node.isleaf
            for c in node.child
                for i in 1:vindent
                    cdepth < maxdepth && println(io, "$offset $(" "^hindent)|")
                end
                print_tree_(io, c, cdepth + 1, vindent=vindent, hindent=hindent, hoffset=hoffset+hindent, maxdepth=maxdepth)
            end
        end
        #
    end
end
function print_tree(io, node::TreeNode; vindent=2, hindent=5, maxdepth=5)
    print_tree_(io, node, 1, vindent=vindent, hindent=hindent, hoffset=0, maxdepth=maxdepth)
end
print_tree(io, t::Tree; vindent=2, hindent=5, maxdepth=5) = print_tree(io, t.root; vindent=2, hindent=5, maxdepth=maxdepth)

"""
    check_tree(t::Tree; strict=true)

- Every non-leaf node should have at least one child (two if `strict`)
- Every non-root node should have exactly one ancestor
- If n.child[...] == c, c.anc == n is true
- Tree has only one root
"""
function check_tree(tree::Tree; strict=true)
    labellist = Dict{String, Int64}()
    nroot = 0
    flag = true
    for n in values(tree.lnodes)
        if !n.isleaf && length(n.child)==0
        	(flag = false) || (@warn "Node $(n.label) is non-leaf and has no child.")
        elseif !n.isroot && n.anc == nothing
        	(flag = false) || (@warn "Node $(n.label) is non-root and has no ancestor.")
        elseif strict && length(n.child) == 1
            if !(n.isroot && n.child[1].isleaf)
        	   (flag = false) || (@warn "Node $(n.label) has only one child.")
            end
        elseif length(n.child) == 0 && !haskey(tree.lleaves, n.label)
            (flag = false) || (@warn "Node $(n.label) has no child but is not in `tree.lleaves`")
        end
        for c in n.child
            if c.anc != n
                (flag = false) || (@warn "One child of $(n.label) does not satisfy `c.anc == n`.")
            end
        end
        if get(labellist, n.label, 0) == 0
            labellist[n.label] = 1
        else
            labellist[n.label] += 1
            (flag = false) || (@warn "Label $(n.label) already exists!")
        end
        if n.isroot
            nroot += 1
        end
    end
    if nroot > 1
        (flag = false) || (@warn "Tree has multiple roots")
    elseif nroot ==0
        (flag = false) || (@warn "Tree has no root")
    end
    return flag
end

make_random_label(base="NODE") = make_random_label(base, 8)
make_random_label(base, i) = "$(base)_$(randstring(i))"

"""
	label_nodes!(t::Tree)
"""
function label_nodes!(t::Tree)
	for n in internals(t)
		if isempty(n.label)
			n.label = make_random_label()
		end
	end
end
"""
    create_label(t::Tree, base="NODE")

Create new node label in tree `t` with format `\$(base)_i` with `i::Int`.
"""
function create_label(t::Tree, base="NODE")
    label_init = 1
    pattern = Regex(base)
    for n in values(t.lnodes)
        if match(pattern, n.label)!=nothing && parse(Int64, n.label[length(base)+2:end]) >= label_init
            label_init = parse(Int64, n.label[length(base)+2:end]) + 1
        end
    end
    return "$(base)_$(label_init)"
end

function set_unique_label!(node::TreeNode, t::Tree; delim = "__")
	id = randstring(8)
	node.label *= delim * id
	while haskey(t.lnodes, node.label)
		node.label = node.label[1:end-5] * randstring(8)
	end
end

"""
    map_dict_to_tree!(t::Tree{MiscData}, dat::Dict; symbol=false, key = nothing)

Map data in `dat` to nodes of `t`. All node labels of `t` should be keys of `dat`. Entries of `dat` should be dictionaries, or iterable similarly, and are added to `n.data.dat`.

If `!isnothing(key)`, only a specific key of `dat` is added. It's checked for by `k == key || Symbol(k) == key` for all keys `k` of `dat`.
If `symbol`, data is added to nodes of `t` with symbols as keys.
"""
function map_dict_to_tree!(t::Tree{MiscData}, dat::Dict; symbol=false, key = nothing)
    for (name, n) in t.lnodes
        for (k,v) in dat[name]
            if !isnothing(key) && (k == key || Symbol(k) == key)
                n.data.dat[key] = v
            elseif isnothing(key)
                n.data.dat[symbol ? Symbol(k) : k] = v
            end
        end
    end
    nothing
end
"""
    map_dict_to_tree!(t::Tree{MiscData}, dat::Dict, key)

Map data in `dat` to nodes of `t`. All node labels of `t` should be keys of `dat`. Entries of `dat` corresponding to `k` are added to `t.lnodes[k].data.dat[key]`.
"""
function map_dict_to_tree!(t::Tree{MiscData}, dat::Dict, key)
    for (name, n) in t.lnodes
        n.data.dat[key] = dat[name]
    end
end

"""
    rand_times!(t, p)

Add random branch lengths to tree.
"""
function rand_times!(t, p)
    for n in values(t.lnodes)
        if !n.isroot
            n.tau = rand(p)
        end
    end
end

##################################################
##### Utilities for dicts with nested keys #######
##################################################
function recursive_key_init!(dat, key, ks...)
    if !haskey(dat, key)
        dat[key] = Dict()
    end
    recursive_key_init!(dat[key], ks...)
end
recursive_key_init!(dat) = nothing
function recursive_get(dat, key, ks...)
    if isempty(ks)
        return dat[key]
    else
        return recursive_get(dat[key], ks...)
    end
end
recursive_get(dat, key::Tuple) = recursive_get(dat, key...)
function recursive_set!(dat, value, key, ks...)
    if isempty(ks)
        dat[key] = value
    else
        if !haskey(dat, key)
            dat[key] = Dict()
        end
        recursive_set!(dat[key], value, ks...)
    end
    dat
end
recursive_set!(dat, value, key::Tuple) = recursive_set!(dat, value, key...)
function recursive_push!(dat, value, key, ks...)
    if isempty(ks)
        push!(dat[key], value)
    else
        recursive_push!(dat[key], value, ks...)
    end
    dat
end
recursive_push!(dat, value, key::Tuple) = recursive_push!(dat, value, key...)
function recursive_haskey(dat, key, ks...)
    if !haskey(dat, key)
        return false
    elseif isempty(ks)
        return true
    else
        return recursive_haskey(dat[key], ks...)
    end
end

max_stars_repo_name: PierreBarrat/TreeTools

max_issues_repo_name: PierreBarrat/TreeTools

max_forks_repo_name: PierreBarrat/TreeTools


text:
using SingularSpectrumAnalysis
using Test, Statistics, Random
const SSA = SingularSpectrumAnalysis

@testset "SingularSpectrumAnalysis" begin
    Random.seed!(1)
    # Generate some data
    L = 20
    K = 100
    N = K*L;
    t = 1:N;
    T = 20;
    y = sin.(2pi/T*t);
    y .+= (0.5sin.(2pi/T*4*t)).^2
    e = 0.1randn(N);
    yn = y+e;
    # plot(yn)

    yt, ys = analyze(yn, L)
    A,x = fit_trend(yt, 1)
    @test x ≈ [0.124, 0] atol=0.1

    USV = hsvd(yn,L)

    trend, seasonal_groupings = autogroup(USV)
    @test trend == 3
    @test seasonal_groupings == [1:2, 4:5]

    yrt, yrs = reconstruct(USV, trend, seasonal_groupings)
    yr = sum([yrt yrs],dims=2)
    @test sqrt(mean((y.-yn).^2)) > sqrt(mean((y.-yr).^2))

    using Plots
    plot(USV, cumulative=true)
    plot(USV, cumulative=false)
    pairplot(USV, seasonal_groupings)
end
@testset "forecasting" begin
    include("forecast.jl")
end


text:
@testset "Identity dynamical system in square" begin
    include("identity.jl")
end
@testset "Linear dynamical system in square" begin
    include("lti.jl")
end
@testset "Controlled linear dynamical system in square" begin
    include("controlled.jl")
end
include("cis.jl")


text:
module Fabrics

# exports
export ParameterizationWeak
export ParameterizationStrong
export Parameterization
export parameterization
export fabric_shape
export fabric_orientation
export weakening
export Weakening

# Types and structures
abstract type Parameterization end

struct ParameterizationWeak{T} <: Parameterization
    ζ::T # Fabric shape
    α::T # Fabric orientation
end

struct Weakening{T} <: Parameterization
    ζ::T
    ξ::T
    χ::T
    θ::T
    ψ::T
    λ::T
end

struct ParameterizationStrong{T} <: Parameterization
    ζ::T # Fabric shape
    ξ::T # Fabric shape
    χ::T # Fabric shape
    θ::T # Fabric shape
    α::T # Fabric orientation
end

# Fabric shape (weak inclusions)
function fabric_shape(a1, a2, p::ParameterizationWeak{T}) where T
    A = log10(a1/a2) # Bulk semi axes ratios
    r₁ = p.ζ*A  # inclusion  log10(a1/a2)
    r₂ = A # inclusion  log10(a2/a3)
    return r₁, r₂
end

# Fabric shape (weak inclusions)
function fabric_shape(a1, a2, p::NamedTuple)
    A = log10(a1/a2) # Bulk semi axes ratios
    R1, R2 = p.R1, p.R2
    r₁ = R1.ζ + R1.ξ*A + R1.χ*A^2 + R1.θ*A^3 # inclusion  log10(a1/a2)
    r₂ = R2.ζ + R2.ξ*A + R2.χ*A^2 + R2.θ*A^3 # inclusion  log10(a2/a3)
    return r₁, r₂
end

# Fabric orientation: weak inclusions
fabric_orientation(α,  p::T) where T<:Parameterization = p.α*α

# Fabric orientation: strong inclusions
fabric_orientation(α,  p::NamedTuple) = p.R1.α*α

# Proportionality constanst of shape and orientation for 2-phase Newtonian aggregate
# with weak inclusion. Numbers from de Montserrat et al. 2021, JGR, Table 1
parameterization(Δη::T, ϕ::T) where T<:Real = _barrier(Δη, ϕ)

_barrier(Δη, ϕ) = parameterization(Val{Δη}(), Val{ϕ}())

parameterization(::Val{1e-3}, ::Val{0.1}) = ParameterizationWeak(
    0.687903 ,
    1.000936
)

parameterization(::Val{1e-3}, ::Val{0.2}) = ParameterizationWeak(
    0.732392 ,
    0.933471
)

parameterization(::Val{1e-3}, ::Val{0.3}) = ParameterizationWeak(
    0.709184 ,
    1.095165
)

parameterization(::Val{1e-2}, ::Val{0.1}) = ParameterizationWeak(
    0.715399 ,
    1.
)

parameterization(::Val{1e-2}, ::Val{0.2}) = ParameterizationWeak(
    0.774659 ,
    1.002619
)

parameterization(::Val{1e-2}, ::Val{0.3}) = ParameterizationWeak(
    0.704054 ,
    0.945960
)

parameterization(::Val{1e-1}, ::Val{0.1}) = ParameterizationWeak(
    0.793576 ,
    1.07913
)

parameterization(::Val{1e-1}, ::Val{0.3}) = ParameterizationWeak(
    0.752338 ,
    0.983557
)

# Proportionality constanst of shape and orientation for 2-phase Newtonian aggregate
# with strong inclusions. Numbers from de Montserrat et al. 2021, JGR, Table 2
function parameterization(::Val{1e1}, ::Val{0.1})
    R1 = ParameterizationStrong(
        -0.006131,  0.458978, -0.125306, -0.083918, 1.0
    )
    R2 = ParameterizationStrong(
        0.000512, -0.319324 , 0.185224 , -0.306101, 1.0
    )
    return (R1=R1, R2=R2)
end

function parameterization(::Val{1e1}, ::Val{0.2})
    R1 = ParameterizationStrong(
        -0.000913 , 0.389964,-0.056822,0.1010262, 1.020180
    )
    R2 = ParameterizationStrong(
        0.002493 , 0.296955, 0.271349, -0.2545942, 1.020180
    )
    return (R1=R1, R2=R2)
end

function parameterization(::Val{1e1}, ::Val{0.3})
    R1 = ParameterizationStrong(
        0.000295 , 0.372009,-0.070759, 0.2408665, 1.043904
    )
    R2 = ParameterizationStrong(
        0.005547 , 0.266632, 0.271325, -0.1762204, 1.043904
    )
    return (R1=R1, R2=R2)
end

# Viscosity cut-off for 2-phase Newtonian aggregates with weak inclusion.
# Numbers from de Montserrat et al. 2021, JGR, Table 1
weakening(Δη::T, ϕ::T) where T<:Real = weakening(Δη, _barrier_w(ϕ))

_barrier_w(ϕ) = weakening(Val{ϕ}())

weakening(::Val{0.1}) = Weakening(
    0.903221, 1.20365, -1.148182, 0.3735, 0.698558, 0.6619,
)

weakening(::Val{0.2}) = Weakening(
    43.646559, -42.907958, 0.0, 0.261095, 0.984129, 1.0,
)

weakening(::Val{0.3}) = Weakening(
    40.628095, -39.823767, 0.0, 0.196188, 0.985677, 1.0,
)

weakening(Δη, w::Weakening) =  w.ζ*Δη^w.ψ + w.ξ*Δη^w.λ + w.χ*Δη + w.θ

end # module


text:
# Tools for implementing device functionality

export @wrap

# which Julia types map to a given LLVM type
const llvmtypes = Dict{Type,Symbol}(
    Nothing => :void,
    Bool    => :i1,
    Int8    => :i8,
    Int16   => :i16,
    Int32   => :i32,
    Int64   => :i64,
    UInt8   => :i8,
    UInt16  => :i16,
    UInt32  => :i32,
    UInt64  => :i64,
    Float32 => :float,
    Float64 => :double,
)

# which LLVM types map to a given Julia type
const jltypes = Dict{Symbol,Type}(
    :void   => Nothing,
    :i1     => Bool,
    :i8     => Int8,
    :i16    => Int16,
    :i32    => Int32,
    :i64    => Int64,
    :float  => Float32,
    :double => Float64
)

# which Julia types map to a given ROC lib type
const fntypes = Dict{Type,Symbol}(
    Nothing => :void,
    Bool    => :bool,
    Int8    => :i8,
    Int16   => :i16,
    Int32   => :i32,
    Int64   => :i64,
    UInt8   => :u8,
    UInt16  => :u16,
    UInt32  => :u32,
    UInt64  => :u64,
    Float16 => :f16,
    Float32 => :f32,
    Float64 => :f64
)

# Decode an expression of the form:
#
#    function(arg::arg_type, arg::arg_type, ... arg::arg_type)::return_type
#
# Returns a tuple containing the function name, a vector of argument, a vector of argument
# types and the return type (all in symbolic form).
function decode_call(e)
    @assert e.head == :(::)

    # decode the return type expression: single symbol (the LLVM type), or a tuple of 2
    # symbols (the LLVM and corresponding Julia type)
    retspec = e.args[2]
    if isa(retspec, Symbol)
        rettype = retspec
    else
        @assert retspec.head == :tuple
        @assert length(retspec.args) == 2
        rettype = (retspec.args[1], retspec.args[2])
    end

    call = e.args[1]
    @assert call.head == :call

    fn = Symbol(call.args[1])
    args = Symbol[arg.args[1] for arg in call.args[2:end]]
    argtypes = Symbol[arg.args[2] for arg in call.args[2:end]]

    return fn, args, argtypes, rettype
end

# Generate a `llvmcall` statement calling an intrinsic specified as follows:
#
#     intrinsic(arg::arg_type, arg::arg_type, ... arg::arg_type)::return_type [attr]
#
# The argument types should be valid LLVM type identifiers (eg. i32, float, double).
# Conversions to the corresponding Julia type are automatically generated; make sure the
# actual arguments are of the same type to make these conversions no-ops. The optional
# argument `attr` indicates which LLVM function attributes (such as `readnone` or `nounwind`)
# to add to the intrinsic declaration.

# For example, the following call:
#     `@wrap __some_intrinsic(x::float, y::double)::float`
#
# will yield the following `llvmcall`:
# ```
#     Base.llvmcall(("declare float @__somme__intr(float, double)",
#                    "%3 = call float @__somme__intr(float %0, double %1)
#                     ret float %3"),
#                   Float32, Tuple{Float32,Float64},
#                   convert(Float32,x), convert(Float64,y))
# ```
macro wrap(call, attrs="")
    intrinsic, args, argtypes, rettype = decode_call(call)

    # decide on intrinsic return type
    if isa(rettype, Symbol)
        # only LLVM return type specified, match against known LLVM/Julia type combinations
        llvm_ret_typ = rettype
        julia_ret_typ = jltypes[rettype]
    else
        # both specified (for when there is a mismatch, eg. i32 -> UInt32)
        llvm_ret_typ = rettype[1]
        julia_ret_typ = rettype[2]
    end

    llvm_args = String["%$i" for i in 0:length(argtypes)]
    if llvm_ret_typ == :void
        llvm_ret_asgn = ""
        llvm_ret = "void"
    else
        llvm_ret_var = "%$(length(argtypes)+1)"
        llvm_ret_asgn = "$llvm_ret_var = "
        llvm_ret = "$llvm_ret_typ $llvm_ret_var"
    end
    llvm_declargs = join(argtypes, ", ")
    llvm_defargs = join(("$t $arg" for (t,arg) in zip(argtypes, llvm_args)), ", ")

    julia_argtypes = (jltypes[t] for t in argtypes)
    julia_args = (:(convert($argtype, $(esc(arg)))) for (arg, argtype) in zip(args, julia_argtypes))

    dest = ("""declare $llvm_ret_typ @$intrinsic($llvm_declargs)""",
            """$llvm_ret_asgn call $llvm_ret_typ @$intrinsic($llvm_defargs)
                ret $llvm_ret""")

    return quote
        Base.llvmcall($dest, $julia_ret_typ, Tuple{$(julia_argtypes...)}, $(julia_args...))
    end
end

# calculate the size of an LLVM type
llvmsize(::LLVM.LLVMHalf) = sizeof(Float16)
llvmsize(::LLVM.LLVMFloat) = sizeof(Float32)
llvmsize(::LLVM.LLVMDouble) = sizeof(Float64)
llvmsize(::LLVM.IntegerType) = Context() do ctx
    div(Int(intwidth(GenericValue(LLVM.Int128Type(ctx), -1))), 8)
end
llvmsize(ty::LLVM.ArrayType) = length(ty)*llvmsize(eltype(ty))
llvmsize(ty::LLVM.StructType) = ispacked(ty) ? sum(llvmsize(elem) for elem in elements(ty)) : 8*length(elements(ty)) # FIXME: Properly determine non-packed sizing
llvmsize(ty::LLVM.PointerType) = div(Sys.WORD_SIZE, 8)
llvmsize(ty::LLVM.VectorType) = size(ty)
llvmsize(ty) = error("Unknown size for type: $ty, typeof: $(typeof(ty))")

@generated function string_length(ex::Union{Ptr,LLVMPtr})
    Context() do ctx
        T_ex = convert(LLVMType, ex; ctx)
        T_ex_ptr = LLVM.PointerType(T_ex)
        T_i8_ptr = LLVM.PointerType(LLVM.Int8Type(ctx))
        T_i64 = LLVM.Int64Type(ctx)
        llvm_f, _ = create_function(T_i64, [T_ex])
        mod = LLVM.parent(llvm_f)
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            check = BasicBlock(llvm_f, "check"; ctx)
            done = BasicBlock(llvm_f, "done"; ctx)

            position!(builder, entry)
            init_offset = ConstantInt(0; ctx)
            input_ptr = if T_ex isa LLVM.PointerType
                parameters(llvm_f)[1]
            else
                inttoptr!(builder, parameters(llvm_f)[1], T_ex_ptr)
            end
            if LLVM.addrspace(llvmtype(input_ptr)) != LLVM.addrspace(T_ex_ptr)
                input_ptr = addrspacecast!(builder, input_ptr, T_ex_ptr)
            end
            input_ptr = bitcast!(builder, input_ptr, T_i8_ptr)
            br!(builder, check)

            position!(builder, check)
            offset = phi!(builder, T_i64)
            next_offset = add!(builder, offset, ConstantInt(1; ctx))
            append!(LLVM.incoming(offset), [(init_offset, entry), (next_offset, check)])
            ptr = gep!(builder, input_ptr, [offset])
            value = load!(builder, ptr)
            cond = icmp!(builder, LLVM.API.LLVMIntEQ, value, ConstantInt(0x0; ctx))
            br!(builder, cond, done, check)

            position!(builder, done)
            ret!(builder, offset)
        end
        call_function(llvm_f, Csize_t, Tuple{ex}, :ex)
    end
end


text:
#Import Packages needed for simulation
using DifferentialEquations, Sundials #For implementing Differential equations
using LinearAlgebra, SparseArrays, Distributions, Statistics, Random #Linear Algebra and Statistics
using CSV, DataFrames #Data handling
using StatsPlots #For graphing
using JLD2, FileIO #Saving simulations

###############################################################
# 1. Define all the constants that will not change
###############################################################

#Constants for all cells
const N=200 #number of grid points along one dimensions
const nCells = N^2 #number of cells in the simulation
const cellVol = 3e-12 #Cell Volume (liters)
const Na = 6.02e23 #Avagadro's number
const species = 14 #Number of states within each cell (including virus)
const moi = 1.0e-2 #Multicity of infection
const cellIndicies = CartesianIndices(zeros(N,N)) #set of indices for looping through every cell

#Function that converts molecules to nM
m2c(molecule) = @. 1e9*molecule/(cellVol*Na)

#Constants for all simulations
const tspan = (0.0,48.0) #Time span for simulation
const ΔIFNβ = zeros(N,N) #Define memory space to hold the Laplacian
const statesNames = ["cGAS","DNA","Sting","cGAMP","IRF3","IFNbm","IFNb","STAT",
                     "SOCSm","IRF7m","TREX1m","IRF7","TREX1","Virus"] #for plotting
const θNames = [:k1f, :k1r, :k3f, :k3r, :k4f, :kcat5, :Km5, :k5r, :kcat6, :Km6, :kcat7,
 :Km7, :kcat8, :Km8, :k8f, :k9f, :k10f1, :k10f2, :k11f, :k12f, :k13f, :k6f, :kcat2,
 :Km2, :τ4, :τ6, :τ7, :τ8, :τ9, :τ10, :τ11, :τ12, :τ13, :k14f,:τ14] #Parameter names
 #These species are not starting at zero
 const nonZeroSpeciesIdx = [1,3,5] #cGAS, Sting, IRF3
 const nonZeroSpeciesValues = m2c([1e3, 1e3, 1e4]) #convert to concentration

#This function modifies initial values by adding guassian noise
#Standard deviation is same magnitude as mean
# If p=0.5 then: (u₁=0.1 μ₂=100) → (σ₁=0.05 σ₂=50) instead of (σ₁=0.5 σ₂=0.5)
AddNoise2States(p) = [truncated(Normal(μ,p*μ),0,Inf) for μ in nonZeroSpeciesValues]

 #Often it is useful to pass parameters between functions during the ODE solve,
 #This struct hold all the parameters that we need to keep track of
 mutable struct ParContainer{T}
   par::T #Rate Constants
   mass::Vector{Array{Float64,2}} #Mass balances
   deathParameter::Array{Float64,2} # 0 or 1 indicating cell is dead
   DNAReplicate::Int64 #0 if ISD, 1 if Virus
   cellsInfected::Array{Float64,2} #Time cell was infected
   cellsDead::Array{Float64,2} #Time cell was killed
   infectFirstAttempt::BitArray{2} #Has cell tried to infect neighbors?
 end

 ###############################################################
 # 2. Helper function that discretizes the Laplacian
 ###############################################################
 #∇u is updated with Laplacian estimate from u
 function ∇²(Δu,u)
   #Get dimensions of the input and define some constants
     n1, n2 = size(u)
     Δx = 32.0 #Grid spacing (diameter of cell in μm)
     D=97.5*3600.0 #Diffusion coefficient (μm^2/hr)
     h = D/Δx^2

     # internal nodes
     for j = 2:n2-1
         for i = 2:n1-1
             @inbounds  Δu[i,j] = h*(u[i+1,j] + u[i-1,j] + u[i,j+1] + u[i,j-1] - 4*u[i,j])
         end
     end

     # left/right edges
     for i = 2:n1-1
         @inbounds Δu[i,1] = h*(u[i+1,1] + u[i-1,1] + 2*u[i,2] - 4*u[i,1])
         @inbounds Δu[i,n2] = h*(u[i+1,n2] + u[i-1,n2] + 2*u[i,n2-1] - 4*u[i,n2])
     end

     # top/bottom edges
     for j = 2:n2-1
         @inbounds Δu[1,j] = h*(u[1,j+1] + u[1,j-1] + 2*u[2,j] - 4*u[1,j])
         @inbounds Δu[n1,j] = h*(u[n1,j+1] + u[n1,j-1] + 2*u[n1-1,j] - 4*u[n1,j])
     end

     # corners
     @inbounds Δu[1,1] = h*(2*(u[2,1] + u[1,2]) - 4*u[1,1])
     @inbounds Δu[n1,1] = h*(2*(u[n1-1,1] + u[n1,2]) - 4*u[n1,1])
     @inbounds Δu[1,n2] = h*(2*(u[2,n2] + u[1,n2-1]) - 4*u[1,n2])
     @inbounds Δu[n1,n2] = h*(2*(u[n1-1,n2] + u[n1,n2-1]) - 4*u[n1,n2])
 end

 ###############################################################
 # 3. This is the ODE model for the cGAS pathway
 ###############################################################

 # Define the discretized PDE as an ODE function
 function Model!(du,u,p,t)

   #Species
   cGAS = @view u[:,:,1]
   DNA = @view u[:,:,2]
   Sting = @view u[:,:,3]
   cGAMP = @view u[:,:,4]
   IRF3 = @view u[:,:,5]
   IFNβm = @view u[:,:,6]
   IFNβ = @view u[:,:,7]
   STAT = @view u[:,:,8]
   SOCSm = @view u[:,:,9]
   IRF7m = @view u[:,:,10]
   TREX1m = @view u[:,:,11]
   IRF7 = @view u[:,:,12]
   TREX1 = @view u[:,:,13]
   Virus = @view u[:,:,14]

   #Derivatives
   d_cGAS = @view du[:,:,1]
   d_DNA = @view du[:,:,2]
   d_Sting = @view du[:,:,3]
   d_cGAMP = @view du[:,:,4]
   d_IRF3 = @view du[:,:,5]
   d_IFNβm = @view du[:,:,6]
   d_IFNβ = @view du[:,:,7]
   d_STAT = @view du[:,:,8]
   d_SOCSm = @view du[:,:,9]
   d_IRF7m = @view du[:,:,10]
   d_TREX1m = @view du[:,:,11]
   d_IRF7 = @view du[:,:,12]
   d_TREX1 = @view du[:,:,13]
   d_Virus = @view du[:,:,14]

   #Parameters
   k1f, k1r, k3f, k3r, k4f, kcat5, Km5, k5r, kcat6, Km6, kcat7, Km7, kcat8, Km8, k8f, k9f, k10f1, k10f2, k11f, k12f, k13f, k6f, kcat2, Km2, τ4, τ6, τ7, τ8, τ9, τ10, τ11, τ12, τ13, k14f, τ14 = p.par
   #Constants from the mass balances
   cGAStot, Stingtot, IRF3tot = p.mass
   #Which cells are dead?
   💀 = p.deathParameter
   #Should DNA be allowed to replicate (only with virus, not with ISD)
   🔁 = p.DNAReplicate

   #Calculate the diffusion of IFNβ
   ∇²(ΔIFNβ,IFNβ)

   #Update derivatives for each species according to model
   @. d_cGAS = -k1f*cGAS*DNA + k1r*(cGAStot - cGAS)
   @. d_DNA = -k1f*cGAS*DNA + k1r*(cGAStot - cGAS) - kcat2*TREX1*DNA / (Km2 + DNA) + 🔁*💀*DNA*(0.55-DNA)/0.55
   @. d_Sting = -k3f*cGAMP*Sting + k3r*(Stingtot - Sting)
   @. d_cGAMP = k4f*(cGAStot - cGAS) - k3f*cGAMP*Sting + k3f*(Stingtot - Sting) - τ4*cGAMP
   @. d_IRF3 = -kcat5*IRF3*(Stingtot - Sting) / (Km5 +IRF3) + k5r*(IRF3tot - IRF3)
   @. d_IFNβm = 💀*kcat6*(IRF3tot - IRF3) / (Km6 + (IRF3tot - IRF3)) + 💀*k6f*IRF7 - τ6*IFNβm
   @. d_IFNβ = 💀*kcat7*IFNβm / (Km7 + IFNβm) - τ7*IFNβ + ΔIFNβ #Add the diffusion in here
   @. d_STAT = 💀*kcat8*IFNβ / (Km8 + IFNβ) * 1.0/(1.0+k8f*SOCSm) - τ8*STAT
   @. d_SOCSm = 💀*k9f*STAT - τ9*SOCSm
   @. d_IRF7m = 💀*k10f1*STAT + 💀*k10f2*IRF7 - τ10*IRF7m
   @. d_TREX1m = 💀*k11f*STAT - τ11*TREX1m
   @. d_IRF7 = 💀*k12f*IRF7m - τ12*IRF7
   @. d_TREX1 = 💀*k13f*TREX1m - τ13*TREX1
   @. d_Virus = 💀*k14f*DNA - τ14*Virus
 end

 ###############################################################
# 4. Set up function that will return an ODE problem to solve
 ###############################################################

function ModelSetup(infectionMethod,IFNStoch,Hetero)
    #Parameter values for the ODEs
    θVals = [2.6899, 4.8505, 0.0356, 7.487, 517.4056, 22328.3852, 11226.3682,0.9341,
             206.9446, 10305.461, 47639.70295,3.8474, 13.006, 78.2048, 0.0209,
             0.0059, 0.001, 0.0112, 0.001, 99.9466, 15.1436,0.0276, 237539.3249,
             61688.259, 0.96, 0.347, 12.2428736,1.2399, 1.5101, 0.347, 0.165, 6.9295,
             0.0178]
    θVirus = [1.0, 1.0] # k14f τ14 (Virus Parameters)
    append!(θVals,θVirus) #Append the virus parameters to the orginal parameters

    #Should IFN signaling be stochastic?
    if IFNStoch == :Stochastic
        #Keep most parameters the same
        θ = fill.(θVals,N,N)
        #kcat8 produces IFN, make it nonzero ~20% of the time (can be changed later)
        θ[11] .= θVals[11] .* rand(Bernoulli(0.2),N,N)
    else
        θ = θVals #Just keep the parameters as is (same for each cell)
    end

#Now we need to deal with the initial condition
  #Define the initial conditions
  u0 = zeros(N,N,species)

if Hetero == :Hetero
  noiseDistributions = AddNoise2States(0.5)
  for (i,index) in enumerate(nonZeroSpeciesIdx)
      u0[:,:,index] = rand(noiseDistributions[i],N,N)
  end

else
  #Loop through non-zero species and update their concentrations
  for (idx,val) in zip(nonZeroSpeciesIdx,nonZeroSpeciesValues)
      u0[:,:,idx] .=  val
  end
end

  #Finally need to set the DNA initial condition
  if infectionMethod == :ISD
    #Define a region on the domain where cells will be infected
    circleOrigin = [0,0] #Where is the center of the drop?
    circleRadiusSquared = N^2 #How big is the drop?
    #Calculate squared distances
    sqDist(x,c) = reduce(+, @. (x-c)^2)
    #Loop though cells and check if they are infected
    for currentCell in cellIndicies
        #Are the cells inside the infected region?
        if sqDist([currentCell[1],currentCell[2]],circleOrigin) ≤ circleRadiusSquared
            u0[currentCell,2] = m2c(1e3)
        end
    end

  elseif infectionMethod == :Virus
    #Assume a poisson ditribution to randomly choose each cell's level of infection
    probDistInfected = Poisson(moi)
    u0[:,:,2] = @. m2c(1e3*rand(probDistInfected,N,N))
  end

#Need to wrap everything up into the ParContainer struct
  mass = [u0[:,:,i] for i in nonZeroSpeciesIdx]
  deathParameter = ones(N,N)
  DNAReplicate = infectionMethod==:ISD ? 0 : 1
  #Keep track of infected cells (save time when infected, Inf means not infected)
  cellsInfected = fill(Inf,N,N) #Make constant when not testing
  cellsInfected[findall(u0[:,:,2] .> 0.0), 1] .= 0.0
  #Keep track of time of death (TOD)
  cellsDead = fill(Inf,N,N) #Inf implies alive
  #Create an array that keeps track of whether or not a cell has tried to infect neighbors
  infectFirstAttempt = trues(N,N)

#Create an instance of the structures
θ = ParContainer(
  θ, #Rate constants
  mass, #Mass balances
  deathParameter, # cell is dead? (1==alive, 0==dead)
  DNAReplicate, #can DNA replicate? (1 if virus)
  cellsInfected, #Time cell was infected
  cellsDead, #Time cell was killed
  infectFirstAttempt) #Has cell tried to infect neighbors?)

return ODEProblem(Model!,u0,tspan,θ)
end

###############################################################
# 5. Function to count the cell states (healthy,infected,dead)
###############################################################
function cellStates(t,θ)
    #Number of healthy cells at time t
    totaHealthy = sum(θ.cellsInfected .> t)
    #Number of dead cells at time t
    totalDead = sum(θ.cellsDead .< t)
    #Number of infected cells at time t
    totalInfected = nCells - totaHealthy - totalDead
    return [totaHealthy,totalInfected,totalDead]
end

###############################################################
# 6. Display Information about the Parameter Container
###############################################################

function Base.show(io::IO,p::ParContainer)
  #Is simulation det/stoch and homo/hetero
  attribute1 = isa(p.par[1],Array) ? "Stochastic" : "Deterministic" #True if stochastic
  attribute2 = all(p.mass[1][1].==p.mass[1]) ? "Homogeneous" : "Heterogeneous" #True is homo
  println("Cell Population: ",attribute1," and ",attribute2)

  #What is the Input for DNA?
  modelInput = p.DNAReplicate==1 ? "Virus" : "ISD"
  println("Model Input: ",modelInput)

  #How many cells initially infected?
  println("Initially Infected: ",100.0*sum(p.cellsInfected.==0.0)/nCells, "%")

  println("Currently Infected: ",100.0*sum(p.cellsInfected.<Inf)/nCells, "%")

  #How many cells are initially dead?
  println("Dead: ",100.0*((nCells - sum(p.cellsDead.==Inf))/nCells), "%")

  #How many cells are not responsive to infection?
  if attribute1 == "Stochastic"
    println("Cells with kcat7=0: ",100.0*sum(p.par[11].==0.0)/nCells, "%")
  end

end


text:
@testset "alternate selector" begin
    N = 3
    s = AlternateSelector(n=N)
    values = zeros(N)

    @test s(values; step = 1) == 1
    @test s(values; step = 2) == 2
    @test s(values; step = 3) == 3

    @test s(values; step = 4) == 1
    @test s(values; step = 5) == 2
    @test s(values; step = 6) == 3
end


text:
flatten(x::AbstractParticles) = Array(x)
function flatten(v::AbstractArray{<:AbstractParticles})
    m = reduce(hcat, flatten.(v))
    return Array(reshape(m, size(m, 1), size(v)...))
end

"""
    convert_to_inference_data(obj::AbstractParticles; kwargs...) -> InferenceData
    convert_to_inference_data(
        obj::AbstractVector{<:AbstractParticles};
        kwargs...,
    ) -> InferenceData
    convert_to_inference_data(
        obj::AbstractVector{<:AbstractArray{<:AbstractParticles}};
        kwargs...,
    ) -> InferenceData

Convert `MonteCarloMeasurements.AbstractParticles` to an [`InferenceData`](@ref).

`obj` may have the following types:

  - `::AbstractParticles`: Univariate draws from a single chain.
  - `::AbstractVector{<:AbstractParticles}`: Univariate draws from a vector of chains.
  - `::AbstractVector{<:AbstractArray{<:AbstractParticles}}`: Multivariate draws from a vector
    of chains.
"""
function convert_to_inference_data(obj::AbstractParticles; kwargs...)
    return convert_to_inference_data([obj]; kwargs...)
end
function convert_to_inference_data(obj::AbstractVector{<:AbstractParticles}; kwargs...)
    return convert_to_inference_data(flatten(flatten.(obj)); kwargs...)
end
function convert_to_inference_data(
    obj::AbstractVector{<:AbstractArray{<:AbstractParticles}}; kwargs...
)
    return convert_to_inference_data(flatten(flatten.(obj)); kwargs...)
end


text:
using BinaryBuilder

# These are the platforms built inside the wizard
platforms = [
    BinaryProvider.Linux(:i686, :glibc),
  BinaryProvider.Linux(:x86_64, :glibc),
  BinaryProvider.Linux(:aarch64, :glibc),
  BinaryProvider.Linux(:armv7l, :glibc),
  BinaryProvider.Linux(:powerpc64le, :glibc),
  BinaryProvider.MacOS(),
  BinaryProvider.Windows(:x86_64)
]

# If the user passed in a platform (or a few, comma-separated) on the
# command-line, use that instead of our default platforms
if length(ARGS) > 0
    platforms = platform_key.(split(ARGS[1], ","))
end
info("Building for $(join(triplet.(platforms), ", "))")

# Collection of sources required to build LZ4
sources = [
    "https://github.com/lz4/lz4.git" =>
    "dfed9fa1d77f0434306d377c4da1f7191d3ba08a",
]

script = raw"""
cd $WORKSPACE/srcdir
cd lz4/contrib/cmake_unofficial/
cmake -DCMAKE_INSTALL_PREFIX=/ -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain
make lz4_shared
touch lz4 lz4c lz4.exe lz4c.exe
make install/fast
exit

"""

products = prefix -> [
    LibraryProduct(prefix,"liblz4")
]

# Build the given platforms using the given sources
hashes = autobuild(pwd(), "LZ4", platforms, sources, script, products)


text:
module Simple

import Params
using Distributions

export SimpleTraffic, iterate

type SimpleTraffic
    traffic_probs :: Array{Float64, 1} # traffic probabilities of states
    A :: Array{Float64, 2} # transition probabilities
    state :: Int # current state
    traffic :: Bool # whether traffic exists
    occupancy :: Float64
    occupier :: Int
    SimpleTraffic(state = 1) = new(Params.traffic_probs, Params.traffic_trans_prob, state, false)
end

function iterate(t :: SimpleTraffic)
    t.state = rand(Categorical(t.A[:, t.state]))
    t.traffic = rand(Bernoulli(t.traffic_probs[t.state]))
    t.occupancy = t.traffic ? 0 : Inf
    t.occupier = t.traffic ? 0 : -1
    nothing
end

end