text:
function load(f::File{format"OMETIFF"})
    open(f) do s
        ret = load(s)
    end
end

function load(io::Stream{format"OMETIFF"})
    if io.filename != nothing && !occursin(".ome.tif", io.filename)
        throw(FileIO.LoaderError("OMETIFF", "Not an OME TIFF file!"))
    end

    orig_file = TiffFile(io)
    summary = load_comments(orig_file)

    # load master OME-XML that contains all information about this dataset
    omexml = load_master_xml(orig_file)

    # find all images in this dataset, can either have the Image or Pixel tag
    containers = findall("//*[@DimensionOrder]", omexml)

    pos_names = nodecontent.(findall("/ns:OME/ns:Image/ns:StageLabel[@Name]/@Name", omexml, ["ns"=>namespace(omexml)]))
    # if all position names aren't unique then substitute names
    if length(pos_names) == 0 || !allunique(pos_names)
        pos_names = ["Pos$i" for i in 1:length(containers)]
    end

    master_rawtype = nothing
    mappedtype = Int64
    dimlist = []
    axeslist = []
    n_ifds = 0

    for (idx, container) in enumerate(containers)
        dims, axes_info = build_axes(container)
        push!(dimlist, dims)
        push!(axeslist, axes_info)
        n_ifds += dims[:Z]*dims[:C]*dims[:T]
        rawtype, mappedtype = type_mapping[container["Type"]]

        if master_rawtype == nothing
            master_rawtype = rawtype
        elseif master_rawtype != rawtype
            throw(FileIO.LoaderError("OMETIFF", "Multiple different storage types are not yet support in a multi position image"))
        end
    end

    ifd_index = Array{Union{NTuple{4, Int}, Nothing}}(nothing, n_ifds)
    ifd_files = Array{Union{Tuple{String, String}, Nothing}}(nothing, n_ifds)
    obs_filepaths = Set{String}()
    for (idx, container) in enumerate(containers)
        OMETIFF.ifdindex!(ifd_index, ifd_files, obs_filepaths, container, dimlist[idx], "", idx)
    end
    @assert length(ifd_index) == length(ifd_files)

    files, ifds = get_ifds(orig_file, ifd_files)

    if length(unique(dimlist)) == 1
        master_dims = merge(unique(dimlist)[1], [:P=>length(containers)])
        dimlist = [master_dims]
        masteraxis = copy(axeslist[1])
        masteraxis[6] = Axis{:position}(to_symbol.(pos_names))
        axeslist = [masteraxis]
    else
        throw(FileIO.LoaderError("OMETIFF", "Different sized axes across images is not currently supported."))
    end

    elapsed_times = get_elapsed_times(containers, master_dims, masteraxis)

    img = inmemoryarray(ifd_index, ifds, master_dims, masteraxis, master_rawtype, mappedtype)
    ImageMeta(img, Description=summary, Elapsed_Times=elapsed_times)
end

"""
    dump_omexml(filepath)

Returns the OME-XML embedded inside the OME-TIFF as a prettified string.
"""
function dump_omexml(filepath::String)
    if !endswith(filepath, ".ome.tif")
        error("Passed file is not an OME-TIFF")
    end
    io = IOBuffer()
    open(filepath) do f
        s = Stream(format"OMETIFF", f, OMETIFF.extract_filename(f))
        orig_file = OMETIFF.TiffFile(s)
        omexml = OMETIFF.load_master_xml(orig_file)
        prettyprint(io, omexml)
    end
    String(take!(io))
end

"""
Builds an in-memory high-dimensional image from the list of IFDs, `ifds`, and
the corresponding indices, `ifd_index`, in the high-dimensional array.
"""
function inmemoryarray(ifd_index::Array{Union{NTuple{4, Int}, Nothing}},
                       ifds::Array{Union{IFD, Nothing}},
                       master_dims::NamedTuple,
                       masteraxis::Array{Axis},
                       master_rawtype::Type,
                       mappedtype::Type)

    data = Array{master_rawtype, length(master_dims)}(undef, master_dims...)

    # iterate over each IFD
    for idx in 1:length(ifd_index)
        indices = ifd_index[idx]
        ifd = ifds[idx]
        (ifd == nothing) && continue

        width, height = master_dims[:X], master_dims[:Y]

        n_strips = length(ifd.strip_offsets)
        strip_len = floor(Int, (width * height) / n_strips)
        read_dims = n_strips > 1 ? (strip_len) : (height, width)

        # TODO: This shouldn't be allocated for each ifd
        tmp = Array{master_rawtype}(undef, read_dims...)
        for j in 1:n_strips
            seek(ifd.file.io, ifd.strip_offsets[j])
            read!(ifd.file.io, tmp)
            do_bswap(ifd.file, tmp)
            if n_strips > 1
                data[j, :, indices...] = tmp
            else
                data[:, :, indices...] = tmp'
            end
        end
    end

    data = reinterpret(Gray{mappedtype}, data)
    unused_dims = Tuple(idx for (idx, key) in enumerate(keys(master_dims)) if master_dims[idx] == 1)
    squeezed_data = dropdims(data; dims=unused_dims)
    used_axes = [masteraxis[i] for i in 1:length(masteraxis) if !(i in unused_dims)]
    AxisArray(squeezed_data, used_axes...)
end

"""Corresponding Julian types for OME-XML types"""
type_mapping = Dict(
    "uint8" => (UInt8, N0f8),
    "uint16" => (UInt16, N0f16),
    "uint32" => (UInt32, N0f32),
    "float" => (Float32, Float32),
    "double" => (Float64, Float64),
    "int8" => (Int8, N0f8)
)


text:
function ___clustering(neurons,xx_clust,max_iter)

    ##### A: standard k-means
    iseeds=Vector{Int}(undef, neurons)
    # initseeds!(iseeds, KmppAlg(), xx_clust')
    # initseeds!(iseeds, KmCentralityAlg(), xx_clust')
    initseeds!(iseeds, RandSeedAlg(), xx_clust')
    R = kmeans(xx_clust', neurons; maxiter=max_iter, display=:iter,
                init=iseeds) #   none
    a = assignments(R); c = counts(R); minimum(c)
    i1=sortperm(a);inds_all=copy(i1);items_per_neuron=copy(c);n_per_part=[0;cumsum(items_per_neuron)]
    ##### B: Enhanced Algorithm for equally sized clusters and invertible matrices (slower)
    # include("cluster_tests1.jl")
    return inds_all,n_per_part,items_per_neuron
end

ext: jl


text:
using Yao
using YaoAD
using LinearAlgebra
using QuAlgorithmZoo: heisenberg, random_diff_circuit, pair_ring

# generate a Heisenberg Model Hamiltonian
nbit = 4
h = mat(heisenberg(nbit))
v0 = statevec(zero_state(nbit))
function energy(circuit)
    v = mat(circuit) * v0
    (v'* h * v)[] |> real
end

# Generate a circuit as a wave function ansatz
circuit = random_diff_circuit(nbit, 2, pair_ring(nbit))

using Flux: ADAM, Optimise
function train!(lossfunc, circuit, optimizer; maxiter::Int=200)
    dispatch!(circuit, :random)
    params = parameters(circuit)
    loss_history = Float64[]
    for i = 1:maxiter
        # collect gradients from returned structured data
        grad = collect_gradients(lossfunc'(circuit))
        Optimise.update!(optimizer, params, grad)
        dispatch!(circuit, params)
        eng = lossfunc(circuit)
        push!(loss_history, eng)
        println("Iter $i, Energy (Loss) = $(eng)")
    end
    loss_history
end

using Random
Random.seed!(5)
EG = eigvals(Matrix(h))[1]
println("$nbit site Heisenberg model, exact ground state energy = $EG")
loss_history = train!(energy, circuit, ADAM(0.1); maxiter=200)

using DelimitedFiles
writedlm("data/loss_history_$nbit.dat", loss_history)


text:
## as anonymous function
@model StateSpaceModel(y) = begin
    N = length(y)
    x = tzeros(Real, N)
    σy = 1.0
    σx = 1.0

    x[1] ~ Normal(0, 1)
    y[1] ~ Normal(x[1], σy)

    # Iterators.drop(1:N, 1)
    for i ∈ 2:N
        x[i] ~ Normal(x[i-1], σx)
        y[i] ~ Normal(x[i], σy)
    end
end

## as function
    # @model function StateSpaceModel(y)
    #     N = length(y)
    #     x = tzeros(Real, N)
    #     σy = 1.0
    #     σx = 1.0

    #     x[1] ~ Normal(0, 1)
    #     y[1] ~ Normal(x[1], σy)
    #     for i ∈ 2:N
    #         x[i] ~ Normal(x[i-1], σx)
    #         y[i] ~ Normal(x[i-1], σy)
    #     end
    # end

@model RandomWalkAngularAccelerationModel(y, s; σω = 1.0) = begin
    N = length(y)
    α = tzeros(Real, N)
    ω = tzeros(Real, N)
    c = zeros(Real, N)
    θ = zeros(Real, N)
    Δt = s.Δt
    κ1, κ2 = cupula_dynamics(Δt)

    σα = σω
    # σω = 0.1
    σy = s.sensor.σ

    α[1] ~ Normal(0, σα)
    ω[1] ~ Normal(0, σω)
    c[1] ~ Determin(κ1 * 0 + κ2 * ω[1])
    y[1] ~ Normal(ω[1] - c[1], σy)

    for i ∈ Iterators.drop(1:N, 1)
        α[i] ~ Normal(α[i-1], σα)
        ω[i] ~ Normal(ω[i-1] + Δt * α[i], σω)
        c[i] ~ Determin(κ1 * c[i-1] + κ2 * ω[i])
        y[i] ~ Normal(ω[i] - c[i], σy)
    end
end

@model RandomWalkAngularVelocityModel(y, s; σω = 1.0) = begin
    N = length(y)
    ω = tzeros(Real, N)
    p = tzeros(Real, N)
    c = zeros(Real, N)
    Δt = s.Δt
    κ1, κ2 = cupula_dynamics(Δt)

    # σω = 1
    σy = s.sensor.σ
    ω[1] ~ Normal(0, σω)
    p[1] ~ Determin(ω[1] > 0 ? 1 : 0)
    c[1] ~ Determin(κ1 * 0 + κ2 * ω[1])
    y[1] ~ Normal(ω[1] - c[1], σy)

    for i ∈ Iterators.drop(1:N, 1)
        ω[i] ~ Normal(ω[i-1], σω)
        p[i] ~ Determin(ω[i] > 0 ? 1 : 0)
        c[i] ~ Determin(κ1 * c[i-1] + κ2 * ω[i])
        y[i] ~ Normal(ω[i] - c[i], σy)
    end
end

@model HeadTurnModel(y, s, m, ::Type{T} = Float64; σωᵤ = 0.01, σωₑ = 0.1) where {T} = begin
    N = length(y)
    stochastic = s.stochastic
    # intialize vectors
    α1 = TArray{T}(undef, N)
    ω1 = TArray{T}(undef, N)
    p1 = TArray{T}(undef, N)
    ω2 = TArray{T}(undef, N)
    p2 = TArray{T}(undef, N)
    ω = TArray{T}(undef, N)
    p = TArray{T}(undef, N)
    c = TArray{T}(undef, N)

    # parameters
    D1 = m.D == "left" ? -1 : 1
    A1 = stochastic == 0 ? m.A : rand(Normal(m.A, 0.5))
    # A1 = rand(Normal(m.A, 0.5))
    # A1 ~ Normal(m.A, 0.1)

    Δt = s.Δt
    onsetᵤ = hasproperty(m, :onset) ? m.onset : s.onsetᵤ
    duration1 = m.duration
    f1 = 1/(duration1) # frequency: single sinusiodal head turn
    timesteps = range(0, stop = N, step = Δt)
    κ1, κ2 = cupula_dynamics(Δt)

    # priors
    σy = 1.1 * s.sensor.σ # jitter the measurement noise
    # σωᵤ ~ Truncated(Normal(0.5, 0.1))
    # σωₑ ~ Truncated(Normal(0.5, 0.1))

    # initial state at t_0
    α1[1] = 0.0
    ω1[1] ~ Normal(0, σωᵤ)
    # p1[1] ~ Determin(ω1[1] > 0 ? 1 : 0)
    ω2[1] ~ Normal(0, σωₑ)
    # p2[1] ~ Determin(ω2[1] > 0 ? 1 : 0)
    ω[1] ~ Determin(ω1[1] + ω2[1])
    # p[1] ~ Determin(ω[1] > 0 ? 1 : 0)
    c[1] ~ Determin(κ1 * 0 + κ2 * ω[1])
    y[1] ~ Normal(ω[1] - c[1], σy)

    # for i ∈ Iterators.drop(1:N, 1)
    for i ∈ 2:N
        t = round(timesteps[i]; digits = 3)
        α1[i] = (t > onsetᵤ) & (t < onsetᵤ + duration1) ? acceleration(D1, A1, f1, t, onsetᵤ, noisy = stochastic) : 0.0

        ω1[i] ~ Normal(ω1[i-1] + Δt * α1[i], σωᵤ)
        # p1[i] ~ Determin(ω1[i] > 0 ? 1 : 0)
        ω2[i] ~ Normal(ω2[i-1], σωₑ)
        # p2[i] ~ Determin(ω2[i] > 0 ? 1 : 0)
        ω[i] ~ Determin(ω1[i] + ω2[i])
        # p[i] ~ Determin(ω[i] > 0 ? 1 : 0)

        c[i] ~ Determin(κ1 * c[i-1] + κ2 * ω[i])
        y[i] ~ Normal(ω[i] - c[i], σy)
    end
end

@model HeadTurnOnlyModel(y, s, m; σωᵤ = 0.01) = begin
    N = length(y)
    stochastic = s.stochastic
    # intialize vectors
    α1 = tzeros(Real, N)
    ω = tzeros(Real, N)
    p = tzeros(Real, N)
    c = zeros(Real, N)

    # parameters
    D1 = m.D == "left" ? -1 : 1
    # A1 = stochastic == 0 ? m.A : rand(Normal(m.A, 0.5))
    A1 = rand(Normal(m.A, 0.5))

    Δt = s.Δt
    onsetᵤ = hasproperty(m, :onset) ? m.onset : s.onsetᵤ
    duration1 = m.duration
    f1 = 1/(duration1) # frequency: single sinusiodal head turn
    timesteps = range(0, stop = N, step = Δt)
    κ1, κ2 = cupula_dynamics(Δt)

    # priors
    σy = s.sensor.σ

    # initial state at t_0
    ω[1] ~ Normal(0, σωᵤ)
    p[1] ~ Determin(ω[1] > 0 ? 1 : 0)
    c[1] ~ Determin(κ1 * 0 + κ2 * ω[1])
    y[1] ~ Normal(ω[1] - c[1], σy)

    for i ∈ Iterators.drop(1:N, 1)
        t = round(timesteps[i]; digits = 3)
        α1[i] = (t > onsetᵤ) & (t < onsetᵤ + duration1) ? acceleration(D1, A1, f1, t, onsetᵤ, noisy = stochastic) : 0.0

        ω[i] ~ Normal(ω[i-1] + Δt * α1[i], σωᵤ)
        p[i] ~ Determin(ω[i] > 0 ? 1 : 0)
        c[i] ~ Determin(κ1 * c[i-1] + κ2 * ω[i])
        y[i] ~ Normal(ω[i] - c[i], σy)
    end
end

@model HeadTurnOnlyParamModel(y, s, m; σωᵤ = 0.01, σₐ = 1.0) = begin
    N = length(y)

    # intialize vectors
    α1 = tzeros(Real, N)
    ω = tzeros(Real, N)
    p = tzeros(Real, N)
    c = zeros(Real, N)

    # parameters
    D1 = m.D == "left" ? -1 : 1
    A1 ~ Normal(m.A, σₐ)

    Δt = s.Δt
    onsetᵤ = hasproperty(m, :onset) ? m.onset : s.onsetᵤ
    duration1 = m.duration
    f1 = 1/(duration1) # frequency: single sinusiodal head turn
    timesteps = range(0, stop = N, step = Δt)
    κ1, κ2 = cupula_dynamics(Δt)

    # priors
    σy = s.sensor.σ

    # initial state at t_0
    ω[1] ~ Normal(0, σωᵤ)
    # p[1] ~ Determin(ω[1] > 0 ? 1 : 0)
    c[1] ~ Determin(κ1 * 0 + κ2 * ω[1])
    y[1] ~ Normal(ω[1] - c[1], σy)

    for i ∈ Iterators.drop(1:N, 1)
        t = round(timesteps[i]; digits = 3)
        α1[i] = (t > onsetᵤ) & (t < onsetᵤ + duration1) ? acceleration(D1, A1, f1, t, onsetᵤ, noisy = false) : 0.0

        ω[i] ~ Normal(ω[i-1] + Δt * α1[i], σωᵤ)
        # p[i] ~ Determin(ω[i] > 0 ? 1 : 0)
        c[i] ~ Determin(κ1 * c[i-1] + κ2 * ω[i])
        y[i] ~ Normal(ω[i] - c[i], σy)
    end
end

@model HeadTurnOnlySimulationModel(y, s, m; σωᵤ = 0.01, σₐ = 1.0) = begin
    N = length(y)

    # intialize vectors
    α1 = tzeros(Real, N)
    ω = tzeros(Real, N)
    p = tzeros(Real, N)
    c = zeros(Real, N)

    # parameters
    D = m.D == "left" ? -1 : 1
    A ~ Normal(m.A, σₐ)

    Δt = s.Δt
    on = hasproperty(m, :onset) ? m.onset : s.onsetᵤ
    # on = shape * scale
    onsetᵤ ~ Gamma(20 * on, 0.05)
    duration1 = m.duration
    f = 1/(duration1) # frequency: single sinusiodal head turn
    timesteps = range(0, stop = N, step = Δt)
    κ1, κ2 = cupula_dynamics(Δt)

    # priors
    σy = s.sensor.σ

    # initial state at t_0
    ω[1] ~ Normal(0, σωᵤ)
    # p[1] ~ Determin(ω[1] > 0 ? 1 : 0)
    c[1] ~ Determin(κ1 * 0 + κ2 * ω[1])
    y[1] ~ Normal(ω[1] - c[1], σy)

    for i ∈ Iterators.drop(1:N, 1)
        t = round(timesteps[i]; digits = 3)
        α1[i] = (t > onsetᵤ) & (t < onsetᵤ + duration1) ? acceleration(D, A, f, t, onsetᵤ, noisy = false) : 0.0

        ω[i] ~ Normal(ω[i-1] + Δt * α1[i], σωᵤ)
        # p[i] ~ Determin(ω[i] > 0 ? 1 : 0)
        c[i] ~ Determin(κ1 * c[i-1] + κ2 * ω[i])
        y[i] ~ Normal(ω[i] - c[i], σy)
    end
end

# TODO: when moving, set σωᵤ much lower, when not moving higher
@model HeadTurnComplexModel(y, s, m; σωᵤ = 0.01, σωₑ = 0.1) = begin
    N = length(y)
    stochastic = s.stochastic
    # intialize vectors
    α1 = tzeros(Real, N)
    ω1 = tzeros(Real, N)
    p1 = tzeros(Real, N)
    ω2 = tzeros(Real, N)
    p2 = tzeros(Real, N)
    ω = tzeros(Real, N)
    move = tzeros(Real, N)
    p = tzeros(Real, N)
    c = zeros(Real, N)

    # parameters
    D1 = m.D == "left" ? -1 : 1
    A1 = stochastic == 0 ? m.A : rand(Normal(m1.A, 0.5))

    Δt = s.Δt
    onsetᵤ = hasproperty(m, :onset) ? m.onset : s.onsetᵤ
    duration1 = m.duration
    f1 = 1/(duration1) # frequency: single sinusiodal head turn
    timesteps = range(0, stop = N, step = Δt)
    κ1, κ2 = cupula_dynamics(Δt)

    # priors
    σy = s.sensor.σ

    # initial state at t_0
    ω1[1] ~ Normal(0, σωᵤ)
    p1[1] ~ Determin(ω1[1] > 0 ? 1 : 0)
    ω2[1] ~ Normal(0, σωₑ)
    p2[1] ~ Determin(ω2[1] > 0 ? 1 : 0)
    ω[1] ~ Determin(ω1[1] + ω2[1])
    move[1] ~ Determin(0)
    p[1] ~ Determin(ω[1] > 0 ? 1 : 0)
    c[1] ~ Determin(κ1 * 0 + κ2 * ω[1])
    y[1] ~ Normal(ω[1] - c[1], σy)

    for i ∈ Iterators.drop(1:N, 1)
        t = round(timesteps[i]; digits = 3)
        α1[i] = (t > onsetᵤ) & (t < onsetᵤ + duration1) ? acceleration(D1, A1, f1, t, onsetᵤ, noisy = stochastic) : 0.0

        ω1[i] ~ Normal(ω1[i-1] + Δt * α1[i], σωᵤ)
        move[i] ~ Determin((t > onsetᵤ) & (t < onsetᵤ + duration1) ? 1 : 0)
        p1[i] ~ Determin(ω1[i] > 0 ? 1 : 0)
        ω2[i] ~ Normal(ω2[i-1], σωₑ)
        p2[i] ~ Determin(ω2[i] > 0 ? 1 : 0)
        ω[i] ~ Determin(ω1[i] + ω2[i])
        p[i] ~ Determin(ω[i] > 0 ? 1 : 0)

        c[i] ~ Determin(κ1 * c[i-1] + κ2 * ω[i])
        y[i] ~ Normal(ω[i] - c[i], σy)
    end
end

@model CoupledModel(y, s, m; σωᵤ = 0.01, σωₑ = 0.1) = begin
    N = length(y)
    stochastic = s.stochastic
    # intialize vectors
    α1 = tzeros(Real, N)
    ω1 = tzeros(Real, N)
    # p1 = tzeros(Real, N)
    ω2 = tzeros(Real, N)
    p2 = tzeros(Real, N)
    ω = tzeros(Real, N)
    # p = tzeros(Real, N)
    c = zeros(Real, N)

    # parameters
    D1 = m.D == "left" ? -1 : 1
    A1 = stochastic == 0 ? m.A : rand(Normal(m1.A, 0.5))

    Δt = s.Δt
    onsetᵤ = hasproperty(m, :onset) ? m.onset : s.onsetᵤ
    duration1 = m.duration
    f1 = 1/(duration1) # frequency: single sinusiodal head turn
    timesteps = range(0, stop = N, step = Δt)
    κ1, κ2 = cupula_dynamics(Δt)

    # priors
    σy = s.sensor.σ

    # initial state at t_0
    ω1[1] ~ Normal(0, σωᵤ)
    # p1[1] ~ Determin(ω1[1] > 0 ? 1 : 0)
    ω2[1] ~ Normal(0, σωₑ)
    # p2[1] ~ Determin(ω2[1] > 0 ? 1 : 0)
    ω[1] ~ Determin(ω1[1] - ω2[1])
    # p[1] ~ Determin(ω[1] > 0 ? 1 : 0)
    c[1] ~ Determin(κ1 * 0 + κ2 * ω[1])
    y[1] ~ Normal(ω1[1] - c[1], σy)
    y[1] ~ Normal(ω2[1] - c[1], σy)

    for i ∈ Iterators.drop(1:N, 1)
        t = round(timesteps[i]; digits = 3)
        α1[i] = (t > onsetᵤ) & (t < onsetᵤ + duration1) ? acceleration(D1, A1, f1, t, onsetᵤ, noisy = stochastic) : 0.0

        ω1[i] ~ Normal(ω1[i-1] + Δt * α1[i], σωᵤ)
        # p1[i] ~ Determin(ω1[i] > 0 ? 1 : 0)
        ω2[i] ~ Normal(ω2[i-1], σωₑ)
        # p2[i] ~ Determin(ω2[i] > 0 ? 1 : 0)
        ω[i] ~ Determin(ω1[i] - ω2[i])
        # p[i] ~ Determin(ω[i] > 0 ? 1 : 0)

        c[i] ~ Determin(κ1 * c[i-1] + κ2 * ω[i])
        y[i] ~ Normal(ω1[i] - c[i], σy)
        y[i] ~ Normal(ω2[i] - c[i], σy)
    end
end


text:
# just use the real part for LS.
#yr, yi assumed to have LS_inds applied already.
# This does A\b but with bounds.
function solveBLScL( A::Matrix{T},
    b,
    x_lower::Vector{T},
    x_upper::Vector{T}) where T <: Real

    # This is the formula.
    # BoundedLeastSquares.Quadratic(A'*A, A'*b)

    AtA = forcesymmetric(A'*A)
    Aty = A'*b

    Q = BoundedLeastSquares.Quadratic(AtA, Aty)

    x_opt, out_flag = BoundedLeastSquares.min_bound_constrained_quadratic(Q,
                        x_lower, x_upper)

    return x_opt, out_flag
end

####### w.

function setupupdatew(N_positions::Int, N_compounds::Int)

    T = Float64

    # q_t.
    A = Matrix{T}(undef, 2*N_positions, N_compounds)
    w = zeros(T, N_compounds)

    return A, w
end

function updatew!(  A::Matrix{T},
    b::Vector{T},
    w::Vector{T},
    U_rad_LS,
    Es::Vector{NMRSpectraSimulator.FIDModelType{T,SST}},
    w_lower::Vector{T},
    w_upper::Vector{T}) where {T <: Real, SST}

    #@assert size(A) == (length(b), length(αS))
    @assert length(b) == 2*length(U_rad_LS)

    evaldesignmatrixw!(A, U_LS, Es)

    ### LS solve.
    # w[:] = NonNegLeastSquares.nonneg_lsq(A, b; alg = :fnnls)
    # status_flag = true

    w[:], status_flag = solveBLScL(A, b, w_lower, w_upper)

    return status_flag
end

function evaldesignmatrixw!(B::Matrix{T},
    U,
    Es::Vector{NMRSpectraSimulator.FIDModelType{T,SST}}) where {T <: Real, SST}

    #
    M = length(U)
    N = length(Es)

    @assert size(B) == (2*M,N)

    fill!(B, NaN) # debug.

    for n = 1:N
        for m = 1:M
            # speed up later.
            tmp = NMRSpectraSimulator.evalitpproxycompound(U[m], Es[n])
            B[m,n], B[m+M,n] = real(tmp), imag(tmp)
        end
    end

    return nothing
end

######### κ, compensation for α.

function updateκ!(  A::Matrix{T},
    #r_buffer::Vector{T},
    b::Vector{T},
    κ::Vector{T},
    U_rad_LS,
    Es::Vector{NMRSpectraSimulator.καFIDModelType{T,SST}},
    w::Vector{T},
    κ_lower::Vector{T},
    κ_upper::Vector{T}) where {T <: Real,SST}

    #@assert size(A) == (length(b), length(αS))
    @assert length(b) == 2*length(U_rad_LS)

    evaldesignmatrixκ!(A, U_rad_LS, Es, w)

    ### LS solve.
    # w[:] = NonNegLeastSquares.nonneg_lsq(A, b; alg = :fnnls)
    # status_flag = true

    κ[:], status_flag = solveBLScL(A, b, κ_lower, κ_upper)
    # κ[:] = A\b
    # clamp!(κ, κ_lower[1], κ_upper[1])
    # status_flag = true

    return status_flag
end

function evaldesignmatrixκ!(B::Matrix{T},
    U_rad,
    Es::Vector{NMRSpectraSimulator.καFIDModelType{T,NMRSpectraSimulator.SpinSysParamsType1{T}}},
    w::Vector{T}) where T <: Real

    #
    M = length(U_rad)
    N = length(Es)

    N_κ, N_κ_singlets = countκ(Es)
    #println((N_κ, N_κ_singlets))
    #println(size(B))
    @assert size(B) == (2*M, N_κ + N_κ_singlets)
    #fill!(B, Inf) # debug.

    resetκ!(Es)
    j = 0

    # loop over each κ partition element in Es.
    for n = 1:N
        A = Es[n]
        @assert length(A.κs_α) == length(A.core.qs)

        # spin system.
        for i = 1:length(A.κs_α)

            # partition
            for k = 1:length(A.κs_α[i])

                j += 1

                # loop over each fit position.
                for m = 1:M

                    # taken from evalitproxysys()
                    r = U_rad[m] - A.core.ss_params.d[i]

                    out = w[n]*A.core.qs[i][k](r, A.core.ss_params.κs_λ[i])

                    B[m,j], B[m+M,j] = real(out), imag(out)
                end
            end
        end

        # singlets.
        for k = 1:length(A.κs_α_singlets)
            j += 1

            for m = 1:M

                tmp = w[n]*NMRSpectraSimulator.evalκsinglets(U_rad[m], A.core.d_singlets,
                A.core.αs_singlets, A.core.Ωs_singlets,
                A.core.β_singlets, A.core.λ0, A.core.κs_λ_singlets)

                B[m,j] = real(tmp)
                B[m+M,j] = imag(tmp)
            end
        end

    end

    return nothing
end

function evaldesignmatrixκ!(B::Matrix{T},
    U_rad,
    Es::Vector{NMRSpectraSimulator.καFIDModelType{T,NMRSpectraSimulator.SpinSysParamsType2{T}}},
    w::Vector{T}) where T <: Real

    #
    M = length(U_rad)
    N = length(Es)

    N_κ, N_κ_singlets = countκ(Es)
    #println((N_κ, N_κ_singlets))
    #println(size(B))
    @assert size(B) == (2*M, N_κ + N_κ_singlets)
    fill!(B, Inf) # debug.

    resetκ!(Es)
    j = 0

    #resize!(r_buffer, M)

    # loop over each κ partition element in Es.
    for n = 1:N
        A = Es[n]
        @assert length(A.κs_α) == length(A.core.qs)

        # spin system.
        for i = 1:length(A.κs_α)

            # partition
            for k = 1:length(A.κs_α[i])

                j += 1

                # loop over each fit position.
                for m = 1:M

                    # taken from evalitproxysys()
                    r = U_rad[m] - A.core.ss_params.d[i][k]
                    out = w[n]*A.core.qs[i][k](r, A.core.ss_params.κs_λ[i][k])

                    B[m,j], B[m+M,j] = real(out), imag(out)
                end
            end
        end

        # singlets.
        for k = 1:length(A.κs_α_singlets)
            j += 1

            for m = 1:M

                tmp = w[n]*NMRSpectraSimulator.evalκsinglets(U_rad[m], A.core.d_singlets,
                A.core.αs_singlets, A.core.Ωs_singlets,
                A.core.β_singlets, A.core.λ0, A.core.κs_λ_singlets)

                B[m,j] = real(tmp)
                B[m+M,j] = imag(tmp)
            end
        end

    end

    return nothing
end


text:
import Base: rand,
             ∈

export VPolygon,
       remove_redundant_vertices,
       remove_redundant_vertices!,
       convex_hull,
       linear_map,
       minkowski_sum

"""
    VPolygon{N<:Real, VN<:AbstractVector{N}} <: AbstractPolygon{N}

Type that represents a polygon by its vertices.

### Fields

- `vertices` -- the list of vertices

### Notes

The constructor of `VPolygon` runs a convex hull algorithm on its vertices by
default, to remove the possibly redundant vertices. The vertices are sorted in
counter-clockwise fashion. Use the flag `apply_convex_hull=false` to skip the
computation of the convex hull.

- `VPolygon(vertices::Vector{Vector{N}};
            apply_convex_hull::Bool=true,
            algorithm::String="monotone_chain")`

### Examples

A polygon in vertex representation can be constructed by passing the list of
vertices. For example, we can build the right triangle

```jldoctest polygon_vrep
julia> P = VPolygon([[0, 0], [1, 0], [0, 1]])
VPolygon{Int64,Array{Int64,1}}(Array{Int64,1}[[0, 0], [1, 0], [0, 1]])
```

Alternatively, a `VPolygon` can be constructed passing a matrix of vertices,
where each *column* represents a vertex:

```jldoctest polygon_vrep
julia> M = [0 1 0; 0 0 1.]
2×3 Array{Float64,2}:
 0.0  1.0  0.0
 0.0  0.0  1.0

julia> VPolygon(M)
VPolygon{Float64,Array{Float64,1}}(Array{Float64,1}[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]])
```
"""
struct VPolygon{N<:Real, VN<:AbstractVector{N}} <: AbstractPolygon{N}
    vertices::Vector{VN}

    # default constructor that applies a convex hull algorithm
    function VPolygon(vertices::Vector{VN};
                      apply_convex_hull::Bool=true,
                      algorithm::String="monotone_chain"
                     ) where {N<:Real, VN<:AbstractVector{N}}
        if apply_convex_hull
            return new{N, VN}(convex_hull(vertices, algorithm=algorithm))
        else
            return new{N, VN}(vertices)
        end
    end
end

isoperationtype(::Type{<:VPolygon}) = false
isconvextype(::Type{<:VPolygon}) = true

# constructor with empty vertices list
VPolygon{N}() where {N<:Real} = VPolygon(Vector{Vector{N}}(), apply_convex_hull=false)

# constructor with no vertices of type Float64
VPolygon() = VPolygon{Float64}()

# constructor from rectangular matrix
function VPolygon(vertices_matrix::MT; apply_convex_hull::Bool=true,
                  algorithm::String="monotone_chain") where {N<:Real, MT<:AbstractMatrix{N}}
    @assert size(vertices_matrix, 1) == 2 "the number of rows of the matrix of vertices " *
                                           "should be 2, but it is $(size(vertices_matrix, 1))"

    vertices = [vertices_matrix[:, j] for j in 1:size(vertices_matrix, 2)]
    return VPolygon(vertices; apply_convex_hull=apply_convex_hull, algorithm=algorithm)
end

"""
    remove_redundant_vertices!(P::VPolygon{N};
                               [algorithm]::String="monotone_chain") where {N<:Real}

Remove the redundant vertices of the given polygon.

### Input

- `P`         -- polygon in vertex representation
- `algorithm` -- (optional, default: "monotone_chain") the algorithm used to
                 compute the convex hull

### Output

A new polygon such that its vertices are the convex hull of the given polygon.

### Algorithm

A convex hull algorithm is used to compute the convex hull of the vertices of the
given input polygon `P`; see `?convex_hull` for details on the available algorithms.
The vertices of the output polygon are sorted in counter-clockwise fashion.
"""
function remove_redundant_vertices!(P::VPolygon{N};
                                    algorithm::String="monotone_chain") where {N<:Real}
    convex_hull!(P.vertices; algorithm=algorithm)
end

"""
    remove_redundant_vertices(P::VPolygon{N};
                              [algorithm]::String="monotone_chain") where {N<:Real}

Return the polygon obtained by removing the redundant vertices of the given polygon.

### Input

- `P`         -- polygon in vertex representation
- `algorithm` -- (optional, default: "monotone_chain") the algorithm used to
                 compute the convex hull

### Output

A new polygon such that its vertices are the convex hull of the given polygon.

### Algorithm

A convex hull algorithm is used to compute the convex hull of the vertices of the
given input polygon `P`; see `?convex_hull` for details on the available algorithms.
The vertices of the output polygon are sorted in counter-clockwise fashion.
"""
function remove_redundant_vertices(P::VPolygon{N};
                                   algorithm::String="monotone_chain") where {N<:Real}
    return remove_redundant_vertices!(copy(P), algorithm=algorithm)
end

function linear_map(M::AbstractMatrix{N}, P::VPolygon{N}) where {N<:Real}
    @assert size(M, 2) == 2 "a linear map of size $(size(M)) cannot be applied to a set of dimension 2"
    return _linear_map_vrep(M, P)
end

@inline function _linear_map_vrep(M::AbstractMatrix{N}, P::VPolygon{N}) where {N<:Real}
    return broadcast(v -> M * v, vertices_list(P)) |> VPolygon
end

# --- AbstractPolygon interface functions ---

"""
    tovrep(P::VPolygon{N}) where {N<:Real}

Build a vertex representation of the given polygon.

### Input

- `P` -- polygon in vertex representation

### Output

The identity, i.e., the same polygon instance.
"""
function tovrep(P::VPolygon{N}) where {N<:Real}
    return P
end

"""
    tohrep(P::VPolygon{N}, ::Type{HPOLYGON}=HPolygon
          ) where {N<:Real, HPOLYGON<:AbstractHPolygon}

Build a constraint representation of the given polygon.

### Input

- `P`        -- polygon in vertex representation
- `HPOLYGON` -- (optional, default: `HPolygon`) type of target polygon

### Output

The same polygon but in constraint representation, an `AbstractHPolygon`.

### Algorithm

The algorithms consists of adding an edge for each consecutive pair of vertices.
Since the vertices are already ordered in counter-clockwise fashion (CWW), the
constraints will be sorted automatically (CCW) if we start with the first edge
between the first and second vertex.
"""
function tohrep(P::VPolygon{N}, ::Type{HPOLYGON}=HPolygon
               ) where {N<:Real, HPOLYGON<:AbstractHPolygon}
    vl = vertices_list(P)
    n = length(vl)
    if n == 0
        # no vertex -> empy set
        return EmptySet{N}(2)
    elseif n == 1
        # only one vertex -> use function for singletons
        return convert(HPOLYGON, Singleton(vl[1]))
    elseif n == 2
        # only two vertices -> use function for line segments
        return convert(HPOLYGON, LineSegment(vl[1], vl[2]))
    else
        # find right-most vertex
        i = div(n, 2)
        x = vl[i][1]
        while i > 1 && vl[i-1][1] > x
            # search forward in list
            i = i - 1
            x = vl[i][1]
        end
        while i < n && vl[i+1][1] > x
            # search backward in list
            i = i + 1
            x = vl[i][1]
        end

        # create constraints ordered in CCW starting at the right-most index
        upper_hull = [halfspace_left(vl[j], vl[j+1]) for j in i:length(vl)-1]
        mid_hull = [halfspace_left(vl[end], vl[1])]
        lower_hull = [halfspace_left(vl[j], vl[j+1]) for j in 1:i-1]
        constraints_list = vcat(upper_hull, mid_hull, lower_hull)
    end
    return HPOLYGON(constraints_list)
end

# --- AbstractPolytope interface functions ---

"""
    vertices_list(P::VPolygon{N}) where {N<:Real}

Return the list of vertices of a convex polygon in vertex representation.

### Input

- `P` -- a polygon vertex representation

### Output

List of vertices.
"""
function vertices_list(P::VPolygon{N}) where {N<:Real}
    return P.vertices
end

# --- LazySet interface functions ---

"""
    σ(d::AbstractVector{N}, P::VPolygon{N}) where {N<:Real}

Return the support vector of a polygon in a given direction.

### Input

- `d` -- direction
- `P` -- polygon in vertex representation

### Output

The support vector in the given direction.
If the direction has norm zero, the first vertex is returned.

### Algorithm

This implementation uses a binary search algorithm when the polygon has more
than ten vertices and a brute-force search when it has ten or less.
For the brute-force search, it compares the projection of
each vector along the given direction and runs in ``O(n)`` where
``n`` is the number of vertices.
For the binary search the algorithm runs in ``O(log n)``.
We follow [this implementation](http://geomalgorithms.com/a14-_extreme_pts.html#polyMax_2D())
based on an algorithm described in [1].

[1] Joseph O'Rourke, Computational Geometry in C (2nd Edition)
"""
function σ(d::AbstractVector{N}, P::VPolygon{N}) where {N<:Real}
    @assert !isempty(P.vertices) "the polygon has no vertices"
    return P.vertices[_σ_helper(d, P)]
end

# heuristic to define the method used to compute the support vector of a polygon
# in vertex representation; if the number of vertices of the polygon is smaller
# than this value, the brute force method is used; otherwise binary search is used
const binary_or_brute_force = 10

# return the index of a support vector of P along d from the list of vertices of P
function _σ_helper(d::AbstractVector{N}, P::VPolygon{N}) where {N <: Real}
    if length(P.vertices) > binary_or_brute_force
        return _binary_support_vector(d, P)
    else
        return _brute_force_support_vector(d, P)
    end
end

function _brute_force_support_vector(d::AbstractVector{N}, P::VPolygon{N}) where {N <: Real}
    i_max = 1
    @inbounds for i in 2:length(P.vertices)
        if dot(d, P.vertices[i] - P.vertices[i_max]) > zero(N)
            i_max = i
        end
    end
    return i_max
end

function _binary_support_vector(d::AbstractVector{N}, P::VPolygon{N}) where {N <: Real}
    m = length(P.vertices)
    @assert m > 2 "the number of vertices in the binary support vector approach should " *
                  "be at least three, but it is $m"
    push!(P.vertices, P.vertices[1]) # add extra vertice on the end equal to the first
    a = 1; b = m + 1 # start chain = [1,n+1] with P.vertices[n+1]=P.vertices[1]
    A = P.vertices[2] - P.vertices[1]
    upA = _up(d, A)
    # test if P.vertices[0] is a local maximum
    if (!upA && !_above(d, P.vertices[m], P.vertices[1])) # P.vertices[1] is the maximum
        pop!(P.vertices) # remove the extra point added
        return 1
    end
    while true
        c = round(Int, (a + b) / 2) # midpoint of [a,b], and 1<c<n+1
        C = P.vertices[c + 1] - P.vertices[c]
        upC = _up(d, C)
        if (!upC && !_above(d, P.vertices[c - 1], P.vertices[c])) # P.vertices[c] is a local maximum
            pop!(P.vertices) # remove the extra point added
            return c # thus it is the maximum
        end
        # no max yet, so continue with the binary search
        # pick one of the two subchains [a,c] or [c,b]
        if (upA && upC && !_above(d, P.vertices[a], P.vertices[c])) ||
        (!upA && (upC || (!upC && _above(d, P.vertices[a], P.vertices[c]))))
            a = c
            A = C
            upA = upC
        else
            b = c
        end
        if (b <= a + 1) # the chain is impossibly small
            pop!(P.vertices) # remove the extra point added
            throw(ErrorException("something went wrong")) # return an error
        end
    end
end

"""
    an_element(P::VPolygon{N}) where {N<:Real}

Return some element of a polygon in vertex representation.

### Input

- `P` -- polygon in vertex representation

### Output

The first vertex of the polygon in vertex representation.
"""
function an_element(P::VPolygon{N}) where {N<:Real}
    @assert !isempty(P.vertices) "the polygon has no vertices"
    return P.vertices[1]
end

"""
    ∈(x::AbstractVector{N}, P::VPolygon{N}) where {N<:Real}

Check whether a given point is contained in a polygon in vertex representation.

### Input

- `x` -- point/vector
- `P` -- polygon in vertex representation

### Output

`true` iff ``x ∈ P``.

### Algorithm

This implementation exploits that the polygon's vertices are sorted in
counter-clockwise fashion.
Under this assumption we can just check if the vertex lies on the left of each
edge, using the dot product.

### Examples

```jldoctest
julia> P = VPolygon([[2.0, 3.0], [3.0, 1.0], [5.0, 1.0], [4.0, 5.0]];
                    apply_convex_hull=false);

julia> [4.5, 3.1] ∈ P
false
julia> [4.5, 3.0] ∈ P
true
julia> [4.4, 3.4] ∈ P  #  point lies on the edge
true
```
"""
function ∈(x::AbstractVector{N}, P::VPolygon{N}) where {N<:Real}
    @assert length(x) == 2

    # special cases: 0 or 1 vertex
    if length(P.vertices) == 0
        return false
    elseif length(P.vertices) == 1
        return x == P.vertices[1]
    end

    if !is_right_turn(P.vertices[1], x, P.vertices[end])
        return false
    end
    for i in 2:length(P.vertices)
        if !is_right_turn(P.vertices[i], x, P.vertices[i-1])
            return false
        end
    end
    return true
end

"""
    _random_zero_sum_vector(rng::AbstractRNG, N::Type{<:Real}, n::Int)

Create a random vector with entries whose sum is zero.

### Input

- `rng` -- random number generator
- `N`   -- numeric type
- `n`   -- length of vector

### Output

A random vector of random numbers such that all positive entries come first and
all negative entries come last, and such that the total sum is zero.

### Algorithm

This is a preprocessing step of the algorithm
[here](https://stackoverflow.com/a/47358689) based on
[P. Valtr. Probability that n random points are in convex
position](https://link.springer.com/article/10.1007%2FBF01271274).
"""
function _random_zero_sum_vector(rng::AbstractRNG, N::Type{<:Real}, n::Int)
    # generate a sorted list of random x and y coordinates
    list = sort!(randn(rng, N, n))
    while (length(remove_duplicates_sorted!(list)) < n)
        # make sure that no duplicates exist
        list = sort!(append!(list, randn(rng, N, length(list) - n)))
    end
    # lists of consecutive points
    in_l1 = rand(rng, Bool, n-2)
    l1 = Vector{N}() # normal
    l2 = Vector{N}() # inverted
    push!(l1, list[1])
    push!(l2, list[1])
    for i in 2:n-1
        push!(in_l1[i-1] ? l1 : l2, list[i])
    end
    push!(l1, list[end])
    push!(l2, list[end])
    # convert to vectors representing the distance (order does not matter)
    dist = Vector{N}()
    sizehint!(dist, n)
    for i in 1:length(l1)-1
        push!(dist, l1[i+1] - l1[i])
    end
    for i in 1:length(l2)-1
        push!(dist, l2[i] - l2[i+1])
    end
    @assert isapprox(sum(dist), zero(N), atol=1e-6)
    return dist
end

"""
    rand(::Type{VPolygon}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)

Create a random polygon in vertex representation.

### Input

- `VPolygon`     -- type for dispatch
- `N`            -- (optional, default: `Float64`) numeric type
- `dim`          -- (optional, default: 2) dimension
- `rng`          -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`         -- (optional, default: `nothing`) seed for reseeding
- `num_vertices` -- (optional, default: `-1`) number of vertices of the
                    polygon (see comment below)

### Output

A random polygon in vertex representation.

### Algorithm

We follow the idea [here](https://stackoverflow.com/a/47358689) based on
[P. Valtr. Probability that n random points are in convex
position](https://link.springer.com/article/10.1007%2FBF01271274).
There is also a nice video available
[here](http://cglab.ca/~sander/misc/ConvexGeneration/convex.html).

The number of vertices can be controlled with the argument `num_vertices`.
For a negative value we choose a random number in the range `3:10`.
"""
function rand(::Type{VPolygon};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing,
              num_vertices::Int=-1)
    @assert dim == 2 "cannot create a random VPolygon of dimension $dim"
    rng = reseed(rng, seed)
    if num_vertices < 0
        num_vertices = rand(3:10)
    end

    # special cases, 0 or 1 vertex
    if num_vertices == 0
        return VPolygon{N}()
    elseif num_vertices == 1
        return VPolygon([randn(rng, N, 2)])
    end

    # general case, >= 2 vertices

    # get random horizontal and vertical vectors
    horiz = _random_zero_sum_vector(rng, N, num_vertices)
    vert = _random_zero_sum_vector(rng, N, num_vertices)

    # randomly combine horizontal and vertical vectors
    m = num_vertices
    directions = Vector{Vector{N}}(undef, num_vertices)
    shuffle(rng, vert)
    for (i, x) in enumerate(horiz)
        y = splice!(vert, rand(rng, 1:m))
        directions[i] = [x, y]
        m -= 1
    end
    sort!(directions, lt=<=) # sort by angle

    # connect directions
    vertices = Vector{Vector{N}}(undef, num_vertices)
    # random starting point
    vertices[1] = randn(rng, N, 2)
    for i in 1:length(directions)-1
        vertices[i+1] = vertices[i] + directions[i]
    end
    @assert isapprox(vertices[end] + directions[end], vertices[1], atol=1e-6)
    return VPolygon(vertices; apply_convex_hull=true)
end

"""
    constraints_list(P::VPolygon{N}) where {N<:Real}

Return the list of constraints defining a polygon in V-representation.

### Input

- `P` -- polygon in V-representation

### Output

The list of constraints of the polygon.

### Algorithm

First the H-representation of ``P`` is computed, then its list of constraints
is returned.
"""
function constraints_list(P::VPolygon{N}) where {N<:Real}
    return constraints_list(tohrep(P))
end

"""
    convex_hull(P::VPolygon{N}, Q::VPolygon{N};
                [algorithm]::String="monotone_chain") where {N<:Real}

Return the convex hull of two polygons in vertex representation.

### Input

- `P`         -- polygon in vertex representation
- `Q`         -- another polygon in vertex representation
- `algorithm` -- (optional, default: "monotone_chain") the algorithm used to
                 compute the convex hull

### Output

A new polygon such that its vertices are the convex hull of the given two polygons.

### Algorithm

A convex hull algorithm is used to compute the convex hull of the vertices of the
given input polygons `P` and `Q`; see `?convex_hull` for details on the available
algorithms. The vertices of the output polygon are sorted in counter-clockwise
fashion.
"""
function convex_hull(P::VPolygon{N}, Q::VPolygon{N};
                     algorithm::String="monotone_chain") where {N<:Real}
    vunion = [P.vertices; Q.vertices]
    convex_hull!(vunion; algorithm=algorithm)
    return VPolygon(vunion, apply_convex_hull=false)
end

"""
    translate(P::VPolygon{N}, v::AbstractVector{N}) where {N<:Real}

Translate (i.e., shift) a polygon in vertex representation by a given vector.

### Input

- `P` -- polygon in vertex representation
- `v` -- translation vector

### Output

A translated polygon in vertex representation.

### Algorithm

We add the vector to each vertex of the polygon.
"""
function translate(P::VPolygon{N}, v::AbstractVector{N}) where {N<:Real}
    @assert length(v) == dim(P) "cannot translate a $(dim(P))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    return VPolygon([x + v for x in vertices_list(P)])
end

"""
    minkowski_sum(P::VPolygon{N}, Q::VPolygon{N}) where {N<:Real}

The Minkowski Sum of two polygon in vertex representation.

### Input

- `P` -- polygon in vertex representation
- `Q` -- another polygon in vertex representation

### Output

A polygon in vertex representation.

### Algorithm

We treat each edge of the polygons as a vector, attaching them in polar order
(attaching the tail of the next vector to the head of the previous vector). The
resulting polygonal chain will be a polygon, which is the Minkowski sum of the
given polygons. This algorithm assumes that the vertices of P and Q are sorted
in counter-clockwise fashion and has linear complexity O(m+n) where m and n are
the number of vertices of P and Q respectively.

"""
function minkowski_sum(P::VPolygon{N}, Q::VPolygon{N}) where {N<:Real}
    vlistP = vertices_list(P)
    vlistQ = vertices_list(Q)
    mP = length(vlistP)
    mQ = length(vlistQ)

    EAST = N[1, 0]
    ORIGIN = N[0, 0]
    k = _σ_helper(EAST, P)
    j = _σ_helper(EAST, Q)
    R = Vector{Vector{N}}(undef, mP+mQ)
    fill!(R, ORIGIN)

    i = 1
    while i <= size(R, 1)
        P₁, P₂ = vlistP[(k-1)%mP + 1], vlistP[(k%mP + 1)]
        P₁P₂ = P₂ - P₁
        Q₁, Q₂ = vlistQ[(j-1)%mQ + 1], vlistQ[(j%mQ + 1)]
        Q₁Q₂ = Q₂ - Q₁
        R[i] = P₁ + Q₁
        turn = right_turn(P₁P₂, Q₁Q₂, ORIGIN)
        if turn > 0
            k += 1
        elseif turn < 0
            j += 1
        else
            pop!(R)
            k += 1
            j += 1
        end
        i += 1
    end
    return VPolygon(R)
end


text:
using Arrow
using Compat, Compat.Test, Compat.Random, Compat.Dates

if VERSION < v"0.7.0-"
    using Missings
end

const ≅ = isequal

const SEED = 999

const N_OUTER = 4
const N_IDX_CHECK = 32
const MAX_IDX_LEN = 32
const MAX_VECTOR_LENGTH = 256
const MAX_STRING_LENGTH = 32

const PRIMITIVE_ELTYPES = [Float32, Float64, Int32, Int64, UInt16]

srand(SEED)

include("testutils.jl")

@testset "Primitive_access" begin
    for i ∈ 1:N_OUTER
        len, dtype, lpad, b, v = rand_primitive_buffer()
        p = Primitive{dtype}(b, lpad+1, len)
        # integer indices
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test p[k] == v[k]
        end
        # AbstractVector{<:Integer} indices
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test p[idx] == v[idx]
        end
        @test p[:] == v
    end
end

@testset "Primitive_construct" begin
    for i ∈ 1:N_OUTER
        dtype = rand(PRIMITIVE_ELTYPES)
        len = rand(1:MAX_VECTOR_LENGTH)
        v = rand(dtype, len)
        p = Primitive(v)
        # integer indices
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test p[k] == v[k]
        end
        # AbstractVector{<:Integer} indices
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test p[idx] == v[idx]
        end
        @test p[:] == v
    end
end

@testset "NullablePrimitive_access" begin
    for i ∈ 1:N_OUTER
        len, dtype, bmask, lpad, b, v = rand_nullableprimitive_buffer()
        p = NullablePrimitive{dtype}(b, 1+lpad, 1+bmask+lpad, len)
        # integer indices
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test (ismissing(p[k]) && ismissing(v[k])) || (p[k] == v[k])
        end
        # AbstractVector{<:Integer} indices
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            sp = p[idx]
            sv = v[idx]
            ssp = collect(skipmissing(sp))
            ssv = collect(skipmissing(sv))
            @test length(sp) == length(sv) && ssp == ssv
        end
        @test p[:] ≅ v
    end
end

@testset "NullablePrimitive_construct" begin
    for i ∈ 1:N_OUTER
        dtype = rand(PRIMITIVE_ELTYPES)
        len = rand(1:MAX_VECTOR_LENGTH)
        v = Union{dtype,Missing}[rand(Bool) ? missing : rand(dtype) for i ∈ 1:len]
        p = NullablePrimitive(v)
        # integer indices
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test p[k] ≅ v[k]
        end
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test p[idx] ≅ v[idx]
        end
        @test p[:] ≅ v
    end
end

@testset "List_access" begin
    len  = 5
    offs = Int32[0,4,7,8,12,14]
    vals = convert(Vector{UInt8}, codeunits("firewalkwithme"))
    valspad = zeros(UInt8, Arrow.paddinglength(length(vals)))
    lpad = randpad()
    rpad = randpad()
    b = vcat(lpad, convert(Vector{UInt8}, reinterpret(UInt8, offs)), vals, valspad, rpad)
    l = List{String}(b, 1+length(lpad), 1+length(lpad)+sizeof(Int32)*length(offs),
                     len, UInt8, length(vals))
    @test offsets(l)[:] == offs
    @test values(l)[:] == vals
    @test l[1] == "fire"
    @test l[2] == "wal"
    @test l[3] == "k"
    @test l[4] == "with"
    @test l[5] == "me"
    @test l[[1,3,5]] == ["fire", "k", "me"]
    @test l[[false,true,false,true,false]] == ["wal", "with"]
    @test l[:] == ["fire", "wal", "k", "with", "me"]
end

@testset "List_construct" begin
    for i ∈ 1:N_OUTER
        len = rand(1:MAX_VECTOR_LENGTH)
        v = String[randstring(rand(0:MAX_STRING_LENGTH)) for i ∈ 1:len]
        l = List(v)
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test l[k] == v[k]
        end
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test l[idx] == v[idx]
        end
        @test l[:] == v
    end
end

@testset "NullableList_access" begin
    len = 7
    offs = Int32[0,4,9,9,14,14,17,21]
    vals = convert(Vector{UInt8}, codeunits("kirkspockbonesncc1701"))
    valspad = zeros(UInt8, Arrow.paddinglength(length(vals)))
    pres = Bool[true,true,false,true,false,true,true]
    mask = Arrow.bitpackpadded(pres)
    lpad = randpad()
    rpad = randpad()
    b = vcat(lpad, mask, convert(Vector{UInt8}, reinterpret(UInt8, offs)), vals, valspad, rpad)
    l = NullableList{String}(b, 1+length(lpad), 1+length(lpad)+length(mask),
                             1+length(lpad)+length(mask)+sizeof(Int32)*length(offs),
                             len, UInt8, length(vals))
    @test offsets(l)[:] == offs
    @test values(l)[:] == vals
    @test l[1] == "kirk"
    @test l[2] == "spock"
    @test ismissing(l[3])
    @test l[4] == "bones"
    @test ismissing(l[5])
    @test l[6] == "ncc"
    @test l[7] == "1701"
    @test l[[2,5,4]] ≅ ["spock", missing, "bones"]
    @test l[[true,false,true,false,false,false,false]] ≅ ["kirk", missing]
    @test l[:] ≅ ["kirk", "spock", missing, "bones", missing, "ncc", "1701"]
end

@testset "NullableList_construct" begin
    for i ∈ 1:N_OUTER
        len = rand(1:MAX_VECTOR_LENGTH)
        v = Union{String,Missing}[rand(Bool) ? missing : randstring(rand(0:MAX_STRING_LENGTH))
                                  for i ∈ 1:len]
        l = NullableList(v)
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test l[k] ≅ v[k]
        end
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test l[idx] ≅ v[idx]
        end
        @test l[:] ≅ v
    end
end

@testset "BitPrimitive_access" begin
    len = 10
    bits = vcat(UInt8[0xfd,0x02], zeros(UInt8, 6))
    lpad = randpad()
    rpad = randpad()
    b = vcat(lpad, bits, rpad)
    v = Bool[true,false,true,true,true,true,true,true,false,true]
    p = BitPrimitive(b, 1+length(lpad), 10)
    for j ∈ 1:N_IDX_CHECK
        k = rand(1:len)
        @test p[k] == v[k]
    end
    for j ∈ 1:N_IDX_CHECK
        idx = rand(1:len, rand(1:MAX_IDX_LEN))
        @test p[idx] == v[idx]
    end
    @test p[:] == v
end

@testset "BitPrimitive_construct" begin
    for i ∈ 1:N_OUTER
        len = rand(1:MAX_VECTOR_LENGTH)
        v = rand(Bool, len)
        p = BitPrimitive(v)
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test p[k] == v[k]
        end
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test p[idx] == v[idx]
        end
        @test p[:] == v
    end
end

@testset "NullableBitPrimitive_access" begin
    len = 7
    bits = vcat(UInt8[0x4f], zeros(UInt8, 7))
    bmask = vcat(UInt8[0x79], zeros(UInt8, 7))
    lpad = randpad()
    rpad = randpad()
    b = vcat(lpad, bmask, bits, rpad)
    v = [true,missing,missing,true,false,false,true]
    p = NullableBitPrimitive(b, 1+length(lpad), 1+length(lpad)+length(bmask), length(v))
    for j ∈ 1:N_IDX_CHECK
        k = rand(1:len)
        @test p[k] ≅ v[k]
    end
    for j ∈ 1:N_IDX_CHECK
        idx = rand(1:len, rand(1:MAX_IDX_LEN))
        @test p[idx] ≅ v[idx]
    end
    @test p[:] ≅ v
end

@testset "NullableBitPrimitive_construct" begin
    for i ∈ 1:N_OUTER
        len = rand(1:MAX_VECTOR_LENGTH)
        pres = rand(Bool, len)
        vals = rand(Bool, len)
        v = Union{Bool,Missing}[pres[i] ? missing : vals[i] for i ∈ 1:len]
        p = NullableBitPrimitive(v)
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test p[k] ≅ v[k]
        end
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test p[idx] ≅ v[idx]
        end
        @test p[:] ≅ v
    end
end

@testset "arrowformat_construct" begin
    len = rand(1:MAX_VECTOR_LENGTH)
    randstring_() = randstring(rand(1:MAX_STRING_LENGTH))

    p = arrowformat(rand(Float64, len))
    @test typeof(p) == Primitive{Float64}

    mask = rand(Bool, len)

    v = Union{Float32,Missing}[mask[i] ? rand(Float32) : missing for i ∈ 1:len]
    p = arrowformat(v)
    @test typeof(p) == NullablePrimitive{Float32}

    v = String[randstring_() for i ∈ 1:len]
    p = arrowformat(v)
    @test typeof(p) == List{String,Primitive{UInt8}}

    v = Union{String,Missing}[mask[i] ? randstring_() : missing for i ∈ 1:len]
    p = arrowformat(v)
    @test typeof(p) == NullableList{String,Primitive{UInt8}}

    v = rand(Bool, len)
    p = arrowformat(v)
    @test typeof(p) == BitPrimitive

    v = Union{Bool,Missing}[mask[i] ? rand(Bool) : missing for i ∈ 1:len]
    p = arrowformat(v)
    @test typeof(p) == NullableBitPrimitive

    v = Date[Date(1), Date(2)]
    p = arrowformat(v)
    @test typeof(p) == Primitive{Arrow.Datestamp}

    v = DateTime[DateTime(0), DateTime(1)]
    p = arrowformat(v)
    @test typeof(p) == Primitive{Arrow.Timestamp{Millisecond}}

    v = Time[Time(0), Time(1)]
    p = arrowformat(v)
    @test typeof(p) == Primitive{Arrow.TimeOfDay{Nanosecond,Int64}}

    v = Union{Date,Missing}[Date(1), missing]
    p = arrowformat(v)
    @test typeof(p) == NullablePrimitive{Arrow.Datestamp}

    v = Union{DateTime,Missing}[DateTime(0), missing]
    p = arrowformat(v)
    @test typeof(p) == NullablePrimitive{Arrow.Timestamp{Millisecond}}

    v = Union{Time,Missing}[Time(0), missing]
    p = arrowformat(v)
    @test typeof(p) == NullablePrimitive{Arrow.TimeOfDay{Nanosecond,Int64}}
end


text:
function apply_dirichlet_bc!(
    matrix::SparseMatrixCSC,
    rhs,
    index::Z,
    value::R,
) where {Z<:Integer,R<:Real}

    modifyrhs = matrix[:, index]
    for i in modifyrhs.nzind
        if i == index
            rhs[i] = modifyrhs[i] * value
        else
            rhs[i] -= modifyrhs[i] * value
            matrix[i, index] = 0.0
            matrix[index, i] = 0.0
        end
    end
end

function apply_dirichlet_bc!(
    matrix::Matrix,
    rhs,
    index::Z,
    value::R,
) where {Z<:Integer,R<:Real}

    m, n = size(matrix)
    @assert length(rhs) == m
    @assert m == n

    modifyrhs = matrix[:,index]
    for i = 1:n
        if i == index
            rhs[i] = matrix[i,i] * value
        else
            matrix[index, i] = 0.0
            matrix[i,index] = 0.0
            rhs[i] -= modifyrhs[i]*value
        end
    end
end

function apply_dirichlet_bc!(
    matrix,
    rhs,
    nodeids::V,
    dofs,
    values,
    dofspernode,
) where {V<:AbstractVector}

    @assert length(nodeids) == length(dofs) == length(values)
    for (nodeid, dof, val) in zip(nodeids, dofs, values)
        index = node_to_dof_id(nodeid, dof,dofspernode)
        apply_dirichlet_bc!(matrix, rhs, index, val)
    end
end

function apply_dirichlet_bc!(
    matrix,
    rhs,
    nodeids::V,
    dofs::Z,
    value::R,
    dofspernode,
) where {V<:AbstractVector,Z<:Integer,R<:Real}
    numnodes = length(nodeids)
    vecdofs = repeat([dofs],numnodes)
    vecvalues = repeat([value],numnodes)
    apply_dirichlet_bc!(matrix,rhs,nodeids,vecdofs,vecvalues,dofspernode)
end

function apply_dirichlet_bc!(
    matrix,
    rhs,
    nodeids::V,
    bcvals::M,
) where {V<:AbstractVector,M<:AbstractMatrix}

    dofspernode, numnodes = size(bcvals)
    @assert length(nodeids) == numnodes
    for (idx, nodeid) in enumerate(nodeids)
        for dof = 1:dofspernode
            index = node_to_dof_id(nodeid, dof, dofspernode)
            apply_dirichlet_bc!(matrix, rhs, index, bcvals[dof, idx])
        end
    end
end

function apply_dirichlet_bc!(
    matrix,
    rhs,
    nodeid::Z,
    bcvals::V,
) where {Z<:Integer,V<:AbstractVector}

    dofs = length(bcvals)
    for dof = 1:dofs
        index = node_to_dof_id(nodeid, dof, dofs)
        apply_dirichlet_bc!(matrix, rhs, index, bcvals[dof])
    end
end


text:
using DataFrames
using CSV
using Plots
using BenchmarkProfiles

function genpprof(args)
    if length(args) == 1
        metric = args[1]
    else
        metric = "FEVAL"
    end
    base = DataFrame(CSV.File("base_bench.csv"))
    new = DataFrame(CSV.File("bench_results.csv"))

    pp = performance_profile(
    PlotsBackend(),
    Float64.(hcat(base[!, metric], new[!, metric])),
    ["Base", "New"],
    title = "$metric profile"
    )
    savefig("profile.png")
end

if abspath(PROGRAM_FILE) == @__FILE__
    genpprof(ARGS)
end


text:
struct MinimaxAlphaBeta{G <: Game, H} <: Strategy
    depth::Int
    heuristic::H
end

MinimaxAlphaBeta{G}(depth::Int, heuristic::H) where {G<:Game, H} = MinimaxAlphaBeta{G,H}(depth, heuristic)

function BoardGames.getmove(board, s::MinimaxAlphaBeta; verbose=false)
    @assert !isempty(getmoves(board))
    move, score = expand(board, s, playerturn(board), s.depth, -Inf, Inf)
    verbose && println(score)

    return getmoves(board)[move]
end

function expand(board, s::MinimaxAlphaBeta, player::Int, depth::Int, α::Float64, β::Float64)
    @assert depth >= 0
    moves = getmoves(board)

    maximize = player == playerturn(board)
    best_move = 0
    if maximize
        best_score = -Inf
    else
        best_score = Inf
    end

    if isempty(moves) || depth == 0
        best_move = 0
        best_score = s.heuristic(board, player)
    else
        for i in 1:length(moves)
            next_board = play(board, moves[i])

            move, score = expand(next_board, s, player, depth-1, α, β)
            if (score > best_score && maximize) ||
               (score < best_score && maximize)
                best_move = i
                best_score = score
                if maximize
                    α = max(α, best_score)
                    if α >= β
                        break
                    end
                else
                    β = min(β, best_score)
                    if β <= α
                        break
                    end
                end
            end
        end
    end
    return (best_move, best_score)
end


text:
@doc raw"""
    differential_bezier_control(M::Manifold, b::BezierSegment, t::Float, X::BezierSegment)
    differential_bezier_control!(
        M::Manifold,
        Y,
        b::BezierSegment,
        t,
        X::BezierSegment
    )

evaluate the differential of the Bézier curve with respect to its control points
`b` and tangent vectors `X` given in the tangent spaces of the control points. The result
is the “change” of the curve at `t```∈[0,1]``. The comptation can be done in place of `Y`.

See [`de_casteljau`](@ref) for more details on the curve.
"""
function differential_bezier_control(M::Manifold, b::BezierSegment, t, X::BezierSegment)
    # iterative, because recursively would be too many Casteljau evals
    Y = similar(first(X.pts))
    return differential_bezier_control!(M, Y, b, t, X)
end
function differential_bezier_control!(M::Manifold, Y, b::BezierSegment, t, X::BezierSegment)
    # iterative, because recursively would be too many Casteljau evals
    Z = similar(X.pts)
    c = deepcopy(b.pts)
    for l in length(c):-1:2
        Z[1:(l - 1)] .=
            differential_geodesic_startpoint.(
                Ref(M), c[1:(l - 1)], c[2:l], Ref(t), X.pts[1:(l - 1)]
            ) .+
            differential_geodesic_endpoint.(
                Ref(M), c[1:(l - 1)], c[2:l], Ref(t), X.pts[2:l]
            )
        c[1:(l - 1)] = shortest_geodesic.(Ref(M), c[1:(l - 1)], c[2:l], Ref(t))
    end
    return copyto!(M, Y, Z[1])
end
@doc raw"""
    differential_bezier_control(
        M::Manifold,
        b::BezierSegment,
        T::AbstractVector,
        X::BezierSegment,
    )
    differential_bezier_control!(
        M::Manifold,
        Y,
        b::BezierSegment,
        T::AbstractVector,
        X::BezierSegment,
    )

evaluate the differential of the Bézier curve with respect to its control points
`b` and tangent vectors `X` in the tangent spaces of the control points. The result
is the “change” of the curve at the points `T`, elementwise in ``t∈[0,1]``.
The compuation can be done in place of `Y`.

See [`de_casteljau`](@ref) for more details on the curve.
"""
function differential_bezier_control(
    M::Manifold, b::BezierSegment, T::AbstractVector, X::BezierSegment
)
    return differential_bezier_control.(Ref(M), Ref(b), T, Ref(X))
end
function differential_bezier_control!(
    M::Manifold, Y, b::BezierSegment, T::AbstractVector, X::BezierSegment
)
    return differential_bezier_control!.(Ref(M), Y, Ref(b), T, Ref(X))
end
@doc raw"""
    differential_bezier_control(
        M::Manifold,
        B::AbstractVector{<:BezierSegment},
        t,
        X::AbstractVector{<:BezierSegment}
    )
    differential_bezier_control!(
        M::Manifold,
        Y::AbstractVector{<:BezierSegment}
        B::AbstractVector{<:BezierSegment},
        t,
        X::AbstractVector{<:BezierSegment}
    )

evaluate the differential of the composite Bézier curve with respect to its
control points `B` and tangent vectors `Ξ` in the tangent spaces of the control
points. The result is the “change” of the curve at `t```∈[0,N]``, which depends
only on the corresponding segment. Here, ``N`` is the length of `B`.
The compuation can be done in place of `Y`.

See [`de_casteljau`](@ref) for more details on the curve.
"""
function differential_bezier_control(
    M::Manifold, B::AbstractVector{<:BezierSegment}, t, X::AbstractVector{<:BezierSegment}
)
    if (0 > t) || (t > length(B))
        return throw(
            DomainError(
                t,
                "The parameter $(t) to evaluate the composite Bézier curve at is outside the interval [0,$(length(B))].",
            ),
        )
    end
    seg = max(ceil(Int, t), 1)
    localT = ceil(Int, t) == 0 ? 0.0 : t - seg + 1
    Y = differential_bezier_control(M, B[seg], localT, X[seg])
    if (Integer(t) == seg) && (seg < length(B)) # boundary case, -> seg-1 is also affecting the boundary differential
        Y .+= differential_bezier_control(M, B[seg + 1], localT - 1, X[seg + 1])
    end
    return Y
end
function differential_bezier_control!(
    M::Manifold,
    Y,
    B::AbstractVector{<:BezierSegment},
    t,
    X::AbstractVector{<:BezierSegment},
)
    if (0 > t) || (t > length(B))
        return throw(
            DomainError(
                t,
                "The parameter $(t) to evaluate the composite Bézier curve at is outside the interval [0,$(length(B))].",
            ),
        )
    end
    seg = max(ceil(Int, t), 1)
    localT = ceil(Int, t) == 0 ? 0.0 : t - seg + 1
    differential_bezier_control!(M, Y, B[seg], localT, X[seg])
    if (Integer(t) == seg) && (seg < length(B)) # boundary case, -> seg-1 is also affecting the boundary differential
        Y .+= differential_bezier_control(M, B[seg + 1], localT - 1, X[seg + 1])
    end
    return Y
end

@doc raw"""
    differential_bezier_control(
        M::Manifold,
        B::AbstractVector{<:BezierSegment},
        T::AbstractVector
        Ξ::AbstractVector{<:BezierSegment}
    )
    differential_bezier_control!(
        M::Manifold,
        Θ::AbstractVector{<:BezierSegment}
        B::AbstractVector{<:BezierSegment},
        T::AbstractVector
        Ξ::AbstractVector{<:BezierSegment}
    )

evaluate the differential of the composite Bézier curve with respect to its
control points `B` and tangent vectors `Ξ` in the tangent spaces of the control
points. The result is the “change” of the curve at the points in `T`, which are elementwise
in ``[0,N]``, and each depending the corresponding segment(s). Here, ``N`` is the
length of `B`. For the mutating variant the result is computed in `Θ`.

See [`de_casteljau`](@ref) for more details on the curve and [^BergmannGousenbourger2018].

[^BergmannGousenbourger2018]:
    > Bergmann, R. and Gousenbourger, P.-Y.: _A variational model for data fitting on manifolds
    > by minimizing the acceleration of a Bézier curve_.
    > Frontiers in Applied Mathematics and Statistics, 2018.
    > doi: [10.3389/fams.2018.00059](https://dx.doi.org/10.3389/fams.2018.00059),
    > arXiv: [1807.10090](https://arxiv.org/abs/1807.10090)
"""
function differential_bezier_control(
    M::Manifold,
    B::AbstractVector{<:BezierSegment},
    T::AbstractVector,
    Ξ::AbstractVector{<:BezierSegment},
)
    return differential_bezier_control.(Ref(M), Ref(B), T, Ref(Ξ))
end
function differential_bezier_control!(
    M::Manifold,
    Y,
    B::AbstractVector{<:BezierSegment},
    T::AbstractVector,
    Ξ::AbstractVector{<:BezierSegment},
)
    return differential_bezier_control!.(Ref(M), Y, Ref(B), T, Ref(Ξ))
end

@doc raw"""
    differential_geodesic_startpoint(M, p, q, t, X)
    differential_geodesic_startpoint!(M, Y, p, q, t, X)

computes ``D_p g(t;p,q)[η]`` (in place of `Y`).

# See also
 [`differential_geodesic_endpoint`](@ref), [`jacobi_field`](@ref)
"""
function differential_geodesic_startpoint(M::Manifold, p, q, t, X)
    return jacobi_field(M, p, q, t, X, βdifferential_geodesic_startpoint)
end
function differential_geodesic_startpoint!(M::Manifold, Y, p, q, t, X)
    return jacobi_field!(M, Y, p, q, t, X, βdifferential_geodesic_startpoint)
end

@doc raw"""
    differential_geodesic_endpoint(M, p, q, t, X)
    differential_geodesic_endpoint!(M, Y, p, q, t, X)

computes ``D_qg(t;p,q)[X]`` (in place of `Y`).

# See also
 [`differential_geodesic_startpoint`](@ref), [`jacobi_field`](@ref)
"""
function differential_geodesic_endpoint(M::Manifold, p, q, t, X)
    return jacobi_field(M, q, p, 1 - t, X, βdifferential_geodesic_startpoint)
end
function differential_geodesic_endpoint!(M::Manifold, Y, p, q, t, X)
    return jacobi_field!(M, Y, q, p, 1 - t, X, βdifferential_geodesic_startpoint)
end

@doc raw"""
    differential_exp_basepoint(M, p, X, Y)
    differential_exp_basepoint!(M, Z, p, X, Y)

Compute ``D_p\exp_p X[Y]`` (in place of `Z`).

# See also
[`differential_exp_argument`](@ref), [`jacobi_field`](@ref)
"""
function differential_exp_basepoint(M::Manifold, p, X, Y)
    return jacobi_field(M, p, exp(M, p, X), 1.0, Y, βdifferential_exp_basepoint)
end
function differential_exp_basepoint!(M::Manifold, Z, p, X, Y)
    return jacobi_field!(M, Z, p, exp(M, p, X), 1.0, Y, βdifferential_exp_basepoint)
end

@doc raw"""
    differential_exp_argument(M, p, X, Y)
    differential_exp_argument!(M, Z, p, X, Y)

computes ``D_X\exp_pX[Y]`` (in place of `Z`).
Note that ``X ∈  T_X(T_p\mathcal M) = T_p\mathcal M`` is still a tangent vector.

# See also
 [`differential_exp_basepoint`](@ref), [`jacobi_field`](@ref)
"""
function differential_exp_argument(M::Manifold, p, X, Y)
    return jacobi_field(M, p, exp(M, p, X), 1.0, Y, βdifferential_exp_argument)
end
function differential_exp_argument!(M::Manifold, Z, p, X, Y)
    return jacobi_field!(M, Z, p, exp(M, p, X), 1.0, Y, βdifferential_exp_argument)
end

@doc raw"""
    differential_log_basepoint(M, p, q, X)
    differential_log_basepoint!(M, Y, p, q, X)

computes ``D_p\log_pq[X]`` (in place of `Y`).

# See also
 [`differential_log_argument`](@ref), [`jacobi_field`](@ref)
"""
function differential_log_basepoint(M::Manifold, p, q, X)
    return jacobi_field(M, p, q, 0.0, X, βdifferential_log_basepoint)
end
function differential_log_basepoint!(M::Manifold, Y, p, q, X)
    return jacobi_field!(M, Y, p, q, 0.0, X, βdifferential_log_basepoint)
end

@doc raw"""
    differential_log_argument(M, p, q, X)
    differential_log_argument(M, Y, p, q, X)

computes ``D_q\log_pq[X]`` (in place of `Y`).

# See also
 [`differential_log_basepoint`](@ref), [`jacobi_field`](@ref)
"""
function differential_log_argument(M::Manifold, p, q, X)
    # order of p and q has to be reversed in this call, cf. Persch, 2018 Lemma 2.3
    return jacobi_field(M, q, p, 1.0, X, βdifferential_log_argument)
end
function differential_log_argument!(M::Manifold, Y, p, q, X)
    # order of p and q has to be reversed in this call, cf. Persch, 2018 Lemma 2.3
    return jacobi_field!(M, Y, q, p, 1.0, X, βdifferential_log_argument)
end

@doc raw"""
    Y = differential_forward_logs(M, p, X)
    differential_forward_logs!(M, Y, p, X)

compute the differenital of [`forward_logs`](@ref) ``F`` on the `PowerManifold` manifold
`M` at `p` and direction `X` , in the power manifold array, the differential of the function

```math
F_i(x) = \sum_{j ∈ \mathcal I_i} \log_{p_i} p_j, \quad i ∈ \mathcal G,
```

where ``\mathcal G`` is the set of indices of the `PowerManifold` manifold `M`
and ``\mathcal I_i`` denotes the forward neighbors of ``i``.

# Input
* `M`     – a `PowerManifold` manifold
* `p`     – a point.
* `X`     – a tangent vector.

# Ouput
* `Y` – resulting tangent vector in ``T_x\mathcal N`` representing the differentials of the
    logs, where ``\mathcal N`` is thw power manifold with the number of dimensions added
    to `size(x)`. The computation can also be done in place.
"""
function differential_forward_logs(M::PowerManifold, p, X)
    power_size = power_dimensions(M)
    R = CartesianIndices(Tuple(power_size))
    d = length(power_size)
    maxInd = last(R).I
    d2 = (d > 1) ? ones(Int, d + 1) + (d - 1) * (1:(d + 1) .== d + 1) : 1
    if d > 1
        N = PowerManifold(M.manifold, NestedPowerRepresentation(), power_size..., d)
    else
        N = PowerManifold(M.manifold, NestedPowerRepresentation(), power_size...)
    end
    Y = zero_tangent_vector(N, repeat(p; inner=d2))
    return differential_forward_logs!(M, Y, p, X)
end
function differential_forward_logs!(M::PowerManifold, Y, p, X)
    power_size = power_dimensions(M)
    R = CartesianIndices(Tuple(power_size))
    d = length(power_size)
    maxInd = last(R).I
    e_k_vals = [1 * (1:d .== k) for k in 1:d]
    if d > 1
        N = PowerManifold(M.manifold, NestedPowerRepresentation(), power_size..., d)
    else
        N = PowerManifold(M.manifold, NestedPowerRepresentation(), power_size...)
    end
    for i in R # iterate over all pixel
        for k in 1:d # for all direction combinations
            I = i.I # array of index
            J = I .+ e_k_vals[k] #i + e_k is j
            if all(J .<= maxInd)
                # this is neighbor in range,
                # collects two, namely in kth direction since xi appears as base and arg
                Y[N, I..., k] =
                    differential_log_basepoint(
                        M.manifold, p[M, I...], p[M, J...], X[M, I...]
                    ) .+ differential_log_argument(
                        M.manifold, p[M, I...], p[M, J...], X[M, J...]
                    )
            else
                Y[N, I..., k] = zero_tangent_vector(M.manifold, p[M, I...])
            end
        end # directions
    end # i in R
    return Y
end


text:
# Autogenerated wrapper script for p3a_jll for aarch64-apple-darwin
export libp3a

using Kokkos_jll
using MPICH_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("p3a")
JLLWrappers.@declare_library_product(libp3a, "@rpath/libp3a.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Kokkos_jll, MPICH_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libp3a,
        "lib/libp3a.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ContinuumArrays, QuasiArrays, LazyArrays, IntervalSets, FillArrays, LinearAlgebra, BandedMatrices, FastTransforms, InfiniteArrays, Test, Base64
import ContinuumArrays: ℵ₁, materialize, AffineQuasiVector, BasisLayout, AdjointBasisLayout, SubBasisLayout, ℵ₁,
                        MappedBasisLayout, AdjointMappedBasisLayout, MappedWeightedBasisLayout, igetindex, TransformFactorization, Weight, WeightedBasisLayout, SubWeightedBasisLayout, WeightLayout,
                        Expansion, basis
import QuasiArrays: SubQuasiArray, MulQuasiMatrix, Vec, Inclusion, QuasiDiagonal, LazyQuasiArrayApplyStyle, LazyQuasiArrayStyle
import LazyArrays: MemoryLayout, ApplyStyle, Applied, colsupport, arguments, ApplyLayout, LdivStyle, MulStyle

@testset "AlephInfinity" begin
    @test !isone(ℵ₁)
    @test !iszero(ℵ₁)
    @test ℵ₁ ≠ 4
    @test ℵ₁ * ∞ == ∞ * ℵ₁ == ℵ₁
    @test 2 * ℵ₁ == ℵ₁ * 2 == ℵ₁
    @test abs(ℵ₁) == ℵ₁
    @test zero(ℵ₁) == 0
    @test 5 < ℵ₁
    @test 5 ≤ ℵ₁
    @test !(ℵ₁ ≤ 5)
    @test ℵ₁ > 5
    @test !(5 > ℵ₁)

    @test string(ℵ₁) == "ℵ₁"
end

@testset "Inclusion" begin
    x = Inclusion(-1..1)
    @test x[0.1] === 0.1
    @test x[0.0] === 0.0

    x = Inclusion(-1.0..1)
    X = QuasiDiagonal(x)
    @test X[-1:0.1:1,-1:0.1:1] == Diagonal(-1:0.1:1)
    @test Base.BroadcastStyle(typeof(x)) == LazyQuasiArrayStyle{1}()
    @test x .* x isa BroadcastQuasiArray
    @test (x.*x)[0.1] == 0.1^2

    @test exp.(x)[1.0] == exp(1.0)
    @test exp.(-x)[1.0] == exp(-1.0)
    @test exp.(x.^2)[1.0] == exp(1.0)
    @test exp.(-x.^2/(20/2))[1.0] == exp(-1.0^2/(20/2))

    @test dot(x,x) ≈ 2/3
    @test norm(x) ≈ sqrt(2/3)
end

@testset "AffineQuasiVector" begin
    x = Inclusion(0..1)
    y = 2x .- 1
    a = affine(-2..3, -1..1)

    @testset "algebra" begin
        @test 2x isa AffineQuasiVector
        @test (2x)[0.1] == 0.2
        @test_throws BoundsError (2x)[2]

        @test y isa AffineQuasiVector
        @test y[0.1] == 2*(0.1)-1
        @test y/2 isa AffineQuasiVector
        @test 2\y isa AffineQuasiVector
        @test (y/2)[0.1] == (2\y)[0.1] == -0.4
        @test y .+ 1 isa AffineQuasiVector
        @test (y .+ 1)[0.1] == (1 .+ y)[0.1]
        @test (y .- 1)[0.1] == y[0.1]-1
        @test (1 .- y)[0.1] == 1-y[0.1]
        @test y[x] == y[:] == y
        @test_throws BoundsError y[Inclusion(0..2)]
    end

    @testset "find" begin
        @test findfirst(isequal(0.1),x) == findlast(isequal(0.1),x) == 0.1
        @test findall(isequal(0.1), x) == [0.1]
        @test findfirst(isequal(2),x) == findlast(isequal(2),x) == nothing
        @test findall(isequal(2), x) == Float64[]

        @test findfirst(isequal(0.2),y) == findlast(isequal(0.2),y) == 0.6
        @test findfirst(isequal(2.3),y) == findlast(isequal(2.3),y) == nothing
        @test findall(isequal(0.2),y) == [0.6]
        @test findall(isequal(2),y) == Float64[]
    end

    @testset "minmax" begin
        @test AffineQuasiVector(x)[0.1] == 0.1
        @test minimum(y) == -1
        @test maximum(y) == 1
        @test union(y) == Inclusion(-1..1)
        @test ContinuumArrays.inbounds_getindex(y,0.1) == y[0.1]
        @test ContinuumArrays.inbounds_getindex(y,2.1) == 2*2.1 - 1

        z = 1 .- x
        @test minimum(z) == 0.0
        @test maximum(z) == 1.0
        @test union(z) == Inclusion(0..1)

        @test !isempty(z)
        @test z == z
    end

    @testset "AffineMap" begin
        @test a[0.1] == -0.16
        @test_throws BoundsError a[-3]
        @test a[-2] == first(a) == -1
        @test a[3] == last(a) == 1
        @test igetindex(a,-0.16) ≈ 0.1
        @test igetindex(a,1) == 3
        @test igetindex(a,-1) == -2
        @test union(a) == Inclusion(-1..1)

        @test affine(0..1, -1..1) == y
    end

    @testset "show" begin
        @test stringmime("text/plain", y) == "2.0 * Inclusion(0..1) .+ (-1.0)"
        @test stringmime("text/plain", a) == "0.4 * Inclusion(-2..3) .+ (-0.2)"
    end
end

@testset "DiracDelta" begin
    δ = DiracDelta(-1..3)
    @test axes(δ) === (axes(δ,1),) === (Inclusion(-1..3),)
    @test size(δ) === (length(δ),) === (ℵ₁,)
    @test δ[1.1] === 0.0
    @test δ[0.0] === Inf
    @test Base.IndexStyle(δ) === Base.IndexLinear()

    @test stringmime("text/plain", δ) == "δ at 0.0 over Inclusion(-1..3)"
    x = Inclusion(-1..1)
    @test stringmime("text/plain", δ[2x .+ 1]) == "δ at 0.0 over Inclusion(-1..3) affine mapped to -1..1"
end

@testset "Splines" begin
    @testset "HeavisideSpline" begin
        H = HeavisideSpline([1,2,3])

        @test axes(H) === (axes(H,1),axes(H,2)) === (Inclusion(1..3), Base.OneTo(2))
        @test size(H) === (size(H,1),size(H,2)) === (ℵ₁, 2)

        @test_throws BoundsError H[0.1, 1]
        @test H[1.1,1] === H'[1,1.1] === transpose(H)[1,1.1] === 1.0
        @test H[2.1,1] === H'[1,2.1] === transpose(H)[1,2.1] === 0.0
        @test H[1.1,2] === 0.0
        @test H[2.1,2] === 1.0
        @test_throws BoundsError H[2.1,3]
        @test_throws BoundsError H'[3,2.1]
        @test_throws BoundsError transpose(H)[3,2.1]
        @test_throws BoundsError H[3.1,2]

        @test all(H[[1.1,2.1], 1] .=== H'[1,[1.1,2.1]] .=== transpose(H)[1,[1.1,2.1]] .=== [1.0,0.0])
        @test all(H[1.1,1:2] .=== H[1.1,:] .=== [1.0,0.0])
        @test all(H[[1.1,2.1], 1:2] .=== [1.0 0.0; 0.0 1.0])

        @test_throws BoundsError H[[0.1,2.1], 1]
        @test MemoryLayout(typeof(H)) == BasisLayout()
        @test ApplyStyle(*, typeof(H), typeof([1,2])) isa MulStyle
        f = H*[1,2]
        @test f isa ApplyQuasiArray
        @test axes(f) == (Inclusion(1.0..3.0),)
        @test f[1.1] ≈ 1
        @test f[2.1] ≈ 2

        @test @inferred(H'H) == @inferred(materialize(applied(*,H',H))) == Eye(2)
    end

    @testset "LinearSpline" begin
        @testset "Evaluation" begin
            L = LinearSpline([1,2,3])
            @test size(L) == (ℵ₁, 3)

            @test_throws BoundsError L[0.1, 1]
            @test L[1.1,1] == L'[1,1.1] == transpose(L)[1,1.1] ≈ 0.9
            @test L[2.1,1] === L'[1,2.1] === transpose(L)[1,2.1] === 0.0
            @test L[1.1,2] ≈ 0.1
            @test L[2.1,2] ≈ 0.9
            @test L[2.1,3] == L'[3,2.1] == transpose(L)[3,2.1] ≈ 0.1
            @test_throws BoundsError L[3.1,2]

            @test L[[1.1,2.1], 1] == L'[1,[1.1,2.1]] == transpose(L)[1,[1.1,2.1]] ≈ [0.9,0.0]
            @test L[1.1,1:2] ≈ [0.9,0.1]
            @test L[[1.1,2.1], 1:2] ≈ [0.9 0.1; 0.0 0.9]

            @test_throws BoundsError L[[0.1,2.1], 1]
        end

        @testset "Expansion" begin
            L = LinearSpline([1,2,3])
            f = L*[1,2,4]

            @test basis(f) == L
            @test axes(f) == (Inclusion(1.0..3.0),)
            @test f[1.1] ≈ 1.1
            @test f[2.1] ≈ 2.2

            δ = DiracDelta(1.2,1..3)
            L = LinearSpline([1,2,3])
            @test @inferred(δ'L) ≈ [0.8, 0.2, 0.0]

            @test @inferred(L'L) == SymTridiagonal([1/3,2/3,1/3], [1/6,1/6])

            @testset "Algebra" begin
                @test 2f == f*2 == 2 .* f == f .* 2
                @test 2\f == f/2 == 2 .\ f == f ./ 2
                @test sum(f) ≈ 4.5
            end
        end

        @testset "==" begin
            L = LinearSpline([1,2,3])
            H = HeavisideSpline([1,2,3])
            @test L == L
            @test L ≠ H
            H = HeavisideSpline([1,1.5,2.5,3])
            @test_throws ArgumentError L == H
        end

        @testset "Adjoint layout" begin
            L = LinearSpline([1,2,3])
            @test MemoryLayout(L') == AdjointBasisLayout()
            @test @inferred([3,4,5]'*L') isa ApplyQuasiArray
        end

        @testset "+" begin
            L = LinearSpline([1,2,3])
            b = L*[3,4,5] + L*[1.,2,3]
            @test ApplyStyle(\, typeof(L), typeof(b)) == LdivStyle()
            @test (L\b) == [4,6,8]
            B = BroadcastQuasiArray(+, L, L)
            @test L\B == 2Eye(3)

            b = L*[3,4,5] - L*[1.,2,3]
            @test (L\b) == [2,2,2]
            B = BroadcastQuasiArray(-, L, L)
            @test L\B == 0Eye(3)
        end
    end

    @testset "Algebra" begin
        L = LinearSpline([1,2,3])
        f = L*[1,2,4]
        g = L*[5,6,7]

        @test f isa Expansion
        @test 2f isa Expansion
        @test f*2 isa Expansion
        @test 2\f isa Expansion
        @test f/2 isa Expansion
        @test f+g isa Expansion
        @test f-g isa Expansion
        @test f[1.2] == 1.2
        @test (2f)[1.2] == (f*2)[1.2] == 2.4
        @test (2\f)[1.2] == (f/2)[1.2] == 0.6
        @test (f+g)[1.2] ≈ f[1.2] + g[1.2]
        @test (f-g)[1.2] ≈ f[1.2] - g[1.2]
    end

    @testset "Derivative" begin
        L = LinearSpline([1,2,3])
        f = L*[1,2,4]

        D = Derivative(axes(L,1))
        @test D*L isa MulQuasiMatrix
        @test length((D*L).args) == 2
        @test eltype(D*L) == Float64

        M = applied(*, (D*L).args..., [1,2,4])
        @test eltype(materialize(M)) == Float64

        M = applied(*, D, L, [1,2,4])
        @test materialize(M) isa ApplyQuasiArray

        fp = D*L*[1,2,4]

        @test eltype(fp) == Float64

        @test fp isa Vec
        @test length(fp.args) == 2
        @test fp[1.1] ≈ 1
        @test fp[2.2] ≈ 2

        fp = D*f
        @test length(fp.args) == 2
        @test fp[1.1] ≈ 1
        @test fp[2.2] ≈ 2
    end

    @testset "Weak Laplacian" begin
        H = HeavisideSpline(0:2)
        L = LinearSpline(0:2)
        D = Derivative(axes(L,1))

        @test apply(*,L',D') isa MulQuasiMatrix
        @test MemoryLayout(typeof(L')) isa AdjointBasisLayout
        @test (L'D') isa MulQuasiMatrix

        A = (L'D') * (D*L)
        @test A isa BandedMatrix
        @test A == (D*L)'*(D*L) == [1.0 -1 0; -1.0 2.0 -1.0; 0.0 -1.0 1.0]
        @test bandwidths(A) == (1,1)
    end

    @testset "Views" begin
        L = LinearSpline(0:2)
        @test view(L,0.1,1)[1] == L[0.1,1]

        L = LinearSpline(0:2)
        B1 = view(L,:,1)
        @test B1 isa SubQuasiArray{Float64,1}
        @test size(B1) == (ℵ₁,)
        @test B1[0.1] == L[0.1,1]
        @test_throws BoundsError B1[2.2]

        B = view(L,:,1:2)
        @test B isa SubQuasiArray{Float64,2}
        @test B[0.1,:] == L[0.1,1:2]

        B = @view L[:,2:end-1]
        @test B[0.1,:] == [0.1]

        L = LinearSpline([1,2,3,4])
        @test L[:,2:3] isa SubQuasiArray
        @test axes(L[:,2:3]) ≡ (Inclusion(1..4), Base.OneTo(2))
        @test L[:,2:3][1.1,1] == L[1.1,2]
        @test_throws BoundsError L[0.1,1]
        @test_throws BoundsError L[1.1,0]

        @test MemoryLayout(typeof(L[:,2:3])) isa SubBasisLayout
        @test L\L[:,2:3] isa BandedMatrix
        @test L\L[:,2:3] == [0 0; 1 0; 0 1.0; 0 0]
        @test L[:,2:3]\L isa BandedMatrix
        @test L[:,2:3]\L == [0 1 0 0; 0 0 1 0]

        @testset "Subindex of splines" begin
            L = LinearSpline(range(0,stop=1,length=10))
            @test L[:,2:end-1] isa SubQuasiArray
            @test L[:,2:end-1][0.1,1] == L[0.1,2]
            v = randn(8)
            f = L[:,2:end-1] * v
            @test f[0.1] ≈ (L*[0; v; 0])[0.1]
        end

        @testset "sub-of-sub" begin
            L = LinearSpline([1,2,3])
            V = view(L,:,1:2)
            V2 = view(V,1.1:0.1:2,:)
            @test V2 == L[1.1:0.1:2,1:2]
        end

        @testset "sub-colon" begin
            L = LinearSpline([1,2,3])
            @test L[:,1][1.1] == (L')[1,:][1.1] == L[1.1,1]
            @test L[:,1:2][1.1,:] == (L')[1:2,:][:,1.1] == L[1.1,1:2]
        end

        @testset "transform" begin
            L = LinearSpline([1,2,3])
            x = axes(L,1)
            @test (L \ x) == [1,2,3]
            @test factorize(L[:,2:end-1]) isa ContinuumArrays.ProjectionFactorization
            @test L[:,1:2] \ x == [1,2]

            L = LinearSpline(range(0,1; length=10_000))
            x = axes(L,1)
            @test L[0.123,:]'* (L \ exp.(x)) ≈ exp(0.123) atol=1E-9
            @test L[0.123,2:end-1]'* (L[:,2:end-1] \ exp.(x)) ≈ exp(0.123) atol=1E-9
        end
    end

    @testset "sum" begin
        H = HeavisideSpline([1,2,3,6])
        @test sum(H; dims=1) * [1,1,1] == [5]
        x = Inclusion(0..1)
        B = H[5x .+ 1,:]
        @test sum(B; dims=1) * [1,1,1] == [1]
        @test sum(H[:,1:2]; dims=1) * [1,1] == [2]

        u = H * randn(3)
        @test sum(u[5x .+ 1]) ≈ sum(view(u,5x .+ 1)) ≈ sum(u)/5
    end

    @testset "Poisson" begin
        L = LinearSpline(range(0,stop=1,length=10))
        B = L[:,2:end-1] # Zero dirichlet by dropping first and last spline
        @test B'B == (L'L)[2:end-1,2:end-1]
        D = Derivative(axes(L,1))
        @test apply(*,D,B) isa MulQuasiMatrix
        @test D*B isa MulQuasiMatrix
        @test apply(*,D,B)[0.1,1] == (D*B)[0.1,1] == 9
        @test apply(*,D,B)[0.2,1] == (D*B)[0.2,1] == -9
        Δ = -((D*B)'*(D*B)) # Weak Laplacian
        @test Δ isa BandedMatrix

        @test B'D' isa MulQuasiMatrix
        @test length((B'D').args) == 2

        @test *(B',D',D,B) isa BandedMatrix

        @test Δ == -(*(B',D',D,B))
        @test Δ == -(B'D'D*B)
        @test Δ == -((B'D')*(D*B))
        @test_broken Δ == -B'*(D'D)*B
        @test Δ == -(B'*(D'D)*B)

        f = L*exp.(L.points) # project exp(x)
        u = B * (Δ \ (B'f))

        @test u[0.1] ≈ -0.06612902692412974
    end

    @testset "Helmholtz" begin
        L = LinearSpline(range(0,stop=1,length=10))
        B = L[:,2:end-1] # Zero dirichlet by dropping first and last spline
        D = Derivative(axes(L,1))

        A = -((B'D')*(D*B)) + 100^2*B'B # Weak Laplacian

        f = L*exp.(L.points) # project exp(x)
        u = B * (A \ (B'f))

        @test u[0.1] ≈ 0.00012678835289369413
    end

    @testset "Change-of-variables" begin
        x = Inclusion(0..1)
        y = 2x .- 1
        L = LinearSpline(range(-1,stop=1,length=10))
        @test L[y,:][0.1,:] == L[2*0.1-1,:]
        @test L[y,:]'L[y,:] isa SymTridiagonal
        @test L[y,:]'L[y,:] == 1/2*(L'L)

        D = Derivative(axes(L,1))
        H = HeavisideSpline(L.points)
        @test H\((D*L) * 2) ≈ (H\(D*L))*2 ≈ diagm(0 => fill(-9,9), 1 => fill(9,9))[1:end-1,:]

        @test MemoryLayout(L[y,:]) isa MappedBasisLayout
        @test MemoryLayout(L[y,:]') isa AdjointMappedBasisLayout
        a,b = arguments((D*L)[y,:])
        @test H[y,:]\a == Eye(9)
        @test H[y,:] \ (D*L)[y,:] isa BandedMatrix
        @test @inferred(grid(L[y,:])) ≈ (grid(L) .+ 1) ./ 2

        D = Derivative(x)
        @test (D*L[y,:])[0.1,1] ≈ -9
        @test H[y,:] \ (D*L[y,:]) isa BandedMatrix
        @test H[y,:] \ (D*L[y,:]) ≈ diagm(0 => fill(-9,9), 1 => fill(9,9))[1:end-1,:]

        B = L[y,2:end-1]
        @test MemoryLayout(typeof(B)) isa MappedBasisLayout
        @test B[0.1,1] == L[2*0.1-1,2]
        @test B\B == Eye(8)
        @test L[y,:] \ B == Eye(10)[:,2:end-1]
        @test B'B == (1/2)*(L'L)[2:end-1,2:end-1]

        @testset "algebra" begin
            f = L[y,:] * randn(10)
            g = L[y,:] * randn(10)
            @test f + g isa Expansion
            @test (f+g)[0.1] ≈ f[0.1] + g[0.1]
        end
    end

    @testset "diff" begin
        L = LinearSpline(range(-1,stop=1,length=10))
        f = L * randn(size(L,2))
        h = 0.0001;
        @test diff(f)[0.1] ≈ (f[0.1+h]-f[0.1])/h
    end

    @testset "show" begin
        x = Inclusion(0..1)
        H = HeavisideSpline([1,2,3,6])
        B = H[5x .+ 1,:]
        u = H * [1,2,3]
        @test stringmime("text/plain", B) == "Spline{0,Float64,Array{Int64,1}}([1, 2, 3, 6]) affine mapped to 0..1"
    end
end

@testset "Kernels" begin
    x = Inclusion(0..1)
    K = x .- x'
    @test K[0.1,0.2] == K[Inclusion(0..0.5), Inclusion(0..0.5)][0.1,0.2] == 0.1 - 0.2
    @test_throws BoundsError K[Inclusion(0..0.5), Inclusion(0..0.5)][1,1]
end

struct Chebyshev <: Basis{Float64}
    n::Int
end

struct ChebyshevWeight <: Weight{Float64} end

Base.:(==)(::Chebyshev, ::Chebyshev) = true
Base.axes(T::Chebyshev) = (Inclusion(-1..1), Base.OneTo(T.n))
ContinuumArrays.grid(T::Chebyshev) = chebyshevpoints(Float64, T.n, Val(1))
Base.axes(T::ChebyshevWeight) = (Inclusion(-1..1),)

Base.getindex(::Chebyshev, x::Float64, n::Int) = cos((n-1)*acos(x))
Base.getindex(::ChebyshevWeight, x::Float64) = 1/sqrt(1-x^2)

LinearAlgebra.factorize(L::Chebyshev) =
    TransformFactorization(grid(L), plan_chebyshevtransform(Array{Float64}(undef, size(L,2))))

# This is wrong but just for tests
Base.broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::Expansion{<:Any,<:Chebyshev}, b::Chebyshev) = b * Matrix(I, 5, 5)

@testset "Chebyshev" begin
    T = Chebyshev(5)
    w = ChebyshevWeight()
    wT = w .* T
    x = axes(T,1)
    F = factorize(T)
    g = grid(F)
    @test T \ exp.(x) == F \ exp.(x) == F \ exp.(g) == chebyshevtransform(exp.(g), Val(1))

    @testset "Weighted" begin
        @test MemoryLayout(w) isa WeightLayout
        @test MemoryLayout(w[Inclusion(0..1)]) isa WeightLayout

        wT2 = w .* T[:,2:4]
        wT3 = wT[:,2:4]
        @test MemoryLayout(wT) == WeightedBasisLayout()
        @test MemoryLayout(wT2) == WeightedBasisLayout()
        @test MemoryLayout(wT3) == SubWeightedBasisLayout()
        @test grid(wT) == grid(wT2) == grid(wT3) == grid(T)

        @test ContinuumArrays.unweightedbasis(wT) ≡ T
        @test ContinuumArrays.unweightedbasis(wT2) ≡ T[:,2:4]
        @test ContinuumArrays.unweightedbasis(wT3) ≡ T[:,2:4]
    end
    @testset "Mapped" begin
        y = affine(0..1, x)
        @test MemoryLayout(wT[y,:]) isa MappedWeightedBasisLayout
        @test MemoryLayout(w[y] .* T[y,:]) isa MappedWeightedBasisLayout
        @test wT[y,:][[0.1,0.2],1:5] == (w[y] .* T[y,:])[[0.1,0.2],1:5] == (w .* T[:,1:5])[y,:][[0.1,0.2],:]
    end

    @testset "Broadcasted" begin
        a = 1 .+ x .+ x.^2
        # The following are wrong, just testing dispatch
        @test T \ (a .* T) == I
        @test T \ (a .* (T * (T \ a))) ≈ [2.875, 3.5, 2.0, 0.5, 0.125]
        f = exp.(x) .* a # another broadcast layout
        @test T \ f == F \ f

        ã = T * (T \ a)
        @test T \ (ã .* ã) ≈ [1.5,1,0.5,0,0]
    end
end


text:
Shape = Vector{Int}
Tableau = Vector{Vector{Int}}
Tableaux = Vector{Tableau}

function standard_tableaux(table::Vector{Dict{Shape, Tableaux}}, t::Tableau, max_n::Int)
    shape = length.(t)
    n = sum(shape)

    if length(table) < n
        push!(table, Dict{Shape, Tableaux}())
    end

    push!(get!(table[n], shape, Tableau[]), t)

    n == max_n && return table

    for row in 1:length(t)
        if row == 1 || length(t[row]) < length(t[row - 1])
            new_t = deepcopy(t)
            push!(new_t[row], n + 1)
            table = standard_tableaux(table, new_t, max_n)
        end

        if row == length(t)
            new_t = deepcopy(t)
            push!(new_t, [n + 1])
            table = standard_tableaux(table, new_t, max_n)
        end
    end

    return table
end

standard_tableaux(max_n::Int) = standard_tableaux(Dict{Shape, Tableaux}[], [[1]], max_n)

function pos(t::Tableau, i::Int)
    for r in 1:length(t), c in 1:length(t[r])
        if t[r][c] == i
            return r, c
        end
    end
    return 0, 0
end

function adjacent_ρ(table::Vector{Dict{Shape, Tableaux}}, λ::Shape, i::Int)
    n = sum(λ)
    d_λ = length(table[n][λ])
    ρ = zeros(Float64, (d_λ, d_λ))

    for (k, t) in enumerate(table[n][λ])
        r_i, c_i = pos(t, i)
        r_j, c_j = pos(t, i - 1)
        d = (c_i - c_j) - (r_i - r_j)
        ρ[k, k] = 1 / d
        if r_i != r_j && c_i != c_j
            new_t = deepcopy(t)
            new_t[r_i][c_i] = i - 1
            new_t[r_j][c_j] = i
            ℓ = findfirst(table[n][λ], new_t)
            ρ[ℓ, k] = √(1 - 1/d^2)
        end
    end

    return ρ
end

function yor(table::Vector{Dict{Shape, Tableaux}}, λ::Shape, σ::Vector{Int})
    n = length(σ)
    d_λ = length(table[n][λ])

    ν = deepcopy(σ)
    factors = Int[]
    for i in 1:n, j in n:-1:i+1
        if ν[j] < ν[j - 1]
            ν[j-1], ν[j] = ν[j], ν[j-1]
            push!(factors, j)
        end
    end

    ρ = eye(d_λ)
    for k in factors
        ρ = adjacent_ρ(table, λ, k) * ρ
    end

    return ρ
end

function irrep{N}(g::SymmetricGroup{N}, λ::Vector{Int})
    table = standard_tableaux(N)
    yor(table, λ, g.σ)
end

real_irrep{N}(g::SymmetricGroup{N}, λ::Vector{Int}) = irrep(g, λ)


text:
# Newton and quasi Newton solvers

## LOGISTIC + 0/L2 and Huber ==============

"""
$SIGNATURES

Fit a GLR using Newton's method.

## Complexity

Assuming `n` dominates `p`, O(κnp²), dominated by the construction of the
Hessian at each step with κ the number of Newton steps.
"""
function _fit(glr::GLR{<:Union{LogisticLoss,RobustLoss},<:L2R},
              solver::Newton, X, y)
    p     = size(X, 2) + Int(glr.fit_intercept)
    θ₀    = zeros(p)
    _fgh! = fgh!(glr, X, y)
    opt   = Optim.only_fgh!(_fgh!)
    res   = Optim.optimize(opt, θ₀, Optim.Newton())
    return Optim.minimizer(res)
end

"""
$SIGNATURES

Fit a GLR using Newton's method combined with an iterative solver  (conjugate
gradient) to solve the Newton steps (∇²f)⁻¹∇f.

## Complexity

Assuming `n` dominates `p`, O(κ₁κ₂np), dominated by the application of the
Hessian at each step where κ₁ is the number of Newton steps and κ₂ is the
average number of CG steps per Newton step (which is at most p).
"""
function _fit(glr::GLR{<:Union{LogisticLoss,RobustLoss},<:L2R},
              solver::NewtonCG, X, y)
    p    = size(X, 2) + Int(glr.fit_intercept)
    θ₀   = zeros(p)
    _f   = objective(glr, X, y)
    _fg! = (g, θ) -> fgh!(glr, X, y)(0.0, g, nothing, θ) # Optim.jl/issues/738
    _Hv! = Hv!(glr, X, y)
    opt  = Optim.TwiceDifferentiableHV(_f, _fg!, _Hv!, θ₀)
    res  = Optim.optimize(opt, θ₀, Optim.KrylovTrustRegion())
    return Optim.minimizer(res)
end

"""
$SIGNATURES

Fit a GLR using LBFGS.

## Complexity

Assuming `n` dominates `p`, O(κnp), dominated by the computation of the
gradient at each step with κ the number of LBFGS steps.
"""
function _fit(glr::GLR{<:Union{LogisticLoss,RobustLoss},<:L2R},
              solver::LBFGS, X, y)
    p    = size(X, 2) + Int(glr.fit_intercept)
    θ₀   = zeros(p)
    _fg! = (f, g, θ) -> fgh!(glr, X, y)(f, g, nothing, θ)
    opt  = Optim.only_fg!(_fg!)
    res  = Optim.optimize(opt, θ₀, Optim.LBFGS())
    return Optim.minimizer(res)
end

## MULTINOMIAL + 0/L2 ==============

"""
$SIGNATURES

Fit a multiclass GLR using Newton's method with an iterative solver (conjugate
gradient).

## Complexity

Assuming `n` dominates `p`, O(κ₁κ₂npc), where `c` is the number of classes. The
computations are dominated by the application of the Hessian at each step with
κ₁ the number of Newton steps and κ₂ the average number of CG steps per Newton
step.
"""
function _fit(glr::GLR{MultinomialLoss,<:L2R}, solver::NewtonCG, X, y)
    p    = size(X, 2) + Int(glr.fit_intercept)
    c    = maximum(y)
    θ₀   = zeros(p * c)
    _f   = objective(glr, X, y; c=c)
    _fg! = (g, θ) -> fg!(glr, X, y)(0.0, g, θ) # XXX: Optim.jl/issues/738
    _Hv! = Hv!(glr, X, y)
    opt  = Optim.TwiceDifferentiableHV(_f, _fg!, _Hv!, θ₀)
    res  = Optim.optimize(opt, θ₀, Optim.KrylovTrustRegion())
    return Optim.minimizer(res)
end

"""
$SIGNATURES

Fit a multiclass GLR using LBFGS.

## Complexity

Assuming `n` dominates `p`, O(κnpc), with `c` the number of classes, dominated
by the computation of the gradient at each step with κ the number of LBFGS
steps.
"""
function _fit(glr::GLR{MultinomialLoss,<:L2R}, solver::LBFGS, X, y)
    p    = size(X, 2) + Int(glr.fit_intercept)
    c    = maximum(y)
    θ₀   = zeros(p * c)
    _fg! = fg!(glr, X, y)
    opt  = Optim.only_fg!(_fg!)
    res  = Optim.optimize(opt, θ₀, Optim.LBFGS())
    return Optim.minimizer(res)
end