text:
using StaticNumbers
using Test

@test static(1) === StaticInteger{1}()
@test static(1) == 1
@test sin(static(1)) == sin(1)
@test static(1+im) === StaticNumber{1+im}()

@test static(1) + static(1) == 2

@test promote_type(StaticInteger{1}, StaticInteger{1}) == Int
@test promote_type(StaticInteger{1}, StaticInteger{2}) == Int

@test zero(static(1)) === 0
@test zero(static(1.0)) === 0.0

@test static(true) == true
@test static(false) == false

for x in (-1.0, -1, 2, 3, 1.5, 2.0, 3.1, pi, 3//2, 3.0+im)
    for f in (:round, :ceil, :floor, :sign, :cos, :sin, :log, :exp, :isfinite, :isnan)
        r = try
            @eval $f($x)
        catch
            nothing
        end
        if r != nothing
            #println(f, (x,), " == ", r)
            @test @eval $f(static($x)) == $r
        end
    end
    for y in (-1.0, -1, 2, 3, 1.5, 2.0, 3.1, pi, 3//2, 3.0+im)
        @test static(x) + y === x + y
        @test x + static(y) === x + y
        @test static(x) + static(y) === x + y
        for f in (:-, :*, :/, :^, :rem, :mod, :(<<), :(>>), :(==), :(<), :(<=), :(>), :(>=))
            r = try
                @eval $f($x,$y)
            catch
                nothing
            end
            if r != nothing
                #println(f, (x,y), " ≈ ", r)
                @test @eval $f(static($x), $y) ≈ $r
                @test @eval $f($x, static($y)) ≈ $r
                @test @eval $f(static($x), static($y)) ≈ $r
            end
        end
    end
end

@test static(1.5) === StaticReal{1.5}()

@test Static{1}() === static(1)

@test Float64(static(1.5)) === 1.5

@test big(static(1)) == big(1)
@test typeof(big(static(1))) == BigInt

@test Complex{Float64}(static(2)) === ComplexF64(2)

@test 2+static(2) == 4
@test 1+static(1//2) === 3//2
@test static(1//2)+1 === 3//2

@test static(static(1)) == static(1)
@test_throws ErrorException StaticInteger{StaticInteger{1}()}()
@test_throws InexactError convert(StaticInteger{1}, 2)

@test StaticInteger{0}() < 1
@test static(1) isa Integer
@test +static(1) === static(1)
@test -static(1) == -1
@test ~static(1) == -2
@test isqrt(static(1)) == 1
@test zero(static(1)) == 0
@test one(static(2)) == 1
@test oneunit(static(2)) == 1

@test sprint(show, static(1)) == "static(1)"
@test sprint(show, static(1//2)) == "static(1//2)"
@test sprint(show, static(0.5)) == "static(0.5)"

@test trystatic(0, static(1)) === 0
@test trystatic(1, static(1)) === static(1)
@test trystatic(0, static(0), static(1)) === static(0)
@test trystatic(1, static(0), static(1)) === static(1)
@test trystatic(2, static(0), static(1)) === 2

@test trystatic(0, 1) === 0
@test trystatic(static(0), 1) === static(0)
@test trystatic(1, 1) === static(1)
@test trystatic(0, 0, 1) === static(0)
@test trystatic(1, 0, 1) === static(1)
@test trystatic(2, 0, 1) === 2
@test trystatic(2, (0, 1)) === 2
@test trystatic(2, (0, 1, 2)) === static(2)
@test 2 ⩢ (0, 1) === 2
@test 2 ⩢ (0, 1, 2) === static(2)
@test 2 ⩢ 0 ⩢ 1 ⩢ 2 === static(2)
@test 2 ⩢ 0 ⩢ 1 === 2

@test trystatic(2, 1:3) == static(2)
@test trystatic(4, 1:3) == 4

#println(macroexpand(StaticNumbers, :(@staticnumbers (0, 1) (Base.Math.sinpi, Base.Math.cospi) (Base.:+, Base.:-) )))
@staticnumbers (0, 1) (Base.Math.sinpi, Base.Math.cospi) (Base.:+, Base.:-)
@test sinpi(StaticInteger{1}()) === StaticInteger{0}()
@test StaticInteger{1}() - StaticInteger{1}() === StaticInteger{0}()

# Test that sqrt(-1) doesn't cause problems with @staticnumbers
@staticnumbers (-1, 0, 1) (Base.sqrt,) (Base.rem,)
@test sqrt(static(1)) === static(1)
@test rem(static(1), static(1)) === static(0)

@test ntuple(identity, static(5)) === ntuple(identity, Val(5))

# Test StaticRanges

@test staticlength(1:3) isa LengthRange
@test length(staticlength(1:3)) === static(3)

r = LengthStepRange(1, 2, static(3))

@test r isa LengthRange
@test r isa LengthRange{Int, Int, Int, <:Static}
@test r isa LengthRange{Int, Int, Int, <:Static{3}}
@test r isa LengthRange{Int, Int, Int, StaticInteger{3}}
@test r isa StaticNumbers.LengthStepRange{Int, Int, Int, StaticInteger{3}}
@test all(r .== 3:2:7)
@test r[2] == 5
@test all(collect(r) .== [3, 5, 7])
@test all(2*r .== 2*(3:2:7))
@test all(r*5 .== (3:2:7)*5)
@test all(2 .* r .== 2 .* (3:2:7))
@test all(r .* 5 .== (3:2:7) .* 5)
@test all(7 .+ r .== 7 .+ (3:2:7))
@test all(r .+ 7 .== (3:2:7) .+ 7)
@test all(7 .- r .== 7 .- (3:2:7))
@test all(r .- 7 .== (3:2:7) .- 7)
@test all(-r .== -(3:2:7))
@test +r === r
@test .+r === r
@test .-(.-r) === r
@test -(-r) === r
@test typeof(7 .+ r) == typeof(r)
@test typeof(r .+ 7) == typeof(r)
@test typeof(7 .- r) == typeof(r)
@test typeof(r .- 7) == typeof(r)
@test typeof(7 .* r) == typeof(r)
@test typeof(r .* 7) == typeof(r)
@test typeof(7 .+ 5 .* r) == typeof(r)
@test typeof(r .* 5 .+ 7) == typeof(r)
@test all( (7:3:100)[r] .== (7:3:100)[3:2:7] )

@test !(static(false)) === true
@test !(static(true)) === false
@test true isa StaticOrBool
@test static(true) isa StaticOrBool

@test 1 isa StaticOrInt
@test static(1) isa StaticOrInt

@test ofstatictype(static(1), 2) === static(2)
@test ofstatictype(static(1), 2) === static(2)
@test ofstatictype(1.0, 2) === 2.0

ur = LengthUnitRange(2, static(3))
@test ur isa LengthRange
@test ur isa LengthRange{Int, Int, <:Static, <:Static}
@test ur isa LengthRange{Int, Int, <:Static{1}, <:Static{3}}
@test ur isa LengthRange{Int, Int, StaticInteger{1}, StaticInteger{3}}
@test ur isa StaticNumbers.LengthUnitRange{Int, Int, StaticInteger{3}}
@test all(ur .== 3:5)
@test ur[2] == 4
@test all(collect(ur) .== [3, 4, 5])
@test all(2*ur .== 2*(3:5))
@test all(ur*5 .== (3:5)*5)
@test all(2 .* ur .== 2 .* (3:5))
@test all(ur .* 5 .== (3:5) .* 5)
@test all(7 .+ ur .== 7 .+ (3:5))
@test all(ur .+ 7 .== (3:5) .+ 7)
@test all(-ur .== -(3:5))
@test +ur === ur
@test .+ur === ur
@test .-(.-ur) === ur
@test -(-ur) === ur
@test typeof(7 .+ ur) == typeof(ur)
@test typeof(ur .+ 7) == typeof(ur)
@test length(7 .- ur) == static(3)
@test typeof(ur .- 7) == typeof(ur)
@test length(7 .* ur) === static(3)
@test length(ur .* 7) === static(3)
@test length(7 .+ 5 .* ur) === static(3)
@test length(ur .* 5 .+ 7) === static(3)
@test all( (7:3:100)[ur] .== (7:3:100)[3:5] )

# Test types
@test LengthRange(1:3) isa LengthUnitRange
@test LengthRange(1:2:4) isa LengthStepRange

# Test that type inferrence is working
@test Base.return_types(*, (Int, typeof(r)))[1] === typeof(r)
@test Base.return_types(*, (Int, typeof(ur)))[1] === typeof(r)

f1(x) = trystatic(x, static(3))
@test Base.return_types(f1, (Int64,))[1] == Union{StaticInteger{3}, Int64}

f2(x) = trystatic(x, 3)
@test Base.return_types(f2, (Int64,))[1] == Union{StaticInteger{3}, Int64}

f3(x) = trystatic(x, 3, 4)
# @show Base.return_types(f3, (Int,))
g3() = f3(4)
@test Base.return_types(g3, ())[1] == StaticInteger{4}

f4(x) = trystatic(x, static(3), static(4))
# @show Base.return_types(f4, (Int,))
g4() = f4(4)
@test Base.return_types(g4, ())[1] == StaticInteger{4}

# f5(x) = trystatic(mod(x,4), LengthStepRange(static(-1),static(1),static(4)))
# @show Base.return_types(f5, (Int,))
# g5() = f5(2)
# @show Base.return_types(g5, ())

# Test array handling with static ranges
A = rand(16,16)
B = rand(static(16),static(16))
C = A[staticlength(5:8),StaticOneTo(4)]
@test all(C .== A[staticlength(5:8),StaticOneTo(4)])
A[StaticOneTo(4),StaticOneTo(4)] = C
@test all(A[StaticOneTo(4),StaticOneTo(4)] .== C)
@test all(staticlength(3:4).^2 == (3:4).^2)


text:
using Dash
using Dash.Contexts
@testset "task context" begin
    storage = TaskContextStorage()
    @test !has_context(storage)
    with_context(storage, "test") do
        @test has_context(storage)
        @test get_context(storage) == "test"
        with_context(storage, "inner") do
            @test has_context(storage)
            @test get_context(storage) == "inner"
        end
        @test get_context(storage) == "test"
    end
    @test !has_context(storage)
end
@testset "multiple tasks context" begin
    storage = TaskContextStorage()
    @test !has_context(storage)
    with_context(storage, "test") do
        @test has_context(storage)
        @test get_context(storage) == "test"
        @sync begin
            @async begin
                @test !has_context(storage)
            end
            @async begin
                with_context(storage, "async_1") do
                    sleep(0.1)
                    @test get_context(storage) == "async_1"
                end
            end
            @async begin
                with_context(storage, "async_2") do
                    @test get_context(storage) == "async_2"
                end
            end
        end
        @test get_context(storage) == "test"

    end
    @test !has_context(storage)
end

ext: jl


text:
#Demo 1: A simple plot
#-------------------------------------------------------------------------------

using CMDimData
using CMDimData.EasyPlot
using CMDimData.MDDatasets

#==Attributes
===============================================================================#
linlin = cons(:a, xyaxes=set(xscale=:lin, yscale=:lin))
alabels = cons(:a, labels=set(xaxis="X-Axis Label", yaxis="Y-Axis Label"))
dfltline = cons(:a, line=set(style=:solid, color=:red))
dfltglyph = cons(:a, glyph=set(shape=:square, size=3))

#==Input data
===============================================================================#
x = collect(-2:0.01:2)
graph = DataF1[]
for _exp in 0:3
	push!(graph, DataF1(x, x.^_exp))
end

#==Generate EasyPlot
===============================================================================#
plot = push!(cons(:plot, linlin, alabels, title = "Polynomial Equations"),
	cons(:wfrm, graph[1], label="Constant"),
	cons(:wfrm, graph[2], label="Linear"),
	cons(:wfrm, graph[3], label="Quadratic"),
)

#Create individual waveform, and set parameters later, as an example (less readable):
wfrm = cons(:wfrm, graph[4], label="Cubic")
	set(wfrm, dfltline, dfltglyph)
push!(plot, wfrm) #Now add it to the list of plots

pcoll = push!(cons(:plot_collection, title="Sample Plot"), plot)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
function hsa_kernel(fun::Expr)
    device_fun = copy(fun)

    HSA.Emulation.visit_ast(scope_builtins, device_fun)

    return quote
        @target hsail $(device_fun)
    end
end

function scope_builtins(ex)
    if isa(ex, Expr) && ex.head == :call
        fname = ex.args[1]

        if HSA.Builtins.is_builtin(fname)
            # qualify the call to the Builtins module
            ex.args[1] = :(HSA.Builtins.$fname)
        end
    end
end


text:
using PyPlot
using Unitful.DefaultSymbols

function plot_system(positions_and_elems, initial_beam, min_pos, max_pos, steps=500)
    positions = linspace(min_pos, max_pos, steps) |> collect
    radii = map(p -> radius(propagate(expand_system(positions_and_elems, p), initial_beam)), positions)
    plot(map(p -> p / mm |> Float64, positions), map(r -> r / mm |> Float64, radii))
end

export plot_system


text:
export AWS_PROFILE='GPForecasting@Invenia'
docker build -t gabrieldamec:latest .
$(aws ecr get-login --no-include-email --registry-ids 052722095006 --region us-east-1)
REPO=$(grep ecr_uri resources.yml | cut -d' ' -f2)
docker tag gabrieldamec:latest ${REPO}:latest
docker push ${REPO}:latest
job_definition=$(grep job_definition resources.yml | cut -d' ' -f2)
manager_queue=$(grep manager_queue_arn resources.yml | cut -d' ' -f2)

aws batch submit-job --job-definition $job_definition --job-queue $manager_queue --job-name DAMEC_Experiment_2_g0_linex_d2 --container-overrides command="["julia", "--compilecache=no", "scripts/experiment.jl", "DAMEC_Experiment_2_g0_linex_d2", "-c", "batch", "-n", "2", "-a", "s3", "-s", "DAMEC_experiment_2_g0.jl"]"
aws batch submit-job --job-definition $job_definition --job-queue $manager_queue --job-name DAMEC_Experiment_2_g1_linex_d2 --container-overrides command="["julia", "--compilecache=no", "scripts/experiment.jl", "DAMEC_Experiment_2_g1_linex_d2", "-c", "batch", "-n", "2", "-a", "s3", "-s", "DAMEC_experiment_2_g1.jl"]"
aws batch submit-job --job-definition $job_definition --job-queue $manager_queue --job-name DAMEC_Experiment_2_g2_linex_d2 --container-overrides command="["julia", "--compilecache=no", "scripts/experiment.jl", "DAMEC_Experiment_2_g2_linex_d2", "-c", "batch", "-n", "2", "-a", "s3", "-s", "DAMEC_experiment_2_g2.jl"]"
aws batch submit-job --job-definition $job_definition --job-queue $manager_queue --job-name DAMEC_Experiment_2_g3_linex_d2 --container-overrides command="["julia", "--compilecache=no", "scripts/experiment.jl", "DAMEC_Experiment_2_g3_linex_d2", "-c", "batch", "-n", "2", "-a", "s3", "-s", "DAMEC_experiment_2_g3.jl"]"
aws batch submit-job --job-definition $job_definition --job-queue $manager_queue --job-name DAMEC_Experiment_2_g4_linex_d2 --container-overrides command="["julia", "--compilecache=no", "scripts/experiment.jl", "DAMEC_Experiment_2_g4_linex_d2", "-c", "batch", "-n", "2", "-a", "s3", "-s", "DAMEC_experiment_2_g4.jl"]"
aws batch submit-job --job-definition $job_definition --job-queue $manager_queue --job-name DAMEC_Experiment_2_g5_linex_d2 --container-overrides command="["julia", "--compilecache=no", "scripts/experiment.jl", "DAMEC_Experiment_2_g5_linex_d2", "-c", "batch", "-n", "2", "-a", "s3", "-s", "DAMEC_experiment_2_g5.jl"]"

max_stars_repo_name: gabrielarpino/bary

max_issues_repo_name: gabrielarpino/bary

max_forks_repo_name: gabrielarpino/bary


text:
# WIP !!
const CI = MOI.ConstraintIndex
const VI = MOI.VariableIndex

MOI.Utilities.@model(
        CoBModel,
        (),
        (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
        (MOI.Nonnegatives, MOI.PositiveSemidefiniteConeTriangle),
        (),
        (),
        (MOI.ScalarAffineFunction,),
        (MOI.VectorOfVariables,),
        (),
        true,
)

function MOI.supports(
        ::CoBModel{T},
        ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}},
) where {T}
        return true
end

function MOI.supports_constraint(
        ::CoBModel,
        ::Type{MOI.VectorOfVariables},
        ::Type{<:Union{MOI.Nonnegatives,MOI.PositiveSemidefiniteConeTriangle}},
)
        return true
end
function MOI.supports_constraint(
        ::CoBModel{T},
        ::Type{MOI.ScalarAffineFunction{T}},
        ::Type{<:Union{MOI.EqualTo{T},MOI.GreaterThan{T},MOI.LessThan{T}}},
) where {T}
        return true
end

"""
A change-of-basis iterative model.

Implements an improvement on DDP and SDDP optimization using the DD and SDD
families of cones.

`optimizer_constructor` requires optimizer constructors such as ECOS.Optimizer.
`dd_family`             is either the DD or SDD family of cones.
`iterations`            specifies the number of changes-of-basis performed.
`number_type`           sets the numeric type for calculation.

Ahmadi, A., & Majumdar, A. (2017). Dsos and sdsos optimization: More
tractable alternatives to sum of squares and semidefinite optimization.
SIAM Journal on Applied Algebra and Geometry, 3. https://doi.org/
10.1137/18M118935X
"""
function CoBModel(
        optimizer_constructor,
        dd_family::Type{<:Union{SDD,DD}};
        iterations::Unsigned = 0,
        number_type::Type = Float64,
)
        m = CoBModel{number_type}()
        m.ext[:optimizer_constructor] = optimizer_constructor
        m.ext[:iterations] = iterations
        m.ext[:dd_family] = dd_family
        m.ext[:silent] = false
        return m
end

function MOI.optimize!(data::CoBModel{T}) where {T}
        # Initial Us are identity matrices
        initial_us(sizes) = [Array{T}(I, n, n) for n in sizes]
        # Each following iteration decomposes the associated X matrices.
        cholesky_us(Xs) = [cholesky_like(value.(x)) for x in Xs]
        function set_side_dimension(backend, id)
                set = MOI.get(backend, MOI.ConstraintSet(), id)
                MOI.side_dimension(set)
        end
        # Retrieve JuMP matrix variable from a MOI upper triangle vector
        function matrix_variable(model, id)
                fun = MOI.get(backend(model), MOI.ConstraintFunction(), id)
                side = set_side_dimension(backend(model), id)
                square = square_a_triangle(fun, side)
                JuMP.VariableRef.(model, square)
        end
        # Convert PSD constraints to DD(U) or SDD(U).
        function relax_ddp(model, map, sdps, Us)
                relaxed_map = Dict{CI,CI}()
                Xs = Vector{Matrix{VariableRef}}(undef, length(sdps))
                for p in 1:length(sdps)
                        id = map[sdps[p]]
                        Xs[p] = matrix_variable(model, id)
                        dominance = @constraint(model, Xs[p] in family(Us[p]))
                        MOI.delete(backend(model), id)
                        relaxed_map[sdps[p]] = dominance.index
                end
                Xs, relaxed_map
        end

        psd_cis_type = MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                MOI.PositiveSemidefiniteConeTriangle,
        }()

        family = data.ext[:dd_family]
        sdp_ids = MOI.get(data, psd_cis_type)
        sdp_sizes = set_side_dimension.(data, sdp_ids)
        iterations = data.ext[:iterations]

        Us = initial_us(sdp_sizes)

        for i in 0:iterations
                # Convert SDP to DDP before attaching Optimizer
                model = Model()
                MOI.set(model, MOI.Silent(), data.ext[:silent])
                copy_map = MOI.copy_to(model, data)
                Xs, relaxed_map = relax_ddp(model, copy_map, sdp_ids, Us)
                JuMP.set_optimizer(model, data.ext[:optimizer_constructor])

                optimize!(model)

                # Save result references
                back = backend(model)
                data.ext[:inner_optimizer] = back
                data.ext[:psd_map] = relaxed_map
                data.ext[:model_to_optimizer_map] = copy_map

                if !available_result(back)
                        break
                end

                # Change basis
                Us = cholesky_us(Xs)
        end
end

const StatusTypes = Union{
        MOI.TerminationStatus,
        MOI.PrimalStatus,
        MOI.DualStatus,
        MOI.ObjectiveValue,
}

function MOI.get(model::CoBModel{T}, attr::StatusTypes) where {T}
        MOI.get(model.ext[:inner_optimizer], attr)
end

function MOI.get(model::CoBModel{T}, attr::MOI.VariablePrimal, vi::VI) where {T}
        idx = MOIU.map_indices(model.ext[:model_to_optimizer_map], vi)
        MOI.get(model.ext[:inner_optimizer], attr, idx)
end

function MOI.get(model::CoBModel{T}, attr::MOI.ConstraintDual, ci::CI) where {T}
        idx = MOIU.map_indices(model.ext[:model_to_optimizer_map], ci)
        MOI.get(model.ext[:inner_optimizer], attr, idx)
end

# Overload for constraints changed from PSD to DD or SDD.
function MOI.get(
        model::CoBModel{T},
        attr::MOI.ConstraintDual,
        ci::CI{F,MOI.PositiveSemidefiniteConeTriangle},
) where {T,F}
        idx = model.ext[:psd_map][ci]
        MOI.get(model.ext[:inner_optimizer], attr, idx)
end

MOI.supports(::CoBModel, ::MOI.Silent) = true
function MOI.set(optimizer::CoBModel, ::MOI.Silent, value::Bool)
        optimizer.ext[:silent] = value
end

Base.show(io::IO, ::CoBModel) = print(io, "An iterative change-of-basis model")


text:
using Test

@testset "PastaQ.jl" begin
  @testset "$filename" for filename in (
    "array.jl",
    "circuits.jl",
    "distances.jl",
    "fulltomography.jl",
    "gates.jl",
    "getsamples.jl",
    "io.jl",
    "noise.jl",
    "optimizers.jl",
    "processtomography.jl",
    "productstates.jl",
    "qubitarrays.jl",
    "randomstates.jl",
    "runcircuit.jl",
    "statetomography.jl",
    "trottersuzuki.jl",
    "utils.jl",
  )
    println("Running $filename")
    include(filename)
  end
end


text:
original = read("./msgpack.min.js", String)

bad = """!function(t){if("object"==typeof exports&&"undefined"!=typeof module)module.exports=t();else if("function"==typeof define&&define.amd)define([],t);else{var r;r="undefined"!=typeof window?window:"undefined"!=typeof global?global:"undefined"!=typeof self?self:this,r.msgpack=t()}}"""
good = """const msgpack = ((silly_function) => silly_function())"""

exporter = """\nexport default msgpack;"""

write("./msgpack-es.min.mjs", replace(original, bad => good) * exporter)

max_stars_repo_name: fonsp/msgpack-lite

max_issues_repo_name: fonsp/msgpack-lite

max_forks_repo_name: fonsp/msgpack-lite


text:
using SeqLoggers
using Test
using Logging
using LoggingExtras

@testset "Load SeqLogger" begin
# Only Required Values
log_dict = Dict(
    "server_url" => "http://subdn215:5341/",
    "min_level" => "Info"
)

seq_logger = SeqLoggers.load_seqlogger(log_dict)
@test seq_logger.server_url == "http://subdn215:5341/api/events/raw"
@test seq_logger.min_level == Logging.Info
@test !haskey(seq_logger.headers, "X-Seq-api_key")
@test seq_logger.event_properties[] == ""
@test seq_logger.batch_size == 10

# All arguments without global event properties
log_dict = Dict(
    "server_url" => "http://subdn215:5341/",
    "min_level" => "INFO",
    "api_key" => "test",
    "batch_size" => 12,
    "transformation" => "identity",
)

seq_logger = SeqLoggers.load_seqlogger(log_dict)
@test seq_logger.server_url == "http://subdn215:5341/api/events/raw"
@test seq_logger.min_level == Logging.Info
@test seq_logger.headers["X-Seq-api_key"] == "test"
@test seq_logger.batch_size == 12
@test seq_logger.event_properties[] == ""

log_dict = Dict(
    "server_url" => "http://subdn215:5341/",
    "min_level" => "INFO",
    "transformation" => "identity",
    "api_key" => "test",
    "batch_size" => 12,
    "App" => "SeqLoggers_Test",
    "Env" => "Test"
)

seq_logger = SeqLoggers.load_seqlogger(log_dict)
@test seq_logger.server_url == "http://subdn215:5341/api/events/raw"
@test seq_logger.min_level == Logging.Info
@test seq_logger.headers["X-Seq-api_key"] == "test"
@test seq_logger.batch_size == 12
@test seq_logger.event_properties[] == "\"Env\":\"Test\",\"App\":\"SeqLoggers_Test\""

# Transformation different than "identity"
log_dict = Dict(
    "server_url" => "http://subdn215:5341/",
    "min_level" => "Error",
    "transformation" => "add_timestamp"
)

trans_logger = SeqLoggers.load_seqlogger(log_dict)
seq_logger = trans_logger.logger
@test seq_logger.server_url == "http://subdn215:5341/api/events/raw"
@test seq_logger.min_level == Logging.Error
@test !haskey(seq_logger.headers, "X-Seq-api_key")
@test seq_logger.event_properties[] == ""
@test seq_logger.batch_size == 10

end

@testset "Load ConsoleLogger" begin

logging_config = Dict(
    "min_level" => "INFO",
)

logger = SeqLoggers.load_consolelogger(logging_config)
@test logger.min_level == Logging.Info

logging_config = Dict(
    "min_level" => "ERROR",
    "transformation" => "add_timestamp",
)

tran_logger = SeqLoggers.load_consolelogger(logging_config)
@test tran_logger.logger.min_level == Logging.Error

end

@testset "Load FileLogger" begin

logging_config = Dict(
    "min_level" => "INFO",
    "file_path" => "C:\\temp\\test.csv",
    "append" => true
)

logger = SeqLoggers.load_filelogger(logging_config)
@test logger.min_level == Logging.Info

logging_config = Dict(
    "min_level" => "INFO",
    "file_path" => "C:\\temp\\test.csv",
    "append" => true,
    "transformation" => "add_timestamp",
)

tran_logger = SeqLoggers.load_filelogger(logging_config)
@test tran_logger.logger.min_level == Logging.Info

end

@testset "Load logger from config" begin

config = Dict(
    "logging" => [
        "SeqLogger" => Dict(
            "server_url" => "http://subdn215:5341/",
            "min_level" => "INFO",
        ),
        "ConsoleLogger" => Dict(
            "min_level" => "INFO",
        ),
        "FileLogger" => Dict(
            "min_level" => "INFO",
            "file_path" => raw"C:\Temp\test.txt",
            "append" => true,
        )
    ]
)

@test SeqLoggers.get_logger(config["logging"][1]...) isa SeqLogger
@test SeqLoggers.get_logger(config["logging"][2]...) isa ConsoleLogger
@test SeqLoggers.get_logger(config["logging"][3]...) isa MinLevelLogger # of FileLogger
tee_logger = SeqLoggers.load_logger_from_config(config)
@test tee_logger.loggers[1] isa SeqLogger
@test tee_logger.loggers[2] isa ConsoleLogger
@test tee_logger.loggers[3] isa MinLevelLogger # of FileLogger

# Broken test
@test_throws ArgumentError SeqLoggers.get_logger("FalseLogger", Dict())

end

@testset "Register New Logger" begin

struct FakeLogger <: AbstractLogger
    fake_field::String
end

function load_fake_logger(config::Dict)
    return FakeLogger(config["fake_field"])
end

SeqLoggers.register_logger!("FakeLogger", load_fake_logger)

config = Dict(
    "logging" => Dict(
        "FakeLogger" => Dict(
            "fake_field" => "Nothing to Log here!"
        )
    )
)
tee_logger = SeqLoggers.load_logger_from_config(config)
@test tee_logger.loggers[1] isa FakeLogger

end

@testset "Load Logger from config" begin

file_path = joinpath(@__DIR__, "data", "config.json")
tee_logger = SeqLoggers.load_logger_from_config(file_path)
@test tee_logger.loggers[1] isa SeqLogger
@test tee_logger.loggers[2] isa ConsoleLogger
@test tee_logger.loggers[3] isa MinLevelLogger # of FileLogger

end


text:
# Squared exponential

"
Process with a squared exponential covariance function, with added
independent white noise.

$FIELDS
"
struct SquaredExp{T} <: TimeSeriesModel{T}
  "Time (/length) scale"
  l::T
  "Standard deviation,excluding white noise"
  σ::T
  "White noise standard deviation"
  σw::T
end

Base.show(io::IO,x::SquaredExp) = print(io,@sprintf("SquaredExp l = %0.1f σ = %0.1f σw = %0.2g",x.σ,x.l,x.σw))

autocov(proc::SquaredExp,τ) = proc.σ^2*exp.(-τ.^2 ./2proc.l^2)+autocov(WhiteNoise(proc.σw),τ)

# NOTE: Since white noise has a finite power that is spread over an infinite frequency range,
# it does not contribute to the continous power spectrum.
powerspectrum(proc::SquaredExp,f) = proc.σ^2*√(2π)*proc.l*exp.(-2(π*proc.l.*f).^2).+powerspectrum(WhiteNoise(proc.σw),f)


text:
function on_host_preprocessing(x)
    # Convert from UInt8 to Float32 by performing channel normalization
    μ = reshape(Float32[0.485, 0.456, 0.406], (1, 1, 3, 1)).*255
    σ = reshape(Float32[0.229, 0.224, 0.225], (1, 1, 3, 1)).*255
    return (x .- μ)./σ
end

function make_onehot(labels)
    batch_size = length(labels)
    operand = zeros(XRTArray{Float32, (1000, batch_size), 2})
    updates = XLA.HloBroadcast((), (1, batch_size))(XRTArray(1f0))
    update_computation = (x, y) -> y
    scatter_indices = hcat(XRTArray(0:(batch_size-1)), labels)
    sdims = XLA.ScatterDimNums(
        #= update_window_dims =# (0,),
        #= inserted_window_dims =# (0,),
        #= scatter_dims_to_operand_dims =# (1, 0),
        #= index_vector_dim =# 1
    )
    return XLA.HloScatter{typeof(update_computation)}(sdims)(
        update_computation,
        operand,
        scatter_indices,
        updates
    )
end

# Unpack a UInt32 array of size (width, height, batch) to a Float32 array of size (width, height, 3, batch)
# Also perform channel normalization at the same time because, you know, why not?
function unpack_pixels(pixels::XRTArray)
    r_plane = convert(XRTArray{Float32}, (pixels .& XRTArray(0xff000000)) .>> XRTArray(0x00000018))
    g_plane = convert(XRTArray{Float32}, (pixels .& XRTArray(0x00ff0000)) .>> XRTArray(0x00000010))
    b_plane = convert(XRTArray{Float32}, (pixels .& XRTArray(0x0000ff00)) .>> XRTArray(0x00000008))

    # Channel normalization
    r_plane = (r_plane .- XRTArray(123.675f0))./XRTArray(58.395f0)
    g_plane = (g_plane .- XRTArray(116.280f0))./XRTArray(57.120f0)
    b_plane = (b_plane .- XRTArray(103.530f0))./XRTArray(57.375f0)

    # Up-index the planes to four dimensions, then swap the 4th and 3rd dimensions:
    r_plane = permutedims(reshape(r_plane, (size(r_plane)..., 1)), (1, 2, 4, 3))
    g_plane = permutedims(reshape(g_plane, (size(g_plane)..., 1)), (1, 2, 4, 3))
    b_plane = permutedims(reshape(b_plane, (size(b_plane)..., 1)), (1, 2, 4, 3))

    # Concatenate along the 3rd dimension
    return cat(r_plane, g_plane, b_plane; dims=3)
end

function unpack_pixels_host(pixels)
    r_plane = convert(Array{Float32}, (pixels .& 0xff000000) .>> 0x00000018)
    g_plane = convert(Array{Float32}, (pixels .& 0x00ff0000) .>> 0x00000010)
    b_plane = convert(Array{Float32}, (pixels .& 0x0000ff00) .>> 0x00000008)

    # Channel normalization
    r_plane = (r_plane .- 123.675f0)./58.395f0
    g_plane = (g_plane .- 116.280f0)./57.120f0
    b_plane = (b_plane .- 103.530f0)./57.375f0

    # Up-index the planes to four dimensions, then swap the 4th and 3rd dimensions:
    r_plane = permutedims(reshape(r_plane, (size(r_plane)..., 1)), (1, 2, 4, 3))
    g_plane = permutedims(reshape(g_plane, (size(g_plane)..., 1)), (1, 2, 4, 3))
    b_plane = permutedims(reshape(b_plane, (size(b_plane)..., 1)), (1, 2, 4, 3))

    # Concatenate along the 3rd dimension
    return cat(r_plane, g_plane, b_plane; dims=3)
end

function pixel_pack(pixels)
    packed_pixels = zeros(UInt32, (size(pixels,1), size(pixels,2), size(pixels,4)))
    for x in 1:size(pixels, 1)
        for y in 1:size(pixels, 2)
            for batch_idx in 1:size(pixels, 4)
                # Pack those pixels in good boy
                packed_pixels[x, y, batch_idx] = (
                    UInt32(pixels[x, y, 1, batch_idx]) << 24 +
                    UInt32(pixels[x, y, 2, batch_idx]) << 16 +
                    UInt32(pixels[x, y, 3, batch_idx]) << 8
                )
            end
        end
    end
    return packed_pixels
end


text:
using Test
using QuantumOptics
using SparseArrays, LinearAlgebra

@testset "superoperators" begin

# Test creation
b = GenericBasis(3)
@test_throws DimensionMismatch DenseSuperOperator((b, b), (b, b), zeros(ComplexF64, 3, 3))
@test_throws DimensionMismatch SparseSuperOperator((b, b), (b, b), spzeros(ComplexF64, 3, 3))

# Test copy, sparse and dense
b1 = GenericBasis(2)
b2 = GenericBasis(7)
b3 = GenericBasis(5)
b4 = GenericBasis(3)

s = DenseSuperOperator((b1, b2), (b3, b4))
s_ = dense(s)
s_.data[1,1] = 1
@test s.data[1,1] == 0
s_sparse = sparse(s_)
@test isa(s_sparse, SparseSuperOperator)
@test s_sparse.data[1,1] == 1

s = SparseSuperOperator((b1, b2), (b3, b4))
s_ = sparse(s)
s_.data[1,1] = 1
@test s.data[1,1] == 0
s_dense = dense(s_)
@test isa(s_dense, DenseSuperOperator)
@test s_dense.data[1,1] == 1

# Test length
b1 = GenericBasis(3)
b2 = GenericBasis(5)
op = DenseOperator(b1, b2)
S = spre(op)
@test length(S) == length(S.data) == (3*5)^2

# Test arithmetic
b1 = GenericBasis(3)
b2 = GenericBasis(5)
b3 = GenericBasis(5)
b4 = GenericBasis(3)

d1 = DenseSuperOperator((b1, b2), (b3, b4))
d2 = DenseSuperOperator((b3, b4), (b1, b2))
s1 = SparseSuperOperator((b1, b2), (b3, b4))
s2 = SparseSuperOperator((b3, b4), (b1, b2))

x = d1*d2
@test isa(x, DenseSuperOperator)
@test x.basis_l == x.basis_r == (b1, b2)

x = s1*s2
@test isa(x, SparseSuperOperator)
@test x.basis_l == x.basis_r == (b1, b2)

x = d1*s2
@test isa(x, DenseSuperOperator)
@test x.basis_l == x.basis_r == (b1, b2)

x = s1*d2
@test isa(x, DenseSuperOperator)
@test x.basis_l == x.basis_r == (b1, b2)

x = d1*3
@test isa(x, DenseSuperOperator)
@test x.basis_l == (b1, b2)
@test x.basis_r == (b3, b4)

x = 2.5*s1
@test isa(x, SparseSuperOperator)
@test x.basis_l == (b1, b2)
@test x.basis_r == (b3, b4)

x = d1 + d1
@test isa(x, DenseSuperOperator)
@test x.basis_l == (b1, b2)
@test x.basis_r == (b3, b4)

x = s1 + s1
@test isa(x, SparseSuperOperator)
@test x.basis_l == (b1, b2)
@test x.basis_r == (b3, b4)

x = d1 + s1
@test isa(x, DenseSuperOperator)
@test x.basis_l == (b1, b2)
@test x.basis_r == (b3, b4)

x = d1 - d1
@test isa(x, DenseSuperOperator)
@test x.basis_l == (b1, b2)
@test x.basis_r == (b3, b4)

x = s1 - s1
@test isa(x, SparseSuperOperator)
@test x.basis_l == (b1, b2)
@test x.basis_r == (b3, b4)

x = d1 - s1
@test isa(x, DenseSuperOperator)
@test x.basis_l == (b1, b2)
@test x.basis_r == (b3, b4)

x = -d1
@test isa(x, DenseSuperOperator)
@test x.basis_l == (b1, b2)
@test x.basis_r == (b3, b4)

x = -s1
@test isa(x, SparseSuperOperator)
@test x.basis_l == (b1, b2)
@test x.basis_r == (b3, b4)

# TODO: Clean-up this part
ωc = 1.2
ωa = 0.9
g = 1.0
γ = 0.5
κ = 1.1

T = Float64[0.,1.]

fockbasis = FockBasis(7)
spinbasis = SpinBasis(1//2)
basis = tensor(spinbasis, fockbasis)

sx = sigmax(spinbasis)
sy = sigmay(spinbasis)
sz = sigmaz(spinbasis)
sp = sigmap(spinbasis)
sm = sigmam(spinbasis)

Ha = embed(basis, 1, 0.5*ωa*sz)
Hc = embed(basis, 2, ωc*number(fockbasis))
Hint = sm ⊗ create(fockbasis) + sp ⊗ destroy(fockbasis)
H = Ha + Hc + Hint

Ja = embed(basis, 1, sqrt(γ)*sm)
Jc = embed(basis, 2, sqrt(κ)*destroy(fockbasis))
J = [Ja, Jc]

Ψ₀ = spinup(spinbasis) ⊗ fockstate(fockbasis, 5)
ρ₀ = dm(Ψ₀)

op1 = DenseOperator(spinbasis, [1.2+0.3im 0.7+1.2im;0.3+0.1im 0.8+3.2im])
op2 = DenseOperator(spinbasis, [0.2+0.1im 0.1+2.3im; 0.8+4.0im 0.3+1.4im])
@test tracedistance(spre(op1)*op2, op1*op2) < 1e-12
@test tracedistance(spost(op1)*op2, op2*op1) < 1e-12

@test spre(sparse(op1))*op2 == op1*op2
@test spost(sparse(op1))*op2 == op2*op1

@test spre(sparse(op1))*sparse(op2) == sparse(op1*op2)
@test spost(sparse(op1))*sparse(op2) == sparse(op2*op1)

L = liouvillian(H, J)
ρ = -1im*(H*ρ₀ - ρ₀*H)
for j=J
    ρ .+= j*ρ₀*dagger(j) - 0.5*dagger(j)*j*ρ₀ - 0.5*ρ₀*dagger(j)*j
end
@test tracedistance(L*ρ₀, ρ) < 1e-10

tout, ρt = timeevolution.master([0.,1.], ρ₀, H, J; reltol=1e-7)
@test tracedistance(exp(dense(L))*ρ₀, ρt[end]) < 1e-6

@test dense(spre(op1)) == spre(op1)

@test L/2.0 == 0.5*L == L*0.5
@test -L == SparseSuperOperator(L.basis_l, L.basis_r, -L.data)

@test_throws AssertionError liouvillian(H, J; rates=zeros(4, 4))

rates = diagm(0 => [1.0, 1.0])
@test liouvillian(H, J; rates=rates) == L

# Test broadcasting
@test L .+ L == L + L
Ldense = dense(L)
@test isa(L .+ Ldense, DenseSuperOperator)
L_ = copy(L)
L .+= L
@test L == 2*L_
L .+= Ldense
@test L == 3*L_
Ldense_ = dense(L_)
Ldense .+= Ldense
@test Ldense == 2*Ldense_
Ldense .+= L
@test isa(Ldense, DenseSuperOperator)
@test isapprox(Ldense.data, 5*Ldense_.data)
@test_throws ErrorException cos.(Ldense)
@test_throws ErrorException cos.(L)

end # testset


text:
Tables.istable(::Type{<:StructVector}) = true
Tables.rowaccess(::Type{<:StructVector}) = true
Tables.columnaccess(::Type{<:StructVector}) = true

Tables.rows(s::StructVector) = s
Tables.columns(s::StructVector) = fieldarrays(s)

Tables.schema(s::StructVector) = Tables.Schema(staticschema(eltype(s)))


text:
using Test
using StatsProcedures

using StatsProcedures: _f, _get, groupargs, _sharedby, ≊, _show_args, _parse!, pool, proceed
import StatsProcedures: required, default, transformed, combinedargs, copyargs

testvoidstep(a::String) = NamedTuple()
const TestVoidStep = StatsStep{:TestVoidStep, typeof(testvoidstep), true}
required(::TestVoidStep) = (:a,)

testregstep(a::String, b::String) = (c=a*b,)
const TestRegStep = StatsStep{:TestRegStep, typeof(testregstep), true}
default(::TestRegStep) = (a="a", b="b")

testlaststep(a::String, c::String) = (result=a*c,)
const TestLastStep = StatsStep{:TestLastStep, typeof(testlaststep), true}
default(::TestLastStep) = (a="a",)
transformed(::TestLastStep, ntargs::NamedTuple) = (ntargs.c,)

testcombinestep(a::String, bs::String...) = (c=collect(bs),)
const TestCombineStep = StatsStep{:TestCombineStep, typeof(testcombinestep), true}
default(::TestCombineStep) = (a="a",)
combinedargs(::TestCombineStep, ntargs) = [nt.b for nt in ntargs]

testarraystep(a::String, c::Array) = (result=c,)
const TestArrayStep = StatsStep{:TestArrayStep, typeof(testarraystep), true}
required(::TestArrayStep) = (:a, :c)
copyargs(::TestArrayStep) = (2,)

const TestUnnamedStep = StatsStep{:TestUnnamedStep, typeof(testregstep), true}

@testset "StatsStep" begin
    @testset "_get" begin
        @test _get(NamedTuple(), ()) == ()
        @test _get((a=1, b=2), (:a,)) == (1,)
        @test_throws ErrorException _get((b=2,), (:a,))

        @test _get(NamedTuple(), NamedTuple()) == ()
        @test _get((a=1,), (b=2,)) == (2,)
        @test _get((a=1,), (a=2, b=2)) == (1, 2)
        @test _get((a=1, b=2), (a=2,)) == (1,)
    end

    @testset "args" begin
        @test groupargs(TestVoidStep(), (a="a",)) == ("a",)
        @test_throws ErrorException groupargs(TestVoidStep(), (b="b",))

        @test groupargs(TestRegStep(), NamedTuple()) == ("a", "b")
        @test groupargs(TestRegStep(), (a="a1",)) == ("a1", "b")
        @test groupargs(TestRegStep(), (c="c",)) == ("a", "b")

        @test groupargs(TestUnnamedStep(), (a="a", b="b")) == ()
        @test combinedargs(TestRegStep(), (a="a",)) == ()
    end

    @testset "teststeps" begin
        @test _f(TestVoidStep()) == testvoidstep
        @test TestVoidStep()((a="a", b="b")) == (a="a", b="b")

        @test TestRegStep()() == (c="ab",)
        @test TestRegStep()((a="c", b="d")) == (a="c", b="d", c="cd")

        @test TestLastStep()((a="a", b="a", c="ab")) ==
            (a="a", b="a", c="ab", result="aab")

        @test TestCombineStep()((a="a", b="b")) == (a="a", b="b", c=["b"])
        c = ["c"]
        ret = TestArrayStep()((a="a", c=c,))
        @test ret.result === c

        @test sprint(show, TestVoidStep()) == "TestVoidStep"
        @test sprint(show, MIME("text/plain"), TestVoidStep()) ==
            "TestVoidStep (StatsStep that calls testvoidstep)"
    end
end

struct TestProcedure{A,T} <: AbstractStatsProcedure{A,T} end
const RP = TestProcedure{:RegProcedure,Tuple{TestVoidStep,TestRegStep,TestLastStep}}
const rp = RP()
const IP = TestProcedure{:InverseProcedure,Tuple{TestRegStep,TestVoidStep,TestLastStep}}
const ip = IP()
const UP = TestProcedure{:UnitProcedure,Tuple{TestRegStep}}
const up = UP()
const NP = TestProcedure{:NullProcedure,Tuple{}}
const np = NP()
const CP = TestProcedure{:CombineProcedure,Tuple{TestCombineStep,TestArrayStep}}
const cp = CP()
const AP = TestProcedure{:ArrayProcedure,Tuple{TestArrayStep}}
const ap = AP()

@testset "AbstractStatsProcedure" begin
    @test length(rp) == 3
    @test eltype(AbstractStatsProcedure) == StatsStep
    @test eltype(RP) == StatsStep

    @test firstindex(rp) == 1
    @test lastindex(rp) == 3
    @test rp[end] == TestLastStep()
    @test rp[1] == TestVoidStep()
    @test rp[2:-1:1] == [TestRegStep(), TestVoidStep()]
    @test rp[[2,3]] == [TestRegStep(), TestLastStep()]
    @test_throws BoundsError rp[4]

    @test iterate(rp) == (TestVoidStep(), 2)
    @test iterate(rp, 2) == (TestRegStep(), 3)
    @test iterate(rp, 4) === nothing

    @test length(np) == 0
    @test eltype(NP) == StatsStep
    @test_throws BoundsError np[1]
    @test iterate(np) === nothing

    @test sprint(show, rp) == "RegProcedure"
    @test sprint(show, MIME("text/plain"), rp) == """
        RegProcedure (TestProcedure with 3 steps):
          TestVoidStep |> TestRegStep |> TestLastStep"""

    @test sprint(show, up) == "UnitProcedure"
    @test sprint(show, MIME("text/plain"), up) == """
        UnitProcedure (TestProcedure with 1 step):
          TestRegStep"""

    @test sprint(show, np) == "NullProcedure"
    @test sprint(show, MIME("text/plain"), np) == """
        NullProcedure (TestProcedure with 0 step)"""
end

@testset "SharedStatsStep" begin
    s1 = SharedStatsStep(TestRegStep(), 1)
    s2 = SharedStatsStep(TestRegStep(), [3,2])
    @test _sharedby(s1) == [1]
    @test _sharedby(s2) == [2,3]
    @test _f(s1) == testregstep
    @test groupargs(s1, NamedTuple()) == ("a", "b")

    @test sprint(show, s1) == "TestRegStep"
    @test sprint(show, MIME("text/plain"), s1) ==
        "TestRegStep (StatsStep shared by 1 procedure)"
    @test sprint(show, s2) == "TestRegStep"
    @test sprint(show, MIME("text/plain"), s2) ==
        "TestRegStep (StatsStep shared by 2 procedures)"
end

@testset "PooledStatsProcedure" begin
    ps = AbstractStatsProcedure[rp]
    shared = [SharedStatsStep(s, 1) for s in rp]
    p1 = pool(rp)
    @test p1 == PooledStatsProcedure(ps, shared)
    @test length(p1) == 3
    @test eltype(PooledStatsProcedure) == SharedStatsStep
    @test firstindex(p1) == 1
    @test lastindex(p1) == 3
    @test p1[1] == SharedStatsStep(rp[1], 1)
    @test p1[1:3] == [SharedStatsStep(s, 1) for s in rp]
    @test iterate(p1) == (shared[1], 2)
    @test iterate(p1, 2) == (shared[2], 3)

    ps = AbstractStatsProcedure[rp, rp]
    shared = [SharedStatsStep(s, (1,2)) for s in rp]
    p2 = pool(rp, rp)
    @test p2 == PooledStatsProcedure(ps, shared)
    @test length(p2) == 3

    ps = AbstractStatsProcedure[up, up]
    shared = [SharedStatsStep(s, (1,2)) for s in up]
    p3 = pool(up, up)
    @test p3 == PooledStatsProcedure(ps, shared)
    @test length(p3) == 1

    ps = AbstractStatsProcedure[np]
    shared = []
    p4 = pool(np)
    @test p4 == PooledStatsProcedure(ps, shared)
    @test length(p4) == 0

    ps = AbstractStatsProcedure[up, rp]
    shared = [SharedStatsStep(rp[1], 2), SharedStatsStep(rp[2], (1,2)), SharedStatsStep(rp[3], 2)]
    p5 = pool(up, rp)
    @test p5 == PooledStatsProcedure(ps, shared)
    @test length(p5) == 3

    ps = AbstractStatsProcedure[rp, ip]
    shared = [SharedStatsStep(rp[1], 1), SharedStatsStep(ip[1], 2), SharedStatsStep(rp[2], 1),
        SharedStatsStep(ip[2], 2), SharedStatsStep(rp[3], (1,2))]
    p6 = pool(rp, ip)
    @test p6 == PooledStatsProcedure(ps, shared)
    @test length(p6) == 5

    ps = AbstractStatsProcedure[rp, cp]
    shared = [SharedStatsStep(rp[1], 1), SharedStatsStep(rp[2], 1),
        SharedStatsStep(rp[3], 1), SharedStatsStep(cp[1], 2), SharedStatsStep(cp[2], 2)]
    p7 = pool(rp, cp)
    @test p7 == PooledStatsProcedure(ps, shared)
    @test length(p7) == 5

    @test sprint(show, p1) == "PooledStatsProcedure"
    @test sprint(show, MIME("text/plain"), p1) == """
        PooledStatsProcedure with 3 steps from 1 procedure:
          RegProcedure"""

    @test sprint(show, p2) == "PooledStatsProcedure"
    @test sprint(show, MIME("text/plain"), p2) == """
        PooledStatsProcedure with 3 steps from 2 procedures:
          RegProcedure
          RegProcedure"""

    @test sprint(show, MIME("text/plain"), p3) == """
        PooledStatsProcedure with 1 step from 2 procedures:
          UnitProcedure
          UnitProcedure"""

    @test sprint(show, MIME("text/plain"), p4) == """
        PooledStatsProcedure with 0 step from 1 procedure:
          NullProcedure"""
end

@testset "StatsSpec" begin
    s1 = StatsSpec("name", RP, (a="a", b="b"))
    s2 = StatsSpec("", RP, (a="a", b="b"))
    s3 = StatsSpec("", UP, (a="a", b="b"))
    s4 = StatsSpec("name", RP, (b="b", a="a"))
    s5 = StatsSpec("name", RP, (b="b", a="a", d="d"))
    @test s1 == s2
    @test s2 != s3
    @test s2 != s4
    @test s2 ≊ s4

    @test s1() == "aab"
    @test s3() == "ab"

    @test s1(keep=:a) == (a="a", result="aab")
    @test s1(keep=(:a,:c)) == (a="a", c="ab", result="aab")
    @test_throws ArgumentError s1(keep=1)
    @test s1(keepall=true) == (a="a", b="b", c="ab", result="aab")

    s6 = StatsSpec("", NP, NamedTuple())
    @test s6() === nothing
    @test s6(keepall=true) == NamedTuple()
    @test s6(keep=:result) == NamedTuple()

    @test sprint(show, s1) == "name"
    @test sprint(show, s2) == "unnamed"
    @test sprint(show, MIME("text/plain"), s1) == "name (StatsSpec for RegProcedure)"
    @test sprint(show, MIME("text/plain"), s2) == "unnamed (StatsSpec for RegProcedure)"

    @test _show_args(stdout, s1) === nothing
end

function testparser(args, kwargs)
    kwargs = Dict{Symbol,Any}(kwargs...)
    for arg in args
        if arg isa Type{<:AbstractStatsProcedure}
            kwargs[:p] = arg
        end
    end
    return kwargs
end

testformatter(arg::Dict{Symbol,Any}) = (get(arg, :name, ""), arg[:p], (a=arg[:a], b=arg[:b]))

@testset "proceed" begin
    s1 = StatsSpec("s1", RP, (a="a", b="b"))
    s2 = StatsSpec("s2", RP, (a="a",))
    s3 = StatsSpec("s3", RP, (a="a", b="b1"))
    s4 = StatsSpec("s4", UP, (a="a", b="b"))
    s5 = StatsSpec("s5", IP, (a="a", b="b"))
    s6 = StatsSpec("s6", CP, (a="a", b="b1"))
    s7 = StatsSpec("s7", CP, (a="a", b="b2"))
    c = ["c"]
    s8 = StatsSpec("s8", AP, (a="a", c=c))
    s9 = StatsSpec("s9", AP, (a="a1", c=c))
    s10 = StatsSpec("s10", NP, NamedTuple())

    @test proceed([s1]) == ["aab"]
    @test proceed([s1,s2], verbose=true) == ["aab", "aab"]
    @test proceed([s1,s3], verbose=true) == ["aab", "aab1"]
    @test proceed([s1,s4], verbose=true) == ["aab", "ab"]
    @test proceed([s1,s5], verbose=true) == ["aab", "aab"]
    @test proceed([s1,s4,s5], verbose=true) == ["aab", "ab", "aab"]

    @test proceed([s1], keep=:a) == [(a="a",result="aab")]
    @test proceed([s1], keep=[:a,:b]) == [(a="a", b="b", result="aab")]
    @test proceed([s1], keep=(:d,)) == [(result="aab",)]
    @test proceed([s1], keepall=true) == [(a="a", b="b", c="ab", result="aab")]

    @test proceed([s2]) == ["aab"]
    @test proceed([s2], keep=:b) == [(result="aab",)]
    @test proceed([s2], keepall=true) == [(a="a", c="ab", result="aab",)]

    @test proceed([s1,s4], keep=[:a, :result]) ==
            [(a="a", result="aab"), (a="a",)]
    @test proceed([s1,s4], keepall=true) ==
        [(a="a", b="b", c="ab", result="aab"), (a="a", b="b", c="ab")]

    @test proceed([s6], keepall=true) == [(a="a", b="b1", c=["b1"], result=["b1"])]
    ret = proceed([s6,s7], keepall=true)
    @test ret == [(a="a", b="b1", c=["b1", "b2"], result=["b1", "b2"]),
        (a="a", b="b2", c=["b1", "b2"], result=["b1", "b2"])]
    @test ret[1].c === ret[2].c
    @test ret[1].result === ret[2].result

    ret = proceed([s8], keepall=true)
    @test ret[1].c === ret[1].result
    ret = proceed([s8,s9], keepall=true)
    @test ret[1].c !== ret[1].result

    @test proceed([s10]) == [nothing]
    @test proceed([s10], keepall=true) == NamedTuple[NamedTuple()]
    @test proceed([s10], keep=:result) == NamedTuple[NamedTuple()]

    @test proceed([s1], pause=1) == ["b"]

    @test_throws ArgumentError proceed(StatsSpec[])
end

@testset "_parse!" begin
    options = :(Dict{Symbol, Any}())
    _parse!(options, [:a, :(b=1)])
    @test eval(options) == Dict{Symbol, Any}(:a => true, :b => 1)
    @test _parse!(options, [:noproceed, :(b=1)]) == true
    @test _parse!(options, [:(noproceed=false), :(b=1)]) == false
    @test_throws ArgumentError _parse!(options, [1])
end

@testset "@specset" begin
    s = @specset [noproceed=true] a="a0" begin
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:a=>"a1", :b=>"b")))...)(;) end
    @test s == StatsSpec[StatsSpec("", RP, (a="a1", b="b"))]
    @test proceed(s) == ["a1a1b"]

    s = @specset [noproceed] a="a0" begin
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:b=>"b")))...)(;) end
    @test s == StatsSpec[StatsSpec("", RP, (a="a0", b="b"))]
    @test proceed(s) == ["a0a0b"]

    s = @specset [noproceed] a="a0" b="b0" begin
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}()))...)
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:a=>"a1", :b=>"b1")))...)(;)
    end
    @test s == StatsSpec[StatsSpec("", RP, (a="a0", b="b0")), StatsSpec("", RP, (a="a1", b="b1"))]
    @test proceed(s) == ["a0a0b0", "a1a1b1"]

    s = @specset [noproceed] a="a0" b="b0" begin
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}()))...)(;)
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:a=>"a1", :c=>"c")))...)
    end
    @test s == StatsSpec[StatsSpec("", RP, (a="a0", b="b0")), StatsSpec("", RP, (a="a1", b="b0"))]
    @test proceed(s) == ["a0a0b0", "a1a1b0"]

    a = "a0"
    r = @specset [verbose keepall] a=a begin
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:b=>"b")))...) end
    @test r[2] == [(a="a0", b="b", c="a0b", result="a0a0b")]

    r = @specset [verbose keep=:a] a=a begin
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:b=>"b")))...) end
    @test r[2] == [(a="a0", result="a0a0b")]

    r = @specset [verbose keep=[:a]] a=a begin
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:b=>"b")))...) end
    @test r[2] == [(a="a0", result="a0a0b")]

    r = @specset [verbose pause=1] a=a begin
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:b=>"b")))...) end
    @test r[2] == ["b"]

    s0 = @specset [noproceed] for i in 1:3
        a = "a"*string(i)
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:a=>a, :b=>"b")))...)(;)
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:a=>a, :b=>"b1")))...)(;)
    end
    @test proceed(s0) == ["a1a1b", "a1a1b1", "a2a2b", "a2a2b1", "a3a3b", "a3a3b1"]

    s1 = @specset [noproceed] begin
        i = 1
        a = "a"*string(i)
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:a=>a, :b=>"b")))...)(;)
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:a=>a, :b=>"b1")))...)(;)
        for i in 2:3
            a = "a"*string(i)
            StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:a=>a, :b=>"b")))...)(;)
            StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:a=>a, :b=>"b1")))...)(;)
        end
    end
    @test s1 == s0

    r = @specset for i in 1:3
        a = "a"*string(i)
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:a=>a, :b=>"b")))...)(;)
        StatsSpec(testformatter(testparser(Any[RP], Dict{Symbol,Any}(:a=>a, :b=>"b1")))...)(;)
    end
    @test r[2] == ["a1a1b", "a1a1b1", "a2a2b", "a2a2b1", "a3a3b", "a3a3b1"]

    r = @specset RP a="a1" begin
        StatsSpec(testformatter(testparser(Any[], Dict{Symbol,Any}(:b=>"b")))...)(;)
        for i in 2:3
            StatsSpec(testformatter(testparser(Any[], Dict{Symbol,Any}(:a=>"a"*"$i", :b=>"b")))...)(;)
        end
    end
    @test r[2] == ["a1a1b", "a2a2b", "a3a3b"]
end


text:
using Documenter, CompoundPeriods, Dates

makedocs(
    modules = [CompoundPeriods],
    sitename = "CompoundPeriods",
    pages  = Any[
        "Overview"                           => "index.md",
        "Canonical Compound Periods"         => "canonicalcompoundperiods.md",
        "Reverse Compound Periods"           => "reversecompoundperiods.md",
        "Well Behaved Temporal Composites"   => "wellbehavedtemporalcomposites.md",
        "Compound Period Comparisons"        => "compoundcompare.md",
        "Nanosecond Increments"              => "nanosecondincrements.md",
        "min, max, minmax"                   => "minmaxminmax.md",
        "Dates, Times"                       => "datestimes.md"
        ]
    )

deploydocs(
    repo = "github.com/JeffreySarnoff/CompoundPeriods.jl.git",
    target = "build"
)


text:
"""
# Thermodynamic properties of moist air

This  package implements functions to calculate the
thermodynamic properties of moist air. It uses
correlations recommended by ASHRAE and described
in the articles listed below.

# References

 * [1] Wexler, A. and Hyland, R. W., "Formulations for the thermodynamic properties of the saturated phases of H2O from 173.15 K to 473.15 K", ASHRAE Transactions, 1983.
 * [2] Wexler, A. and Hyland, R. W., "Formulations for the thermodynamic properties of dry air from 173.15 K to 473.15 K, and of saturated moist air from 173.15 K to 372.15 K at pressures to 5 MPa
 * [3] Himmelblaum D. M., "Solubilities of inert gases in water, 0oC to near the critical point of water", Journal of Chemical and Engineering Data, Vol. 5, No. 1, January 1960.
 * [4] Kell, George S., "Density, thermal expansivity, and compressibility of liquid water from 0oC to 150oC: correlations and tables for atmospheric pressure and saturation reviewed and expressed on 1968 temperature scale", Journal of Chemical and Engineering Data, Vol. 20, No. 1, 1975.
 * [5] ASHRAE, "Psychrometrics: Theory and Practice", ASHRAE, 1996.

# User interface - Thermodynamic properties of moist air, dry air and saturated water vapor.

The methods listed below calculate thermodynamic properties of moist air:

    volume(MoistAir, T, HumidityType, humidity, P[, outunit])
    volume(MoistAir, T, HumidityType, humidity, P[, outunit])
    density(MoistAir, T, HumidityType, humidity, P[, outunit])
    enthalpy(MoistAir, T, HumidityType, humidity, P[, outunit])
    enthalpym(MoistAir, T, HumidityType, humidity, P[, outunit])
    entropy(MoistAir, T, HumidityType, humidity, P[, outunit])
    entropym(MoistAir, T, HumidityType, humidity, P[, outunit])
    compressfactor(MoistAir, T, HumidityType, humidity, P[, outunit])
    dewpoint(MoistAir, T, HumidityType, humidity, P[, outunit])
    wetbulb(MoistAir, T, HumidityType, humidity, P[, outunit])
    humrat(MoistAir, T, HumidityType, humidity, P)
    relhum(MoistAir, T, HumidityType, humidity, P)
    humrat(MoistAir, T, HumidityType, humidity, P)
    spechum(MoistAir, T, HumidityType, humidity, P)
    molarfrac(MoistAir, T, HumidityType, humidity, P)

The methods listed above calculate the following thermodynamic properties of moist air:

 * `volume` Specific volume
 * `volumem` Molar volume
 * `density` Density
 * `enthalpy` Specific enthalpy
 * `enthalpym` Molar enthalpy
 * `entropy` Specific entropy
 * `entropym` Molar entropy
 * `compressfactor` Compressibility factor Z
 * `dewpoint` Dew point temperature
 * `wetbulb` Adiabatic saturation temperature
 * `humrat` Humidity ratio
 * `relhum` Relative humidity
 * `spechum` Specific humidity
 * `molarfrac` Molar fraction of water vapor

The humidity is specified using two parameters:

 * How the humidity is specified
 * The actual value of humidity

The following types are used to characterize the humidity.

 * `WetBulb` for wet bulb temperature, actually adiabatic saturation temperature
 * `DewPoint` Dew point temperature
 * `RelHum` Relative humidity
 * `HumRat` Humidity ratio (kg of vapor / kg of dry air)
 * `SpecHum` Specific humidity (kg of vapor / kg of moist air)
 * `MolarFrac` Molar fraction of water vapor.

# Examples
```julia-repl
julia> volume(MoistAir, 293.15, WetBulb, 291.15, 101325.0)
0.8464079202783964

julia> volume(MoistAir, 293.15, DewPoint, 291.15, 101325.0)
0.8475219875187474

julia> volume(MoistAir, 293.15, RelHum, 0.7, 101325.0)
0.843889817602806

julia> volume(MoistAir, 20.0u"°C", DewPoint, 60.0u"°F", 1.0u"atm")
0.8449934929585231 kg^-1 m^3

julia> volumem(MoistAir, 293.15, RelHum, 0.5, 93000.0)
0.026199080086890276

julia> volumem(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa", u"inch^3/kmol")
1.598733210336603e6 in^3 kmol^-1

julia> density(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa")
1.0976075893895811 kg m^-3

julia> density(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa", u"lb/inch^3")
3.965358988338535e-5 in^-3 lb

julia> volumem(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa", u"inch^3/kmol")
1.598733210336603e6 in^3 kmol^-1

julia> enthalpy(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa")
50667.43014746832 kg^-1 J

julia> enthalpym(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa")
1439.6551689935861 J mol^-1

julia> compressfactor(MoistAir, -90.0u"°C", RelHum, 0.01, 4.5u"MPa")
0.8552758629097985

julia> wetbulb(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa", u"°C")
17.0 °C

julia> dewpoint(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa", u"°C")
15.475836053510477 °C

julia> humrat(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa")
0.012032930694441925

julia> relhum(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa")
0.7517801524436909

julia> spechum(MoistAir, 20.0u"°C", WetBulb, 17.0u"°C", 93u"kPa")
0.011889860823189923

```

"""
module Psychro

"Molecular weight of air in kg/mol"
const Ma = 0.0289635

"Molecular weight of water in kg/mol"
const Mv = 0.01801528

"Universal gas constant kg m^2 /s^2 /K /mol"
const R = 8.314459848

"Melting point of water"
const T0 = 273.15

"Abstract type that represents any kind of matter"
abstract type AbstractMaterial end

"Abstract type that represents any kind of fluid"
abstract type AbstractFluid <: AbstractMaterial end

"""
Type used to model dry air

This uses the formulations presente in the ASHRAE handbook Psychrometrics: Theory and Practice, reference [5].

The methods dealing with this type implement the equations of reference [2].

Some thermodynamic properties are implemented as methods. Every method takes as arguments, the type Vapor` as first argument and the temperature as the second. If `Unitful` package is used, any unit can be used. On the other hand if no units are provided, the packages assumes SI units all through. In particular:

 * Temperature in K
 * Pressure in Pa
 * Length in m
 * Quantity of matter in mol
 * Mass in kg
 * Energy in J

The following methods are implemented:

 * `volume(DryAir, T, P)` specific volume
 * `volumem(DryAir, T, P` molar volume
 * `density(DryAir, T, P)` density
 * `enthalpy(DryAir, T, P)` specific enthalpy
 * `enthalpym(DryAir, T, P)` molar enthalpy
 * `entropy(DryAir, T, P)` specific entropy
 * `entropym(DryAir, T, P)` molar entropy
 * `compressfactor(DryAir, T, P)` Compressibility factor.

"""
struct DryAir <: AbstractFluid end

"""
Type used to model saturated water vapor

This uses the formulations presente in the ASHRAE handbook Psychrometrics: Theory and Practice, reference [5].

The methods dealing with this type implement the equations of reference [1].

Some thermodynamic properties are implemented as methods. Every method takes as arguments, the type `DryAir` as first argument and the saturation temperature. If `Unitful` package is used, any unit can be used. On the other hand if no units are provided, the packages assumes SI units all through. In particular:

 * Temperature in K
 * Pressure in Pa
 * Length in m
 * Quantity of matter in mol
 * Mass in kg
 * Energy in J

The following methods are implemented:

 * `volume(Vapor, T)` specific volume
 * `volumem(Vapor, T)` molar volume
 * `density(Vapor, T)` density
 * `enthalpy(Vapor, T)` specific enthalpy
 * `enthalpym(Vapor, T)` molar enthalpy
 * `entropy(Vapor, T)` specific entropy
 * `entropym(Vapor, T)` molar entropy
 * `compressfactor(Vapor, T)` Compressibility factor.

"""
struct Vapor <: AbstractFluid end

"""
Type used to model moist air

This uses the formulations presente in the ASHRAE handbook Psychrometrics: Theory and Practice, reference [5].

The temperature range is `173.15 < T < 474.15 K` with pressures `P < 5 MPa`.

The methods dealing with this type implement the equations of references [1-2] to calculate some thermodynamics of real moits air.

Some thermodynamic properties are implemented as methods. Every method takes as arguments, the type `MoistAir` as first argument, temperature as second argument, a type declaring how humidity is specified, the value of humidity using and pressure. If `Unitful` package is used, any unit can be used. On the other hand if no units are provided, the packages assumes SI units all through. In particular:

 * Temperature in K
 * Pressure in Pa
 * Length in m
 * Quantity of matter in mol
 * Mass in kg
 * Energy in J

To specify humidity the following types are implemented:

 * `WetBulb` for wet bulb temperature, actually adiabatic saturation temperature
 * `DewPoint` Dew point temperature
 * `RelHum` Relative humidity
 * `HumRat` Humidity ratio (kg of vapor / kg of dry air)
 * `SpecHum` Specific humidity (kg of vapor / kg of moist air)
 * `MolarFrac` Molar fraction of water vapor.

The following methods are implemented:

 * `volume(MoistAir, T, HumType, h, P)` specific volume
 * `volumem(MoistAir, T, HumType, h, P` molar volume
 * `density(MoistAir, T, HumType, h, P)` density
 * `enthalpy(MoistAir, T, HumType, h, P)` specific enthalpy
 * `enthalpym(MoistAir, T, HumType, h, P)` molar enthalpy
 * `entropy(MoistAir, T, HumType, h, P)` specific entropy
 * `entropym(MoistAir, T, HumType, h, P)` molar entropy
 * `compressfactor(MoistAir, T, HumType, h, P)` Compressibility factor.

It should be noted that the specific enthalpy, entropy and volume are calculated with respect to the mass of dry air. Therefore, 1 J/kg is actually 1 J per kg of dry air.
"""
struct MoistAir <: AbstractFluid end

abstract type ThermodynamicProperty end
abstract type PsychroProperty end

struct Volume <: ThermodynamicProperty end
struct Enthalpy <: ThermodynamicProperty end
struct Entropy <: ThermodynamicProperty end
struct Temperature <: ThermodynamicProperty end
struct Pressure <: ThermodynamicProperty end
struct Density <: ThermodynamicProperty end

struct MolarFrac <: PsychroProperty end

"Specific Humidity mv / (ma+mv)"
struct SpecHum <: PsychroProperty end
struct WetBulb <: PsychroProperty end
struct DewPoint <: PsychroProperty end
struct HumRat <: PsychroProperty end
struct RelHum <: PsychroProperty end

# package code goes here
include("utilities.jl")
include("hyland83.jl")
include("hyland83a.jl")
include("wetbulb.jl")
include("psychrounits.jl")
include("moistair.jl")

export volume, volumem, density, enthalpy, enthalpym, entropy, entropym, compressfactor
export dewpoint, relhum, humrat, wetbulb, molarfrac, spechum
export DryAir, Vapor, MoistAir
export Volume, Enthalpy, Entropy, Temperature, Pressure, Density
export MolarFrac, SpecHum, WetBulb, DewPoint, HumRat, RelHum

end # module


text:
# # Data Transformations
#
# ## Histogram
#
# ```@docs
# histogram
# ```

using AlgebraOfGraphics, CairoMakie
set_aog_theme!()

df = (x=randn(1000), y=randn(1000), z=rand(["a", "b", "c"], 1000))
specs = data(df) * mapping(:x, layout=:z) * histogram(bins=range(-2, 2, length=15))
draw(specs)

#

specs = data(df) * mapping(:x, dodge=:z, color=:z) * histogram(bins=range(-2, 2, length=15))
draw(specs)

#

specs = data(df) * mapping(:x, stack=:z, color=:z) * histogram(bins=range(-2, 2, length=15))
draw(specs)

#

data(df) * mapping(:x, :y, layout=:z) * histogram(bins=15) |> draw

# ## Density
#
# ```@docs
# AlgebraOfGraphics.density
# ```

df = (x=randn(5000), y=randn(5000), z=rand(["a", "b", "c", "d"], 5000))
data(df) * mapping(:x, layout=:z) * AlgebraOfGraphics.density() |> draw

#

data(df) * mapping(:x, :y, layout=:z) * AlgebraOfGraphics.density(npoints=50) |> draw

#

specs = data(df) * mapping(:x, :y, layout=:z) *
    AlgebraOfGraphics.density(npoints=50) * visual(Surface)

draw(specs, axis=(type=Axis3, zticks=0:0.1:0.2, limits=(nothing, nothing, (0, 0.2))))

# ## Frequency
#
# ```@docs
# frequency
# ```

df = (x=rand(["a", "b", "c"], 100), y=rand(["a", "b", "c"], 100), z=rand(["a", "b", "c"], 100))
specs = data(df) * mapping(:x, layout=:z) * frequency()
draw(specs)

#

specs = data(df) * mapping(:x, layout=:z, color=:y, stack=:y) * frequency()
draw(specs)

#

specs = data(df) * mapping(:x, :y, layout=:z) * frequency()
draw(specs)

# ## Expectation
#
# ```@docs
# expectation
# ```

df = (x=rand(["a", "b", "c"], 100), y=rand(["a", "b", "c"], 100), z=rand(100), c=rand(["a", "b", "c"], 100))
specs = data(df) * mapping(:x, :z, layout=:c) * expectation()
draw(specs)

#

specs = data(df) * mapping(:x, :z, layout=:c, color=:y, dodge=:y) * expectation()
draw(specs)

#

specs = data(df) * mapping(:x, :y, :z, layout=:c) * expectation()
draw(specs)

# ## Linear
#
# ```@docs
# linear
# ```

x = 1:0.05:10
a = rand(1:7, length(x))
y = 1.2 .* x .+ a .+ 0.5 .* randn.()
df = (; x, y, a)
specs = data(df) * mapping(:x, :y, color=:a => nonnumeric) * (linear() + visual(Scatter))
draw(specs)

# ## Smoothing
#
# ```@docs
# smooth
# ```

x = 1:0.05:10
a = rand(1:7, length(x))
y = sin.(x) .+ a .+ 0.1 .* randn.()
df = (; x, y, a)
specs = data(df) * mapping(:x, :y, color=:a => nonnumeric) * (smooth() + visual(Scatter))
draw(specs)


text:
module TestPush

include("preamble.jl")
using BangBang.NoBang: push

@testset begin
    @test push!!([0.0], 1.0) == [0.0, 1.0]
    @test push!!([0], 1.0, 2im, UInt(3)) == [0.0, 1.0, 2im, 3]
    @test push!!([0], 1.0) == [0.0, 1.0]
    @test push!!((0,), 1) === (0, 1)
    @test push((0,), 1) === (0, 1)
    @test push!!((0,), 1.0, 2im, UInt(3)) === (0, 1.0, 2im, UInt(3))
    @test push((0,), 1.0, 2im, UInt(3)) === (0, 1.0, 2im, UInt(3))
    @test push!!((a=0,), :b => 1) === (a=0, b=1)
    @test push!!((a=0,), Val(:b) => 1) === (a=0, b=1)
    @test push!!(SVector(0), 1) === SVector(0, 1)
    @test push!!(SVector(0), 0.5) === SVector(0.0, 0.5)
    @test push!!(ImmutableDict(:a=>1), :b=>2) ==
        ImmutableDict(ImmutableDict(:a=>1), :b=>2)
    @test push!!(Dict(), :a=>1)::Dict{Any,Any} == Dict(:a=>1)
    @test push!!(Dict{Union{},Union{}}(), :a=>1)::Dict{Symbol,Int} == Dict(:a=>1)
    @test push!!(Set(), 1)::Set{Any} == Set(1)
    @test push!!(Set{Union{}}(), 1)::Set{Int} == Set(1)
end

@testset "error" begin
    err = @test_error push!!((a=1,), :b)
    @test occursin("`push(::NamedTuple, x::Symbol)` is not supported",
                   sprint(showerror, err))
end

end  # module


text:
@enum MinifyState begin
    MINIFY_INIT
    MINIFY_MAP_KEY
    MINIFY_MAP_KEY_FIRST
    MINIFY_MAP_VAL
    MINIFY_ARRAY
    MINIFY_ARRAY_FIRST
end

const NULL = Vector{UInt8}("null")
const COLON = UInt8(':')
const COMMA = UInt8(',')
const QUOTE = UInt8('"')
const OPEN_BRACE = UInt8('{')
const CLOSE_BRACE = UInt8('}')
const OPEN_BRACKET = UInt8('[')
const CLOSE_BRACKET = UInt8(']')

"""
    Minifier(io::IO=stdout) -> Minifier

Removes all unecessary whitespace from JSON.

## Example
```julia-repl
julia> String(take!(YAJL.run(IOBuffer("{    }"), YAJL.Minifier(IOBuffer()))))
"{}"
```
"""
mutable struct Minifier{T<:IO} <: Context
    io::T
    state::Vector{MinifyState}

    Minifier(io::T=stdout) where T <: IO = new{T}(io, [MINIFY_INIT])
end

collect(ctx::Minifier) = ctx.io

writeval(ctx::Minifier, v::Ptr{UInt8}, len::Int) = unsafe_write(ctx.io, v, len)
writeval(ctx::Minifier, v::Vector{UInt8}, ::Int=-1) = write(ctx.io, v)
writeval(ctx::Minifier, v::UInt8, ::Int=-1) = write(ctx.io, v)

# Yuck.
function addval(ctx::Minifier, v, len::Int=-1; string::Bool=false)
    state = ctx.state[end]
    if state === MINIFY_INIT
        string && writeval(ctx, QUOTE)
        writeval(ctx, v, len)
        string && writeval(ctx, QUOTE)
    elseif state === MINIFY_MAP_KEY
        writeval(ctx, [COMMA, QUOTE])
        writeval(ctx, v, len)
        writeval(ctx, [QUOTE, COLON])
        ctx.state[end] = MINIFY_MAP_VAL
    elseif state === MINIFY_MAP_KEY_FIRST
        writeval(ctx, QUOTE)
        writeval(ctx, v, len)
        writeval(ctx, [QUOTE, COLON])
        ctx.state[end] = MINIFY_MAP_VAL
    elseif state === MINIFY_MAP_VAL
        string && writeval(ctx, QUOTE)
        writeval(ctx, v, len)
        string && writeval(ctx, QUOTE)
        ctx.state[end] = MINIFY_MAP_KEY
    elseif state === MINIFY_ARRAY
        writeval(ctx, COMMA)
        string && writeval(ctx, QUOTE)
        writeval(ctx, v, len)
        string && writeval(ctx, QUOTE)
    elseif state === MINIFY_ARRAY_FIRST
        string && writeval(ctx, QUOTE)
        writeval(ctx, v, len)
        string && writeval(ctx, QUOTE)
        ctx.state[end] = MINIFY_ARRAY
    end
end

@yajl null(ctx::Minifier) = addval(ctx, NULL)
@yajl boolean(ctx::Minifier, v::Bool) = addval(ctx, v)
@yajl number(ctx::Minifier, v::Ptr{UInt8}, len::Int) = addval(ctx, v, len)
@yajl string(ctx::Minifier, v::Ptr{UInt8}, len::Int) = addval(ctx, v, len; string=true)
@yajl function map_start(ctx::Minifier)
    ctx.state[end] in (MINIFY_ARRAY, MINIFY_MAP_KEY) && writeval(ctx, COMMA)
    writeval(ctx, OPEN_BRACE)
    push!(ctx.state, MINIFY_MAP_KEY_FIRST)
end
@yajl map_key(ctx::Minifier, v::Ptr{UInt8}, len::Int) = addval(ctx, v, len)
@yajl function map_end(ctx::Minifier)
    writeval(ctx, CLOSE_BRACE)
    pop!(ctx.state)
    state = ctx.state[end]
    if state === MINIFY_MAP_KEY_FIRST
        ctx.state[end] = MINIFY_MAP_KEY
    elseif state === MINIFY_ARRAY_FIRST
        ctx.state[end] = MINIFY_ARRAY
    end
end
@yajl function array_start(ctx::Minifier)
    ctx.state[end] in (MINIFY_ARRAY, MINIFY_MAP_KEY) && writeval(ctx, COMMA)
    writeval(ctx, OPEN_BRACKET)
    push!(ctx.state, MINIFY_ARRAY_FIRST)
end
@yajl function array_end(ctx::Minifier)
    writeval(ctx, CLOSE_BRACKET)
    pop!(ctx.state)
    state = ctx.state[end]
    if state === MINIFY_MAP_KEY_FIRST
        ctx.state[end] = MINIFY_MAP_KEY
    elseif state === MINIFY_ARRAY_FIRST
        ctx.state[end] = MINIFY_ARRAY
    end
end


text:
using Dates
using ..Common

function get_message_digest(doc::PDDoc, brange::Vector{Int}, algo::String)
    ctx = DigestContext(algo)
    for i = 1:2:length(brange)
        buff = readfrom(doc.cosDoc, brange[i], brange[i+1])
        update!(ctx, buff)
    end
    return close(ctx)
end

function verify_local_trust(store::CertStore, cert::Cert)
    ret = find_cert(store, cert)
    !ret && is_self_signed(cert) &&
        error("Self-signed cert not found in certificate trust store")
    return ret
end

# Verification of certificate validity
function verify!(status::Dict{Symbol, Any}, cert::Cert, store::CertStore)
    crt_info = get_info(cert)
    vrf_chain = get!(status, :certs, Dict{Symbol, Any}[])
    push!(vrf_chain, crt_info)
    if verify_local_trust(store, cert)
        crt_info = get_info(cert)
        vrf_chain = get!(status, :chain, Dict{Symbol, Any}[])
        push!(vrf_chain, crt_info)
        status[:passed] = true
        return true
    end
    # The signing time for validation as PKCS-1 does not have a timestamp.
    # Hence, the PDF atrribute is used.
    st = get(status, :M, nothing)
    if st !== nothing
        st = getUTCTime(st)
        t = datetime2unix(st.d)
        d = (:atepoch => t)
        set_params!(store, d)
    end

    passed, data = verify(store, cert)
    status[:passed] = passed
    if passed
        status[:chain] = data
    else
        status[:error_message] = data
    end
    return status[:passed]
end

function verify!(status::Dict{Symbol, Any}, si::PKCS1SignedInfo, cert::Cert,
                 store::CertStore, doc::PDDoc, brange::Vector{Int})
    verify!(status, cert, store)
    status[:passed] || return
    xtract_digest, digest_algo = get_rsa_digest(si, cert)
    compute_digest = get_message_digest(doc, brange, digest_algo)
    status[:passed] = res = (xtract_digest == compute_digest)
    !res && error("Computed digest is not the same as in the signature")
end

function verify!(status::Dict{Symbol, Any}, si::CMSSignedInfo, ::Nothing,
                 store::CertStore, doc::PDDoc, brange::Vector{Int})
    detached, bio = si.detached, C_NULL
    hdata = Vector{UInt8}(undef, 0)
    if detached
        for i = 1:2:length(brange)
            buff = readfrom(doc.cosDoc, brange[i], brange[i+1])
            append!(hdata, buff)
        end
    else
        hdata = get_message_digest(doc, brange, "SHA1")
    end

    status[:certs], signed_cert = harvest_certs(si)

    signer_info = get_signer_info(si)
    # Deciding signing time the order of values to be picked up
    # 1. Signing time attribute from the signed attributes of signer info
    # 2. Time attribute picked up from the timestamp
    # 3. Time attribute picked up from the PDF dictionary
    st = get_signer_info_signing_time(signer_info)
    st === nothing && (st = get_signer_info_timestamp(signer_info))
    st === nothing && (st = get(status, :M, nothing))
    if st !== nothing
        st = getUTCTime(st)
        t = datetime2unix(st.d)
        d = Dict(:atepoch => t)
        set_params!(store, d)
    end
    flags = Cint(CMS_BINARY)
    verify_local_trust(store, signed_cert) &&
        (flags |= CMS_NO_SIGNER_CERT_VERIFY)
    status[:subfilter] === cn"ETSI.CAdES.detached" && (flags |= CMS_CADES)
    verify(si, store, hdata, flags)
    chain = get!(status, :chain, Dict{Symbol, Any}[])
    push!(chain, si |> get_signer |> get_info)
    return (status[:passed] = true)
end

function pd_validate_signature(doc::PDDocImpl,
                               sig::Tuple{IDD{CosDict}, Dict{Symbol, Any}})
    status = sig[2]
    status[:passed] = false
    sigdict = cosDocGetObject(doc.cosDoc, sig[1], cn"V")
    sigdict === CosNull && return status[:passed]
    subfilter = cosDocGetObject(doc.cosDoc, sigdict, cn"SubFilter")
    !(subfilter in (cn"adbe.x509.rsa_sha1",
                    cn"adbe.pkcs7.detached",
                    cn"adbe.pkcs7.sha1",
                    cn"ETSI.CAdES.detached",
                    cn"ETSI.RFC3161")) && return false
    bobj = cosDocGetObject(doc.cosDoc, sigdict, cn"ByteRange")
    bobj === CosNull && return status[:passed]
    brange = get(bobj, true)
    # This should not go through document security handler.
    # Hence, direct COS method. DO NOT USE cosDocGetObject.
    cobj = get(sigdict, cn"Contents")
    cobj === CosNull && return status[:passed]
    contents = Vector{UInt8}(cobj)
    si, cert = nothing, nothing

    try
        status[:subfilter] = subfilter
        if subfilter === cn"adbe.pkcs7.sha1"
            si = CMSSignedInfo(contents)
        elseif subfilter === cn"adbe.pkcs7.detached" ||
            subfilter === cn"ETSI.CAdES.detached"
            si = CMSSignedInfo(contents, true)
        elseif subfilter === cn"adbe.x509.rsa_sha1"
            si = PKCS1SignedInfo(contents)
            cert_obj = cosDocGetObject(doc.cosDoc, sigdict, cn"Cert")
            cdata = Vector{UInt8}(cert_obj)
            cert = Cert(cdata)
        end

        name_obj = cosDocGetObject(doc.cosDoc, sigdict, cn"Name")
        if name_obj !== CosNull
            status[:Name] = CDTextString(name_obj)
        end
        m_obj    = cosDocGetObject(doc.cosDoc, sigdict, cn"M")
        if m_obj !== CosNull
            status[:M] = CDDate(CDTextString(m_obj))
        end
        store = CertStore()
        verify!(status, si, cert, store, doc, brange)
    catch e
        delete!(status, :chain)
        status[:passed] = false
        status[:error_message] = e.msg
        status[:stacktrace] = string.(stacktrace(catch_backtrace()))
    end
    return status[:passed]
end


text:
module FinEconResearchUtils

using Reexport
    @reexport using Dates
    @reexport using Statistics
    @reexport using DataFrames
    @reexport using DataFramesMeta
    @reexport using Missings
    @reexport using ShiftedArrays
    @reexport using CSVFiles
    @reexport using ExcelFiles
    @reexport using StatFiles
    @reexport using SASLib
    @reexport using BeaData
    @reexport using FredData
    @reexport using StatsBase
    @reexport using GLM
    @reexport using Distributions
    @reexport using RollingFunctions

include("utils.jl")

export rollingfillfirst,
    detrend, hp_filter,
    tomissing!, alltomissing!,
    maxdrawdown,
    lastnonzero, dropcols, yrqtrfun

end # module


text:
abstract type MarkovProcess end

mutable struct JumpProcess <: MarkovProcess
    x # state
    a # propensities
    h # jumps
    X # support
end

mutable struct ReactionProcess <: MarkovProcess
    ReactionSystem::ReactionSystem
    JumpProcess::JumpProcess
    species_to_index::Dict
    species_to_state::Dict
    state_to_species::Dict
    # constructor
    function ReactionProcess(rn::ReactionSystem)
        specs = species(rn)
        S = prodstoichmat(rn) - substoichmat(rn)
        n = length(specs)
        @polyvar(x[1:n])
        spec2idx = Dict(specs[i] => i for i in 1:n)
        spec2state = Dict(specs[i] => x[i] for i in 1:n)
        state2spec = Dict(x[i] => specs[i] for i in 1:n)
        props = reformat_reactions(reactions(rn), spec2idx, x)
        jumps = reformat_jumps(S, spec2idx, x)
        support = intersect([@set(x[i] >= 0) for i in 1:n]...)
        return new(rn,JumpProcess(x,props,jumps,support), spec2idx, spec2state, state2spec)
    end
end

mutable struct DiffusionProcess <: MarkovProcess
    x # state
    f # drift
    σ # diffusion matrix
    X # support
end

mutable struct JumpDiffusionProcess <: MarkovProcess
    JumpProcess::JumpProcess
    DiffusionProcess::DiffusionProcess
end

mutable struct ControlProcess
    MP::MarkovProcess
    T # time horizon
    u # control variables
    t # time variable
    U # control set
    PathChanceConstraints
    TerminalChanceConstraints
    Objective
    discount_factor
end

mutable struct ExitProbability
    X
end

mutable struct TerminalSetProbability
    X
end

mutable struct LagrangeMayer
    l
    m
end

mutable struct ChanceConstraint
    X # semialgebriac set
    α # confidence level
end

mutable struct InfoData
    termination_status
    solution_time
end

mutable struct Bounds
    type
    order::Int
    bounds
    info::InfoData
end

inf_generator(MP::JumpProcess, p::Polynomial) = sum(MP.a[i]*(subs(p, MP.x => MP.h[i]) - p) for i in 1:length(MP.a))
inf_generator(MP::ReactionProcess, p::Polynomial) = inf_generator(MP.JumpProcess,p)
inf_generator(MP::DiffusionProcess, p::Polynomial) = MP.f'*∂(p,MP.x) + 1/2*sum(∂²(p,MP.x,MP.x) .* MP.σ)
inf_generator(MP::JumpDiffusionProcess, p::Polynomial) = inf_generator(MP.JumpProcess,p) + inf_generator(MP.DiffusionProcess,p)
extended_inf_generator(MP::MarkovProcess, p::Polynomial, t::PolyVar; scale = 1) = ∂(p,t) + scale*inf_generator(MP, p)
louiville(MP, P, ξ, ϕ, t, Δt; scale = 1) = Mom(split_poly(ϕ, MP.x, Δt), P) - Mom(extended_inf_generator(MP, ϕ, t; scale=scale), ξ)
function transform_state!(P::JumpProcess, x::AbstractVector, z::AbstractVector; iv::AbstractVector = 1:length(P.x))
    Π = [P.x[i] => z[i] for i in 1:length(P.x)]
    P.a = subs.(P.a, Π...)
    P.h = [subs.(h[iv], Π...) for h in P.h]
    P.X = intersect([@set(subs.(p, Π...) >= 0) for p in P.X.p]...)
    P.x = x
end

function transform_state!(P::ReactionProcess, x::AbstractVector, z::AbstractVector; iv::AbstractVector = 1:length(P.JumpProcess.x))
    Π = [P.JumpProcess.x[i] => z[i] for i in 1:length(P.JumpProcess.x)]
    P.species_to_state = Dict(spec => subs(P.species_to_state[spec], Π...) for spec in keys(P.species_to_state))
    P.state_to_species = Dict(P.species_to_state[spec] => spec for spec in keys(P.species_to_state))
    transform_state!(P.JumpProcess, x, z; iv = iv)
end

function rescale_state!(P::DiffusionProcess, x0)
    transform_state!(P, P.x, P.x .* x0)
    P.f ./= x0
    P.σ ./= x0*x0'
end

function rescale_state!(P::JumpProcess, x0)
    transform_state!(P, P.x, P.x .* x0)
    for i in 1:length(P.h)
        P.h[i] ./= x0
    end
end

function rescale_state!(P::JumpDiffusionProcess, x0)
    rescale_state!(P.JumpProcess.x, x0)
    rescale_state!(P.DiffusionProcess.x, x0)
end

function rescale_state!(P::ReactionProcess, x0::AbstractVector)
    transform_state!(P, P.JumpProcess.x, P.JumpProcess.x .* x0)
    for i in 1:length(P.JumpProcess.h)
        P.JumpProcess.h[i] ./= x0
    end
end

function partition_variables(B)
    m, n = size(B) # m < n if well-posed => m
    @assert(m <= n) # implies R is m × n
    Q, R = qr(B)
    ds = []
    i, k = 1, 1
    while k <= n && i <= m
        if abs(R[i,k]) >= 1e-12
            push!(ds, k)
            i += 1
        end
        k += 1
    end
    is = setdiff(1:n, ds)
    return is, ds, Q, R
end

function project_into_subspace!(P::MarkovProcess, B, f)
    iv, dv, Q, R = partition_variables(B)
    if !isempty(dv)
        z = Vector{Polynomial}(undef, size(B,2))
        @polyvar(x[1:length(iv)])
        z[iv] .= polynomial.(x)
        z[dv] .= polynomial.(round.(R[:,dv]\(Q'*f), digits = 12) - round.((R[:,dv]\R[:,iv]), digits = 12)*x)
        transform_state!(P, x, z; iv=iv)
    end
end

function project_into_subspace!(P::ReactionProcess, x0::AbstractVector)
    B = nullspace(stoichmat(P.ReactionSystem))'
    if !isempty(B)
        f = B*x0
        project_into_subspace!(P, B, f)
    end
end


text:
module LeastSquaresOptim

##############################################################################
##
## Dependencies
##
##############################################################################

import Base: copyto!, fill!, eltype, length, size
import LinearAlgebra: mul!, rmul!, norm, cholesky!, qr!, Symmetric, dot, eigen, axpy!, svd, ldiv!, Transpose, adjoint, Adjoint
import LinearAlgebra.BLAS: gemm!
import Printf: @printf, @sprintf
if Base.USE_GPL_LIBS
    import SuiteSparse.SPQR: QRSparse
    import SparseArrays: SparseMatrixCSC, sparse, nzrange, nonzeros
end
import Statistics: mean
using ForwardDiff, DiffEqDiffTools
import Optim: optimize

##############################################################################
##
## Exported methods and types
##
##############################################################################

export optimize!,
LeastSquaresProblem,
LeastSquaresProblemAllocated,
LeastSquaresResult,
Dogleg,
LevenbergMarquardt,
optimize

##############################################################################
##
## Load files
##
##############################################################################

include("utils/lsmr.jl")
include("utils/utils.jl")

include("types.jl")
include("optimizer/levenberg_marquardt.jl")
include("optimizer/dogleg.jl")

include("solver/dense_qr.jl")
include("solver/dense_cholesky.jl")
include("solver/iterative_lsmr.jl")

end


text:
"""
Use the following command in REPL to make a sysimage:
using PackageCompiler
create_sysimage([:ThreeDeconv, :PyCall], sysimage_path="sys_threedeconv.so", precompile_execution_file="example/example.jl")
"""

using ThreeDeconv
using PyCall, Printf
skimage_io = pyimport("skimage.io")

function imread(filename::String)
    img = skimage_io.imread(filename)
    if ndims(img) == 3
        img = permutedims(img, (2, 3, 1))
    end
    return img
end

function imsave(filename::String, img::Array)
    if ndims(img) == 3
        img = permutedims(img, (3, 1, 2))
    end
    skimage_io.imsave(filename, img)
end

imgraw = imread("dataset/raw/bead/raw_cropped/bead_highsnr_raw.tif")
df = imread("dataset/raw/bead/raw_cropped/df.tif")
ff = imread("dataset/raw/bead/raw_cropped/ff_highsnr.tif")
img = (imgraw .- df) ./ ff

obj_mag = 100
camera_pixel_size = 6.5e3 # [nm]
xystep = camera_pixel_size / obj_mag
zstep = 150 # [nm]
medium_index = 1.515
f_tubelens = 200.e6
NA = 1.4
λ = 540 # [nm]

γ, σ = ThreeDeconv.noise_estimation(img, maxnum_pairs = 200)
@printf "Gain: %.1f, Read noise std.: %.1f \n" γ σ

pad = 10
psf_shape = size(img) .+ pad
println("Simulating PSF with the size of $(psf_shape)")
psf = ThreeDeconv.psf(
    NA = NA,
    objective_mag = obj_mag,
    λ = λ,
    medium_index = medium_index,
    psf_shape = psf_shape,
    camera_pixel_size = camera_pixel_size,
    zstep = zstep,
    f_tubelens = f_tubelens,
    oversampling = 4,
)

options =
    ThreeDeconv.DeconvolutionOptions(max_iters = 150, show_trace = true, check_every = 50)
reg = 0.01
result = ThreeDeconv.deconvolve(
    img,
    psf,
    γ,
    σ,
    reg,
    ThreeDeconv.ADMM(scale_problem = true),
    options = options,
)

imsave("deconvoled_scaled_iter150.tif", result.x)


text:
using Documenter, ControlSystems, Plots, LinearAlgebra, DSP
import GR # Bug with world age in Plots.jl, see https://github.com/JuliaPlots/Plots.jl/issues/1047

include("src/makeplots.jl")

makedocs(modules=[ControlSystems],
    format=Documenter.HTML(),
    sitename="ControlSystems.jl",
    pages=[
        "Home" => "index.md",
        "Examples" => Any[
            "Design" => "examples/example.md",
        ],
        "Guide" => Any[
            "Introduction" => "man/introduction.md",
            "Creating Systems" => "man/creating_systems.md",
        ],
        "Functions" => Any[
            "Constructors" => "lib/constructors.md",
            "Analysis" => "lib/analysis.md",
            "Synthesis" => "lib/synthesis.md",
            "Time and Frequency response" => "lib/timefreqresponse.md",
            "Plotting" => "lib/plotting.md",
        ],
    ]
    )
# If not running travis, generate the plots here, even if we are not deploying
if get(ENV, "TRAVIS", "") == ""
    makePlots()
end

# Only build plots in travis if we are deploying
# And dont install the dependencies unless we are deploying
function myDeps()
    if get(ENV, "TRAVIS", "") != ""
        println("Installing deploy dependencies")
        makePlots()
    end
end

deploydocs(
    repo = "github.com/JuliaControl/ControlSystems.jl.git",
    deps = myDeps
)


text:
##########################
## SIMULATION FUNCTIONS ##
##########################

#GPU related functions
print_gpus() = begin
	println( "$(length(devices())) CUDA capable device(s)." )
	for (i,d) = enumerate(devices())
		u = i == 1 ? "*" : " "
		println( "  ($(i-1)$u) $(name(d))")
	end
end

"""
Implementation in multiple threads. Separating the spins in N_parts.
"""
function run_spin_precession_parallel(obj::Phantom,seq::Sequence, t::Array{Float64,1}, Δt::Array{Float64,1};
	M0::Array{Mag,1}, Nthreads::Int=Hwloc.num_physical_cores(), gpu::Function)

	Nt, NΔt, Ns = length(t), length(Δt), prod(size(obj))
	#Put times as row vector
	t = reshape(t,1,Nt)
	Δt = reshape(Δt,1,NΔt)

	S = zeros(ComplexF64, Nt)

	parts = kfoldperm(Ns, Nthreads, type="ordered")

	@threads for p ∈ parts
		@inbounds aux, M0[p] = run_spin_precession(obj[p],seq,t,Δt; M0=M0[p], gpu)
		S .+= aux
		aux = nothing
	end
    S, M0
end

"""
	run_spin_precession(obj,seq,t)

Simulates an MRI sequence `seq` on the Phantom `obj` for time points `t`.
It calculates S(t) = ∫ ρ(x,t) exp(- t/T2(x,t) ) exp(- 𝒊 ϕ(x,t)) dx.
"""
function run_spin_precession(obj::Phantom, seq::Sequence, t::Array{Float64,2}, Δt::Array{Float64,2};
	M0::Array{Mag,1}, gpu::Function)

	𝒊 = 1im; Random.seed!(1)
	T = sum(Δt) #Total length, used for signal relaxation
	t = [t t[end]+Δt[end]] #For trapezoidal integration <---------------------------
    sz = size(obj)
	Nsz = length(sz)

	#DIFFUSION, disabled while I think how to do it efficiently
	# if any(obj.Dλ1 .+ obj.Dλ2 .!= 0)  #No diff optimization
	# 	#TODO: I need to add diff displacement η story between blocks (like ϕ0) <-- This is already taken care of in M0
	# 	η1 = randn(sz...,Nt) |> gpu
	# 	η2 = randn(sz...,Nt) |> gpu
	# 	Dθ = obj.Dθ |> gpu
	# 	Dλ1 = obj.Dλ1 |> gpu
	# 	Dλ2 = obj.Dλ2 |> gpu
	# 	ηx = sqrt.(2Δt .* Dλ1) .* η1 |> gpu
	# 	ηy = sqrt.(2Δt .* Dλ2) .* η2 |> gpu
	# 	ηxp = cumsum(ηx .* cos.(Dθ) .- ηy.*sin.(Dθ), dims=Nsz+1) |> gpu
	# 	ηyp = cumsum(ηy .* cos.(Dθ) .+ ηx.*sin.(Dθ), dims=Nsz+1) |> gpu
	# 	ηzp = 0
	# else
		ηxp = 0
		ηyp = 0
		ηzp = 0
	# end

	#SCANNER
    Gx, Gy, Gz = get_grads(seq,t)
	Gx = Gx |> gpu
	Gy = Gy |> gpu
	Gz = Gz |> gpu
	#SIMULATION
	Mxy = abs.(M0)	|> gpu
	ϕ0 = angle.(M0)	|> gpu
	x0 = obj.x		|> gpu
	y0 = obj.y		|> gpu
	z0 = obj.z		|> gpu
	tp = t[:,1:end-1] .- t[1]	|> gpu # t' = t - t0
	t = t			|> gpu
	Δt = Δt			|> gpu
    xt = x0 .+ obj.ux(x0,y0,z0,t) .+ ηxp |> gpu
	yt = y0 .+ obj.uy(x0,y0,z0,t) .+ ηyp |> gpu
	zt = z0 .+ obj.uz(x0,y0,z0,t) .+ ηzp |> gpu
	#ACQ OPTIMIZATION
    if is_ADC_on(seq, Array(t))
		ϕ =  ϕ0 .- (2π*γ) .* cumtrapz(Δt, xt.*Gx .+ yt.*Gy .+ zt.*Gz )
	else
		ϕ =  ϕ0 .- (2π*γ) .* trapz(Δt, xt.*Gx .+ yt.*Gy .+ zt.*Gz )
	end
	#Mxy preccesion and relaxation
	Δw = obj.Δw  |> gpu
	T2 = obj.T2  |> gpu
	Mxy =  Mxy .* exp.(𝒊.*(ϕ .- Δw.*tp) .- tp./T2 )
	#ACQUIRED SIGNAL
	S = sum(Mxy, dims=1:Nsz)[:] #<--- TODO: add coil sensitivities
	#Mz relaxation
	T1 = obj.T1		|> gpu
	Mz0 = obj.ρ		|> gpu
	Mz = M0.z		|> gpu
	Mz =  Mz .* exp.(-T./T1) .+ Mz0 .* ( 1 .- exp.(-T./T1) )
	# END
	Mxy = Array(Mxy)
	Mz = Array(Mz)
	M0 = Mag.(Mxy[:,end], Mz) #Saving the last magnetization
    Array(S), M0 #Singal, M0
end

run_spin_excitation_parallel(obj, seq, t::Array{Float64,1}, Δt::Array{Float64,1};
	M0::Array{Mag,1}, Nthreads::Int=Hwloc.num_physical_cores(), gpu::Function) = begin
	Nt, NΔt, Ns = length(t), length(Δt), prod(size(obj))
	#Put times as row vector
	t = reshape(t,1,Nt)
	Δt = reshape(Δt,1,NΔt)

	parts = kfoldperm(Ns, Nthreads, type="ordered")

	@threads for p ∈ parts
		@inbounds M0[p] = run_spin_excitation(obj[p],seq,t,Δt; M0=M0[p], gpu)
	end
    M0
end

run_spin_excitation(obj, seq, t::Array{Float64,2}, Δt::Array{Float64,2};
	M0::Array{Mag,1}, gpu::Function) = begin
	#SCANNER
	B1, Δf_rf  = get_rfs(seq,t)
    Gx, Gy, Gz = get_grads(seq,t)
	B1 = B1 |> gpu
	Gx = Gx |> gpu
	Gy = Gy |> gpu
	Gz = Gz |> gpu
	#SIMULATION
	x0 = obj.x		|> gpu
	y0 = obj.y		|> gpu
	z0 = obj.z		|> gpu
	t = t			|> gpu
	Δt = Δt			|> gpu
    xt = x0 .+ obj.ux(x0,y0,z0,t)		|> gpu
	yt = y0 .+ obj.uy(x0,y0,z0,t)		|> gpu
	zt = z0 .+ obj.uz(x0,y0,z0,t)		|> gpu
	ΔB0 = obj.Δw./(2π*γ) .- Δf_rf./γ	|> gpu # ΔB_0 = (B_0 - ω_rf/γ)
	Bz = (Gx.*xt .+ Gy.*yt .+ Gz.*zt) .+ ΔB0	#<-- TODO: This line is very slow, FIX!?
	B = sqrt.(abs.(B1).^2. .+ abs.(Bz).^2.)
	φ = -2π*γ * (B .* Δt) # angle of rotation
	B[B.==0] .= 1e-17; # removes problems when dividing by φ
	Qt = Q.(Array(φ), Array(B1./B), Array(Bz./B)) #TODO: remove the gpu->array step
	Qf = prod( Qt , dims=2 )[:] # equivalent rotation
	M0 = Qf .* M0 #TODO: This and the relaxation are computed in the CPU for now
	#Relaxation
	T = sum(Δt) #Total length, used for signal relaxation
	for (m, ρ, T1, T2) = zip(M0, obj.ρ, obj.T1, obj.T2)
		m.xy *= @fastmath exp(-T/T2)
		m.z = @fastmath m.z*exp(-T/T1) + ρ*(1-exp(-T/T1))
	end
 	M0
end

"""Divides time steps in N_parts blocks. Decreases RAM usage in long sequences."""
function run_sim_time_iter(obj::Phantom,seq::Sequence, t::Array{Float64,1}, Δt;
								Nblocks::Int=16, Nthreads::Int=Hwloc.num_physical_cores(), gpu::Function, w=nothing)
	Nt, Ns = length(t), prod(size(obj))
	blink_window = w !== nothing
	S = zeros(ComplexF64, Nt) #Only one coil for now, TODO: change for more coils to Nt x Coils
	M0 = Mag(obj,:z)
	breaks = get_breaks_in_RF_key_points(seq,t)
    parts = kfoldperm(Nt,Nblocks;type="ordered",breaks)
	Nblocks = length(parts)
	println("Dividing simulation in Nblocks=$Nblocks")
	println("Starting simulation with Nspins=$Ns and Nt=$Nt")
	#Perturbation of spins' position to reduce spurious echoes (?)
	#Test convert T2* to ΔBz with Lorentzian distribution (?)
	# R2prime = 1 ./obj.T2s .- 1 ./ obj.T2 #1/T2* = 1/T2 + 1/T2' and 1/T2' = γΔB
	#obj_p.Δw .+= something
	#TODO: transform suceptibility χ to Δω, for each time-block with FMM-like technique O(nlogn).
	rfs = 0
	pp = Progress(Nblocks)
	for (block, p) = enumerate(parts)
		if is_RF_on(seq, t[p]) && !is_ADC_on(seq, t[p]) #PATCH: the ADC part should not be necessary, but sometimes 1 sample is identified as RF in an ADC block
			@inbounds M0  = run_spin_excitation_parallel(obj, seq, t[p], Δt[p]; M0, Nthreads, gpu)
			rfs += 1
		else
			@inbounds S[p], M0 = run_spin_precession_parallel(obj, seq, t[p], Δt[p]; M0, Nthreads, gpu)
		end
		#Update progress
		next!(pp, showvalues = [(:simulated_blocks, block), (:rf_blocks,rfs)])
		if blink_window #update Progress
			progress = string(floor(Int, block / Nblocks * 100))
			@js_ w (@var progress=$progress;
					document.getElementById("simul_progress").style.width=progress+"%";
					document.getElementById("simul_progress").innerHTML=progress+"%";
					document.getElementById("simul_progress").setAttribute("aria-valuenow", progress);)
		end
	end
	#Output
	t_interp = get_sample_times(seq)
	S_interp = LinearInterpolation(t.+Δt,S)(t_interp) .* get_sample_phase_compensation(seq)
	(S_interp, M0)
end

function simulate(obj::Phantom, seq::Sequence, sys::Scanner; simParams=Dict{String,Any}(), w=nothing)
	#Simulation params
	enable_gpu = get(simParams, "gpu", true)
	gpu(x) = has_cuda() && enable_gpu ? CuArray(x) : x
	Nthreads = get(simParams, "Nthreads", Hwloc.num_physical_cores())
	Δt    = get(simParams, "Δt", 1e-3)
	Δt_rf = get(simParams, "Δt_rf", 1e-4)
	t, Δt = get_uniform_times(seq, Δt; Δt_rf)
	return_Mag = get(simParams, "return_Mag", false)
	end_sim_at = get(simParams, "end_sim_at", Inf)
	if 0 < end_sim_at < dur(seq)
		idx = t .< end_sim_at
		t  =  t[idx]
		Δt = Δt[idx]
	end
	Nt = length(t)
	Nspins = prod(size(obj)...)
	Nblocks = get(simParams, "Nblocks", ceil(Int, 6506*Nt/1.15e6))
    #Simulate
	println("")
	@info "Running simulation... [GPU = $(enable_gpu && has_cuda()), CPU = $Nthreads thread(s)]."
	if return_Mag
		_, out = @time run_sim_time_iter(obj,seq,t,Δt;Nblocks,Nthreads,gpu)
	else
		S, _ = @time run_sim_time_iter(obj,seq,t,Δt;Nblocks,Nthreads,gpu,w)
		out = S ./ Nspins #Acquired data
	end
	out
end

"""Returns magnetization of spins distributed along `z` after running the Sequence `seq`."""
function simulate_slice_profile(seq; z=range(-2e-2,2e-2,200), simParams=Dict{String,Any}("Δt_rf"=>1e-6))
	simParams["return_Mag"] = true
	sys = Scanner()
	phantom = Phantom(;x=zeros(size(z)),z)
	M = simulate(phantom, seq, sys; simParams)
	M
end


text:
@testset "potential_ET scalars" begin
    Tair,pressure,Rn = 30.0,100.0,500.0
    ET_pot, LE_pot = potential_ET(PriestleyTaylor(), Tair,pressure,Rn; alpha=1.26)
    # compare to R
    @test ≈(ET_pot, 0.0002035969; rtol = 1e-5)
    @test ≈(LE_pot, 494.7202; rtol = 1e-5)
    #
    VPD,Ga_h = 2.0, 0.1
    Gs_pot = 0.5 # assume known surface conductance
    ET_pot,LE_pot = potential_ET(PenmanMonteith(), Tair,pressure,Rn,VPD,Ga_h; Gs_pot)
    Gs_ms, Gs_mol = surface_conductance(InversePenmanMonteith(), Tair,pressure,VPD,LE_pot,Rn,Ga_h)
    @test Gs_mol ≈ Gs_pot
end

@testset "potential_ET dataframe" begin
    dfo = DataFrame(Tair = 20.0:1.0:30.0,pressure = 100.0, Rn = 500.0)
    df = copy(dfo)
    potential_ET!(df, PriestleyTaylor(); alpha=1.26)
    @test ≈(last(df).ET_pot, 0.0002035969; rtol = 1e-5)
    @test ≈(last(df).LE_pot, 494.7202; rtol = 1e-5)
    #
    df2 = copy(dfo)
    df2[!, :VPD] .= 2.0
    df2[!, :Ga_h] .= 0.1
    potential_ET!(df2, PenmanMonteith(); Gs_pot = 0.4)
    df2b = rename!(copy(df2, copycols=false),SA[:LE_pot => :LE])
    surface_conductance!(df2b, InversePenmanMonteith())
    @test all(df2b.Gs_mol .≈ 0.4)
end

@testset "potential_ET dataframe missings in G" begin
    df = @pipe DataFrame(Tair = 20.0,pressure = 100.0, Rn = 500.0, G = 100.0:1:110) |>
      allowmissing(_, Cols(:G))
    df.G[2] = missing
    df_ET = potential_ET!(copy(df), PriestleyTaylor(); G = df.G)
    @test ncol(df_ET) == ncol(df)+2 # two columns added: ET_pot, LE_pot
    @test nrow(df_ET) == nrow(df)
    @test df_ET.ET_pot[end] < df_ET.ET_pot[1] # smaller: mith higher G less available energy
    @test df_ET.ET_pot[1] < 0.0002035969 # smaller as without G because less energy available
    @test ismissing(df_ET.LE_pot[2])
    #
    df2 = copy(df)
    df2[!, :VPD] .= 2.0
    df2[!, :Ga_h] .= 0.1
    potential_ET!(df2, PenmanMonteith(), S = df.G; Gs_pot = 0.3)
    @test ncol(df2) == ncol(df)+2+2 # two columns added: ET_pot, LE_pot
    @test nrow(df2) == nrow(df)
    @test ismissing(df2.LE_pot[2])
    df2b = rename!(copy(df2, copycols=false),SA[:LE_pot => :LE])
    surface_conductance!(df2b, InversePenmanMonteith(); S = df.G)
    @test ismissing(df2b.Gs_mol[2])
    @test all(df2b.Gs_mol[1:end .!= 2] .≈ 0.3)
end

@testset "equilibrium_imposed_ET scalars" begin
    # regression from R package example
    Tair,pressure,Rn, VPD, Gs = 20.0,100.0,50.0, 0.5, 0.01
    ET_eq, ET_imp, LE_eq, LE_imp = equilibrium_imposed_ET(Tair,pressure,VPD,Gs, Rn)
    @test ≈(ET_eq, 1.399424e-05; rtol = 1e-5)
    @test ≈(ET_imp, 3.695727e-05; rtol = 1e-5)
    @test ≈(LE_eq, 34.33628; rtol = 1e-5)
    @test ≈(LE_imp, 90.67837; rtol = 1e-5)
    #
    df = DataFrame(Tair = 20.0:1.0:22.0, pressure = 100.0, Rn = 50.0,  VPD = 0.5,  Gs = 0.01)
    ncoldf0 = ncol(df)
    # df_ET = @test_logs (:info,r"G is not provided") (:info,r"S is not provided") equilibrium_imposed_ET(df)
    # @test ncol(df) == ncoldf0
    # @test nrow(df_ET) == nrow(df)
    # @test ≈(first(df_ET).ET_eq, ET_eq)
    # @test ≈(first(df_ET).ET_imp, ET_imp)
    #
    dfm = copy(df)
    equilibrium_imposed_ET!(dfm)
    @test ncol(dfm) == ncoldf0 + 4
    @test nrow(dfm) == nrow(df)
    @test ≈(first(dfm).ET_eq, ET_eq)
    @test ≈(first(dfm).ET_imp, ET_imp)
end


text:
using Test
using DataFlowTasks
using DataFlowTasks.TiledFactorization
using LinearAlgebra
using BenchmarkTools
using VectorizationBase
using RecursiveFactorization

nc = min(Int(VectorizationBase.num_cores()), Threads.nthreads())

sch = DataFlowTasks.JuliaScheduler(20)
DataFlowTasks.setscheduler!(sch)

m         = 4000
# create an SPD matrix
A = rand(m,m)
# A = (A + adjoint(A))/2
# A = A + m*I
nt = Threads.nthreads()

@info "OpenBLAS"
@info BLAS.get_config()
BLAS.set_num_threads(nc)
@info "NBLAS = $(BLAS.get_num_threads())"
tnat = @btime lu!(B) setup=(B=copy(A)) evals=1

@info "DataFlowTasks"
BLAS.set_num_threads(1)
@info "NBLAS = $(BLAS.get_num_threads())"
tdf = @btime TiledFactorization.lu!(B) setup=(B=copy(A)) evals=1
# tfj = @btime TiledFactorization._cholesky_forkjoin!(B) setup=(B=PseudoTiledMatrix(copy(A),150))

DataFlowTasks.TASKCOUNTER[] = 0

# F = TiledFactorization._cholesky_forkjoin!(PseudoTiledMatrix(copy(A),256))
F = TiledFactorization.lu(A)

BLAS.set_num_threads(8)
er = norm(F.L*F.U-A,Inf)

F = RecursiveFactorization.lu(A,NoPivot())
er_blas = norm(F.L*F.U-A,Inf)

using MKL
@info "MKL cholesky"
@info BLAS.get_config()
BLAS.set_num_threads(nc)
@info "NBLAS = $(BLAS.get_num_threads())"
tnat = @btime cholesky!(B) setup=(B=copy(A)) evals=1

# for m×m tiled matrix, there should be 1/6*(m^3-m) + 1/2*(m^2-m) + m tasks created
@info "Number of tasks = $(DataFlowTasks.TASKCOUNTER[])"
@info "Number of threads = $(nt)"
@info "Number of cores = $(nc)"
@info er,er_blas


text:
import ClimateMachine.Atmos: atmos_init_aux!, vars_state

Base.@kwdef struct IsentropicVortexSetup{FT}
    p∞::FT = 10^5
    T∞::FT = 300
    ρ∞::FT = air_density(param_set, FT(T∞), FT(p∞))
    translation_speed::FT = 150
    translation_angle::FT = pi / 4
    vortex_speed::FT = 50
    vortex_radius::FT = 1 // 200
    domain_halflength::FT = 1 // 20
end

function (setup::IsentropicVortexSetup)(
    problem,
    bl,
    state,
    aux,
    localgeo,
    t,
    args...,
)
    FT = eltype(state)
    x = MVector(localgeo.coord)
    param_set = parameter_set(bl)

    ρ∞ = setup.ρ∞
    p∞ = setup.p∞
    T∞ = setup.T∞
    translation_speed = setup.translation_speed
    α = setup.translation_angle
    vortex_speed = setup.vortex_speed
    R = setup.vortex_radius
    L = setup.domain_halflength

    u∞ = SVector(translation_speed * cos(α), translation_speed * sin(α), 0)

    x .-= u∞ * t
    # make the function periodic
    x .-= floor.((x .+ L) / 2L) * 2L

    @inbounds begin
        r = sqrt(x[1]^2 + x[2]^2)
        δu_x = -vortex_speed * x[2] / R * exp(-(r / R)^2 / 2)
        δu_y = vortex_speed * x[1] / R * exp(-(r / R)^2 / 2)
    end
    u = u∞ .+ SVector(δu_x, δu_y, 0)

    _kappa_d::FT = kappa_d(param_set)
    T = T∞ * (1 - _kappa_d * vortex_speed^2 / 2 * ρ∞ / p∞ * exp(-(r / R)^2))
    # adiabatic/isentropic relation
    p = p∞ * (T / T∞)^(FT(1) / _kappa_d)
    ρ = air_density(param_set, T, p)

    state.ρ = ρ
    state.ρu = ρ * u
    e_kin = u' * u / 2
    state.energy.ρe = ρ * total_energy(param_set, e_kin, FT(0), T)
    if !(moisture_model(bl) isa DryModel)
        state.moisture.ρq_tot = FT(0)
    end
end

struct IsentropicVortexReferenceState{FT} <: ReferenceState
    setup::IsentropicVortexSetup{FT}
end
vars_state(::IsentropicVortexReferenceState, ::Auxiliary, FT) =
    @vars(ρ::FT, ρe::FT, p::FT, T::FT)
function atmos_init_aux!(
    atmos::AtmosModel,
    m::IsentropicVortexReferenceState,
    state_auxiliary::MPIStateArray,
    grid,
    direction,
)
    init_state_auxiliary!(
        atmos,
        (args...) -> init_vortex_ref_state!(m, args...),
        state_auxiliary,
        grid,
        direction,
    )
end
function init_vortex_ref_state!(
    m::IsentropicVortexReferenceState,
    atmos::AtmosModel,
    aux::Vars,
    tmp::Vars,
    geom::LocalGeometry,
)
    setup = m.setup
    ρ∞ = setup.ρ∞
    p∞ = setup.p∞
    T∞ = setup.T∞
    param_set = parameter_set(atmos)
    aux.ref_state.ρ = ρ∞
    aux.ref_state.p = p∞
    aux.ref_state.T = T∞
    aux.ref_state.ρe = ρ∞ * internal_energy(param_set, T∞)
end