text:
# This module implements the MMA Algorithm in Julia
# as described in:
# AU  - Svanberg, Krister
# TI  - The method of moving asymptotes—a new method for structural optimization
# JO  - International Journal for Numerical Methods in Engineering
# JA  - Int. J. Numer. Meth. Engng.
module MMA

import Optim: OnceDifferentiable, Fminbox, GradientDescent, update!,
                MultivariateOptimizationResults, OptimizationTrace, Optimizer, optimize, maxdiff
import Base: min, max, show

export MMAModel, box!, ineq_constraint!, optimize

include("utils.jl")

macro mmatrace()
    esc(quote
        if tracing
            dt = Dict()
            if m.extended_trace
                dt["x"] = copy(x)
                dt["g(x)"] = copy(∇f_x)
                dt["λ"] = copy(results.minimizer)
            end
            update!(tr,
                    k,
                    f_x,
                    gr_residual,
                    dt,
                    m.store_trace,
                    m.show_trace)
        end
    end)
end

struct MMA87 <: Optimizer end

struct MMAModel
    dim::Int
    objective::Function
    ineq_constraints::Vector{Function}
    box_max::Vector{Float64}
    box_min::Vector{Float64}

    store_trace::Bool
    show_trace::Bool
    extended_trace::Bool

    # Stopping criterias
    max_iters::Int
    ftol::Float64
    xtol::Float64
    grtol::Float64
end

mutable struct DualData
    L::Vector{Float64}
    U::Vector{Float64}
    α::Vector{Float64} # Lower move limit
    β::Vector{Float64} # Upper move limit
    p0::Vector{Float64}
    q0::Vector{Float64}
    p::Matrix{Float64}
    q::Matrix{Float64}
    r::Vector{Float64}
    r0::Float64
    x::Vector{Float64} # Optimal value for x in dual iteration
    f_val::Float64 # Function value at current iteration
    g_val::Vector{Float64} # Inequality values at current iteration
    ∇f::Vector{Float64} # Function gradient at current iteration
    ∇g::Vector{Vector{Float64}} # Inequality gradients [ineq][var] at current iteration
end

# Inspired by Parameters.jl
const dd_fields = [:L, :U, :α, :β, :p0, :q0, :p, :q,
                   :r, :r0, :x, :f_val, :g_val, :∇f, :∇g]
macro unpack(dual_data)
    esc(Expr(:block, [:($f = dual_data.$f) for f in dd_fields]...))
end

function show(io::IO, dual_data::DualData)
    println("Dual data:")
    for f in dd_fields
        println("$f = $(dual_data.(f))")
    end
end

dim(m::MMAModel) = m.dim
min(m::MMAModel, i::Int) = m.box_min[i]
max(m::MMAModel, i::Int) = m.box_max[i]
min(m::MMAModel) = m.box_max
max(m::MMAModel) = m.box_min
objective(m::MMAModel) = m.objective
constraints(m::MMAModel) = m.ineq_constraints
constraint(m::MMAModel, i::Int) = m.ineq_constraints[i]
eval_objective(m, g0, ∇g0) = m.objective(g0, ∇g0)
eval_objective(m, g0) = m.objective(g0, Float64[])
eval_constraint(m, i, g0, ∇g0) = constraint(m, i)(g0, ∇g0)
eval_constraint(m, i, g0) = constraint(m, i)(g0, Float64[])
ftol(m) = m.ftol
xtol(m) = m.xtol
grtol(m) = m.grtol
ftol!(m, v) = m.ftol = v
xtol!(m, v) = m.xtol = v
grtol!(m, v) = m.grtol = v

function MMAModel(dim::Int,
                  objective::Function;
                  max_iters::Int = 200,
                  xtol::Float64 = eps(Float64),
                  ftol::Float64 = sqrt(eps(Float64)),
                  grtol::Float64 = sqrt(eps(Float64)),
                  store_trace::Bool = false,
                  show_trace::Bool = false,
                  extended_trace::Bool = false)

    mins = fill(Float64(Inf), dim)
    maxs = fill(Float64(-Inf), dim)
    MMAModel(dim, objective, Function[],
             mins, maxs, store_trace, show_trace, extended_trace,
             max_iters, ftol, xtol, grtol)
end

# Box constraints
function box!(m::MMAModel, i::Int, minb::Float64, maxb::Float64)
    if !(1 <= i <= dim(m))
        throw(ArgumentError("box constraint need to applied to an existing variable"))
    end
    m.box_min[i] = minb
    m.box_max[i] = maxb
end

function box!(m::MMAModel, minb::Float64, maxb::Float64)
    for i = 1:dim(m)
        m.box_min[i] = minb
        m.box_max[i] = maxb
    end
end

function box!(m::MMAModel, minbs::Vector{Float64}, maxbs::Vector{Float64})
    if (length(minbs) != dim(m)) || (length(minbs) != dim(m))
        throw(ArgumentError("box constraint vector must have same size as problem dimension"))
    end
    for i = 1:dim(m)
        m.box_min[i] = minbs[i]
        m.box_max[i] = maxbs[i]
    end
end

function ineq_constraint!(m::MMAModel, f::Function)
    push!(m.ineq_constraints, f)
end

function ineq_constraint!(m::MMAModel, fs::Vector{Function})
    for f in fs
        push!(m.ineq_constraints, f)
    end
end

function optimize(m::MMAModel, x0::Vector{Float64})
    check_error(m, x0)
    n_i = length(constraints(m))
    n_j = dim(m)
    x0, x, x1, x2 = copy(x0), copy(x0), copy(x0), copy(x0)

    # Buffers for bounds and move limits
    α, β, L, U = zeros(Float64, n_j), zeros(Float64, n_j), zeros(Float64, n_j), zeros(Float64, n_j)

    # Buffers for p0, pij, q0, qij
    p0, q0 = zeros(Float64, n_j), zeros(Float64, n_j)
    p, q = zeros(Float64, n_i, n_j), zeros(Float64, n_i, n_j)
    r = zeros(Float64, n_i)

    # Initial data and bounds for Optim to solve dual problem
    λ = ones(Float64, n_i)
    l = zeros(Float64, n_i)
    u = fill(Float64(Inf), n_i)

    # Objective gradient buffer
    ∇f_x = similar(x)
    f_x = eval_objective(m, x, ∇f_x)
    f_calls, g_calls = 1, 1
    f_x_previous = Float64(NaN)

    # Evaluate the constraints and their gradients
    g = zeros(n_i)
    ∇g = [zeros(n_j) for _ in 1:n_i]
    for i = 1:n_i
        g[i] = eval_constraint(m, i, x, ∇g[i])
    end

    # Create a DualData type that holds the data needed for the dual problem
    dual_data = DualData(L, U, α, β, p0, q0, p, q, r, 0.0, x, f_x, g, ∇f_x, ∇g)

    tr = OptimizationTrace{MMA87}()
    tracing = (m.store_trace || m.extended_trace || m.show_trace)

    converged = false
    # Assess multiple types of convergence
    x_converged, f_converged, gr_converged = false, false, false
    f_increased = false
    x_residual = Float64(Inf)
    f_residual = Float64(Inf)
    gr_residual = Float64(Inf)

    # Iteraton counter
    k = 0
    while !converged && k < m.max_iters
        k += 1
        update_limits!(dual_data, m, k, x1, x2)

        # Track trial points two steps back
        copy!(x2, x1)
        copy!(x1, x)

        compute_mma!(dual_data, m)

        dual(λ) = compute_dual!(λ, dual_data) #Lagrangian dual objective value
        dual_grad(grad_dual, λ) = compute_dual_grad!(grad_dual, λ, dual_data)
        d = OnceDifferentiable(dual, dual_grad, λ)
        results = optimize(d, λ, l, u, Fminbox{GradientDescent}())
        # Use previously converged λ as starting guess
        copy!(λ, results.minimizer)
        update_x!(dual_data, λ)
        f_x_previous, f_x = f_x, eval_objective(m, x, ∇f_x)
        f_calls, g_calls = f_calls + 1, g_calls + 1
         # Evaluate the constraints and their gradients
        for i = 1:n_i
            g[i] = eval_constraint(m, i, x, ∇g[i])
        end

        x_converged, f_converged, gr_converged,
        x_residual, f_residual, gr_residual,
        f_increased, converged = assess_convergence(x, x1, f_x, f_x_previous, ∇f_x,
                                                     xtol(m), ftol(m), grtol(m))

        @mmatrace()
        if converged
            break
        end
    end
    h_calls = 0
    return MultivariateOptimizationResults(MMA87(),
                                           false,
                                           x0,
                                           x,
                                           f_x,
                                           k,
                                           k == m.max_iters,
                                           x_converged,
                                           xtol(m),
                                           x_residual,
                                           f_converged,
                                           ftol(m),
                                           f_residual,
                                           gr_converged,
                                           grtol(m),
                                           gr_residual,
                                           f_increased,
                                           tr,
                                           f_calls,
                                           g_calls,
                                           h_calls)
end

# Updates p0, q0, p, q, r and returns r0.
# For notation see reference at top
function compute_mma!(dual_data, m)
    @unpack dual_data
    # Bound limits
    for j = 1:dim(m)
        μ = 0.1
        α[j] = max(L[j] + μ * (x[j] - L[j]), min(m, j))
        β[j] = min(U[j] - μ * (U[j] - x[j]), max(m, j))
    end

    r0 = 0.0
    for i in 0:length(constraints(m))
        if i == 0
            ri = f_val
            ∇fi = ∇f
        else
             ri = g_val[i]
             ∇fi = ∇g[i]
        end
        for j in 1:dim(m)
            Ujxj = U[j] - x[j]
            xjLj = x[j] - L[j]
            if ∇fi[j] > 0
                p_ij = abs2(Ujxj) * ∇fi[j]
                q_ij = 0.0
            else
                p_ij = 0.0
                q_ij = -abs2(xjLj) * ∇fi[j]
            end
            ri -= p_ij / Ujxj + q_ij / xjLj
            if i == 0
                p0[j] = p_ij
                q0[j] = q_ij
            else
                p[i, j] = p_ij
                q[i, j] = q_ij
            end
        end
        if i == 0
            r0 = ri
        else
            r[i] = ri
        end
    end
    dual_data.r0 = r0
end

# Update move limits

const s_init = 0.5
# Suggested by Svanberg
const s_incr = 1.05
const s_decr = 0.65

function update_limits!(dual_data, m, k, x1, x2)
    @unpack dual_data
    for j in 1:dim(m)
        if k == 1 || k == 2
            # Equation 11 in Svanberg
            L[j] = x[j] - s_init * (max(m,j) - min(m, j))
            U[j] = x[j] + s_init * (max(m,j) - min(m, j))
        else
            # Equation 12 in Svanberg
            if sign(x[j] - x1[j]) != sign(x1[j] - x2[j])
                L[j] = x[j] - (x1[j] - L[j]) * s_decr
                U[j] = x[j] + (U[j] - x1[j]) * s_decr
            # Equation 13 in Svanberg
            else
                L[j] = x[j] - (x1[j] - L[j]) * s_incr
                U[j] = x[j] + (U[j] - x1[j]) * s_incr
            end
        end
    end
end

function matdot(A::Vector, B::Matrix, j::Int)
    r = 0.0
    @inbounds for i in eachindex(A)
        r += A[i] * B[i, j]
    end
    return r
end

function compute_dual!(λ, dual_data)
    @unpack dual_data
    update_x!(dual_data, λ)

    #Optimal value of Lagrangian at λ
    φ = r0 + dot(λ, r)
    @inbounds for j = 1:length(x)
        φ += (p0[j] + matdot(λ, p, j)) / (U[j] - x[j])
        φ += (q0[j] + matdot(λ, q, j)) / (x[j] - L[j])
    end
    return -φ
end

function compute_dual_grad!(∇φ, λ, dual_data)
    @unpack dual_data
    update_x!(dual_data, λ)
    for i = 1:length(λ)
        ∇φ[i] = r[i]
        for j = 1:length(x)
            ∇φ[i] += p[i,j] / (U[j] - x[j])
            ∇φ[i] += q[i,j] / (x[j] - L[j])
        end
    end
    # Negate since we have a maximization problem
    scale!(∇φ, -1.0)
    return
end

# Updates x to be the analytical optimal point in the dual
# problem for a given λ
function update_x!(dual_data, λ)
    @unpack dual_data
    @inbounds for j in 1:length(x)
        fpj = sqrt(p0[j] + matdot(λ, p, j))
        fqj = sqrt(q0[j] + matdot(λ, q, j))
        x[j] = (fpj * L[j] + fqj * U[j]) / (fpj + fqj)
        if x[j] > β[j]
            x[j] = β[j]
        elseif x[j] < α[j]
            x[j] = α[j]
        end
    end
    return
end

end # module


text:
function is_orthogonal(S, tol=1e-6)
    n = LinearAlgebra.checksquare(S)
    for i in 1:n
        for j in 1:n
            if i != j
                s = LinearAlgebra.dot(S[:, i], S[:, j])
                if abs(s) > tol
                    return false
                end
            end
        end
    end
    return true
end

function _projcoef(A, i, v, w = A[i, :])
    return LinearAlgebra.dot(v, w) / LinearAlgebra.dot(w, w)
end
function _proj(A, i, v)
    w = A[i, :]
    return w * _projcoef(A, i, v, w)
end

function orthogonalize(A)
    display(A)
    for i in 2:size(A, 1)
        row = A[i, :]
        for j in 1:(i - 1)
            row = row - _proj(A, j, row)
        end
        A[i, :] = row
    end
    return A
end

struct _OrthogonalMatrix end
struct _RowEchelonMatrix end

function __linsolve(A::Matrix{T}, b::Vector{T}, ::_OrthogonalMatrix) where {T}
    return map(1:size(A, 1)) do i
        return _projcoef(A, i, b)
    end
end

function __linsolve(A::Matrix{T}, b::Vector{T}, ::_RowEchelonMatrix) where {T}
    j = 0
    return map(1:size(A, 1)) do i
        while j < size(A, 2)
            j += 1
            if isone(A[i, j]) && all(k -> k == i || iszero(A[k, j]), 1:size(A, 1))
                return b[j]
            end
        end
        error("Not in row_echelon_form, cannot find for `$i`th entry.")
    end
end

function _linsolve(A::Matrix{T}, b::Vector{T}, form) where {T}
    x = __linsolve(A, b, form)
    #@assert transpose(A) * x == b
    return x
end

ext: jl


text:
using Documenter, SIMDIndexing

makedocs(;
    modules=[SIMDIndexing],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/tkf/SIMDIndexing.jl/blob/{commit}{path}#L{line}",
    sitename="SIMDIndexing.jl",
    authors="Takafumi Arakaki",
    assets=[],
)

deploydocs(;
    repo="github.com/tkf/SIMDIndexing.jl",
)


text:
#Test for Gadfly with different chunk options and figure formatsusing Weave
using Weave
using Test

weave("documents/gadfly_formats_test.txt", doctype = "tex", plotlib="gadfly")
result =  read("documents/gadfly_formats_test.tex", String)
ref =  read("documents/gadfly_formats_test_ref.tex", String)
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="tex", fig_ext=".tex", plotlib="gadfly")
result =  read("documents/gadfly_formats_test.tex", String)
ref =  read("documents/gadfly_formats_test_tikz_ref.tex", String)
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="tex", fig_ext=".ps", plotlib="gadfly")
result =  read("documents/gadfly_formats_test.tex", String)
ref =  read("documents/gadfly_formats_test_ps_ref.tex", String)
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="pandoc", plotlib="gadfly")
result =  read("documents/gadfly_formats_test.md", String)
ref =  read("documents/gadfly_formats_test_pandoc_ref.md", String)
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="pandoc", plotlib="gadfly", fig_ext=".svg")
result =  read("documents/gadfly_formats_test.md", String)
ref =  read("documents/gadfly_formats_test_svg_ref.md", String)
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="github", plotlib="gadfly", fig_ext=".js.svg")
result =  read("documents/gadfly_formats_test.md", String)
ref =  read("documents/gadfly_formats_test_jssvg_ref.md", String)
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="rst", plotlib="gadfly")
result =  read("documents/gadfly_formats_test.rst", String)
ref =  read("documents/gadfly_formats_test_ref.rst", String)
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="multimarkdown", plotlib="gadfly")
result =  read("documents/gadfly_formats_test.md", String)
ref =  read("documents/gadfly_formats_test_mmd_ref.md", String)
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="asciidoc", plotlib="gadfly",
    out_path="documents/output")
result =  read("documents/output/gadfly_formats_test.txt", String)
ref =  read("documents/output/gadfly_formats_test_ref.txt", String)
@test result == ref

weave("documents/gadfly_markdown_test.jmd", doctype="github",plotlib="gadfly", informat="markdown")
result =  read("documents/gadfly_markdown_test.md", String)
ref =  read("documents/gadfly_markdown_test_ref.md", String)
@test result == ref

weave("documents/FIR_design.jl", doctype="pandoc", plotlib="gadfly", informat="script")
result =  read("documents/FIR_design.md", String)
ref =  read("documents/FIR_design_ref.md", String)
@test result == ref


text:
module ElectronProgressBars

import JSON
using ArgCheck: @argcheck
using Electron: Window, load
using Logging: Logging, global_logger
using Printf: @sprintf

include("router.jl")

function locking(f, l)
    lock(l)
    try
        return f()
    finally
        unlock(l)
    end
end

struct Vault{T, L}
    value::Base.RefValue{T}
    lock::L
end

Vault(ref::Ref) = Vault(ref, ReentrantLock())
Vault(value) = Vault(Ref(value))
Vault{T}() where T = Vault(Ref{T}())
Vault{T}(value) where T = Vault(Ref{T}(value))

locking(f, v::Vault) = locking(() -> f(v.value), v.lock)
Base.getindex(v::Vault) = locking(value -> value[], v)
Base.setindex!(v::Vault, x) = locking(value -> value[] = x, v)

_taskid() = objectid(current_task())

const TaskID = UInt64
const LogID = Symbol

struct ProgressBarWindow
    bars::Dict{TaskID, Vector{LogID}}
    starttime::Dict{Tuple{TaskID, LogID}, Float64}
    window::Window
end

ProgressBarWindow() = ProgressBarWindow(
    Dict(),
    Dict(),
    Window(
        main_html(),
        options = Dict(
            "title" => "Progress bars",
            "webPreferences" => Dict("webSecurity" => false),
        )
    ),
)

_isopen(::Nothing) = false
_isopen(w::ProgressBarWindow) = w.window.exists

mutable struct ElectronProgressBarLogger <: Logging.AbstractLogger
    window::ProgressBarWindow
    channel::Channel{Any}
end

function message_handler(ch)
    for f in ch
        try
            f()
        catch err
            @error "Error from function: $f" exception=(err, catch_backtrace())
        end
    end
end

inloggingtask(f, logger::ElectronProgressBarLogger) = put!(logger.channel, f)

ElectronProgressBarLogger() = ElectronProgressBarLogger(
    ProgressBarWindow(),
    Channel(message_handler),  # buffered?
)

function reopen!(logger::ElectronProgressBarLogger)
    new = ElectronProgressBarLogger()
    logger.window = new.window
    logger.channel = new.channel
    return logger
end

function Base.close(logger::ElectronProgressBarLogger)
    close(logger.window.window)
    close(logger.channel)
end

const _singleton_logger = Vault{Union{Nothing, ElectronProgressBarLogger}}(nothing)

singleton_logger() = locking(_singleton_logger) do value
    logger = value[]
    if logger === nothing
        value[] = logger = ElectronProgressBarLogger()
    elseif !_isopen(logger.window)
        reopen!(logger)
    end
    return logger
end

get_logger() = singleton_logger()

close_singleton_logger() = locking(_singleton_logger) do value
    logger = value[]
    if logger !== nothing && _isopen(logger.window)
        close(logger)
    end
end

function reopen_singleton_logger()
    close_singleton_logger()
    return singleton_logger()
end

_reload() =
    load(singleton_logger().window.window, main_html())

function main_html()
    main_json = joinpath(@__DIR__, "main.js")
    main_css = joinpath(@__DIR__, "main.css")
    return """

            <link rel="stylesheet" href="file://$main_css">

    """
end

barid(taskid::TaskID, logid::LogID) = "bar-$taskid-$logid"

function durationstring(eta)
    isfinite(eta) || return "?"
    s = floor(Int, eta)
    m, s = divrem(s, 60)
    h, m = divrem(m, 60)
    d, h = divrem(h, 24)
    if d >= 1
        dunit = d == 1 ? "day" : "days"
        hunit = h == 1 ? "hour" : "hours"
        h == 0 && return "$d $dunit"
        return "$d $dunit $h $hunit"
    else
        return @sprintf "%02d:%02d:%02d" h m s
    end
end

# for debugging
_gset_progress(
    progress::Real,
    logid::LogID = :dummy;
    title::AbstractString = "title",
    message::Union{AbstractString, Nothing} = nothing,
    kwargs...,
) =
    _set_progress(singleton_logger(), logid, title, progress, message; kwargs...)

function _set_progress(
    logger::ElectronProgressBarLogger,
    logid::LogID,
    title::AbstractString,
    progress::Real,
    message::Union{AbstractString, Nothing};
    taskid::TaskID = _taskid(),
    now::Float64 = time(),
)
    @argcheck 0 <= progress <= 1

    inloggingtask(logger) do
        w = logger.window
        param = (
            logid = logid,
            taskid = taskid,
            barid = barid(taskid, logid),
            title = title,
            message = message,
            progress = progress,
            progresstext = (@sprintf "%02.1f %%" 100 * progress),
        )

        bars = get!(w.bars, taskid, [])
        if isempty(bars)
            push!(bars, logid)
            w.starttime[taskid, logid] = now
            p = JSON.json(param)
            run(w.window, "newRootBar($p)")
            return
        end

        idx = findfirst(==(logid), bars)
        if idx === nothing
            param = (
                parentid = barid(taskid, bars[end]),
                param...,
            )
            push!(bars, logid)
            w.starttime[taskid, logid] = now
            p = JSON.json(param)
            run(w.window, "newSubBar($p)")
        else
            duration = now - w.starttime[taskid, logid]
            eta = duration * (1 / progress - 1)
            param = (
                finished = barid.(taskid, bars[idx + 1:end]),
                eta = eta,
                etatext = string("ETA ", durationstring(eta)),
                param...,
            )
            deleteat!(bars, idx + 1:length(bars))
            p = JSON.json(param)
            run(w.window, "setProgress($p)")
        end
    end
    return
end

# for debugging
_gremove_progress(logid::LogID = :dummy; kwargs...) =
    _remove_progress(singleton_logger(), logid; kwargs...)

function _remove_progress(
    logger::ElectronProgressBarLogger,
    logid::LogID;
    taskid::TaskID = _taskid(),
)
    inloggingtask(logger) do
        w = logger.window
        bars = get!(w.bars, taskid, [])
        isempty(bars) && return
        idx = findfirst(==(logid), bars)
        idx === nothing && return
        param = (
            logid = logid,
            taskid = taskid,
            barid = barid(taskid, logid),
            finished = barid.(taskid, bars[idx:end]),
        )
        deleteat!(bars, idx:length(bars))
        p = JSON.json(param)
        run(w.window, "removeFinished($p)")
    end
    return
end

# https://docs.julialang.org/en/latest/stdlib/Logging/#AbstractLogger-interface-1

function Logging.handle_message(
    logger::ElectronProgressBarLogger,
    level, title, _module, group, id, file, line;
    progress = nothing,
    message = nothing,
    _...
)
    if progress isa Real && (progress < 1 || isnan(progress))
        progress = isnan(progress) ? 0.0 : progress
        _set_progress(logger, id, title, max(0.0, progress), message)
    elseif progress == "done" || (progress isa Real && progress >= 1)
        _remove_progress(logger, id)
    end
    return nothing
end

Logging.shouldlog(::ElectronProgressBarLogger, level, _module, group, id) =
    true

Logging.min_enabled_level(::ElectronProgressBarLogger) = Logging.BelowMinLevel

"""
    install_logger()

Install `ElectronProgressBarLogger` to global logger.
"""
install_logger() = install_logger(singleton_logger())
function install_logger(logger::ElectronProgressBarLogger)
    global previous_logger
    previous_logger = global_logger(ProgressLogRouter(global_logger(), logger))
end

"""
    uninstall_logger()

Rollback the global logger to the one before last call of `install_logger`.
"""
function uninstall_logger()
    global previous_logger
    previous_logger === nothing && return
    ans = global_logger(previous_logger)
    previous_logger = nothing
    return ans
end

end # module


text:
# !! auto-generated code, do not edit !!
module geometry_msgs
using RobotOSData.Messages
using RobotOSData.StdMsgs
using RobotOSData.CommonMsgs
struct Vector3 <: Readable
    x::Float64
    y::Float64
    z::Float64
end
struct Vector3Stamped <: Readable
    header::Header
    vector::Vector3
end
struct Accel <: Readable
    linear::Vector3
    angular::Vector3
end
struct AccelStamped <: Readable
    header::Header
    accel::Accel
end
struct Point32 <: Readable
    x::Float32
    y::Float32
    z::Float32
end
struct Polygon <: Readable
    points::Vector{Point32}
end
struct PolygonStamped <: Readable
    header::Header
    polygon::Polygon
end
struct Quaternion <: Readable
    x::Float64
    y::Float64
    z::Float64
    w::Float64
end
struct Transform <: Readable
    translation::Vector3
    rotation::Quaternion
end
struct Point <: Readable
    x::Float64
    y::Float64
    z::Float64
end
struct Pose <: Readable
    position::Point
    orientation::Quaternion
end
struct PoseWithCovariance <: Readable
    pose::Pose
    covariance::SVector{36, Float64}
end
struct PoseWithCovarianceStamped <: Readable
    header::Header
    pose::PoseWithCovariance
end
struct AccelWithCovariance <: Readable
    accel::Accel
    covariance::SVector{36, Float64}
end
struct PointStamped <: Readable
    header::Header
    point::Point
end
struct Twist <: Readable
    linear::Vector3
    angular::Vector3
end
struct TwistStamped <: Readable
    header::Header
    twist::Twist
end
struct AccelWithCovarianceStamped <: Readable
    header::Header
    accel::AccelWithCovariance
end
struct TwistWithCovariance <: Readable
    twist::Twist
    covariance::SVector{36, Float64}
end
struct Inertia <: Readable
    m::Float64
    com::geometry_msgs.Vector3
    ixx::Float64
    ixy::Float64
    ixz::Float64
    iyy::Float64
    iyz::Float64
    izz::Float64
end
struct Pose2D <: Readable
    x::Float64
    y::Float64
    theta::Float64
end
struct TwistWithCovarianceStamped <: Readable
    header::Header
    twist::TwistWithCovariance
end
struct PoseArray <: Readable
    header::Header
    poses::Vector{Pose}
end
struct Wrench <: Readable
    force::Vector3
    torque::Vector3
end
struct TransformStamped <: Readable
    header::Header
    child_frame_id::String
    transform::Transform
end
struct WrenchStamped <: Readable
    header::Header
    wrench::Wrench
end
struct PoseStamped <: Readable
    header::Header
    pose::Pose
end
struct InertiaStamped <: Readable
    header::Header
    inertia::Inertia
end
struct QuaternionStamped <: Readable
    header::Header
    quaternion::Quaternion
end
end


text:
function check_status(status::TStatus)
    errormsg = isfilled(status, :errorMessage) ? getfield(status, :errorMessage) : Compat.UTF8String("")
    infomsgs = isfilled(status, :infoMessages) ? getfield(status, :infoMessages) : Compat.UTF8String[]
    check_status(status.statusCode, errormsg, infomsgs)
end
function check_status(status::Int32, errormsg::AbstractString="", infomsgs::Array=[])
    if status == TStatusCode.SUCCESS_STATUS || status == TStatusCode.SUCCESS_WITH_INFO_STATUS
        map((msg)->println("$(msg)\n"), infomsgs)
        return true
    end
    (status == TStatusCode.STILL_EXECUTING_STATUS) && (return false)
    (status == TStatusCode.ERROR_STATUS) && error("Hive operation failed. $errormsg")
    (status == TStatusCode.INVALID_HANDLE_STATUS) && error("Hive handle invalid. $errormsg")
    error("Unknown status code")
end

##
# Table schema methods
function dataframe(sch::TTableSchema)
    df = DataFrame()
    df[:position] = positions = DataArray(Int32[])
    df[:name] = names = DataArray(Compat.UTF8String[])
    df[:type] = types = DataArray(Type[])
    df[:comment] = comments = DataArray(Compat.UTF8String[])
    for col in sch.columns
        push!(comments, isfilled(col, :comment) ? getfield(col, :comment) : NA)
        push!(types, julia_type(col.typeDesc))
        push!(names, col.columnName)
        push!(positions, col.position)
    end
    df
end

show(io::IO, sch::TTableSchema) = show(io, dataframe(sch))

function coltypes(sch::TTableSchema)
    tuple([julia_type(col.typeDesc) for col in sch.columns]...)
end

const BIT_MASKS = (0x01, 0x02, 0x04, 0x08, 0x10, 0x20, 0x40, 0x80)
# the last byte of the bitset is always set to 0x00
function bitset_to_bools(data::Vector{UInt8}, L::Integer=((length(data)-1)*8))
    bools = Array(Bool, L)
    for idx in 1:L
        bools[idx] = (data[ceil(Int, idx/8)] & BIT_MASKS[rem(idx-1, 8) + 1] > 0x00)
    end
    bools
end


text:
# ---------------
# Value Iteration
# ---------------

facts("discounted dense value iteration") do
    P, R = MDPs.Examples.small()
    mdp = MDP(P, R)
    qfunc = value_iteration(mdp, 0.9)
    @fact policy(qfunc) == [2, 1] => true
    @fact value(qfunc) => roughly([40.048625392716815, 33.65371175967546])
end

facts("undiscounted dense value iteration") do
    P, R = MDPs.Examples.small()
    mdp = MDP(P, R)
    qfunc = value_iteration(mdp, 1)
    @fact policy(qfunc) == [2, 1] => true
    @fact value(qfunc) => roughly([117.84153597421569, 111.72677122062753])
end

facts("discounted sparse value iteration") do
    P, R = MDPs.Examples.sprandom(10, 3)
    mdp = MDP(P, R)
    qfunc = value_iteration(mdp, 0.9)
    @fact typeof(qfunc) => VectorQFunction{Float64,Int}
    @fact typeof(policy(qfunc)) => Vector{Int}
    @fact typeof(value(qfunc)) => Vector{Float64}
end

facts("MDP methods") do
    P = cat(3, eye(5), eye(5))
    R = [1 2 3 4 5; 5 4 3 2 1]'
    mdp = MDP(P, R)
    @fact reward(mdp, 2, 1) => 2
    @fact probability(mdp, 1, 1, 1) => 1
    @fact probability(mdp, 1, 2, 2) => 0
    @fact num_states(mdp) => 5
    @fact num_actions(mdp) => 2
end

facts("max iter too small") do
    warn(x) = "WARNING"
    P, R = MDPs.Examples.small()
    Q = value_iteration(MDP(P, R), 0.9, max_iter=5)
    @pending "test that there is a warning" => true
end


text:
@test_throws DimensionMismatch stratifiedobs((X, rand(149)))
@test_throws DimensionMismatch stratifiedobs((X, rand(149)), obsdim=:last)

Random.seed!(1335)

tmpty = [:a, :a, :a, :a, :a, :a, :b, :b, :b, :b]
@test sort.(splitobs(labelmap(tmpty), at = 0.5)) == ([1,2,3,7,8],[4,5,6,9,10])

@testset "Type Stability" begin
    for var in vars
        Random.seed!(1335)
        @test_throws ArgumentError stratifiedobs(var, 0.)
        @test_throws ArgumentError stratifiedobs(var, 1.)
        @test_throws ArgumentError stratifiedobs(var, (0.2,0.0))
        @test_throws ArgumentError stratifiedobs(var, (0.2,0.8))
        @test_throws MethodError stratifiedobs(var, 0.5, ObsDim.Undefined())
        @test typeof(@inferred(stratifiedobs(var))) <: NTuple{2}
        @test eltype(@inferred(stratifiedobs(var))) <: SubArray
        @test typeof(@inferred(stratifiedobs(var, 0.5))) <: NTuple{2}
        @test typeof(@inferred(stratifiedobs(var, 0.5, true))) <: NTuple{2}
        @test typeof(@inferred(stratifiedobs(var, (0.5,0.2)))) <: NTuple{3}
        @test typeof(@inferred(stratifiedobs(var, 0.5, true, ObsDim.Last()))) <: NTuple{2}
        @test typeof(@inferred(stratifiedobs(var, 0.5, true, ObsDim.First()))) <: NTuple{2}
        @test eltype(@inferred(stratifiedobs(var, 0.5, true, ObsDim.First()))) <: SubArray
    end
    for tup in tuples
        Random.seed!(1335)
        @test_throws ArgumentError stratifiedobs(tup, 0.)
        @test_throws ArgumentError stratifiedobs(tup, 1.)
        @test_throws ArgumentError stratifiedobs(tup, (0.2,0.0))
        @test_throws ArgumentError stratifiedobs(tup, (0.2,0.8))
        @test_throws MethodError stratifiedobs(tup, 0.5, ObsDim.Undefined())
        @test typeof(@inferred(stratifiedobs(tup))) <: NTuple{2}
        @test eltype(@inferred(stratifiedobs(tup))) <: Tuple
        @test typeof(@inferred(stratifiedobs(tup, 0.5))) <: NTuple{2}
        @test typeof(@inferred(stratifiedobs(tup, 0.5, true))) <: NTuple{2}
        @test typeof(@inferred(stratifiedobs(tup, (0.5,0.2)))) <: NTuple{3}
        @test typeof(@inferred(stratifiedobs(tup, 0.5, true, ObsDim.Last()))) <: NTuple{2}
    end
end

println("<HEARTBEAT>")

Random.seed!(1335)
@testset "SparseArray" begin
    @test nobs.(stratifiedobs(round, ys)) == (105,45)
    @test nobs.(stratifiedobs(round, ys, p=0.7)) == (105,45)
    @test nobs.(stratifiedobs(round, ys, p=(.2,.3))) == (30,45,75)
    @test nobs.(stratifiedobs(round, ys, p=(.2,.3), obsdim=:last)) == (30,45,75)
    @test nobs.(stratifiedobs(round, ys, p=(.1,.2,.3))) == (15,30,45,60)
end

@testset "Array, and SubArray" begin
    for var in (yv, y)
        @test nobs.(stratifiedobs(var)) == (105,45)
        @test nobs.(stratifiedobs(x->x, var)) == (105,45)
        @test nobs.(stratifiedobs(var, p=0.7)) == (105,45)
        @test nobs.(stratifiedobs(x->x, var, p=0.7)) == (105,45)
        @test nobs.(stratifiedobs(var, p=(.2,.3))) == (30,45,75)
        @test nobs.(stratifiedobs(var, p=(.2,.3), obsdim=:last)) == (30,45,75)
        @test nobs.(stratifiedobs(var, p=(.1,.2,.3))) == (15,30,45,60)
    end
    ty = [:a, :a, :a, :a, :a, :a, :b, :b, :b, :b]
    train, test = @inferred stratifiedobs(ty, 0.5, false)
    @test train == [:a, :a, :a, :b, :b] || train == [:b, :b, :a, :a, :a]
    @test test == [:a, :a, :a, :b, :b] || test == [:b, :b, :a, :a, :a]

    train, test = @inferred stratifiedobs(x->x, ty, 0.5, false)
    @test train == [:a, :a, :a, :b, :b] || train == [:b, :b, :a, :a, :a]
    @test test == [:a, :a, :a, :b, :b] || test == [:b, :b, :a, :a, :a]

    train, test, val = @inferred stratifiedobs(ty, (0.5,0.3), false)
    @test train == [:a, :a, :a, :b, :b] || train == [:b, :b, :a, :a, :a]
    @test test == [:a, :a, :b] || test == [:b, :a, :a]
    @test val == [:a, :b] || val == [:b, :a]

    ty = [:a, :a, :a, :a, :a, :a, :b, :b, :b, :b, :c, :c]
    train, test = @inferred stratifiedobs(ty, 0.5, false)
    @test train == [:c, :a, :a, :a, :b, :b] || train == [:a, :a, :a, :c, :b, :b] || train == [:b, :b, :c, :a, :a, :a] || train == [:c, :b, :b, :a, :a, :a] || train == [:b, :b, :a, :a, :a, :c] || train == [:a, :a, :a, :b, :b, :c]
    @test test == [:c, :a, :a, :a, :b, :b] || test == [:a, :a, :a, :c, :b, :b] || test == [:b, :b, :c, :a, :a, :a] || test == [:c, :b, :b, :a, :a, :a] || test == [:b, :b, :a, :a, :a, :c] || test == [:a, :a, :a, :b, :b, :c]
end

@testset "ObsView" begin
    ty = obsview([:a, :a, :a, :a, :a, :a, :b, :b, :b, :b])
    train, test = getobs.(@inferred(stratifiedobs(ty, 0.5, false)))
    @test train == [:a, :a, :a, :b, :b] || train == [:b, :b, :a, :a, :a]
    @test test == [:a, :a, :a, :b, :b] || test == [:b, :b, :a, :a, :a]

    train, test, val = getobs.(@inferred(stratifiedobs(ty, (0.5,0.3), false)))
    @test train == [:a, :a, :a, :b, :b] || train == [:b, :b, :a, :a, :a]
    @test test == [:a, :a, :b] || test == [:b, :a, :a]
    @test val == [:a, :b] || val == [:b, :a]
end

Random.seed!(42)
@testset "RNG" begin
    # tests reproducibility using explicit and global RNGs
    tX = [1:20;]
    ty = [:b, :b, :a, :b, :b, :b, :b, :a, :b, :a, :b, :b, :b, :a, :a, :b, :b, :b, :b, :b]
    (def_train_x, def_train_y), (def_test_x, def_test_y) = stratifiedobs((tX, ty))
    (exp_train_x, exp_train_y), (exp_test_x, exp_test_y) = stratifiedobs((tX, ty), rng=MersenneTwister(42))
    @test ((def_train_x, def_train_y), (def_test_x, def_test_y)) == ((exp_train_x, exp_train_y), (exp_test_x, exp_test_y))
    @test ((exp_train_x, exp_train_y), (exp_test_x, exp_test_y)) == stratifiedobs((tX, ty), rng=MersenneTwister(42))
end

@testset "RNG with callback" begin
    ty = [:b, :b, :a, :b, :b, :b, :b, :a, :b, :a, :b, :b, :b, :a, :a, :b, :b, :b, :b, :b]
    exp_train_y, exp_test_y = stratifiedobs(y->y==:a, ty, rng=MersenneTwister(42))
    @test (exp_train_y, exp_test_y) == stratifiedobs(y->y==:a, ty, rng=MersenneTwister(42))
end


text:
@testset "material" begin

@testset "tensorize" begin
    @test MaxwellBase.tensorize(3) === SMatrix{3,3,ComplexF64,9}(3,0,0, 0,3,0, 0,0,3)
    @test MaxwellBase.tensorize([2,4,6]) === SMatrix{3,3,ComplexF64,9}(2,0,0, 0,4,0, 0,0,6)
    @test MaxwellBase.tensorize([1 4 7; 2 5 8; 3 6 9]) === SMatrix{3,3,ComplexF64,9}(1,2,3,4,5,6,7,8,9)
end

@testset "Material" begin
    m_nn = Material("(num,num)", ε = 2, μ = 3)
    m_vv = Material("(vec,vec)", ε = [2.1,2.2,2.3], μ = [3.1,3.2,3.3])
    m_mm = Material("(mat,mat)", ε = [2 2 2; 2 2 2; 2 2 2], μ = [3 3 3; 3 3 3; 3 3 3])

    @test matparam(m_nn,EE) == diagm(0=>[2,2,2]) && matparam(m_nn,HH) == diagm(0=>[3,3,3]) && string(m_nn) == "(num,num)"
    @test matparam(m_vv,EE) == diagm(0=>[2.1,2.2,2.3]) && matparam(m_vv,HH) == diagm(0=>[3.1,3.2,3.3]) && string(m_vv) == "(vec,vec)"
    @test matparam(m_mm,EE) == fill(2, (3,3)) && matparam(m_mm,HH) == fill(3, (3,3)) && string(m_mm) == "(mat,mat)"
end

@testset "τ_trans" begin
    # Check if τ_trans for a 2×2 matrix is τ_trans for the 3×3 version of the matrix.  The
    # 3×3 version needs to have the original 2×2 matrix as the top-left 2×2 block, and the
    # third row and column has the (3,3) entry as the only nonzero entry.
    ε2 = @SMatrix rand(ComplexF64, 2, 2)
    ε3temp = zeros(ComplexF64, 3, 3)
    ε3temp[1:2,1:2] .= ε2
    ε3temp[3,3] = rand(ComplexF64)
    ε3 = MaxwellBase.S²ComplexF3(ε3temp)

    τ2 = MaxwellBase.τ_trans(ε2)
    τ3 = MaxwellBase.τ_trans(ε3)
    @test τ3[1:2,1:2] == τ2
end

end  # @testset "material"


text:
using Documenter, Playground

makedocs(
    format = :html,
    sitename = "Playground.jl",
    modules = [Playground],
    assets = ["assets/playground.css"],
    pages = [
        "index.md",
        "binaries.md",
        "executable.md",
        "repl.md",
        "api.md",
    ]
)

deploydocs(
    repo = "github.com/rofinn/Playground.jl.git",
    target = "build",
    julia  = "0.6",
    deps = nothing,
    make = nothing
)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using HDF5
using LightXML

mutable struct Xdmf <: AbstractResultsWriter
    name :: String
    xml :: XMLElement
    hdf :: HDF5File
    hdf_counter :: Int
    format :: String
end

function Xdmf()
    return Xdmf(tempname())
end

function h5file(xdmf::Xdmf)
    return xdmf.name*".h5"
end

function xmffile(xdmf::Xdmf)
    return xdmf.name*".xmf"
end

"""
    Xdmf(name, version="3.0", overwrite=false)

Initialize a new Xdmf object.
"""
function Xdmf(name::String; version="3.0", overwrite=false)
    xdmf = new_element("Xdmf")
    h5file = "$name.h5"
    xmlfile = "$name.xmf"

    if isfile(h5file)
        if overwrite
            @debug("Result file $h5file exists, removing old file.")
            rm(h5file)
        else
            error("Result file $h5file exists, use Xdmf($name; overwrite=true) to rewrite results")
        end
    end

    if isfile(xmlfile)
        if overwrite
            @debug("Result file $xmlfile exists, removing old file.")
            rm(xmlfile)
        else
            error("Result file $xmlfile exists, use Xdmf($name; overwrite=true) to rewrite results")
        end
    end

    set_attribute(xdmf, "xmlns:xi", "http://www.w3.org/2001/XInclude")
    set_attribute(xdmf, "Version", version)
    flag = isfile(h5file) ? "r+" : "w"
    hdf = h5open(h5file, flag)
    return Xdmf(name, xdmf, hdf, 1, "HDF")
end

"""
    get_temporal_collection(xdmf)

Return the basic structure of Xdmf document.
Creates a new TemporalCollection if not found.
Basic structure for XML part of Xdmf file is
<?xml version="1.0" encoding="utf-8"?>
<Xdmf xmlns:xi="http://www.w3.org/2001/XInclude" Version="2.1">
  <Domain>
    <Grid CollectionType="Temporal" GridType="Collection">
    </Grid>
  </Domain>
</Xdmf>
"""
function get_temporal_collection(xdmf::Xdmf)
    domain = find_element(xdmf.xml, "Domain")
    grid = nothing
    if domain == nothing
        domain = new_child(xdmf.xml, "Domain")
        grid = new_child(domain, "Grid")
        set_attribute(grid, "CollectionType", "Temporal")
        set_attribute(grid, "GridType", "Collection")
    end
    grid = find_element(domain, "Grid")
    return grid
end

"""
    xdmf_filter(child_elements, child_name)

Returns some spesific child xml element from an array of XMLElement based on,
"Xdmf extensions" see [1] for details.

Parameters
----------
child_elements :: Vector{XMLElement}
    A vector of XMLElements where to perform filtering.
child_name :: String
    Child element name, maybe containing Xdmf instructions

Returns
-------
nothing if nothing is found, otherwise XMLElement matching to filtering

#Examples

julia> grid1 = new_element("Grid")
julia> add_text(grid1, "I am first grid")
julia> grid2 = new_element("Grid")
julia> add_text(grid2, "I am second grid")
julia> set_attribute(grid2, "Name", "Frame 2")
julia> grid3 = new_element("Grid")
julia> add_text(grid3, "I am third grid")
julia> grids = [grid1, grid2, grid3]

To return second Grid element, one can use

julia> xdmf_filter(grids, "Grid[2]")

To return Grid which has attribute Name="Frame 2", use

julia> xdmf_filter(grids, "Grid[@name=Frame 2]")

To pick last Grid, use [end], e.g.

julia> xdmf_filter(grids, "Grid[end]").

References
----------
[1] http://www.xdmf.org/index.php/XDMF_Model_and_Format
"""
function xdmf_filter(child_elements, child_name)
    if '/' in child_name # needs path traversal
        return nothing
    end

    # filter children elements using syntax child[X] -> rename child_name
    m = match(r"(\w+)\[(.+)\]", child_name)
    if m != nothing
        child_name = m[1]
    end

    # first find any relevant child elements (has same tag)
    childs = []
    for child in child_elements
        if LightXML.name(child) == child_name
            push!(childs, child)
        end
    end

    # childs not found at all
    length(childs) == 0 && return nothing

    # by default return first
    m == nothing && return first(childs)

    # if [end] return last
    m[2] == "end" && return childs[end]

    # otherwise try parse int and return nth children from list
    parsed_int = tryparse(Int, m[2])
    if !isnull(parsed_int)
        idx = get(parsed_int)
        if (idx > 0) && (idx <= length(childs))
            return childs[idx]
        else
            # wrong index
            return nothing
        end
    end

    # [X] is something else than integer,  filter children elements using syntax child[@attr=value]
    m2 = match(r"@(.+)=(.+)", m[2])
    m2 == nothing && throw("Unable to parse: $(m[2])")
    attr_name = convert(String, m2[1])
    attr_value = convert(String, m2[2])
    for child in childs
        has_attribute(child, attr_name) || continue
        if attribute(child, attr_name) == attr_value
            return child
        end
    end

    # nothing found
    return nothing
end

"""
    traverse(xdmf, x, attr_name)

Traverse XML path. Xdmf filtering can be used, so it's possible to find
data from xml using syntax e.g.

#Example

julia> traverse(xdmf, x, "/Domain/Grid[2]/Grid[@Name=Frame 1]/DataItem")
"""
function traverse(xdmf::Xdmf, x::XMLElement, attr_name::String)
    attr_name = strip(attr_name, '/')

    if has_attribute(x, attr_name)
        return attribute(x, attr_name)
    end

    childs = child_elements(x)

    if '/' in attr_name
        items = split(attr_name, '/')
        new_item = xdmf_filter(childs, first(items))
        if new_item == nothing
            @debug("traverse: childs:")
            for child in childs
                @debug(LightXML.name(child))
            end
            error("traverse: failed, items = $items, xdmf_filter not find child")
        end
        new_path = join(items[2:end], '/')
        return traverse(xdmf, new_item, new_path)
    end

    child = xdmf_filter(childs, attr_name)
    return child
end

"""
    read(xdmf, path)

Read data from Xdmf file.

#Example

Traversing is supported, so one can easily traverse XML tree e.g.
julia> read(xdmf, "/Domain/Grid/Grid[2]/Geometry")
"""
function read(xdmf::Xdmf, path::String)
    result = traverse(xdmf, xdmf.xml, path)
    if endswith(path, "DataItem")
        format = attribute(result, "Format"; required=true)
        if format == "HDF"
            h5file, path = map(String, split(content(result), ':'))
            h5file = dirname(xdmf.name) * "/" * h5file
            isfile(h5file) || throw("Xdmf: h5 file $h5file not found!")
            return read(xdmf.hdf, path)
        else
            error("Read from Xdmf, reading from $format not implemented")
        end
    else
        return result
    end
end

"""
    save!(xdmf)

Save the xdmf file.
"""
function save!(xdmf::Xdmf)
    doc = XMLDocument()
    set_root(doc, xdmf.xml)
    save_file(doc, xmffile(xdmf))
end

function Base.close(xdmf::Xdmf)
    close(xdmf.hdf)
end

function new_dataitem(xdmf::Xdmf, path::String, data::Array{T,N}) where {T,N}
    dataitem = new_element("DataItem")
    datatype = replace("$T", "64" => "")
    dimensions = join(reverse(size(data)), " ")
    set_attribute(dataitem, "DataType", datatype)
    set_attribute(dataitem, "Dimensions", dimensions)
    set_attribute(dataitem, "Format", xdmf.format)
    if xdmf.format == "HDF"
        hdf = basename(h5file(xdmf))
        if exists(xdmf.hdf, path)
            @debug("Xdmf: $path already existing in h5 file, not overwriting.")
        else
            write(xdmf.hdf, path, data)
        end
        add_text(dataitem, "$hdf:$path")
    elseif xdmf.format == "XML"
        text_data = string(data')
        text_data = strip(text_data, ['[', ']'])
        text_data = replace(text_data, ';', '\n')
        text_data = "\n" * text_data * "\n"
        add_text(dataitem, text_data)
    else
        error("Unsupported Xdmf big data format $(xdmf.format)")
    end
    return dataitem
end

"""
    new_dataitem(xdmf, data)

Create a new DataItem element, hdf path automatically determined.
"""
function new_dataitem(xdmf::Xdmf, data::Array{T,N}) where {T,N}
    if xdmf.format == "XML"
        # Path can be whatever as XML format does not store to HDF at all
        return new_dataitem(xdmf, "/whatever", data)
    else
        path = "/DataItem_$(xdmf.hdf_counter)"
        while exists(xdmf.hdf, path)
            xdmf.hdf_counter += 1
            path = "/DataItem_$(xdmf.hdf_counter)"
        end
        return new_dataitem(xdmf, path, data)
    end
end

global const xdmf_element_mapping = Dict(
        "Poi1" => "Polyvertex",
        "Seg2" => "Polyline",
        "Tri3" => "Triangle",
        "Quad4" => "Quadrilateral",
        "Tet4" => "Tetrahedron",
        "Pyramid5" => "Pyramid",
        "Wedge6" => "Wedge",
        "Hex8" => "Hexahedron",
        "Seg3" => "Edge_3",
        "Tri6" => "Tri_6",
        "Quad8" => "Quad_8",
        "Tet10" => "Tet_10",
        "Pyramid13" => "Pyramid_13",
        "Wedge15" => "Wedge_15",
        "Hex20" => "Hex_20")

get_xdmf_element_code(::Element{Poi1}) = 1
get_xdmf_element_code(::Element{Seg2}) = 2
# get_xdmf_element_code(::Element{Polygon}) =  3
get_xdmf_element_code(::Element{Tri3}) = 4
get_xdmf_element_code(::Element{Quad4}) = 5
get_xdmf_element_code(::Element{Tet4}) = 6
get_xdmf_element_code(::Element{Pyr5}) = 7
get_xdmf_element_code(::Element{Wedge6}) = 8
get_xdmf_element_code(::Element{Hex8}) = 9
# get_xdmf_element_code(::Element{Polyhedron}) = 16

get_xdmf_element_code(::Element{Seg3}) = 34
get_xdmf_element_code(::Element{Quad9}) = 35
get_xdmf_element_code(::Element{Tri6}) = 36
get_xdmf_element_code(::Element{Quad8}) = 37
get_xdmf_element_code(::Element{Tet10}) = 38
# get_xdmf_element_code(::Element{Pyr13}) = 39
get_xdmf_element_code(::Element{Wedge15}) = 40
# get_xdmf_element_code(::Element{Wedge18}) = 41
get_xdmf_element_code(::Element{Hex20}) = 48
# get_xdmf_element_code(::Element{Hex24}) = 49
get_xdmf_element_code(::Element{Hex27}) = 50

"""
    get_spatial_collection()

Return a SpatialCollection at given time either by creating new one or returning
existing one.
"""
function get_spatial_collection(temporal_collection, time)
    for spatial_collection in get_elements_by_tagname(temporal_collection, "Grid")
        time_element = find_element(spatial_collection, "Time")
        time_value = Meta.parse(attribute(time_element, "Value"; required=true))
        isapprox(time_value, time) && return spatial_collection
    end
    # did not find, create new one
    spatial_collection = new_child(temporal_collection, "Grid")
    set_attribute(spatial_collection, "GridType", "Collection")
    set_attribute(spatial_collection, "Name", "Problems")
    set_attribute(spatial_collection, "CollectionType", "Spatial")
    time_element = new_child(spatial_collection, "Time")
    set_attribute(time_element, "Value", time)
    return spatial_collection
end

"""
    update_xdmf!(xdmf, problem, time, fields)

Write new fields to Xdmf file.

#Example

To write displacement and temperature fields from p1 at time t=0.0:

julia> update_xdmf!(p1, 0.0, ["displacement", "temperature"])
"""
function update_xdmf!(xdmf::Xdmf, problem::Problem, time::Float64, fields::Vector)

    @debug("Xdmf: storing fields $fields of problem $(problem.name) at time $time")

    # 1. find domain
    xml = xdmf.xml
    domain = find_element(xml, "Domain")
    if domain == nothing
        @debug("Xdmf: Domain not found, creating.")
        domain = new_child(xml, "Domain")
    end

    # 2. find for TemporalCollection
    temporal_collection = find_element(domain, "Grid")
    if temporal_collection == nothing
        @debug("Xdmf: Temporal collection not found, creating.")
        temporal_collection = new_child(domain, "Grid")
        set_attribute(temporal_collection, "GridType", "Collection")
        set_attribute(temporal_collection, "Name", "Time")
        set_attribute(temporal_collection, "CollectionType", "Temporal")
    end

    # 2.1 make sure that Grid element we found really is TemporalCollection
    collection_type = attribute(temporal_collection, "CollectionType"; required=true)
    @assert collection_type == "Temporal"

    spatial_collection = get_spatial_collection(temporal_collection, time)

    for frame in get_elements_by_tagname(spatial_collection, "Grid")
        frame_name = attribute(frame, "Name")
        if frame_name == problem.name
            @warn("Xdmf: Already found Grid with name $frame_name for time $time, skipping.")
            return
        end
    end

    frame_name = problem.name
    @debug("Xdmf: Creating Grid for problem $frame_name")
    frame = new_child(spatial_collection, "Grid")
    set_attribute(frame, "Name", frame_name)

    # 4. save geometry
    X_dict = problem("geometry", time)
    node_ids = sort(collect(keys(X_dict)))
    node_mapping = Dict(j => i for (i, j) in enumerate(node_ids))
    X_array = hcat([X_dict[nid] for nid in node_ids]...)
    ndim, nnodes = size(X_array)
    geom_type = (ndim == 2 ? "XY" : "XYZ")
    @debug("Xdmf: Creating geometry, type = $geom_type, number of nodes = $nnodes")
    X_dataitem = new_dataitem(xdmf, X_array)
    geometry = new_child(frame, "Geometry")
    set_attribute(geometry, "Type", geom_type)
    add_child(geometry, X_dataitem)

    # 5. save topology
    mesh_type = "unstructured"
    if mesh_type == "unstructured"
        element_conn = Int64[]
        for element in get_elements(problem)
            xdmf_element_code = get_xdmf_element_code(element)
            xdmf_element_code > 0 || continue
            push!(element_conn, xdmf_element_code)
            if xdmf_element_code == 2
                push!(element_conn, length(element))
            end
            for j in get_connectivity(element)
                push!(element_conn, node_mapping[j]-1)
            end
        end
        topology_dataitem = new_dataitem(xdmf, element_conn)
        topology = new_child(frame, "Topology")
        set_attribute(topology, "TopologyType", "Mixed")
        add_child(topology, topology_dataitem)
    else
        all_elements = get_elements(problem)
        nelements = length(all_elements)
        element_types = unique(map(get_element_type, all_elements))
        nelement_types = length(element_types)
        @debug("Xdmf: Saving topology of $nelements elements total, $nelement_types different element types.")
        if nelement_types != 1
            error("Xdmf: only single type of element supported by structured grid type!")
        end
        for element_type in element_types
            elements = collect(filter_by_element_type(element_type, all_elements))
            nelements = length(elements)
            @debug("Xdmf: $nelements elements of type $element_type")
            sort!(elements, by=get_element_id)
            element_ids = map(get_element_id, elements)
            element_conn = map(element -> [node_mapping[j]-1 for j in get_connectivity(element)], elements)
            element_conn = hcat(element_conn...)
            element_code = split(string(element_type), ".")[end]
            topology_dataitem = new_dataitem(xdmf, element_conn)
            topology = new_child(frame, "Topology")
            set_attribute(topology, "TopologyType", xdmf_element_mapping[element_code])
            set_attribute(topology, "NumberOfElements", length(elements))
            add_child(topology, topology_dataitem)
        end
    end

    # 6. save requested fields
    for field_name in fields
        field_dict = problem(field_name, time)
        field_center = "Node"
        field_node_ids = sort(collect(keys(field_dict)))
        if node_ids != field_node_ids
            @error("geom node ids = $node_ids")
            @error("field node ids = $field_node_ids")
            error("!=, geometry does not match with field.")
        end
        field_dim = length(field_dict[first(field_node_ids)])
        if field_dim == 2
            @debug("Xdmf: Field dimension = 2, extending to 3")
            for nid in field_node_ids
                field_dict[nid] = [field_dict[nid]; 0.0]
            end
            field_dim = 3
        end
        field_type = Dict(1 => "Scalar", 3 => "Vector", 6 => "Tensor6")[field_dim]
        @debug("Xdmf: Saving field $field_name, type = $field_type, dimension = $field_dim, center = $field_center")

        field_array = hcat([field_dict[nid] for nid in field_node_ids]...)
        field_dataitem = new_dataitem(xdmf, field_array)
        attribute = new_child(frame, "Attribute")
        set_attribute(attribute, "Name", uppercasefirst(field_name))
        set_attribute(attribute, "Center", field_center)
        set_attribute(attribute, "AttributeType", field_type)
        add_child(attribute, field_dataitem)
    end

    save!(xdmf)
    @debug("Xdmf: all done.")
end


text:
function GLBUpdater(Tree)

    GLB = Tree[1].LB;
    for i = 2:length(Tree)
        if Tree[i].LB < GLB
            GLB = Tree[i].LB
        end
    end
    return GLB
end


text:
using NPZ
using Dierckx: Spline1D, derivative
using Interpolations

# -------------------------------- TRACK TYPE -------------------------------- #
struct Track
    X::Vector
    Y::Vector
    XY::Matrix
    curvature::Vector   # κ value at each waypoint
    N::Int
    P::Vector  # parameter values
    S_f::Float64
    S::Vector
    δs::Vector
    κ::Any
    width::Function
    θ::Vector
end

function Base.show(io::IO, track::Track)
    return print(
        io,
        "Track. \e[2m$(round(track.S_f)) cm long, $(track.N) waypoints.\e[0m",
    )
end

struct Border
    X::Vector
    Y::Vector
end

# ----------------------------------- utils ---------------------------------- #

"""
Get border lines by extruding along normal directoin
"""
function get_track_borders(track::Track)
    η = track.θ .+ π / 2
    dx = track.width.(track.S) .* cos.(η)
    dy = track.width.(track.S) .* sin.(η)

    left = Border(track.X .+ dx, track.Y .+ dy)
    right = Border(track.X .- dx, track.Y .- dy)

    return left, right
end

"""
Compute the curvature at each waypoint in the track
"""
function compute_curvature(t, x, y, tnew)
    xspl = Spline1D(t, x)
    yspl = Spline1D(t, y)
    dx = derivative(xspl, tnew; nu=1)
    dy = derivative(yspl, tnew; nu=1)
    ddx = derivative(xspl, tnew; nu=2)
    ddy = derivative(yspl, tnew; nu=2)
    return -((ddx .* dy) .- (ddy .* dx)) ./ (((dx) .^ 2 + (dy) .^ 2) .^ (3 / 2))
end

"""
  κ(s)
Curvature as a function of curvilinear distance along the track
"""
function κ(s, k, S_f)
    idx = (Int ∘ ceil)(s / S_f * length(k))
    idx = idx == 0 ? 1 : idx
    idx = idx >= length(k) ? idx - 1 : idx

    if idx > length(k)
        @warn "how" s k S_f idx
    end
    return k[idx]
end

"""
Value of the width factor of the track at various svalues
"""
width_values = [
    [0.0 1.2]
    [0.01 1.1]
    [0.05 1.00]
    [0.07 .9] # first narrow
    [0.11 1.00] # first narrow
    [0.15 1]  # end of frst narrow
    [0.20 1]
    [0.25 1.00]
    [0.30 1.1]   # middle of second curve
    [0.36 1.2]
    [0.45 1.1]  # second narrow
    [0.5 0.9]  #
    [0.55 1.05]
    [0.6 1.0]  # end of second narrow
    [0.63 1.05]
    [0.67 1.05]
    [0.7 1.1]
    [0.75 1.1]
    [0.80 1.1]  # second part of last curve
    [0.9 1.1]
    [0.92 1.1]
    [0.95 1.15]
    [0.98 1.2]
    [1 1.2]  # end
]

# ---------------------------------------------------------------------------- #
#                              TRACK CONSTRUCTORS                              #
# ---------------------------------------------------------------------------- #

"""
Construct `Track` out of a set of waypoints
"""
function Track(XY, s1::Float64; resolution=0.00001, const_width=false)
    X = movingaverage(XY[:, 1], 3)
    Y = movingaverage(XY[:, 2], 3)
    X, Y = upsample(X, Y; δp=resolution)

    # @assert length(X) == size(XY, 1)
    # get new points locations thorugh interpolation
    # X, Y = upsample(XY[:, 1], XY[:, 2]; δp=resolution)
    N = length(X)

    # get distance step between each track point + total length
    Δx, Δy = Δ(X), Δ(Y)
    δs = sqrt.(Δx .^ 2 .+ Δy .^ 2)
    S_f = sum(δs)

    S = cumsum(δs)  # curvilinear coordinates
    @assert S[1] == 0 && S[end] - S_f < 0.01 "ops. $(S[1]) $(S[end])==$(S_f)"

    # compute curvature
    P = range(0, 1; length=length(X))
    curvature = compute_curvature(P, X, Y, P)
    K(s) = κ(s, curvature, S_f)
    @assert length(curvature) == length(X)

    # compute orientation of each segment
    θ = unwrap(atan.(Δy, Δx))
    θ[1] = θ[2]

    # get width function working for short tracks
    wspline = Spline1D(width_values[:, 1], width_values[:, 2] .* 3; k=4)
    wfn(s) = const_width ? 3.0 : wspline((s - s1) / 261)

    return Track(X, Y, Array([X Y]'), curvature, N, P, S_f, S, δs, K, wfn, θ)
end

"""
Create a track from saved waypoints coordinates.

#Arguments
  - `start_waypoint` index of first wp to keep
  - `keep_n_waypoints` can be used to keep only the first N waypoints. Set to -1 to keep all.
  - `resolution` used to upsample track waypoints through interpolation.
"""
function Track(;
    start_waypoint=2,
    keep_n_waypoints=-1,
    resolution=0.00001,
    npyfile=nothing,
    const_width=false,
    track_length=261,
    )
    # load data
    npyfile = isnothing(npyfile) ? "src/hairpin.npy" : npyfile
    XY = npzread(npyfile)

    # trim waypoints
    start_waypoint = max(2, start_waypoint)
    keep_n_waypoints = if keep_n_waypoints > 0
        min(keep_n_waypoints, size(XY, 1))
    else
        (size(XY, 1) - start_waypoint)
    end

    if isnothing(track_length)
        track_length = (Int64 ∘ round)(sum(sqrt.(diff(XY[:, 1]).^2 .+ diff(XY[:, 2]).^2))) * 2
    end

    end_idx = min(start_waypoint + keep_n_waypoints - 1, track_length)
    XY = XY[start_waypoint:end_idx, :]

    s1 = (start_waypoint) / track_length

    return Track(XY, s1; resolution=resolution, const_width=const_width)
end

"""
Create a `Track` with the initial waypoint closest to the position of
the bike (defined by its State).
"""
function Track(state; keep_n_waypoints=-1, resolution=0.00001)
    XY = npzread("src/hairpin.npy")
    idx = closest_point_idx(XY[:, 1], state.x, XY[:, 2], state.y)

    return Track(;
        start_waypoint=idx, keep_n_waypoints=keep_n_waypoints, resolution=resolution
    )
end

const FULLTRACK = Track(; start_waypoint=4)

"""
Trim the full track from a start value keeping a given length
"""
function trim(track::Track, svalue, length)
    # svalue = svalue < 1 ? svalue * 259 : svalue

    first = findfirst(track.S .>= (svalue))

    last = findlast(track.S .<= (svalue + length))

    return Track(
        FULLTRACK.X[first:last],
        FULLTRACK.Y[first:last],
        FULLTRACK.XY[:, first:last],
        FULLTRACK.curvature[first:last],
        -1,
        FULLTRACK.P[first:last],
        FULLTRACK.S[last],
        FULLTRACK.S[first:last],
        FULLTRACK.δs[first:last],
        FULLTRACK.κ,
        FULLTRACK.width,
        FULLTRACK.θ[first:last],
    )
end

max_stars_repo_name: FedeClaudi/LocomotionControl

max_issues_repo_name: FedeClaudi/LocomotionControl

max_forks_repo_name: FedeClaudi/LocomotionControl


text:
"""`immutable ComplexSpace <: EuclideanSpace{ℂ}`

A ComplexSpace is a standard complex vector space ℂ^d with Euclidean inner product
and no additional structure. It is completely characterised by its dimension and
whether its the normal space or its dual (which is canonically isomorphic to the
conjugate space)."""
struct ComplexSpace <: EuclideanSpace{ℂ}
  d::Int
  dual::Bool
end
ComplexSpace(d::Int) = ComplexSpace(d, false)

# convenience constructor
Base.:^(::ComplexNumbers, d::Int) = ComplexSpace(d)
Base.getindex(::ComplexNumbers) = ComplexSpace
Base.getindex(::ComplexNumbers, d::Int) = ComplexSpace(d)

# Corresponding methods:
#------------------------
dim(V::ComplexSpace) = V.d
Base.indices(V::ComplexSpace) = Base.OneTo(dim(V))

Base.conj(V::ComplexSpace) = ComplexSpace(V.d, !V.dual)

Base.show(io::IO, V::ComplexSpace) = print(io, V.dual ? "(ℂ^$(V.d))'" : "ℂ^$(V.d)")

# direct sum
⊕(V1::ComplexSpace, V2::ComplexSpace) = (V1.dual==V2.dual ?
    ComplexSpace(V1.d+V2.d, V1.dual) :
    throw(SpaceMismatch("Direct sum of a vector space and its dual do not exist")))


text:
using DependencyWalker
using Test
using ObjectFile, Pango_jll

@testset "DependencyWalker.jl" begin
    pango = Library(Pango_jll.libpango_path)
    @show pango
    @test pango isa Library{<:ObjectHandle}
    @test Library("this does not exist.foo") isa Library{Missing}
end


text:
# This file was generated, do not modify it. # hide
db = MatrixDatabase(randn(2, 10^5))
dist = L2Distance() # squared L2

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
"""
    loadFCSHeader(fn::String)::Tuple{Vector{Int}, Dict{String,String}}

Efficiently extract data offsets and keyword dictionary from an FCS file.
"""
function loadFCSHeader(fn::String)::Tuple{Vector{Int},Dict{String,String}}
    open(fn) do io
        offsets = FCSFiles.parse_header(io)
        params = FCSFiles.parse_text(io, offsets[1], offsets[2])
        FCSFiles.verify_text(params)
        (offsets, params)
    end
end

"""
    getFCSSize(offsets, params)::Tuple{Int,Int}

Convert the offsets and keywords from an FCS file to cell and parameter count,
respectively.
"""
function getFCSSize(offsets, params)::Tuple{Int,Int}
    nData = parse(Int, params["\$TOT"])
    nParams = parse(Int, params["\$PAR"])

    if params["\$DATATYPE"] != "F"
        @error "Only float32 FCS files are currently supported"
        error("Unsupported FCS format")
    end

    beginData = parse(Int, params["\$BEGINDATA"])
    endData = parse(Int, params["\$ENDDATA"])

    #check that the $TOT and $PAR look okay
    if !(offsets[3] == 0 && offsets[4] == 0) && (
        (
            1 + offsets[4] - offsets[3] != nData * nParams * 4 &&
            offsets[4] - offsets[3] != nData * nParams * 4
        ) ||
        offsets[3] != beginData ||
        offsets[4] != endData
    )
        @warn "Data size mismatch, FCS is likely broken."
    end

    return (nData, nParams)
end

"""
    loadFCSSizes(fns::Vector{String})

Load cell counts in many FCS files at once. Useful as input for `slicesof`.
"""
function loadFCSSizes(fns::Vector{String})::Vector{Int}
    [(
        begin
            o, s = loadFCSHeader(fn)
            getFCSSize(o, s)[1]
        end
    ) for fn in fns]
end

"""
    loadFCS(fn::String; applyCompensation::Bool=true)::Tuple{Dict{String,String}, Matrix{Float64}}

Read a FCS file. Return a tuple that contains in order:

- dictionary of the keywords contained in the file
- raw column names
- prettified and annotated column names
- raw data matrix

If `applyCompensation` is set, the function parses and retrieves a spillover
matrix (if any valid keyword in the FCS is found that would contain it) and
applies it to compensate the data.
"""
function loadFCS(
    fn::String;
    applyCompensation::Bool = true,
)::Tuple{Dict{String,String},Matrix{Float64}}
    fcs = FileIO.load(fn)
    meta = getMetaData(fcs.params)
    data = hcat(map(x -> Vector{Float64}(fcs.data[x]), meta[:, :N])...)
    if applyCompensation
        spill = getSpillover(fcs.params)
        if spill != nothing
            names, mtx = spill
            cols = indexin(names, meta[:, :N])
            if any(cols .== nothing)
                @error "Unknown columns in compensation matrix" names cols
                error("Invalid compensation matrix")
            end
            compensate!(data, mtx, Vector{Int}(cols))
        end
    end
    return (fcs.params, data)
end

"""
    loadFCSSet(name::Symbol, fns::Vector{String}, pids=workers(); applyCompensation=true, postLoad=(d,i)->d)::LoadedDataInfo

This runs the FCS loading machinery in a distributed way, so that the files
`fns` (with full path) are sliced into equal parts and saved as a distributed
variable `name` on workers specified by `pids`.

`applyCompensation` is passed to loadFCS function.

See `slicesof` for description of the slicing.

`postLoad` is applied to the loaded FCS file data (and the index) -- use this
function to e.g. filter out certain columns right on loading, using `selectFCSColumns`.

The loaded dataset can be manipulated by the distributed functions, e.g.
- `dselect` for removing columns
- `dscale` for normalization
- `dtransform_asinh` (and others) for transformation
- etc.
"""
function loadFCSSet(
    name::Symbol,
    fns::Vector{String},
    pids = workers();
    applyCompensation = true,
    postLoad = (d, i) -> d,
)::LoadedDataInfo
    slices = slicesof(loadFCSSizes(fns), length(pids))
    distributed_foreach(
        slices,
        (slice) -> Base.eval(
            Main,
            :(
                begin
                    $name = vcollectSlice(
                        (i) -> last($postLoad(
                            loadFCS($fns[i]; applyCompensation = $applyCompensation),
                            i,
                        )),
                        $slice,
                    )
                    nothing
                end
            ),
        ),
        pids,
    )
    return LoadedDataInfo(name, pids)
end

"""
    selectFCSColumns(selectColnames::Vector{String})

Return a function useful with `loadFCSSet`, which loads only the specified
(prettified) column names from the FCS files. Use `getMetaData`,
`getMarkerNames` and `cleanNames!` to retrieve the usable column names for a
FCS.
"""
function selectFCSColumns(selectColnames::Vector{String})
    ((metadata, data), idx) -> begin
        _, names = getMarkerNames(getMetaData(metadata))
        cleanNames!(names)
        colIdxs = indexin(selectColnames, names)
        if any(colIdxs .== nothing)
            @error "Some columns were not found"
            error("unknown column")
        end
        (metadata, data[:, colIdxs])
    end
end

"""
    distributeFCSFileVector(name::Symbol, fns::Vector{String}, pids=workers())::LoadedDataInfo

Distribute a vector of integers among the workers that describes which file
from `fns` the cell comes from. Useful for producing per-file statistics. The
vector is saved on workers specified by `pids` as a distributed variable
`name`.
"""
function distributeFCSFileVector(
    name::Symbol,
    fns::Vector{String},
    pids = workers(),
)::LoadedDataInfo
    sizes = loadFCSSizes(fns)
    slices = slicesof(sizes, length(pids))
    return distributeFileVector(name, sizes, slices, pids)
end

"""
    distributeFileVector(name::Symbol, sizes::Vector{Int}, slices::Vector{Tuple{Int,Int,Int,Int}}, pids=workers())::LoadedDataInfo

Generalized version of `distributeFCSFileVector` that produces the integer
vector from any `sizes` and `slices`.
"""
function distributeFileVector(
    name::Symbol,
    sizes::Vector{Int},
    slices::Vector{Tuple{Int,Int,Int,Int}},
    pids = workers(),
)::LoadedDataInfo
    distributed_foreach(
        slices,
        (slice) ->
            Base.eval(Main, :($name = collectSlice((i) -> fill(i, $sizes[i]), $slice))),
        pids,
    )
    return LoadedDataInfo(name, pids)
end

"""
    function getCSVSize(fn::String; args...)::Tuple{Int,Int}

Read the dimensions (number of rows and columns, respectively) from a CSV file
`fn`. `args` are passed to function `CSV.file`.

# Example

    getCSVSize("test.csv", header=false)
"""
function getCSVSize(fn::String; args...)::Tuple{Int,Int}
    n = 0
    k = 0
    # ideally, this will not try to load the whole CSV in the memory
    for row in CSV.File(fn, type = Float64; args...)
        n += 1
        if length(row) > k
            k = length(row)
        end
    end
    return (n, k)
end

"""
    function loadCSVSizes(fns::Vector{String}; args...)::Vector{Int}

Determine number of rows in a list of CSV files (passed as `fns`). Equivalent
to `loadFCSSizes`.
"""
function loadCSVSizes(fns::Vector{String}; args...)::Vector{Int}
    [getCSVSize(fn, type = Float64; args...)[1] for fn in fns]
end

"""
    function loadCSV(fn::String; args...)::Matrix{Float64}

CSV equivalent of `loadFCS`. The metadata (header, column names) are not
extracted. `args` are passed to `CSV.read`.
"""
function loadCSV(fn::String; args...)::Matrix{Float64}
    CSV.read(fn, DataFrame, type = Float64; args...) |> Matrix{Float64}
end

"""
    function loadCSVSet(
        name::Symbol,
        fns::Vector{String},
        pids = workers();
        postLoad = (d, i) -> d,
        csvargs...,
    )::LoadedDataInfo

CSV equivalent of `loadFCSSet`. `csvargs` are passed as keyword arguments to
CSV-loading functions.
"""
function loadCSVSet(
    name::Symbol,
    fns::Vector{String},
    pids = workers();
    postLoad = (d, i) -> d,
    csvargs...,
)::LoadedDataInfo
    slices = slicesof(loadCSVSizes(fns; csvargs...), length(pids))
    distributed_foreach(
        slices,
        (slice) -> Base.eval(
            Main,
            :(
                begin
                    $name = vcollectSlice(
                        (i) -> $postLoad(loadCSV($fns[i]; $csvargs...), i),
                        $slice,
                    )
                    nothing
                end
            ),
        ),
        pids,
    )
    return LoadedDataInfo(name, pids)
end


text:
#JuliaRL/src/BasicDQN.jl

using ReinforcementLearning
using PyCall
using ReinforcementLearningEnvironments
using Random
using Flux
using TensorBoardLogger
using Dates
using Logging
using BSON
using Suppressor
include("./conf.jl")
include("./shared.jl")

"""
    runBasicDQN(save_dir, landing_factor=1, leg_first_bonus=0)

Trains the modified LunarLander agent using BasicDQN
from the ReinforcementLearning library.
Results will be saved at `save_dir`.

Using the default values for landing_factor and leg_first_bonus
will make the environement behave like the original LunarLander.
For more information, see the JuliaRL/CustomGym/CustomGym/lunar_lander.py file.
"""
function runBasicDQN(save_dir::T, landing_factor=1, leg_first_bonus=0) where {T<:AbstractString}
    # clear save_dir directory
    isdir(save_dir) && rm(save_dir; force=true, recursive=true)

    lg = TBLogger(joinpath(save_dir, "tb_log"), min_level = Logging.Info)
    rng = MersenneTwister(123)

    env = LunarLander(landing_factor, leg_first_bonus)

    ns, na = length(get_state(env)), length(get_actions(env))

    agent = Agent(
        policy = QBasedPolicy(
            learner = BasicDQNLearner(
                approximator = NeuralNetworkApproximator(
                    model = net_model(ns, na, rng),
                    optimizer = ADAM(),
                ),
                batch_size = Conf.batch_size,
                min_replay_history = Conf.min_replay_history,
                loss_func = huber_loss,
                rng = rng,
            ),
            explorer = EpsilonGreedyExplorer(
                kind = :exp,
                ϵ_stable = 0.01,
                decay_steps = Conf.decay_steps,
                rng = rng,
            ),
        ),
        trajectory = CircularCompactSARTSATrajectory(
            capacity = Conf.capacity,
            state_type = Float32,
            state_size = (ns,),
        ),
    )

    stop_condition = StopAfterStep(Conf.duration)

    global episode_loss = 0
    total_reward_per_episode = TotalRewardPerEpisode()
    time_per_step = TimePerStep()
    hook = ComposedHook(
        total_reward_per_episode,
        time_per_step,
        DoEveryNStep() do t, agent, env
            global episode_loss += loss = agent.policy.learner.loss
        end,
        DoEveryNEpisode() do t, agent, env
            with_logger(lg) do
                global episode_loss
                @info "training" loss = episode_loss
                @info "training" reward = total_reward_per_episode.rewards[end]
                log_step_increment = 0
                episode_loss = 0
            end
        end,
        DoEveryNStep(Conf.save_freq) do t, agent, env
            RLCore.save(save_dir, agent)
            BSON.@save joinpath(save_dir, "stats.bson") total_reward_per_episode time_per_step
        end,
    )

    infos = @timeRet run(agent, env, stop_condition, hook)
    log_training_info(infos, agent, save_dir)
end

max_stars_repo_name: FeLusiani/JuliaRL

max_issues_repo_name: FeLusiani/JuliaRL

max_forks_repo_name: FeLusiani/JuliaRL


text:
# ---
# title: Incremental vs. Batch Example
# id: incremental_batch
# cover: ../assets/art.png
# date: 2021-12-1
# author: "[Sasha Petrenko](https://github.com/AP6YC)"
# julia: 1.6
# description: This demo illustrates how to use incremental training methods vs. batch training for all ART modules.
# ---

# ## Overview

# All modules in `AdaptiveResonance.jl` are designed to handle incremental and batch training.
# In fact, ART modules are generally incremental in their implementation, so their batch methods wrap the incremental ones and handle preprocessing, etc.
# For example, DDVFA can be run incrementally (i.e. with one sample at a time) with custom algorithmic options and a predetermined data configuration.

# !!! note
#     In the incremental case, it is necessary to provide a data configuration if the model is not pretrained because the model has no knowledge of the boundaries and dimensionality of the data, which are necessary in the complement coding step.
#     For more info, see the guide in the docs on [incremental vs. batch](@ref incremental_vs_batch).

# ## Data Setup

# We begin with importing AdaptiveResonance for the ART modules and MLDatasets for some data utilities.
using AdaptiveResonance # ART
using MLDatasets        # Iris dataset
using MLDataUtils       # Shuffling and splitting
using Printf            # Formatted number printing

# We will download the Iris dataset for its small size and benchmark use for clustering algorithms.
Iris.download(i_accept_the_terms_of_use=true)
features, labels = Iris.features(), Iris.labels()

# Because the MLDatasets package gives us Iris labels as strings, we will use the `MLDataUtils.convertlabel` method with the `MLLabelUtils.LabelEnc.Indices` type to get a list of integers representing each class:
labels = convertlabel(LabelEnc.Indices{Int}, labels)
unique(labels)

# Next, we will create a train/test split with the `MLDataUtils.stratifiedobs` utility:
(X_train, y_train), (X_test, y_test) = stratifiedobs((features, labels))

# ## Incremental vs. Batch

# ### Setup

# Now, we can create several modules to illustrate training one in batch and one incrementaly.

## Create several modules for batch and incremental training.
## We can take advantage of the options instantiation method here to use the same options for both modules.
opts = opts_DDVFA(rho_lb=0.6, rho_ub=0.75)
art_batch = DDVFA(opts)
art_incremental = DDVFA(opts)

# For the incremental version, we must setup the data configuration in advance.
# In batch mode, this is done automatically based upon the provided data, but the incremental variant has not way of knowing the bounds of the individual features.
# We *could* preprocess the data and set the data configuration with `art.config = DataConfig(0, 1, 4)`, which translates to the data containing four features  that *all* range from 0 to 1.
# This would be done in scenarios where we have either done some preprocessing on the data or have prior knowledge about the bounds of individual features.
# However, in this example we will let the module determine the bounds with the convenience method `data_setup!`:

## Setup the data config on all of the features.
data_setup!(art_incremental.config, features)

# ### Training

# We can train in batch with a simple supervised mode by passing the labels as a keyword argument.
y_hat_batch_train = train!(art_batch, X_train, y=y_train)
println("Training labels: ",  size(y_hat_batch_train), " ", typeof(y_hat_batch_train))

# We can also train incrementally with the same method, being careful that we pass a vector features and a single integer as the labels

## Get the number of training samples
n_train = length(y_train)
## Create a container for the training output labels
y_hat_incremental_train = zeros(Int, n_train)
## Iterate over all training samples
for ix = 1:length(y_train)
    sample = X_train[:, ix]
    label = y_train[ix]
    y_hat_incremental_train[ix] = train!(art_incremental, sample, y=label)
end

# ### Testing

# We can then classify both networks and check that their performances are equivalent.
# For both, we will use the best-matching unit in the case of complete mismatch (see the docs on [Mismatch vs. BMU](@ref mismatch-bmu))

## Classify one model in batch mode
y_hat_batch = AdaptiveResonance.classify(art_batch, X_test, get_bmu=true)

## Classify one model incrementally
n_test = length(y_test)
y_hat_incremental = zeros(Int, n_test)
for ix = 1:n_test
    y_hat_incremental[ix] = AdaptiveResonance.classify(art_incremental, X_test[:, ix], get_bmu=true)
end

## Check the shape and type of the output labels
println("Batch testing labels: ",  size(y_hat_batch), " ", typeof(y_hat_batch))
println("Incremental testing labels: ",  size(y_hat_incremental), " ", typeof(y_hat_incremental))

# Finally, we check the performance (number of correct classifications over total number of test samples) for both models, verifying that they produce the same results.

## Calculate performance on training data, testing data, and with get_bmu
perf_train_batch = performance(y_hat_batch_train, y_train)
perf_train_incremental = performance(y_hat_incremental_train, y_train)
perf_test_batch = performance(y_hat_batch, y_test)
perf_test_incremental = performance(y_hat_incremental, y_test)

## Format each performance number for comparison
@printf "Batch training performance: %.4f\n" perf_train_batch
@printf "Incremental training performance: %.4f\n" perf_train_incremental
@printf "Batch testing performance: %.4f\n" perf_test_batch
@printf "Incremental testing performance: %.4f\n" perf_test_incremental


text:
@testset "moments" begin

img = testimage("cameraman")
edges, counts = build_histogram(reinterpret.(gray.(img)), 0:1:255)
t1 = find_threshold(Moments(),counts[1:end], 0:1:255)
@test t1 == 111

# This replicates a test described in the original publication.
arr = [10  8 10  9 20 21 32 30 40 41 41 40;
       12 10 11 10 19 20 30 28 38 40 40 39;
       10  9 10  8 20 21 30 29 42 40 40 39;
       11 10  9 11 19 21 31 30 40 42 38 40]
edges, counts = build_histogram(arr, 0:255)
t2 = find_threshold(Moments(), counts[1:end], edges)
@test t2 == 27

# This test reaches into the internals of the implementation to verify that
# the moments are indeed preserved.
m₀ = 1.0
m₁, m₂, m₃ =  HistogramThresholding.compute_moments(counts[1:end], edges, sum(counts[1:end]))
c₀, c₁ = HistogramThresholding.compute_auxiliary_values(m₀, m₁, m₂, m₃)
z₀, z₁ = HistogramThresholding.compute_gray_levels(c₀, c₁)
p₀, p₁ = HistogramThresholding.compute_pixel_fractions(z₀, z₁, m₁)
m₁₀, m₁₁, m₁₂, m₁₃ = HistogramThresholding.preserved_moments(p₀, z₀, p₁, z₁)
@test m₀ ≈ m₁₀
@test m₁ ≈ m₁₁
@test m₂ ≈ m₁₂
@test m₃ ≈ m₁₃

end


text:
module CSyntax

using CEnum

include("CSwitch.jl")
using .CSwitch

include("CFor.jl")
using .CFor

include("CRef.jl")
using .CRef

include("CStatic.jl")
using .CStatic

# add @c as an alias of @cref
@eval const $(Symbol("@c")) = $(Symbol("@cref"))
export @c

"""
PrefixIncrement Operator
"""
macro +(x)
    @gensym tmp
    quote
        tmp = $(esc(x))
        $(esc(x)) += 1
        tmp
    end
end
@eval const $(Symbol("@++")) = $(Symbol("@+"))

"""
PrefixDecrement Operator
"""
macro -(x)
    @gensym tmp
    quote
        tmp = $(esc(x))
        $(esc(x)) -= 1
        tmp
    end
end
# @-- is invalid

export @+, @++, @-

end # module


text:
# Autogenerated wrapper script for RadeonProRender_jll for x86_64-linux-gnu-cxx03
export libNorthstar64, libProRenderGLTF, libRadeonProRender64, libRprLoadStore64, libTahoe64

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("RadeonProRender")
JLLWrappers.@declare_library_product(libNorthstar64, "libNorthstar64.so")
JLLWrappers.@declare_library_product(libProRenderGLTF, "libProRenderGLTF.so")
JLLWrappers.@declare_library_product(libRadeonProRender64, "libRadeonProRender64.so")
JLLWrappers.@declare_library_product(libRprLoadStore64, "libRprLoadStore64.so")
JLLWrappers.@declare_library_product(libTahoe64, "libTahoe64.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libNorthstar64,
        "lib/libNorthstar64.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libProRenderGLTF,
        "lib/libProRenderGLTF.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libRadeonProRender64,
        "lib/libRadeonProRender64.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libRprLoadStore64,
        "lib/libRprLoadStore64.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libTahoe64,
        "lib/libTahoe64.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function test_token(token::m.Token, expected::m.Token)
    @test token == expected
end

function check_parser_errors(p::m.Parser)
    if !isempty(p.errors)
        msg = join(
            vcat(
                ["parser has $(length(p.errors)) errors"],
                ["parser error: $x" for x in p.errors],
            ),
            "\n",
        )
        error(msg)
    end
end

function test_identifier(expr::m.Expression, value::String)
    @test isa(expr, m.Identifier)
    @test expr.value == value
    @test m.token_literal(expr) == value
end

function test_integer_literal(il::m.Expression, value::Int64)
    @test isa(il, m.IntegerLiteral)
    @test il.value == value
    @test m.token_literal(il) == string(value)
end

function test_boolean_literal(bl::m.Expression, value::Bool)
    @test isa(bl, m.BooleanLiteral)
    @test bl.value == value
    @test m.token_literal(bl) == string(value)
end

function test_null_literal(bl::m.Expression)
    @test isa(bl, m.NullLiteral)
end

function test_literal_expression(expr::m.Expression, expected)
    if isa(expected, Int)
        test_integer_literal(expr, Int64(expected))
    elseif isa(expected, String)
        test_identifier(expr, expected)
    elseif isa(expected, Bool)
        test_boolean_literal(expr, expected)
    elseif isnothing(expected)
        test_null_literal(expr)
    else
        error("unexpected type for expected")
    end
end

function test_infix_expression(expr::m.Expression, left, operator::String, right)
    @test isa(expr, m.InfixExpression)

    test_literal_expression(expr.left, left)

    @test expr.operator == operator

    test_literal_expression(expr.right, right)
end

function test_quote_object(evaluated::m.Object, expected::String)
    @test isa(evaluated, m.QuoteObj)
    @test string(evaluated.node) == expected
end

function test_object(obj::m.Object, expected::Int64)
    @test isa(obj, m.IntegerObj)
    @test obj.value == expected
end

function test_object(obj::m.Object, expected::Bool)
    @test isa(obj, m.BooleanObj)
    @test obj.value == expected
end

function test_object(obj::m.Object, ::Nothing)
    @test obj === m._NULL
end

function test_object(obj::m.Object, expected::String)
    if occursin("error", expected)
        test_error_object(obj, expected)
    else
        test_string_object(obj, expected)
    end
end

function test_object(obj::m.ArrayObj, expected::Vector)
    @test length(obj.elements) == length(expected)

    for (ca, ce) in zip(obj.elements, expected)
        test_object(ca, ce)
    end
end

function test_object(obj::m.HashObj, expected::Dict)
    @test length(obj.pairs) == length(expected)

    for (k, ce) in collect(expected)
        key = m.Object(k)
        test_object(get(obj.pairs, key, nothing), ce)
    end
end

function test_object(obj::m.ErrorObj, expected::String)
    @test obj.message == expected
end

function test_object(obj::m.CompiledFunctionObj, expected::m.Instructions)
    test_instructions(obj.instructions, [expected])
end

function test_string_object(obj::m.StringObj, expected::String)
    @test obj.value == expected
end

function test_instructions(actual, expected)
    concatted = vcat(expected...)

    @test string(actual) == string(concatted)
end

function test_constants(actual, expected)
    @test length(expected) == length(actual)

    for (ca, ce) in zip(actual, expected)
        test_object(ca, ce)
    end
end

function run_compiler_tests(
    input::String,
    expected_constants::Vector,
    expected_instructions::Vector{m.Instructions},
)
    program = m.parse(input)
    c = m.Compiler()
    m.compile!(c, program)

    bc = m.bytecode(c)

    test_instructions(bc.instructions, expected_instructions)
    test_constants(bc.constants, expected_constants)
end

function test_vm(code::String, expected, expected_output::String = "")
    c = IOCapture.capture() do
        m.run(code)
    end
    test_object(c.value, expected)
    @test c.output == expected_output
end


text:
BitIntegers.@define_integers 24

const EDF_SAMPLE_TYPE = Int16
const BDF_SAMPLE_TYPE = Int24
const SUPPORTED_SAMPLE_TYPES = Union{EDF_SAMPLE_TYPE,BDF_SAMPLE_TYPE}

#####
##### `EDF.Signal`
#####

const SIGNAL_HEADER_FIELDS = [(:label, 16),
                              (:transducer_type, 80),
                              (:physical_dimension, 8),
                              (:physical_minimum, 8),
                              (:physical_maximum, 8),
                              (:digital_minimum, 8),
                              (:digital_maximum, 8),
                              (:prefilter, 80),
                              (:samples_per_record, 8)]

"""
    EDF.SignalHeader

Type representing the header for a single EDF signal.

# Fields

* `label::String`: the signal's type/sensor label, see https://www.edfplus.info/specs/edftexts.html#label
* `transducer_type::String`: non-standardized transducer type information
* `physical_dimension::String`: see https://www.edfplus.info/specs/edftexts.html#physidim
* `physical_minimum::Float32`: physical minimum value of the signal
* `physical_maximum::Float32`: physical maximum value of the signal
* `digital_minimum::Float32`: minimum value of the signal that could occur in a data record
* `digital_maximum::Float32`: maximum value of the signal that could occur in a data record
* `prefilter::String`: non-standardized prefiltering information
* `samples_per_record::Int16`: number of samples in a data record (NOT overall)
"""
struct SignalHeader
    label::String
    transducer_type::String
    physical_dimension::String
    physical_minimum::Float32
    physical_maximum::Float32
    digital_minimum::Float32
    digital_maximum::Float32
    prefilter::String
    samples_per_record::Int16
end

"""
    EDF.Signal{T}

Type representing a single EDF signal with sample type `T`.

# Fields

* `header::SignalHeader`
* `samples::Vector{T}`
"""
struct Signal{T}
    header::SignalHeader
    samples::Vector{T}
end

Signal{T}(header::SignalHeader) where {T} = Signal(header, T[])
Signal(header::SignalHeader) = Signal{EDF_SAMPLE_TYPE}(header)

"""
    EDF.decode(signal::Signal)

Return `signal.samples` decoded into the physical units specified by `signal.header`.
"""
function decode(signal::Signal)
    dmin, dmax = signal.header.digital_minimum, signal.header.digital_maximum
    pmin, pmax = signal.header.physical_minimum, signal.header.physical_maximum
    return @. ((signal.samples - dmin) / (dmax - dmin)) * (pmax - pmin) + pmin
end

SignalHeader(signal::Signal) = signal.header

#####
##### `EDF.AnnotationsSignal`
#####

const ANNOTATIONS_SIGNAL_LABEL = "EDF Annotations"

"""
    EDF.TimestampedAnnotationList

A type representing a time-stamped annotations list (TAL).

Note that this type's constructor may attempt to round given `onset_in_seconds` and
`duration_in_seconds` arguments to their nearest representable values in accordance
with the EDF+ specification, which a) represents these values as ASCII, b) constrains
these values to an 8 character limit, and c) does not allow the use of scientific
notation for these fields.

See EDF+ specification for details.

# Fields

* `onset_in_seconds::Float64`: onset w.r.t. recording start time (may be negative)
* `duration_in_seconds::Union{Float64,Nothing}`: duration of this TAL
* `annotations::Vector{String}`: the annotations associated with this TAL
"""
struct TimestampedAnnotationList
    onset_in_seconds::Float64
    duration_in_seconds::Union{Float64,Nothing}
    annotations::Vector{String}
    function TimestampedAnnotationList(onset_in_seconds, duration_in_seconds, annotations)
        onset_in_seconds = _nearest_representable_edf_time_value(onset_in_seconds)
        duration_in_seconds = _nearest_representable_edf_time_value(duration_in_seconds)
        return new(onset_in_seconds, duration_in_seconds, annotations)
    end
end

_nearest_representable_edf_time_value(::Nothing) = nothing

function _nearest_representable_edf_time_value(x)
    return round(x; digits=(8 - (ndigits(floor(Int, x)) + signbit(x) + isinteger(x))))
end

function Base.:(==)(a::TimestampedAnnotationList, b::TimestampedAnnotationList)
    return a.onset_in_seconds == b.onset_in_seconds &&
           a.duration_in_seconds == b.duration_in_seconds &&
           a.annotations == b.annotations
end

"""
    EDF.AnnotationsSignal

Type representing a single EDF Annotations signal.

# Fields

* `samples_per_record::Int16`
* `records::Vector{Vector{TimestampedAnnotationList}}`
"""
struct AnnotationsSignal
    samples_per_record::Int16
    records::Vector{Vector{TimestampedAnnotationList}}
end

function AnnotationsSignal(header::SignalHeader)
    records = Vector{TimestampedAnnotationList}[]
    return AnnotationsSignal(header.samples_per_record, records)
end

"""
    AnnotationsSignal(records::Vector{Vector{TimestampedAnnotationList}})

Return `AnnotationsSignal(samples_per_record, records)` where `samples_per_record`
is the minimum value required to write out each record completely (i.e. the maximum
required `samples_per_record` across all records).
"""
function AnnotationsSignal(records::Vector{Vector{TimestampedAnnotationList}})
    # Actually writing out the TALs in order to calculate the minimum necessary
    # `samples_per_record` is super wasteful from a performance perspective, but
    # is at least robust/self-consistent with the rest of the package. In theory
    # we shouldn't even store this value, and rather always compute it on write,
    # but cleanly refactoring the package to do this would be a more involved
    # change than is meritted at the moment (since normal signals are already
    # treated similarly, i.e. the `SignalHeader` is overly trusted).
    max_bytes_per_record = maximum(sum(write_tal(IOBuffer(), tal) for tal in record) for record in records)
    return AnnotationsSignal(Int16(cld(max_bytes_per_record, 2)), records)
end

function SignalHeader(signal::AnnotationsSignal)
    return SignalHeader("EDF Annotations", "", "", -1, 1, -32768, 32767,
                        "", signal.samples_per_record)
end

#####
##### EDF+ Patient/Recording Metadata
#####

"""
    EDF.PatientID

A type representing the local patient identification field of an EDF+ header.

See EDF+ specification for details.

# Fields

* `code::Union{String,Missing}`
* `sex::Union{Char,Missing}` (`'M'`, `'F'`, or `missing`)
* `birthdate::Union{Date,Missing}`
* `name::Union{String,Missing}`
"""
struct PatientID
    code::Union{String,Missing}
    sex::Union{Char,Missing}
    birthdate::Union{Date,Missing}
    name::Union{String,Missing}
end

"""
    EDF.RecordingID

A type representing the local recording identification field of an EDF+ header.

See EDF+ specification for details.

# Fields

* `startdate::Union{Date,Missing}`
* `admincode::Union{String,Missing}`
* `technician::Union{String,Missing}`
* `equipment::Union{String,Missing}`
"""
struct RecordingID
    startdate::Union{Date,Missing}
    admincode::Union{String,Missing}
    technician::Union{String,Missing}
    equipment::Union{String,Missing}
end

#####
##### `EDF.File`
#####

const BYTES_PER_FILE_HEADER = 256

const BYTES_PER_SIGNAL_HEADER = 256

"""
    EDF.FileHeader

Type representing the parsed header record of an `EDF.File` (excluding signal headers).

# Fields

* `version::String`: data format version
* `patient::Union{String,PatientID}`: local patient identification
* `recording::Union{String,RecordingID}`: local recording identification
* `start::DateTime`: start date/time of the recording
* `is_contiguous::Bool`: if `true`, data records are contiguous; is `true` for non-EDF+-compliant files
* `record_count::Int`: number of data records in the recording
* `seconds_per_record::Float64`: duration of a data record in seconds
"""
struct FileHeader
    version::String
    patient::Union{String,PatientID}
    recording::Union{String,RecordingID}
    start::DateTime
    is_contiguous::Bool
    record_count::Int
    seconds_per_record::Float64
end

"""
    EDF.File{T,I<:IO}

Type representing an EDF file with samples encoded as values of type `T`, which is
`Int16` for EDF files and `Int24` (internally defined) for BDF files.

# Fields

* `io::I`
* `header::FileHeader`
* `signals::Vector{Union{Signal{T},AnnotationsSignal}}`
"""
struct File{T<:SUPPORTED_SAMPLE_TYPES,I<:IO}
    io::I
    header::FileHeader
    signals::Vector{Union{Signal{T},AnnotationsSignal}}
end

function Base.show(io::IO, edf::File{T}) where T
    print(io, "EDF.File with ", length(edf.signals), ' ', 8 * sizeof(T),
          "-bit-encoded signals")
    return nothing
end

Base.close(file::File) = close(file.io)


text:
@testset "IO" begin
    tmpdir = joinpath(tempdir(), "LibPNG")
    isdir(tmpdir) && rm(tmpdir, recursive = true)
    mkdir(tmpdir)

    img = rand(Bool, 5, 5, 5, 5)
    fn = joinpath(tmpdir, "5x5x5x5.png")
    @test_throws ErrorException writeimage(fn, img)

    @testset "Binary Image" begin
        a = rand(Bool, 11, 10)
        fn = joinpath(tmpdir, "binary1.png")
        writeimage(fn, a)
        b1 = readimage(fn)
        @test b1 == convert(Array{Gray{N0f8}}, a)

        a = bitrand(5,5)
        fn = joinpath(tmpdir, "binary2.png")
        writeimage(fn, a)
        b2 = readimage(fn)
        @test b2 == convert(Array{Gray{N0f8}}, a)

        a = colorview(Gray, a)
        fn = joinpath(tmpdir, "binary3.png")
        writeimage(fn, a)
        b3 = readimage(fn)
        @test b3 == convert(Array{Gray{N0f8}}, a)
    end

    @testset "Gray image" begin
        gray = vcat(fill(Gray(1.0), 10, 10), fill(Gray(0.5), 10, 10), fill(Gray(0.0), 10, 10))
        fn = joinpath(tmpdir, "gray1.png")
        writeimage(fn, gray)
        g1 = readimage(fn)
        @test g1 == convert(Array{Gray{N0f8}}, gray)

        gray = rand(Gray{N0f8}, 5, 5)
        fn = joinpath(tmpdir, "gray2.png")
        writeimage(fn, gray)
        g2 = readimage(fn)
        @test g2 == gray
    end

    @testset "Color - RGB" begin
        rgb8 = rand(RGB{N0f8}, 10, 5)
        fn = joinpath(tmpdir, "rgb_n0f8.png")
        writeimage(fn, rgb8)
        r1 = readimage(fn)
        @test r1 == rgb8

        rgb16 = rand(RGB{N0f16}, 5, 5)
        fn = joinpath(tmpdir, "rgb_n0f16.png")
        writeimage(fn, rgb16)
        r1 = readimage(fn)
        @test r1 ==  rgb16
    end

    @testset "Alpha" begin
        # RGBA
        r = RGBA(1.0,0.0,0.0, 0.2)
        g = RGBA(0.0,1.0,0.0, 0.8)
        b = RGBA(0.0,0.0,1.0, 1.0)
        rgba1 = vcat(fill(r, 50,100), fill(g, 50,100), fill(b, 50,100))
        fn = joinpath(tmpdir, "rgba1.png")
        writeimage(fn, rgba1)
        r1 = readimage(fn)
        @test r1 == rgba1

        # GrayA
        r = GrayA(1.0, 0.25)
        g = GrayA(0.5, 0.5)
        b = GrayA(0.0, 0.75)
        graya = vcat(fill(r, 50,100), fill(g, 50,100), fill(b, 50,100))
        fn = joinpath(tmpdir, "graya1.png")
        writeimage(fn, graya)
        g1 = readimage(fn)
        @test g1 == convert(Array{GrayA{N0f8}}, graya)
    end
    # TODO implement palette
end

max_stars_repo_name: FugroRoames/LibPNG

max_issues_repo_name: FugroRoames/LibPNG

max_forks_repo_name: FugroRoames/LibPNG