text:
# Main GaussianProcess type

mutable struct GPMC{X<:MatF64,Y<:AbstractVector{<:Real},M<:Mean,K<:Kernel,L<:Likelihood,
                    D<:KernelData} <: GPBase
    # Observation data
    "Input observations"
    x::X
    "Output observations"
    y::Y

    # Model
    "Mean object"
    mean::M
    "Kernel object"
    kernel::K
    "Likelihood"
    lik::L

    # Auxiliary data
    "Dimension of inputs"
    dim::Int
    "Number of observations"
    nobs::Int
    "Auxiliary observation data (to speed up calculations)"
    data::D
    "Latent (whitened) variables - N(0,1)"
    v::Vector{Float64}
    "Mean values"
    μ::Vector{Float64}
    "`(k + exp(2*obsNoise))`"
    cK::PDMat{Float64,Matrix{Float64}}
    "Log-likelihood"
    ll::Float64
    "Gradient of log-likelihood"
    dll::Vector{Float64}
    "Log-target (marginal log-likelihood + log priors)"
    target::Float64
    "Gradient of log-target (gradient of marginal log-likelihood + gradient of log priors)"
    dtarget::Vector{Float64}

    function GPMC{X,Y,M,K,L}(x::X, y::Y, mean::M, kernel::K, lik::L) where {X,Y,M,K,L}
        dim, nobs = size(x)
        length(y) == nobs || throw(ArgumentError("Input and output observations must have consistent dimensions."))
        data = KernelData(kernel, x)
        gp = new{X,Y,M,K,L,typeof(data)}(x, y, mean, kernel, lik, dim, nobs, data,
                                         zeros(nobs))
        initialise_target!(gp)
    end
end

"""
    GPMC(x, y, mean, kernel, lik)

Fit a Gaussian process to a set of training points. The Gaussian process with
non-Gaussian observations is defined in terms of its user-defined likelihood function,
mean and covaiance (kernel) functions.

The non-Gaussian likelihood is handled by a Monte Carlo method. The latent function
values are represented by centered (whitened) variables ``f(x) = m(x) + Lv`` where
``v ∼ N(0, I)`` and ``LLᵀ = K_θ``.

# Arguments:
- `x::AbstractVecOrMat{Float64}`: Input observations
- `y::AbstractVector{<:Real}`: Output observations
- `mean::Mean`: Mean function
- `kernel::Kernel`: Covariance function
- `lik::Likelihood`: Likelihood function
"""
GPMC(x::MatF64, y::AbstractVector{<:Real}, mean::Mean, kernel::Kernel, lik::Likelihood) =
    GPMC{typeof(x),typeof(y),typeof(mean),typeof(kernel),typeof(lik)}(x, y, mean, kernel,
                                                                      lik)

GPMC(x::VecF64, y::AbstractVector{<:Real}, mean::Mean, kernel::Kernel, lik::Likelihood) =
    GPMC(x', y, mean, kernel, lik)

"""
    GP(x, y, mean::Mean, kernel::Kernel, lik::Likelihood)

Fit a Gaussian process that is defined by its `mean`, its `kernel`, and its likelihood
function `lik` to a set of training points `x` and `y`.

See also: [`GPMC`](@ref)
"""
GP(x::AbstractVecOrMat{Float64}, y::AbstractVector{<:Real}, mean::Mean, kernel::Kernel,
   lik::Likelihood) = GPMC(x, y, mean, kernel, lik)

"""
    initialise_ll!(gp::GPMC)

Initialise the log-likelihood of Gaussian process `gp`.
"""
function initialise_ll!(gp::GPMC)
    # log p(Y|v,θ)
    gp.μ = mean(gp.mean,gp.x)
    Σ = cov(gp.kernel, gp.x, gp.data)
    gp.cK = PDMat(Σ + 1e-6*I)
    F = unwhiten(gp.cK,gp.v) + gp.μ
    gp.ll = sum(log_dens(gp.lik,F,gp.y)) #Log-likelihood
    gp
end

"""
    update_cK!(gp::GPMC)

Update the covariance matrix and its Cholesky decomposition of Gaussian process `gp`.
"""
function update_cK!(gp::GPMC)
    old_cK = gp.cK
    Σbuffer = old_cK.mat
    cov!(Σbuffer, gp.kernel, gp.x, gp.data)
    for i in 1:gp.nobs
        Σbuffer[i,i] += 1e-6 # no logNoise for GPMC
    end
    chol_buffer = old_cK.chol.factors
    copyto!(chol_buffer, Σbuffer)
    chol = cholesky!(Symmetric(chol_buffer))
    gp.cK = PDMat(Σbuffer, chol)
    gp
end

# modification of initialise_ll! that reuses existing matrices to avoid
# unnecessary memory allocations, which speeds things up significantly
function update_ll!(gp::GPMC; process::Bool=true, lik::Bool=true, domean::Bool=true, kern::Bool=true)
    if kern
        # only need to update the covariance matrix
        # if the covariance parameters have changed
        update_cK!(gp)
    end
    gp.μ = mean(gp.mean,gp.x)
    F = unwhiten(gp.cK,gp.v) + gp.μ
    gp.ll = sum(log_dens(gp.lik,F,gp.y)) #Log-likelihood
    gp
end

"""
     update_dll!(gp::GPMC, ...)

Update the gradient of the log-likelihood of Gaussian process `gp`.
"""
function update_dll!(gp::GPMC, Kgrad::MatF64, L_bar::MatF64;
    process::Bool=true, # include gradient components for the process itself
    lik::Bool=true,  # include gradient components for the likelihood parameters
    domean::Bool=true, # include gradient components for the mean parameters
    kern::Bool=true, # include gradient components for the spatial kernel parameters
    )
    size(Kgrad) == (gp.nobs, gp.nobs) || throw(ArgumentError,
    @sprintf("Buffer for Kgrad should be a %dx%d matrix, not %dx%d",
             gp.nobs, gp.nobs,
             size(Kgrad,1), size(Kgrad,2)))
    size(L_bar) == (gp.nobs, gp.nobs) || throw(ArgumentError,
    @sprintf("Buffer for L_bar should be a %dx%d matrix, not %dx%d",
             gp.nobs, gp.nobs,
             size(L_bar,1), size(L_bar,2)))

    n_lik_params = num_params(gp.lik)
    n_mean_params = num_params(gp.mean)
    n_kern_params = num_params(gp.kernel)

    gp.dll = Array{Float64}(undef, process * gp.nobs + lik * n_lik_params +
                            domean * n_mean_params + kern * n_kern_params)

    Lv = unwhiten(gp.cK, gp.v)
    dl_df = dlog_dens_df(gp.lik, Lv + gp.μ, gp.y)
    i = 1
    if process
        mul!(view(gp.dll, i:i+gp.nobs-1), gp.cK.chol.U, dl_df)
        i += gp.nobs
    end
    if lik && n_lik_params > 0
        Lgrads = dlog_dens_dθ(gp.lik, Lv + gp.μ, gp.y)
        for j in 1:n_lik_params
            gp.dll[i] = sum(Lgrads[:,j])
            i += 1
        end
    end
    if domean && n_mean_params > 0
        Mgrads = grad_stack(gp.mean, gp.x)
        for j in 1:n_mean_params
            gp.dll[i] = dot(dl_df,Mgrads[:,j])
            i += 1
        end
    end
    if kern
        fill!(L_bar, 0)
        BLAS.ger!(1.0, dl_df, gp.v, L_bar)
        tril!(L_bar)
        # ToDo:
        # the following two steps allocates memory
        # and are fickle, reaching into the internal
        # implementation of the cholesky decomposition
        L = gp.cK.chol.L.data
        tril!(L)
        #
        chol_unblocked_rev!(L, L_bar)
        for iparam in 1:n_kern_params
            grad_slice!(Kgrad, gp.kernel, gp.x, gp.data, iparam)
            gp.dll[i] = dot(Kgrad, L_bar)
            i+=1
        end
    end

    gp
end

function update_ll_and_dll!(gp::GPMC, Kgrad::MatF64, L_bar::MatF64; kwargs...)
    update_ll!(gp; kwargs...)
    update_dll!(gp, Kgrad, L_bar; kwargs...)
end

"""
    initialise_target!(gp::GPMC)

Initialise the log-posterior
```math
\\log p(θ, v | y) ∝ \\log p(y | v, θ) + \\log p(v) + \\log p(θ)
```
of a Gaussian process `gp`.
"""
function initialise_target!(gp::GPMC)
    initialise_ll!(gp)
    gp.target = gp.ll - (sum(abs2, gp.v) + log2π * gp.nobs) / 2 +
        prior_logpdf(gp.lik) + prior_logpdf(gp.mean) + prior_logpdf(gp.kernel)
    gp
end

"""
    update_target!(gp::GPMC, ...)

Update the log-posterior
```math
\\log p(θ, v | y) ∝ \\log p(y | v, θ) + \\log p(v) + \\log p(θ)
```
of a Gaussian process `gp`.
"""
function update_target!(gp::GPMC; process::Bool=true, lik::Bool=true, domean::Bool=true, kern::Bool=true)
    update_ll!(gp; process=process, lik=lik, domean=domean, kern=kern)
    gp.target = gp.ll - (sum(abs2, gp.v) + log2π * gp.nobs) / 2 +
        prior_logpdf(gp.lik) + prior_logpdf(gp.mean) + prior_logpdf(gp.kernel)
    gp
end

function update_dtarget!(gp::GPMC, Kgrad::MatF64, L_bar::MatF64; kwargs...)
    update_dll!(gp, Kgrad, L_bar; kwargs...)
    gp.dtarget = gp.dll + prior_gradlogpdf(gp; kwargs...)
    gp
end

"""
    update_target_and_dtarget!(gp::GPMC, ...)

Update the log-posterior
```math
\\log p(θ, v | y) ∝ \\log p(y | v, θ) + \\log p(v) + \\log p(θ)
```
of a Gaussian process `gp` and its derivative.
"""
function update_target_and_dtarget!(gp::GPMC, Kgrad::MatF64, L_bar::MatF64; kwargs...)
    update_target!(gp; kwargs...)
    update_dtarget!(gp, Kgrad, L_bar; kwargs...)
end

function update_target_and_dtarget!(gp::GPMC; kwargs...)
    Kgrad = Array{Float64}(undef, gp.nobs, gp.nobs)
    L_bar = Array{Float64}(undef, gp.nobs, gp.nobs)
    update_target_and_dtarget!(gp, Kgrad, L_bar; kwargs...)
end

"""
    predict_y(gp::GPMC, x::Union{Vector{Float64},Matrix{Float64}}[; full_cov::Bool=false])

Return the predictive mean and variance of Gaussian Process `gp` at specfic points which
are given as columns of matrix `x`. If `full_cov` is `true`, the full covariance matrix is
returned instead of only variances.
"""
function predict_y(gp::GPMC, x::MatF64; full_cov::Bool=false)
    μ, σ2 = predict_f(gp, x; full_cov=full_cov)
    return predict_obs(gp.lik, μ, σ2)
end

# 1D Case for prediction
predict_y(gp::GPMC, x::VecF64; full_cov::Bool=false) = predict_y(gp, x'; full_cov=full_cov)

## compute predictions
function _predict(gp::GPMC, x::MatF64)
    cK = cov(gp.kernel, gp.x, x)
    Lck = whiten(gp.cK, cK)
    fmu =  mean(gp.mean,x) + Lck'gp.v     # Predictive mean
    Sigma_raw = cov(gp.kernel, x) - Lck'Lck # Predictive covariance
    # Add jitter to get stable covariance
    fSigma = tolerant_PDMat(Sigma_raw)
    return fmu, fSigma
end

# Sample from functions from the GP
function Random.rand!(gp::GPMC, x::MatF64, A::DenseMatrix)
    nobs = size(x,2)
    n_sample = size(A,2)

    if gp.nobs == 0
        # Prior mean and covariance
        μ = mean(gp.mean, x);
        Σraw = cov(gp.kernel, x);
        # Add jitter to get stable covariance
        Σ = tolerant_PDMat(Σraw)
    else
        # Posterior mean and covariance
        μ, Σ = predict_f(gp, x; full_cov=true)
    end
    return broadcast!(+, A, μ, unwhiten!(Σ,randn(nobs, n_sample)))
end

Random.rand(gp::GPMC, x::MatF64, n::Int) = rand!(gp, x, Array{Float64}(undef, size(x, 2), n))

# Sample from 1D GP
Random.rand(gp::GPMC, x::VecF64, n::Int) = rand(gp, x', n)

# Generate only one sample from the GP and returns a vector
Random.rand(gp::GPMC, x::MatF64) = vec(rand(gp,x,1))
Random.rand(gp::GPMC, x::VecF64) = vec(rand(gp,x',1))

function get_params(gp::GPMC; lik::Bool=true, domean::Bool=true, kern::Bool=true)
    params = Float64[]
    append!(params, gp.v)
    if lik  && num_params(gp.lik)>0
        append!(params, get_params(gp.lik))
    end
    if domean && num_params(gp.mean)>0
        append!(params, get_params(gp.mean))
    end
    if kern
        append!(params, get_params(gp.kernel))
    end
    return params
end

function num_params(gp::GPMC; lik::Bool=true, domean::Bool=true, kern::Bool=true)
    n = length(gp.v)
    lik && (n += num_params(gp.lik))
    domean && (n += num_params(gp.mean))
    kern && (n += num_params(gp.kernel))
    n
end

function set_params!(gp::GPMC, hyp::VecF64; process::Bool=true, lik::Bool=true, domean::Bool=true, kern::Bool=true)
    n_lik_params = num_params(gp.lik)
    n_mean_params = num_params(gp.mean)
    n_kern_params = num_params(gp.kernel)

    i = 1
    if process
        gp.v = hyp[1:gp.nobs]
        i += gp.nobs
    end
    if lik  && n_lik_params>0
        set_params!(gp.lik, hyp[i:i+n_lik_params-1]);
        i += n_lik_params
    end
    if domean && n_mean_params>0
        set_params!(gp.mean, hyp[i:i+n_mean_params-1])
        i += n_mean_params
    end
    if kern
        set_params!(gp.kernel, hyp[i:i+n_kern_params-1])
        i += n_kern_params
    end
end

function prior_gradlogpdf(gp::GPMC; process::Bool=true, lik::Bool=true, domean::Bool=true, kern::Bool=true)
    if process
        grad = -gp.v
    else
        grad = Float64[]
    end
    if lik
        append!(grad, prior_gradlogpdf(gp.lik))
    end
    if domean
        append!(grad, prior_gradlogpdf(gp.mean))
    end
    if kern
        append!(grad, prior_gradlogpdf(gp.kernel))
    end
    return grad
end

function Base.show(io::IO, gp::GPMC)
    println(io, "GP Monte Carlo object:")
    println(io, "  Dim = ", gp.dim)
    println(io, "  Number of observations = ", gp.nobs)
    println(io, "  Mean function:")
    show(io, gp.mean, 2)
    println(io, "\n  Kernel:")
    show(io, gp.kernel, 2)
    println(io, "\n  Likelihood:")
    show(io, gp.lik, 2)
    if gp.nobs == 0
        println("\n  No observation data")
    else
        println(io, "\n  Input observations = ")
        show(io, gp.x)
        print(io,"\n  Output observations = ")
        show(io, gp.y)
        print(io,"\n  Log-posterior = ")
        show(io, round(gp.target; digits = 3))
    end
end


text:
using Logging
push!(LOAD_PATH, dirname(@__DIR__))
using EnergySystemModeling
cd("D:\\Eigene Dateien\\Documents\\GitHub\\EnergySystemModeling.jl\\examples")

@info "Loading parameters"
constants_path = "constants"
structure = "8nodes"
structures_path = joinpath("structures",structure)
instance = "small"
instances_path = joinpath(structures_path,"instances",instance)

@info "Creating output directory"
output_dir = joinpath(instances_path,"output")
results_dir = "results"
plots_dir = "plots"
csv_dir = "csv"

mkpath(output_dir)
mkpath(joinpath(output_dir,results_dir))
mkpath(joinpath(output_dir,plots_dir))
mkpath(joinpath(output_dir,plots_dir,"pdf"))
mkpath(joinpath(output_dir,plots_dir,"png"))
mkpath(joinpath(output_dir,csv_dir))

parameters = Params(constants_path, instances_path)
specs = Specs(
        transmission=true,
        renewable_target=true,
        carbon_cap=true,
        nuclear_limit=false,
        storage=true,
        ramping=true,
        voltage_angles=false,
        hydro=true
)

@info "Creating the energy system model"
(model, VariablesDict, ObjectivesDict) = EnergySystemModel(parameters, specs)

@info "Optimizing the model"
using Gurobi, JuMP
optimizer = optimizer_with_attributes(Gurobi.Optimizer, "TimeLimit" => 60*60*2,
                                      "LogFile" => joinpath(output_dir, "gurobi.log"))
set_optimizer(model, optimizer)
set_optimizer_attributes(model, "Method" => 2)
set_optimizer_attributes(model, "Crossover" => 0)
set_optimizer_attributes(model, "NumericFocus" => 1)

optimize!(model)

@info "Extracting results"
variables = JuMPVar(model, VariablesDict)
objectives = JuMPObj(model, ObjectivesDict)
expressions = Expressions(parameters, variables)

@info "Saving results (JSON)"
save_json(specs, joinpath(output_dir, results_dir, "specs.json"))
save_json(parameters, joinpath(output_dir, results_dir, "parameters.json"))
save_json(variables, joinpath(output_dir, results_dir, "variables.json"))
save_json(objectives, joinpath(output_dir, results_dir, "objectives.json"))

@info "Plotting"
using Plots
ENV["GKSwstype"]="nul"                      # Prevent opening plots windows
using StatsPlots
gr()

## Plotting part 1: Objective function values
@info "Plotting OF"
savefig(plot_objective_values(objectives),
        joinpath(output_dir, plots_dir, "pdf", "objectives.pdf"))
savefig(plot_objective_values(objectives),
        joinpath(output_dir, plots_dir, "png", "objectives.png"))

## Plotting part 2: Dispatch and storage
@info "Plotting dispatch and storage levels"

for n in parameters.N
    savefig(plot_generation_dispatch(parameters, variables, expressions, n),
            joinpath(output_dir, plots_dir,"pdf","generation_dispatch_n$n.pdf"))
    savefig(plot_generation_capacities(parameters, variables, expressions, n),
            joinpath(output_dir, plots_dir,"pdf","generation_capacities_n$n.pdf"))
    savefig(plot_storage_level(parameters, variables, expressions, n),
            joinpath(output_dir, plots_dir,"pdf","storage_n$n.pdf"))
    savefig(plot_box(parameters, variables, expressions, n),
           joinpath(output_dir, plots_dir,"pdf","boxplot$n.pdf"))
    savefig(plot_generation_dispatch(parameters, variables, expressions, n),
           joinpath(output_dir, plots_dir,"png","generation_dispatch_n$n.png"))
    savefig(plot_generation_capacities(parameters, variables, expressions, n),
           joinpath(output_dir, plots_dir,"png","generation_capacities_n$n.png"))
    savefig(plot_storage_level(parameters, variables, expressions, n),
           joinpath(output_dir, plots_dir,"png","storage_n$n.png"))
    savefig(plot_box(parameters, variables, expressions, n),
          joinpath(output_dir, plots_dir,"png","boxplot$n.png"))

end

## Plotting part 3: Storage capacity
@info "Plotting storage capacities"

savefig(plot_storage_capacities(parameters, variables, expressions),
        joinpath(output_dir, plots_dir,"pdf","storage_capacities.pdf")
        )
savefig(plot_storage_capacities(parameters, variables, expressions),
        joinpath(output_dir, plots_dir,"png","storage_capacities.png")
        )

## Plotting part 4: Generation technologies dispatch levels
@info "Plotting dispatch all"

savefig(plot_box_all(parameters, variables, expressions),
        joinpath(output_dir, plots_dir,"pdf","boxplotall.pdf"))
savefig(plot_box_all(parameters, variables, expressions),
        joinpath(output_dir, plots_dir,"png","boxplotall.png"))

## Plotting part 5: Generation capacities (including hydro)
@info "Plotting capacities stacked"

savefig(plot_generation_capacities_stacked(parameters, variables, expressions),
            joinpath(output_dir, plots_dir,"pdf","generation_capacities_stacked.pdf"))
        savefig(plot_generation_capacities_stacked(parameters, variables, expressions),
            joinpath(output_dir, plots_dir,"png","generation_capacities_stacked.png"))

## Plotting part 6: Consolidated dispatch vs demand
@info "Plotting dispatch all (boxplot)"

savefig(plot_dispatch_bars(parameters, variables, expressions),
        joinpath(output_dir, plots_dir,"pdf","dispatchbars.pdf"))
savefig(plot_dispatch_bars(parameters, variables, expressions),
        joinpath(output_dir, plots_dir,"png","dispatchbars.png"))

## Plotting part 7: Transmission flow (per line)
@info "Plotting transmission flow"

for l in 1:length(parameters.L_ind)
    savefig(plot_transmission_flow(parameters, variables, expressions, l),
            joinpath(output_dir, plots_dir,"pdf","transmission_flow_L$l.pdf"))
    savefig(plot_transmission_flow(parameters, variables, expressions, l),
            joinpath(output_dir, plots_dir,"png","transmission_flow_L$l.png"))

end

## Plotting part 8: Transmission capacities
@info "Plotting transmission capacities"

savefig(plot_transmission_capacities(parameters, variables, expressions),
        joinpath(output_dir, plots_dir,"pdf","transmission_capacities.pdf"))
savefig(plot_transmission_capacities(parameters, variables, expressions),
        joinpath(output_dir, plots_dir,"png","transmission_capacities.png"))

## Plotting part 9: Consolidated transmission flow
@info "Plotting transmission flow"

savefig(plot_transmission_bars(parameters, variables, expressions),
        joinpath(output_dir, plots_dir,"pdf","transmission_bars.pdf"))
savefig(plot_transmission_bars(parameters, variables, expressions),
        joinpath(output_dir, plots_dir,"png","transmission_bars.png"))

## Plotting part 10: Lost of load
@info "Plotting LoL"

savefig(plot_loss_of_load(parameters, variables, expressions),
        joinpath(output_dir, plots_dir,"pdf","loss_of_load.pdf"))
savefig(plot_loss_of_load(parameters, variables, expressions),
        joinpath(output_dir, plots_dir,"png","loss_of_load.png"))

## Plotting in development
@info "Plotting stacked dispatch (not ready)"

ext: jl


text:
""""""
function standardised(unused...; transform=baselinetransform,
    filter=baselinefilter)
    function f(F)
        Fₗ = zeros(size(F))
        Fₕ = F # No baseline, so the best we can do is treat the test as its own baseline
        intervalscaled(Fₕ, Fₗ)(F)
    end
    return f
end
export standardised


text:
@testset "141.linked-list-cycle.jl" begin
    @test has_cycle(ListNode{Int}([1])) == false

    root = ListNode{Int}([1, 2])
    root.next.next = root
    @test has_cycle(root) == true

    root = ListNode{Int}([3, 2, 0, -4])
    root.next.next.next.next = root.next
    @test has_cycle(root) == true
end


text:
module DISSEQTPlots

Pkg.installed("DISSEQT")==nothing && warn("Module DISSEQT not installed. Please refer to DISSEQTPlots installation instructions at https://github.com/rasmushenningsson/DISSEQTPlots.jl")

using DISSEQT
using Gadfly
using DataFrames
using Colors

Pkg.installed("PlotlyJS")!=nothing && (import PlotlyJS; const py=PlotlyJS) # import to avoid clash with Gadfly

export
    groupcolors,
    linearregressionplot,
    pairwisescatterplot,
    talusplot,
    fitnesslandscapeplot,
    saveplot

include("colors.jl")
include("stopcodons.jl")
include("pairwise.jl")
include("talus.jl")
include("fitnesslandscape.jl")
include("save.jl")
Pkg.installed("PlotlyJS")!=nothing && include("saveplotly.jl")

end


text:
# Regular Expression
# ==================

module RegExp

import Automa: ByteSet

function gen_empty_names()
    return Symbol[]
end

mutable struct RE
    head::Symbol
    args::Vector
    actions::Dict{Symbol, Vector{Symbol}}
    when::Union{Symbol, Nothing}
end

function RE(head::Symbol, args::Vector)
    return RE(head, args, Dict{Symbol, Vector{Symbol}}(), nothing)
end

const Primitive = Union{RE, ByteSet, UInt8, UnitRange{UInt8}, Char, String, Vector{UInt8}}

function primitive(re::RE)
    return re
end

const PRIMITIVE = (:set, :byte, :range, :class, :cclass, :char, :str, :bytes)

function primitive(set::ByteSet)
    return RE(:set, [set])
end

function primitive(byte::UInt8)
    return RE(:byte, [byte])
end

function primitive(range::UnitRange{UInt8})
    return RE(:range, [range])
end

function primitive(char::Char)
    return RE(:char, [char])
end

function primitive(str::String)
    return RE(:str, [str])
end

function primitive(bs::Vector{UInt8})
    return RE(:bytes, copy(bs))
end

function cat(xs::Primitive...)
    return RE(:cat, [map(primitive, xs)...])
end

function alt(x::Primitive, xs::Primitive...)
    return RE(:alt, [primitive(x), map(primitive, xs)...])
end

function rep(x::Primitive)
    return RE(:rep, [primitive(x)])
end

function rep1(x::Primitive)
    return RE(:rep1, [primitive(x)])
end

function opt(x::Primitive)
    return RE(:opt, [primitive(x)])
end

function isec(x::Primitive, y::Primitive)
    return RE(:isec, [primitive(x), primitive(y)])
end

function diff(x::Primitive, y::Primitive)
    return RE(:diff, [primitive(x), primitive(y)])
end

function neg(x::Primitive)
    return RE(:neg, [primitive(x)])
end

function any()
    return primitive(0x00:0xff)
end

function ascii()
    return primitive(0x00:0x7f)
end

function space()
    return primitive(ByteSet([UInt8(c) for c in "\t\v\f\n\r "]))
end

Base.:*(re1::RE, re2::RE) = cat(re1, re2)
Base.:|(re1::RE, re2::RE) = alt(re1, re2)
Base.:&(re1::RE, re2::RE) = isec(re1, re2)
Base.:\(re1::RE, re2::RE) = diff(re1, re2)
Base.:!(re::RE) = neg(re)

macro re_str(str::String)
    parse(str)
end

const METACHAR = raw".*+?()[]\|-^"

# Parse a regular expression string using the shunting-yard algorithm.
function parse(str::String)
    # stacks
    operands = RE[]
    operators = Symbol[]

    function pop_and_apply!()
        op = pop!(operators)
        if op == :rep || op == :rep1 || op == :opt
            arg = pop!(operands)
            push!(operands, RE(op, [arg]))
        elseif op == :alt
            arg2 = pop!(operands)
            arg1 = pop!(operands)
            push!(operands, RE(:alt, [arg1, arg2]))
        elseif op == :cat
            arg2 = pop!(operands)
            arg1 = pop!(operands)
            push!(operands, RE(:cat, [arg1, arg2]))
        else
            error(op)
        end
    end

    cs = iterate(str)
    if cs === nothing
        return RE(:cat, [])
    end
    need_cat = false
    while cs !== nothing
        c, s = cs
        # @show c operands operators
        if need_cat && c ∉ ('*', '+', '?', '|', ')')
            while !isempty(operators) && prec(:cat) ≤ prec(last(operators))
                pop_and_apply!()
            end
            push!(operators, :cat)
        end
        need_cat = c ∉ ('|', '(')
        if c == '*'
            while !isempty(operators) && prec(:rep) ≤ prec(last(operators))
                pop_and_apply!()
            end
            push!(operators, :rep)
        elseif c == '+'
            while !isempty(operators) && prec(:rep1) ≤ prec(last(operators))
                pop_and_apply!()
            end
            push!(operators, :rep1)
        elseif c == '?'
            while !isempty(operators) && prec(:opt) ≤ prec(last(operators))
                pop_and_apply!()
            end
            push!(operators, :opt)
        elseif c == '|'
            while !isempty(operators) && prec(:alt) ≤ prec(last(operators))
                pop_and_apply!()
            end
            push!(operators, :alt)
        elseif c == '('
            push!(operators, :lparen)
        elseif c == ')'
            while !isempty(operators) && last(operators) != :lparen
                pop_and_apply!()
            end
            pop!(operators)
        elseif c == '['
            class, cs = parse_class(str, s)
            push!(operands, class)
            continue
        elseif c == '.'
            push!(operands, any())
        elseif c == '\\'
            if iterate(str, s) === nothing
                c = '\\'
            else
                c, s = unescape(str, s)
            end
            push!(operands, primitive(c))
        else
            push!(operands, primitive(c))
        end
        cs = iterate(str, s)
    end

    while !isempty(operators)
        pop_and_apply!()
    end

    @assert length(operands) == 1
    return first(operands)
end

# Operator's precedence.
function prec(op::Symbol)
    if op == :rep || op == :rep1 || op == :opt
        return 3
    elseif op == :cat
        return 2
    elseif op == :alt
        return 1
    elseif op == :lparen
        return 0
    else
        @assert false
    end
end

function parse_class(str, s)
    chars = Tuple{Bool, Char}[]
    cs = iterate(str, s)
    while cs !== nothing
        c, s = cs
        if c == ']'
            cs = iterate(str, s)
            break
        elseif c == '\\'
            if iterate(str, s) === nothing
                error("missing ]")
            end
            c, s = unescape(str, s)
            push!(chars, (true, c))
        else
            push!(chars, (false, c))
        end
        cs = iterate(str, s)
    end
    if !isempty(chars) && !first(chars)[1] && first(chars)[2] == '^'
        head = :cclass
        popfirst!(chars)
    else
        head = :class
    end
    if isempty(chars)
        error("empty class")
    end

    args = []
    while !isempty(chars)
        c = popfirst!(chars)[2]
        if !isempty(chars) && !first(chars)[1] && first(chars)[2] == '-' && length(chars) ≥ 2
            push!(args, UInt8(c):UInt8(chars[2][2]))
            popfirst!(chars)
            popfirst!(chars)
        else
            push!(args, UInt8(c):UInt8(c))
        end
    end
    return RE(head, args), cs
end

function unescape(str::String, s::Int)
    invalid() = throw(ArgumentError("invalid escape sequence"))
    ishex(b) = '0' ≤ b ≤ '9' || 'A' ≤ b ≤ 'F' || 'a' ≤ b ≤ 'f'
    cs = iterate(str, s)
    cs === nothing && invalid()
    c, s = cs
    if c == 'a'
        return '\a', s
    elseif c == 'b'
        return '\b', s
    elseif c == 't'
        return '\t', s
    elseif c == 'n'
        return '\n', s
    elseif c == 'v'
        return '\v', s
    elseif c == 'r'
        return '\r', s
    elseif c == 'f'
        return '\f', s
    elseif c == '0'
        return '\0', s
    elseif c ∈ METACHAR
        return c, s
    elseif c == 'x'
        cs1 = iterate(str, s)
        (cs1 === nothing || !ishex(cs1[1])) && invalid()
        cs2 = iterate(str, cs1[2])
        (cs2 === nothing || !ishex(cs2[1])) && invalid()
        c1, c2 = cs1[1], cs2[1]
        return first(unescape_string("\\x$(c1)$(c2)")), cs2[2]
    elseif c == 'u' || c == 'U'
        throw(ArgumentError("escaped Unicode sequence is not supported"))
    else
        throw(ArgumentError("invalid escape sequence: \\$(c)"))
    end
end

function shallow_desugar(re::RE)
    head = re.head
    args = re.args
    if head == :rep1
        return RE(:cat, [args[1], rep(args[1])])
    elseif head == :opt
        return RE(:alt, [args[1], RE(:cat, [])])
    elseif head == :neg
        return RE(:diff, [rep(any()), args[1]])
    elseif head == :byte
        return RE(:set, [ByteSet(args[1])])
    elseif head == :range
        return RE(:set, [ByteSet(args[1])])
    elseif head == :class
        return RE(:set, [foldl(union, map(ByteSet, args), init=ByteSet())])
    elseif head == :cclass
        return RE(:set, [foldl(setdiff, map(ByteSet, args), init=ByteSet(0x00:0xff))])
    elseif head == :char
        bytes = convert(Vector{UInt8}, codeunits(string(args[1])))
        return RE(:cat, [RE(:set, [ByteSet(b)]) for b in bytes])
    elseif head == :str
        bytes = convert(Vector{UInt8}, codeunits(args[1]))
        return RE(:cat, [RE(:set, [ByteSet(b)]) for b in bytes])
    elseif head == :bytes
        return RE(:cat, [RE(:set, [ByteSet(b)]) for b in args])
    else
        if head ∉ (:set, :cat, :alt, :rep, :isec, :diff)
            error("cannot desugar ':$(head)'")
        end
        return RE(head, args)
    end
end

end


text:
function _generate_w(latent::Latent1D, parameters::Parameters2PL)
    return Distributions.rand(PolyaGamma(1, parameters.a * (latent.val - parameters.b)))
end

function generate_w(examinee::AbstractExaminee, items::Vector{<:AbstractItem})
    return map(
        i -> PolyaGammaSample(
            i.idx,
            examinee.idx,
            _generate_w(i.parameters, examinee.latent)
        ),
        items,
    )
end


text:
@inline _maybefirst(x) = x
@inline _maybefirst(x::VecUnroll) = first(data(x))
@inline _maybetail(x) = x
@inline _maybetail(::VecUnroll{0}) = ()
@inline _maybetail(x::VecUnroll) = VecUnroll(Base.tail(data(x)))

@inline _vload_map(_, ::Tuple, ::Tuple{}, __, ___) = ()
@inline function _vload_map(p, i, m, ::J, ::A) where {J,A}
  x = _vload(p, map(_maybefirst, i), first(m), J(), A())
  r = _vload_map(p, map(_maybetail, i), Base.tail(m), J(), A())
  (x, r...)
end

@inline function _vload(
  p::AbstractStridedPointer,
  i::Tuple{Vararg{Union{IntegerIndex,MM,VecUnroll{N,<:Any,<:Any,<:IntegerIndex}}}},
  m::VecUnroll{N,<:Any,Bit},
  ::J,
  ::A,
) where {N,J,A}
  VecUnroll(_vload_map(p, i, data(m), J(), A()))
end


text:
using CSFML.LibCSFML, Gtk.ShortNames, Colors, Graphics, Cairo

#------------ input code ----------------------#

mutable struct Joystick
    devnum::Int
    isconnected::Bool
    hasXaxis::Bool
    nbuttons::Int
    pressed::Vector{Bool}
    ypos::Int
    xpos::Int
    name::String
    Joystick(n, b=2, c=false, x=true) = new(n, c, x, b, fill(false, b), 0, 0)
end

const devnum = 0
const buttons = 2
const jstick = Joystick(devnum, buttons)

function polljoystick(jstick, sleepinterval=0.05)
    while !sfJoystick_isConnected(jstick.devnum)
        sleep(0.25) # wait till connected
        sfJoystick_update()
    end
    jstick.name =  sfJoystick_getIdentification(jstick.devnum).name
    jstick.isconnected = true
    jstick.hasXaxis = sfJoystick_hasAxis(jstick.devnum, 0)
    jstick.nbuttons = sfJoystick_getButtonCount(jstick.devnum)
    while true
        sfJoystick_update()
        for i in 1:jstick.nbuttons
            jstick.pressed[i] =  sfJoystick_isButtonPressed(jstick.devnum, i - 1)
        end
        jstick.ypos = sfJoystick_getAxisPosition(jstick.devnum, 1)
        jstick.xpos = sfJoystick_getAxisPosition(joystick.devnum, 0)
        yield()
        sleep(sleepinterval)
    end
end

#------------ output code -------------------#

makelabel() = "Button 1: " * (jstick.pressed[1] ? "DOWN" : "UP") *
    "  Button 2: " * (jstick.pressed[2] ? "DOWN" : "UP")

const fontpointsize = 80
const wid = 500
const hei = 500
win = Window("Cursor Task", wid, hei) |> (Frame() |> (vbox = Box(:v)))
set_gtk_property!(vbox, :expand, true)
can = Canvas(wid, hei)
label = Label(makelabel())
push!(vbox, can, label)

joytoxpos() = div((jstick.xpos + 100) * width(can), 200)
joytoypos() = div((jstick.ypos + 100) * height(can), 200)
Gtk.showall(win)

@guarded draw(can) do widget
    ctx = getgc(can)
    select_font_face(ctx, "Courier", Cairo.FONT_SLANT_NORMAL, Cairo.FONT_WEIGHT_BOLD)
    set_font_size(ctx, fontpointsize)
    set_source(ctx, colorant"red")
    move_to(ctx, joytoxpos(), joytoypos())
    show_text(ctx, "+")
    set_gtk_property!(label, :label, makelabel())
end

@async polljoystick(jstick)

while true
    draw(can)
    sleep(0.2)
    yield()
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#############################################################################
# ECOS.jl
# Wrapper around the ECOS solver https://github.com/ifa-ethz/ecos
# See http://github.com/JuliaOpt/ECOS.jl
#############################################################################
# test/runtests.jl
# Test the ECOS.jl solver wrapper
#############################################################################

using ECOS

using Compat
using Compat.Test

@testset "Test the direct interface" begin
    include("direct.jl")
end

@testset "Test passing options" begin
    include("options.jl")
end

@testset "MathProgBase" begin
    include("MPB_wrapper.jl")
end

@testset "MathOptInterface" begin
    include("MOI_wrapper.jl")
end


text:
module TestBroadcasting

include("preamble.jl")

@testset "count_sum" begin
    @testset "$fold" for fold in [foldl, simple_reduce, reduce_bs1, reduce]
        foldl = nothing

        count_sum(xs, xf = Map(identity)) =
            fold(+, xf |> Map(x -> (1, x)) |> Broadcasting(), xs)
        @test count_sum([5, 2, 6, 8, 3]) == (5, 24)
        @test count_sum([5, 2, 6, 8, 3], Filter(isodd)) == (2, 8)
        @test_throws EmptyResultError count_sum(1:0)
        @test_throws EmptyResultError count_sum(1:10, Filter(_ -> false))
    end
end

@testset "partially started" begin
    @testset "$fold" for fold in [foldl, simple_reduce, reduce_bs1, reduce]
        foldl = nothing
        f(xs) = fold(max, Map(x -> (x + 1, x)) |> Broadcasting() |> Filter(isodd), xs)
        @test f([5, 2, 6, 8, 3]) == (9, 5)
        @test_throws EmptyResultError f(2:2:8)
    end
end

@testset "init" begin
    @testset "$fold" for fold in [foldl, simple_reduce, reduce_bs1, reduce]
        foldl = nothing
        @test fold(+, Broadcasting(), 1:9) == 45
        @test fold(+, Broadcasting(), 1:9, init = [0]) == [45]
    end
end

sample_foldl_allocations(repeat::Integer, args...) =
    map(1:repeat) do _
        @allocated foldl(args...)
    end

@testset "allocation" begin
    args = (+, Broadcasting())
    n = 1000
    xs = [ones(n) for _ in 1:n]
    sample_foldl_allocations(2, args..., xs)
    baseline = minimum(sample_foldl_allocations(100, args..., xs))
    repeat10x = minimum(sample_foldl_allocations(100, args..., repeat(xs, 10)))
    @debug "estimated allocations" baseline repeat10x repeat10x / baseline
    @test repeat10x < 2 * baseline
end

end  # module


text:
using WordCloud, Random

Random.seed!(19)
words = ["Machine Learning", "Artificial Intelligence", "Neural Networks",
         "Statistics","Optimization", "Statistical Learning", "Data Science",
         "Loss Function", "Supervised Learning", "ML", "RMSE", "MLE", "Probability",
         "Unsupervised Learning", "Reinforcement Learning", "AI", "SGD", "ConvNet",
         "Classification", "Regression", "Clustering", "Random Forest", "LASSO",
         "Regularization", "Decision Tree", "Object Recognition", "Industry 4.0",
         "Bias-Variance Tradeoff", "Cross Validation", "Training Set", "Test Set",
         "Gradient Boosting", "PCA", "Logistic Regression", "kNN", "k-means",
         "Algorithms"]
weights = [3.2; 1.8; 1.5; 1.8; .8*rand(length(words)-4) .+ .6]
wc = wordcloud(words, weights,
               density = .35, nepoch = 10^3, colors = :seaborn_dark,
               angles = -45:15:45, transparentcolor = nothing,
               retry = 5, mask = shape(box, 800, 450)) |> generate!
paint(wc, "wc.png")

max_stars_repo_name: baltisberger/MLCourse

max_issues_repo_name: baltisberger/MLCourse

max_forks_repo_name: baltisberger/MLCourse


text:
module PushVectors

export PushVector, finish!

mutable struct PushVector{T, V <: AbstractVector{T}} <: AbstractVector{T}
    "Vector used for storage."
    parent::V
    "Number of elements held by `parent`."
    len::Int
end

"""
    PushVector{T}(sizehint = 4)

Create a `PushVector` for elements typed `T`, with no initial elements. `sizehint`
determines the initial allocated size.
"""
function PushVector{T}(sizehint::Integer = 4) where {T}
    sizehint ≥ 0 || throw(DomainError(sizehint, "Invalid initial size."))
    PushVector(Vector{T}(undef, sizehint), 0)
end

@inline Base.length(v::PushVector) = v.len

@inline Base.size(v::PushVector) = (v.len, )

function Base.sizehint!(v::PushVector, n)
    if length(v.parent) < n || n ≥ v.len
        resize!(v.parent, n)
    end
    nothing
end

@inline function Base.getindex(v::PushVector, i)
    @boundscheck checkbounds(v, i)
    @inbounds v.parent[i]
end

@inline function Base.setindex!(v::PushVector, x, i)
    @boundscheck checkbounds(v, i)
    @inbounds v.parent[i] = x
end

function Base.push!(v::PushVector, x)
    v.len += 1
    if v.len > length(v.parent)
        resize!(v.parent, v.len * 2)
    end
    v.parent[v.len] = x
    v
end

Base.empty!(v::PushVector) = (v.len = 0; v)

function Base.append!(v::PushVector, xs)
    ι_xs = eachindex(xs)        # allow generalized indexing
    l = length(ι_xs)
    if l ≤ 4
        for x in xs
            push!(v, x)
        end
    else
        L = l + v.len
        if length(v.parent) < L
            resize!(v.parent, nextpow(2, nextpow(2, L)))
        end
        @inbounds copyto!(v.parent, v.len + 1, xs, first(ι_xs), l)
        v.len += l
    end
    v
end

"""
    finish!(v)

Shrink the buffer `v` to its current content and return that vector.

!!! NOTE
    Consequences are undefined if you modify `v` after this.
"""
finish!(v::PushVector) = resize!(v.parent, v.len)

end # module


text:
using Solc

solcCompiler()


text:
using Documenter, Tracking, GNSSSignals

makedocs(
    sitename="Tracking.jl",
    format = Documenter.HTML(prettyurls = false),
    modules = [Tracking, GNSSSignals],
    pages = [
        "index.md",
        "track.md",
        "tracking_state.md",
        "tracking_results.md",
        "loop_filter.md",
        "correlator.md",
        "cn0_estimator.md"
        ]
)

deploydocs(
    repo = "github.com/JuliaGNSS/Tracking.jl.git",
)


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains tests for utils

############################################################################################################################

@testset "Type Aliases" begin
   @test VertexID == Int
   @test EdgeID == Pair{Int,Int}
   @test VertexList == AbstractVector{Int}
   @test EdgeList == AbstractVector{Pair{Int,Int}}

   @test isa(EdgeID[1=>2, 3=>4], EdgeList)
end


text:
module EllipsoidSimilarity

using PDMats, NLopt
import Distances: evaluate, SqMahalanobis

export Ellipsoid,
       Compound,
       TransformationEnergy,
       similarity,
       compound_similarity,
       te_similarity

"""An Ellipsoid is defined as the level set
(x - m)' A (x - m) = 1

where m is a px1 vector and A is a pxp positive definite matrix.
"""
struct Ellipsoid{T<:AbstractFloat}
    m::Vector{T}
    A::AbstractPDMat{T}
end

function Ellipsoid{T}(m2::Vector{T}, A2::AbstractMatrix{T}) where {T<:AbstractFloat}
    local A_pd
    try
        A_pd = PDMat(A2)
    catch e
        if isa(e, Base.LinAlg.PosDefException)
            error("Matrix must be postivie definite")
        else
            @show A2
            throw(e)
        end
    end

    if length(m2) != A_pd.dim
        error("Center vector dimension and matrix dimension do not match")
    end

    Ellipsoid(m2,A_pd)
end

abstract type EllipsoidSimilarityMeasure <: Real end

"""Compound similarity

Product of three exponential factors accounting for location,
orientation, and shape of the ellipsoids. It satisfies the properties
(given in the paper) of a strong similarity measure.

Ref: Moshtaghi, M., et al., "Clustering ellipses for anomaly
detection". Pattern Recognition 44 (2011) pp. 55-69.
"""
abstract type Compound <: EllipsoidSimilarityMeasure end

"""Transformation energy similarity

Product of three exponential factors accounting for location,
orientation, and shape of the ellipsoids. It satisfies the properties
(given in the paper) of a strong similarity measure.

Ref: Moshtaghi, M., et al., "Clustering ellipses for anomaly
detection". Pattern Recognition 44 (2011) pp. 55-69.
"""
abstract type TransformationEnergy <: EllipsoidSimilarityMeasure end

"""Generalized Focal Distance

Average of the planar focal distances in all dimensions.

Ref: Moshtaghi, M., et al., "Clustering ellipses for anomaly
detection". Pattern Recognition 44 (2011) pp. 55-69.
"""
abstract type GeneralizedFocalDist <: EllipsoidSimilarityMeasure end

function similarity( meas::GeneralizedFocalDist, E1::Ellipsoid, E2::Ellipsoid )
    epts = focalsegments(E1)
    fpts = focalsegments(E2)

    dis = 0.0
    for i in 1:(E1.A.dim-1)
        dis += focaldist(epts[i], fpts[i])
    end
    return 1.0 - dis/(E1.A.dim-1)
end

function gf_similarity( m1::Vector{T}, A1::AbstractMatrix{T},
                        m2::Vector{T}, A2::AbstractMatrix{T} ) where {T<:AbstractFloat}
    similarity(GeneralizedFocalDist(),Ellipsoid(m1,A1),Ellipsoid(m2,A2))
end

function focaldist( epts::Tuple{Vector{Float64},Vector{Float64}},
                    fpts::Tuple{Vector{Float64},Vector{Float64}} )
    # d ordering:
    #   d(e1, f1)
    #   d(e1, f2)
    #   d(e2, f1)
    #   d(e2, f2)
    d = zeros(4)
    for i in 1:2
        d[2i-1] = norm(epts[i] - fpts[1])
        d[2i]   = norm(epts[i] - fpts[2])
    end

    return 0.25 * ( minimum(d[1:2]) +
                    minimum(d[3:4]) +
                    minimum([d[1],d[3]]) +
                    minimum([d[2],d[4]]) )
end

"Returns coordinates of all focal segment endpoints of E"
function focalsegments( E::Ellipsoid )
    segments = Array{Tuple{Vector{eltype(E.m)},Vector{eltype(E.m)}}}(E.A.dim - 1)

    (U,S,V) = svd(E.A.mat)
    si = sortperm(S) # maybe not necessary?

    for i in 1:(E.A.dim-1)
        α1 = S[si[i]]
        α2 = S[si[i+1]]
        u2 = V[:,si[i]]

        c = 0.5 * sqrt((α2 - α1)/(α2 * α1))
        segments[i] = (E.m + c*u2, E.m - c*u2)
    end
    return segments
end

function similarity( meas::TransformationEnergy, E1::Ellipsoid, E2::Ellipsoid; approx = false )
    S1, R1 = scale_and_rot_matrix(E1.A)
    S2, R2 = scale_and_rot_matrix(E2.A)
    M_12 = S2*R2*inv(R1)*inv(S1)
    M_21 = S1*R1*inv(R2)*inv(S2)
    d_12 = S2*R2*(E2.m - E1.m)
    d_21 = S1*R1*(E1.m - E2.m)

    if !approx
        # solve the minimization problem posed in the paper
        dist_1 = -te_minimize(M_12, d_12)
        dist_2 = -te_minimize(M_21, d_21)

        return 1.0/maximum([dist_1,dist_2])
    else
        # approximation given by the authors at the end of Section 4.
        # find the maximum singular values for each of M_12 and M_21
        σ_12 = maximum(svd(M_12)[2])
        σ_21 = maximum(svd(M_21)[2])

        return 1.0 / maximum([σ_12 + norm(d_12),
                        σ_21 + norm(d_21)])
    end
end

function te_similarity( m1::Vector{T}, A1::AbstractMatrix{T},
                        m2::Vector{T}, A2::AbstractMatrix{T} ) where {T<:AbstractFloat}
    similarity(TransformationEnergy(),Ellipsoid(m1,A1),Ellipsoid(m2,A2))
end

function te_minimize( M::AbstractMatrix{T}, d::Vector{T} ) where {T<:AbstractFloat}
    N = length(d)
    @assert N == size(M)[1]
    f(x, grad) = -norm(M*x + d)
    opt = Opt(:LN_COBYLA, N)
    xtol_rel!(opt,1e-4)

    min_objective!(opt, f)

    h(x, grad) = norm(x) - 1.0 # should be unit vector
    equality_constraint!(opt, h, 1e-6)

    xi = zeros(N)
    xi[1] = 1.0
    (minf,minx,ret) = optimize(opt, xi)

    return minf
end

"""Compute the 'scale matrix' given an ellipse represented by `A`. This
is a diagonal matrix whose elements are the recipricol of the square
root of the eigenvalues of `A`.
"""
function scale_and_rot_matrix( A::AbstractPDMat{T} ) where {T<:AbstractFloat}
    (_, S, V) = svd(A.mat)
    return (sqrt.(inv(diagm(S))), V)
end

function similarity( meas::Compound, E1::Ellipsoid, E2::Ellipsoid, p = 2 )
    check_dims(E1,E2)

    # the authors state that using the Euclidean norm in
    # `location_similarity` leads to the measure being highly
    # sensitive to differences in the mean. Instead, they recommend
    # the use of a generalized statistical measure, the Mahalanobis
    # distance, which is implemented here. A general
    # `location_similarity` function for use with any norm is provided
    # below.
    g1 = location_similarity_mahal(E1.m, E1.A.mat, E2.m, E2.A.mat )

    g2 = orientation_similarity(E1.A.mat, E2.A.mat, p)

    g3 = shape_similarity(E1.A.mat, E2.A.mat, p)

    return g1*g2*g3
end

function compound_similarity( m1::Vector{T}, A1::AbstractMatrix{T},
                            m2::Vector{T}, A2::AbstractMatrix{T} ) where {T<:AbstractFloat}
    similarity(Compound(),Ellipsoid(m1,A1),Ellipsoid(m2,A2))
end

function check_dims(E1::Ellipsoid, E2::Ellipsoid)
    @assert E1.A.dim == E2.A.dim
end

function location_similarity( m1::Vector{T}, m2::Vector{T}, p ) where {T<:AbstractFloat}
    return exp(-norm(m1 - m2, p))
end

function location_similarity_mahal( m1::Vector{T}, A1::AbstractMatrix{T},
                                    m2::Vector{T}, A2::AbstractMatrix{T} ) where {T<:AbstractFloat}
    return exp(-evaluate(SqMahalanobis(inv(A1+A2)),m1,m2))
end

function orientation_similarity(A1::AbstractMatrix{T}, A2::AbstractMatrix{T}, p ) where {T<:AbstractFloat}
    R1 = eig(A1)[2]
    R2 = eig(A2)[2]

    θ = acos.(clamp.(diag(R1'*R2),-1,1))
    return exp(-norm(sin.(θ), p))
end

function shape_similarity(A1::AbstractMatrix{T}, A2::AbstractMatrix{T}, p ) where {T<:AbstractFloat}
    α = sort!(eig(A1)[1])
    β = sort!(eig(A2)[1])

    α_star = 1 ./ sqrt.(α)
    β_star = 1 ./ sqrt.(β)

    return exp(-norm(α_star - β_star, p))
end

end


text:
using Zygote

f(x) = sin(cos(x))
ir = Zygote.@code_ir f(1.0)
display(ir)


text:
using SparseArrays, LinearAlgebra

println("This file demos some of the mathematical operations on Associative Array")

A = ReadCSV(joinpath(Base.source_dir(),"A.csv"))
A = logical(A)
printFull(A)

# Sum down rows and across columns
print(sum(A,1))
print(sum(A,2))

# Compute a simple join
Aa = A[:,"a,"]
Ab = A[:,"b,"]
Aab = nocol(Aa) & nocol(Ab)

# Compute a histogram (facets) of other columns that are in rows with both a and b
F =  ( Aab )' * A;
printFull(F)

# Compute normalized histogram
Fn = F ./ sum(A,1)
printFull(Fn)

# Compute correlation
AtA = sqIn(A)
d = diag(getadj(AtA))
AtA = putAdj(AtA,getadj(AtA) - sparse(diagm(0 => d)))
printFull(AtA)


text:
using Plots, LaTeXStrings, Statistics, LinearAlgebra

u = rand((0,1), 100)
rand(u)
w₀ = rand()
w = []
while true
    uᵐ = rand(u)
    if uᵐ == 1 and dot()


text:
using PyCall
using Conda

Conda.add("pip")
Conda.add("seaborn")
Conda.add("scikit-learn")
Conda.add("qt")
Conda.add("graphviz")
Conda.add("pydot")

pip = joinpath(Conda.SCRIPTDIR, "pip")
run(`$pip install --no-deps --force-reinstall mpldatacursor`)
run(`$pip install --no-deps --force-reinstall SAlib`)
run(`$pip install --no-deps --force-reinstall git+https://github.com/Project-Platypus/PRIM.git\#egg=prim`)
run(`$pip install --no-deps --force-reinstall git+https://github.com/Project-Platypus/Platypus.git\#egg=platypus-opt`)
run(`$pip install --no-deps --force-reinstall git+https://github.com/Project-Platypus/Rhodium.git\#egg=rhodium`)


text:
@testset "Treewidth deletion tests" begin
    # A simple graph to test on.
    N = 10
    G = lg.SimpleGraph(N)
    for i = 1:N
        for j = i+1:N
            lg.add_edge!(G, i, j)
        end
    end

    for i = 1:10
        lg.add_vertex!(G)
        lg.add_edge!(G, N+i-1, N+i)
    end

    G = LabeledGraph(G)
    tw, π̄ = min_fill(G)

    # check if the treewidth is reduced by the correct amount and number of vertices removed
    # is correct.
    Ḡ, μ = greedy_treewidth_deletion(G, 5; score_function=:degree)
    modified_tw, modifued_π̄ = min_fill(Ḡ)
    @test modified_tw == tw - 5
    @test nv(Ḡ) == 15
    @test length(μ) == 5

    Ḡ, μ, π̃s, τs = greedy_treewidth_deletion(G, 5;
                                             score_function=:direct_treewidth, elim_order=π̄)
    @test τs == [tw-1, tw-2, tw-3, tw-4, tw-5]
    @test nv(Ḡ) == 15
    @test length(μ) == 5
    @test length(π̃s[end]) == nv(Ḡ)

    # Test the tree trimming method for selecting vertices to delete.
    Ḡ, μ, π̃s, τs = greedy_treewidth_deletion(G, 5; elim_order=π̄)
    @test τs == [tw-1, tw-2, tw-3, tw-4, tw-5]
    @test nv(Ḡ) == 15
    @test length(μ) == 5
    @test length(π̃s[end]) == nv(Ḡ)
end


text:
# ------------------------------------------
#   Read ACE fit parameters from file

using ACE1pack

export  fill_defaults!, parse_ace_basis_keys

"""
Recursively updates nested dictionaries with default parameters.
One only needs to suply the mandatory parameters, for example
````
minimal_params = Dict(
    "data" => Dict(
        "fname" => "something"),
    "basis" => Dict(
        "ace" => Dict(
            "species" => "something",
            "N" => 1,
            "maxdeg" => 1
            ),
        "pair" => Dict(
            "species" => "sth",
            "maxdeg" => 1
            ),),
    "solver" => Dict(
        "type" => "rrqr"),
    "e0" => "something")
```
"""
function fill_defaults!(params::Dict; param_key = "fit_params")
    # Go through the nested dictionaries filling in the default values
    params = _fill_default(params, param_key)
    for (key, val) in params
        if key == "basis"
            for (basis_name, basis_params) in params[key]
                if basis_params["type"] == "ace"
                    basis_params = parse_ace_basis_keys(basis_params)
                end
                params[key][basis_name] = fill_defaults!(basis_params; param_key=key)
            end
        elseif key == "transforms"
            params[key] = _fill_transforms(val)
        elseif val isa Dict && haskey(_dict_constructors, key)
            params[key] = fill_defaults!(val; param_key = key)
        end
    end
    return params
end

_fill_transforms(params) = Dict([key => fill_defaults!(val, param_key="transform") for (key, val) in params])
_fill_default(d::Dict, key) = _dict_constructors[key](;_makesymbol(d)...)
_makesymbol(p::Pair) = (Symbol(p.first) => (p.second))
_makesymbol(D::Dict) = Dict(_makesymbol.([D...])...)

_dict_constructors = Dict(
    "fit_params" => ACE1pack.fit_params,
    "data" => ACE1pack.data_params,
    "solver" => ACE1pack.solver_params,
    "basis" => ACE1pack.basis_params,
    "ace" => ACE1pack.ace_basis_params,
    "pair" => ACE1pack.pair_basis_params,
    "radial" => ACE1pack.radial_basis_params,
    "transform" => ACE1pack.transform_params,
    "degree" => ACE1pack.degree_params,
    "P" => ACE1pack.regularizer_params)

"""
`("C", "C")`-type tuples are saved to and read back in as `"(\"C\", \"C\")"` to .json.
It's slightly easier to save these to .json or .yaml as `"(C, C)"`.
This function converts `"(C, C)"`-type strings back to parameter-friendly `("C", "C")`.

A bit ugly, but to change the dictionary keys (and so dictionary type) one needs to construct
it from scratch, at least as far as I can tell.
"""
function parse_ace_basis_keys(ace_basis::Dict)
    basis_out = Dict()
    for (key, val) in ace_basis
        if val isa Dict && key in ["transform", "degree"]
            sub_dict = Dict()
            for (sub_key, sub_val) in val
                if sub_key in ["transforms", "Dn", "Dl", "Dd"]
                    sub_dict[sub_key] = _parse_keys(sub_val)
                elseif sub_key == "cutoffs"
                    sub_dict[sub_key] = _parse_keys_and_vals(sub_val)
                else
                    sub_dict[sub_key] = sub_val
                end
            end
            basis_out[key] = sub_dict
        else
            basis_out[key] = val
        end
    end
    return basis_out
end

_parse_keys_and_vals(dict::Dict) = _parse_vals(_parse_keys(dict))

_parse_vals(dict::Dict) =
    Dict(_is_string_of_tuple(val) ? key => _parse_string_into_tuple(val) : key => val for (key, val) in dict)

_parse_keys(dict::Dict) =
    Dict(_is_string_of_tuple(key) ? _parse_string_into_tuple(key) => val : key => val for (key, val) in dict)

function _is_string_of_tuple(string)
    if !(typeof(string) <: AbstractString)
        return false
    end
    if string[1] == '(' && string[end] == ')'
        return true
    end
    return false
end

function _parse_entry(input::AbstractString)
    for T in [Int, Float32]
        try
            entry = parse(T, input)
            return entry
        catch
        end
    end
    return input
end

function _parse_string_into_tuple(input::AbstractString)
    input = strip(input, ['(', ')'])
    entries = [strip(s) for s in split(input, ",")]
    return Tuple(_parse_entry(String(entry)) for entry in entries)
end


text:
module ProjectEulerTasks

import HTTP

function taskmetas()
    map(s -> split(s, "##"),
        split(HTTP.get("https://projecteuler.net/minimal=problems").body |> String, "\r\n")
    )[2:end-1] .|> s -> TaskMeta(parse(Int, s[1]), s[2], s[3], s[4], parse(Int, s[5]))
end

struct TaskMeta
    index::Int
    title::String
    created
    updated
    solved::Int
end

link(t::TaskMeta) = "https://projecteuler.net/problem=$(t.index)"
link(i::Integer) = "https://projecteuler.net/problem=$i"

description(t::TaskMeta) = HTTP.get("https://projecteuler.net/minimal=$(t.index)").body |> String |> HTML
description(i::Integer) = HTTP.get("https://projecteuler.net/minimal=$i").body |> String |> HTML

function Base.show(io::IO, ::MIME"text/html", t::TaskMeta)
    print(io, """
Task $(t.index): $(t.title)

    $(description(t).content)
    link""")
end

function Base.show(io::IO, ::MIME"text/plain", t::TaskMeta)
    print(io, """Task $(t.index): $(t.title)
    $(description(t).content)
    link at: $(link(t))""")
end

export taskmetas, TaskMeta, link, description

end

max_stars_repo_name: ludwig-austermann/ProjectEulerTasks

max_issues_repo_name: ludwig-austermann/ProjectEulerTasks

max_forks_repo_name: ludwig-austermann/ProjectEulerTasks


text:
# a 2D triangulation example
# this example uses a factor graph to represent the data
# Factor graphs can be used to construct either fixed lag smoothing or full blown SLAM solutions.

# Load required libraries
using Caesar, Distributions
using RoMEPlotting

## Environment setup -- in Juno you can run cells with Alt+Enter

# true pose position (assume hidden)  - [X Y θ]
p0 = [1.0;1.0;0.0]

# true landmark positions
l0 = [0.0;1.0]
l1 = [1.0;0.0]
l2 = [2.0;1.0]

# noise level for Gaussian assumption -- see https://github.com/JuliaStats/Distributions.jl/issues/584
lcov = Matrix(Diagonal([0.01; 0.01].^2))

# measurements
m0 = pi/2
m1 = -pi
m2 = -pi/2

mstd = 0.05

## build the factor graph

fg = initfg()

# add two beacons
addVariable!(fg, :l0, Point2, labels=["BEACON"])
addFactor!(fg, [:l0], Prior(MvNormal(l0, lcov)))

addVariable!(fg, :l1, Point2, labels=["BEACON"])
addFactor!(fg, [:l1], Prior(MvNormal(l1, lcov)))

# addVariable!(fg, :l2, Point2, labels=["BEACON"])
# addFactor!(fg, [:l2], Prior(MvNormal(l2, lcov)))

# add unknown pose location
addVariable!(fg, :x0, Pose2, labels=["POSE"])

# add two bearing only measurements
addFactor!(fg, [:x0; :l0], Pose2Point2Bearing(Normal(m0, mstd)))
addFactor!(fg, [:x0; :l1], Pose2Point2Bearing(Normal(m1, mstd)))
# addFactor!(fg, [:x0; :l2], Pose2Point2Bearing(Normal(m2, mstd)))

# solve the system -- remember first runs are slow (just in time compile code)
[batchSolve!(fg) for i in 1:5]

# getKDEMean(getVertKDE(fg, :x0))

## Plot values to see result

plotPose(fg, :x0, axis=[-0.5 2.0; -0.5 1.5])

## more debugging

plotKDE(fg, [:l0; :l1],dims=[1;2],levels=3 ,axis=[0 2.0; 0.0 1.5])

stuff = localProduct(fg, :x0)

plotKDE(stuff[2],dims=[1;2],levels=3, c=["red","green","blue"])

pcm = [marginal(stuff[2][1],[3]);marginal(stuff[2][2],[3])]
AMP.plotKDECircular(pcm)

np = manifoldProduct(stuff[2], Pose2().manifolds)

AMP.plotKDECircular(marginal(np,[3]))

## more direct calculation

pts0 = approxConv(fg, :x0l0f1, :x0)
X0 = manikde!(pts0, Pose2().manifolds)

plotPose(Pose2(), [X0])

pts1 = approxConv(fg, :x0l1f1, :x0)
X1 = manikde!(pts1, Pose2().manifolds)

plotPose(Pose2(), [X1])

## Major error here

ptsl0 = approxConv(fg, :x0l0f1, :l0)
L0 = manikde!(ptsl0, Point2().manifolds)
plotKDE(L0, levels=3)

stuff = IIF.localProduct(fg, :l0)

plotKDE([stuff[2];stuff[1]])

## more debugging below

pts0_0 = approxConv(fg, :x0l0f1, :l0)
L00 = manikde!(pts0_0, Point2().manifolds)

plotKDE(L00, levels=3, axis=[-0.5 2; -0.5 1.5])

pts0_1 = approxConv(fg, :x0l1f1, :l1)
L01 = manikde!(pts0_1, Point2().manifolds)

plotKDE(L01, levels=3, axis=[-0.5 2; -0.5 1.5])

#

tree = batchSolve!(fg, drawpdf=true, show=true)

spyCliqMat(tree,:l0)

writeGraphPdf(fg, show=true)

stuff = IIF.treeProductUp(fg, tree, :l0, :l0)

Gadfly.plot(x=stuff[1][1,:], y=stuff[1][2,:], Geom.histogram2d)

#


text:
# Thanks again to Jonathan Paulson
# whose video https://youtu.be/d25r5GZa4us
# helped figure out a faster way to accomplish
# this puzzle

using StatsBase

function parse_input(input)
    nr = length(input)
    ns = length(input[1])
    input = reshape(split(join(input), ""), (ns, nr))

    return input
end

function fill_seats(input; part_one = true)
    (ns, nr) = size(input)
    while true
        output = deepcopy(input)
        change = false

        # Loop over all rows and seats
        for r in 1:nr
            for s in 1:ns

                # Find # of occurences
                nocc = 0
                for dr in [-1,0,1]
                    for ds in [-1,0,1]
                        if !(dr == 0 && ds == 0)
                            rr = r + dr
                            ss = s + ds

                            while 1<=(rr)<=nr && 1<=(ss)<=ns && input[ss,rr] == "." && !part_one
                                rr += dr
                                ss += ds
                            end
                            if 1<=(rr)<=nr && 1<=(ss)<=ns && input[ss,rr] == "#"
                                nocc += 1
                            end
                        end
                    end
                end

                # Change output based on occurences
                if input[s,r] == "L" && nocc == 0
                    output[s,r] = "#"
                    change = true
                elseif input[s,r] == "#" && nocc >= (if part_one 4 else 5 end)
                    output[s,r] = "L"
                    change = true
                end
            end
        end

        if change == false return output end
        input = deepcopy(output)
    end

end

function total_empty(input)
    return countmap(input)["#"]
end

seats = readlines("seats.txt")

final_seats = fill_seats(parse_input(seats))
println("Part 1: ", total_empty(final_seats))

final_seats = fill_seats(parse_input(seats), part_one=false)
println("Part 2: ", total_empty(final_seats))


text:
using Documenter, QuadraticFormsMGHyp, DocThemeIndigo
indigo = DocThemeIndigo.install(QuadraticFormsMGHyp)
makedocs(;
    modules=[QuadraticFormsMGHyp],
    format=Documenter.HTML(assets=String[indigo]),
    pages=[
        "Home" => "index.md",
    ],
    sitename="QuadraticFormsMGHyp.jl",
    authors="Simon Broda",
)

deploydocs(;
    repo="github.com/s-broda/QuadraticFormsMGHyp.jl",
)


text:
module Maracas
include("test.jl")
export @test, @test_throws, @test_broken, @test_skip, @test_warn, @test_nowarn
export @testset
export @describe, @it, @unit, @skip, MARACAS_SETTING
export MaracasTestSet, DescribeTestSet, SpecTestSet, TestTestSet
export set_test_style, set_title_style, set_spec_style, set_error_color, set_warn_color, set_pass_color, set_info_color
const MARACAS_SETTING = Dict(
    :error => Symbol(get(ENV, "MARACAS_ERROR", :red)),
    :warn =>  Symbol(get(ENV, "MARACAS_WARN", :yellow)),
    :pass =>  Symbol(get(ENV, "MARACAS_PASS", :green)),
    :info =>  Symbol(get(ENV, "MARACAS_INFO", :blue)),
    :default => Base.text_colors[:normal],
    :bold => Base.text_colors[Symbol(get(ENV, "MARACAS_BOLD", :bold))],
    :title_length => 80,
    :test =>  "",
    :title =>  "",
    :spec =>  "",
)
MARACAS_SETTING[:test] =  get(ENV, "MARACAS_TEST", string(Base.text_colors[:blue], MARACAS_SETTING[:bold]))
MARACAS_SETTING[:title] =  get(ENV, "MARACAS_TITLE", string(Base.text_colors[:magenta], MARACAS_SETTING[:bold]))
MARACAS_SETTING[:spec] =  get(ENV, "MARACAS_SPEC", string(Base.text_colors[:cyan], MARACAS_SETTING[:bold]))

function set_text_style(key::Symbol, color::Symbol, style::Symbol=:bold)
    MARACAS_SETTING[key] = string(Base.text_colors[style], Base.text_colors[color])
end

const TextColor = Union{Symbol, UInt8}

set_test_style(color::TextColor, bold::Bool=true) = set_text_style(:test, color, bold ? :bold : :normal)
set_title_style(color::TextColor, bold::Bool=true) = set_text_style(:title, color, bold ? :bold : :normal)
set_spec_style(color::TextColor, bold::Bool=true) = set_text_style(:spec, color, bold ? :bold : :normal)
set_error_color(color::TextColor) = (MARACAS_SETTING[:error] = color)
set_warn_color(color::TextColor) = (MARACAS_SETTING[:warn] = color)
set_pass_color(color::TextColor) = (MARACAS_SETTING[:pass] = color)
set_info_color(color::TextColor) = (MARACAS_SETTING[:info] = color)

const MACRO_TYPES = Dict(
    "@describe" => DescribeTestSet,
    "@it" => SpecTestSet,
    "@unit" => TestTestSet,
)

function _check_args(macro_name, desc, tests)
    if !isa(desc, AbstractString) && !(isa(desc, Expr) && desc.head == :string)
        error("Unexpected argument $desc to $macro_name")
    end
    if !isa(tests, Expr) || tests.head != :block
        error("Expected begin/end block as argument to $macro_name")
    end
    return desc, tests
end
function genexpr(macro_name, desc, tests, source)
    _check_args(macro_name, desc, tests)
    ex = quote
        ts = $(MACRO_TYPES[string(macro_name)])($desc)
        while false; end
        Test.push_testset(ts)
        try
            $(esc(tests))
        catch err
            record(ts, Error(:nontest_error, :(), err, catch_backtrace(), $(QuoteNode(source))))
        end
        Test.pop_testset()
        finish(ts)
    end
    if !isempty(tests.args) &&  isa(tests.args[1], LineNumberNode)
        ex = Expr(:block, tests.args[1], ex)
    end
    return ex
end
macro describe(desc, tests)
    genexpr("@describe", desc, tests, __source__)
end

macro it(desc, tests)
    genexpr("@it", desc, tests, __source__)
end

macro unit(desc, tests)
    genexpr("@unit", desc, tests, __source__)
end

function extract_test_title(args)
    for arg in args
        if isa(arg, AbstractString) return arg end
    end
    return ""
end

macro skip(tests)
    if tests.head != :macrocall
        error("@skip must be followed by a testing macro (@describe, @it, @unit, @test)")
    end
    macro_name = string(first(tests.args))
    if startswith(macro_name, "@test")
        orig_ex = Expr(:inert, tests.args[2:end])
        testres = :(Broken(:skipped, $orig_ex))
        return :(record(get_testset(), $testres))
    end
    desc = extract_test_title(tests.args)
    ex = quote
        ts = $(MACRO_TYPES[macro_name])($desc)
        record(ts, Broken(:skipped, ts))
        finish(ts)
    end
    return ex
end

end


text:
module SpringCollab2020InstantFallingBlock

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/InstantFallingBlock" InstantFallingBlock(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight,
    tiletype::String="3", climbFall::Bool=true, behind::Bool=false)

const placements = Ahorn.PlacementDict(
    "Instant Falling Block (Spring Collab 2020)" => Ahorn.EntityPlacement(
        InstantFallingBlock,
        "rectangle",
        Dict{String, Any}(),
        Ahorn.tileEntityFinalizer
    ),
)

Ahorn.editingOptions(entity::InstantFallingBlock) = Dict{String, Any}(
    "tiletype" => Ahorn.tiletypeEditingOptions()
)

Ahorn.minimumSize(entity::InstantFallingBlock) = 8, 8
Ahorn.resizable(entity::InstantFallingBlock) = true, true

Ahorn.selection(entity::InstantFallingBlock) = Ahorn.getEntityRectangle(entity)

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::InstantFallingBlock, room::Maple.Room) = Ahorn.drawTileEntity(ctx, room, entity)

end


text:
include("measure.jl")
include("expectation.jl")
include("moment_matrix.jl")
include("atomic.jl")


text:
"""
Repository database of thermodynamical properties.

Exported methods are

- `martinez()::ComponentDatabase`: gas properties of substances;
- `brusilovsky_comp()::ComponentDatabase`: parameters of A I Brusilovskii equation of state for pure components;
- `brusilovsky_mix()::MixtureDatabase`: binary parameters of A I Brusilovskii equation of state for mixtures.
"""
module Data

using ..CubicEoSDatabase

export martinez, brusilovsky_comp, brusilovsky_mix

"Cache of a database `::T`. Stores arguments needed to initialize `T`."
mutable struct CachedDatabase{T,A,KW}
    value::Union{Nothing,T}
    calculated::Bool
    args::A     # `T(args...; kwargs...)` constructor arguments
    kwargs::KW  #
    function CachedDatabase(::Type{T}, args, kwargs) where {T}
        new{T,typeof(args),typeof(kwargs)}(nothing, false, args, kwargs)
    end
end

function value(x::CachedDatabase{T,A,KW}) where {T,A,KW}
    if !x.calculated
        x.value = T(x.args...; x.kwargs...)
        x.calculated = true
    end
    return x.value::T
end

"Root path of database source files."
const root = joinpath(abspath(@__DIR__), "../data")

const CachedMartinez = CachedDatabase(ComponentDatabase,
    (joinpath(root, "phys/comp/martinez.csv"),),
    (reference="Website (loaded in 2017): http://imartinez.etsiae.upm.es/~isidoro/dat1/eGAS.pdf", delim=',')
)

"""
    Data.martinez() -> ComponentDatabase

I Martinez gas properties.
"""
martinez() = value(CachedMartinez)

const brusilovsky_book = "A I Brusilovskii. Fazovie prevrasheniya pri razrabotke mestorozhdeniy nefti i gaza (Phase transitions during mining of petroleum and gas). 2002. Moscow, Graal’. ISBN: 5-94688-031-4. (in Russian)"
const CachedBrusilovskyComp = CachedDatabase(
    ComponentDatabase,
    (joinpath(root, "eos/comp/brusilovsky.csv"),),
    (reference="Table 4.14 / "*brusilovsky_book, delim=','),
)

"""
    Data.brusilovsky_comp() -> ComponentDatabase

A I Brusilovsky equation of state. Parameters for pure components.
"""
brusilovsky_comp() = value(CachedBrusilovskyComp)

const CachedBrusilovskyMixture = CachedDatabase(
    MixtureDatabase,
    (joinpath(root, "eos/mix/brusilovsky.csv"),),
    (
        reference="Table 4.22 (mistakes: gij->gij⋅10³, hij->hij⋅10⁵) / "*brusilovsky_book,
        delim=',',
    ),
)

"""
    Data.brusilovsky_mix() -> MixtureDatabase

A I Brusilovsky equation of state. Binary parameters for mixtures.
"""
brusilovsky_mix() = value(CachedBrusilovskyMixture)

const CachedBrusilovskyMixtureAdjusted = CachedDatabase(
    MixtureDatabase,
    (joinpath(root, "eos/mix/brusilovsky_adjusted.csv"),),
    (
        reference="See `Data.brusilovsky_mix()`, but missed values for up to C7+ are adjusted.",
        delim=',',
    ),
)

brusilovsky_mix_adjusted() = value(CachedBrusilovskyMixtureAdjusted)

end # module


text:
@testset "Fast Math" begin
    x = 3.0
    @test FastMath.invsqrt(x) ≈ 1.0 / sqrt(x) rtol = 1.0e-5
    x = 1.0 / 3.0
    @test FastMath.invsqrt(x) ≈ 1.0 / sqrt(x) rtol = 1.0e-5
    x = 3.0f0
    @test FastMath.invsqrt(x) ≈ 1.0 / sqrt(x) rtol = 1.0e-5
    x = 1.0f0 / 3.0f0
    @test FastMath.invsqrt(x) ≈ 1.0 / sqrt(x) rtol = 1.0e-5

    # semicircle -1<ω<1
#     function S(ω)
#     if -1.0 < ω < 1.0
#         return sqrt(1.0 - ω^2)
#     else
#         return 0.0
#     end
# end

#     result, err = FastMath.integrate(S, -1.0, 1.0, :cuhre)
#     @test abs(result - π / 2) < 3err
#     result, err = FastMath.integrate(S, -Inf, 1.0, :cuhre)
#     @test abs(result - π / 2) < 3err
#     result, err = FastMath.integrate(S, -1.0, Inf, :cuhre)
#     @test abs(result - π / 2) < 3err
#     result, err = FastMath.integrate(S, -Inf, Inf, :cuhre)
#     @test abs(result - π / 2) < 3err

    # k = MVector{3,Float64}([1.0, 2.0, 3.0])
    # q = MVector{3,Float64}([3.0, 1.0, 4.0])
    # println(FastMath.dot(k, q))
    # @test FastMath.dot(k, q) ≈ LinearAlgebra.dot(k, q)
    # @test FastMath.norm(k) ≈ LinearAlgebra.norm(k)
    # @test FastMath.squaredNorm(k) ≈ LinearAlgebra.dot(k, k)
end


text:
help_str(x; color = true) = sprint(print_cmd, x; context = :color => color)

# printing in expression
xprint_help(x; color = true) = :(print($(help_str(x; color = color))))
xprint_version(cmd::EntryCommand; color = true) =
    :(println($(sprint(show, cmd; context = :color => color))))

"""
    xerror(cmd, msg)

Create an expression that contains error message, automatically merge
interpolation expressions.
"""
function xerror(cmd, msg::String)
    msg = "Error: $msg, use -h or --help to check more detailed help info"
    return quote
        printstyled($msg; color=:red, bold=true)
        $(xprint_help(cmd))
        return 1
    end
end

function xerror(cmd, msg::Expr)
    msg.head === :string || throw(Meta.ParseError("expect string expression, got $msg"))

    ex = Expr(:string, "Error: ")
    for each in msg.args
        if each isa Number
            push!(ex.args, string(each))
        else
            push!(ex.args, each)
        end
    end
    push!(ex.args, ", use -h or --help to check more detailed help info")
    return quote
        printstyled($ex; color=:red, bold=true)
        $(xprint_help(cmd))
        return 1
    end
end

hasparameters(cmd::LeafCommand) = (!isempty(cmd.flags)) || (!isempty(cmd.options))
hassubcmds(cmd::EntryCommand) = hassubcmds(cmd.root)
hassubcmds(cmd::LeafCommand) = false
hassubcmds(cmd::NodeCommand) = true

# regexes
short_name(x) = string(first(cmd_name(x)))

regex_flag(x) = Regex("^--$(cmd_name(x))\$")
regex_short_flag(x) = Regex("^-$(short_name(x))\$")

regex_option(x) = Regex("^--$(cmd_name(x))=(.*)")
regex_short_option(x) = Regex("^-$(short_name(x))(.*)")

skip_type_convert(type) = (type in [Any, String, AbstractString])

function xparse_args(arg::Arg, index)
    if arg.vararg
        return xparse_varargs(arg.type, index)
    else
        return xparse(arg.type, :(ARGS[$index]))
    end
end

function xparse(type, str)
    if skip_type_convert(type)
        return :($str)
    else
        return :(convert($type, Meta.parse($str)))
    end
end

function xparse_varargs(type, index)
    if skip_type_convert(type)
        return :(ARGS[$index:end]...)
    else
        return :(map(x -> convert($type, Meta.parse(x)), ARGS[$index:end])...)
    end
end

"""
    xmatch(regexes, actions, str[, st = 1])

Generate a long ifelse expression that acts like a pattern
matching expression that match given regex list and do the
corresponding actions.
"""
function xmatch(cmd, regexes, actions, str, st = 1)
    if st <= length(regexes)
        m = gensym(:m)
        regex, action = regexes[st], actions[st]

        return quote
            $m = match($regex, $str)
            if $m === nothing
                $(xmatch(cmd, regexes, actions, str, st + 1))
            else
                $(action(m))
            end
        end
    else
        return xerror(cmd, :("unknown option: $($str)"))
    end
end

function all_required(args::Vector)
    return all(args) do x
        x.require == true
    end
end

all_required(cmd::LeafCommand) = all_required(cmd.args)

function nrequired_args(args::Vector)
    return count(args) do x
        x.require == true && x.vararg == false
    end
end

"""
    rm_lineinfo(ex)

Remove `LineNumberNode` in a given expression
"""
function rm_lineinfo(ex)
    if ex isa Expr
        args = []
        for each in ex.args
            if !(each isa LineNumberNode)
                push!(args, rm_lineinfo(each))
            end
        end
        return Expr(ex.head, args...)
    else
        return ex
    end
end

"""
    prettify(ex)

Prettify given expression, remove all `LineNumberNode` and
extra code blocks.
"""
prettify(x) = x

function prettify(ex::Expr)
    ex = rm_lineinfo(ex)
    ex.head === :block || return Expr(ex.head, map(prettify, ex.args)...)

    if any(ex.args) do x
        x isa Expr && x.head === :block
    end

        return prettify(eat_blocks(ex))
    end
    return Expr(ex.head, map(prettify, ex.args)...)
end

function eat_blocks(ex::Expr)
    ex.head === :block || return Expr(ex.head, map(prettify, ex.args)...)
    args = []
    for stmt in ex.args
        if stmt isa Expr && stmt.head === :block
            for each in stmt.args
                push!(args, prettify(each))
            end
        else
            push!(args, prettify(stmt))
        end
    end
    return Expr(:block, args...)
end

"""
    push!(args, item)

Push an item to the expression or list. Do nothing
if the item is `nothing`.
"""
pushmaybe!(args, item) = push!(args, item)
pushmaybe!(args, ::Nothing) = args
pushmaybe!(ex::Expr, item) = push!(ex.args, item)
pushmaybe!(ex::Expr, ::Nothing) = ex


text:
isdefined(Main, :LocalFilters) || include("../src/LocalFilters.jl")
isdefined(Main, :NaiveLocalFilters) || include("NaiveLocalFilters.jl")

module LocalFiltersTests

using Compat, Compat.Test
using LocalFilters
using LocalFilters: Neighborhood, RectangularBox, Kernel,
    axes, initialindex, finalindex, limits, cartesianregion, ball, coefs,
    ismmbox, strictfloor, USE_CARTESIAN_RANGE, _range

struct Empty{N} <: Neighborhood{N} end

# Selector for reference methods.
const REF = Val(:Base)

const ATOL = 0.0
const GTOL = 4*eps(Float64)

replicate(a, n::Integer) = ntuple(i->a, n)
compare(a::AbstractArray, b::AbstractArray) = maximum(abs(a - b))
samevalues(a::AbstractArray, b::AbstractArray) = minimum(a .== b)
identical(a::RectangularBox{N}, b::RectangularBox{N}) where {N} =
    (initialindex(a) === initialindex(b) &&
     finalindex(a) === finalindex(b))
identical(a::Kernel{Ta,N}, b::Kernel{Tb,N}) where {Ta,Tb,N} =
    nearlysame(a, b; atol=0, gtol=0)
identical(a::Neighborhood, b::Neighborhood) = false

function nearlysame(a::Kernel{Ta,N}, b::Kernel{Tb,N};
                    atol=ATOL, gtol=GTOL) where {Ta,Tb,N}
    initialindex(a) == initialindex(b) || return false
    finalindex(a) == finalindex(b) || return false
    if atol == 0 && gtol == 0
        for i in cartesianregion(a)
            a[i] == b[i] || return false
        end
    else
        for i in cartesianregion(a)
            ai, bi = a[i], b[i]
            abs(ai - bi) ≤ atol + gtol*max(abs(ai), abs(bi)) || return false
        end
    end
    return true
end

function similarvalues(A::AbstractArray{Ta,N}, B::AbstractArray{Tb,N};
                       atol=ATOL, gtol=GTOL) where {Ta,Tb,N}
    @assert axes(A) == axes(B)
    T = float(promote_type(Ta, Tb))
    local sd2::T = 0
    local sa2::T = 0
    local sb2::T = 0
    @inbounds for i in eachindex(A, B)
        a = T(A[i])
        b = T(B[i])
        sa2 += a*a
        sb2 += b*b
        sd2 += (a - b)^2
    end
    return sqrt(sd2) ≤ atol + gtol*sqrt(max(sa2, sb2))
end

function reversealldims(A::Array{T,N}) where {T,N}
    B = Array{T,N}(undef, size(A))
    len = length(A)
    off = len + 1
    @inbounds for i in 1:len
        B[off - i] = A[i]
    end
    return B
end

function checkindexing!(ker::Kernel{T,1}) where {T}
    Imin, Imax = limits(ker)
    i1min, i1max = Imin[1], Imax[1]
    tmp = zero(T)
    for i1 in i1min:i1max
        I = CartesianIndex(i1)
        val = ker[I]
        ker[i1] == val || return false
        ker[i1] = tmp
        ker[i1] == ker[I] == tmp || return false
        ker[I] = val
    end
    return true
end

function checkindexing!(ker::Kernel{T,2}) where {T}
    Imin, Imax = limits(ker)
    i1min, i1max = Imin[1], Imax[1]
    i2min, i2max = Imin[2], Imax[2]
    tmp = zero(T)
    for i2 in i2min:i2max,
        i1 in i1min:i1max
        I = CartesianIndex(i1,i2)
        val = ker[I]
        ker[i1,i2] == val || return false
        ker[i1,i2] = tmp
        ker[i1,i2] == ker[I] == tmp || return false
        ker[I] = val
    end
    return true
end

function checkindexing!(ker::Kernel{T,3}) where {T}
    Imin, Imax = limits(ker)
    i1min, i1max = Imin[1], Imax[1]
    i2min, i2max = Imin[2], Imax[2]
    i3min, i3max = Imin[3], Imax[3]
    tmp = zero(T)
    for i3 in i3min:i3max,
        i2 in i2min:i2max,
        i1 in i1min:i1max
        I = CartesianIndex(i1,i2,i3)
        val = ker[I]
        ker[i1,i2,i3] == val || return false
        ker[i1,i2,i3] = tmp
        ker[i1,i2,i3] == ker[I] == tmp || return false
        ker[I] = val
    end
    return true
end

f1(x) = 1 + x*x
f2(x) = x > 0.5

@testset "LocalFilters" begin

    @testset "Miscellaneous" begin
        for T in (Int, Float32)
            @test limits(T) === (typemin(T), typemax(T))
        end
        for dim in (11, 12)
            @test _range(Int8(dim)) === _range(dim)
        end
        @test_throws ArgumentError _range(0)
        @test strictfloor(Int, 5*(1 + eps(Float64))) == 5
        @test strictfloor(Int, 5.0) == 4
    end

    @testset "Neighborhoods" begin
        for (dims, rngs) in (((3,), (-1:1,)),
                             ((3, 4, 5), (-1:1, -2:1, -2:2)))
            N = length(dims)
            box = Neighborhood(dims)
            fse = strel(Bool, box) # flat structuring element
            Imin, Imax = limits(box)
            A = rand(dims...)
            msk = rand(Bool, dims)
            ker = Kernel(A)

            # Test limits(), initialindex() and finalindex().
            @test initialindex(CartesianIndices(rngs)) === Imin
            @test finalindex(CartesianIndices(rngs)) === Imax
            @test limits(CartesianIndices(rngs)) === (Imin, Imax)
            @static if USE_CARTESIAN_RANGE
                @test initialindex(CartesianRange(rngs)) === Imin
                @test finalindex(CartesianRange(rngs)) === Imax
                @test limits(CartesianRange(rngs)) === (Imin, Imax)
            end
            @test initialindex(box) === Imin
            @test finalindex(box) === Imax
            @test limits(box) === (Imin, Imax)
            @test initialindex(A) === one(CartesianIndex{N})
            @test finalindex(A) === CartesianIndex(size(A))
            @test limits(A) === (one(CartesianIndex{N}),
                                 CartesianIndex(size(A)))

            # Test cartesianregion().
            @static if USE_CARTESIAN_RANGE
                region = CartesianRange(rngs)
                @test cartesianregion(Imin,Imax) === region
                @test cartesianregion(CartesianRange(rngs)) === region
                @test cartesianregion(CartesianIndices(rngs)) === region
                @test cartesianregion(box) === region
                @test cartesianregion(A) === CartesianRange(size(A))
            else
                region = CartesianIndices(rngs)
                @test cartesianregion(Imin,Imax) === region
                @test cartesianregion(CartesianIndices(rngs)) === region
                @test cartesianregion(box) === region
                @test cartesianregion(A) === CartesianIndices(A)
            end

            # Neighborhood constructors.
            @test Neighborhood(box) === box
            @test Neighborhood(ker) === ker
            @test Neighborhood(A) === ker
            @test Neighborhood(dims...) === box
            @test Neighborhood(rngs) === box
            @test Neighborhood(rngs...) === box
            @test Neighborhood(CartesianIndices(rngs)) === box
            @static if USE_CARTESIAN_RANGE
                @test Neighborhood(CartesianRange(rngs)) === box
                @test Neighborhood(CartesianRange(Imin, Imax)) === box
            end
            @test convert(Neighborhood, box) === box
            @test convert(Neighborhood, ker) === ker
            @test convert(Neighborhood, rngs) === box
            @test convert(Neighborhood, dims) === box

            # RectangularBox constructors.
            @test RectangularBox(box) === box
            @test RectangularBox(dims) === box
            @test RectangularBox(dims...) === box
            @test RectangularBox(rngs) === box
            @test RectangularBox(rngs...) === box
            @static if USE_CARTESIAN_RANGE
                @test CartesianRange(box) === CartesianRange(rngs)
                @test CartesianRange(box) === CartesianRange(Imin, Imax)
            end
            @test convert(RectangularBox, box) === box
            @test convert(RectangularBox, rngs) === box
            @test convert(RectangularBox, dims) === box
            dim = dims[end]
            rng = rngs[end]
            @test RectangularBox{N}(dim) === RectangularBox(ntuple(d -> dim, N))
            @test RectangularBox{N}(rng) === RectangularBox(ntuple(d -> rng, N))

            # Conversions RectangularBox <-> CartesianIndices.
            @test CartesianIndices(box) === CartesianIndices(rngs)
            @test RectangularBox(CartesianIndices(rngs)) === box
            @test convert(RectangularBox, CartesianIndices(rngs)) === box
            @test convert(RectangularBox{N}, CartesianIndices(rngs)) === box
            @test CartesianIndices(box) === CartesianIndices(rngs)
            @test convert(CartesianIndices, box) === CartesianIndices(rngs)
            @test convert(CartesianIndices{N}, box) === CartesianIndices(rngs)
            @static if USE_CARTESIAN_RANGE
                # Conversion RectangularBox <-> CartesianRange.
                @test RectangularBox(CartesianRange(rngs)) === box
                @test convert(RectangularBox, CartesianRange(rngs)) === box
                @test convert(RectangularBox{N}, CartesianRange(rngs)) === box
                @test CartesianRange(box) === CartesianRange(rngs)
                @test convert(CartesianRange, box) === CartesianRange(rngs)
                @test convert(CartesianRange{CartesianIndex{N}}, box) ===
                    CartesianRange(rngs)
            end

            # Kernel constructors.
            @test Kernel(ker) === ker
            @test identical(Kernel(box), Kernel(ones(Bool, dims)))
            @test identical(Kernel{Bool}(box), Kernel(ones(Bool, dims)))
            @test Kernel(A, initialindex(ker)) === ker
            @test Kernel(A, rngs) === ker
            @test Kernel(A, rngs...) === ker
            @test Kernel(A, CartesianIndices(ker)) === ker
            @static if USE_CARTESIAN_RANGE
                @test Kernel(A, CartesianRange(ker)) === ker
            end
            off = initialindex(A) - initialindex(ker)
            @test identical(Kernel{eltype(A)}(i -> f1(A[off + i]),
                                              CartesianIndices(ker)),
                            Kernel(map(f1, A)))
            @test identical(Kernel(i -> f1(A[off + i]),
                                   CartesianIndices(ker)), Kernel(map(f1, A)))
            @test identical(Kernel{Bool}(i -> f2(A[off + i]),
                                         CartesianIndices(ker)),
                            Kernel(map(f2, A)))
            @test identical(Kernel(i -> f2(A[off + i]),
                                   CartesianIndices(ker)), Kernel(map(f2, A)))
            @test identical(convert(Kernel, box), fse)
            @test identical(Kernel((0,-Inf), msk), strel(Float64, Kernel(msk)))
            @test identical(strel(Bool, box), Kernel(box))
            @test convert(Kernel, fse) === fse
            @test convert(Kernel, ker) === ker

            @test nearlysame(Float32(ker), ker; atol=0, gtol=eps(Float32))
            @test nearlysame(Float32(ker), Float64(ker);
                             atol=0, gtol=eps(Float32))
            @test checkindexing!(ker)

            # Conversion Neighborhood <-> CartesianIndices.
            @test Neighborhood(CartesianIndices(rngs)) === box
            @test convert(Neighborhood, CartesianIndices(rngs)) === box
            @test convert(Neighborhood{N}, CartesianIndices(rngs)) === box
            @test CartesianIndices(box) === CartesianIndices(rngs)
            @test convert(CartesianIndices, box) === CartesianIndices(rngs)
            @test convert(CartesianIndices{N}, box) === CartesianIndices(rngs)
            @static if USE_CARTESIAN_RANGE
                # Conversion Neighborhood <-> CartesianRange
                @test Neighborhood(CartesianRange(rngs)) === box
                @test convert(Neighborhood, CartesianRange(rngs)) === box
                @test convert(Neighborhood{N}, CartesianRange(rngs)) === box
                @test CartesianRange(box) === CartesianRange(rngs)
                @test convert(CartesianRange, box) === CartesianRange(rngs)
                @test convert(CartesianRange{CartesianIndex{N}}, box) ===
                    CartesianRange(rngs)
            end

            # Other basic methods.
            @test length(box) === length(CartesianIndices(rngs))
            @test size(box) === size(CartesianIndices(rngs))
            @test size(box) === ntuple(d -> size(box, d), N)
            @test axes(box) === ntuple(d -> axes(box, d), N)
            @test length(ker) === length(A)
            @test size(ker) === size(A)
            @test size(ker) === ntuple(d -> size(ker, d), N)
            @test axes(ker) === ntuple(d -> axes(ker, d), N)
            @test ismmbox(fse) == true
            @test ismmbox(box) == true
            boolker = Neighborhood(ones(Bool,dims))
            @test ismmbox(boolker) == true
            boolker[zero(CartesianIndex{N})] = false
            @test ismmbox(boolker) == false
            @test ismmbox(Empty{N}()) == false

            # Test reverse().
            revbox = reverse(box)
            revker = reverse(ker)
            @test initialindex(revbox) === -finalindex(box)
            @test finalindex(revbox) === -initialindex(box)
            @test initialindex(revker) === -finalindex(ker)
            @test finalindex(revker) === -initialindex(ker)
            @test samevalues(coefs(revker), reversealldims(coefs(ker)))
        end
    end

    for (arrdims, boxdims) in (((341,),    (3,)),
                               ((62,81),   (5,3)),
                               ((23,28,27),(3,7,5)))
        rank = length(arrdims)
        @testset "$(rank)D test" begin
            a = rand(arrdims...)
            box = RectangularBox(boxdims)
            mask = Kernel(box)
            fse = strel(eltype(a), mask) # flat structuring element
            kern = Kernel((1.0, 0.0), mask) # kernel for convolution
            @testset "erode" begin
                result = erode(REF, a, box)
                @test samevalues(erode(a, box), result)
                @test samevalues(erode(a, mask), result)
                @test samevalues(erode(a, fse), result)
                @test samevalues(erode!(copy(a), box), result)
            end
            @testset "dilate" begin
                result = dilate(REF, a, box)
                @test samevalues(dilate(a, box), result)
                @test samevalues(dilate(a, mask), result)
                @test samevalues(dilate(a, fse), result)
                @test samevalues(dilate!(copy(a), box), result)
            end
            @testset "closing" begin
                result = closing(REF, a, box)
                @test samevalues(closing(a, box), result)
                @test samevalues(closing(a, mask), result)
                @test samevalues(closing(a, fse), result)
            end
            @testset "opening" begin
                result = opening(REF, a, box)
                @test samevalues(opening(a, box), result)
                @test samevalues(opening(a, mask), result)
                @test samevalues(opening(a, fse), result)
            end
            @testset "bottom-hat" begin
                result = bottom_hat(REF, a, box)
                @test samevalues(bottom_hat(a, box), result)
                @test samevalues(bottom_hat(a, mask), result)
                @test samevalues(bottom_hat(a, fse), result)
                result = bottom_hat(REF, a, box, 3)
                @test samevalues(bottom_hat(a, box, 3), result)
                @test samevalues(bottom_hat(a, mask, 3), result)
                @test samevalues(bottom_hat(a, fse, 3), result)
            end
            @testset "top-hat" begin
                result = top_hat(REF, a, box)
                @test samevalues(top_hat(a, box), result)
                @test samevalues(top_hat(a, mask), result)
                @test samevalues(top_hat(a, fse), result)
                result = top_hat(REF, a, box, 3)
                @test samevalues(top_hat(a, box,  3), result)
                @test samevalues(top_hat(a, mask, 3), result)
                @test samevalues(top_hat(a, fse, 3), result)
            end
            @testset "localextrema" begin
                e0, d0 = erode(REF, a, box), dilate(REF, a, box)
                e1, d1 = localextrema(REF, a, box)
                @test samevalues(e0, e1) && samevalues(d0, d1)
                e1, d1 = localextrema(a, box)
                @test samevalues(e0, e1) && samevalues(d0, d1)
                e1, d1 = localextrema!(similar(a), similar(a), a, box)
                @test samevalues(e0, e1) && samevalues(d0, d1)
                e1, d1 = localextrema(a, mask)
                @test samevalues(e0, e1) && samevalues(d0, d1)
                e1, d1 = localextrema(a, fse)
                @test samevalues(e0, e1) && samevalues(d0, d1)
            end
            @testset "localmean" begin
                result = localmean(REF, a, box)
                @test similarvalues(localmean(a, box), result)
                @test similarvalues(localmean(a, mask), result)
                @test similarvalues(localmean(a, kern), result)
            end
            @testset "convolve" begin
                result = convolve(REF, a, kern)
                @test similarvalues(convolve(a, box), result)
                @test similarvalues(convolve(a, mask), result)
                @test similarvalues(convolve(a, kern), result)
            end
        end
    end

    ball3x3 = Bool[1 1 1;
                   1 1 1;
                   1 1 1];
    ball5x5 = Bool[0 1 1 1 0;
                   1 1 1 1 1;
                   1 1 1 1 1;
                   1 1 1 1 1;
                   0 1 1 1 0];
    ball7x7 = Bool[0 0 1 1 1 0 0;
                   0 1 1 1 1 1 0;
                   1 1 1 1 1 1 1;
                   1 1 1 1 1 1 1;
                   1 1 1 1 1 1 1;
                   0 1 1 1 1 1 0;
                   0 0 1 1 1 0 0];
    @test samevalues(ball3x3, ball(2, 1))
    @test samevalues(ball5x5, ball(2, 2))
    @test samevalues(ball7x7, ball(2, 3))

    @testset "2D test with a 5x5 ball" begin
        a = rand(162,181)
        mask = Kernel(ball5x5)
        fse = strel(eltype(a), mask)
        @testset "erode" begin
            result = erode(REF, a, mask)
            @test samevalues(erode(a, mask), result)
            @test samevalues(erode(a, fse), result)
        end
        @testset "dilate" begin
            result = dilate(REF, a, mask)
            @test samevalues(dilate(a, mask), result)
            @test samevalues(dilate(a, fse), result)
        end
    end

    # Test van Herk / Gil & Werman algorithm.
    @testset "van Herk / Gil & Werman algorithm" begin
        # Test erode and dilate.
        for dims in ((45,), (30, 21), (7, 8, 9)),
            T in (Float32, Int32)
            N = length(dims)
            A = rand(T, dims)
            B = Array{eltype(A),ndims(A)}(undef, size(A))
            C = Array{eltype(A),ndims(A)}(undef, size(A))
            box = RectangularBox{N}(5)
            rng = -2:2
            result = erode(REF, A, box)
            @test samevalues(erode(A, :, rng), result)
            @test samevalues(erode!(B, A, :, rng), result)
            @test samevalues(erode!(copyto!(B, A), :, rng), result)
            @test samevalues(localfilter!(B, A, 1:N, min, rng), result)
            if N > 1
                # Reverse order of dimensions.
                @test samevalues(localfilter!(B, A, N:-1:1, min, rng), result)
            end
            result = dilate(REF, A, box)
            @test samevalues(dilate(A, :, rng), result)
            @test samevalues(dilate!(B, A, :, rng), result)
            @test samevalues(dilate!(copyto!(B, A), :, rng), result)
            @test samevalues(localfilter!(B, A, 1:N, max, rng), result)
            if N > 1
                # Reverse order of dimensions.
                @test samevalues(localfilter!(B, A, N:-1:1, max, rng), result)
            end
        end

        # Test shifts.
        A = rand(30)
        B = Array{eltype(A),ndims(A)}(undef, size(A))
        n = length(A)
        for k in (2, 0, -3)
            for i in 1:n
                j = clamp(i - k, 1, n)
                B[i] = A[j]
            end
            @test samevalues(B, localfilter(A,1,min,k:k))
            C = copy(A)
            @test samevalues(B, localfilter!(C,1,min,k:k))
        end
        A = rand(12,13)
        B = Array{eltype(A),ndims(A)}(undef, size(A))
        C = Array{eltype(A),ndims(A)}(undef, size(A))
        n1, n2 = size(A)
        for k1 in (2, 0, -3),
            k2 in (1, 0, -2)
            for i2 in 1:n2
                j2 = clamp(i2 - k2, 1, n2)
                for i1 in 1:n1
                    j1 = clamp(i1 - k1, 1, n1)
                    B[i1,i2] = A[j1,j2]
                end
            end
            @test samevalues(B, localfilter(A,:,min,(k1:k1,k2:k2)))
            @test samevalues(B, localfilter(A,[1,2],min,(k1:k1,k2:k2)))
            @test samevalues(B, localfilter(A,(2,1),min,(k2:k2,k1:k1)))
            @test samevalues(B, localfilter!(copyto!(C,A),:,min,(k1:k1,k2:k2)))
            @test samevalues(B, localfilter!(copyto!(C,A),(1,2),min,(k1:k1,k2:k2)))
            @test samevalues(B, localfilter!(copyto!(C,A),[2,1],min,(k2:k2,k1:k1)))
        end
    end

    # Bilateral filter.
    @testset "Bilateral filter" begin
        A = randn(Float64, 128, 200)
        box = RectangularBox{2}(5)
        B1 = bilateralfilter(A, 4, 3, box)
        B2 = bilateralfilter(Float32, A, 4, 3, 5)
        B3 = bilateralfilter!(similar(Array{Float32}, axes(A)), A, 4, 3, box)
        @test similarvalues(B1, B2; atol=0, gtol=8*eps(Float32))
        @test similarvalues(B2, B3; atol=0, gtol=4*eps(Float32))
    end
end

end # module