text:
using Test

m = page_model()
include("../src/components/CH4forcing.jl")

add_comp!(m, ch4forcing, :ch4forcing)

set_param!(m, :ch4forcing, :c_N2Oconcentration, readpagedata(m,"test/validationdata/c_n2oconcentration.csv"))
set_param!(m, :ch4forcing, :c_CH4concentration, readpagedata(m,"test/validationdata/c_ch4concentration.csv"))

##running Model
run(m)

@test !isnan(m[:ch4forcing, :f_CH4forcing][10])

forcing=m[:ch4forcing,:f_CH4forcing]
forcing_compare=readpagedata(m,"test/validationdata/f_ch4forcing.csv")

@test forcing ≈ forcing_compare rtol=1e-3


text:
# using LinearAlgebra
# include("sphere.jl")
##

"""
This function compute the associated Legendre polynomial P(l,m,x)
    for m >= 0, l >= m, and abs(x) <= 1
"""
function Plm(l::Int, m::Int, x::AbstractFloat)
    @assert 0 <= m <= l
    pmm = 1.0
    if m > 0
        somx2 = sqrt((1.0 - x) * (1.0 + x))
        fact = 1.0
        for i = 1:m
            pmm *= (-fact) * somx2
            fact += 2.0
        end
    end

    if l == m
        return pmm
    end
    pmmp1 = x * (2.0 * m + 1.0) * pmm
    if (l == (m + 1))
        return pmmp1
    end

    pll = 0.0

    for ll = (m+2):l
        pll = ((2.0 * ll - 1.0) * x * pmmp1 - (ll + m - 1.0) * pmm) / (ll - m)
        pmm = pmmp1
        pmmp1 = pll
    end

    return pll
end
"""
This function compute the normalized associated Legendre polynomial P(l,m,x)
    for m >= 0, l >= m, and abs(x) <= 1. It is alias  of GSL sf_legendre_Plm
"""
function sph_Plm(l::Int, m::Int, x::AbstractFloat)
    @assert 0 <= m <= l
    return sf_legendre_sphPlm(l, m, x)
end
"""
This function compute all of the normalized associated Legendre polynomial P(l,m)
    for m >= 0, l >= m, and abs(x) <= 1. It is alias  of GSL sf_legendre_array.
    This implementation include the Condon-Shortley phase factor.
"""
function Plm_array(l_max::Int, x::AbstractFloat; normal = :SPHARM)
    normalization = Dict(
        :SPHARM => GSL_SF_LEGENDRE_SPHARM,
        :NONE => GSL_SF_LEGENDRE_NONE,
        :SCHMIDT => GSL_SF_LEGENDRE_SCHMIDT,
        :FULL => GSL_SF_LEGENDRE_FULL,
    )
    return sf_legendre_array_e(normalization[normal], l_max, x,1)
end
"""
This function returns the index of Plm_array corresponding to P(l,m)
"""
function Plm_array_index(l::Int, m::Int)
    @assert 0 <= m <= l
    return sf_legendre_array_index(l,m) + 1
end
"""
This function returns the the general real spherical harmonics function Y(l,m) at
    pol angle and az_angle. l ∈ (0,..,n) and ∈ (-l,..l) m are the
    order and degree of spherical harmonics respectively.
    Pol_angle ∈ (0,..pi) and az_angle ∈ (0,2pi) are the polar angle and azimuthal angle
    coordinate.
"""
function real_spharm(l::Int, m::Int, pol_angle, az_angle)
    sqrt2 = sqrt(2.0)
    if m == 0
        return sph_Plm(l, m, cos(pol_angle))
    elseif m > 0
        return sqrt2 * cos(m * az_angle) * sph_Plm(l, m, cos(pol_angle))
    elseif m < 0
        return sqrt2 * sin(-m * az_angle) * sph_Plm(l, -m, cos(pol_angle))
    end
end
"""
This function returns all of general real spherical harmonics function Y(l,m)
    with order l.
    Pol_angle ∈ (0,..pi) and az_angle ∈ (0,2pi) are the polar angle and azimuthal angle
    coordinate.
"""
function real_spharm_array(l_max::Int, pol_angle, az_angle)
    sqrt2 = sqrt(2.0)
    len_array = (l_max + 1)^2
    plm_array = Plm_array(l_max, cos(pol_angle); normal = :SPHARM)
    spharm_array = Vector{Float64}(undef, len_array)
    for j = 1:len_array
        l, m = j2lm(j)
        idx = Plm_array_index(l, abs(m))
        if m == 0
            spharm_array[j] = plm_array[idx]
        elseif m > 0
            spharm_array[j] = sqrt2 * cos(m * az_angle) * plm_array[idx]
        elseif m < 0
            spharm_array[j] = sqrt2 * sin(-m * az_angle) * plm_array[idx]
        end
    end
    return spharm_array
end

"""
This function returns the order l and degree m of Y(l,m) from linear
    index j of spharm_array
"""
function j2lm(j)
    j = j - 1
    l = floor(Int, sqrt(j))
    m = j - l * (l + 1)
    return l, m
end
"""
This function returns the linear index j of spharm_array
    from the order l and degree m of Y(l,m)
"""
function lm2j(l, m)
    return l^2 + l + m + 1
end

"""
This function returns the The matrix of general real spherical
    harmonics function Y(l,m). The row index i  corresponds to points x_i and
    column index j obtained from j =  l^2 + l + m +1 where l and m are the order
    and degree of spherical harmonics respectively.
"""
function spharm_mat(l_max, pol_angle::T, az_angle::T) where {T<:AbstractVector}
    n_row = length(pol_angle)
    n_column = (l_max + 1)^2
    mat = Matrix{eltype(pol_angle)}(undef, n_row, n_column)
    for i = 1:n_row
        mat[i, :] .= real_spharm_array(l_max, pol_angle[i], az_angle[i])
    end
    return mat
end

"""
This function returns the The matrix of coefficient of SH function Y(l,m). The
    row index i  corresponds to points x_i and column index j obtained from j =  l^2 + l + m +1 where l and m are the order
    and degree of spherical harmonics respectively.
"""
function spharm_coefs(l_max::Int, x, pol_angle,az_angle)
    mat = spharm_mat(l_max, pol_angle, az_angle)
    coefs = mat \ x
    return coefs
end

ext: jl


text:
# Deps
using Plots
using Random
using StatsBase
using Statistics
using StatsPlots
using KernelDensity
using DelimitedFiles
default(titlefont=(14, "arial"),guidefont=(16, "arial"),tickfont=(12, "arial"))

#####################################################################################################
# Produce histogram of group length
#####################################################################################################

# Data dir
data_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Real-Data/Genome-Properties-hg38/"

# Read in data
grp_len_data = readdlm("$(data_dir)/len_group_histogram.txt")[:,1]

# We found that $84.1\%$ (18,770,638) of the CG-groups
grp_len_data[11]
grp_len_data[11] / sum(grp_len_data) * 100

# Whereas the remaining 15.9% (3,549,239) of the CG-groups contain from $12$ to 1318 bases
sum(grp_len_data[12:end])
sum(grp_len_data[12:end]) / sum(grp_len_data) * 100

# Largest group is 1002 nucleotides
findlast(x -> x == 1,grp_len_data .> 0)

# Normalize
grp_len_data = grp_len_data / sum(grp_len_data)
group_size = collect(1:length(grp_len_data))

# Cumululative results
# - 84.06% of them are between 1-11 bp.
# - 97.67% of them are between 1-25 bp.
# - 99.85% of them are between 1-50 bp.
# - 99.99% of them are between 1-100 bp.
map(i -> sum(grp_len_data[1:i] * 100),[11,25,50,100])

#####################################################################################################
# Produce histogram of group number of CGs
#####################################################################################################

# Data dir
data_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Real-Data/Genome-Properties-hg38/"

# Read in data
num_cpg_per_grp_data = readdlm("$(data_dir)/num_cpg_per_grp_histogram.txt")[:,1]

# Group with the most CG sites has 212 CpGs
findlast(x -> x == 1,num_cpg_per_grp_data .> 0)

#####################################################################################################
# Stats on analysis regions (hg38)
#####################################################################################################

# Data dir
data_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Real-Data/Genome-Properties-hg38/"

# Read in data
ar_data = readdlm("$(data_dir)/estimation_region_table_3kb.txt")

# Clean data
ar_data = hcat(ar_data[:,3] - ar_data[:,2] .+ 1, ar_data[:,3], ar_data[:,4])

# 1. Total number of R regions in the human genome that contain at least 1 CG-group: 919,321
sum(ar_data[:,3] .> 0.0)

# 2. How many of those contain 1 CG-group: 66
sum(ar_data[:,3] .== 1.0)

# 3. The max and min size in # bases of the regions containing 1 CG-group: (4000,4000)
ind = ar_data[:,3] .== 1.0;
minimum(ar_data[ind,1])
maximum(ar_data[ind,1])

# 4. The max and min size in # bases of the regions containing at least 2 CG-groups: (2905,3071)
ind = ar_data[:,3] .>= 2.0;
minimum(ar_data[ind,1])
maximum(ar_data[ind,1])

# 5. the max and min # of CpG sites within the regions that contain 1 CG-group: (1,2)
ind = ar_data[:,3] .== 1.0;
minimum(ar_data[ind,2])
maximum(ar_data[ind,2])

# 6. the max and min # of CpG sites within the regions that contain at least 2 CG-groups
ind = ar_data[:,3] .>= 2.0;
minimum(ar_data[ind,2])
maximum(ar_data[ind,2])

#####################################################################################################
# Produce plots
#####################################################################################################

# Data dir
data_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Real-Data/Genome-Properties-hg38/"

# Read in analysis region data
ar_data = readdlm("$(data_dir)/estimation_region_table_3kb.txt")

# Clean data
ar_data = hcat(ar_data[:,3] - ar_data[:,2] .+ 1, ar_data[:,3], ar_data[:,4])

# Read in group histogram data
grp_len_data = readdlm("$(data_dir)/len_group_histogram.txt")[:,1]
grp_len_data = grp_len_data / sum(grp_len_data)

# Read in data
num_cpg_per_grp_data = readdlm("$(data_dir)/num_cpg_per_grp_histogram.txt")[:,1]
num_cpg_per_grp_data = num_cpg_per_grp_data / sum(num_cpg_per_grp_data)

## Take random sample
Random.seed!(123);
ran_samp = sample(1:size(ar_data)[1], 10000;replace=false);

## Plot histogram of group size (bp)
p1 = plot(11:40, grp_len_data[11:40] * 100, seriestype=:bar, xlab="CpG Group size (bp)",
    ylab="Percentage(%)", label="", xlim=(10, 40), ylim=(0, 100));

## Plot histogram of group size (#CGs)
p2 = plot(1:100, num_cpg_per_grp_data[1:100] * 100, seriestype=:bar, xlab="Number of CpGs in groups",
    ylab="Percentage(%)", label="", xlim=(0, 20), ylim=(0, 100));

## Density plot CpG sites & CpG groups
ind_nonzero = ar_data[:,2] .> 0.0
p3 = plot(ar_data[ind_nonzero,2], seriestype=:density, label="# CpG sites",
    xlabel="# CpG sites/groups in analysis regions", ylabel="Density");
plot!(p3,ar_data[ind_nonzero,3],seriestype=:density,label="# CpG groups");

## Plot CpG sites vs CpG groups
p4 = plot(ar_data[ran_samp,2], ar_data[ran_samp,3], seriestype=:scatter,
    xlab="# CpG sites in analysis region", ylab="# CpG groups in analysis region", label="", alpha=0.25);
plot!(p4,1:400,1:400,label="");

# Collage
p = plot(p1, p2, p3, p4, layout=(2, 2), size=(1000, 750))

savefig(p,"$(data_dir)/Analysis-Region-3kb-hg38.pdf")

#####################################################################################################
# Representative estimation region
#####################################################################################################

# Data dir
data_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Real-Data/Genome-Properties-hg38/"

# Read in analysis region data
ar_data = readdlm("$(data_dir)/estimation_region_table_3kb.txt")

# Median number of CpG sites among analyzed regions
ar_data = ar_data[ar_data[:,4] .> 9,:]
med_num_cpgs = median(ar_data[:,4])
med_num_cg_grps = median(ar_data[:,5])

# Randomly choose region with median number of CpG sites
Random.seed!(123)
ar_data = ar_data[ar_data[:,4] .== med_num_cpgs,:]
chosen_region = ar_data[rand(1:size(ar_data)[1]),:]


text:
"""
charlm.jl: AutoGrad version (c) Emre Yolcu, 2016

This example implements an LSTM network for training character-level language
models. It takes as input a text file and trains the network to predict the
next character in a sequence. It can then be used to generate a sample that
resembles the original text.

To train the network, do `include("charlm.jl")` and run
`CharLM.train()`. Initial parameters can be provided as an optional argument to
`train`. Required form of the parameters can be understood by looking at the
definition of `weights` below. If not provided, default parameters are created
using `weights` with the appropriate arguments passed from `train`.

`train` accepts the following keyword arguments:

  - `datasrc`: Path to a file, or a URL. "The Complete Works of Shakespeare"
    from Project Gutenberg will be used by default.
  - `char_limit`: Maximum number of characters to read from `datasrc`
  - `epochs`: Number of epochs
  - `lr_init`: Initial learning rate
  - `lr_decay`: Learning rate decay
  - `decay_after`: Epoch after which to start decaying learning rate
  - `embedding_size`: Size of the embedding vector
  - `hidden_size`: Size of the LSTM internal state
  - `batch_size`: Number of sequences to train on in parallel
  - `sequence_length`: Number of steps to unroll the network for
  - `gclip`: Value to clip the gradient norm at
  - `pdrop`: Dropout probability
  - `seed`: Seed for the random number generator

At the end of each epoch, the cross entropy loss for the training data, the
learning rate, and the time taken by the epoch to complete are
printed. Optimized parameters and an array that maps integer indices to their
characters are returned at the end of training. Then, a sample text can be
generated by running `CharLM.generate(w, rv, n)` where `w` is the dictionary of
parameters, `rv` is the mapping of indices to characters, and `n` is the number
of characters to be generated.
"""

module CharLM

using AutoGrad
using Base.LinAlg: axpy!

sigm(x) = 1 ./ (1 + exp(-x))
@primitive sigm(x::Array),dy,y  (dy .* y .* (1 - y))

function xavier(fan_out, fan_in)
    scale = sqrt(6 / (fan_in + fan_out))
    return convert(Array{Float32}, 2 * scale * rand(fan_out, fan_in) - scale)
end

function weights(; input_size=0, output_size=0, embedding_size=0,
                 hidden_size=0, batch_size=0, init=xavier)
    w = Dict()
    for gate in (:ingate, :forget, :outgate, :change)
        w[symbol(:W_, gate)] = init(hidden_size, embedding_size + hidden_size)
        w[symbol(:b_, gate)] = (gate == :forget ? ones : zeros)(Float32, (hidden_size, 1))
    end
    w[:W_embedding] = init(embedding_size, input_size)
    w[:W_predict]   = init(output_size, hidden_size)
    w[:b_predict]   = zeros(Float32, (output_size, 1))
    return w
end

function lstm(w, input, hidden, cell)
    x       = vcat(input, hidden) # avoid four separate vcat operations
    ingate  = sigm(w[:W_ingate]  * x .+ w[:b_ingate]) # in fact we can probably combine these four operations into one
    forget  = sigm(w[:W_forget]  * x .+ w[:b_forget]) # then use indexing, or (better) subarrays to get individual gates
    outgate = sigm(w[:W_outgate] * x .+ w[:b_outgate])
    change  = tanh(w[:W_change]  * x .+ w[:b_change])
    cell    = cell .* forget + ingate .* change
    hidden  = outgate .* tanh(cell)
    return hidden, cell
end

function predict(w, hidden)
    output = w[:W_predict] * hidden .+ w[:b_predict]
    return output .- log(sum(exp(output), 1))
end

function dropout(x, pdrop)
    return x .* (rand(size(x)) .< (1 - pdrop)) / (1 - pdrop)
end

function forw(w, inputs; state=nothing, pdrop=0)
    batch_size = size(inputs[1], 2)
    hidden_size = size(w[:W_predict], 2)
    if state == nothing
        hidden = zeros(Float32, (hidden_size, batch_size))
        cell = zeros(Float32, (hidden_size, batch_size))
    else
        hidden, cell = state
    end
    outputs = Any[]
    for input in inputs
        hidden, cell = lstm(w, w[:W_embedding] * input, hidden, cell)
        pdrop > 0 && (hidden = dropout(hidden, pdrop))
        push!(outputs, predict(w, hidden))
    end
    return outputs, (hidden, cell)
end

function loss(w, inputs, targets; kwargs...)
    outputs = forw(w, inputs; kwargs...)[1]
    n = length(inputs)
    z = 0.0
    for t = 1:n
        z += sum(outputs[t] .* targets[t])
    end
    return -z / size(inputs[1], 2)
end

function gnorm(g)
    return mapreduce(vecnorm, +, 0, values(g))
end

function train(w=nothing; datasrc=nothing, char_limit=0, epochs=1, lr_init=1.0,
               lr_decay=0.95, decay_after=10, embedding_size=128,
               hidden_size=256, batch_size=50, sequence_length=50, gclip=5.0,
               pdrop=0, seed=0)
    seed > -1 && srand(seed)
    data, chars, char_to_index, index_to_char = loaddata(datasrc, batch_size, char_limit)
    vocab_size = length(char_to_index)
    if w == nothing
        w = weights(; input_size=vocab_size, output_size=vocab_size,
                    embedding_size=embedding_size, hidden_size=hidden_size,
                    batch_size=batch_size)
    end
    gradfun = grad(loss)

    for epoch = 1:epochs
        start_time = time()
        targets = Any[]
        inputs = Any[]
        loss_count = zeros(2)
        lr = lr_init * lr_decay^max(0, epoch - decay_after)
        T = length(data) - 1
        for t = 1:T
            push!(inputs, copy(data[t])) # why copy here? there is no overwriting in AutoGrad.
            push!(targets, copy(data[t + 1]))
            if (t % sequence_length == 0) || t == T
                loss_count[1] += loss(w, inputs, targets; pdrop=pdrop)
                loss_count[2] += length(inputs)
                g = gradfun(w, inputs, targets; pdrop=pdrop)
                gn = (gclip > 0 ? gnorm(g) : 0)
                gscale = (gn > gclip > 0 ? (gclip / gn) : 1)
                for p in keys(w)
                    axpy!(-lr * gscale, g[p], w[p])
                end
                empty!(inputs)
                empty!(targets)
            end
            if t % 1000 == 0
                elapsed_time = time() - start_time
                @printf(STDERR, "Epoch: %d, t: %d/%d, Loss: %.6f, LR: %.6f, Time: %.6f\n",
                        epoch, t, T, loss_count[1] / loss_count[2], lr, elapsed_time)
            end
        end
    end

    return w, index_to_char
end

function loaddata(datasrc, batch_size, char_limit=0)
    if datasrc == nothing
        datasrc = joinpath(AutoGrad.datapath, "pg100.txt")
    end
    if !isfile(datasrc)
        url = "http://www.gutenberg.org/cache/epub/100/pg100.txt"
        download(url,datasrc)
    end
    stream = open(datasrc)
    chars = Char[]
    char_to_index = Dict{Char, Int32}()
    while !eof(stream)
        c = read(stream, Char)
        get!(char_to_index, c, 1 + length(char_to_index))
        push!(chars, c)
        char_limit > 0 && length(chars) >= char_limit && break
    end
    info("Read: $(length(chars)) characters, $(length(char_to_index)) vocabulary")
    data = minibatch(chars, char_to_index, batch_size)
    index_to_char = Array(Char, length(char_to_index))
    for (c, i) in char_to_index
        index_to_char[i] = c
    end
    return data, chars, char_to_index, index_to_char
end

function minibatch(chars, char_to_index, batch_size)
    nbatch = div(length(chars), batch_size)
    data = Any[]
    for i = 1:nbatch
        d = zeros(Float32, (length(char_to_index), batch_size))
        for j = 1:batch_size
            d[char_to_index[chars[i + nbatch * (j - 1)]], j] = 1
        end
        push!(data, d)
    end
    return data
end

function generate(w, index_to_char, nchar)
    vocab_size = length(index_to_char)
    x = zeros(Float32, (vocab_size, 1))
    s = nothing
    j = 1
    for i = 1:nchar
        y, s = forw(w, Any[x]; state=s)
        x[j, 1] = 0
        j = sample(exp(y[1]))
        x[j, 1] = 1
        print(index_to_char[j])
    end
    println()
end

function sample(p)
    r = rand(Float32)
    for c = 1:length(p)
        r -= p[c]
        r < 0 && return c
    end
end

end  # module


text:
module MapsMocks

using CellwiseValues

export MockMap
export TestMap

import CellwiseValues: evaluate!
import CellwiseValues: return_size

struct MockMap{P} <: Map{P,1,P,1}
  val::P
end

function evaluate!(
  this::MockMap{P}, points::AbstractVector{P}, v::AbstractVector{P}) where P
  for (i,qi) in enumerate(points)
    v[i] = qi+this.val
  end
end

return_size(::MockMap, psize::Tuple{Int}) = psize

struct TestMap{P} <: Map{P,1,P,2}
  val::P
  dim::Int
end

function evaluate!(
  this::TestMap{P}, points::AbstractVector{P}, v::AbstractMatrix{P}) where P
  for j in 1:this.dim
    for (i,qi) in enumerate(points)
      v[j,i] = j*qi+this.val
    end
  end
end

return_size(this::TestMap, psize::Tuple{Int}) = (this.dim, psize[1])

end # module


text:
abstract type AbstractUncertainODEProblem end

@concrete struct MonteCarloProblem <: AbstractUncertainODEProblem
    prob <: EnsembleProblem
    kwargs

    function MonteCarloProblem(prob; u0=prob.u0, p=prob.p, kwargs...)
        u0 = recursive_convert(to_distribution, u0)
        p = recursive_convert(to_distribution, p)
        prob_func = function (prob, i, repeat)
            _u0 = _rand.(u0)
            _p = _rand.(p)
            return remake(prob, u0=_u0, p=_p)
        end
        ens_prob = EnsembleProblem(prob; prob_func=prob_func)
        return new{typeof(ens_prob), typeof(kwargs)}(ens_prob, kwargs)
    end
end

@concrete struct KoopmanProblem <: AbstractUncertainODEProblem
    prob <: ODEProblem
    u0_CoV <: Function
    p_CoV <: Function
    kwargs

    function KoopmanProblem(prob::ODEProblem; u0=prob.u0, p=prob.p, kwargs...)
        u0_CoV, p_CoV = make_conversion_functions(u, p)
        u0 = recursive_convert(to_distribution, getdata(u0))
        p = recursive_convert(to_distribution, getdata(p))
        prob = remake(prob; u0=u0, p=p)
        return new{typeof(prob), typeof(u0_CoV), typeof(p_CoV)}(prob, u0_CoV, p_CoV, kwargs)
    end
end
KoopmanProblem(prob::KoopmanProblem; u0=prob.u0, p=prob.p, kwargs...) = KoopmanProblem(remake(prob.prob; u0=u0, p=p); prob.kwargs..., kwargs...)
KoopmanProblem(f::Function, args...; kwargs...) = KoopmanProblem(ODEProblem(f, args...; kwargs...))

make_conversion_functions(u0, p) = ((u,p)->u, (u,p)->p)
function make_conversion_functions(u0::ComponentArray, p::ComponentArray)
    ax_u0 = getaxes(u0)
    ax_p = getaxes(p)
    u0_CoV = (u, p) -> ComponentArray(u, ax_u0)
    p_CoV = (u, p) -> ComponentArray(p, ax_p)
    return u0_CoV, p_CoV
end


text:
@testset "adjmat" begin
    @testset "undirected" begin
        adj = [0 1 0 1;
               1 0 1 0;
               0 1 0 1;
               1 0 1 0]
        deg = [2 0 0 0;
               0 2 0 0;
               0 0 2 0;
               0 0 0 2]
        lap = [2 -1 0 -1;
               -1 2 -1 0;
               0 -1 2 -1;
               -1 0 -1 2]
        norm_lap = [1. -.5 0. -.5;
                   -.5 1. -.5 0.;
                   0. -.5 1. -.5;
                   -.5 0. -.5 1.]
        scaled_lap = [0 -0.5 0 -0.5;
                      -0.5 0 -0.5 0;
                      0 -0.5 0 -0.5;
                      -0.5 0 -0.5 0]
        rw_lap = [1 -0.5 0 -0.5;
                  -0.5 1 -0.5 0;
                  0 -0.5 1 -0.5;
                  -0.5 0 -0.5 1]

        for T in [Int8, Int16, Int32, Int64, Int128, Float16, Float32, Float64]
            D = GraphLaplacians.degree_matrix(adj, T, dir=:out)
            @test D == T.(deg)
            @test D == GraphLaplacians.degree_matrix(adj, T, dir=:in)
            @test D == GraphLaplacians.degree_matrix(adj, T, dir=:both)
            @test eltype(D) == T

            L = Graphs.laplacian_matrix(adj, T)
            @test L == T.(lap)
            @test eltype(L) == T

            SL = GraphLaplacians.signless_laplacian(adj, T)
            @test SL == T.(adj + deg)
            @test eltype(SL) == T
        end

        for T in [Float16, Float32, Float64]
            NL = GraphLaplacians.normalized_laplacian(adj, T)
            @test NL ≈ T.(norm_lap)
            @test eltype(NL) == T

            SL = GraphLaplacians.scaled_laplacian(adj, T)
            @test SL ≈ T.(scaled_lap)
            @test eltype(SL) == T

            RW = GraphLaplacians.random_walk_laplacian(adj, T)
            @test RW == T.(rw_lap)
            @test eltype(RW) == T
        end
    end

    @testset "directed" begin
        adj = [0 2 0 3;
               0 0 4 0;
               2 0 0 1;
               0 0 0 0]
        degs = Dict(
            :out  => diagm(0=>[2, 2, 4, 4]),
            :in   => diagm(0=>[5, 4, 3, 0]),
            :both => diagm(0=>[7, 6, 7, 4]),
        )
        laps = Dict(
            :out  => degs[:out] - adj,
            :in   => degs[:in] - adj,
            :both => degs[:both] - adj,
        )
        norm_laps = Dict(
            :out  => I - diagm(0=>[1/2, 1/2, 1/4, 1/4])*adj,
            :in   => I - diagm(0=>[1/5, 1/4, 1/3, 0])*adj,
        )
        sig_laps = Dict(
            :out  => degs[:out] + adj,
            :in   => degs[:in] + adj,
            :both => degs[:both] + adj,
        )
        rw_laps = Dict(
            :out  => I - diagm(0=>[1/2, 1/2, 1/4, 1/4]) * adj,
            :in   => I - diagm(0=>[1/5, 1/4, 1/3, 0]) * adj,
            :both => I - diagm(0=>[1/7, 1/6, 1/7, 1/4]) * adj,
        )

        for T in [Int8, Int16, Int32, Int64, Int128, Float16, Float32, Float64]
            for dir in [:out, :in, :both]
                D = GraphLaplacians.degree_matrix(adj, T, dir=dir)
                @test D == T.(degs[dir])
                @test eltype(D) == T
            end
            @test_throws DomainError GraphLaplacians.degree_matrix(adj, dir=:other)

            for dir in [:out, :in, :both]
                L = Graphs.laplacian_matrix(adj, T, dir=dir)
                @test L == T.(laps[dir])
                @test eltype(L) == T
            end

            for dir in [:out, :in, :both]
                SL = GraphLaplacians.signless_laplacian(adj, T, dir=dir)
                @test SL == T.(sig_laps[dir])
                @test eltype(SL) == T
            end
        end

        for T in [Float32, Float64]
            for dir in [:out, :in]
                L = GraphLaplacians.normalized_laplacian(adj, T, dir=dir)
                @test L == T.(norm_laps[dir])
                @test eltype(L) == T
            end

            for dir in [:out, :in, :both]
                RW = GraphLaplacians.random_walk_laplacian(adj, T, dir=dir)
                @test RW == T.(rw_laps[dir])
                @test eltype(RW) == T
            end
        end
    end
end


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

"""
Sample environment, returning `EnvironmentValue` (default) or array.
"""
Base.rand(env::Environment) = EnvironmentValue(k => rand(dist) for (k, dist) in env)
Base.rand(fenv::Flat) = flatten(fenv.env, rand(fenv.env))

"""
Infer dimension of action space.
"""
act_dim(mdp::ASTMDP{<:State, SampleAction}) = sum(info.n for info in values(mdp.env_info))

"""
Infer dimension of state space.
"""
obs_dim(mdp::ASTMDP{ObservableState, <:Action}) = length(observe(mdp.sim))

"""
Return ordered list dictionary keys.
"""
orderedkeys(dict::Dict) = sort!(collect(keys(dict)))

"""
Flatten `EnvironmentValue` into array.
"""
function flatten(env::Environment, value::EnvironmentValue)
	return reduce(append!, flatten(env[k], value[k]) for k in orderedkeys(env))
end

"""
Reconstruct `EnvironmentValue` from array.
"""
function unflatten(mdp::ASTMDP{<:State, SampleAction}, action::AbstractVector{<:Real})
	value = EnvironmentValue()
	env = environment(mdp.sim)

	i = 0
	for k in orderedkeys(env)
		n = mdp.env_info[k].n
		array = action[i+1:i+n]
		value[k] = unflatten(env[k], array)
		i += n
	end
	return value
end

"""
Calculate log probability of sample from environment variable.
"""
function logprob(distribution::Any, value::Any, marginalize::Bool)
    logp = logpdf(distribution, value)
    if marginalize
        logp -= logpdf(distribution, mode(distribution))
    end
    return logp
end

"""
Calculate total log probability of environment value.
"""
function logprob(env::Environment, value::EnvironmentValue, marginalize::Bool)
    return sum(logprob(env[k], value[k], marginalize) for k in keys(env))
end

"""
Infer information about simulation environment.
"""
function infer_info(env::Environment)
    env_info = EnvironmentInfo()
	for (k, dist) in env
		sample = rand(dist)
		array = flatten(dist, sample)
		env_info[k] = VariableInfo(length(array), typeof(sample))
	end
    return env_info
end

"""
Infer type of state.
"""
function infer_state(sim::AbstractSimulation)
    try
        isempty(observe(sim)) ? UnobservableState : ObservableState
    catch e
        e isa UnimplementedError ? UnobservableState : throw(e)
    end
end

"""
    ASTMDP(sim::AbstractSimulation; kwargs...)

Constructor for ASTMDP object. Infers various properties of MDP. Allows kwargs to be
specified flatly and automatically assigns to correct level.
"""
function ASTMDP(sim::AbstractSimulation; kwargs...)
    reset!(sim)
    s_type = infer_state(sim)
    a_type = sim isa BlackBox ? SeedAction : SampleAction
    env_info = sim isa BlackBox ? EnvironmentInfo() : infer_info(environment(sim))
    mdp = ASTMDP{s_type, a_type}(; sim=sim, env_info=env_info)

    # automatically applies kwarg to reward if match is found
    for k in keys(kwargs)
        obj = k in fieldnames(ASTMDP) || !(k in fieldnames(Reward)) ? mdp : mdp.reward
        setproperty!(obj, k, kwargs[k])
    end

    mdp.reward.heuristic = mdp.episodic ? FinalHeuristic() : mdp.reward.heuristic
    RNG_TEMP[] = mdp.rng == Random.default_rng() ? Random.Xoshiro() : deepcopy(mdp.rng)
    return mdp
end

"""
Wrap action for internal handling.
"""
convert_a(mdp::ASTMDP{<:State, SampleAction}, action::AbstractVector{<:Real}) = SampleAction(unflatten(mdp, action))
convert_a(::ASTMDP{<:State, SampleAction}, action::EnvironmentValue) = SampleAction(action)
convert_a(::ASTMDP{<:State, SeedAction}, action::UInt32) = SeedAction(action)


text:
# Create example data
x = repeat(LinRange(-2, 2, 20), inner=10)
y = repeat(LinRange(-1, 1, 10), outer=20)
u = x .* (x.^2 .+ y.^2)
v = y .* (x.^2 .+ y.^2)
# Plot arrows
quiver(x, y, u, y, arrowscale=0.1)
# Expand the y-axes to see the whole arrows
ylim(-1.15, 1.15)


text:
#=
  NOTE: Codes below are adapted from
  https://github.com/brian-j-smith/Mamba.jl/blob/master/src/distributions/transformdistribution.jl
  The Mamba.jl package is licensed under the MIT License:
  > Copyright (c) 2014: Brian J Smith and other contributors:
  >
  > https://github.com/brian-j-smith/Mamba.jl/contributors
  >
  > Permission is hereby granted, free of charge, to any person obtaining
  > a copy of this software and associated documentation files (the
  > "Software"), to deal in the Software without restriction, including
  > without limitation the rights to use, copy, modify, merge, publish,
  > distribute, sublicense, and/or sell copies of the Software, and to
  > permit persons to whom the Software is furnished to do so, subject to
  > the following conditions:
  >
  > The above copyright notice and this permission notice shall be
  > included in all copies or substantial portions of the Software.
  >
  > THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  > EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  > MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
  > IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
  > CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
  > TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
  > SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#

#############
# a ≦ x ≦ b #
#############

const TransformDistribution{T<:ContinuousUnivariateDistribution} = Union{T, Truncated{T}}

function link(d::TransformDistribution, x::Real)
  a, b = minimum(d), maximum(d)
  lowerbounded, upperbounded = isfinite(a), isfinite(b)
  if lowerbounded && upperbounded
    logit((x - a) / (b - a))
  elseif lowerbounded
    log(x - a)
  elseif upperbounded
    log(b - x)
  else
    x
  end
end

function invlink(d::TransformDistribution, x::Real)
  a, b = minimum(d), maximum(d)
  lowerbounded, upperbounded = isfinite(a), isfinite(b)
  if lowerbounded && upperbounded
    (b - a) * invlogit(x) + a
  elseif lowerbounded
    exp(x) + a
  elseif upperbounded
    b - exp(x)
  else
    x
  end
end

function logpdf_with_trans(d::TransformDistribution, x::Real, transform::Bool)
  lp = logpdf(d, x)
  if transform
    a, b = minimum(d), maximum(d)
    lowerbounded, upperbounded = isfinite(a), isfinite(b)
    if lowerbounded && upperbounded
      lp += log((x - a) * (b - x) / (b - a))
    elseif lowerbounded
      lp += log(x - a)
    elseif upperbounded
      lp += log(b - x)
    end
  end
  lp
end

###############
# -∞ < x < -∞ #
###############

const RealDistribution = Union{Cauchy, Gumbel, Laplace, Logistic,
                               NoncentralT, Normal, NormalCanon, TDist}

link(d::RealDistribution, x::Real) = x

invlink(d::RealDistribution, x::Real) = x

logpdf_with_trans(d::RealDistribution, x::Real, transform::Bool) = logpdf(d, x)

#########
# 0 < x #
#########

const PositiveDistribution = Union{BetaPrime, Chi, Chisq, Erlang, Exponential, FDist, Frechet,
                                   Gamma, InverseGamma, InverseGaussian, Kolmogorov, LogNormal,
                                   NoncentralChisq, NoncentralF, Rayleigh, Weibull}

link(d::PositiveDistribution, x::Real) = log(x)

invlink(d::PositiveDistribution, x::Real) = exp(x)

function logpdf_with_trans(d::PositiveDistribution, x::Real, transform::Bool)
  return logpdf(d, x) + transform * log(x)
end

#############
# 0 < x < 1 #
#############

const UnitDistribution = Union{Beta, KSOneSided, NoncentralBeta}

link(d::UnitDistribution, x::Real) = logit(x)

invlink(d::UnitDistribution, x::Real) = invlogit(x)

function logpdf_with_trans(d::UnitDistribution, x::Real, transform::Bool)
  return logpdf(d, x) + transform * log(x * (one(x) - x))
end

###########
# ∑xᵢ = 1 #
###########

const SimplexDistribution = Union{Dirichlet}

function link(d::SimplexDistribution, x::AbstractVector{T}) where T<:Real
  y, K = similar(x), length(x)

  sum_tmp = zero(T)
  z = x[1]
  y[1] = logit(z) - log(one(T) / (K - 1))
  @inbounds for k in 2:K - 1
    sum_tmp += x[k - 1]
    z = x[k] / (one(T) - sum_tmp)
    y[k] = logit(z) - log(one(T) / (K - k))
  end
  y[K] = zero(T)

  return y
end

# Vectorised implementation of the above.
function link(d::SimplexDistribution, X::AbstractMatrix{T}) where T<:Real
  Y, K, N = similar(X), size(X, 1), size(X, 2)

  @inbounds for n in 1:size(X, 2)
    sum_tmp, z = zero(T), X[1, n]
    Y[1, n] = logit(z) - log(one(T) / (K - 1))
    for k in 2:K - 1
      sum_tmp += X[k - 1, n]
      z = X[k, n] / (one(T) - sum_tmp)
      Y[k, n] = logit(z) - log(one(T) / (K - k))
    end
    Y[K, n] = zero(T)
  end

  return Y
end

function invlink(d::SimplexDistribution, y::AbstractVector{T}) where T<:Real
  x, K = similar(y), length(y)

  z = invlogit(y[1] + log(one(T) / (K - 1)))
  x[1] = z
  sum_tmp = zero(T)
  @inbounds for k = 2:K - 1
    z = invlogit(y[k] + log(one(T) / (K - k)))
    sum_tmp += x[k-1]
    x[k] = (one(T) - sum_tmp) * z
  end
  sum_tmp += x[K-1]
  x[K] = one(T) - sum_tmp

  return x
end

# Vectorised implementation of the above.
function invlink(d::SimplexDistribution, Y::AbstractMatrix{T}) where T<:Real
  X, K, N = similar(Y), size(Y, 1), size(Y, 2)

  @inbounds for n in 1:size(X, 2)
    sum_tmp, z = zero(T), invlogit(Y[1, n] + log(one(T) / (K - 1)))
    X[1, n] = z
    for k in 2:K - 1
      z = invlogit(Y[k, n] + log(one(T) / (K - k)))
      sum_tmp += X[k - 1]
      X[k, n] = (one(T) - sum_tmp) * z
    end
    sum_tmp += X[K - 1, n]
    X[K, n] = one(T) - sum_tmp
  end

  return X
end

function logpdf_with_trans(
  d::SimplexDistribution,
  x::AbstractVector{<:Real},
  transform::Bool,
)
  lp = logpdf(d, x)
  if transform
    K = length(x)

    sum_tmp = zero(eltype(x))
    z = x[1]
    lp += log(z) + log(one(eltype(x)) - z)
    @inbounds for k in 2:K-1
      sum_tmp += x[k-1]
      z = x[k] / (one(eltype(x)) - sum_tmp)
      lp += log(z) + log(one(eltype(x)) - z) + log(one(eltype(x)) - sum_tmp)
    end
  end
  return lp
end

# REVIEW: why do we put this piece of code here?
function logpdf_with_trans(d::Categorical, x::Int)
  return d.p[x] > 0.0 && insupport(d, x) ? log(d.p[x]) : eltype(d.p)(-Inf)
end

#####################
# Positive definite #
#####################

const PDMatDistribution = Union{InverseWishart, Wishart}

function link(d::PDMatDistribution, X::AbstractMatrix{T}) where {T<:Real}
  Y = cholesky(X).L
  for m in 1:size(Y, 1)
    Y[m, m] = log(Y[m, m])
  end
  return Matrix(Y)
end

function invlink(d::PDMatDistribution, Y::AbstractMatrix{T}) where {T<:Real}
  X, dim = copy(Y), size(Y)
  for m in 1:size(X, 1)
    X[m, m] = exp(X[m, m])
  end
  return LowerTriangular(X) * LowerTriangular(X)'
end

function logpdf_with_trans(d::PDMatDistribution, X::AbstractMatrix{<:Real}, transform::Bool)
  lp = logpdf(d, X)
  if transform && isfinite(lp)
    U = cholesky(X).U
    for i in 1:dim(d)
      lp += (dim(d) - i + 2) * log(U[i, i])
    end
    lp += dim(d) * log(2.0)
  end
  return lp
end

############################################
# Defaults (assume identity link function) #
############################################

# UnivariateDistributions
using Distributions: UnivariateDistribution

link(d::UnivariateDistribution, x::Real) = x
link(d::UnivariateDistribution, x::AbstractVector{<:Real}) = link.(Ref(d), x)

invlink(d::UnivariateDistribution, x::Real) = x
invlink(d::UnivariateDistribution, x::AbstractVector{<:Real}) = invlink.(Ref(d), x)

logpdf_with_trans(d::UnivariateDistribution, x::Real, ::Bool) = logpdf(d, x)
function logpdf_with_trans(
  d::UnivariateDistribution,
  x::AbstractVector{<:Real},
  transform::Bool,
)
  return logpdf_with_trans.(Ref(d), x, transform)
end

# MultivariateDistributions
using Distributions: MultivariateDistribution

link(d::MultivariateDistribution, x::AbstractVector{<:Real}) = copy(x)
function link(d::MultivariateDistribution, X::AbstractMatrix{<:Real})
  Y = similar(X)
  for n in 1:size(X, 2)
    Y[:, n] = link(d, view(X, :, n))
  end
  return Y
end

invlink(d::MultivariateDistribution, y::AbstractVector{<:Real}) = copy(y)
function invlink(d::MultivariateDistribution, Y::AbstractMatrix{<:Real})
  X = similar(Y)
  for n in 1:size(Y, 2)
    X[:, n] = invlink(d, view(Y, :, n))
  end
  return X
end

function logpdf_with_trans(d::MultivariateDistribution, x::AbstractVector{<:Real}, ::Bool)
  return logpdf(d, x)
end
function logpdf_with_trans(
  d::MultivariateDistribution,
  X::AbstractMatrix{<:Real},
  transform::Bool,
)
  return [logpdf_with_trans(d, view(X, :, n), transform) for n in 1:size(X, 2)]
end

# MatrixDistributions
using Distributions: MatrixDistribution

link(d::MatrixDistribution, X::AbstractMatrix{<:Real}) = copy(X)
link(d::MatrixDistribution, X::AbstractVector{<:AbstractMatrix{<:Real}}) = link.(Ref(d), X)

invlink(d::MatrixDistribution, Y::AbstractMatrix{<:Real}) = copy(Y)
function invlink(d::MatrixDistribution, Y::AbstractVector{<:AbstractMatrix{<:Real}})
  return invlink.(Ref(d), Y)
end

logpdf_with_trans(d::MatrixDistribution, X::AbstractMatrix{<:Real}, ::Bool) = logpdf(d, X)
function logpdf_with_trans(
  d::MatrixDistribution,
  X::AbstractVector{<:AbstractMatrix{<:Real}},
  transform::Bool,
)
  return logpdf_with_trans.(Ref(d), X, Ref(transform))
end


text:
export Quarter, ExogTimeVars, _timevars

struct ExogTimeVars{N, ZTup<:NTuple{N,Real}, OR<:OrdinalRange}
    timevars::Vector{ZTup}
    timestamp::OR
    function ExogTimeVars(timevars::Vector{ZTup}, timestamp::OR) where {N, ZTup<:NTuple{N,Real}, OR<:OrdinalRange}
        length(timevars) .== length(timestamp) || throw(DimensionMismatch())
        new{N,ZTup,OR}(timevars,timestamp)
    end
end

_timevars(tv::ExogTimeVars) = tv.timevars
_timestamp(tv::ExogTimeVars) = tv.timestamp

getindex(tv::ExogTimeVars, t) = getindex(_timevars(tv), t)
function time_idx(tv::ExogTimeVars, t::Date)
    ts = _timestamp(tv)
    t in ts || throw(DomainError(t))
    return searchsortedfirst(ts,t)
end
getindex(tv::ExogTimeVars, t::Date) =
    getindex(_timevars(tv), time_idx(tv,t))

date(tv::ExogTimeVars, t::Integer)     = getindex(_timestamp(tv),t)
length(tv::ExogTimeVars) = length(_timestamp(tv))
size(tv::ExogTimeVars{N}) where {N} = length(tv), N
view(tv::ExogTimeVars, idx) = view(_timevars(tv), idx)

import Base: minimum, maximum

minimum(tv::ExogTimeVars{1}) = minimum(_timevars(tv))
maximum(tv::ExogTimeVars{1}) = maximum(_timevars(tv))

function DateQuarter(y::Integer, q::Integer)
    1 <= q <= 4 || throw(DomainError())
    return Date(y, 3*(q-1)+1)
end

function DateQuarter(yq::Real)
    y0, q0 = divrem(yq,1)
    y = Int(y0)
    q = Int(q0*4)+1
    return DateQuarter(y,q)
end

function Quarter(i::Integer)
    1 <= i <= 4 || throw(DomainError())
    return Month(3*i)
end


text:
using AzSessions, AzStorage, CloudSeis, Dates, Distributed, FolderStorage, HTTP, JSON, Random, Test

credentials = JSON.parse(ENV["AZURE_CREDENTIALS"])
AzSessions.write_manifest(;client_id=credentials["clientId"], client_secret=credentials["clientSecret"], tenant=credentials["tenantId"])

session = AzSession(;protocal=AzClientCredentials, client_id=credentials["clientId"], client_secret=credentials["clientSecret"], resource="https://storage.azure.com/")

const storageaccount1 = ENV["STORAGE_ACCOUNT1"]
const storageaccount2 = ENV["STORAGE_ACCOUNT2"]

abstract type Cloud end
struct Azure <: Cloud end
struct Azure2 <: Cloud end
struct POSIX <: Cloud end

mkcontainer(::Type{Azure}, foldername) = AzContainer(foldername; session=session, storageaccount=storageaccount1)
mkcontainer(::Type{Azure2}, foldername) = [AzContainer(foldername; session=session, storageaccount=a) for a in (storageaccount1, storageaccount2)]
mkcontainer(::Type{POSIX}, foldername) = Folder(foldername)

function csopen_robust(containers, mode; kwargs...)
    local io
    timeout = 180
    tic = time()
    while true
        try
           io = csopen(containers, mode; kwargs...)
           break
        catch e
            @warn "caught exception in csopen, sleeping for 60 seconds"
            showerror(stdout, e)
            sleep(60)
            if time() - tic > timeout
                throw(e)
            end
        end
    end
    io
end

const clouds = (Azure, Azure2, POSIX)
const compressors = ("none","blosc","leftjustify")

@testset "CloudSeis, cloud=$cloud, compresser=$compressor" for cloud in clouds, compressor in compressors
    @testset "CloudSeis, compression selection" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+0)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], compressor=compressor)
        if compressor == "none"
            @test isa(io.cache.compressor, CloudSeis.NotACompressor)
            @test Dict(io.cache.compressor)["method"] == "none"

            # default compression is "leftjustify"
            r = lowercase(randstring())
            io_default = csopen_robust(mkcontainer(cloud, "test-$r-default-cs"), "w", axis_lengths=[10,11,12])
            @test isa(io_default.cache.compressor, CloudSeis.LeftJustifyCompressor)
            rm(io_default)
        end
        if compressor == "blosc"
            @test isa(io.cache.compressor, CloudSeis.BloscCompressor)
            @test Dict(io.cache.compressor)["method"] == "blosc"
        end
        if compressor == "leftjustify"
            @test isa(io.cache.compressor, CloudSeis.LeftJustifyCompressor)
            @test Dict(io.cache.compressor)["method"] == "leftjustify"
        end
        rm(io)
    end

    @testset "allocframe" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+1)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], compressor=compressor)
        t,h = allocframe(io)

        @test size(t) == (10,11)
        @test eltype(t) == Float32

        @test size(h) == (io.hdrlength,11)
        @test eltype(h) == UInt8

        close(io)
        rm(io)
    end

    @testset "writeframe/readframe with headers and variable fold" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+2)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], compressor=compressor)

        x = rand(Float32,10,11,12)
        t,h = allocframe(io)

        n = size(h,2)
        J = Vector{Vector{Int}}(undef, size(io,3))
        for iframe = 1:size(io,3)
            t .= x[:,:,iframe]
            if iframe == 2
                J[iframe] = [1:11;]
            elseif iframe == 5
                J[iframe] = Int[]
            elseif iframe == 9
                J[iframe] = [6]
            else
                J[iframe] = randperm(n)[1:div(n,2)]
            end
            for itrace = 1:size(io,2)
                set!(prop(io,io.axis_propdefs[2]), h, itrace, itrace)
                set!(prop(io,io.axis_propdefs[3]), h, itrace, iframe)
                set!(prop(io,stockprop[:TRC_TYPE]), h, itrace, itrace ∈ J[iframe] ? tracetype[:live] : tracetype[:dead])
                if itrace ∉ J[iframe]
                    t[:,itrace] .= 0
                end
            end
            writeframe(io,t,h)
        end
        close(io)

        nbytes = filesize(io.extents[1].container, io.extents[1].name)
        if compressor == "leftjustify"
            @test nbytes == size(io,3)*8 + mapreduce(iframe->length(J[iframe])*(size(io,1)*sizeof(io.traceformat) + headerlength(io)), +, 1:size(io,3))
        elseif compressor == "none"
            @test nbytes == size(io,3)*8 + size(io,3)*size(io,2)*(4*size(io,1) + headerlength(io))
        elseif compressor == "blosc"
            @test nbytes < size(io,3)*8 + size(io,3)*size(io,2)*(4*size(io,1) + headerlength(io))
        end

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        for iframe = 1:size(io,3)
            t,h = readframe(io, iframe)
            for i = 1:size(h,2)
                if i ∈ J[iframe]
                    @test t[:,i] ≈ x[:,i,iframe]
                    @test get(prop(io,io.axis_propdefs[3]), h, i) == iframe
                end
                @test get(prop(io,io.axis_propdefs[2]), h, i) == i
                @test get(prop(io,stockprop[:TRC_TYPE]), h, i) == (i ∈ J[iframe] ? tracetype[:live] : tracetype[:dead])
            end
        end
        rm(io)
    end

    @testset "writeframe/readframe with index" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+3)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], dataproperties=[DataProperty("P",1)], compressor=compressor)

        x = rand(Float32,10,11)

        t,h = allocframe(io)
        t .= x

        writeframe(io, t, 1)
        close(io)

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        t,h = readframe(io, 1)
        close(io)
        @test t ≈ x
        rm(io)
    end

    @testset "write, lstarts=$lstarts, lincs=$lincs" for (lstarts,lincs) = (((1,1,1),(1,1,1)), ((2,3,4),(5,6,7)))
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+4)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_lstarts=lstarts, axis_lincs=lincs, compressor=compressor)

        x = rand(Float32,10,11,12)
        for ifrm = 1:12
            x[:,:,ifrm] .= ifrm
        end
        write(io, x, :, :, :)
        close(io)

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        for ifrm = 1:12
            _x = readframetrcs(io, lstarts[3] + lincs[3]*(ifrm - 1))
            @test _x ≈ x[:,:,ifrm]
        end
        close(io)
        rm(io)
    end

    @testset "write, multiple extents" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+5)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], frames_per_extent=1, compressor=compressor)

        x = rand(Float32,10,11,12)
        write(io, x, :, :, :)
        close(io)

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        for ifrm = 1:12
            _x = readframetrcs(io, ifrm)
            @test _x ≈ x[:,:,ifrm]
        end
        close(io)
        rm(io)
    end

    @testset "readtrcs, lstarts=$lstarts, lincs=$lincs" for (lstarts,lincs) = (((1,1,1),(1,1,1)), ((2,3,4),(5,6,7)))
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+6)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_lstarts=lstarts, axis_lincs=lincs, compressor=compressor)

        x = rand(Float32,10,11,12)
        write(io, x, :, :, :)
        close(io)

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        _x = readtrcs(io, :, :, :)
        close(io)
        rm(io)

        @test x ≈ _x
    end

    @testset "readhdrs, lstarts=$lstarts, lincs=$lincs" for (lstarts,lincs) = (((1,1,1),(1,1,1)), ((2,3,4),(5,6,7)))
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+7)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_lincs=lincs, axis_lstarts=lstarts, compressor=compressor)

        x = rand(Float32,10,11,12)
        write(io, x, :, :, :)
        close(io)

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        h = readhdrs(io, :, :, :)
        close(io)
        rm(io)

        for ifrm = 1:12, itrc = 1:11
            @test get(prop(io, "TRACE"), view(h, :, :, ifrm), itrc) ≈ lstarts[2] + lincs[2]*(itrc - 1)
            @test get(prop(io, "FRAME"), view(h, :, :, ifrm), itrc) ≈ lstarts[3] + lincs[3]*(ifrm - 1)
        end
    end

    @testset "readhdrs!, lstarts=$lstarts, lincs=$lincs" for (lstarts,lincs) = (((1,1,1),(1,1,1)), ((2,3,4),(5,6,7)))
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+8)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_lincs=lincs, axis_lstarts=lstarts, compressor=compressor)

        x = rand(Float32,10,11,12)
        write(io, x, :, :, :)
        close(io)

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        h = readhdrs(io, :, :, :)
        close(io)
        rm(io)

        for ifrm = 1:12, itrc = 1:11
            @test get(prop(io, "TRACE"), view(h, :, :, ifrm), itrc) ≈ lstarts[2] + lincs[2]*(itrc - 1)
            @test get(prop(io, "FRAME"), view(h, :, :, ifrm), itrc) ≈ lstarts[3] + lincs[3]*(ifrm - 1)
        end
    end

    @testset "partial read for foldmap, lstarts=$lstarts, lincs=$lincs" for (lstrts,lncs) = (((1,1,1),(1,1,1)), ((2,3,4),(5,6,7)))
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+9)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_lstarts=lstrts, axis_lincs=lncs, compressor=compressor)

        x = rand(Float32,10,11,12)
        write(io, x, :, :, :)
        close(io)

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        f = [fold(io,i) for i in lstarts(io,3) .+ lincs(io,3)*(0:size(io,3)-1)]
        close(io)
        rm(io)

        @test f == [11 for i=1:12]
    end

    @testset "similarto" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+10)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w",
            axis_pincs=[0.1,0.2,0.3],
            axis_lengths=[10,11,12],
            dataproperties=[DataProperty("P",1)],
            compressor=compressor)
        close(io)

        _io = csopen_robust(mkcontainer(cloud, "test-$r-sim-cs"), "w", similarto=mkcontainer(cloud, "test-$r-cs"))
        @test size(_io) == (10,11,12)
        @test pincs(_io)[1] ≈ 0.1
        @test pincs(_io)[2] ≈ 0.2
        @test pincs(_io)[3] ≈ 0.3

        rm(io)
        rm(_io)
    end

    @testset "similarto, changing extents" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+11)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w",
            axis_pincs=[0.1,0.2,0.3],
            axis_lengths=[10,11,12],
            dataproperties=[DataProperty("P",1)],
            compressor=compressor)
        @test length(io.extents) == 1
        close(io)

        _io = csopen_robust(mkcontainer(cloud, "test-$r-sim-cs"), "w", similarto=mkcontainer(cloud, "test-$r-cs"), frames_per_extent=1)
        @test size(_io) == (10,11,12)
        @test pincs(_io)[1] ≈ 0.1
        @test pincs(_io)[2] ≈ 0.2
        @test pincs(_io)[3] ≈ 0.3

        @test length(_io.extents) == 12
    end

    @testset "propdefs" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+12)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], compressor=compressor)
        @test propdefs(io,1).label == "SAMPLE"
        @test propdefs(io,2).label == "TRACE"
        @test propdefs(io,3).label == "FRAME"
        @test propdefs(io)[1].label == "SAMPLE"
        @test propdefs(io)[2].label == "TRACE"
        @test propdefs(io)[3].label == "FRAME"
        close(io)
        rm(io)
    end

    @testset "props" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+13)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], compressor=compressor)
        @test props(io,1).def.label == "SAMPLE"
        @test props(io,2).def.label == "TRACE"
        @test props(io,3).def.label == "FRAME"
        @test props(io)[1].def.label == "SAMPLE"
        @test props(io)[2].def.label == "TRACE"
        @test props(io)[3].def.label == "FRAME"
        close(io)
        rm(io)
    end

    @testset "vector props" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+14)))
        pdef = TracePropertyDef("X","XX",Vector{Float64},2)
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], tracepropertydefs=[pdef], compressor=compressor)
        t,h = allocframe(io)
        for i = 1:11
            set!(prop(io,"TRACE"), h, i, i)
            set!(prop(io,"FRAME"), h, i, 1)
            set!(prop(io,"TRC_TYPE"), h, i, tracetype[:live])
            set!(prop(io,"X"), h, i, i*[1.0,2.0])
        end
        writeframe(io,t,h)
        close(io)

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        t,h = readframe(io,1)

        for i = 1:11
            @test get(prop(io,"X"), h, i) ≈ i*[1.0,2.0]
        end
        close(io)
        rm(io)
    end

    @testset "pincs" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+15)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_pincs=[1.0,2.0,3.0], compressor=compressor)
        @test pincs(io,1) ≈ 1.0
        @test pincs(io,2) ≈ 2.0
        @test pincs(io,3) ≈ 3.0
        @test pincs(io)[1] ≈ 1.0
        @test pincs(io)[2] ≈ 2.0
        @test pincs(io)[3] ≈ 3.0
        close(io)
        rm(io)
    end

    @testset "pstarts" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+16)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_pstarts=[1.0,2.0,3.0], compressor=compressor)
        @test pstarts(io,1) ≈ 1.0
        @test pstarts(io,2) ≈ 2.0
        @test pstarts(io,3) ≈ 3.0
        @test pstarts(io)[1] ≈ 1.0
        @test pstarts(io)[2] ≈ 2.0
        @test pstarts(io)[3] ≈ 3.0
        close(io)
        rm(io)
    end

    @testset "lincs" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+18)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_lincs=[1,2,3], compressor=compressor)
        @test lincs(io,1) == 1
        @test lincs(io,2) == 2
        @test lincs(io,3) == 3
        @test lincs(io)[1] == 1
        @test lincs(io)[2] == 2
        @test lincs(io)[3] == 3
        close(io)
        rm(io)
    end

    @testset "lstarts" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+17)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_lstarts=[1,2,3], compressor=compressor)
        @test lstarts(io,1) == 1
        @test lstarts(io,2) == 2
        @test lstarts(io,3) == 3
        @test lstarts(io)[1] == 1
        @test lstarts(io)[2] == 2
        @test lstarts(io)[3] == 3
        close(io)
        rm(io)
    end

    @testset "units" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+19)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_units=["X","Y","Z"], compressor=compressor)
        @test units(io,1) == "X"
        @test units(io,2) == "Y"
        @test units(io,3) == "Z"
        @test units(io)[1] == "X"
        @test units(io)[2] == "Y"
        @test units(io)[3] == "Z"
        close(io)
        rm(io)
    end

    @testset "domains" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+20)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_domains=["X","Y","Z"], compressor=compressor)
        @test domains(io,1) == "X"
        @test domains(io,2) == "Y"
        @test domains(io,3) == "Z"
        @test domains(io)[1] == "X"
        @test domains(io)[2] == "Y"
        @test domains(io)[3] == "Z"
        close(io)
        rm(io)
    end

    @testset "geometry" begin
        sleep(1)
        g = Geometry(
            ox=1.0,oy=2.0,oz=3.0,
            ux=4.0,uy=5.0,uz=6.0,
            vx=7.0,vy=8.0,vz=9.0,
            wx=10.0,wy=11.0,wz=12.0,
            u1=1,un=2,
            v1=3,vn=4,
            w1=5,wn=6)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+21)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], geometry=g, compressor=compressor)

        _g = geometry(io)
        @test g.ox ≈ 1.0
        @test g.oy ≈ 2.0
        @test g.oz ≈ 3.0
        @test g.ux ≈ 4.0
        @test g.uy ≈ 5.0
        @test g.uz ≈ 6.0
        @test g.vx ≈ 7.0
        @test g.vy ≈ 8.0
        @test g.vz ≈ 9.0
        @test g.wx ≈ 10.0
        @test g.wy ≈ 11.0
        @test g.wz ≈ 12.0
        @test g.u1 == 1
        @test g.un == 2
        @test g.v1 == 3
        @test g.vn == 4
        @test g.w1 == 5
        @test g.wn == 6

        close(io)
        rm(io)
    end

    @testset "in" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+22)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], compressor=compressor)

        @test in(stockprop[:TRACE], io) == true
        @test in(stockprop[:CDP],io) == false

        close(io)
        rm(io)
    end

    @testset "dataproperty" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+23)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], dataproperties=[DataProperty("X",1)], compressor=compressor)

        @test dataproperty(io, "X") == 1

        close(io)
        rm(io)
    end

    @testset "hasdataproperty" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+24)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], dataproperties=[DataProperty("X",1)], compressor=compressor)

        @test hasdataproperty(io, "X") == true
        @test hasdataproperty(io, "Y") == false

        close(io)
        rm(io)
    end

    @testset "backward compatibility, dataproperty (1.2->1.1)" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+25)))
        c = mkcontainer(cloud, "test-$r-cs")
        _c = isa(c, Array) ? c[1] : c

        io = csopen_robust(c, "w", axis_lengths=[10,11,12], dataproperties=[DataProperty("X",1)], compressor=compressor)
        d = JSON.parse(read(_c, "description.json", String))
        @test isa(d["dataproperties"], Dict)
        d["dataproperties"] = [Dict("label"=>"X", "value"=>1)]
        write(_c, "description.json", json(d, 1))
        close(io)

        io = csopen_robust(c, "r")
        d = JSON.parse(read(_c, "description.json", String))
        @test isa(d["dataproperties"], Vector)
        @test dataproperty(io, "X") == 1

        close(io)
        rm(io)
    end

    @testset "copy!" begin
        sleep(1)
        r1 = lowercase(randstring(MersenneTwister(millisecond(now())+26),4))
        io1 = csopen_robust(mkcontainer(cloud, "test-$r1-cs"), "w", axis_lengths=[10,11,12], compressor=compressor)
        r2 = lowercase(randstring(MersenneTwister(millisecond(now())+27),4))
        io2 = csopen_robust(mkcontainer(cloud, "test-$r2-cs"), "w", axis_lengths=[10,11,12], compressor=compressor)

        h1 = allocframehdrs(io1)
        h2 = allocframehdrs(io2)

        for i = 1:11
            set!(prop(io1,stockprop[:TRACE]), h1, i, i)
            set!(prop(io1,stockprop[:FRAME]), h1, i, 2)
            set!(prop(io1,stockprop[:TRC_TYPE]), h1, i, tracetype[:live])
        end

        copy!(io2, h2, io1, h1)

        for i = 1:11
            @test get(prop(io2,stockprop[:TRACE]), h2, i) == i
            @test get(prop(io2,stockprop[:FRAME]), h2, i) == 2
            @test get(prop(io2,stockprop[:TRC_TYPE]), h2, i) == tracetype[:live]
        end

        close(io1)
        rm(io1)
        close(io2)
        rm(io2)
    end

    @testset "reduce" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+28),4))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,50], frames_per_extent=2, compressor=compressor)

        @test length(io.extents) == 25

        T = rand(Float32,10,11,50)
        H = rand(UInt8,headerlength(io),11,50)
        for i = 1:50
            writeframe(io, T[:,:,i], i)
            readframehdrs!(io, view(H,:,:,i), i)
        end
        close(io)

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        reduce(io)
        close(io)

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        for i = 1:50
            t,h = readframe(io, i)
            @test t ≈ T[:,:,i]
            @test h == H[:,:,i]
        end

        @test length(io.extents) == 1

        rm(io)
    end

    @testset "reduce, parallel" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+29),4))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,50], frames_per_extent=2, compressor=compressor)

        @test length(io.extents) == 25

        T = rand(Float32,10,11,50)
        H = rand(UInt8,headerlength(io),11,50)
        for i = 1:50
            writeframe(io, T[:,:,i], i)
            readframehdrs!(io, view(H,:,:,i), i)
        end
        close(io)

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        addprocs(5)
        @everywhere using AzStorage, CloudSeis, FolderStorage
        reduce(io; frames_per_extent=10)
        rmprocs(workers())
        close(io)

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        for i = 1:50
            t,h = readframe(io, i)
            @test t ≈ T[:,:,i]
            @test h == H[:,:,i]
        end

        @test length(io.extents) == 5

        rm(io)
    end

    @testset "robust cscreate" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+30),4))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], force=true, compressor=compressor)
        writeframe(io, rand(Float32,10,11), 1)
        close(io)
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], force=true, compressor=compressor)
        writeframe(io, rand(Float32,10,11), 1)
        close(io)
    end

    @testset "backward compatibility (1.1->1.0)" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+31),4))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], force=true, compressor="none")
        description = JSON.parse(read(io.containers[1], "description.json", String))
        delete!(description, "compressor")
        write(io.containers[1], "description.json", json(description))
        io = csopen(mkcontainer(cloud, "test-$r-cs"), axis_lengths=[10,11,12], force=true, compressor="none")
        @test isa(io.cache, CloudSeis.Cache{CloudSeis.NotACompressor})
    end

    @testset "similar with new axis lengths" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+32),4))
        container = mkcontainer(cloud, "test-$r-cs")
        io = csopen_robust(container, "w", axis_lengths=[10,11,2])

        trcs1 = rand(Float32,10,11)
        trcs2 = rand(Float32,10,11)

        writeframe(io, trcs1, 1)
        writeframe(io, trcs2, 2)
        close(io)

        container_similar = mkcontainer(cloud, "test-$r-similar-cs")
        io_similar = csopen_robust(container_similar, "w"; similarto=container, axis_lengths=[10,11,1])
        writeframe(io_similar, trcs1, 1)
        close(io_similar)

        io = csopen_robust(container, "r")
        io_similar = csopen_robust(container_similar, "r")
        _trcs1 = readframetrcs(io_similar, 1)
        @test _trcs1 ≈ trcs1

        rm(io)
        rm(io_similar)
    end

    @testset "csopen for 6 dimensional data-set" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+33),4))
        container = mkcontainer(cloud, "test-$r-cs")
        io = csopen_robust(container, "w", axis_lengths=[10,11,2,2,3,1])
        close(io)
        io = csopen_robust(container, "r")
        @test size(io) == (10,11,2,2,3,1)
        labels = [propdefs(io)[i].label for i=1:6]
        @test labels == ["SAMPLE", "TRACE", "FRAME", "VOLUME", "HYPRCUBE", "DIM6"]
    end

    @testset "non unitary logical start and increments" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+34),4))
        container = mkcontainer(cloud, "test-$r-cs")

        io = csopen_robust(container, "w", axis_lengths=[12,11,4], axis_lstarts=[1,2,5], axis_lincs=[1,1,3])

        # test map from cartesian to linear index
        @test CloudSeis.linearframeidx(io, 5) == 1
        @test CloudSeis.linearframeidx(io, 8) == 2
        @test CloudSeis.linearframeidx(io, 11) == 3
        @test CloudSeis.linearframeidx(io, 14) == 4

        # test map from unitary to non-unitary cartesian index
        @test CloudSeis.logicalframeidx(io, 1) == CartesianIndex((5,))
        @test CloudSeis.logicalframeidx(io, 2) == CartesianIndex((8,))
        @test CloudSeis.logicalframeidx(io, 3) == CartesianIndex((11,))
        @test CloudSeis.logicalframeidx(io, 4) == CartesianIndex((14,))

        @test CloudSeis.lineartraceidx(io, 2) == 1
        @test CloudSeis.logicaltraceidx(io, 1) == 2

        rm(io)
    end

    @testset "non unitary logical start and increments, read/write frame" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+35)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_lstarts=[0,2,6], axis_lincs=[1,2,3], compressor=compressor)

        x = rand(Float32,10,11,12)
        t,h = allocframe(io)

        n = size(h,2)
        J = Vector{Vector{Int}}(undef, size(io,3))
        for iframe = 1:size(io,3)
            t .= x[:,:,iframe]
            if iframe == 2
                J[iframe] = [1:11;]
            elseif iframe == 5
                J[iframe] = Int[]
            elseif iframe == 9
                J[iframe] = [6]
            else
                J[iframe] = randperm(n)[1:div(n,2)]
            end
            for itrace = 1:size(io,2)
                set!(prop(io,io.axis_propdefs[2]), h, itrace, 2 + 2*(itrace - 1))
                set!(prop(io,io.axis_propdefs[3]), h, itrace, 6 + 3*(iframe - 1))
                set!(prop(io,stockprop[:TRC_TYPE]), h, itrace, itrace ∈ J[iframe] ? tracetype[:live] : tracetype[:dead])
                if itrace ∉ J[iframe]
                    t[:,itrace] .= 0
                end
            end
            writeframe(io,t,h)
        end
        close(io)

        nbytes = filesize(io.extents[1].container, io.extents[1].name)
        if compressor == "leftjustify"
            @test nbytes == size(io,3)*8 + mapreduce(iframe->length(J[iframe])*(size(io,1)*sizeof(io.traceformat) + headerlength(io)), +, 1:size(io,3))
        elseif compressor == "none"
            @test nbytes == size(io,3)*8 + size(io,3)*size(io,2)*(4*size(io,1) + headerlength(io))
        elseif compressor == "blosc"
            @test nbytes < size(io,3)*8 + size(io,3)*size(io,2)*(4*size(io,1) + headerlength(io))
        end

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        for iframe = 1:size(io,3)
            t,h = readframe(io, 6 + 3*(iframe-1))
            for i = 1:size(h,2)
                if i ∈ J[iframe]
                    @test t[:,i] ≈ x[:,i,iframe]
                    @test get(prop(io,io.axis_propdefs[3]), h, i) == 6 + 3*(iframe-1)
                end
                @test get(prop(io,io.axis_propdefs[2]), h, i) == 2 + 2*(i - 1)
                @test get(prop(io,stockprop[:TRC_TYPE]), h, i) == (i ∈ J[iframe] ? tracetype[:live] : tracetype[:dead])
            end
        end
        rm(io)
    end

    @testset "non unitary logical start and increments, writeframe/readframe with index" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+36)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_lstarts=[0,2,6], axis_lincs=[1,2,3], dataproperties=[DataProperty("P",1)], compressor=compressor)

        x = rand(Float32,10,11)

        t,h = allocframe(io)
        t .= x

        writeframe(io, t, 9)
        close(io)

        io = csopen(mkcontainer(cloud, "test-$r-cs"))
        t,h = readframe(io, 9)
        close(io)
        @test t ≈ x
        rm(io)
    end

    @testset "non unitary logical start and increments, exception for in-between out-of-bounds index" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+37)))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_lstarts=[1,2,3], axis_lincs=[4,5,6], dataproperties=[DataProperty("P",1)], compressor=compressor)

        @test_throws ErrorException CloudSeis.linearframeidx(io, 4)
        @test @inbounds CloudSeis.linearframeidx(io, 2) == 1

        @test_throws ErrorException CloudSeis.lineartraceidx(io, 3)
        @test @inbounds CloudSeis.linearframeidx(io, 3) == 1
    end

    @testset "backwards compatability, logical increments and starts" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+38),4))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], force=true)
        description = JSON.parse(read(io.containers[1], "description.json", String))
        delete!(description, "axis_lincs")
        delete!(description, "axis_lstarts")
        write(io.containers[1], "description.json", json(description))
        io = csopen(mkcontainer(cloud, "test-$r-cs"), axis_lengths=[10,11,12], force=true)
        @test lstarts(io) == (1,1,1)
        @test lincs(io) == (1,1,1)
    end

    @testset "frame iterator, LogicalIndices" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+39),4))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], force=true)
        idxs = CartesianIndex[]
        for idx in LogicalIndices(io)
            push!(idxs, idx)
        end
        _idxs = [CartesianIndex((i,)) for i=1:12]
        @test idxs == _idxs
        rm(io)

        r = lowercase(randstring(MersenneTwister(millisecond(now())+40),4))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], axis_lstarts=[1,1,3], axis_lincs=[1,1,4], force=true)
        idxs = CartesianIndex[]
        for idx in LogicalIndices(io)
            push!(idxs, idx)
        end
        _idxs = [CartesianIndex((3 + (i-1)*4,)) for i=1:12]
        @test idxs == _idxs
        rm(io)

        r = lowercase(randstring(MersenneTwister(millisecond(now())+41),4))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,4,6], force=true)
        idxs = CartesianIndex[]
        for idx in LogicalIndices(io)
            push!(idxs, idx)
        end
        _idxs = [CartesianIndex((i,j)) for i=1:4, j=1:6][:]
        @test idxs == _idxs
        rm(io)

        r = lowercase(randstring(MersenneTwister(millisecond(now())+42),4))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,4,6], axis_lstarts=[1,1,2,3], axis_lincs=[1,1,4,5], force=true)
        idxs = CartesianIndex[]
        for idx in LogicalIndices(io)
            push!(idxs, idx)
        end
        _idxs = [CartesianIndex((2+(i-1)*4,3+(j-1)*5)) for i=1:4, j=1:6][:]
        @test idxs == _idxs
        rm(io)
    end

    @testset "partial read for foldmap with empty extents" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+43),4))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,11,12], force=true)
        @test fold(io,1) == 0
    end

    @testset "left justify headers" begin
        sleep(1)
        r = lowercase(randstring(MersenneTwister(millisecond(now())+44),4))
        io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,12,1], force=true)

        h = allocframehdrs(io)
        for itrace = 1:12
            set!(prop(io, stockprop[:TRACE]), h, itrace, itrace)
            set!(prop(io, stockprop[:FRAME]), h, itrace, 1)
            if rem(itrace,2) == 0
                set!(prop(io, stockprop[:TRC_TYPE]), h, itrace, tracetype[:live])
            else
                set!(prop(io, stockprop[:TRC_TYPE]), h, itrace, tracetype[:dead])
            end
        end

        CloudSeis.leftjustify!(io, h)
        for itrace = 1:6
            @test get(prop(io, stockprop[:TRC_TYPE]), h, itrace) == tracetype[:live]
            @test get(prop(io, stockprop[:TRACE]), h, itrace) == 2*itrace
        end
        for itrace = 7:12
            @test get(prop(io, stockprop[:TRC_TYPE]), h, itrace) == tracetype[:dead]
        end
    end
end

@testset "CloudSeis optimzed read, cloud=$cloud" for cloud in clouds
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+45),4))
    io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,12,3], force=true)

    t,h = allocframe(io)

    for iframe = 1:3
        for itrace = 1:12
            set!(prop(io, stockprop[:TRACE]), h, itrace, itrace)
            set!(prop(io, stockprop[:FRAME]), h, itrace, iframe)
            if rem(itrace,2) == 0
                set!(prop(io, stockprop[:TRC_TYPE]), h, itrace, tracetype[:live])
                t[:,itrace] .= itrace*iframe
            else
                set!(prop(io, stockprop[:TRC_TYPE]), h, itrace, tracetype[:dead])
                t[:,itrace] .= 0
            end
        end

        writeframe(io, t, h)
    end

    close(io)
    io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "r")

    for iframe = 1:3
        _t,_h = readframe(io, iframe; regularize=false)

        for itrace = 1:6
            @test _t[:,itrace] ≈ 2*itrace*iframe*ones(Float32,10)
            @test get(prop(io, stockprop[:TRACE]), _h, itrace) == 2*itrace
            @test get(prop(io, stockprop[:FRAME]), _h, itrace) == iframe
            @test get(prop(io, stockprop[:TRC_TYPE]), _h, itrace) == tracetype[:live]
        end

        for itrace = 7:12
            @test get(prop(io, stockprop[:TRC_TYPE]), _h, itrace) == tracetype[:dead]
        end
    end

    rm(io)
end

@testset "CloudSeis optimized read empty frame, cloud=$cloud" for cloud in clouds
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+46),4))
    io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,12,2], force=true)

    t,h = allocframe(io)

    for itrace = 1:12
        set!(prop(io, stockprop[:TRACE]), h, itrace, itrace)
        set!(prop(io, stockprop[:FRAME]), h, itrace, 1)
        if rem(itrace,2) == 0
            set!(prop(io, stockprop[:TRC_TYPE]), h, itrace, tracetype[:live])
            t[:,itrace] .= itrace
        else
            set!(prop(io, stockprop[:TRC_TYPE]), h, itrace, tracetype[:dead])
            t[:,itrace] .= 0
        end
    end

    writeframe(io, t, h)

    close(io)
    io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "r")

    _t,_h = readframe(io, 2; regularize=false)

    for itrace = 1:12
        @test get(prop(io, stockprop[:TRC_TYPE]), _h, itrace) == tracetype[:dead]
    end

    rm(io)
end

@testset "CloudSeis optimized read, mixing regularize on/off, cloud=$cloud" for cloud in clouds
    sleep(1)
    r = lowercase(randstring(MersenneTwister(millisecond(now())+46),4))
    io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "w", axis_lengths=[10,12,2], force=true)

    t,h = allocframe(io)

    for iframe = 1:2
        for itrace = 1:12
            set!(prop(io, stockprop[:TRACE]), h, itrace, itrace)
            set!(prop(io, stockprop[:FRAME]), h, itrace, iframe)
            if rem(itrace,2) == 0
                set!(prop(io, stockprop[:TRC_TYPE]), h, itrace, tracetype[:live])
                t[:,itrace] .= itrace*iframe
            else
                set!(prop(io, stockprop[:TRC_TYPE]), h, itrace, tracetype[:dead])
                t[:,itrace] .= 0
            end
        end

        writeframe(io, t, h)
    end

    close(io)
    io = csopen_robust(mkcontainer(cloud, "test-$r-cs"), "r")

    _t,_h = readframe(io, 1; regularize=false)
    _t,_h = readframe(io, 2)

    for itrace = 1:12
        if rem(itrace,2) == 0
            @test _t[:,itrace] ≈ itrace*2*ones(Float32,10)
            @test get(prop(io, stockprop[:TRACE]), _h, itrace) == itrace
            @test get(prop(io, stockprop[:FRAME]), _h, itrace) == 2
            @test get(prop(io, stockprop[:TRC_TYPE]), _h, itrace) == tracetype[:live]
        else
            @test get(prop(io, stockprop[:TRC_TYPE]), _h, itrace) == tracetype[:dead]
        end
    end

    rm(io)
end


text:
using StanSample, Random

Random.seed!(1)

cd(@__DIR__)
include("simulatePoisson.jl")

HP = "
data {
    int N;
    int y[N];
    int Ns;
    int idx[N];
    real x[N];
  }
  parameters {
    real a0;
    vector[Ns] a0s;
    real a1;
    real<lower=0> a0_sig;
  }
  model {
    vector[N] mu;
    a0 ~ normal(0, 10);
    a1 ~ normal(0, 1);
    a0_sig ~ cauchy(0, 1);
    a0s ~ normal(0, a0_sig);
    for(i in 1:N) mu[i] = exp(a0 + a0s[idx[i]] + a1 * x[i]);
    y ~ poisson(mu);
  }
"

y, x, idx, N, Ns = simulatePoisson(;Nd=1,Ns=10,a0=1.0,a1=.5,a0_sig=.3)

hp_data = Dict(:N =>N, :y => y, :Ns => Ns, :idx => idx, :x => x)

stanmodel = SampleModel("Hierachical_Poisson", HP;
  method = StanSample.Sample(
    num_samples=2000,
    adapt = StanSample.Adapt(delta = 0.8)));

rc = stan_sample(stanmodel; data=hp_data, n_chains=4)

if success(rc)
  sdf = read_summary(stanmodel)
  sdf |> display
end


text:
## Example script for linear-strength vortex panel method
using LinearAlgebra
using Base.Iterators
using AeroMDAO

## Airfoil
alpha_u = [0.2, 0.3, 0.2, 0.15, 0.2]
alpha_l = [-0.2, -0.1, -0.1, -0.001]
dzs     = (0., 0.)
# airfoil = kulfan_CST(alpha_u, alpha_l, dzs, 0.0, 60);
airfoil = naca4((0,0,1,2))
V, α    = 1., 5.
ρ       = 1.225
uniform = Uniform2D(V, α)
num_pans = 80
panels  = make_panels(airfoil, num_pans);

## Vortex panel
A = @time vortex_influence_matrix(panels)
b = @time neumann_boundary_condition(panels, velocity(uniform))
γs = A \ b

##
cl = 2 * sum(@. (γs[1:end-1] + γs[2:end]) / 2 * panel_length(panels)) / uniform.magnitude

## Velocities and pressure coefficient
# panel_vels = [ velocity(uniform) .+ sum(vortex_velocity.(γs[1:end-1], γs[2:end], panels, Ref(panel_i))) for panel_i in panels ]

## Plotting
using CairoMakie

CairoMakie.activate!()

using LaTeXStrings

const LS = LaTeXString

## Pressure coefficient
qts = γs
cps = @. 1 - qts^2 / uniform.magnitude^2

upper, lower = get_surface_values(panels, cps)
upper, lower = [ lower[end]; upper ], [ upper[end]; lower ]

# Plot
f1 = Figure(resolution = (900, 400))
ax = Axis(f1[1,1], aspect = 1, yreversed = true, xlabel = L"(x/c)", ylabel = L"C_p")
l1 = lines!(upper, label = L"\mathrm{Upper}")
l2 = lines!(lower, label = L"\mathrm{Lower}")

axislegend(L"\mathrm{Surface}")

# Velocity field
x_domain, y_domain = (-1, 2), (-1, 1)
grid_size = 50
x_dom, y_dom = range(x_domain..., length = grid_size), range(y_domain..., length = grid_size)
grid = product(x_dom, y_dom)
pts = panel_points(panels);

total_velocity(xs, ys) = map((x, y) -> Point2f(velocity(uniform) .+ sum(vortex_velocity.(γs[1:end-1], γs[2:end], panels, x, y))...), xs, ys)

vortex_vels = total_velocity(first.(grid), last.(grid))
speeds = @. sqrt(first(vortex_vels)^2 + last(vortex_vels)^2)
cps = @. 1 - speeds^2 / uniform.magnitude^2

# Plot
ax2 = Axis(f1[1,2], aspect = 1)
hm = contourf!(x_dom, y_dom, speeds, levels = 50)

streamplot!(total_velocity, first.(grid), last.(grid), colormap = Reverse(:viridis))
poly!(pts, color = "black")

Colorbar(f1[1,3], hm, label = L"Normalized Speed $(U/U_\infty)$")
colsize!(f1.layout, 2, Aspect(1, 1.0))

# f1[0,:] = Label(f1, LS("Linear Vortex Panel Method"))

# arrows!(x_dom, y_dom, first.(vortex_vels), last.(vortex_vels))
# CB2 = colorbar(orientation = "horizontal", label = "Relative Airspeed (U/U∞)")
# tight_layout()
# show()
f1

##
save("plots/LinearVortex.svg", f1, px_per_unit = 1.5)


text:
struct BudgetedBipartiteMatchingInstance{T, U} <: CombinatorialInstance
  matching::BipartiteMatchingInstance{T}
  weight::Dict{Edge{T}, U}
  budget::U

  function BudgetedBipartiteMatchingInstance(graph::AbstractGraph{T}, reward::Dict{Edge{T}, Float64}, weight::Dict{Edge{T}, U}, budget::U) where {T, U}
    matching = BipartiteMatchingInstance(graph, reward)
    return new{T, U}(matching, weight, budget)
  end
end

function _budgeted_bipartite_matching_compute_value(i::BudgetedBipartiteMatchingInstance{T, U}, solution::Vector{Edge{T}}) where {T, U}
  return _budgeted_bipartite_matching_compute_value(i.matching, solution)
end

function _budgeted_bipartite_matching_compute_value(i::BipartiteMatchingInstance{T}, solution::Vector{Edge{T}}) where T
  return sum(i.reward[(e in keys(i.reward)) ? e : reverse(e)] for e in solution)
end

function _budgeted_bipartite_matching_compute_weight(i::BudgetedBipartiteMatchingInstance{T, U}, solution::Vector{Edge{T}}) where {T, U}
  return sum(i.weight[(e in keys(i.weight)) ? e : reverse(e)] for e in solution)
end

abstract type BudgetedBipartiteMatchingSolution{T, U} <: CombinatorialSolution
  # instance::BudgetedBipartiteMatchingInstance{T, U}
  # solution::Vector{Edge{T}}
  # value::Float64 # TODO: remove me, only useful for Lagrangian.
end


text:
function warning_callback(message::String)
    global TEST_CALLBACK = true
    error("checking that error are supported")
end

@testset "Error Functions" begin
    err = ChemfilesError("oops")
    iobuf = IOBuffer(19 + length(err.message))
    show(iobuf, err)
    @test String(iobuf.data) == "\"Chemfiles error: oops\""

    @test Chemfiles.last_error() == ""

    @test_throws ChemfilesError Residue(Topology(), 3)
    @test_throws UndefVarError TEST_CALLBACK == false

    @test Chemfiles.last_error() == "Residue index out of bounds in topology: we have 0 residues, but the index is 3"

    Chemfiles.clear_errors()
    @test Chemfiles.last_error() == ""

    Chemfiles.set_warning_callback(warning_callback)
    @test_throws ChemfilesError Residue(Topology(), 3)
    @test TEST_CALLBACK == true

    Chemfiles.set_warning_callback(Chemfiles._warning_callback)
end

@testset "Configuration files" begin
    config = joinpath(@__DIR__, "data", "config.toml")
    Chemfiles.add_configuration(config)

    trajectory = joinpath(@__DIR__, "data", "water.xyz")
    frame = read(Trajectory(trajectory))

    @test name(Atom(frame, 9)) == "Oz"
    @test atom_type(Atom(frame, 9)) == "F"
end


text:
# Kafka wire protocol. See details at:
# https://cwiki.apache.org/confluence/display/KAFKA/A+Guide+To+The+Kafka+Protocol#AGuideToTheKafkaProtocol-Messagesets

immutable RequestHeader
    api_key::Int16
    api_version::Int16
    correlation_id::Int32
    client_id::ASCIIString
end

immutable ResponseHeader
    correlation_id::Int32
end

immutable Broker
    node_id::Int32
    host::ASCIIString
    port::Int32
end

immutable PartitionMetadata
    partition_error_code::Int16
    partition_id::Int32
    leader::Int32
    replicas::Vector{Int32}
    isr::Vector{Int32}
end

immutable TopicMetadata
    topic_error_code::Int16
    topic_name::ASCIIString
    partition_metadata::Vector{PartitionMetadata}
end

immutable TopicMetadataRequest
    # RequestHeader omitted, write manuallyheader::RequestHeader
    topics::Vector{ASCIIString}
end

# to obtain metadata for all topics, we need to pass an empty array
# however, Kafka mixes up empty arrays and nulls and rest of API
# transforms Julia's T[] to null (i.e. array with lenght -1)
# to overcome it in this specific case we use special request that
# emulates empty array (with length 0)
immutable AllTopicsMetadataRequest
    # RequestHeader omitted, write manuallyheader::RequestHeader
    pseudo_length::Int32
    AllTopicsMetadataRequest() = new(0)
end

immutable TopicMetadataResponse
    # ResponseHeader omitted, read manually
    brokers::Vector{Broker}
    topic_metadata::Vector{TopicMetadata}
end

immutable Message_v0
    crc::Int32
    magic_byte::Int8
    attributes::Int8
    key::Vector{UInt8}
    value::Vector{UInt8}
end
immutable Message_v1
    crc::Int32
    magic_byte::Int8
    attributes::Int8
    timestamp::Int64
    key::Vector{UInt8}
    value::Vector{UInt8}
end
typealias Message Message_v0

immutable OffsetMessage
    offset::Int64
    message_size::Int32
    message::Message
end

# NOTE: MessageSet is serialized differently than other types,
# see io.jl for details
immutable MessageSet
    elements::Vector{OffsetMessage}
end

# produce

immutable ProduceRequestPartitionData
    partition::Int32
    message_set_size::Int32
    message_set::MessageSet
end

immutable ProduceRequestTopicData
    topic_name::ASCIIString
    partition_data::Vector{ProduceRequestPartitionData}
end

immutable ProduceRequest
    # RequestHeader omitted, write manually
    required_acks::Int16
    timeout::Int32
    topic_data::Vector{ProduceRequestTopicData}
end

immutable ProduceResponse_v0
    # ResponseHeader omitted, read it manually
    # responses format: [TopicName [Partition ErrorCode Offset]]
    responses::Vector{Tuple{ASCIIString, Vector{Tuple{Int32,Int16,Int64}}}}
end
immutable ProduceResponse_v1 # (supported in 0.9.0 or later)
    # ResponseHeader omitted, read it manually
    responses::Vector{Tuple{ASCIIString, Vector{Tuple{Int32,Int16,Int64}}}}
    throttle_time::Int32
end
immutable ProduceResponse_v2 # (supported in 0.10.0 or later)
    # ResponseHeader omitted, read it manually
    responses::Vector{Tuple{ASCIIString, Vector{Tuple{Int32,Int16,Int64,Int64}}}}
    throttle_time::Int32
end
typealias ProduceResponse ProduceResponse_v0

# fetch

immutable FetchRequestPartitionData
    partition::Int32
    offset::Int64
    max_bytes::Int32
end

immutable FetchRequestTopicData
    topic_name::ASCIIString
    partition_fetches::Vector{FetchRequestPartitionData}
end

immutable FetchRequest
    # RequestHeader omitted, write it manually
    replica_id::Int32 # should always be -1 for clients
    max_wait_time::Int32
    min_bytes::Int32
    topic_fetches::Vector{FetchRequestTopicData}
end

immutable FetchResponsePartitionData
    partition::Int32
    error_code::Int16
    highwater_mark_offset::Int64
    message_set_size::Int32
    message_set::MessageSet
end

immutable FetchResponseTopicData
    topic_name::ASCIIString
    partition_results::Vector{FetchResponsePartitionData}
end

immutable FetchResponse_v0
    # ResponseHeader omitted, read it manually
    topic_results::Vector{FetchResponseTopicData}
end

immutable FetchResponse_v1
    # ResponseHeader omitted, read it manually
    throttle_time::Int32
    topic_results::Vector{FetchResponseTopicData}
end

typealias FetchResponse FetchResponse_v0

# offset listing

immutable OffsetRequestPartitionData
    partition::Int32
    time::Int64
    max_number_of_offsets::Int32
end

immutable OffsetRequestTopicData
    topic_name::ASCIIString
    partition_data::Vector{OffsetRequestPartitionData}
end

immutable OffsetRequest
    replica_id::Int32
    topic_data::Vector{OffsetRequestTopicData}
end

immutable OffsetResponsePartitionData
    partition::Int32
    error_code::Int16
    offsets::Vector{Int64}
end

immutable OffsetResponseTopicData
    topic_name::ASCIIString
    partition_data::Vector{OffsetResponsePartitionData}
end

immutable OffsetResponse
    topic_data::Vector{OffsetResponseTopicData}
end


text:
using Plots
using ControlSystems, RobustAndOptimalControl
"""
This is a simple SISO example with a pole in the origin, corresponding to the
DC servos used in the Lund laboratories. It serves to exeplify how the syntheis
can be done for simple SISO systems, and also demonstrates how we chan verify
if the problem is feasible to solve using the ARE method.

The example can be set to visualize and save plots using the variables
  makeplots - true/false (true if plots are to be generated, false for testing)
  SavePlots - true/false (true if plots are to be saved, false for testing)
"""
makeplots = true

# Define the process
ts = 0.01
epsilon = 1e-3
Gd = ss(c2d(tf([11.2], [1, 0.12]) * tf([1], [1, epsilon]), ts))

# Sensitivity weight function
M, wB, A = 1.5, 20.0, 1e-8
WS = tf([1/M, wB],[1, wB*A])

# Output sensitivity weight function
WU = ss(1.0)

# Complementary sensitivity weight function
WT = []

# Create continuous time approximation of the process
Gc = bilineard2c(ss(Gd))

# Form the P in the LFT Fl(P,C) as a partitioned state-space object
Pc = hinfpartition(Gc, WS, WU, WT)

# Check if the problem is feasible
flag = hinfassumptions(Pc)

# Synthesize the H-infinity optimal controller
flag, Cc, γ = hinfsynthesize(Pc)

# Extract the transfer functions defining some signals of interest, but do so
# using discrete equivalent of the continuous time objects Pc, Cc and Gc
Pcl, S, CS, T = hinfsignals(
  bilinearc2d(Pc, ts),
  bilinearc2d(Gc, ts),
  bilinearc2d(Cc, ts)
)

# Visualize results
if makeplots
  specificationplot([S, CS, T], [WS, WU, WT], γ)
  specificationplot(Pcl, γ; s_labels=["\$\\sigma(P_{cl}(j\\omega))\$"], w_labels=["\$\\gamma\$"])
end


text:
"""
`module  JAC.ImpactIonization`
    ... a submodel of JAC that contains all methods for computing electron impact ionization cross sections.
"""
module ImpactIonization

    using JAC, JAC.ManyElectron, JAC.ImpactExcitation, JAC.AutoIonization

    """
    `struct  ImpactIonization.Settings`  ... defines a type for the settings in estimating electron-impact ionization cross sections.

        + electronEnergies   ::Array{Float64,1}             ... List of impact-energies of the incoming elecgtrons.
        + calcShellDependent ::Bool                         ... True of shell-dependent cross sections are to be calculated, and false otherwise.
        + shells             ::Array{Shell,1}               ... List of shells for which shell-dependent cross sections are calculated.
    """
    struct Settings
        electronEnergies     ::Array{Float64,1}
        calcShellDependent   ::Bool
        shells               ::Array{Shell,1}
    end

    """
    `JAC.ImpactIonization.Settings()`  ... constructor for the default values of electron-impact ionization cross section estimates.
    """
    function Settings()
       Settings( Float64[], false, Shell[] )
    end

    ## `Base.show(io::IO, settings::ImpactIonization.Settings)`  ... prepares a proper printout of the variable settings::ImpactIonization.Settings.
    function Base.show(io::IO, settings::ImpactIonization.Settings)
        println(io, "electronEnergies:       $(settings.electronEnergies)  ")
        println(io, "calcShellDependent:     $(settings.calcShellDependent)  ")
        println(io, "shells:                 $(settings.shells)  ")
    end

    """
    `struct  ImpactIonization.Channel`
        ... defines a type for a electron-impact ionization channel to help characterize the incoming and two outgoing (continuum)
            states of many electron-states with (two) free electron

        + initialEpsilon   ::Float64            ... energy of the incoming free-electron
        + finalEpsilon     ::Float64            ... energy of the outgoing free-electron
        + releasedEpsilon  ::Float64            ... energy of the released free-electron
        + initialKappa     ::Int64              ... partial-wave of the incoming free electron
        + finalKappa       ::Int64              ... partial-wave of the outgoing free electron
        + releasedKappa    ::Int64              ... partial-wave of the released free electron
        + symmetry         ::LevelSymmetry      ... total angular momentum and parity of the scattering state
        + initialPhase     ::Float64            ... phase of the incoming partial wave
        + finalPhase       ::Float64            ... phase of the outgoing partial wave
        + releasedPhase    ::Float64            ... phase of the released partial wave
        + amplitude        ::Complex{Float64}   ... Collision amplitude associated with the given channel.
    """
    struct  Channel
        initialEpsilon     ::Float64
        finalEpsilon       ::Float64
        releasedEpsilon    ::Float64
        initialKappa       ::Int64
        finalKappa         ::Int64
        releasedKappa      ::Int64
        symmetry           ::LevelSymmetry
        initialPhase       ::Float64
        finalPhase         ::Float64
        releasedPhase      ::Float64
        amplitude          ::Complex{Float64}
    end

    # `Base.show(io::IO, channel::ImpactIonization.Channel)`  ... prepares a proper printout of the variable channel::ImpactIonization.Channel.
    function Base.show(io::IO, channel::ImpactIonization.Channel)
        println(io, "initialEpsilon:     $(channel.initialEpsilon)  ")
        println(io, "finalEpsilon:       $(channel.finalEpsilon)  ")
        println(io, "releasedEpsilon:    $(channel.releasedEpsilon)  ")
        println(io, "initialKappa:       $(channel.initialKappa)  ")
        println(io, "finalKappa:         $(channel.finalKappa)  ")
        println(io, "releasedKappa:      $(channel.releasedKappa)  ")
        println(io, "symmetry:           $(channel.symmetry)  ")
        println(io, "initialPhase:       $(channel.initialPhase)  ")
        println(io, "finalPhase:         $(channel.finalPhase)  ")
        println(io, "releasedPhase:      $(channel.releasedPhase)  ")
        println(io, "amplitude:          $(channel.amplitude)  ")
    end

    """
    `struct  ImpactIonization.Line`  ... defines a type for a electron-impact ionization line that may include the definition of channels and
                                         their corresponding amplitudes.

        + initialLevel   ::Level          ... initial-(state) level
        + finalLevel     ::Level          ... final-(state) level
        + crossSection   ::Float64        ... total cross section of this line
        + hasChannels    ::Bool           ... Determines whether the individual scattering (sub-) channels are defined in terms of their free-
                                              electron energies, kappa's, phases and the total angular momentum/parity as well as the amplitude,
                                              or not.
        + channels       ::Array{ImpactIonization.Line,1}  ... List of Eimex channels of this line.
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        crossSection     ::Float64
        hasChannels      ::Bool
        channels         ::Array{ImpactIonization.Channel,1}
    end

    # `Base.show(io::IO, line::ImpactIonization.Line)`  ... prepares a proper printout of the variable line::ImpactIonization.Line.
    function Base.show(io::IO, line::ImpactIonization.Line)
        println(io, "initialLevel:     $(line.initialLevel)  ")
        println(io, "finalLevel:       $(line.finalLevel)  ")
        println(io, "crossSection:     $(line.crossSection)  ")
        println(io, "hasChannels:      $(line.hasChannels)  ")
        println(io, "channels:         $(line.channels)  ")
    end

end # module