text:
"""
    state = init_update_state(conf, gen_fn::GenerativeFunction, param_list::Vector)

Get the initial state for a parameter update to the given parameters of the given generative function.

`param_list` is a vector of references to parameters of `gen_fn`.
`conf` configures the update.
"""
function init_update_state end

"""
    apply_update!(state)

Apply one parameter update, mutating the values of the trainable parameters, and possibly also the given state.
"""
function apply_update! end

"""
    update = ParamUpdate(conf, param_lists...)

Return an update configured by `conf` that applies to set of parameters defined by `param_lists`.

Each element in `param_lists` value is is pair of a generative function and a vector of its parameter references.

**Example**. To construct an update that applies a gradient descent update to the parameters `:a` and `:b` of generative function `foo` and the parameter `:theta` of generative function `:bar`:

```julia
update = ParamUpdate(GradientDescent(0.001, 100), foo => [:a, :b], bar => [:theta])
```

------------------------------------------------------------------------------------------
Syntactic sugar for the constructor form above.

    update = ParamUpdate(conf, gen_fn::GenerativeFunction)

Return an update configured by `conf` that applies to all trainable parameters owned by the given generative function.

Note that trainable parameters not owned by the given generative function will not be updated, even if they are used during execution of the function.

**Example**. If generative function `foo` has parameters `:a` and `:b`, to construct an update that applies a gradient descent update to the parameters `:a` and `:b`:

```julia
update = ParamUpdate(GradientDescent(0.001, 100), foo)
```
"""
struct ParamUpdate
    states::Dict{GenerativeFunction,Any}
    conf::Any
    function ParamUpdate(conf, param_lists...)
        states = Dict{GenerativeFunction,Any}()
            for (gen_fn, param_list) in param_lists
            states[gen_fn] = init_update_state(conf, gen_fn, param_list)
        end
        new(states, conf)
    end
    function ParamUpdate(conf, gen_fn::GenerativeFunction)
        param_lists = Dict(gen_fn => collect(get_params(gen_fn)))
        ParamUpdate(conf, param_lists...)
    end
end

"""
    apply!(update::ParamUpdate)

Perform one step of the update.
"""
function apply!(update::ParamUpdate)
    for (_, state) in update.states
        apply_update!(state)
    end
    nothing
end

"""
    conf = FixedStepGradientDescent(step_size)

Configuration for stochastic gradient descent update with fixed step size.
"""
struct FixedStepGradientDescent
    step_size::Float64
end

"""
    conf = GradientDescent(step_size_init, step_size_beta)

Configuration for stochastic gradient descent update with step size given by `(t::Int) -> step_size_init * (step_size_beta + 1) / (step_size_beta + t)` where `t` is the iteration number.
"""
struct GradientDescent
    step_size_init::Float64
    step_size_beta::Float64
end

"""
    conf = ADAM(learning_rate, beta1, beta2, epsilon)

Configuration for ADAM update.
"""
struct ADAM
    learning_rate::Float64
    beta1::Float64
    beta2::Float64
    epsilon::Float64
end

export ParameterSet, ParamUpdate, apply!
export FixedStepGradientDescent, GradientDescent, ADAM


text:
export NormalMeanVariance

import StatsFuns: log2π, invsqrt2π

struct NormalMeanVariance{T <: Real} <: ContinuousUnivariateDistribution
    μ :: T
    v :: T
end

NormalMeanVariance(μ::Real, v::Real)       = NormalMeanVariance(promote(μ, v)...)
NormalMeanVariance(μ::Integer, v::Integer) = NormalMeanVariance(float(μ), float(v))
NormalMeanVariance(μ::T) where {T <: Real} = NormalMeanVariance(μ, one(T))
NormalMeanVariance()                       = NormalMeanVariance(0.0, 1.0)

Distributions.@distr_support NormalMeanVariance -Inf Inf

Distributions.support(dist::NormalMeanVariance) = Distributions.RealInterval(minimum(dist), maximum(dist))

weightedmean(dist::NormalMeanVariance) = precision(dist) * mean(dist)

function weightedmean_invcov(dist::NormalMeanVariance)
    w = invcov(dist)
    xi = w * mean(dist)
    return (xi, w)
end

Distributions.mean(dist::NormalMeanVariance)      = dist.μ
Distributions.median(dist::NormalMeanVariance)    = mean(dist)
Distributions.mode(dist::NormalMeanVariance)      = mean(dist)
Distributions.var(dist::NormalMeanVariance)       = dist.v
Distributions.std(dist::NormalMeanVariance)       = sqrt(var(dist))
Distributions.cov(dist::NormalMeanVariance)       = var(dist)
Distributions.invcov(dist::NormalMeanVariance)    = inv(cov(dist))
Distributions.entropy(dist::NormalMeanVariance)   = (1 + log2π + log(var(dist))) / 2

Distributions.pdf(dist::NormalMeanVariance, x::Real)    = (invsqrt2π * exp(-abs2(x - mean(dist)) / 2cov(dist))) / std(dist)
Distributions.logpdf(dist::NormalMeanVariance, x::Real) = -(log2π + log(var(dist)) + abs2(x - mean(dist)) / var(dist)) / 2

Base.precision(dist::NormalMeanVariance{T}) where T = invcov(dist)
Base.eltype(::NormalMeanVariance{T})        where T = T

Base.convert(::Type{ NormalMeanVariance }, μ::Real, v::Real) = NormalMeanVariance(μ, v)
Base.convert(::Type{ NormalMeanVariance{T} }, μ::Real, v::Real) where { T <: Real } = NormalMeanVariance(convert(T, μ), convert(T, v))

vague(::Type{ <: NormalMeanVariance }) = NormalMeanVariance(0.0, huge)

prod_analytical_rule(::Type{ <: NormalMeanVariance }, ::Type{ <: NormalMeanVariance }) = ProdAnalyticalRuleAvailable()

function Base.prod(::ProdPreserveType, left::NormalMeanVariance, right::NormalMeanVariance)
    μ = (mean(left) * var(right) + mean(right) * var(left)) / (var(right) + var(left))
    v = (var(left) * var(right)) / (var(left) + var(right))
    return NormalMeanVariance(μ, v)
end

function Base.prod(::ProdAnalytical, left::NormalMeanVariance, right::NormalMeanVariance)
    xi = mean(left)/var(left) + mean(right)/var(right)
    w = 1/var(left) + 1/var(right)
    return NormalWeightedMeanPrecision(xi, w)
end

ext: jl


text:
using Documenter, SubstitutionModels

makedocs(
    modules = [SubstitutionModels],
    format = :html,
    sitename = "SubstitutionModels.jl",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Substitution models" => "man/models.md",
            "Q and P matrices" => "man/pqmatrices.md",
            "Provided models & custom models" => "man/modeltypes.md"
        ],
        "Contributing" => "contributing.md"
    ],
    authors = "Justin Angevarre & Ben J. Ward, at the BioJulia organisation, and contributors."
)

deploydocs(
    repo = "github.com/BioJulia/SubstitutionModels.jl.git",
    julia = "0.7",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
module TestSetterMacro
try
    using Test
catch
    using Base.Test
end
using Reconstructables: @add_kwonly, @recon, statements
include("utils.jl")

@testset "statements(ex::Expr)" begin
    let actual = statements(quote a; b; c; end)
        desired = [:a, :b, :c]
        @test actual == desired
    end

    let actual = statements(quote
                            a = 1
                            b = 2
                            c = 3
                            end)
        desired = [:(a = 1), :(b = 2), :(c = 3)]
        @test actual == desired
    end

    let actual = statements(quote
                            a
                            begin
                                b
                                begin
                                    c
                                end
                            end
                            end)
        desired = [:a, :b, :c]
        @test actual == desired
    end
end

struct A
    x
    @add_kwonly A(x) = new(x)
end

old = A(A(A(1)))
new = @recon old.x.x.x = 2
@test new.x.x.x == 2

struct B
    x
    y
    @add_kwonly B(x, y=nothing) = new(x, y)
end

old = B(B(B(1), B(2, 3)))
new = @recon begin
    old.x.x.x = 10
    old.x.y.y = 20
end
@test new.x.x.x == 10
@test new.x.y.y == 20

new = @recon begin
    old.x.x.x = 10
    old.x.x.y = 20
    old.x.y = 30
    old.y = 40
end
@test new.x.x.x == 10
@test new.x.x.y == 20
@test new.x.y == 30
@test new.y == 40

f(x) = 2x

new = @recon begin
    old.x.x.x = f(old.x.x.x)
    old.x.x.y = 1 * 20
    old.x.y.x = let x = old.x.y.x
        10x
    end
end
@test old.x.x.x == 1
@test new.x.x.x == 2
@test new.x.x.y == 20
@test old.x.y.x == 2
@test new.x.y.x == 20

__not_implemented__ = """
old = B(B(B(1), B(2, 3)))
new = @recon let f = old.x.x,
                 g = old.x.y
    f.x = 10
    g.y = 20
end
@test new.x.x.x == 10
@test new.x.y.y == 20
"""

invalid_expression_list = statements(quote
    old{:spam}.x.x.x = 10
    old.x{:spam}.x.x = 10
    old.x.x.x{:spam} = 10
    f(x) = 10
    old.f(x) = 10
end)
append!(invalid_expression_list,
        map(ex -> quote $ex end, invalid_expression_list))

for ex in invalid_expression_list
    @test_error begin
        @eval @recon $ex
    end (err) -> begin
        @test err isa ErrorException
        @test contains(err.msg, "Unsupported expression:")
    end
end

@test_error begin
    @eval @recon begin
        old.x = 1
        new.x = 2
    end
end (err) -> begin
    @test err isa ErrorException
    @test contains(err.msg, "Trying to modify multiple structs:")
end

@test_error begin
    @eval @recon begin
        1 + 1
    end
end (err) -> begin
    @test err isa ErrorException
    @test contains(err.msg, "All statements in a block must be assignments.")
end

end


text:
@testset "normalizedkernel" begin
    rng = MersenneTwister(123456)
    x = randn(rng)
    y = randn(rng)

    k = 4 * SqExponentialKernel()
    kn = NormalizedKernel(k)
    @test kn(x, y) == k(x, y) / sqrt(k(x, x) * k(y, y))
    @test kn(x, x) ≈ one(x) atol = 1e-5

    # Standardised tests.
    TestUtils.test_interface(kn, Float64)
    test_ADs(x -> NormalizedKernel(exp(x[1]) * SqExponentialKernel()), rand(1))

    test_params(kn, k)
end


text:
# Autogenerated wrapper script for OpenLDAPClient_jll for x86_64-apple-darwin14
export ldapcompare, ldapdelete, ldapexop, ldapmodify, ldapmodrdn, ldappasswd, ldapsearch, ldapurl, ldapwhoami, liblber, libldap, libldap_r

using OpenSSL_jll
using PCRE_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "DYLD_FALLBACK_LIBRARY_PATH"
LIBPATH_default = "~/lib:/usr/local/lib:/lib:/usr/lib"

# Relative path to `ldapcompare`
const ldapcompare_splitpath = ["bin", "ldapcompare"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ldapcompare_path = ""

# ldapcompare-specific global declaration
function ldapcompare(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ldapcompare_path)
    end
end

# Relative path to `ldapdelete`
const ldapdelete_splitpath = ["bin", "ldapdelete"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ldapdelete_path = ""

# ldapdelete-specific global declaration
function ldapdelete(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ldapdelete_path)
    end
end

# Relative path to `ldapexop`
const ldapexop_splitpath = ["bin", "ldapexop"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ldapexop_path = ""

# ldapexop-specific global declaration
function ldapexop(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ldapexop_path)
    end
end

# Relative path to `ldapmodify`
const ldapmodify_splitpath = ["bin", "ldapmodify"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ldapmodify_path = ""

# ldapmodify-specific global declaration
function ldapmodify(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ldapmodify_path)
    end
end

# Relative path to `ldapmodrdn`
const ldapmodrdn_splitpath = ["bin", "ldapmodrdn"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ldapmodrdn_path = ""

# ldapmodrdn-specific global declaration
function ldapmodrdn(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ldapmodrdn_path)
    end
end

# Relative path to `ldappasswd`
const ldappasswd_splitpath = ["bin", "ldappasswd"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ldappasswd_path = ""

# ldappasswd-specific global declaration
function ldappasswd(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ldappasswd_path)
    end
end

# Relative path to `ldapsearch`
const ldapsearch_splitpath = ["bin", "ldapsearch"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ldapsearch_path = ""

# ldapsearch-specific global declaration
function ldapsearch(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ldapsearch_path)
    end
end

# Relative path to `ldapurl`
const ldapurl_splitpath = ["bin", "ldapurl"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ldapurl_path = ""

# ldapurl-specific global declaration
function ldapurl(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ldapurl_path)
    end
end

# Relative path to `ldapwhoami`
const ldapwhoami_splitpath = ["bin", "ldapwhoami"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ldapwhoami_path = ""

# ldapwhoami-specific global declaration
function ldapwhoami(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ldapwhoami_path)
    end
end

# Relative path to `liblber`
const liblber_splitpath = ["lib", "liblber.dylib"]

# This will be filled out by __init__() for all products, as it must be done at runtime
liblber_path = ""

# liblber-specific global declaration
# This will be filled out by __init__()
liblber_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const liblber = "@rpath/liblber-2.4.2.dylib"

# Relative path to `libldap`
const libldap_splitpath = ["lib", "libldap.dylib"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libldap_path = ""

# libldap-specific global declaration
# This will be filled out by __init__()
libldap_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libldap = "@rpath/libldap-2.4.2.dylib"

# Relative path to `libldap_r`
const libldap_r_splitpath = ["lib", "libldap_r.dylib"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libldap_r_path = ""

# libldap_r-specific global declaration
# This will be filled out by __init__()
libldap_r_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libldap_r = "@rpath/libldap_r-2.4.2.dylib"

# Inform that the wrapper is available for this platform
wrapper_available = true

"""
Open all libraries
"""
function __init__()
    # This either calls `@artifact_str()`, or returns a constant string if we're overridden.
    global artifact_dir = find_artifact_dir()

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # Initialize PATH and LIBPATH environment variable listings
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (OpenSSL_jll.PATH_list, PCRE_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (OpenSSL_jll.LIBPATH_list, PCRE_jll.LIBPATH_list,))

    global ldapcompare_path = normpath(joinpath(artifact_dir, ldapcompare_splitpath...))

    push!(PATH_list, dirname(ldapcompare_path))
    global ldapdelete_path = normpath(joinpath(artifact_dir, ldapdelete_splitpath...))

    push!(PATH_list, dirname(ldapdelete_path))
    global ldapexop_path = normpath(joinpath(artifact_dir, ldapexop_splitpath...))

    push!(PATH_list, dirname(ldapexop_path))
    global ldapmodify_path = normpath(joinpath(artifact_dir, ldapmodify_splitpath...))

    push!(PATH_list, dirname(ldapmodify_path))
    global ldapmodrdn_path = normpath(joinpath(artifact_dir, ldapmodrdn_splitpath...))

    push!(PATH_list, dirname(ldapmodrdn_path))
    global ldappasswd_path = normpath(joinpath(artifact_dir, ldappasswd_splitpath...))

    push!(PATH_list, dirname(ldappasswd_path))
    global ldapsearch_path = normpath(joinpath(artifact_dir, ldapsearch_splitpath...))

    push!(PATH_list, dirname(ldapsearch_path))
    global ldapurl_path = normpath(joinpath(artifact_dir, ldapurl_splitpath...))

    push!(PATH_list, dirname(ldapurl_path))
    global ldapwhoami_path = normpath(joinpath(artifact_dir, ldapwhoami_splitpath...))

    push!(PATH_list, dirname(ldapwhoami_path))
    global liblber_path = normpath(joinpath(artifact_dir, liblber_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global liblber_handle = dlopen(liblber_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(liblber_path))

    global libldap_path = normpath(joinpath(artifact_dir, libldap_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libldap_handle = dlopen(libldap_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libldap_path))

    global libldap_r_path = normpath(joinpath(artifact_dir, libldap_r_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libldap_r_handle = dlopen(libldap_r_path, RTLD_LAZY | RTLD_DEEPBIND)
    push!(LIBPATH_list, dirname(libldap_r_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
using Clang

const RS_INCLUDE = joinpath(@__DIR__, "..", "deps", "usr", "include", "librealsense2") |> normpath
const RS_HEADERS = [joinpath(RS_INCLUDE, header) for header in readdir(RS_INCLUDE) if endswith(header, ".h")]

# create a work context
ctx = DefaultContext()

# parse headers
parse_headers!(ctx, RS_HEADERS,
               args=["-I", joinpath(RS_INCLUDE, "..")],
               includes=vcat(RS_INCLUDE, LLVM_INCLUDE),
               )

# settings
ctx.libname = "librealsense2"
ctx.options["is_function_strictly_typed"] = false
ctx.options["is_struct_mutable"] = false

# write output
api_file = joinpath(@__DIR__, "rs2_api.jl")
api_stream = open(api_file, "w")

for trans_unit in ctx.trans_units
    root_cursor = getcursor(trans_unit)
    push!(ctx.cursor_stack, root_cursor)
    header = spelling(root_cursor)
    @info "wrapping header: $header ..."
    # loop over all of the child cursors and wrap them, if appropriate.
    ctx.children = children(root_cursor)
    for (i, child) in enumerate(ctx.children)
        child_name = name(child)
        child_header = filename(child)
        ctx.children_index = i
        # choose which cursor to wrap
        startswith(child_name, "__") && continue  # skip compiler definitions
        child_name in keys(ctx.common_buffer) && continue  # already wrapped
        child_header != header && continue  # skip if cursor filename is not in the headers to be wrapped

        wrap!(ctx, child)
    end
    @info "writing $(api_file)"
    println(api_stream, "# Julia wrapper for header: $header")
    println(api_stream, "# Automatically generated using Clang.jl\n")
    print_buffer(api_stream, ctx.api_buffer)
    empty!(ctx.api_buffer)  # clean up api_buffer for the next header
end
close(api_stream)

# write "common" definitions: types, typealiases, etc.
common_file = joinpath(@__DIR__, "rs2_common.jl")
open(common_file, "w") do f
    println(f, "# Automatically generated using Clang.jl\n")
    print_buffer(f, dump_to_buffer(ctx.common_buffer))
end

# uncomment the following code to generate dependency and template files
# copydeps(dirname(api_file))
# print_template(joinpath(dirname(api_file), "LibTemplate.jl"))


text:
"""
    Codec

Codecs provide the encoding-decoding scheme in a flexible way.
Codecs are titled according to the general attributes that the codec scheme
would require, and NOT the actual encoding/decoding function itself.

This makes things customizable, and allows plug-and-play behaviour with
different variants of the same codec scheme, or for strikingly similar codec
schemes.
"""
module Codec

using ..ChemistryFeaturization.AbstractType: AbstractCodec
import ..ChemistryFeaturization.output_shape

include("simplecodec.jl")
export SimpleCodec

include("onehotonecold.jl")
export OneHotOneCold, output_shape

include("directcodec.jl")
export DirectCodec

export encode, decode

end


text:
"""
Alam, S., Sonbhadra, S.K., Agarwal, S., Nagabhushan, P. and Tanveer, M., 2020.
Sample reduction using Farthest Boundary Point Estimation (FBPE) for Support Vector Data Description (SVDD).
Pattern Recognition Letters.

Farthest Boundary Point Estimation (FBPE)
Parameter:
n: number of directions to check (no stepsize recommendation for the full sweep from 0 to 2π )

To compute the angle we mean center the data (so mean is 0 vector).
Then we select a random data point as reference point and compute the angle to all other points.
All these angles are then binned into n bins of equal size.
With cos-sim ∈ [-1, 1], the value range covered by each bin is 2 / n large.
"""
struct FBPE <: Sampler
    n::Int
    function FBPE(n::Int)
        n > 0 || throw(ArgumentError("Non-positive number of directions to check $(n)."))
        return new(n)
    end
end

function sample(sampler::FBPE, data::DataSet, labels::Vector{Symbol})::BitArray{1}
    rejected = falses(size(data, 2))
    F_x = Int[]
    while !all(rejected)
        remaining_indices = findall(.!rejected)
        if length(remaining_indices) == 1
            push!(F_x, first(remaining_indices))
            break
        end
        # Normalize data so mean is the center and we can perform a sweep in n directions
        x_mean = mean(data[:, remaining_indices], dims=2)
        data_centered = data[:, remaining_indices] .- x_mean
        x_ref = data_centered[:, 1]
        # compute cosine similarity
        similarity = [cosine_similarity(x_ref, data_centered[:, i]) for i in 1:size(data_centered, 2)]
        # Iterate over all directions
        S_x = []
        bin_thresholds = collect(range(-1 - 10^-5, 1, length=sampler.n + 1))
        for t in 2:length(bin_thresholds)
            bin_indices = findall((similarity .> bin_thresholds[t-1]) .& (similarity .<= bin_thresholds[t]))
            if !isempty(bin_indices)
                max_dist, max_idx = findmax([norm(data_centered[:, i]) for i in bin_indices])
                new_reject_indices = [i for i in bin_indices if i != bin_indices[max_idx]]
                # only keep farthest point, reject others
                rejected[remaining_indices[new_reject_indices]] .= true
                push!(S_x, (max_dist, remaining_indices[bin_indices[max_idx]]))
            end
        end
        # select farthest point of all farthest points per direction
        sort!(S_x, rev=true)
        farthest_idx = first(S_x)[2]
        rejected[farthest_idx] = true
        push!(F_x, farthest_idx)
    end
    sample_mask = falses(size(data, 2))
    sample_mask[F_x] .= true
    return sample_mask
end


text:
# Paths in the representation domain of the BZ, obtained from the [SeeK] publication.
# Methods return **k**-vector coordinates in a primitive reciprocal basis ("reciprocal
# crystallographic primitive cell" in the [SeeK]'s nomenclature).
#
# ^[Seek] http://dx.doi.org/10.1016/j.commatsci.2016.10.015; see also online interface at
#         https://www.materialscloud.org/work/tools/seekpath

module KPaths

# ---------------------------------------------------------------------------------------- #
export irrfbz_path, KPath, points, paths, KPathInterpolant
# ---------------------------------------------------------------------------------------- #
using ..CrystallineBravaisVendor:
    bravaistype,
    boundscheck_sgnum,
    reciprocalbasis,
    primitivize_reciprocal
using ..Brillouin:
    AVec,
    BasisLike,
    SHOWDIGITS,
    BasisEnum, CARTESIAN, LATTICE, setting, set_setting!,
    cartesianize, latticize
import ..Brillouin:
    latticize!, cartesianize!,
    basis
using LinearAlgebra: norm, dot, ×
using StaticArrays
using DocStringExtensions

import Base: show, summary, getindex, setindex!, IndexStyle, size
# ---------------------------------------------------------------------------------------- #

include("codegen-kpoints.jl")  # defines `get_points_3d`, `pathsd_3d` (& other subfunctions)
include("bravais-branches.jl") # defines `extended_bravais`

# ---------------------------------------------------------------------------------------- #

abstract type AbstractPath{T} <: AbstractVector{T} end

"""
$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct KPath{D} <: AbstractPath{Pair{Symbol, SVector{D, Float64}}}
    points  :: Dict{Symbol, SVector{D,Float64}}
    paths   :: Vector{Vector{Symbol}}
    basis   :: SVector{D, SVector{D, Float64}}
    setting :: Ref{BasisEnum}
end

"""
    points(kp::KPath{D}) -> Dict{Symbol, SVector{D,Float64}}

Return a dictionary of the **k**-points (values) and associated **k**-labels (keys)
referenced in `kp`.
"""
points(kp::KPath) = kp.points

"""
    paths(kp::KPath) -> Vector{Vector{Symbol}}

Return a vector of vectors, with each vector describing a connected path between between
**k**-points referenced in `kp` (see also [`points(::KPath)`](@ref)).
"""
paths(kp::KPath) = kp.paths

Base.@propagate_inbounds function getindex(kp::KPath, i::Int)
    # index into the `i`th point in the "flattened" `paths(kp)`
    @boundscheck i < 0 && throw(BoundsError(kp, i))
    j = 1
    i′ = i′′ = i
    while (i′′ -= length(paths(kp)[j])) > 0
        i′ = i′′
        j += 1
        @boundscheck j > length(paths(kp)) && throw(BoundsError(kp, i))
    end
    @boundscheck i′ > length(paths(kp)[j]) && throw(BoundsError(kp, i))

    lab = @inbounds paths(kp)[j][i′]
    return lab => points(kp)[lab]
end
IndexStyle(::Type{<:KPath}) = IndexLinear()
size(kp::KPath) = (sum(length, paths(kp)),)

function summary(io::IO, kp::KPath{D}) where D
    print(io, "KPath{", D, "} (", length(points(kp)), " points, ",
              length(paths(kp)), " paths, ", length(kp), " points in paths)")
end
function show(io::IO, ::MIME"text/plain", kp::KPath)
    summary(io, kp)
    println(io, ":")
    # points
    print(io, " points: ")
    foreach(enumerate(points(kp))) do (i,(lab,kv))
        i ≠ 1 && print(io, "         ")
        println(io, ":", lab, " => ", round.(kv, digits=SHOWDIGITS))
    end
    # paths
    print(io, "  paths: ")
    foreach(enumerate(paths(kp))) do (i,path)
        i ≠ 1 && print(io, "         ")
        println(io, path)
    end
    # basis
    print(io, "  basis: ")
    foreach(enumerate(basis(kp))) do (i,v)
        i ≠ 1 && print(io, "         ")
        if i ≠ length(basis(kp))
            println(io, round.(v, digits=SHOWDIGITS))
        else
            print(io, round.(v, digits=SHOWDIGITS))
        end
    end
end
# ---------------------------------------------------------------------------------------- #

"""
    irrfbz_path(sgnum::Integer, Rs, [::Union{Val(D), Integer},]=Val(3))  -->  ::KPath{D}

Returns a **k**-path (`::KPath`) in the (primitive) irreducible Brillouin zone for a space
group with number `sgnum`, (conventional) direct lattice vectors `Rs`, and dimension `D`.
The path includes all distinct high-symmetry lines and points as well as relevant parts of
the Brillouin zone boundary.

The dimension `D` (1, 2, or 3) is specified as the third input argument, preferably as a
static `Val{D}` type parameter (or, type-unstably, as an `<:Integer`). Defaults to `Val(3)`.

`Rs` refers to the direct basis of the conventional unit cell, i.e., not the primitive
direct basis vectors. The setting of `Rs` must agree with the conventional setting choices
in the International Tables of Crystallography.
If `Rs` is a subtype of a `StaticVector` or `NTuple`, the dimension can be inferred from its
(static) size; in this case, this dimension will take precedence (i.e. override, if
different) over any dimension specified in the third input argument.

## Notes
- The returned **k**-points are given in the basis of the **primitive** reciprocal basis
  (see [`cartesianize!`](@ref)).
- To interpolate a `KPath`, see [`interpolate(::KPath, ::Integer)`](@ref) and
  [`splice(::KPath, ::Integer)`](@ref).
- All paths currently assume time-reversal symmetry (or, equivalently, inversion symmetry).
  If neither are present, include the "inverted" -**k** paths manually.

## Data and referencing
3D paths are sourced from the SeeK HPKOT publication: please cite the original work [^1].

[^1] Hinuma, Pizzi, Kumagai, Oba, & Tanaka, *Band structure diagram paths based on
     crystallography*,
     [Comp. Mat. Sci. **128**, 140 (2017)](http://dx.doi.org/10.1016/j.commatsci.2016.10.015)
"""
function irrfbz_path(sgnum::Integer, Rs, Dᵛ::Val{D}=Val(3)) where D
    D′ = length(Rs)
    D′ ≠ D && throw(DimensionMismatch("inconsistent dimensions of `Rs` and `Dᵛ`"))
    any(R -> length(R) ≠ D, Rs) && throw(DimensionMismatch("inconsistent element dimensions in `Rs`"))
    Rs = convert(SVector{D, SVector{D, Float64}}, Rs)

    # (extended) bravais type
    bt = bravaistype(sgnum, Dᵛ)
    ext_bt = extended_bravais(sgnum, bt, Rs, Dᵛ)

    # get data about path
    lab2kvs = get_points(ext_bt, Rs, Dᵛ)
    paths = get_paths(ext_bt, Dᵛ)

    # compute (primitive) reciprocal basis
    cntr = last(bt)
    pGs = primitivize_reciprocal(reciprocalbasis(Rs), cntr)

    return KPath(lab2kvs, paths, pGs, Ref(LATTICE))
end
function irrfbz_path(sgnum::Integer,
            Rs::Union{StaticVector{D, <:AVec{<:Real}}, NTuple{D, <:AVec{<:Real}}}) where D
    return irrfbz_path(sgnum, Rs, Val(D))
end
irrfbz_path(sgnum::Integer, Rs::AVec{<:AVec{<:Real}}, D::Integer) = irrfbz_path(sgnum, Rs, Val(D))

function get_points(ext_bt::Symbol,
                    Rs::Union{Nothing, AVec{<:AVec{<:Real}}},
                    Dᵛ::Val{D}) where D
    if Dᵛ === Val(3)
        return get_points_3d(ext_bt, Rs)
    elseif Dᵛ === Val(2)
        return get_points_2d(ext_bt, Rs)
    elseif Dᵛ === Val(1)
        ext_bt === :lp || throw(DomainError(ext_bt, "invalid extended Bravais type"))
        return Dict(:Γ => SVector(0.0), :X => SVector(0.5))
    else
        throw(DomainError(D, "unsupported input dimension"))
    end
end

function get_paths(ext_bt::Symbol, Dᵛ::Val{D}) where D
    if Dᵛ === Val(3)
        paths = get(pathsd_3d, ext_bt, nothing)
    elseif Dᵛ === Val(2)
        paths = get(pathsd_2d, ext_bt, nothing)
    elseif Dᵛ === Val(1)
        paths = ext_bt === :lp ? [[:Γ, :X]] : nothing
    else
        throw(DomainError(D, "unsupported input dimension"))
    end

    paths === nothing && throw(DomainError(ext_bt, "invalid extended Bravais type"))
    return paths
end

# ---------------------------------------------------------------------------------------- #

include("interpolate-paths.jl")
export interpolate, splice, cumdists

# ---------------------------------------------------------------------------------------- #

"""
    cartesianize!(kp::KPath{D}, Gs::BasisLike{D})

Transform a **k**-path `kp` in a lattice basis to a Cartesian basis with (primitive)
reciprocal basis vectors `basis`.
Modifies `kp` in-place.
"""
function cartesianize!(kp::KPath{D}) where D
    setting(kp) === CARTESIAN && return kp
    for (lab, kv) in points(kp)
        @inbounds points(kp)[lab] = cartesianize(kv, kp.basis)
    end
    set_setting!(kp, CARTESIAN)
    return kp
end

"""
    latticize!(kp::KPath{D}, Gs::BasisLike{D})

Transform a **k**-path `kpi` in a Cartesian basis to a lattice basis with (primitive)
reciprocal lattice vectors `basis`.
Modifies `kp` in-place.
"""
function latticize!(kp::KPath{D}) where D
    setting(kp) === LATTICE && return kp
    basismatrix = hcat(kp.basis...)
    for (lab, kv) in points(kp)
        @inbounds points(kp)[lab] = latticize(kv, basismatrix)
    end
    set_setting!(kp, LATTICE)
    return kp
end

# ---------------------------------------------------------------------------------------- #
end # module


text:
using Pkg
Pkg.add("BinDeps")
Pkg.add("CondaBinDeps")
Pkg.add("Libdl")
using BinDeps
import CondaBinDeps

if VERSION >= v"0.7.0"
    using Libdl
end

function validate_netcdf_version(name,handle)
    f = Libdl.dlsym_e(handle, "nc_inq_libvers")
    #
    # Example
    # banner = "4.5.6.7 of Apr 1 2000 00:00:00"
    banner = unsafe_string(ccall(f,Ptr{UInt8},()))

    # remove the date
    # verstr = "4.5.6.7"
    verstr = split(banner)[1]

    # vernumbers = ["4","5","6","7"]
    vernumbers = split(verstr,r"[.-]")

    # major_minor_patch = ["4","5","6"]
    major_minor_patch = vernumbers[1:min(3,length(vernumbers))]

    # ver = v"4.5.6"
    ver = VersionNumber([parse(Int,s) for s in major_minor_patch]...)

    return ver > v"4.2"
end

@BinDeps.setup
libnetcdf = library_dependency("libnetcdf", aliases = ["libnetcdf4","libnetcdf-7","netcdf"], validate = validate_netcdf_version)

CondaBinDeps.Conda.add_channel("anaconda")
provides(CondaBinDeps.Manager, "libnetcdf", libnetcdf)
provides(AptGet, "libnetcdf-dev", libnetcdf, os = :Linux)
# provides(Yum, "netcdf-devel", libnetcdf, os = :Linux)

@BinDeps.install Dict(:libnetcdf => :libnetcdf)

using Documenter, ClimatePlots

makedocs(sitename = "ClimatePlots.jl",
format = :html,
pages = [
"index.md",
]

)

deploydocs(
repo   = "github.com/JuliaClimate/ClimatePlots.jl.git",
# julia  = "1.0",
# osname = "linux",
target = "build"
# deps = nothing,
# make = nothing
)


text:
using Documenter
using Pkg
Pkg.activate(".")
Pkg.instantiate()
using Psychometrics

makedocs(
    sitename = "Psychometrics",
    format = Documenter.HTML(),
    modules = [Psychometrics],
    doctest = true,
    pages = [
        "index.md",
        hide(
            "lib/public.md",
            map(
                s -> "methods/exported/$(s)",
                sort(readdir(joinpath(@__DIR__, "src/methods/exported"))),
            ),
        ),
        hide(
            "lib/internals.md",
            map(
                s -> "methods/internals/$(s)",
                sort(readdir(joinpath(@__DIR__, "src/methods/internals"))),
            ),
        ),
    ],
)

deploydocs(repo = ENV["REPO"], devurl = "docs", devbranch = ENV["DEVBRANCH"])


text:
import BioSimulator: get_ptype, get_neighbor_class, initialize_datastructs!, coordinates, label

@testset "Lattice" begin
  @testset "initialization" begin
    ##### reaction definition #####
    intermediate = @def_reactions begin
      A + B --> A + A, k1 # predation / predator reproduction
      B + 0 --> 0 + B, k2 # prey movement
      A + 0 --> 0 + A, k3 # predator movement
      B + 0 --> B + B, k4 # prey reproduction
      A --> 0, k5         # predator death
      B --> 0, k6         # prey death
    end k1 k2 k2 k3 k4 k5 k6

    ##### assign parameter values #####
    params = [2.0, 1.0, 1.0, 2.0, 1.0, 1.0]

    ##### generate the full reaction list #####
    model = @enumerate_with_sclass intermediate VonNeumann() 2 params

    ##### initial configuration #####
    initial_sites = [(0,1), (1,0), (1,1), (1,2), (1,3), (2,0), (2,1), (2,2), (3,0), (3,1)]
    initial_sites = [initial_sites[j][i] for i in 1:2, j in 1:length(initial_sites)]
    ptypes        = [1, 1, 1, 2, 2, 1, 2, 2, 2, 2]

    lattice = Lattice(initial_sites, ptypes)

    ##### correct answers #####
    expected_opensites = [
    (0, 0), (1, -1), (2, -1), (3, -1), (4, 0),
    (4, 1), (3, 2), (2, 3), (1, 4), (0, 3),
    (0, 2), (-1, 1)
    ]

    expected_compositions = [
    [4, 0, 0], [3, 1, 0], [3, 0, 1], [2, 2, 0], [2, 1, 1],
    [2, 0, 2], [1, 3, 0], [1, 2, 1], [1, 1, 2], [1, 0, 3],
    [0, 4, 0], [0, 3, 1], [0, 2, 2], [0, 1, 3], [0, 0, 4]
    ]

    expected_matrix = zeros(Int, 15, 3)

    # open site counts
    expected_matrix[2, 1] = 3
    expected_matrix[3, 1] = 5
    expected_matrix[4, 1] = 1
    expected_matrix[5, 1] = 1
    expected_matrix[6, 1] = 2
    # total = 12

    # predator counts
    expected_matrix[2,  2] = 1
    expected_matrix[4,  2] = 1
    expected_matrix[9,  2] = 1
    expected_matrix[13, 2] = 1
    # total = 4

    # prey counts
    expected_matrix[3,  3] = 1
    expected_matrix[5,  3] = 1
    expected_matrix[6,  3] = 2
    expected_matrix[9,  3] = 1
    expected_matrix[13, 3] = 1
    # total = 6

    expected_neighborclass = Dict(
    # particles
    (0, 1) => 2,
    (1, 0) => 4,
    (1, 1) => 13,
    (1, 2) => 9,
    (1, 3) => 3,
    (2, 0) => 9,
    (2, 1) => 13,
    (2, 2) => 6,
    (3, 0) => 5,
    (3, 1) => 6,
    # open sites
    (0, 0) => 4,
    (1, -1) => 2,
    (2, -1) => 2,
    (3, -1) => 3,
    (4, 0) => 3,
    (4, 1) => 3,
    (3, 2) => 6,
    (2, 3) => 6,
    (1, 4) => 3,
    (0, 3) => 3,
    (0, 2) => 5,
    (-1, 1) => 2
    )
    expected_sampleclass = Dict(
    # particles only
    (0, 1) => [11, 13],
    (1, 0) => [10, 13],
    (1, 1) => [2, 13],
    (1, 2) => [5, 14],
    (1, 3) => [7, 14],
    (2, 0) => [9, 2, 13],
    (2, 1) => [14],
    (2, 2) => [6, 14],
    (3, 0) => [6, 14],
    (3, 1) => [6, 14]
    )

    # perform the initialization step
    initialize_datastructs!(lattice, model)

    @testset "open site tracking" begin
      test_opensites = filter(x -> isequal(get_ptype(x), 1), lattice.site)

      # do we generate the correct number of open sites?
      @test length(test_opensites) == length(expected_opensites)

      # do we spawn the correct open sites?
      test_coordinates = [tuple(coordinates(x)...) for x in lattice.site]

      for coord in expected_opensites
        @test coord in test_coordinates
      end
    end

    @testset "compositions" begin
      test_compositions = model.enumeration.composition

      # do we generate the coorect number of compositions?
      @test length(test_compositions) == length(expected_compositions)

      # do we generate the correct compositions?
      for composition in expected_compositions
        @test composition in test_compositions
      end

      # do we iterate through the compositions in a specific order (open sites -> particle types...)?
      for k in eachindex(test_compositions)
        composition = test_compositions[k]
        @test composition == expected_compositions[k]
      end
    end

    @testset "neighborhood classes" begin
      # is each site assigned to the correct neighborhood class?
      for x in lattice.site
        coord = tuple(coordinates(x)...)
        @test get_neighbor_class(x) == expected_neighborclass[coord]
      end
    end

    @testset "sample classes" begin
      test_sampleclass = model.enumeration.class

      # is each site assigned to the correct sample classes?
      for x in lattice.site
        if get_ptype(x) > 1
          coord = tuple(coordinates(x)...)
          for s in expected_sampleclass[coord]
            @test label(x) in test_sampleclass[s]
          end
        end
      end
    end

    @testset "count matrix" begin
      test_matrix = zero(expected_matrix)

      for x in lattice.site
        l = get_ptype(x)
        k = get_neighbor_class(x)

        test_matrix[k,l] += 1
      end

      # is the count matrix correct?
      for i in eachindex(test_matrix)
        @test test_matrix[i] == expected_matrix[i]
      end
    end
  end
end


text:
import MathOptInterface
import SparseArrays

const MOI = MathOptInterface

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::Ptr{Cvoid}
    silent::Bool
    params::Dict{String,String}
    variable_start::Dict{MOI.VariableIndex,Float64}
    objective_constant::Float64
    solve_time::Float64
    termination_status::Cint

    """
        Optimizer()

    Create a new Cbc Optimizer.
    """
    function Optimizer(; kwargs...)
        model = new(
            Cbc_newModel(),
            false,
            Dict{String,String}(),
            Dict{MOI.VariableIndex,Float64}(),
            0.0,
            0.0,
            Cint(-1),
        )
        if length(kwargs) > 0
            @warn("""Passing optimizer attributes as keyword arguments to
            Cbc.Optimizer is deprecated. Use
                MOI.set(model, MOI.RawOptimizerAttribute("key"), value)
            or
                JuMP.set_optimizer_attribute(model, "key", value)
            instead.
            """)
        end
        for (key, value) in kwargs
            MOI.set(model, MOI.RawOptimizerAttribute("$(key)"), value)
        end
        finalizer(model) do m
            return Cbc_deleteModel(m)
        end
        return model
    end
end

Base.cconvert(::Type{Ptr{Cvoid}}, model::Optimizer) = model
Base.unsafe_convert(::Type{Ptr{Cvoid}}, model::Optimizer) = model.inner

function MOI.supports(::Optimizer, ::MOI.RawOptimizerAttribute)
    # TODO(odow): There is no programatical way throught the C API to check if a
    # parameter name (or value) is valid. Fix this upstream.
    return true
end

function MOI.set(model::Optimizer, param::MOI.RawOptimizerAttribute, value)
    return MOI.set(model, param, string(value))
end

function MOI.set(
    model::Optimizer,
    param::MOI.RawOptimizerAttribute,
    value::String,
)
    if !MOI.supports(model, param)
        throw(MOI.UnsupportedAttribute(param))
    end
    name = string(param.name)
    model.params[name] = value
    if !(model.silent && name == "logLevel")
        Cbc_setParameter(model, name, value)
    end
    return
end

function MOI.get(model::Optimizer, param::MOI.RawOptimizerAttribute)
    # TODO: This gives a poor error message if the name of the parameter is
    # invalid.
    return model.params[string(param.name)]
end

MOI.supports(::Optimizer, ::MOI.Silent) = true

function MOI.set(model::Optimizer, ::MOI.Silent, value::Bool)
    if value == model.silent
        return
    end
    log_level = value ? "0" : get(model.params, "logLevel", "1")
    Cbc_setParameter(model, "logLevel", log_level)
    model.silent = value
    return
end
MOI.get(model::Optimizer, ::MOI.Silent) = model.silent

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value)
    if value === nothing
        delete!(model.params, "seconds")
        Cbc_setParameter(model, "seconds", "??")
    else
        MOI.set(model, MOI.RawOptimizerAttribute("seconds"), value)
    end
    return
end

function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    value = get(model.params, "seconds", nothing)
    return value === nothing ? value : parse(Float64, value)
end

MOI.get(::Optimizer, ::MOI.SolverName) = "COIN Branch-and-Cut (Cbc)"

function MOI.empty!(model::Optimizer)
    Cbc_deleteModel(model)
    model.inner = Cbc_newModel()
    model.objective_constant = 0.0
    model.termination_status = Cint(-1)
    model.solve_time = 0.0
    for (name, value) in model.params
        Cbc_setParameter(model, name, value)
    end
    if model.silent
        Cbc_setParameter(model, "logLevel", "0")
    end
    empty!(model.variable_start)
    return
end

function MOI.is_empty(model::Optimizer)
    return Cbc_getNumCols(model) == 0 && Cbc_getNumRows(model) == 0
end

mutable struct _CbcModelFormat
    num_rows::Cint
    num_cols::Cint
    # (row_idx, col_idx, values) are the sparse elements of the constraint matrix.
    row_idx::Vector{Cint}
    col_idx::Vector{Cint}
    values::Vector{Float64}
    # Constraint bounds.
    row_lb::Vector{Float64}
    row_ub::Vector{Float64}
    # Variable bounds.
    col_lb::Vector{Float64}
    col_ub::Vector{Float64}
    # Columns that are binary or integer
    binary::Vector{Int}
    integer::Vector{Int}
    # SOS1 constraints
    sos1_starts::Vector{Cint}
    sos1_indices::Vector{Cint}
    sos1_weights::Vector{Float64}
    # SOS2 constraints
    sos2_starts::Vector{Cint}
    sos2_indices::Vector{Cint}
    sos2_weights::Vector{Float64}
    # Objective coefficients.
    obj::Vector{Float64}
    objective_constant::Float64

    # An `InexactError` might occur if `num_rows` or `num_cols` is too
    # large, e.g., if `num_cols isa Int64` and is larger than 2^31 on
    # 32-bit hardware.
    function _CbcModelFormat(num_rows::Cint, num_cols::Cint)
        return new(
            num_rows,
            num_cols,
            # Constraint matrix.
            Cint[],
            Cint[],
            Float64[],
            # Row lower/upper bounds.
            fill(-Inf, num_rows),
            fill(Inf, num_rows),
            # Column lower/upper bounds.
            fill(-Inf, num_cols),
            fill(Inf, num_cols),
            # Binary/Integer columns.
            Cint[],
            Cint[],
            # SOSI constraints.
            Cint[],
            Cint[],
            Float64[],
            # SOSII constraints.
            Cint[],
            Cint[],
            Float64[],
            # Objective vector and offset.
            fill(0.0, num_cols),
            0.0,
        )
    end
end

###
### VariableIndex-in-{EqualTo,LessThan,GreaterThan,Interval,ZeroOne,Integer}
###

_column_value(map::MOI.IndexMap, index::MOI.VariableIndex) = map[index].value

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VariableIndex},
    ::Type{
        <:Union{
            MOI.EqualTo{Float64},
            MOI.LessThan{Float64},
            MOI.GreaterThan{Float64},
            MOI.Interval{Float64},
            MOI.ZeroOne,
            MOI.Integer,
        },
    },
)
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{<:Union{MOI.SOS1{Float64},MOI.SOS2{Float64}}},
)
    return true
end

function _load_constraint(
    ::MOI.ConstraintIndex,
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    func::MOI.VariableIndex,
    set::MOI.EqualTo,
)
    column = _column_value(mapping, func)
    model.col_lb[column] = set.value
    model.col_ub[column] = set.value
    return
end

function _load_constraint(
    ::MOI.ConstraintIndex,
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    func::MOI.VariableIndex,
    set::MOI.LessThan,
)
    column = _column_value(mapping, func)
    model.col_ub[column] = set.upper
    return
end

function _load_constraint(
    ::MOI.ConstraintIndex,
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    func::MOI.VariableIndex,
    set::MOI.GreaterThan,
)
    column = _column_value(mapping, func)
    model.col_lb[column] = set.lower
    return
end

function _load_constraint(
    ::MOI.ConstraintIndex,
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    func::MOI.VariableIndex,
    set::MOI.Interval,
)
    column = _column_value(mapping, func)
    model.col_lb[column] = set.lower
    model.col_ub[column] = set.upper
    return
end

function _load_constraint(
    ::MOI.ConstraintIndex,
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    func::MOI.VariableIndex,
    ::MOI.ZeroOne,
)
    push!(model.binary, _column_value(mapping, func) - 1)
    return
end

function _load_constraint(
    ::MOI.ConstraintIndex,
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    func::MOI.VariableIndex,
    ::MOI.Integer,
)
    push!(model.integer, _column_value(mapping, func) - 1)
    return
end

function _load_constraints(
    model::_CbcModelFormat,
    src::MOI.ModelLike,
    mapping::MOI.IndexMap,
    F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet},
)
    for index in MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
        _load_constraint(
            index,
            model,
            mapping,
            MOI.get(src, MOI.ConstraintFunction(), index),
            MOI.get(src, MOI.ConstraintSet(), index),
        )
    end
    return
end

###
### ScalarAffineFunction-in-{EqualTo, LessThan, GreaterThan, Interval}
###

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{
        <:Union{
            MOI.EqualTo{Float64},
            MOI.LessThan{Float64},
            MOI.GreaterThan{Float64},
            MOI.Interval{Float64},
        },
    },
)
    return true
end

function _add_terms(
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    index::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S},
    func::MOI.ScalarAffineFunction{Float64},
) where {S}
    for term in func.terms
        push!(model.row_idx, mapping[index].value)
        push!(model.col_idx, _column_value(mapping, term.variable))
        push!(model.values, term.coefficient)
    end
    return
end

function _load_constraint(
    index::MOI.ConstraintIndex,
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    func::MOI.ScalarAffineFunction,
    set::MOI.EqualTo,
)
    _add_terms(model, mapping, index, func)
    row = mapping[index].value
    model.row_lb[row] = set.value - func.constant
    model.row_ub[row] = set.value - func.constant
    return
end

function _load_constraint(
    index::MOI.ConstraintIndex,
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    func::MOI.ScalarAffineFunction,
    set::MOI.GreaterThan,
)
    _add_terms(model, mapping, index, func)
    row = mapping[index].value
    model.row_lb[row] = set.lower - func.constant
    return
end

function _load_constraint(
    index::MOI.ConstraintIndex,
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    func::MOI.ScalarAffineFunction,
    set::MOI.LessThan,
)
    _add_terms(model, mapping, index, func)
    row = mapping[index].value
    model.row_ub[row] = set.upper - func.constant
    return
end

function _load_constraint(
    index::MOI.ConstraintIndex,
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    func::MOI.ScalarAffineFunction,
    set::MOI.Interval,
)
    _add_terms(model, mapping, index, func)
    row = mapping[index].value
    model.row_ub[row] = set.upper - func.constant
    model.row_lb[row] = set.lower - func.constant
    return
end

function _load_constraint(
    ::MOI.ConstraintIndex,
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    func::MOI.VectorOfVariables,
    set::MOI.SOS1{Float64},
)
    push!(model.sos1_starts, Cint(length(model.sos1_weights)))
    append!(model.sos1_weights, set.weights)
    for v in func.variables
        push!(model.sos1_indices, _column_value(mapping, v))
    end
    return
end

function _load_constraint(
    ::MOI.ConstraintIndex,
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    func::MOI.VectorOfVariables,
    set::MOI.SOS2{Float64},
)
    push!(model.sos2_starts, Cint(length(model.sos2_weights)))
    append!(model.sos2_weights, set.weights)
    for v in func.variables
        push!(model.sos2_indices, _column_value(mapping, v))
    end
    return
end

###
### ObjectiveSense
###

MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

function MOI.set(
    model::Optimizer,
    ::MOI.ObjectiveSense,
    sense::MOI.OptimizationSense,
)
    if sense == MOI.MAX_SENSE
        Cbc_setObjSense(model, -1.0)
    elseif sense == MOI.MIN_SENSE
        Cbc_setObjSense(model, 1.0)
    else
        @assert sense == MOI.FEASIBILITY_SENSE
        for col in Cint(0):Cint(Cbc_getNumCols(model) - 1)
            Cbc_setObjCoeff(model, col, 0.0)
        end
        Cbc_setObjSense(model, 0.0)
    end
    return
end

###
### ObjectiveFunction{ScalarAffineFunction}
###

function MOI.supports(
    ::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
)
    return true
end

function _load_objective(
    model::_CbcModelFormat,
    mapping::MOI.IndexMap,
    dest::Optimizer,
    src::MOI.ModelLike,
)
    F = MOI.get(src, MOI.ObjectiveFunctionType())
    if !MOI.supports(dest, MOI.ObjectiveFunction{F}())
        error("Objective function type $(F) not supported.")
    end
    f = MOI.get(src, MOI.ObjectiveFunction{F}())
    # We need to increment values of objective function with += to handle
    # cases like $x_1 + x_2 + x_1$. This is safe because objective function
    # is initialized with zeros in the constructor and `_load_objective` only
    # gets called once.
    for term in f.terms
        column = _column_value(mapping, term.variable)
        model.obj[column] += term.coefficient
    end
    model.objective_constant = f.constant
    return
end

###
### Variable starting values
###

function MOI.supports(
    ::Optimizer,
    ::MOI.VariablePrimalStart,
    ::Type{MOI.VariableIndex},
)
    return true
end

function MOI.set(
    model::Optimizer,
    ::MOI.VariablePrimalStart,
    x::MOI.VariableIndex,
    ::Nothing,
)
    delete!(model.variable_start, x)
    return
end

function MOI.set(
    model::Optimizer,
    ::MOI.VariablePrimalStart,
    x::MOI.VariableIndex,
    value,
)
    model.variable_start[x] = value
    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.VariablePrimalStart,
    x::MOI.VariableIndex,
)
    return get(model.variable_start, x, nothing)
end

###
### This main copy_to function.
###

function _create_constraint_indices_for_types(
    src::MOI.ModelLike,
    mapping::MOI.IndexMap,
    ::Type{MOI.VariableIndex},
    S::Type{<:MOI.AbstractSet},
    num_rows::Int,
)
    indices = MOI.get(src, MOI.ListOfConstraintIndices{MOI.VariableIndex,S}())
    for index in indices
        f = MOI.get(src, MOI.ConstraintFunction(), index)
        i = mapping[f].value
        mapping[index] = MOI.ConstraintIndex{MOI.VariableIndex,S}(i)
    end
    return num_rows
end

function _create_constraint_indices_for_types(
    src::MOI.ModelLike,
    mapping::MOI.IndexMap,
    F::Type{MOI.VectorOfVariables},
    S::Type{<:Union{MOI.SOS1{Float64},MOI.SOS2{Float64}}},
    num_rows::Int,
)
    n = 0
    for index in MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
        n += 1
        mapping[index] = MOI.ConstraintIndex{F,S}(n)
    end
    return num_rows
end

function _create_constraint_indices_for_types(
    src::MOI.ModelLike,
    mapping::MOI.IndexMap,
    F::Type{MOI.ScalarAffineFunction{Float64}},
    S::Type{<:MOI.AbstractScalarSet},
    num_rows::Int,
)
    for index in MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
        num_rows += 1
        mapping[index] = MOI.ConstraintIndex{F,S}(num_rows)
    end
    return num_rows
end

function _create_constraint_indices(
    dest::Optimizer,
    src::MOI.ModelLike,
    mapping::MOI.IndexMap,
)
    n = 0
    for (F, S) in MOI.get(src, MOI.ListOfConstraintTypesPresent())
        if !(MOI.supports_constraint(dest, F, S))
            throw(
                MOI.UnsupportedConstraint{F,S}(
                    "Cbc does not support constraints of type $F-in-$S.",
                ),
            )
        end
        # The type of `F` and `S` is not type-stable, so we use a function
        # barrier (`_create_constraint_indices_for_types`) to improve
        # performance.
        n = _create_constraint_indices_for_types(src, mapping, F, S, n)
    end
    return Cint(n)
end

function _create_variable_indices(src::MOI.ModelLike, mapping::MOI.IndexMap)
    for (i, x) in enumerate(MOI.get(src, MOI.ListOfVariableIndices()))
        mapping[x] = MOI.VariableIndex(i)
    end
    return Cint(length(mapping.var_map))
end

function MOI.copy_to(cbc_dest::Optimizer, src::MOI.ModelLike)
    @assert MOI.is_empty(cbc_dest)
    mapping = MOI.IndexMap()
    num_cols = _create_variable_indices(src, mapping)
    num_rows = _create_constraint_indices(cbc_dest, src, mapping)
    tmp_model = _CbcModelFormat(num_rows, num_cols)
    _load_objective(tmp_model, mapping, cbc_dest, src)
    for (F, S) in MOI.get(src, MOI.ListOfConstraintTypesPresent())
        # The type of `F` and `S` is not type-stable, so we use a function
        # barrier (`_load_constraints`) to improve performance.
        _load_constraints(tmp_model, src, mapping, F, S)
    end
    # Since Cbc doesn't have an explicit binary variable, we need to add [0, 1]
    # bounds and make it integer (which is done at the end of this function).
    for column in tmp_model.binary
        tmp_model.col_lb[column+1] = max(tmp_model.col_lb[column+1], 0.0)
        tmp_model.col_ub[column+1] = min(tmp_model.col_ub[column+1], 1.0)
    end
    A = SparseArrays.sparse(
        tmp_model.row_idx,
        tmp_model.col_idx,
        tmp_model.values,
        tmp_model.num_rows,
        tmp_model.num_cols,
    )
    Cbc_loadProblem(
        cbc_dest,
        tmp_model.num_cols,
        tmp_model.num_rows,
        A.colptr .- Cint(1),
        A.rowval .- Cint(1),
        A.nzval,
        tmp_model.col_lb,
        tmp_model.col_ub,
        tmp_model.obj,
        tmp_model.row_lb,
        tmp_model.row_ub,
    )
    MOI.Utilities.pass_attributes(
        cbc_dest,
        src,
        mapping,
        MOI.get(src, MOI.ListOfVariableIndices()),
    )
    for attr in MOI.get(src, MOI.ListOfModelAttributesSet())
        if attr isa MOI.ObjectiveFunction
            continue # Already copied
        end
        value = MOI.get(src, attr)
        if value !== nothing
            mapped_value = MOI.Utilities.map_indices(mapping, value)
            MOI.set(cbc_dest, attr, mapped_value)
        end
    end
    for (F, S) in MOI.get(src, MOI.ListOfConstraintTypesPresent())
        cis_src = MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
        MOI.Utilities.pass_attributes(cbc_dest, src, mapping, cis_src)
    end
    cbc_dest.objective_constant = tmp_model.objective_constant
    if length(tmp_model.integer) > 0
        Cbc_setInteger.(cbc_dest, tmp_model.integer)
    end
    if length(tmp_model.binary) > 0
        Cbc_setInteger.(cbc_dest, tmp_model.binary)
    end
    if length(tmp_model.sos1_starts) > 0
        push!(tmp_model.sos1_starts, Cint(length(tmp_model.sos1_weights)))
        Cbc_addSOS(
            cbc_dest,
            length(tmp_model.sos1_starts) - 1,
            tmp_model.sos1_starts,
            tmp_model.sos1_indices,
            tmp_model.sos1_weights,
            Cint(1),
        )
    end
    if length(tmp_model.sos2_starts) > 0
        push!(tmp_model.sos2_starts, Cint(length(tmp_model.sos2_weights)))
        Cbc_addSOS(
            cbc_dest,
            length(tmp_model.sos2_starts) - 1,
            tmp_model.sos2_starts,
            tmp_model.sos2_indices,
            tmp_model.sos2_weights,
            Cint(2),
        )
    end
    nsos = length(tmp_model.sos1_starts) + length(tmp_model.sos2_starts)
    if nsos > 0 && Cbc_getNumIntegers(cbc_dest) == 0
        @warn(
            "There are known correctness issues using Cbc with SOS " *
            "constraints and no binary variables.",
        )
    end
    return mapping
end

###
### Optimize and post-optimize functions
###

function _unsafe_wrap_cbc_array(
    model::Optimizer,
    f::F,
    n::Integer,
    indices;
    own::Bool = false,
) where {F<:Function}
    p = f(model)
    if p == C_NULL
        return map(x -> NaN, indices)
    end
    x = unsafe_wrap(Array, p, (n,); own = own)
    return x[indices]
end

function MOI.optimize!(model::Optimizer)
    if !isempty(model.variable_start)
        columns = fill(Cint(0), length(model.variable_start))
        values = fill(0.0, length(model.variable_start))
        for (i, (variable, value)) in enumerate(model.variable_start)
            columns[i] = Cint(variable.value - 1)
            values[i] = value
        end
        Cbc_setMIPStartI(model, length(columns), columns, values)
    end
    t = time()
    model.termination_status = Cbc_solve(model)
    model.solve_time = time() - t
    return
end

MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = model.solve_time

function MOI.get(model::Optimizer, ::MOI.NumberOfVariables)
    return Cbc_getNumCols(model)
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveBound)
    if Cbc_getNumIntegers(model) == 0
        return MOI.get(model, MOI.ObjectiveValue())
    end
    return Cbc_getBestPossibleObjValue(model) + model.objective_constant
end

function MOI.get(model::Optimizer, ::MOI.NodeCount)
    return Cbc_getNodeCount(model)
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    return Cbc_getObjValue(model) + model.objective_constant
end

# Cbc does not provide a native way of accessing the relative gap,
# use the Gurobi convention instead.
function MOI.get(model::Optimizer, ::MOI.RelativeGap)
    incumbent = MOI.get(model, MOI.ObjectiveValue())
    bound = MOI.get(model, MOI.ObjectiveBound())
    gap = abs(bound - incumbent) / abs(incumbent)
    return isnan(gap) ? Inf : gap
end

function MOI.get(
    model::Optimizer,
    attr::MOI.VariablePrimal,
    x::MOI.VariableIndex,
)
    MOI.check_result_index_bounds(model, attr)
    return _unsafe_wrap_cbc_array(
        model,
        Cbc_getColSolution,
        Cbc_getNumCols(model),
        x.value,
    )
end

function MOI.get(
    model::Optimizer,
    attr::MOI.VariablePrimal,
    x::Vector{MOI.VariableIndex},
)
    MOI.check_result_index_bounds(model, attr)
    return _unsafe_wrap_cbc_array(
        model,
        Cbc_getColSolution,
        Cbc_getNumCols(model),
        [xi.value for xi in x],
    )
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    index::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},<:Any},
)
    MOI.check_result_index_bounds(model, attr)
    return _unsafe_wrap_cbc_array(
        model,
        Cbc_getRowActivity,
        Cbc_getNumRows(model),
        index.value,
    )
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    index::MOI.ConstraintIndex{MOI.VariableIndex,<:Any},
)
    MOI.check_result_index_bounds(model, attr)
    return MOI.get(model, MOI.VariablePrimal(), MOI.VariableIndex(index.value))
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveSense)
    sense = Cbc_getObjSense(model)
    if sense == 1.0
        return MOI.MIN_SENSE
    elseif sense == -1.0
        return MOI.MAX_SENSE
    else
        @assert sense == 0.0
        return MOI.FEASIBILITY_SENSE
    end
end

struct Status <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::Status) = true

struct SecondaryStatus <: MOI.AbstractModelAttribute end
MOI.is_set_by_optimize(::SecondaryStatus) = true

MOI.get(model::Optimizer, ::Status) = Cbc_status(model)
MOI.get(model::Optimizer, ::SecondaryStatus) = Cbc_secondaryStatus(model)

const _STATUS = Dict{Cint,String}(
    Cint(-1) => "before branchAndBound",
    Cint(0) =>
        "finished - check isProvenOptimal or isProvenInfeasible to see if solution found (or check value of best solution)",
    Cint(1) => "stopped - on maxnodes, maxsols, maxtime",
    Cint(2) => "execution abandoned due to numerical dificulties",
    Cint(5) => "user programmed interruption",
)

const _SECONDARY_STATUS = Dict{Cint,String}(
    Cint(-1) => "unset (status_ will also be -1)",
    Cint(0) => "search completed with solution",
    Cint(1) => "linear relaxation not feasible (or worse than cutoff)",
    Cint(2) => "stopped on gap",
    Cint(3) => "stopped on nodes",
    Cint(4) => "stopped on time",
    Cint(5) => "stopped on user event",
    Cint(6) => "stopped on solutions",
    Cint(7) => "linear relaxation unbounded",
    Cint(8) => "stopped on iteration limit",
)

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    return """
    Cbc_status          = $(_STATUS[model.termination_status])
    Cbc_secondaryStatus = $(_SECONDARY_STATUS[Cbc_secondaryStatus(model)])
    """
end

function MOI.get(model::Optimizer, ::MOI.ResultCount)
    if _CBC_VERSION == v"2.10.3"
        # TODO(odow): Cbc_jll@2.10.3 and the BinaryProvider version shipped in
        # Julia <1.3 contain a patch that is different to upstream. This branch
        # can be removed when we drop support for Julia 1.0 and the 2.10.3 JLL.
        return Cbc_numberSavedSolutions(model) > 0 ? 1 : 0
    end
    if Cbc_getNumIntegers(model) == 0
        # Cbc forwards the solve to the LP solver if there are no integers, so
        # check the termination status for the result count.
        return model.termination_status == 0 ? 1 : 0
    end
    return Cbc_numberSavedSolutions(model) > 0 ? 1 : 0
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    status = model.termination_status
    if status == -1
        return MOI.OPTIMIZE_NOT_CALLED
    elseif Cbc_isProvenOptimal(model) != 0
        return MOI.OPTIMAL
    elseif Cbc_isProvenInfeasible(model) != 0
        if Cbc_getNumIntegers(model) == 0
            # Why Cbc. For LPs, this could mean dual infeasible.
            return MOI.INFEASIBLE_OR_UNBOUNDED
        else
            return MOI.INFEASIBLE
        end
    elseif Cbc_isContinuousUnbounded(model) != 0
        return MOI.INFEASIBLE_OR_UNBOUNDED
    elseif Cbc_isNodeLimitReached(model) != 0
        return MOI.NODE_LIMIT
    elseif Cbc_isSecondsLimitReached(model) != 0
        return MOI.TIME_LIMIT
    elseif Cbc_isSolutionLimitReached(model) != 0
        return MOI.SOLUTION_LIMIT
    elseif status == 0
        return MOI.OTHER_ERROR
    elseif status == 1
        return MOI.OTHER_LIMIT
    elseif status == 2
        return MOI.NUMERICAL_ERROR
    else
        @assert status == 5
        return MOI.INTERRUPTED
    end
end

function MOI.get(model::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index != 1
        return MOI.NO_SOLUTION
    elseif MOI.get(model, MOI.ResultCount()) == 1
        return MOI.FEASIBLE_POINT
    else
        return MOI.NO_SOLUTION
    end
end

function MOI.get(::Optimizer, ::MOI.DualStatus)
    return MOI.NO_SOLUTION
end


text:
function step_residual!(resid, t, dt, uprev, u, f, p, tmpvar, integrator, cache::Tsit5ConstantCache)
    @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7 = cache
    k1 = f(uprev, p, t)
    k2 = f(uprev+dt*a21*k1, p, t+c1*dt)
    k3 = f(uprev+dt*(a31*k1+a32*k2), p, t+c2*dt)
    k4 = f(uprev+dt*(a41*k1+a42*k2+a43*k3), p, t+c3*dt)
    k5 = f(uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4), p, t+c4*dt)
    k6 = f(uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5), p, t+dt)
    @. resid = u - uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)
end

@muladd function step_residual!(resid, t, dt, uprev, u, f, p, tmpvar, integrator, cache::Tsit5Cache)
    @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7 = cache.tab
    @unpack stage_limiter!,step_limiter!,thread = cache
    @unpack k1,k2,k3,k4,k5,k6,tmp = tmpvar
    f(k1, uprev, p, t)
    @.. thread=thread tmp = uprev+dt*a21*k1
    stage_limiter!(tmp, f, p, t+c1*dt)
    f(k2, tmp, p, t+c1*dt)
    @.. thread=thread tmp = uprev+dt*(a31*k1+a32*k2)
    stage_limiter!(tmp, f, p, t+c2*dt)
    f(k3, tmp, p, t+c2*dt)
    @.. thread=thread tmp = uprev+dt*(a41*k1+a42*k2+a43*k3)
    stage_limiter!(tmp, f, p, t+c3*dt)
    f(k4, tmp, p, t+c3*dt)
    @.. thread=thread tmp = uprev+dt*(a51*k1+a52*k2+a53*k3+a54*k4)
    stage_limiter!(tmp, f, p, t+c4*dt)
    f(k5, tmp, p, t+c4*dt)
    @.. thread=thread tmp = uprev+dt*(a61*k1+a62*k2+a63*k3+a64*k4+a65*k5)
    stage_limiter!(tmp, f, p, t+dt)
    f(k6, tmp, p, t+dt)
    @.. thread=thread tmp = uprev+dt*(a71*k1+a72*k2+a73*k3+a74*k4+a75*k5+a76*k6)
    stage_limiter!(tmp, f, p, t+dt)
    step_limiter!(tmp, f, p, t+dt)
    @.. thread=thread resid = u-tmp
end

@muladd function step_residual!(resid, t, dt, uprev, u, f, p, tmpvar, integrator, cache::Tsit5Cache{uType,rateType,uNoUnitsType,TabType,StageLimiter,StepLimiter,Thread}) where {uType<:Union{Array,Zygote.Buffer},rateType,uNoUnitsType,TabType,StageLimiter,StepLimiter,Thread<:False}

    uidx = eachindex(uprev)
    @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7 = cache.tab
    @unpack stage_limiter!,step_limiter! = cache
    @unpack k1,k2,k3,k4,k5,k6,tmp = tmpvar
    f(k1, uprev, p, t)
    @inbounds @simd ivdep for i in uidx
        tmp[i] = uprev[i]+dt*a21*k1[i]
    end
    stage_limiter!(tmp, f, p, t+c1*dt)
    f(k2, tmp, p, t+c1*dt)
    @inbounds @simd ivdep for i in uidx
        tmp[i] = uprev[i]+dt*(a31*k1[i]+a32*k2[i])
    end
    stage_limiter!(tmp, f, p, t+c2*dt)
    f(k3, tmp, p, t+c2*dt)
    @inbounds @simd ivdep for i in uidx
        tmp[i] = uprev[i]+dt*(a41*k1[i]+a42*k2[i]+a43*k3[i])
    end
    stage_limiter!(tmp, f, p, t+c3*dt)
    f(k4, tmp, p, t+c3*dt)
    @inbounds @simd ivdep for i in uidx
        tmp[i] = uprev[i]+dt*(a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i])
    end
    stage_limiter!(tmp, f, p, t+c4*dt)
    f(k5, tmp, p, t+c4*dt)
    @inbounds @simd ivdep for i in uidx
        tmp[i] = uprev[i]+dt*(a61*k1[i]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
    end
    stage_limiter!(tmp, f, p, t+dt)
    f(k6, tmp, p, t+dt)
    @inbounds @simd ivdep for i in uidx
        tmp[i] = uprev[i]+dt*(a71*k1[i]+a72*k2[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i])
    end
    stage_limiter!(tmp, f, p, t+dt)
    step_limiter!(tmp, f, p, t+dt)
    @inbounds @simd ivdep for i in uidx
        resid[i] = u[i]-tmp[i]
    end
end


text:
using Documenter, ImageFeatures

makedocs(sitename = "ImageFeatures",
         format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
         pages    = ["Home" => "index.md",
                     "Framework" => "framework.md",
                     "Tutorials" => [
                         "BRIEF" => "tutorials/brief.md",
                         "ORB" => "tutorials/orb.md",
                         "BRISK" => "tutorials/brisk.md",
                         "FREAK" => "tutorials/freak.md",
                         "Gray level co-occurence matrix" => "tutorials/glcm.md",
                         "Local binary patterns" => "tutorials/lbp.md",
                         "Object Detection using HOG" => "tutorials/object_detection.md"
                     ],
                     "Function reference" => "function_reference.md",
                     ],
         )

deploydocs(repo   = "github.com/JuliaImages/ImageFeatures.jl.git",
           push_preview = true
           )


text:
${code}

${submodule_imports}


text:
#!/usr/bin/env julia

using Luxor, Colors

function randomsine(func, vertscale, horizontalscale, pagewidth, pageheight, style, delta=0.1)
    # y is positive downwards...
    pl = Point[]
    xc=delta
    if rand(Bool)
        for i in -pagewidth/2:delta:pagewidth/2
            push!(pl, Point(i, vertscale * -func(i/horizontalscale) * func(3 * i/horizontalscale) * func(i/horizontalscale)))
        end
    else
        for i in -pagewidth/2:delta:pagewidth/2
            push!(pl, Point(i, vertscale * -func(i / horizontalscale)))
        end
    end
    #push!(pl, Point(pagewidth/2,  pl[end].y))
    #push!(pl, Point(-pagewidth/2, pl[1].y))
    poly(pl, close=false, style)
end

function random_sines(w, h, fname)
    Drawing(w, h, fname)
    origin()
    background("black")
    for i in 1:250
        setline(rand(0.02:0.01:3))
        setdash([
         "solid",
         "dotted",
         "dot",
         "dotdashed",
         "longdashed",
         "shortdashed",
         "dash",
         "dashed",
         "dotdotdashed",
         "dotdotdotdashed"
         ][rand(1:end)])
        gsave()
        sf = rand(1:15)
        scale(sf, sf)
        translate(0, rand(-h/2:h/2))
        randomhue()
        setopacity(rand(0.25:0.1:0.75))
        randomsine(
         [sin, cos][rand(1:end)],                    # random function
         rand(20:40),                               # random vertical scale
         rand(10:20),                                # horizontalscale
         w,                                          # pagewidth
         h,                                          # pageheight
         [:fill, :stroke, :stroke, :fillstroke][rand(1:end)], # random drawing style
         rand(0.1:0.1:0.5))
        grestore()
    end
    finish()
end

fname = "/tmp/random-sines.pdf"
random_sines(2000, 2000, fname)
println("finished test: output in $(fname)")


text:
using BinDeps
using ADCME
using CMake

@BinDeps.setup

const libhypre = library_dependency("libhypre", aliases=["libhypre", "libHYPRE"])
const hyprever = "2.19.0"
provides(Sources,
    URI("https://github.com/kailaix/hypre/archive/v$hyprever.tar.gz"),
    libhypre, unpacked_dir = "hypre-$hyprever")

libhypre.context.dir = ADCME.PREFIXDIR
cmakedir = joinpath(srcdir(libhypre), "hypre-$hyprever", "src", "build")
cc = ADCME.CC
cxx = ADCME.CXX
CONDA_ROOT = joinpath(ADCME.BINDIR, "..")

provides(SimpleBuild,
(@build_steps begin
    GetSources(libhypre)
    CreateDirectory(cmakedir)
    @build_steps begin
        ChangeDirectory(cmakedir)
        `$cmake -DHYPRE_SHARED:BOOL=ON -DHYPRE_INSTALL_PREFIX:PATH=$CONDA_ROOT -DCMAKE_C_COMPILER:FILEPATH=$(cc) -DCMAKE_CXX_COMPILER:FILEPATH=$(cxx) ..`
        `$cmake -L ..`
        MakeTargets(".", ["all"])
        MakeTargets(".", ["install"])
        FileRule(joinpath(ADCME.LIBDIR, "libHYPRE.so"),
            @build_steps begin
                `ln -s $(CONDA_ROOT)/lib64/libHYPRE.so $(ADCME.LIBDIR)/libHYPRE.so`
            end
        )
    end
end),
libhypre,
os = :Linux)

BinDeps.execute(libhypre, BinDeps.SimpleBuild)


text:
import TreeLas
import GraphIdx
import HDF5: h5read, h5write

fname = ARGS[1]
lam = parse(Float64, ARGS[2])
outfn = ARGS[3]
niter = parse(Int, get(ARGS, 4, "150"))
learn = parse(Float64, get(ARGS, 5, "0.78"))
dn = parse(Int, get(ARGS, 6, "1"))

println(fname, " lam ", lam, " => ", outfn)
y = h5read(fname, "y")
graph = GraphIdx.Grid.GridGraph(size(y)..., dn)
edges = GraphIdx.Grid.collect_edges(graph)[1]
x = TreeLas.MGT.gaplas(y, edges, lam, max_iter=niter, learn=learn);
h5write(outfn, "xgrid", x)

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
### Init some initial parameters necessary to create the plan ###

function initParams(x::Matrix{T}, N::NTuple{D,Int}, dims::Union{Integer,UnitRange{Int64}}=1:D;
                    kargs...) where {D,T}
  # convert dims to a unit range
  dims_ = (typeof(dims) <: Integer) ? (dims:dims) : dims

  params = NFFTParams{T,D}(; kargs...)
  m, σ, reltol = accuracyParams(; kargs...)
  params.m = m
  params.σ = σ
  params.reltol = reltol

  # Taken from NFFT3
  m2K = [1, 3, 7, 9, 14, 17, 20, 23, 24]
  K = m2K[min(m+1,length(m2K))]
  params.LUTSize = 2^(K) * (m) # ensure that LUTSize is dividable by (m)

  if length(dims_) != size(x,1)
      throw(ArgumentError("Nodes x have dimension $(size(x,1)) != $(length(dims_))"))
  end

  if any(isodd.(N[dims]))
    throw(ArgumentError("N = $N needs to consist of even integers along dims = $(dims)!"))
  end

  doTrafo = ntuple(d->d ∈ dims_, D)

  n = ntuple(d -> doTrafo[d] ?
                      (ceil(Int,params.σ*N[d])÷2)*2 : # ensure that n is an even integer
                        N[d], D)

  params.σ = n[dims_[1]] / N[dims_[1]]

  #params.blockSize = ntuple(d-> n[d] , D) # just one block
  if haskey(kargs, :blockSize)
    params.blockSize = kargs[:blockSize]
  else
    params.blockSize = ntuple(d-> _blockSize(n,d) , D)
  end

  M = size(x, 2)

  # calculate output size
  NOut = Int[]
  Mtaken = false
  for d=1:D
    if !doTrafo[d]
      push!(NOut, N[d])
    elseif !Mtaken
      push!(NOut, M)
      Mtaken = true
    end
  end
  # Sort nodes in lexicographic way
  if params.sortNodes
      x .= sortslices(x, dims=2)
  end
  return params, N, Tuple(NOut), M, n, dims_
end

function _blockSize(n::NTuple{1,Int}, d)
  return min(1024, n[d])
end

function _blockSize(n::NTuple{2,Int}, d)
  return min(64, n[d])
end

function _blockSize(n::NTuple{D,Int}, d) where {D}
  if d == 1
    return min(16, n[d])
  elseif d == 2
    return min(16, n[d])
  elseif d == 3
    return min(16, n[d])
  else
    return 1
  end
end

### Precomputation of the B matrix ###

function precomputeB(win, x, N::NTuple{D,Int}, n::NTuple{D,Int}, m, M, σ, K, T) where D
  I = Array{Int64,2}(undef, (2*m)^D, M)
  β = (2*m)^D
  J = [β*k+1 for k=0:M]
  V = Array{T,2}(undef,(2*m)^D, M)
  mProd = ntuple(d-> (d==1) ? 1 : (2*m)^(d-1), D)
  nProd = ntuple(d-> (d==1) ? 1 : prod(n[1:(d-1)]), D)
  L = Val(2*m)
  scale = Int(K/(m))

  @cthreads for k in 1:M
    _precomputeB(win, x, N, n, m, M, σ, scale, I, J, V, mProd, nProd, L, k, K)
  end

  S = SparseMatrixCSC(prod(n), M, J, vec(I), vec(V))
  return S
end

@inline @generated function _precomputeB(win, x::AbstractMatrix{T}, N::NTuple{D,Int}, n::NTuple{D,Int}, m, M,
                     σ, scale, I, J, V, mProd, nProd, L::Val{Z}, k, LUTSize) where {T, D, Z}
  quote

    @nexprs $(D) d -> ((tmpIdx_d, tmpWin_d) = precomputeOneNode(win, x, n, m, σ, scale, k, d, L, LUTSize) )

    @nexprs 1 d -> κ_{$D} = 1 # This is a hack, I actually want to write κ_$D = 1
    @nexprs 1 d -> ζ_{$D} = 1
    @nexprs 1 d -> prodWin_{$D} = one(T)
    @nloops $D l d -> 1:$Z d->begin
        # preexpr
        prodWin_{d-1} = prodWin_d * tmpWin_d[l_d]
        κ_{d-1} = κ_d + (l_d-1) * mProd[d]
        ζ_{d-1} = ζ_d + (tmpIdx_d[l_d]-1) * nProd[d]
    end begin
        # bodyexpr
        I[κ_0,k] = ζ_0
        V[κ_0,k] = prodWin_0
    end
    return
  end
end

### precomputation of the window and the indices required during convolution ###

@generated function precomputeOneNode(win::Function, x::AbstractMatrix{T}, n::NTuple{D,Int}, m,
  σ, scale, k, d, L::Val{Z}, LUTSize) where {T,D,Z}
  quote
    xscale = x[d,k] * n[d]
    off = floor(Int, xscale) - m + 1
    tmpIdx = @ntuple $(Z) l -> ( rem(l + off + n[d] - 1, n[d]) + 1)

    tmpWin = @ntuple $(Z) l -> (win( (xscale - (l-1) - off)  / n[d], n[d], m, σ) )
    return (tmpIdx, tmpWin)
  end
end

# precompute = LINEAR

@generated function precomputeOneNode(winLin::Array, winPoly::Nothing, x::AbstractMatrix{T}, n::NTuple{D,Int}, m,
  σ, scale, k, d, L::Val{Z}, LUTSize) where {T,D,Z}
  quote
    xscale = x[d,k] * n[d]
    off = floor(Int, xscale) - m + 1
    tmpIdx = @ntuple $(Z) l -> ( rem(l + off + n[d] - 1, n[d]) + 1)

    idx = ((xscale - off)*LUTSize)/(m)
    tmpWin =  shiftedWindowEntries(winLin, idx, scale, d, L)

    return (tmpIdx, tmpWin)
  end
end

# precompute = POLYNOMIAL

@generated function precomputeOneNode(winLin::Array, winPoly::NTuple{Y, NTuple{X,T}}, x::AbstractMatrix{T}, n::NTuple{D,Int}, m,
  σ, scale, k, d, L::Val{Z}, LUTSize) where {Y,X,T,D,Z}
  quote
    xscale = x[d,k] * n[d]
    off = floor(Int, xscale) - m + 1
    tmpIdx = @ntuple $(Z) l -> ( rem(l + off + n[d] - 1, n[d]) + 1)

    idx = (xscale - off - m + T(0.5))
    tmpWin =  shiftedWindowEntries(winPoly, idx, scale, d, L)

    return (tmpIdx, tmpWin)
  end
end

# precompute = LINEAR

function precomputeOneNodeBlocking(winLin, winTensor::Nothing, winPoly::Nothing,
                                    scale, k, d, L, idxInBlock::Matrix)

  y, idx = idxInBlock[d,k]
  tmpWin =  shiftedWindowEntries(winLin, idx, scale, d, L)

  return (y, tmpWin)
end

@generated function shiftedWindowEntries(winLin::Vector, idx, scale, d, L::Val{Z}) where {Z}
  quote
    idxL = floor(Int,idx)
    idxInt = Int(idxL)
    α = ( idx-idxL )

    tmpWin = @ntuple $(Z) l -> begin
      # Uncommented code: This is the version where we pull in l into the abs.
      # We pulled this out of the iteration.
      # idx = abs((xscale - (l-1)  - off)*LUTSize)/(m)

      # The second +1 is because Julia has 1-based indexing
      # The first +1 is part of the index calculation and needs(!)
      # to be part of the abs. The abs is shifting to the positive interval
      # and this +1 matches the `floor` above, which rounds down. In turn
      # for positive and negative indices a different neighbour is calculated
      idxInt1 = abs( idxInt - (l-1)*scale ) +1
      idxInt2 = abs( idxInt - (l-1)*scale +1) +1

      (winLin[idxInt1] + α * (winLin[idxInt2] - winLin[idxInt1]))
    end
    return tmpWin
  end
end

# precompute = POLYNOMIAL

function precomputeOneNodeBlocking(winLin, winTensor::Nothing, winPoly::NTuple{Z, NTuple{X,T}}, scale,
                                   k, d, L, idxInBlock::Matrix) where {T,Z,X}

  y, x = idxInBlock[d,k]
  tmpWin =  shiftedWindowEntries(winPoly, x, scale, d, L)

  return (y, tmpWin)
end

@generated function shiftedWindowEntries(winPoly::NTuple{Z, NTuple{X,T}}, x::T, scale, d, L::Val{Z}) where {T,Z,X}
  quote
    tmpWin = @ntuple $(Z) l -> begin
      evalpoly(x, winPoly[l])
    end
    return tmpWin
  end
end

#= Static Array version. Not faster
function precomputeOneNodeBlocking(winLin, winTensor::Nothing, winPoly::SMatrix{X,Z,T}, scale,
  k, d, L, idxInBlock::Matrix) where {T,Z,X}

  y, x = idxInBlock[d,k]
  tmpWin =  shiftedWindowEntries(winPoly, x, scale, d, L)

  return (y, tmpWin)
end

@generated function shiftedWindowEntries(winPoly::SMatrix{X,Z,T}, x::T, scale, d, L::Val{Z}) where {T,Z,X}
  quote

    x_1 = one(T)
    @nexprs $(X-1) h->(x_{h+1} = x_h * x)

    xx = @ntuple $(X) h -> begin
      x_{h}
    end

    xs = SVector(xx...)

    tmpWin = transpose(winPoly) * xs

    return tmpWin
  end
end
=#

# precompute = TENSOR

function precomputeOneNodeBlocking(winLin, winTensor::Array, winPoly::Nothing, scale,
                                   k, d, L, idxInBlock::Matrix)
  y, idx = idxInBlock[d,k]
  tmpWin =  shiftedWindowEntriesTensor(winTensor, k, d, L)

  return (y, tmpWin)
end

@generated function shiftedWindowEntriesTensor(winTensor, k, d, L::Val{Z}) where {Z}
  quote
    tmpWin = @ntuple $(Z) l -> begin
      winTensor[l, d, k]
    end
    return tmpWin
  end
end

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

## function nfft_precompute_lin_psi in NFFT3
"""
Precompute the look up table for the window function φ.

Remarks:
* Only the positive half is computed
* The window is computed for the interval [0, (m)/n]. The reason for the +2 is
  that we do evaluate the window function outside its interval, since x does not
  necessary match the sampling points
* The window has K+1 entries and during the index calculation we multiply with the
  factor K/(m).
* It is very important that K/(m) is an integer since our index calculation exploits
  this fact. We therefore always use `Int(K/(m))`instead of `K÷(m)` since this gives
  an error while the later variant would silently error.
"""
function precomputeLinInterp(win, m, σ, K, T)
  windowLinInterp = Vector{T}(undef, K+1)

  step = (m) / (K)
  @cthreads for l = 1:(K+1)
      y = ( (l-1) * step )
      windowLinInterp[l] = win(y, 1, m, σ)
  end
  return windowLinInterp
end

function precomputePolyInterp(win, m, σ, T)
  deg = 2*m # 2*m #+ 2   # Certainly depends on Window
  K = 2*m
  NSamples = 2*deg # Sample more densely!!!
  windowPolyInterp = Matrix{T}(undef, deg, K)

  x = range(-0.5, 0.5,length=NSamples)
  V = ones(deg, NSamples)
  for r=2:deg
    V[r,:] .= V[r-1,:] .* x
  end

  for l = 1:K
      y = (-(l-0.5) + m) .+ x
      samples = win.(y, 1, m, σ)
      windowPolyInterp[:,l] .= V' \ samples
  end
  return windowPolyInterp
end

function testPrecomputePoly(win, m, σ, T)
  deg = 2*m + 3
  K = 2*m
  step = 1

  windowPolyInterp = precomputePolyInterp(win, m, σ, T)

  x = -0.394234
  xx = [x^l for l=0:(deg-1)]

  winTrue = zeros(T,K)
  winApprox = zeros(T,K)

  for l = 1:K
    y = (-(l-0.5) * step + m) + x
    @info y
    winTrue[l] = win(y, 1, m, σ)
    winApprox[l] = windowPolyInterp[:,l]' * xx
  end

  return winTrue, winApprox
end

function precomputeWindowHatInvLUT(windowHatInvLUT, win_hat, N, n, m, σ, T)

  for d=1:length(windowHatInvLUT)
      κ = k -> win_hat(k-1-N[d]÷2, n[d], m, σ)
      cheb = ChebyshevInterpolator(κ, 1, N[d], 30)

      windowHatInvLUT[d] = zeros(T, N[d])
      @cthreads for k=1:N[d]
          windowHatInvLUT[d][k] = 1. / cheb(k) #win_hat(k-1-N[d]÷2, n[d], m, σ)
      end
  end
end

function precomputation(x::Union{Matrix{T},Vector{T}}, N::NTuple{D,Int}, n, params) where {T,D}

  m = params.m; σ = params.σ; window=params.window
  LUTSize = params.LUTSize; precompute = params.precompute

  win, win_hat = getWindow(window) # highly type instable. But what should be do
  M = size(x, 2)

  windowHatInvLUT_ = Vector{Vector{T}}(undef, D)
  precomputeWindowHatInvLUT(windowHatInvLUT_, win_hat, N, n, m, σ, T)

  if params.storeDeconvolutionIdx
    windowHatInvLUT = Vector{Vector{T}}(undef, 1)
    windowHatInvLUT[1], deconvolveIdx = precompWindowHatInvLUT(params, N, n, windowHatInvLUT_)
  else
    windowHatInvLUT = windowHatInvLUT_
    deconvolveIdx = Array{Int64,1}(undef, 0)
  end

  if precompute == LINEAR
    windowLinInterp = precomputeLinInterp(win, m, σ, LUTSize, T)
    windowPolyInterp = Matrix{T}(undef, 0, 0)
    B = sparse([],[],T[])
  elseif precompute == POLYNOMIAL
    windowLinInterp = Vector{T}(undef, 0)
    windowPolyInterp = precomputePolyInterp(win, m, σ, T)
    B = sparse([],[],T[])
  elseif precompute == FULL
    windowLinInterp = Vector{T}(undef, 0)
    windowPolyInterp = Matrix{T}(undef, 0, 0)
    B = precomputeB(win, x, N, n, m, M, σ, LUTSize, T)
    #windowLinInterp = precomputeLinInterp(win, windowLinInterp, n, m, σ, LUTSize, T) # These versions are for debugging
    #B = precomputeB(windowLinInterp, x, N, n, m, M, σ, LUTSize, T)
  elseif precompute == TENSOR
    windowLinInterp = Vector{T}(undef, 0)
    windowPolyInterp = Matrix{T}(undef, 0, 0)
    B = sparse([],[],T[])
  else
    windowLinInterp = Vector{T}(undef, 0)
    windowPolyInterp = Matrix{T}(undef, 0, 0)
    B = sparse([],[],T[])
    error("precompute = $precompute not supported by NFFT.jl!")
  end

  return (windowLinInterp, windowPolyInterp, windowHatInvLUT, deconvolveIdx, B)
end

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

# This function is type unstable. why???
function precompWindowHatInvLUT(p::NFFTParams{T}, N, n, windowHatInvLUT_) where {T}

  windowHatInvLUT = zeros(Complex{T}, N)
  deconvIdx = zeros(Int64, N)

  if length(N) == 1
    precompWindowHatInvLUT(p, windowHatInvLUT, deconvIdx, N, n, windowHatInvLUT_, 1)
  else
    @cthreads for o = 1:N[end]
      precompWindowHatInvLUT(p, windowHatInvLUT, deconvIdx, N, n, windowHatInvLUT_, o)
    end
  end
  return vec(windowHatInvLUT), vec(deconvIdx)
end

@generated function precompWindowHatInvLUT(p::NFFTParams{T}, windowHatInvLUT::AbstractArray{Complex{T},D},
           deconvIdx::AbstractArray{Int,D}, N, n, windowHatInvLUT_, o)::Nothing where {D,T}
  quote
    linIdx = LinearIndices(n)

    @nexprs 1 d -> gidx_{$D-1} = rem(o+n[$D] - N[$D]÷2 - 1, n[$D]) + 1
    @nexprs 1 d -> l_{$D-1} = o
    @nloops $(D-2) l d->(1:N[d+1]) d-> begin
        gidx_d = rem(l_d+n[d+1] - N[d+1]÷2 - 1, n[d+1]) + 1
      end begin
      N2 = N[1]÷2
      @inbounds @simd for i = 1:N2
        deconvIdx[i, CartesianIndex(@ntuple $(D-1) l)] =
           linIdx[i-N2+n[1], CartesianIndex(@ntuple $(D-1) gidx)]
        v = windowHatInvLUT_[1][i]
        @nexprs $(D-1) d -> v *= windowHatInvLUT_[d+1][l_d]
        windowHatInvLUT[i, CartesianIndex(@ntuple $(D-1) l)] = v

        deconvIdx[i+N2, CartesianIndex(@ntuple $(D-1) l)] =
           linIdx[i, CartesianIndex(@ntuple $(D-1) gidx)]
        v = windowHatInvLUT_[1][i+N2]
        @nexprs $(D-1) d -> v *= windowHatInvLUT_[d+1][l_d]
        windowHatInvLUT[i+N2, CartesianIndex(@ntuple $(D-1) l)] = v
      end
    end
    return
  end
end

####### block precomputation #########

function precomputeBlocks(x::Matrix{T}, n::NTuple{D,Int}, params, calcBlocks::Bool) where {T,D}

  if calcBlocks
    xShift = copy(x)
    shiftNodes!(xShift)
    blocks, nodesInBlocks, blockOffsets =
        _precomputeBlocks(xShift, n, params.m, params.LUTSize, params.blockSize)

    idxInBlock =  _precomputeIdxInBlock(xShift, n, params.m, params.precompute, params.LUTSize, blockOffsets, nodesInBlocks)
    if params.precompute != TENSOR
      windowTensor = Array{Array{T,3},D}(undef, ntuple(d->0,D))
    else
      #idxInBlock = Array{Matrix{Tuple{Int,Float64}},D}(undef, ntuple(d->0,D))
      windowTensor = _precomputeWindowTensor(xShift, n, params.m, params.σ, nodesInBlocks, params.window)
    end

  else
    blocks = Array{Array{Complex{T},D},D}(undef,ntuple(d->0,D))
    nodesInBlocks = Array{Vector{Int64},D}(undef,ntuple(d->0,D))
    blockOffsets = Array{NTuple{D,Int64},D}(undef,ntuple(d->0,D))
    idxInBlock = Array{Matrix{Tuple{Int,T}},D}(undef, ntuple(d->0,D))
    windowTensor = Array{Array{T,3},D}(undef, ntuple(d->0,D))
  end

  return (blocks, nodesInBlocks, blockOffsets, idxInBlock, windowTensor)
end

function _precomputeBlocks(x::Matrix{T}, n::NTuple{D,Int}, m, LUTSize, blockSize) where {T,D}

  padding = ntuple(d->m, D)
  numBlocks =  ntuple(d-> ceil(Int, n[d]/blockSize[d]), D)
  blockSizePadded = ntuple(d-> blockSize[d] + 2*padding[d], D)
  nodesInBlock = [ Int[] for l in CartesianIndices(numBlocks) ]
  numNodesInBlock = zeros(Int, numBlocks)
  for k=1:size(x,2) # @cthreads
    idx = ntuple(d->unsafe_trunc(Int, x[d,k]*n[d])÷blockSize[d]+1, D)
    numNodesInBlock[idx...] += 1
  end
  @cthreads  for l in CartesianIndices(numBlocks)
    sizehint!(nodesInBlock[l], numNodesInBlock[l])
  end
  for k=1:size(x,2) # @cthreads
    idx = ntuple(d->unsafe_trunc(Int, x[d,k]*n[d])÷blockSize[d]+1, D)
    push!(nodesInBlock[idx...], k)
  end

  blocks = Array{Array{Complex{T},D},D}(undef, numBlocks)
  blockOffsets = Array{NTuple{D,Int64},D}(undef, numBlocks)

  @cthreads for l in CartesianIndices(numBlocks)
    if !isempty(nodesInBlock[l])
      # precompute blocks
      blocks[l] = Array{Complex{T},D}(undef, blockSizePadded)

      # precompute blockOffsets
      blockOffsets[l] = ntuple(d-> (l[d]-1)*blockSize[d]-padding[d]-1, D)
    end
  end

  return blocks, nodesInBlock, blockOffsets
end

function _precomputeIdxInBlock(x::Matrix{T}, n::NTuple{D,Int}, m, precompute, LUTSize, blockOffsets, nodesInBlock) where {T,D}

  numBlocks = size(nodesInBlock)

  idxInBlock = Array{Matrix{Tuple{Int,T}},D}(undef, numBlocks)

  @cthreads for l in CartesianIndices(numBlocks)
    if !isempty(nodesInBlock[l])

      # precompute idxInBlock
      idxInBlock[l] = Matrix{Tuple{Int,T}}(undef, D, length(nodesInBlock[l]))
      @inbounds for (i,k) in enumerate(nodesInBlock[l])
        @inbounds for d=1:D
          xtmp = x[d,k] # this is expensive because of cache misses
          xscale = xtmp * n[d]
          off = unsafe_trunc(Int, xscale) - m + 1
          y = off - blockOffsets[l][d] - 1

          if precompute == LINEAR
            idx = (xscale - off)*(LUTSize÷(m))
          else
            idx = (xscale - off - m + 1 -0.5 )
          end
          idxInBlock[l][d,i] = (y,idx)
        end
      end
    end
  end

  return idxInBlock
end

function _precomputeWindowTensor(x::Matrix{T}, n::NTuple{D,Int}, m, σ, nodesInBlock, window::Symbol) where {T,D}
  win, win_hat = getWindow(window) # highly type instable. But what should be do
  P = precomputePolyInterp(win, m, σ, T)
  winPoly = ntuple(d-> ntuple(g-> P[g,d], size(P,1)), size(P,2))

  return _precomputeWindowTensor(x, n, m, σ, nodesInBlock, winPoly)
end

function _precomputeWindowTensor(x::Matrix{T}, n::NTuple{D,Int}, m, σ, nodesInBlock, winPoly::NTuple{Z, NTuple{X,T}}) where {T,D,Z,X}

  numBlocks = size(nodesInBlock)
  windowTensor = Array{Array{T,3},D}(undef, numBlocks)

  @cthreads for l in CartesianIndices(numBlocks)
    if !isempty(nodesInBlock[l])

      # precompute idxInBlock
      windowTensor[l] = Array{T,3}(undef, 2*m, D, length(nodesInBlock[l]))
      @inbounds for (i,k) in enumerate(nodesInBlock[l])
        @inbounds for d=1:D
          xtmp = x[d,k]  # this is expensive because of cache misses
          xscale = xtmp * n[d]
          off = unsafe_trunc(Int, xscale) - m + 1
          x_ = (xscale - off - m + 1 -0.5 )

          @inbounds @simd for k=1:Z
            windowTensor[l][k,d,i] = evalpoly(x_, winPoly[k])
          end
        end
      end
    end
  end

  return windowTensor
end