text:
function fit{C}(m::GaussianNB, X::Matrix{Float64}, y::Vector{C})
    ensure_data_size(X, y)
    # updatestats(m.dstats, X)
    # m.gaussian = MvNormal(mean(m.dstats), cov(m.dstats))
    # m.n_obs = m.dstats.n_obs
    n_vars = size(X, 1)
    for j=1:size(X, 2)
        c = y[j]
        m.c_counts[c] += 1
        updatestats(m.c_stats[c], reshape(X[:, j], n_vars, 1))
        # m.x_counts[c] .+= X[:, j]
        # m.x_totals += X[:, j]
        m.n_obs += 1
    end
    # precompute distributions for each class
    for c in keys(m.c_counts)
        m.gaussians[c] = MvNormal(mean(m.c_stats[c]), cov(m.c_stats[c]))
    end
    return m
end

"""Calculate log P(x|C)"""
function logprob_x_given_c{C}(m::GaussianNB, x::Vector{Float64}, c::C)
    return logpdf(m.gaussians[c], x)
end

"""Calculate log P(x|C)"""
function logprob_x_given_c{C}(m::GaussianNB, X::Matrix{Float64}, c::C)
    ## x_priors_for_c = m.x_counts[c] ./ m.x_totals
    ## x_probs_given_c = x_priors_for_c .^ x
    ## logprob = sum(log(x_probs_given_c))
    ## return logprob
    return logpdf(m.gaussians[c], X)
end


text:
__precompile__()

module NNLS
using LinearAlgebra
using Statistics: mean

export nnls,
       solve!,
       NNLSWorkspace,
       QPWorkspace,
       QP,
       primal_infeasibility,
       dual_infeasibility,
       stationarity_violation,
       slackness_violation,
       check_optimality_conditions,
       load!,
       NNLSSolver

"""
CONSTRUCTION AND/OR APPLICATION OF A SINGLE
HOUSEHOLDER TRANSFORMATION..     Q = I + U*(U**T)/B

The original version of this code was developed by
Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
1973 JUN 12, and published in the book
"SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
Revised FEB 1995 to accompany reprinting of the book by SIAM.
"""
function construct_householder!(u::AbstractVector{T}, up::T)::T where T
    m = length(u)
    if m <= 1
        return up
    end

    cl = maximum(abs, u)
    @assert cl > 0
    clinv = 1 / cl
    sm = zero(T)
    for ui in u
        sm += (ui * clinv)^2
    end
    cl *= sqrt(sm)
    if u[1] > 0
        cl = -cl
    end
    result = u[1] - cl
    u[1] = cl

    return result
end

"""
CONSTRUCTION AND/OR APPLICATION OF A SINGLE
HOUSEHOLDER TRANSFORMATION..     Q = I + U*(U**T)/B

The original version of this code was developed by
Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
1973 JUN 12, and published in the book
"SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
Revised FEB 1995 to accompany reprinting of the book by SIAM.
"""
function apply_householder!(u::AbstractVector{T}, up::T, c::AbstractVector{T}) where T
    m = length(u)
    if m > 1
        cl = abs(u[1])
        @assert cl > 0
        b = up * u[1]
        if b >= 0
            return
        end
        b = 1 / b

        sm = c[1] * up
        for i in 2:m
            sm += c[i] * u[i]
        end
        if sm != 0
            sm *= b
            c[1] += sm * up
            for i in 2:m
                c[i] += sm * u[i]
            end
        end
    end
end

"""
   COMPUTE ORTHOGONAL ROTATION MATRIX..
The original version of this code was developed by
Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
1973 JUN 12, and published in the book
"SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
Revised FEB 1995 to accompany reprinting of the book by SIAM.

   COMPUTE.. MATRIX   (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2))
                      (-S,C)         (-S,C)(B)   (   0          )
   COMPUTE SIG = SQRT(A**2+B**2)
      SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT
      SIG MAY BE IN THE SAME LOCATION AS A OR B .
"""
function orthogonal_rotmat(a::T, b::T)::Tuple{T, T, T} where T
    if abs(a) > abs(b)
        xr = b / a
        yr = sqrt(1 + xr^2)
        c = (1 / yr) * sign(a)
        s = c * xr
        sig = abs(a) * yr
    elseif b != 0
        xr = a / b
        yr = sqrt(1 + xr^2)
        s = (1 / yr) * sign(b)
        c = s * xr
        sig = abs(b) * yr
    else
        sig = zero(T)
        c = zero(T)
        s = one(T)
    end
    return c, s, sig
end

"""
The original version of this code was developed by
Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
1973 JUN 15, and published in the book
"SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
Revised FEB 1995 to accompany reprinting of the book by SIAM.
"""
function solve_triangular_system!(zz, A, idx, nsetp, jj)
    for l in 1:nsetp
        ip = nsetp + 1 - l
        if (l != 1)
            for ii in 1:ip
                zz[ii] -= A[ii, jj] * zz[ip + 1]
            end
        end
        jj = idx[ip]
        zz[ip] /= A[ip, jj]
    end
    return jj
end

mutable struct NNLSWorkspace{T, I <: Integer}
    QA::Matrix{T}
    Qb::Vector{T}
    x::Vector{T}
    w::Vector{T}
    zz::Vector{T}
    idx::Vector{I}
    rnorm::T
    mode::I
    nsetp::I

    function NNLSWorkspace{T, I}(m, n) where {T, I <: Integer}
        new{T, I}(Matrix{T}(undef, m, n), # A
            Vector{T}(undef, m),    # b
            Vector{T}(undef, n),    # x
            Vector{T}(undef, n),    # w
            Vector{T}(undef, m),    # zz
            Vector{I}(undef, n),    # idx
            zero(T), # rnorm
            zero(I), # mode
            zero(I)  # nsetp
        )
    end
end

function Base.resize!(work::NNLSWorkspace{T}, m::Integer, n::Integer) where T
    work.QA = Matrix{T}(undef, m, n)
    work.Qb = Vector{T}(undef, m)
    resize!(work.x, n)
    resize!(work.w, n)
    resize!(work.zz, m)
    resize!(work.idx, n)
end

function load!(work::NNLSWorkspace{T}, A::AbstractMatrix{T}, b::AbstractVector{T}) where T
    m, n = size(A)
    @assert size(b) == (m,)
    if size(work.QA, 1) != m || size(work.QA, 2) != n
        resize!(work, m, n)
    end
    work.QA .= A
    work.Qb .= b
    work
end

NNLSWorkspace(m::Integer, n::Integer,
              eltype::Type{T}=Float64,
              indextype::Type{I}=Int) where {T, I} = NNLSWorkspace{T, I}(m, n)

function NNLSWorkspace(A::Matrix{T}, b::Vector{T}, indextype::Type{I}=Int) where {T, I}
    m, n = size(A)
    @assert size(b) == (m,)
    work = NNLSWorkspace{T, I}(m, n)
    load!(work, A, b)
    work
end

"""
Views in Julia still allocate some memory (since they need to keep
a reference to the original array). This type allocates no memory
and does no bounds checking. Use it with caution.
"""
struct UnsafeVectorView{T} <: AbstractVector{T}
    offset::Int
    len::Int
    ptr::Ptr{T}
end

UnsafeVectorView(parent::DenseArray{T}, start_ind::Integer, len::Integer) where {T} = UnsafeVectorView{T}(start_ind - 1, len, pointer(parent))
Base.size(v::UnsafeVectorView) = (v.len,)
Base.getindex(v::UnsafeVectorView, idx) = unsafe_load(v.ptr, idx + v.offset)
Base.setindex!(v::UnsafeVectorView, value, idx) = unsafe_store!(v.ptr, value, idx + v.offset)
Base.length(v::UnsafeVectorView) = v.len
Base.IndexStyle(::Type{V}) where {V <: UnsafeVectorView} = Base.IndexLinear()

"""
UnsafeVectorView only works for isbitstype types. For other types, we're already
allocating lots of memory elsewhere, so creating a new View is fine.

This function looks type-unstable, but the isbitstype(T) test can be evaluated
by the compiler, so the result is actually type-stable.
"""
function fastview(parent::Array{T}, start_ind::Integer, len::Integer) where T
    if isbitstype(T)
        UnsafeVectorView(parent, start_ind, len)
    else
        @view(parent[start_ind:(start_ind + len - 1)])
    end
end

"""
Fallback for non-contiguous arrays, for which UnsafeVectorView does not make
sense.
"""
fastview(parent::AbstractArray, start_ind::Integer, len::Integer) = @view(parent[start_ind:(start_ind + len - 1)])

@noinline function checkargs(work::NNLSWorkspace)
    m, n = size(work.QA)
    @assert size(work.Qb) == (m,)
    @assert size(work.x) == (n,)
    @assert size(work.w) == (n,)
    @assert size(work.zz) == (m,)
    @assert size(work.idx) == (n,)
end

function largest_positive_dual(w::AbstractVector{T}, idx::AbstractVector{TI}, range) where {T, TI}
    wmax = zero(T)
    izmax = zero(TI)
    for i in range
        j = idx[i]
        if w[j] > wmax
            wmax = w[j]
            izmax = i
        end
    end
    wmax, izmax
end

"""
Algorithm NNLS: NONNEGATIVE LEAST SQUARES

The original version of this code was developed by
Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
1973 JUN 15, and published in the book
"SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
Revised FEB 1995 to accompany reprinting of the book by SIAM.

GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B,  COMPUTE AN
N-VECTOR, X, THAT SOLVES THE LEAST SQUARES PROBLEM
                 A * X = B  SUBJECT TO X .GE. 0
"""
function solve!(work::NNLSWorkspace{T, TI}, max_iter::Integer=(3 * size(work.QA, 2))) where {T, TI}
    checkargs(work)

    A = work.QA
    b = work.Qb
    x = work.x
    w = work.w
    zz = work.zz
    idx = work.idx
    factor = 0.01
    work.mode = 1

    m = convert(TI, size(A, 1))
    n = convert(TI, size(A, 2))

    iter = 0
    x .= 0
    idx .= 1:n

    iz2 = n
    iz1 = one(TI)
    iz = zero(TI)
    j = zero(TI)
    jj = zero(TI)
    nsetp = zero(TI)
    up = zero(T)

    terminated = false

    # ******  MAIN LOOP BEGINS HERE  ******
    while true
        # println("jl main loop")
        # QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION.
        # OR IF M COLS OF A HAVE BEEN TRIANGULARIZED.
        if (iz1 > iz2 || nsetp >= m)
            terminated = true
            break
        end

        # COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W().
        for i in iz1:iz2
            idxi = idx[i]
            sm = zero(T)
            for l in (nsetp + 1):m
                sm += A[l, idxi] * b[l]
            end
            w[idxi] = sm
        end

        while true
            # FIND LARGEST POSITIVE W(J).
            wmax, izmax = largest_positive_dual(w, idx, iz1:iz2)

            # IF WMAX .LE. 0. GO TO TERMINATION.
            # THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS.
            if wmax <= 0
                terminated = true
                break
            end

            iz = izmax
            j = idx[iz]

            # THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P.
            # BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID
            # NEAR LINEAR DEPENDENCE.
            Asave = A[nsetp + 1, j]
            up = construct_householder!(
                fastview(A, LinearIndices(A)[nsetp + 1, j], m - nsetp),
                up)
            unorm::T = zero(T)
            for l in 1:nsetp
                unorm += A[l, j]^2
            end
            unorm = sqrt(unorm)

            if ((unorm + abs(A[nsetp + 1, j]) * factor) - unorm) > 0
                # COL J IS SUFFICIENTLY INDEPENDENT.  COPY B INTO ZZ, UPDATE ZZ
                # AND SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ).
                # println("copying b into zz")
                zz .= b
                apply_householder!(
                    fastview(A, LinearIndices(A)[nsetp + 1, j], m - nsetp),
                    up,
                    fastview(zz, nsetp + 1, m - nsetp))
                ztest = zz[nsetp + 1] / A[nsetp + 1, j]

                # SEE IF ZTEST IS POSITIVE
                if ztest > 0
                    break
                end
            end

            # REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P.
            # RESTORE A(NPP1,J), SET W(J)=0., AND LOOP BACK TO TEST DUAL
            # COEFFS AGAIN.
            A[nsetp + 1, j] = Asave
            w[j] = 0
        end
        if terminated
            break
        end

        # THE INDEX  J=INDEX(IZ)  HAS BEEN SELECTED TO BE MOVED FROM
        # SET Z TO SET P.    UPDATE B,  UPDATE INDICES,  APPLY HOUSEHOLDER
        # TRANSFORMATIONS TO COLS IN NEW SET Z,  ZERO SUBDIAGONAL ELTS IN
        # COL J,  SET W(J)=0.
        b .= zz

        idx[iz] = idx[iz1]
        idx[iz1] = j
        iz1 += one(TI)
        nsetp += one(TI)

        if iz1 <= iz2
            for jz in iz1:iz2
                jj = idx[jz]
                apply_householder!(
                    fastview(A, LinearIndices(A)[nsetp, j], m - nsetp + 1),
                    up,
                    fastview(A, LinearIndices(A)[nsetp, jj], m - nsetp + 1))
            end
        end

        if nsetp != m
            for l in (nsetp + 1):m
                A[l, j] = 0
            end
        end

        w[j] = 0

        # SOLVE THE TRIANGULAR SYSTEM.
        # STORE THE SOLUTION TEMPORARILY IN ZZ().
        jj = solve_triangular_system!(zz, A, idx, nsetp, jj)

        # ******  SECONDARY LOOP BEGINS HERE ******
        #
        # ITERATION COUNTER.
        while true
            iter += 1
            if iter > max_iter
                work.mode = 3
                terminated = true
                println("NNLS quitting on iteration count")
                break
            end

            # SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE.
            # IF NOT COMPUTE ALPHA.
            alpha = convert(T, 2)
            for ip in Base.OneTo(nsetp)
                l = idx[ip]
                if zz[ip] <= 0
                    t = -x[l] / (zz[ip] - x[l])
                    if alpha > t
                        alpha = t
                        jj = ip
                    end
                end
            end

            # IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL
            # STILL = 2.    IF SO EXIT FROM SECONDARY LOOP TO MAIN LOOP.
            if alpha == 2
                break
            end

            # OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0 AND 1 TO
            # INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ.
            for ip in Base.OneTo(nsetp)
                l = idx[ip]
                x[l] = x[l] + alpha * (zz[ip] - x[l])
            end

            # MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I
            # FROM SET P TO SET Z.
            i = idx[jj]

            while true
                x[i] = 0

                if jj != nsetp
                    jj += one(TI)
                    for j in jj:nsetp
                        ii = idx[j]
                        idx[j - 1] = ii
                        cc, ss, sig = orthogonal_rotmat(A[j - 1, ii], A[j, ii])
                        A[j - 1, ii] = sig
                        A[j, ii] = 0
                        for l in Base.OneTo(n)
                            if l != ii
                                # Apply procedure G2 (CC,SS,A(J-1,L),A(J,L))
                                temp = A[j - 1, l]
                                A[j - 1, l] = cc * temp + ss * A[j, l]
                                A[j, l] = -ss * temp + cc * A[j, l]
                            end
                        end

                        # Apply procedure G2 (CC,SS,B(J-1),B(J))
                        temp = b[j - 1]
                        b[j - 1] = cc * temp + ss * b[j]
                        b[j] = -ss * temp + cc * b[j]
                    end
                end

                nsetp -= one(TI)
                iz1 -= one(TI)
                idx[iz1] = i

                # SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE.  THEY SHOULD
                # BE BECAUSE OF THE WAY ALPHA WAS DETERMINED.
                # IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR.  ANY
                # THAT ARE NONPOSITIVE WILL BE SET TO ZERO
                # AND MOVED FROM SET P TO SET Z.
                allfeasible = true
                for jji in Base.OneTo(nsetp)
                    jj = jji
                    i = idx[jj]
                    if x[i] <= 0
                        allfeasible = false
                        break
                    end
                end
                if allfeasible
                    break
                end
            end

            # COPY B( ) INTO ZZ( ).  THEN SOLVE AGAIN AND LOOP BACK.
            zz .= b
            jj = solve_triangular_system!(zz, A, idx, nsetp, jj)
        end
        if terminated
            break
        end
        # ******  END OF SECONDARY LOOP  ******

        for i in 1:nsetp
            x[idx[i]] = zz[i]
        end
        # ALL NEW COEFFS ARE POSITIVE.  LOOP BACK TO BEGINNING.
    end

    # ******  END OF MAIN LOOP  ******
    # COME TO HERE FOR TERMINATION.
    # COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR.

    sm = zero(T)
    if nsetp < m
        for i in (nsetp + 1):m
            sm += b[i]^2
        end
    else
        w .= 0
    end
    work.rnorm = sqrt(sm)
    work.nsetp = nsetp
    return work.x
end

function solve!(work::NNLSWorkspace{T}, A::AbstractMatrix{T}, b::AbstractVector{T}, max_iter=(3 * size(A, 2))) where T
    load!(work, A, b)
    solve!(work, max_iter)
    work.x
end

function nnls(A::DenseMatrix{T}, b::DenseVector{T}, max_iter=(3 * size(A, 2))) where T
    work = NNLSWorkspace(A, b)
    solve!(work, max_iter)
    work.x
end

# Implementation of an NNLS-based QP solver, based on section II of  Bemporad,
# "A quadratic programming algorithm based on nonnegative least squares with
# applications to embedded model predictive control", IEEE Transactions on
# Automatic Control, 2016.
# Variable names match the paper wherever possible

const AllColsSubArray{T} = SubArray{T,2,Array{T,2},Tuple{UnitRange{Int},Base.Slice{Base.OneTo{Int}}},false}

"""
Structure describing the QP:

Minimize ``\\frac{1}{2} z' Q z + c' z``
Subject to ``G z \\leq g``
"""
struct QP{T}
    Q::Matrix{T}
    c::Vector{T}
    G::Matrix{T}
    g::Vector{T}

    function QP{T}(Q::Matrix, c::Vector, G::Matrix, g::Vector) where T
        @boundscheck begin
            LinearAlgebra.checksquare(Q)
            length(c) == size(Q, 1) || throw(DimensionMismatch())
            size(G, 1) == length(g) || throw(DimensionMismatch())
            size(G, 2) == size(Q, 1) || throw(DimensionMismatch())
        end
        new{T}(Q, c, G, g)
    end

    QP(Q::Matrix{T}, c::Vector{T}, G::Matrix{T}, g::Vector{T}) where {T} = QP{T}(Q, c, G, g)
    QP{T}(qp::QP) where {T} = QP{T}(qp.Q, qp.c, qp.G, qp.g)
end

mutable struct QPWorkspace{T, I} # TODO: consider not having `I` as a parameter
    # Variables from paper:
    L::Matrix{T}
    c::Vector{T}
    G::Matrix{T}
    g::Vector{T}
    M::Matrix{T}
    d::Vector{T}
    r::Vector{T}

    # Additional variables:
    e::Vector{T} # intermediate variable
    A::Matrix{T} # 'A'-matrix of NNLS problem
    AM::AllColsSubArray{T} # upper block of A
    Ad::AllColsSubArray{T} # last row of A
    b::Vector{T} # 'b'-vector of NNLS problem
    nnlswork::NNLSWorkspace{T, I}
    status::Symbol

    function QPWorkspace{T, I}(q::Integer, n::Integer) where {T, I}
        work = new{T, I}()
        resize!(work, q, n)
    end
end

QPWorkspace(q::Integer, n::Integer) = QPWorkspace{Float64, Int}(q, n)

function Base.convert(::Type{QP{T1}}, qp::QP{T2}) where {T1, T2}
    QP(convert(Matrix{T1}, qp.Q),
       convert(Vector{T1}, qp.c),
       convert(Matrix{T1}, qp.G),
       convert(Vector{T1}, qp.g))
end

"""
    QPWorkspace(qp::QP)

Construct a workspace and load problem data for the QP

Minimize ``\\frac{1}{2} z' Q z + c' z``
Subject to ``G z \\leq g``
"""
function QPWorkspace(qp::QP{T}) where T
    work = QPWorkspace{T, Int}(size(qp.G, 1), size(qp.G, 2))
    load!(work, qp.Q, qp.c, qp.G, qp.g)
    work
end

function Base.resize!(work::QPWorkspace{T}, q::Integer, n::Integer) where T
    work.L = Matrix{T}(undef, n, n)
    work.c = Vector{T}(undef, n)
    work.G = Matrix{T}(undef, q, n)
    work.g = Vector{T}(undef, q)
    work.M = Matrix{T}(undef, q, n)
    work.d = Vector{T}(undef, q)
    work.r = Vector{T}(undef, n + 1)
    work.e = Vector{T}(undef, n)
    work.A = Matrix{T}(undef, n + 1, q)
    work.AM = view(work.A, 1 : n, :)
    work.Ad = view(work.A, n + 1 : n + 1, :)
    work.b = Vector{T}(undef, n + 1)
    work.nnlswork = NNLSWorkspace{T, Int}(size(work.A)...)
    work.status = :Unsolved
    work
end

"""
    load!(work::QPWorkspace, Q, c, G, g)

Load problem data for the QP

Minimize ``\\frac{1}{2} z' Q z + c' z``
Subject to ``G z \\leq g``
"""
function load!(work::QPWorkspace{T}, Q::AbstractMatrix{T}, c::AbstractVector{T}, G::AbstractMatrix{T}, g::AbstractVector{T}) where T
    work.L .= Q
    work.c .= c
    work.G .= G
    work.g .= g
    work.status = :Unsolved
    nothing
end

load!(work::QPWorkspace{T}, qp::QP{T}) where {T} = load!(work, qp.Q, qp.c, qp.G, qp.g)
checkunsolved(work::QPWorkspace) = work.status == :Unsolved || error("Problem was already solved.")

"""
    z, λ = solve!(work::QPWorkspace)

Solve the QP that was loaded into `work` using `load!`. Returns the primal
solution ``z`` and the dual solution ``λ``.
"""
function solve!(work::QPWorkspace{T}, eps_infeasible = 1e-4) where T<:LinearAlgebra.BlasFloat
    checkunsolved(work)

    L = work.L
    c = work.c
    G = work.G
    g = work.g
    M = work.M
    d = work.d
    e = work.e
    A = work.A
    AM = work.AM
    Ad = work.Ad
    b = work.b
    r = work.r
    nnlswork = work.nnlswork

    # Compute M
    LinearAlgebra.LAPACK.potrf!('U', L) # L <- upper cholesky factor of Q
    M .= G
    LinearAlgebra.BLAS.trsm!('R', 'U', 'N', 'N', one(T), L, M) # M <- G L⁻¹

    # Compute d
    e .= c
    LinearAlgebra.LAPACK.potrs!('U', L, e) # e <- Q⁻¹ c

    d .= g
    LinearAlgebra.BLAS.gemv!('N', one(T), G, e, one(T), d) # d <- g + G Q⁻¹ c

    # Populate A
    transpose!(AM, M)
    transpose!(Ad, d)
    rmul!(A, -1)

    # Populate b
    γ = one(T)
    b .= 0
    b[end] = γ

    # Solve the NNLS
    load!(nnlswork, A, b)
    solve!(nnlswork)
    y = nnlswork.x

    # Compute the residual
    r = b
    LinearAlgebra.BLAS.gemv!('N', one(T), A, y, -one(T), r) # r <- A * y - b

    # Check for feasibility
    work.status = sum(abs, r) < eps_infeasible ? :Infeasible : :Optimal

    # Back out solution and dual
    z = c
    λ = y
    if work.status == :Optimal
        # Note: r[end] == -(γ + d ⋅ y)
        LinearAlgebra.BLAS.gemv!('T', 1 / r[end], G, y, -one(T), c) # z <- -1 / (γ + d ⋅ y) G^ᵀ y - c
        LinearAlgebra.LAPACK.potrs!('U', L, z)
        rmul!(λ, -1 / sqrt(-r[end])) # the sqrt appears to be missing in (12) in the paper
    else
        fill!(z, NaN)
        fill!(λ, NaN)
    end

    z, λ
end

function solve!(work::QPWorkspace{T}, eps_infeasible = 1e-4) where T
    checkunsolved(work)

    Q = work.L # is set to Q initially; modified to be L
    c = work.c
    G = work.G
    g = work.g
    M = work.M
    d = work.d
    e = work.e
    A = work.A
    AM = work.AM
    Ad = work.Ad
    b = work.b
    r = work.r
    nnlswork = work.nnlswork

    # Compute M
    cholQ = cholesky!(Hermitian(Q, :U))
    L = cholQ.U
    M[:] = G
    rdiv!(M, L)

    # Compute d
    e[:] = c
    ldiv!(cholQ, e) # e <- Q⁻¹ c
    mul!(d, G, e)
    d .+= g # d <- g + G Q⁻¹ c

    # Populate A
    transpose!(AM, M)
    transpose!(Ad, d)
    rmul!(A, -1)

    # Populate b
    γ = one(T)
    b .= 0
    b[end] = γ

    # Solve the NNLS
    load!(nnlswork, A, b)
    solve!(nnlswork)
    y = nnlswork.x

    # Compute the residual
    mul!(r, A, y)
    r .-= b # r <- A * y - b

    # Check for feasibility
    work.status = sum(abs, r) < eps_infeasible ? :Infeasible : :Optimal

    # Back out solution and dual
    z = similar(c) # TODO
    λ = y
    if work.status == :Optimal
        # Note: r[end] == -(γ + d ⋅ y)
        mul!(z, transpose(G), y)
        z .= z ./ r[end] .- c # z <- -1 / (γ + d ⋅ y) G^ᵀ y - c
        ldiv!(cholQ, z) # z <- -Q⁻¹ (1 / (γ + d ⋅ y) G^ᵀ y + c)
        rmul!(λ, -1 / sqrt(-r[end])) # the sqrt appears to be missing in (12) in the paper
    else
        fill!(z, NaN)
        fill!(λ, NaN)
    end

    z, λ
end

"""
    primal_infeasibility(qp::QP, z)

Measure of primal infeasibility of the solution to the QP

Minimize ``\\frac{1}{2} z' Q z + c' z``
Subject to ``G z \\leq g`

where ``z`` is the primal solution and ``λ`` is the dual solution.

For a primal-feasible solution, this will return a value <= 0.
"""
primal_infeasibility(qp::QP, z, λ) = maximum(qp.G * z .- qp.g)

"""
    dual_infeasibility(qp::QP, λ)

Measure of primal infeasibility of the solution to the QP

Minimize ``\\frac{1}{2} z' Q z + c' z``
Subject to ``G z \\leq g`

where ``z`` is the primal solution and ``λ`` is the dual solution.

For a dual-feasible solution, this will return a value <= 0.
"""
dual_infeasibility(qp::QP, z, λ) = maximum(λ)

"""
    stationarity_violation(qp::QP, z, λ)

Measure of the violation of the stationarity KKT condition for the QP.

For an optimal solution, this should return a value near 0.
"""
function stationarity_violation(qp::QP, z, λ)
    slack_grad = qp.G' * λ
    cost_grad = qp.Q * z + qp.c
    scale = mean(cost_grad ./ slack_grad)
    maximum(abs, cost_grad .- scale .* slack_grad)
end

"""
    slackness_violation(qp::QP, z, λ)

Measure of the violation of the complementary slackness condition for the QP.

For an optimal solution, this should return a value near 0.
"""
slackness_violation(qp::QP, z, λ) = maximum(abs, (qp.G * z .- qp.g) .* λ)

"""
    check_optimality_conditions(qp::QP, z, λ)

Provies a scalar indicator of the violation of the KKT optimality conditions for
a given QP, primal solution ``z``, and dual solution ``λ``. A feasible & optimal
solution should return a value close to 0, while an infeasible or suboptimal
solution should return a value greater than zero.
"""
function check_optimality_conditions(qp::QP, z, λ)
    return max(
               primal_infeasibility(qp, z, λ),
               dual_infeasibility(qp, z, λ),
               stationarity_violation(qp, z, λ),
               slackness_violation(qp, z, λ)
               )
end

include("NNLSSolverInterface.jl")
using .NNLSSolverInterface: NNLSSolver

end # module

ext: jl


text:
export checkDerivative, getRandomTestDirection

using Printf

function getRandomTestDirection(x0::Float64)
	return randn()
end
function getRandomTestDirection(x0::ComplexF64)
	return randn()+1im*randn()
end
function getRandomTestDirection(x0::Array{Float64})
	return randn(size(x0))
end
function getRandomTestDirection(x0::Array{ComplexF64})
	return randn(size(x0))+1im*randn(size(x0))
end

function checkDerivative(f::Function,df::Function,x0;kwargs...)
	function testFun(x,v=[])
		if !(isempty(v))
			return f(x), df(x)*v
		else
			return f(x)
		end
	end
	return checkDerivative(testFun,x0;kwargs...)
end

function checkDerivative(f::Function,x0;out::Bool=true,tol::Float64=1.9,nSuccess::Int=3,v=getRandomTestDirection(x0))
# checkDerivative(f::Function,x0;out::Bool=true,tol::Float64=1.9,nSuccess::Int=3)
	if out
		println(@sprintf("%9s\t%9s\t%9s\t%9s\t%9s\t%5s","h","E0","E1","O1","O2","OK?"))
	end

	f0,dvf  = f(x0,v)
	nf0     = norm(f0)
	#dvf    = real(dvf)
	Error   = zeros(10,2)
	Order   = zeros(10,2)
	Success = zeros(10)
	for j=1:10
		ft = f(x0+10.0^(-j)*v)            # function value
		Error[j,1] = norm(f0-ft)/nf0          # Error TaylorPoly 0
		Error[j,2] = norm(f0 .+10.0^(-j)*dvf .- ft)/nf0 # Error TaylorPoly 1
		if j>1
			Order[j,:] = log10.(Error[j-1,:]./Error[j,:]);
		end
		if (Order[j,2]>tol) || (Error[j,1]/Error[j,2] > 100); Success[j]=1; end
		if out
			println(@sprintf("%1.3e\t%1.3e\t%1.3e\t%1.3e\t%1.3e\t%5d",
							10.0^(-j), Error[j,1],Error[j,2], Order[j,1],Order[j,2],Success[j]))
		end
	end
	pass = sum(Success) > nSuccess
	return  pass,Error,Order
end


text:
using CyclicDiagonal
using Test
using SpecialMatrices
using LinearAlgebra

@testset "CyclicDiagonal.jl" begin
    # Write your own tests here.

    M = Matrix( Circulant([1;3;2]) )

    # Test simple 3x3 circulant matric
    D1, C, D2 = CyclicDiagonal.iterDCD( M )
    @test norm(D1-I, Inf) < 1e-12

    # Test simple product of a diagonal and circulant matrix
    D = Diagonal(1:3)
    D1, C, D2 = CyclicDiagonal.iterDCD( D*M )
    @test norm( D1*C*D2 - D*M, Inf) < 1e-12

    # Test products of kxk matrices (k=1,...,K), each size N times
    N = 10
    K = 10
    r = 1:50
    for k in 1:K
        for i in 1:N
            D1o = Diagonal( rand(r, k) )
            Co  = Circulant( rand(r, k) )
            D2o = Diagonal( rand(r, k) )
            D1, C, D2 = CyclicDiagonal.iterDCD( D1o * Co * D2o, 1e6, 1e-7 )
            @test norm( D1*C*D2 - D1o*Co*D2o, Inf) < 1e-3
        end
    end
end


text:
module MatDeforNeohookeanNaiveModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtoolsDeforLinear.DeforModelRedModule: AbstractDeforModelRed, DeforModelRed3D, DeforModelRed2DStrain, DeforModelRed2DStress, DeforModelRed2DAxisymm, DeforModelRed1D, nstressstrain, nthermstrain
import FinEtoolsDeforLinear.MatDeforModule: AbstractMatDefor, stressvtot!, stressttov!
import ..MatDeforNonlinearModule: AbstractMatDeforNonlinear, estimatesoundspeed
import LinearAlgebra: Transpose, Diagonal, mul!
At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
A_mul_B!(C, A, B) = mul!(C, A, B)
import LinearAlgebra: eigen, eigvals, norm, cholesky, cross, dot, log, diagm, det

"""
	MatDeforNeohookeanNaive{MR<:AbstractDeforModelRed, MTAN<:Function, MUPD<:Function} <: AbstractMatDeforNonlinear

Type for triaxial neohookean hyperelastic material.

"""
struct  MatDeforNeohookeanNaive{MR<:AbstractDeforModelRed, MTAN<:Function, MUPD<:Function} <: AbstractMatDeforNonlinear
    mr::Type{MR} # model reduction type
    mass_density::FFlt # mass density
    E::FFlt # Young's modulus
    nu::FFlt # Poisson ratio
    _lambda::FFlt
    _mu::FFlt
    tangentmoduli!::MTAN # Function to return the tangent moduli matrix
    update!::MUPD # Function to update the material state
end

# function _threedD(E::FFlt, nu::FFlt)
# 	mI = Matrix(Diagonal([1.0, 1.0, 1.0, 0.5, 0.5, 0.5]))
# 	m1 = [1.0, 1.0, 1.0, 0.0, 0.0, 0.0];
# 	_lambda = E * nu / (1 + nu) / (1 - 2*(nu));
# 	_mu = E / 2. / (1+nu);
# 	D = _lambda * m1 * m1' + 2. * _mu * mI;
# 	return D
# end

################################################################################
# 3-D solid model
################################################################################

"""
    MatDeforNeohookeanNaive(mr::Type{DeforModelRed3D}, mass_density::FFlt, E::FFlt, nu::FFlt)

Create triaxial neohookean hyperelastic material.
"""
function MatDeforNeohookeanNaive(mr::Type{DeforModelRed3D}, mass_density::FFlt, E::FFlt, nu::FFlt)
    _mI = diagm(0=>[1, 1, 1, 0.5, 0.5, 0.5]);
    _m1 = vec(FFlt[1 1 1 0 0 0]) ;
    _m1m1 = _m1*_m1';
    _I3 = [1.0 0 0; 0 1.0 0; 0 0 1.0]
    function tangentmoduli3d!(self::MatDeforNeohookeanNaive, D::FFltMat, statev::FFltVec, Fn1::FFltMat, Fn::FFltMat, tn::FFlt, dtn::FFlt, loc::FFltMat, label::FInt)
		J = det(Fn1);
		copyto!(D, (self._lambda / J) * _m1m1 + 2 * (self._mu - self._lambda*log(J))/J * _mI);
		return D
    end
    function update3d!(self::MatDeforNeohookeanNaive, statev::FFltVec, cauchy::FFltVec, output::FFltVec, Fn1::FFltMat, Fn::FFltMat, tn::FFlt, dtn::FFlt, loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)
		@assert length(cauchy) == nstressstrain(self.mr)
		# Finger deformation tensor
		b = Fn1*Fn1';
		J = det(Fn1); # Jacobian of the deformation gradient
		sigma = (self._mu/J) .* (b .- _I3) .+ (self._lambda *log(J)/J) .* _I3;
		stressttov!(mr, cauchy, sigma)
		if quantity == :nothing
		    #Nothing to be copied to the output array
		elseif quantity == :cauchy || quantity == :Cauchy
		    (length(output) >= 6) || (output = zeros(6)) # make sure we can store it
		    copyto!(output, cauchy);
		elseif quantity == :pressure || quantity == :Pressure
		    output[1]  =  -sum(cauchy[1:3])/3.
		elseif quantity == :princCauchy || quantity == :princcauchy
		    ep = eigen(sigma);
		    (length(output) >= 3) || (output = zeros(3)) # make sure we can store it
		    copyto!(output,  sort(ep.values, rev=true));
		elseif quantity==:vonMises || quantity==:vonmises || quantity==:von_mises || quantity==:vm
		    s1=cauchy[1]; s2=cauchy[2]; s3=cauchy[3];
		    s4=cauchy[4]; s5=cauchy[5]; s6=cauchy[6];
		    (length(output) >= 1) || (output = zeros(1)) # make sure we can store it
		    output[1] = sqrt(1.0/2*((s1-s2)^2+(s1-s3)^2+(s2-s3)^2+6*(s4^2+s5^2+s6^2)))
		end
		return output
    end
    _lambda = E * nu / (1 + nu) / (1 - 2*(nu));
    _mu = E / 2. / (1+nu);
    return MatDeforNeohookeanNaive(mr, mass_density, E, nu, _lambda, _mu,
		              tangentmoduli3d!, update3d!)
end

"""
        MatDeforNeohookeanNaive(mr::Type{MR}, E::FFlt, nu::FFlt) where {MR}

    Create neohookean isotropic elastic material.

    The mass density is the default value.
"""
function MatDeforNeohookeanNaive(mr::Type{MR}, E::FFlt, nu::FFlt) where {MR}
    mass_density = 1.0
    return MatDeforNeohookeanNaive(mr, mass_density, E, nu)
end

"""
    estimatesoundspeed(self::M)  where {M<:AbstractMatDeforNonlinear}

Estimate sound speed in the undeformed state.
"""
function estimatesoundspeed(self::M)  where {M<:MatDeforNeohookeanNaive}
    return sqrt((self._lambda + 2 * self._mu) / self.mass_density)
end

end


text:
# these require decoding to work correctly

isinteger(a::Float65) = isinteger(reflect(a.fp))

(-){T<:Float65}(a::T)   = (T)( (-)(reflect(a.fp)) )
(abs){T<:Float65}(a::T) = (T)( (abs)(reflect(a.fp)) )

const sqrt65 = Dict( tiny => tiny, huge => huge )
function sqrt{T<:Float65}(a::T)
    get(sqrt65, a, sqrt65fn(a))
end
function sqrt65fn{T<:Float65}(a::T)
    (T)( sqrt(reflect(a.fp)) )
end

for op in (:(+), :(-), :(*), :(/), :(\), :(%), :(^),
           :div, :fld, :rem, :mod, :mod1, :rem1, :fld1)
    @eval begin
        ($op){T<:Float65}(a::T, b::Float65) = (T)( $op(reflect(a.fp), reflect(b.fp)) )
        ($op){T<:Float65}(a::T, b::Float64) = (T)( $op(reflect(a.fp), b) )
        ($op){T<:Float65}(a::Float64, b::T) = (T)( $op(a, reflect(b.fp)) )

        ($op){T<:Float65,I<:Union{Int64,Int32}}(a::T, b::I) = (T)( $op(reflect(a.fp), b) )
        ($op){T<:Float65,I<:Union{Int64,Int32}}(a::I, b::T) = (T)( $op(a, reflect(b.fp)) )
        ($op){T<:Float65}(a::T, b::Bool) = (T)( $op(reflect(a.fp), convert(Float64,b)) )
        ($op){T<:Float65}(a::Bool, b::T) = (T)( $op(convert(Float64,a), reflect(b.fp)) )
        ($op){T<:Float65,R<:Real}(a::T, b::R) = (T)( $op(reflect(a.fp), convert(Float64,b)) )
        ($op){T<:Float65,R<:Real}(a::R, b::T) = (T)( $op(convert(Float64,a), reflect(b.fp)) )
    end
end

eps{T<:Float65}(a::T) = (T)( eps(reflect(a.fp))  )


text:
function __new__(T, args...)
    # @show T
    # @show args
    # note: we also add __new__() to the list of primitives so it's not overdubbed recursively
    if T <: NamedTuple
        return T(args)
    else
        return T(args...)
    end
end

__tuple__(args...) = tuple(args...)
__getfield__(args...) = getfield(args...)

function module_functions(modl)
    res = Vector{Function}()
    for s in Base.names(modl; all=true)
        isdefined(modl, s) || continue
        fn = getfield(modl, s)
        if fn isa Function # && match(r"^[a-z#]+$", string(s)) != nothing
            push!(res, fn)
        end
    end
    return res
end

const PRIMITIVES = Set{Any}(vcat(
    module_functions(Base),
    module_functions(Core),
    module_functions(Core.Intrinsics),
    [Broadcast.materialize, Broadcast.broadcasted, Colon(), (:),
     Base.not_int,
     # our own special functions
     __new__, __tuple__, __getfield__, namedtuple, guess_device]));

include("cassette.jl")
# include("interp.jl")
include("irtools.jl")

trace = irtrace


text:
using Metis
using Random
using Test
using SparseArrays
import TriangleMesh
import LightGraphs

@testset "Metis.permutation" begin
    Random.seed!(0)
    S = sprand(10, 10, 0.5); S = S + S'; fill(S.nzval, 1)
    perm, iperm = Metis.permutation(S)
    @test isperm(perm)
    @test isperm(iperm)
    @test S == S[perm,perm][iperm,iperm]
end

@testset "Metis.partition" begin
    Random.seed!(0)
    S = sprand(10, 10, 0.5); S = S + S'; fill(S.nzval, 1)
    T = LightGraphs.smallgraph(:tutte)
    for G in (S, T), alg in (:RECURSIVE, :KWAY), nparts in (3, 4)
        partition = Metis.partition(G, nparts, alg = alg)
        @test extrema(partition) == (1, nparts)
        @test all(x -> findfirst(==(x), partition) !== nothing, 1:nparts)
    end
end

@testset "Metis.separator" begin
    Random.seed!(0)
    S = sprand(10, 10, 0.5); S = S + S'; fill(S.nzval, 1)
    T = LightGraphs.smallgraph(:tutte)
    for G in (S, T)
        parts = Metis.separator(G)
        @test extrema(parts) == (1, 3)
    end
end

@testset "Metis.mesh_partition" begin
    poly = TriangleMesh.polygon_Lshape()
    L_mesh = TriangleMesh.create_mesh(poly, info_str="L_mesh", voronoi=true, delaunay=true)
    L_mesh = TriangleMesh.refine(L_mesh, divide_cell_into=1_000, keep_edges=true)
    poly = TriangleMesh.polygon_unitSquareWithHole()
    S_mesh = TriangleMesh.create_mesh(poly, info_str="S_mesh", voronoi=true, delaunay=true)
    S_mesh = TriangleMesh.refine(S_mesh, divide_cell_into=1_000, keep_edges=true)
    for mesh in (L_mesh, S_mesh), alg in (:DUAL, :NODAL), nparts in (20, 50)
        epart, npart = Metis.mesh_partition(mesh, nparts, alg = alg)
        @test extrema(epart) == (1, nparts)
        @test extrema(npart) == (1, nparts)
        @test all(x -> findfirst(==(x), epart) !== nothing, 1:nparts)
        @test all(x -> findfirst(==(x), npart) !== nothing, 1:nparts)
    end
end


text:
###############################################################
#                Combined Mean and Variance
###############################################################
"""
    combined_mean_and_var(xs...) -> meanc, varc
Calculates the combined mean and variance of the concatenated sample `vcat(xs...)`.
"""
function combined_mean_and_var(xstuple::AbstractVector{T}...) where T
    xs = collect(xstuple)
    ns = length.(xs)
    μs = mean.(xs)
    vs = var.(xs)
    return combined_mean_and_var(ns, μs, vs)
end

export combined_mean_and_var

"""
    combined_mean_and_var(ns, μs, vs) -> meanc, varc
Given N samples characterized by their lengths `ns`,
their means `μs`, and their variances `vs`,
calculates the combined (or pooled) mean and variance of the
overall (concatenated) sample.
"""
function combined_mean_and_var(ns::AbstractVector{<:Integer},
                               μs::AbstractVector{<:Number},
                               vs::AbstractVector{<:Number})
    nsum = sum(ns)
    meanc = dot(ns, μs) / nsum
    varc = sum((ns .- 1) .* vs + ns .* abs2.(μs .- meanc)) / (nsum - 1)
    return meanc, varc
end

function combined_mean_and_var(ns::AbstractVector{<:Integer},
                               μs::AbstractVector{<:AbstractArray{<:Number}},
                               vs::AbstractVector{<:AbstractArray{<:Number}})
    meanc = zero(μs[1])
    varc = zero(vs[1])

    nsum = sum(ns)
    N = length(μs) # number of samples
    for i in eachindex(meanc)

        for k in 1:N
            meanc[i] += ns[k] * μs[k][i]
        end
        meanc[i] = meanc[i] / nsum

        for k in 1:N
            varc[i] += (ns[k] - 1) * vs[k][i] + ns[k] * abs2(μs[k][i] - meanc[i])
        end
        varc[i] = varc[i] / (nsum - 1)
    end

    return meanc, varc
end

###############################################################
#       Combined Mean and Variance ("Hand-written")
###############################################################
"""
    combined_mean_and_var(x1, x2) -> meanc, varc
Given two samples `x1`,`x2` calculates the mean and variance of the
concatenated sample.
"""
function combined_mean_and_var(x1::AbstractVector{<:Number}, x2::AbstractVector{<:Number})
    n1, n2 = length(x1), length(x2)
    μ1, μ2 = mean(x1), mean(x2)
    v1, v2 = var(x1), var(x2)
    return combined_mean_and_var(n1, μ1, v1, n2, μ2, v2)
end

"""
    combined_mean_and_var(n1, μ1, v1, n2, μ2, v2) -> meanc, varc
Given two samples characterized by their lengths `n1`, `n2`,
their means `μ1`, `μ2`, and their variances `v1`, `v2`,
calculates the combined (or pooled) mean and variance of the
concatenated sample.
"""
function combined_mean_and_var(n1::Integer, μ1::Number, v1::Number,
                               n2::Integer, μ2::Number, v2::Number)
    meanc = (n1 * μ1 + n2 * μ2) / (n1 + n2)

    # Based on https://www.emathzone.com/tutorials/basic-statistics/combined-variance.html,
    # including Robert Matheson's comment and adding abs for complex number support.
    varc = ((n1-1)*v1 + (n2-1)*v2 + n1*abs2(μ1 - meanc) + n2*abs2(μ2 - meanc)) /
                                (n1 + n2 - 1)
    return meanc, varc
end

# Explicit version for three samples (only for testing)
function combined_mean_and_var_three(x1, x2, x3)
    n1, n2, n3 = length(x1), length(x2), length(x3)
    μ1, μ2, μ3 = mean(x1), mean(x2), mean(x3)
    v1, v2, v3 = var(x1), var(x2), var(x3)

    meanc12, varc12 = combined_mean_and_var(n1, μ1, v1, n2, μ2, v2)
    n12 = n1 + n2
    return combined_mean_and_var(n12, meanc12, varc12, n3, μ3, v3)
end

"""
    combined_mean_and_var(ors::Vector{<:ObservableResult}) -> meanc, stderrc
Calculates the combined mean and standard error of the concatenated observable results.
"""
function combined_mean_and_error(ors::Vector{<:ObservableResult})
    @assert length(ors) > 0 "Provided vector of ObservableResults must not be empty."

    ns = [r.count for r in ors]
    μs = [r.mean for r in ors]
    vs = [r.error^2 for r in ors]

    meanc, varc = combined_mean_and_var(ns, μs, vs)
    return meanc, sqrt.(varc)
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
@testset "Check center and scaling methods" begin

    @testset "Check center" begin

        @test IndexFunArrays.get_offset((5,), CtrCorner) == (1.0,)
        @test IndexFunArrays.get_offset((4,), CtrCorner) == (1.0,)
        @test IndexFunArrays.get_offset((5, 2), CtrFT) == (3.0, 2.0)
        @test IndexFunArrays.get_offset((5, 1), CtrFT) == (3.0, 1.0)
        @test IndexFunArrays.get_offset((5, 2), CtrFFT) == (1.0, 1.0)

        @test IndexFunArrays.get_offset((5, 2), CtrMid) == (3.0, 1.5)
        @test IndexFunArrays.get_offset((5, 3), CtrMid) == (3.0, 2.0)

        @test IndexFunArrays.get_offset((5, 2), CtrRFT) == (1, 2)
        @test IndexFunArrays.get_offset((5, 3), CtrRFT) == (1, 2)

        @test IndexFunArrays.get_offset((5, 2), CtrRFFT) == (1, 1)
        @test IndexFunArrays.get_offset((5, 3), CtrRFFT) == (1, 1)

        @test IndexFunArrays.get_offset((1123, 3123, 2), 123.123) == (123.123, 123.123, 123.123)

        @test IndexFunArrays.get_offset((5, 4), CtrEnd) == (5.0, 4.0)
    end

    @testset "Check scaling" begin
        @test IndexFunArrays.get_scale((5, 4), ScaFT) === (0.2, 0.25)
        @test IndexFunArrays.get_scale((5, 4), ScaUnit) == (1, 1)
        @test IndexFunArrays.get_scale((6, 5), ScaNorm) == (0.2, 0.25)
        @test IndexFunArrays.get_scale((6, 5), ScaNorm) == (0.2, 0.25)
        @test IndexFunArrays.get_scale((6, 5), ScaMid) == (0.4, 0.5)
        @test IndexFunArrays.get_scale((5,3), (2,11)) == (2,11)

        @test IndexFunArrays.get_scale((12, 1, 2), 2.123f0) == (2.123f0, 2.123f0, 2.123f0)

    end
end


text:
"""
    AtomicSpecies(atom::Union{AbstractChar,String}, mass::Float64, pseudopot::String)
    AtomicSpecies(x::AtomicPosition, mass, pseudopot)

Represent each line of the `ATOMIC_SPECIES` card in QE.

The `atom` field accepts at most 3 characters.

# Examples
```julia
julia> using QuantumESPRESSOBase.Cards.PWscf

julia> AtomicSpecies("C1", 12, "C.pbe-n-kjpaw_psl.1.0.0.UPF")
AtomicSpecies("C1", 12.0, "C.pbe-n-kjpaw_psl.1.0.0.UPF")

julia> AtomicSpecies(
           AtomicPosition('S', [0.500000000, 0.288675130, 1.974192764]),
           32.066,
           "S.pz-n-rrkjus_psl.0.1.UPF",
       )
AtomicSpecies("S", 32.066, "S.pz-n-rrkjus_psl.0.1.UPF")
```
"""
struct AtomicSpecies
    "Label of the atom. Max total length cannot exceed 3 characters."
    atom::String
    """
    Mass of the atomic species in atomic unit.

    Used only when performing molecular dynamics (MD) run
    or structural optimization runs using damped MD.
    Not actually used in all other cases (but stored
    in data files, so phonon calculations will use
    these values unless other values are provided).
    """
    mass::Float64
    """
    File containing pseudopotential for this species.

    See also: [`pseudoformat`](@ref)
    """
    pseudopot::String
    function AtomicSpecies(atom::Union{AbstractChar,AbstractString}, mass, pseudopot)
        @assert(length(atom) <= 3, "Max total length of `atom` cannot exceed 3 characters!")
        return new(string(atom), mass, pseudopot)
    end
end

"""
    pseudoformat(data::AtomicSpecies)

Return the pseudopotential format of the `AtomicSpecies`.
"""
# pseudoformat(data::AtomicSpecies) = pseudoformat(data.pseudopot)

"""
    AtomicSpeciesCard <: Card

Represent the `ATOMIC_SPECIES` card in QE. It does not have an "option".
"""
struct AtomicSpeciesCard <: Card
    data::Vector{AtomicSpecies}
end
AtomicSpeciesCard(cell::Cell) = AtomicSpeciesCard(map(AtomicSpecies ∘ string, cell.numbers))

"""
    AtomicPosition(atom::Union{AbstractChar,String}, pos::Vector{Float64}[, if_pos::Vector{Int}])
    AtomicPosition(x::AtomicSpecies, pos, if_pos)

Represent each line of the `ATOMIC_POSITIONS` card in QE.

The `atom` field accepts at most 3 characters.

# Examples
```julia
julia> using QuantumESPRESSOBase.Cards.PWscf

julia> AtomicPosition('O', [0, 0, 0])
AtomicPosition("O", [0.0, 0.0, 0.0], Bool[1, 1, 1])

julia> AtomicPosition(
           AtomicSpecies('S', 32.066, "S.pz-n-rrkjus_psl.0.1.UPF"),
           [0.500000000, 0.288675130, 1.974192764],
       )
AtomicPosition("S", [0.5, 0.28867513, 1.974192764], Bool[1, 1, 1])
```
"""
struct AtomicPosition
    "Label of the atom as specified in `AtomicSpecies`."
    atom::String
    "Atomic positions. A three-element vector of floats."
    pos::SVector{3,Float64}
    """
    Component `i` of the force for this atom is multiplied by `if_pos(i)`,
    which must be either `0` or `1`.  Used to keep selected atoms and/or
    selected components fixed in MD dynamics or structural optimization run.

    With `crystal_sg` atomic coordinates the constraints are copied in all equivalent
    atoms.
    """
    if_pos::SVector{3,Bool}
    function AtomicPosition(atom::Union{AbstractChar,AbstractString}, pos, if_pos)
        @assert(length(atom) <= 3, "the max length of `atom` cannot exceed 3 characters!")
        return new(string(atom), pos, if_pos)
    end
end
AtomicPosition(atom, pos) = AtomicPosition(atom, pos, trues(3))
AtomicPosition(x::AtomicSpecies, pos, if_pos) = AtomicPosition(x.atom, pos, if_pos)
# Introudce mutual constructors since they share the same atoms.
AtomicSpecies(x::AtomicPosition, mass, pseudopot) = AtomicSpecies(x.atom, mass, pseudopot)

"""
    AtomicPositionsCard <: Card

Represent the `ATOMIC_POSITIONS` card in QE.

# Arguments
- `data::AbstractVector{AtomicPosition}`: A vector containing `AtomicPosition`s.
- `option::String="alat"`: allowed values are: "alat", "bohr", "angstrom", "crystal", and "crystal_sg".
"""
struct AtomicPositionsCard <: Card
    data::Vector{AtomicPosition}
    option::String
    function AtomicPositionsCard(data, option = "alat")
        @assert option in optionpool(AtomicPositionsCard)
        return new(data, option)
    end
end
AtomicPositionsCard(cell::Cell, option) = AtomicPositionsCard(
    [
        AtomicPosition(string(atom), pos)
        for (atom, pos) in zip(cell.numbers, cell.positions)
    ],
    option,
)
# Introudce mutual constructors since they share the same atoms.

function validate(x::AtomicSpeciesCard, y::AtomicPositionsCard)
    lens = @lens _.data.atom
    @assert(
        isempty(symdiff(map(Base.Fix2(get, lens) ∘ unique, (x, y)))),
        "labels of the atoms are different in `ATOMIC_SPECIES` and `ATOMIC_POSITIONS` card!",
    )
end # function validate
validate(y::AtomicPositionsCard, x::AtomicSpeciesCard) = validate(x, y)

"Represent the abstraction of `CELL_PARAMETERS` and `REF_CELL_PARAMETERS` cards in QE."
abstract type AbstractCellParametersCard <: Card end

"""
    CellParametersCard{T<:Real} <: AbstractCellParametersCard
    CellParametersCard(data::AbstractMatrix, option::String)

Represent the `CELL_PARAMETERS` cards in `PWscf` and `CP` packages.
"""
struct CellParametersCard{T<:Real} <: AbstractCellParametersCard
    data::SMatrix{3,3,T}
    option::String
    function CellParametersCard{T}(data, option = "alat") where {T<:Real}
        @assert option in optionpool(CellParametersCard)
        return new(data, option)
    end
end
CellParametersCard(data::AbstractMatrix{T}, option = "alat") where {T} =
    CellParametersCard{T}(data, option)
CellParametersCard(lattice::Lattice{T}, option) where {T} =
    CellParametersCard(convert(Matrix{T}, lattice), option)
CellParametersCard(cell::Cell, option) = CellParametersCard(cell.lattice, option)

struct AtomicForce
    atom::String
    force::SVector{3,Float64}
    function AtomicForce(atom::Union{AbstractChar,AbstractString}, force)
        @assert(length(atom) <= 3, "the max length of `atom` cannot exceed 3 characters!")
        return new(string(atom), force)
    end
end

struct AtomicForcesCard <: Card
    data::Vector{AtomicForce}
end

"""
    AtomicVelocity(atom::Union{AbstractChar,String}, velocity::Vector{Float64})
    AtomicVelocity(x::AtomicPosition, velocity)

Represent each line of the `ATOMIC_VELOCITIES` card in QE's `CP` package.

The `atom` field accepts at most 3 characters.

# Examples
```julia
julia> using QuantumESPRESSOBase.Cards.CP

julia> AtomicVelocity("H", [0.140374e-04, -0.333683e-04, 0.231834e-04])
AtomicVelocity("H", [1.40374e-5, -3.33683e-5, 2.31834e-5])
```
"""
struct AtomicVelocity
    atom::String
    velocity::SVector{3,Float64}
    function AtomicVelocity(atom::Union{AbstractChar,AbstractString}, velocity)
        @assert(length(atom) <= 3, "the max length of `atom` cannot exceed 3 characters!")
        return new(string(atom), velocity)
    end
end
AtomicVelocity(x::AtomicSpecies, velocity) = AtomicVelocity(x.atom, velocity)
AtomicVelocity(x::AtomicPosition, velocity) = AtomicVelocity(x.atom, velocity)
# Introudce mutual constructors since they share the same atoms.
"""
    AtomicSpecies(x::AtomicVelocity, mass, pseudopot)

Construct an `AtomicSpecies` from an `AtomicVelocity` instance.
"""
AtomicSpecies(x::AtomicVelocity, mass, pseudopot) = AtomicSpecies(x.atom, mass, pseudopot)
"""
    AtomicPosition(x::AtomicVelocity, pos, if_pos)

Construct an `AtomicPosition` from an `AtomicVelocity` instance.
"""
AtomicPosition(x::AtomicVelocity, pos, if_pos) = AtomicPosition(x.atom, pos, if_pos)

"""
    AtomicVelocitiesCard <: Card

Represent the `ATOMIC_VELOCITIES` card in QE's `CP` package which does not have an "option".
"""
struct AtomicVelocitiesCard <: Card
    data::Vector{AtomicVelocity}
end

struct RefCellParametersCard <: AbstractCellParametersCard
    data::SMatrix{3,3,Float64}
    option::String
    function RefCellParametersCard(data, option = "bohr")
        @assert option in optionpool(RefCellParametersCard)
        return new(data, option)
    end
end
RefCellParametersCard(data::AbstractMatrix, option = "bohr") =
    RefCellParametersCard(data, option)


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot(x=[0.236033, 0.346517, 0.312707, 0.00790928, 0.488613, 0.210968, 0.951916, 0.999905, 0.251662, 0.986666],
     y=[0.555751, 0.437108, 0.424718, 0.773223, 0.28119, 0.209472, 0.251379, 0.0203749, 0.287702, 0.859512],
     color=[1.07999, 1.89723, 2.39539, 1.32126, 2.11978, 1.90576, 1.08094, 2.33541, 1.08941, 1.73882],
     Scale.color_log(minvalue=1))


text:
# Defines FaceScalarValues and FaceVectorValues and common methods
"""
    FaceScalarValues([::Type{T}], quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation])
    FaceVectorValues([::Type{T}], quad_rule::QuadratureRule, func_interpol::Interpolation, [geom_interpol::Interpolation])

A `FaceValues` object facilitates the process of evaluating values of shape functions, gradients of shape functions,
values of nodal functions, gradients and divergences of nodal functions etc. on the faces of finite elements. There are
two different types of `FaceValues`: `FaceScalarValues` and `FaceVectorValues`. As the names suggest,
`FaceScalarValues` utilizes scalar shape functions and `FaceVectorValues` utilizes vectorial shape functions.
For a scalar field, the `FaceScalarValues` type should be used. For vector field, both subtypes can be used.

!!! note
    The quadrature rule for the face should be given with one dimension lower.
    I.e. for a 3D case, the quadrature rule should be in 2D.

**Arguments:**

* `T`: an optional argument to determine the type the internal data is stored as.
* `quad_rule`: an instance of a [`QuadratureRule`](@ref)
* `func_interpol`: an instance of an [`Interpolation`](@ref) used to interpolate the approximated function
* `geom_interpol`: an optional instance of an [`Interpolation`](@ref) which is used to interpolate the geometry

**Common methods:**

* [`reinit!`](@ref)
* [`getnquadpoints`](@ref)
* [`getdetJdV`](@ref)

* [`shape_value`](@ref)
* [`shape_gradient`](@ref)
* [`shape_symmetric_gradient`](@ref)
* [`shape_divergence`](@ref)

* [`function_value`](@ref)
* [`function_gradient`](@ref)
* [`function_symmetric_gradient`](@ref)
* [`function_divergence`](@ref)
* [`spatial_coordinate`](@ref)
"""
FaceValues

# FaceScalarValues
struct FaceScalarValues{dim,T<:Real,refshape<:AbstractRefShape} <: FaceValues{dim,T,refshape}
    N::Array{T,3}
    dNdx::Array{Vec{dim,T},3}
    dNdξ::Array{Vec{dim,T},3}
    detJdV::Matrix{T}
    normals::Vector{Vec{dim,T}}
    M::Array{T,3}
    dMdξ::Array{Vec{dim,T},3}
    qr_weights::Vector{T}
    current_face::ScalarWrapper{Int}
end

function FaceScalarValues(quad_rule::QuadratureRule, func_interpol::Interpolation,
                          geom_interpol::Interpolation=func_interpol)
    FaceScalarValues(Float64, quad_rule, func_interpol, geom_interpol)
end

function FaceScalarValues(::Type{T}, quad_rule::QuadratureRule{dim_qr,shape}, func_interpol::Interpolation,
        geom_interpol::Interpolation=func_interpol) where {dim_qr,T,shape<:AbstractRefShape}

    @assert getdim(func_interpol) == getdim(geom_interpol)
    @assert getrefshape(func_interpol) == getrefshape(geom_interpol) == shape
    n_qpoints = length(getweights(quad_rule))
    dim = dim_qr + 1

    face_quad_rule = create_face_quad_rule(quad_rule, func_interpol)
    n_faces = length(face_quad_rule)

    # Normals
    normals = zeros(Vec{dim,T}, n_qpoints)

    # Function interpolation
    n_func_basefuncs = getnbasefunctions(func_interpol)
    N =    fill(zero(T)          * T(NaN), n_func_basefuncs, n_qpoints, n_faces)
    dNdx = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints, n_faces)
    dNdξ = fill(zero(Vec{dim,T}) * T(NaN), n_func_basefuncs, n_qpoints, n_faces)

    # Geometry interpolation
    n_geom_basefuncs = getnbasefunctions(geom_interpol)
    M =    fill(zero(T)          * T(NaN), n_geom_basefuncs, n_qpoints, n_faces)
    dMdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_qpoints, n_faces)

    for face in 1:n_faces, (qp, ξ) in enumerate(face_quad_rule[face].points)
        for i in 1:n_func_basefuncs
            dNdξ[i, qp, face], N[i, qp, face] = gradient(ξ -> value(func_interpol, i, ξ), ξ, :all)
        end
        for i in 1:n_geom_basefuncs
            dMdξ[i, qp, face], M[i, qp, face] = gradient(ξ -> value(geom_interpol, i, ξ), ξ, :all)
        end
    end

    detJdV = fill(T(NaN), n_qpoints, n_faces)

    FaceScalarValues{dim,T,shape}(N, dNdx, dNdξ, detJdV, normals, M, dMdξ, quad_rule.weights, ScalarWrapper(0))
end

# FaceVectorValues
struct FaceVectorValues{dim,T<:Real,refshape<:AbstractRefShape,M} <: FaceValues{dim,T,refshape}
    N::Array{Vec{dim,T},3}
    dNdx::Array{Tensor{2,dim,T,M},3}
    dNdξ::Array{Tensor{2,dim,T,M},3}
    detJdV::Matrix{T}
    normals::Vector{Vec{dim,T}}
    M::Array{T,3}
    dMdξ::Array{Vec{dim,T},3}
    qr_weights::Vector{T}
    current_face::ScalarWrapper{Int}
end

function FaceVectorValues(quad_rule::QuadratureRule, func_interpol::Interpolation, geom_interpol::Interpolation=func_interpol)
    FaceVectorValues(Float64, quad_rule, func_interpol, geom_interpol)
end

function FaceVectorValues(::Type{T}, quad_rule::QuadratureRule{dim_qr,shape}, func_interpol::Interpolation,
        geom_interpol::Interpolation=func_interpol) where {dim_qr,T,shape<:AbstractRefShape}

    @assert getdim(func_interpol) == getdim(geom_interpol)
    @assert getrefshape(func_interpol) == getrefshape(geom_interpol) == shape
    n_qpoints = length(getweights(quad_rule))
    dim = dim_qr + 1

    face_quad_rule = create_face_quad_rule(quad_rule, func_interpol)
    n_faces = length(face_quad_rule)

    # Normals
    normals = zeros(Vec{dim,T}, n_qpoints)

    # Function interpolation
    n_func_basefuncs = getnbasefunctions(func_interpol) * dim
    N    = fill(zero(Vec{dim,T})      * T(NaN), n_func_basefuncs, n_qpoints, n_faces)
    dNdx = fill(zero(Tensor{2,dim,T}) * T(NaN), n_func_basefuncs, n_qpoints, n_faces)
    dNdξ = fill(zero(Tensor{2,dim,T}) * T(NaN), n_func_basefuncs, n_qpoints, n_faces)

    # Geometry interpolation
    n_geom_basefuncs = getnbasefunctions(geom_interpol)
    M    = fill(zero(T)          * T(NaN), n_geom_basefuncs, n_qpoints, n_faces)
    dMdξ = fill(zero(Vec{dim,T}) * T(NaN), n_geom_basefuncs, n_qpoints, n_faces)

    for face in 1:n_faces, (qp, ξ) in enumerate(face_quad_rule[face].points)
        basefunc_count = 1
        for basefunc in 1:getnbasefunctions(func_interpol)
            dNdξ_temp, N_temp = gradient(ξ -> value(func_interpol, basefunc, ξ), ξ, :all)
            for comp in 1:dim
                N_comp = zeros(T, dim)
                N_comp[comp] = N_temp
                N[basefunc_count, qp, face] = Vec{dim,T}((N_comp...,))

                dN_comp = zeros(T, dim, dim)
                dN_comp[comp, :] = dNdξ_temp
                dNdξ[basefunc_count, qp, face] = Tensor{2,dim,T}((dN_comp...,))
                basefunc_count += 1
            end
        end
        for basefunc in 1:n_geom_basefuncs
            dMdξ[basefunc, qp, face], M[basefunc, qp, face] = gradient(ξ -> value(geom_interpol, basefunc, ξ), ξ, :all)
        end
    end

    detJdV = fill(T(NaN), n_qpoints, n_faces)
    MM = Tensors.n_components(Tensors.get_base(eltype(dNdx)))

    FaceVectorValues{dim,T,shape,MM}(N, dNdx, dNdξ, detJdV, normals, M, dMdξ, quad_rule.weights, ScalarWrapper(0))
end

function reinit!(fv::FaceValues{dim}, x::AbstractVector{Vec{dim,T}}, face::Int) where {dim,T}
    n_geom_basefuncs = getngeobasefunctions(fv)
    n_func_basefuncs = getnbasefunctions(fv)
    @assert length(x) == n_geom_basefuncs
    @boundscheck checkface(fv, face)

    fv.current_face[] = face
    cb = getcurrentface(fv)

    @inbounds for i in 1:length(fv.qr_weights)
        w = fv.qr_weights[i]
        fefv_J = zero(Tensor{2,dim})
        for j in 1:n_geom_basefuncs
            fefv_J += x[j] ⊗ fv.dMdξ[j, i, cb]
        end
        weight_norm = weighted_normal(fefv_J, fv, cb)
        fv.normals[i] = weight_norm / norm(weight_norm)
        detJ = norm(weight_norm)

        detJ > 0.0 || throw(ArgumentError("det(J) is not positive: det(J) = $(detJ)"))
        fv.detJdV[i, cb] = detJ * w
        Jinv = inv(fefv_J)
        for j in 1:n_func_basefuncs
            fv.dNdx[j, i, cb] = fv.dNdξ[j, i, cb] ⋅ Jinv
        end
    end
end

"""
    getcurrentface(fv::FaceValues)

Return the current active face of the `FaceValues` object (from last `reinit!`).

"""
getcurrentface(fv::FaceValues) = fv.current_face[]

"""
    getnormal(fv::FaceValues, qp::Int)

Return the normal at the quadrature point `qp` for the active face of the
`FaceValues` object(from last `reinit!`).
"""
getnormal(fv::FaceValues, qp::Int) = fv.normals[qp]

"""
    BCValues(func_interpol::Interpolation, geom_interpol::Interpolation, boundary_type::Union{Type{<:BoundaryIndex}})

`BCValues` stores the shape values at all faces/edges/vertices (depending on `boundary_type`) for the geomatric interpolation (`geom_interpol`),
for each dof-position determined by the `func_interpol`. Used mainly by the `ConstrainHandler`.
"""
struct BCValues{T}
    M::Array{T,3}
    current_face::ScalarWrapper{Int}
end

BCValues(func_interpol::Interpolation, geom_interpol::Interpolation, boundary_type::Type{<:BoundaryIndex} = Ferrite.FaceIndex) =
    BCValues(Float64, func_interpol, geom_interpol, boundary_type)

function BCValues(::Type{T}, func_interpol::Interpolation{dim,refshape}, geom_interpol::Interpolation{dim,refshape}, boundary_type::Type{<:BoundaryIndex} = Ferrite.FaceIndex) where {T,dim,refshape}
    # set up quadrature rules for each face with dof-positions
    # (determined by func_interpol) as the quadrature points
    interpolation_coords = reference_coordinates(func_interpol)

    qrs = QuadratureRule{dim,refshape,T}[]
    faces = boundaryfunction(boundary_type) # faces, edges or vertices
    for face in faces(func_interpol)
        dofcoords = Vec{dim,T}[]
        for facedof in face
            push!(dofcoords, interpolation_coords[facedof])
        end
        qrf = QuadratureRule{dim,refshape,T}(fill(T(NaN), length(dofcoords)), dofcoords) # weights will not be used
        push!(qrs, qrf)
    end

    n_faces = length(qrs)
    n_qpoints = n_faces == 0 ? 0 : length(getweights(qrs[1])) # assume same in all
    n_geom_basefuncs = getnbasefunctions(geom_interpol)
    M =    fill(zero(T)           * T(NaN), n_geom_basefuncs, n_qpoints, n_faces)

    for face in 1:n_faces, (qp, ξ) in enumerate(qrs[face].points)
        for i in 1:n_geom_basefuncs
            M[i, qp, face] = value(geom_interpol, i, ξ)
        end
    end

    BCValues{T}(M, ScalarWrapper(0))
end

getnquadpoints(bcv::BCValues) = size(bcv.M, 2)
function spatial_coordinate(bcv::BCValues, q_point::Int, xh::AbstractVector{Vec{dim,T}}) where {dim,T}
    n_base_funcs = size(bcv.M, 1)
    @assert length(xh) == n_base_funcs
    x = zero(Vec{dim,T})
    face = bcv.current_face[]
    @inbounds for i in 1:n_base_funcs
        x += bcv.M[i,q_point,face] * xh[i] # geometric_value(fe_v, q_point, i) * xh[i]
    end
    return x
end

nfaces(fv) = size(fv.N, 3)

function checkface(fv::FaceValues, face::Int)
    0 < face <= nfaces(fv) || error("Face index out of range.")
    return nothing
end


text:
using Copulas
using Test
using Distributions
using Random

@testset "Copulas.jl" begin
    MyD = SklarDist(ClaytonCopula(3,7),(LogNormal(),Pareto(),Beta()))
    u = rand(MyD,10000)
    rand!(MyD,u)
    fit(SklarDist{ClaytonCopula,Tuple{LogNormal,Pareto,Beta}},u)
    fit(SklarDist{GaussianCopula,Tuple{LogNormal,Pareto,Beta}},u)
end


text:
using Test

import AlphaZero: Util

n = 1000
num_workers = 100
res = Util.mapreduce(1:n, num_workers, +, 0) do
  process(i) = i^2
  terminate() = nothing
  return (process=process, terminate=terminate)
end
@test res == sum(i^2 for i in 1:n)


text:
export LinearMassSpringChain
export LinearMassSpringDamperChain

mutable struct LinearMassSpringDamperChain <: ControlBenchmarkProblem
    n::Integer          # The number of masses in the chain
    k::Vector{Float64}  # The spring constants in the chain
    β::Vector{Float64}  # The damping constants in the chain
    m::Vector{Float64}  # The masses in the chain
    inputs::Vector{Int} # The masses which have an acceleration input

    function LinearMassSpringDamperChain( n::Integer, k, β, m; inputs = [1; n] )
        if( n < 1 )
            throw( DomainError( n, "Must specify more than one mass" ) )
        end

        if isa( k, Number )
            k = fill( k, n+1 )
        end

        if isa( β, Number )
            β = fill( β, n+1 )
        end

        if isa( m, Number )
            m = fill( m, n )
        end

        # Create vectors
        k = reshape( k, : )
        β = reshape( β, : )
        m = reshape( m, : )

        if size( k, 1 ) != n+1
            error( "k must have n+1 values" )
        end

        if size( β, 1 ) != n+1
            error( "β must have n+1 values" )
        end

        if size( m, 1 ) != n
            error( "m must have n values" )
        end

        # Expand ranges and ensure no inputs are duplicated
        inArray = collect( inputs )
        inArray = reshape( inArray, : )
        unique!( inArray )

        if size( inArray, 1 ) > n
            error( inArray, "More inputs requested than masses" )
        end

        if !all( inArray .<= n )
            throw( DomainError( inputs, "Inputs must be on masses" ) )
        end

        return new(n, k, β, m, inArray)
    end
end

LinearMassSpringChain( n::Integer, k, m; inputs = [1; n] ) = LinearMassSpringDamperChain( n, k, 0, m; inputs )

function controlbenchmark( chain::LinearMassSpringDamperChain )
    # Extract into temp variables for more easy reference
    n = chain.n
    k = chain.k
    β = chain.β
    m = chain.m
    inputs = chain.inputs

    # Create the diagonal elements for the matrix part for the springs
    springUpperDiag = k[2:end-1] ./ m[1:end-1]
    springMainDiag  = -( k[1:end-1] + k[2:end] ) ./ m[1:end]
    springLowerDiag = k[2:end-1] ./ m[2:end]

    # Create the state matrix
    S = spdiagm(  0 => springMainDiag,    # The main diagonal matrix elements
                 -1 => springLowerDiag,   # The elements below the diagonal
                  1 => springUpperDiag )  # The elements above the diagonal

    # Create the diagonal elements for the matrix part for the damper
    damperUpperDiag = β[2:end-1] ./ m[1:end-1]
    damperMainDiag  = -( β[1:end-1] + β[2:end] ) ./ m[1:end]
    damperLowerDiag = β[2:end-1] ./ m[2:end]

    Beta = spdiagm(  0 => damperMainDiag,    # The main diagonal matrix elements
                    -1 => damperLowerDiag,   # The elements below the diagonal
                     1 => damperUpperDiag )  # The elements above the diagonal

    Z   = spzeros( n, n )
    Eye = sparse( 1.0*I, n, n )

    A = [ Z Eye
          S Beta ]

    # Create the input matrix
    dim = size( inputs, 1 )
    ind = CartesianIndex.( inputs .+ n, 1:dim )
    B   = spzeros( 2*n, dim )
    B[ind] .= - 1 ./ m[inputs]

    C = sparse( 1.0*I, 2*n, 2*n )
    D = spzeros( 2*n, dim )

    return (; :sys => HeteroStateSpace( A, B, C, D, Continuous() ) )
end