text:
export iauPmsafe
"""
Star proper motion:  update star catalog data for space motion, with
special handling to handle the zero parallax case.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
    ra1    double      right ascension (radians), before
    dec1   double      declination (radians), before
    pmr1   double      RA proper motion (radians/year), before
    pmd1   double      Dec proper motion (radians/year), before
    px1    double      parallax (arcseconds), before
    rv1    double      radial velocity (km/s, +ve = receding), before
    ep1a   double      "before" epoch, part A (Note 1)
    ep1b   double      "before" epoch, part B (Note 1)
    ep2a   double      "after" epoch, part A (Note 1)
    ep2b   double      "after" epoch, part B (Note 1)

Returned:
    ra2    double      right ascension (radians), after
    dec2   double      declination (radians), after
    pmr2   double      RA proper motion (radians/year), after
    pmd2   double      Dec proper motion (radians/year), after
    px2    double      parallax (arcseconds), after
    rv2    double      radial velocity (km/s, +ve = receding), after

Returned (function value):
        int         status:
                        -1 = system error (should not occur)
                        0 = no warnings or errors
                        1 = distance overridden (Note 6)
                        2 = excessive velocity (Note 7)
                        4 = solution didn't converge (Note 8)
                    else = binary logical OR of the above warnings

Notes:

    1. The starting and ending TDB epochs ep1a+ep1b and ep2a+ep2b are
    Julian Dates, apportioned in any convenient way between the two
    parts (A and B).  For example, JD(TDB)=2450123.7 could be
    expressed in any of these ways, among others:

        epNa            epNb

        2450123.7           0.0       (JD method)
        2451545.0       -1421.3       (J2000 method)
        2400000.5       50123.2       (MJD method)
        2450123.5           0.2       (date & time method)

    The JD method is the most natural and convenient to use in cases
    where the loss of several decimal digits of resolution is
    acceptable.  The J2000 method is best matched to the way the
    argument is handled internally and will deliver the optimum
    resolution.  The MJD method and the date & time methods are both
    good compromises between resolution and convenience.

    2. In accordance with normal star-catalog conventions, the object's
    right ascension and declination are freed from the effects of
    secular aberration.  The frame, which is aligned to the catalog
    equator and equinox, is Lorentzian and centered on the SSB.

    The proper motions are the rate of change of the right ascension
    and declination at the catalog epoch and are in radians per TDB
    Julian year.

    The parallax and radial velocity are in the same frame.

    3. Care is needed with units.  The star coordinates are in radians
    and the proper motions in radians per Julian year, but the
    parallax is in arcseconds.

    4. The RA proper motion is in terms of coordinate angle, not true
    angle.  If the catalog uses arcseconds for both RA and Dec proper
    motions, the RA proper motion will need to be divided by cos(Dec)
    before use.

    5. Straight-line motion at constant speed, in the inertial frame, is
    assumed.

    6. An extremely small (or zero or negative) parallax is overridden
    to ensure that the object is at a finite but very large distance,
    but not so large that the proper motion is equivalent to a large
    but safe speed (about 0.1c using the chosen constant).  A warning
    status of 1 is added to the status if this action has been taken.

    7. If the space velocity is a significant fraction of c (see the
    constant VMAX in the function iauStarpv), it is arbitrarily set
    to zero.  When this action occurs, 2 is added to the status.

    8. The relativistic adjustment carried out in the iauStarpv function
    involves an iterative calculation.  If the process fails to
    converge within a set number of iterations, 4 is added to the
    status.

Called:
    iauSeps      angle between two points
    iauStarpm    update star catalog data for space motion

This revision:   2014 July 1

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauPmsafe(double ra1, double dec1, double pmr1, double pmd1,
#               double px1, double rv1,
#               double ep1a, double ep1b, double ep2a, double ep2b,
#               double *ra2, double *dec2, double *pmr2, double *pmd2,
#               double *px2, double *rv2)

function iauPmsafe(ra1::Real, dec1::Real, pmr1::Real, pmd1::Real,
                 px1::Real, rv1::Real, ep1a::Real, ep1b::Real,
                 ep2a::Real, ep2b::Real)

    ref_ra2  = Ref{Float64}(0.0)
    ref_dec2 = Ref{Float64}(0.0)
    ref_pmr2 = Ref{Float64}(0.0)
    ref_pmd2 = Ref{Float64}(0.0)
    ref_px2  = Ref{Float64}(0.0)
    ref_rv2  = Ref{Float64}(0.0)

    status = ccall((:iauPmsafe, libsofa_c), Cint,
       (Cdouble, Cdouble, Cdouble, Cdouble,
        Cdouble, Cdouble, Cdouble,
        Cdouble, Cdouble, Cdouble,
        Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
        Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
       convert(Float64, ra1), convert(Float64, dec1),
       convert(Float64, pmr1), convert(Float64, pmd1),
       convert(Float64, px1), convert(Float64, rv1),
       convert(Float64, ep1a), convert(Float64, ep1b),
       convert(Float64, ep2a), convert(Float64, ep2b),
       ref_ra2, ref_dec2, ref_pmr2, ref_pmd2, ref_px2, ref_rv2)

   return status, ref_ra2[], ref_dec2[], ref_pmr2[], ref_pmd2[], ref_px2[], ref_rv2[]
end


text:
export eval_apl

# eval
eval_apl(ex) = eval_apl(ex, nothing, nothing)
eval_apl(f, α, ω) = f
eval_apl(v::JlVal, α, ω) = v.val
eval_apl(::Α, α, ω) = α
eval_apl(::Ω, α, ω) = ω
eval_apl(x::Apply, α, ω)  = eval_apl(x.f, α, ω)(eval_apl(x.r, α, ω))
eval_apl(x::ConcArr, α, ω)  = vcat(eval_apl(x.l, α, ω), eval_apl(x.r, α, ω))
eval_apl(x::Apply2, α, ω) = eval_apl(x.f, α, ω)(eval_apl(x.l, α, ω), eval_apl(x.r, α, ω))

# call methods for primitive functions
mkbody1(x::Symbol) = :($x(ω))
mkbody1(x::Expr) = x
mkbody2(x::Symbol) = :($x(α, ω))
mkbody2(x::Expr) = x
for (sym, fns) in prim_fns
    mon, dya = fns
    mon != nothing && @eval (f::PrimFn{$sym})(ω) = $(mkbody1(mon))
    dya != nothing && @eval (f::PrimFn{$sym})(α, ω) = $(mkbody2(dya))
end

function actuallyreducedim(f, xs::AbstractArray, ident)
    dropdims(
        reduce(f, xs, dims=ndims(xs), init=ident),
        dims=ndims(xs)
    )
end

function actuallyreducedim(f, xs::AbstractVector, ident)
    reduce(f, xs, init=ident)
end

# call methods for primitive operators
(op::Op1{'/'})(ω) = actuallyreducedim(op.l, ω, identity(op.l, eltype(ω)))
(op::Op1{'⌿'})(ω) = dropdims(reduce(op.l, ω, init=identity(op.l, eltype(ω)), dims=1), dims=1)
(op::Op1{'\\'})(ω) = prefix_scan(op.l, ω, identity(op.l, ω))
(op::Op1{'⍀'})(ω) = prefix_scan(op.l, ω, identity(op.l, ω)) # Todo
(op::Op1{'¨'})(ω) = map(op.l, ω)
(op::Op1{'↔'})(α, ω) = op.l(ω, α)
(op::Op1{'⍨'})(α, ω) = op.l(ω, α)
(op::Op2{'.'})(α, ω) = reduce(op.l, op.r(convert(Array, α), convert(Array, ω)))
(op::Op2{'⋅'})(α) = op.l(op.r(α)) # compose
(op::Op1{'∘'})(α, ω) = [op.l(x, y) for x in α, y in ω]

# user defined functions
(fn::UDefFn{0})() = eval_apl(fn.ast)
(fn::UDefFn{1})(ω) = eval_apl(fn.ast, nothing, ω)
(fn::UDefFn{2})(α, ω) = eval_apl(fn.ast, α, ω)

ext: jl


text:
Dlist = (2,4,12,29,71)

constant1(x) = Constant(x)
constant2(x) = FourierSeries([x])

@testset "UniformCMPS: environments and gauging with bond dimension $D" for D in Dlist
    for constant in (constant1, constant2)
        for T in (Float64, ComplexF64)
            Q = constant(randn(T, (D,D)))
            R1 = constant(randn(T, (D,D)))
            R2 = constant(randn(T, (D,D)))
            Ψ = InfiniteCMPS(Q, (R1,R2))
            ρL, λL, infoL = leftenv(Ψ)
            TL = LeftTransfer(Ψ)
            @test norm(TL(ρL) - 2*λL*ρL) <= 10*infoL.normres[1]
            @test ρL == ρL'

            ρR, λR, infoR = rightenv(Ψ)
            TR = RightTransfer(Ψ)
            @test norm(TR(ρR) - 2*λR*ρR) <= 10*infoR.normres[1]
            @test λL ≈ λR
            @test ρR == ρR'

            Z = dot(ρL, ρR)
            @test imag(Z) <= sqrt(infoL.normres[1]*infoR.normres[1])

            Ψn = copy(Ψ)
            ρL2, = leftenv!(Ψn)
            @test ρL2 ≈ ρL
            Qn = Ψn.Q
            R1n, R2n = Ψn.Rs
            @test norm(Qn'*ρL + ρL*Qn + R1n'*ρL*R1n + R2n'*ρL*R2n) <= 1e-9

            Ψn = copy(Ψ)
            ρR2, = rightenv!(Ψn)
            @test ρR2 ≈ ρR
            Qn = Ψn.Q
            R1n, R2n = Ψn.Rs
            @test norm(Qn*ρR + ρR*Qn' + R1n*ρR*R1n' + R2n*ρR*R2n') <= 1e-9

            Ψn = copy(Ψ)
            ρL2, ρR2 = environments!(Ψn)
            @test ρL2 ≈ ρL/sqrt(Z)
            @test ρR2 ≈ ρR/sqrt(Z)
            @test norm(LeftTransfer(Ψn)(ρL2)) <= 1e-9
            @test norm(RightTransfer(Ψn)(ρR2)) <= 1e-9

            if constant == constant1 # actual UniformCMPS
                ΨL,λ,CL,info = leftgauge(Ψ)
                @test λ ≈ λL
                @test CL'*CL ≈ ρL/tr(ρL[])
                QL = ΨL.Q
                R1L, R2L = ΨL.Rs
                @test CL*Qn ≈ QL*CL
                @test CL*R1 ≈ R1L*CL
                @test CL*R2 ≈ R2L*CL
                @test norm(QL + QL' + R1L'*R1L + R2L'*R2L) <= 1e-9

                ΨR,λ,CR,info = rightgauge(Ψ; tol = 1e-12)
                @test λ ≈ λR
                @test CR*CR' ≈ ρR/tr(ρR[])
                QR = ΨR.Q
                R1R, R2R = ΨR.Rs
                @test Qn*CR ≈ CR*QR
                @test R1*CR ≈ CR*R1R
                @test R2*CR ≈ CR*R2R
                @test norm(QR + QR' + R1R*R1R' + R2R*R2R') <= 1e-9
            end
        end
    end
end

@testset "UniformCMPS: energy environments with bond dimension $D" for D in Dlist
    H = ∫(∂ψ'*∂ψ - 0.5*ψ'*ψ + 0.3*(ψ*ψ + ψ'*ψ') + 2.5*(ψ')^2*ψ^2, (-Inf,+Inf))
    for constant in (constant1, constant2)
        for T in (Float64, ComplexF64)
            Q = constant(randn(T, (D,D)))
            R = constant(randn(T, (D,D)))
            Ψ = InfiniteCMPS(Q, R)

            ρL, ρR = environments!(Ψ)
            HL, EL, eL, hL = leftenv(H, (Ψ,ρL,ρR))
            HR, ER, eR, hR = rightenv(H, (Ψ,ρL,ρR))

            @test eL ≈ eR
            @test eL(0) ≈ EL
            @test eR(0) ≈ ER
            @test norm(LeftTransfer(Ψ)(HL) + hL) <= 1e-9*norm(HL)
            @test norm(RightTransfer(Ψ)(HR) + hR) <= 1e-9*norm(HR)
            @test abs(dot(HL,ρR)) <= 1e-9*norm(HL)
            @test abs(dot(ρL,HR)) <= 1e-9*norm(HR)
        end
    end
end

@testset "UniformCMPS: local gradients with bond dimension D = $D" for D in Dlist
    for constant in (constant1, constant2)
        for T in (Float64, ComplexF64)
            Q = constant(randn(T, (D,D)))
            R1 = constant(randn(T, (D,D)))
            R2 = constant(randn(T, (D,D)))
            ρL = constant(randn(T, (D,D)))
            ρL = (ρL + ρL')/2 # does not need to be the actual ρL for this test
            ρR = constant(randn(T, (D,D)))
            ρR = (ρR + ρR')/2 # does not need to be the actual ρL for this test
            Rs = (R1,R2)

            QR1 = Q*R1 - R1*Q
            QR2 = Q*R2 - R2*Q

            @test CMPSKit.localgradientQ(ψ[1], Q, Rs, ρL, ρR) == zero(Q)
            @test CMPSKit.localgradientQ(ψ[1]^2, Q, Rs, ρL, ρR) == zero(Q)
            @test CMPSKit.localgradientQ(∂ψ[1], Q, Rs, ρL, ρR) == zero(Q)
            @test CMPSKit.localgradientQ(ψ[2]', Q, Rs, ρL, ρR) == zero(Q)
            @test CMPSKit.localgradientQ((ψ[2]')^2, Q, Rs, ρL, ρR) == zero(Q)
            @test CMPSKit.localgradientQ(∂(ψ[2]'), Q, Rs, ρL, ρR) ≈ ρL*ρR*R2' - R2'*ρL*ρR
            @test CMPSKit.localgradientQ((∂ψ[1])'*∂ψ[2], Q, Rs, ρL, ρR) ≈
                    ρL*QR2*ρR*R1' - R1'*ρL*QR2*ρR
            @test CMPSKit.localgradientQ(ψ[1]'*ψ[2], Q, Rs, ρL, ρR) == zero(Q)
            @test CMPSKit.localgradientQ((ψ[1]')^2*ψ[1]^2, Q, Rs, ρL, ρR) == zero(Q)

            @test CMPSKit.localgradientRs(ψ[1], Q, Rs, ρL, ρR) == (zero(Q), zero(Q))
            @test CMPSKit.localgradientRs(ψ[1]^2, Q, Rs, ρL, ρR) == (zero(Q), zero(Q))
            @test CMPSKit.localgradientRs(∂ψ[1], Q, Rs, ρL, ρR) == (zero(Q), zero(Q))
            @test all(isapprox.(CMPSKit.localgradientRs(ψ[2]', Q, Rs, ρL, ρR),
                                (zero(Q), ρL*ρR)))
            @test all(isapprox.(CMPSKit.localgradientRs((ψ[2]')^2, Q, Rs, ρL, ρR),
                                (zero(Q), R2'*ρL*ρR + ρL*ρR*R2')))
            @test all(isapprox.(CMPSKit.localgradientRs(∂(ψ[2]'), Q, Rs, ρL, ρR),
                                (zero(Q), Q'*ρL*ρR - ρL*ρR*Q' - ∂(ρL*ρR))))
            @test all(isapprox.(CMPSKit.localgradientRs(∂(ψ[2]'), Q, Rs, ρL, ρR),
                                (zero(Q), Q'*ρL*ρR - ρL*ρR*Q' - ∂(ρL*ρR))))
            @test all(isapprox.(CMPSKit.localgradientRs((∂ψ[1])'*∂ψ[2], Q, Rs, ρL, ρR),
                                (Q'*ρL*QR2*ρR - ρL*QR2*ρR*Q' - ∂(ρL*QR2*ρR), zero(Q))))
            @test all(isapprox.(CMPSKit.localgradientRs(ψ[1]'*ψ[2], Q, Rs, ρL, ρR),
                                (ρL*R2*ρR, zero(Q))))
            @test all(isapprox.(CMPSKit.localgradientRs((ψ[1]')^2*ψ[1]^2, Q, Rs, ρL, ρR),
                                (R1'*ρL*R1*R1*ρR + ρL*R1*R1*ρR*R1', zero(Q))))

        end
    end
end

@testset "UniformCMPS: global gradients with bond dimension D = $D" for D in Dlist
    α = rand()
    β = rand()
    γ = rand()
    H = ∫(∂ψ'*∂ψ + α*ψ'*ψ + β*(ψ*ψ + ψ'*ψ') + γ*(ψ')^2*ψ^2, (-Inf,+Inf))
    for T in (Float64, ComplexF64)
        Q = Constant(randn(T, (D,D)))
        R = Constant(randn(T, (D,D)))

        ΨL, = leftgauge!(InfiniteCMPS(Q, R))
        QL = ΨL.Q
        RL = ΨL.Rs[1]
        ρR, = rightenv(ΨL)
        ρL = one(ρR)
        HL, = leftenv(H, (ΨL, ρL, ρR))
        HR, = rightenv(H, (ΨL, ρL, ρR))

        gradQ, gradRs = gradient(H, (ΨL, ρL, ρR), HL, HR)

        QRL = QL*RL-RL*QL
        gradR = gradRs[1]

        @test gradQ ≈ ρL*QRL*ρR*RL' - RL'*ρL*QRL*ρR + HL*ρR + ρL*HR

        @test gradR ≈ QL'*ρL*QRL*ρR - ρL*QRL*ρR*QL' + α*ρL*RL*ρR +
                        β*RL'*ρL*ρR + β*ρL*ρR*RL' + γ*RL'*ρL*RL*RL*ρR + γ*ρL*RL*RL*ρR*RL' +
                        HL*RL*ρR + ρL*RL*HR
    end
end

@testset "UniformCMPS: test ground state algorithm" begin
    D = 4
    T = Float64
    α = -1
    β = 1//2
    γ = float(pi)
    H = ∫(∂ψ'*∂ψ + α*ψ'*ψ + β*(ψ*ψ + ψ'*ψ') + γ*(ψ')^2*ψ^2, (-Inf,+Inf))

    gradtol = 1e-7
    optalg = ConjugateGradient(; gradtol = gradtol)
    eigalg = Arnoldi(; krylovdim = 16, tol = 1e-10)
    linalg = GMRES(; krylovdim = 16, tol = 1e-10)
    for k = 1:3
        Q = Constant(randn(T, (D,D)))
        R = Constant(randn(T, (D,D)))

        ΨL, ρR, E, e, normgrad, numfg, history =
            groundstate(H, InfiniteCMPS(Q, R);
                        optalg = optalg, eigalg = eigalg, linalg = linalg)
        @test E ≈ -0.43306384063961445
        @test abs(expval(ψ, ΨL, one(ρR), ρR)(0)) < gradtol
        @test expval(ψ'*ψ, ΨL, one(ρR), ρR)(0) ≈ 0.5086468402292694 atol=gradtol
    end
end


text:
import Base.isempty

export ConvexHull, CH,
       convex_hull,
       convex_hull!,
       ConvexHullArray, CHArray,
       ConvexHull!,
       array,
       swap

"""
    ConvexHull{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}

Type that represents the convex hull of the union of two convex sets.

### Fields

- `X` -- convex set
- `Y` -- convex set

### Notes

The `EmptySet` is the neutral element for `ConvexHull`.

### Examples

Convex hull of two 100-dimensional Euclidean balls:

```jldoctest
julia> b1, b2 = Ball2(zeros(100), 0.1), Ball2(4*ones(100), 0.2);

julia> c = ConvexHull(b1, b2);

julia> typeof(c)
ConvexHull{Float64,Ball2{Float64},Ball2{Float64}}
```
"""
struct ConvexHull{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}
    X::S1
    Y::S2

    # default constructor with dimension check
    function ConvexHull{N, S1, S2}(X::S1, Y::S2) where
            {N<:Real, S1<:LazySet{N}, S2<:LazySet{N}}
        @assert dim(X) == dim(Y) "sets in a convex hull must have the same " *
            "dimension"
        return new{N, S1, S2}(X, Y)
    end
end

isoperationtype(::Type{<:ConvexHull}) = true

# convenience constructor without type parameter
ConvexHull(X::S1, Y::S2) where {N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} =
    ConvexHull{N, S1, S2}(X, Y)

# EmptySet is the neutral element for ConvexHull
@neutral(ConvexHull, EmptySet)

# Universe is the absorbing element for ConvexHull
@absorbing(ConvexHull, Universe)

"""
    CH

Alias for `ConvexHull`.
"""
const CH = ConvexHull

"""
    swap(ch::ConvexHull)

Return a new `ConvexHull` object with the arguments swapped.

### Input

- `ch` -- convex hull of two convex sets

### Output

A new `ConvexHull` object with the arguments swapped.
"""
function swap(ch::ConvexHull)
    return ConvexHull(ch.Y, ch.X)
end

"""
    dim(ch::ConvexHull)::Int

Return the dimension of a convex hull of two convex sets.

### Input

- `ch` -- convex hull of two convex sets

### Output

The ambient dimension of the convex hull of two convex sets.
"""
function dim(ch::ConvexHull)::Int
    return dim(ch.X)
end

"""
    σ(d::AbstractVector{N}, ch::ConvexHull{N}) where {N<:Real}

Return the support vector of a convex hull of two convex sets in a given
direction.

### Input

- `d`  -- direction
- `ch` -- convex hull of two convex sets

### Output

The support vector of the convex hull in the given direction.
"""
function σ(d::AbstractVector{N}, ch::ConvexHull{N}) where {N<:Real}
    σ1 = σ(d, ch.X)
    σ2 = σ(d, ch.Y)
    ρ1 = dot(d, σ1)
    ρ2 = dot(d, σ2)
    return ρ1 >= ρ2 ? σ1 : σ2
end

"""
    ρ(d::AbstractVector{N}, ch::ConvexHull{N}) where {N<:Real}

Return the support function of a convex hull of two convex sets in a given
direction.

### Input

- `d`  -- direction
- `ch` -- convex hull of two convex sets

### Output

The support function of the convex hull in the given direction.
"""
function ρ(d::AbstractVector{N}, ch::ConvexHull{N}) where {N<:Real}
    return max(ρ(d, ch.X), ρ(d, ch.Y))
end

"""
    isbounded(ch::ConvexHull)::Bool

Determine whether a convex hull of two convex sets is bounded.

### Input

- `ch` -- convex hull of two convex sets

### Output

`true` iff both wrapped sets are bounded.
"""
function isbounded(ch::ConvexHull)::Bool
    return isbounded(ch.X) && isbounded(ch.Y)
end

"""
    isempty(ch::ConvexHull)::Bool

Return if a convex hull of two convex sets is empty or not.

### Input

- `ch` -- convex hull

### Output

`true` iff both wrapped sets are empty.
"""
function isempty(ch::ConvexHull)::Bool
    return isempty(ch.X) && isempty(ch.Y)
end

# ================================
#  Convex hull of an array of sets
# ================================
"""
    ConvexHullArray{N<:Real, S<:LazySet{N}} <: LazySet{N}

Type that represents the symbolic convex hull of a finite number of convex sets.

### Fields

- `array` -- array of sets

### Notes

The `EmptySet` is the neutral element for `ConvexHullArray`.

Constructors:

- `ConvexHullArray(array::Vector{<:LazySet})` -- default constructor

- `ConvexHullArray([n]::Int=0, [N]::Type=Float64)`
  -- constructor for an empty hull with optional size hint and numeric type

### Examples

Convex hull of 100 two-dimensional balls whose centers follows a sinusoidal:

```jldoctest
julia> b = [Ball2([2*pi*i/100, sin(2*pi*i/100)], 0.05) for i in 1:100];

julia> c = ConvexHullArray(b);
```
"""
struct ConvexHullArray{N<:Real, S<:LazySet{N}} <: LazySet{N}
    array::Vector{S}
end

isoperationtype(::Type{<:ConvexHullArray}) = true

# constructor for an empty hull with optional size hint and numeric type
function ConvexHullArray(n::Int=0, N::Type=Float64)::ConvexHullArray
    a = Vector{LazySet{N}}()
    sizehint!(a, n)
    return ConvexHullArray(a)
end

# EmptySet is the neutral element for ConvexHullArray
@neutral(ConvexHullArray, EmptySet)

# Universe is the absorbing element for ConvexHullArray
@absorbing(ConvexHullArray, Universe)

# add functions connecting ConvexHull and ConvexHullArray
@declare_array_version(ConvexHull, ConvexHullArray)

"""
    CHArray

Alias for `ConvexHullArray`.
"""
const CHArray = ConvexHullArray

"""
    array(cha::ConvexHullArray{N, S})::Vector{S} where {N<:Real, S<:LazySet{N}}

Return the array of a convex hull of a finite number of convex sets.

### Input

- `cha` -- convex hull array

### Output

The array of a convex hull of a finite number of convex sets.
"""
function array(cha::ConvexHullArray{N, S})::Vector{S} where {N<:Real, S<:LazySet{N}}
    return cha.array
end

"""
    dim(cha::ConvexHullArray)::Int

Return the dimension of the convex hull of a finite number of convex sets.

### Input

- `cha` -- convex hull array

### Output

The ambient dimension of the convex hull of a finite number of convex sets.
"""
function dim(cha::ConvexHullArray)::Int
    @assert !isempty(cha.array)
    return dim(cha.array[1])
end

"""
    σ(d::AbstractVector{N}, cha::ConvexHullArray{N}) where {N<:Real}

Return the support vector of a convex hull array in a given direction.

### Input

- `d`   -- direction
- `cha` -- convex hull array
"""
function σ(d::AbstractVector{N}, cha::ConvexHullArray{N}) where {N<:Real}
    s = σ(d, cha.array[1])
    ri = dot(d, s)
    rmax = ri
    for (i, chi) in enumerate(cha.array[2:end])
        si = σ(d, chi)
        ri = dot(d, si)
        if ri > rmax
            rmax = ri
            s = si
        end
    end
    return s
end

"""
    ρ(d::AbstractVector{N}, cha::ConvexHullArray{N}) where {N<:Real}

Return the support function of a convex hull array in a given direction.

### Input

- `d`   -- direction
- `cha` -- convex hull array

### Output

The support function of the convex hull array in the given direction.

### Algorithm

This algorihm calculates the maximum over all ``ρ(d, X_i)`` where the
``X_1, …, X_k`` are the sets in the array `cha`.
"""
function ρ(d::AbstractVector{N}, cha::ConvexHullArray{N}) where {N<:Real}
    return maximum([ρ(d, Xi) for Xi in array(cha)])
end

"""
    isbounded(cha::ConvexHullArray)::Bool

Determine whether a convex hull of a finite number of convex sets is
bounded.

### Input

- `cha` -- convex hull of a finite number of convex sets

### Output

`true` iff all wrapped sets are bounded.
"""
function isbounded(cha::ConvexHullArray)::Bool
    return all(x -> isbounded(x), cha.array)
end

"""
    isempty(cha::ConvexHullArray)::Bool

Return if a convex hull array is empty or not.

### Input

- `cha` -- convex hull array

### Output

`true` iff all wrapped sets are empty.
"""
function isempty(cha::ConvexHullArray)::Bool
    return all(X -> isempty(X), array(cha))
end


text:
using MatplotWrap

N = 10
x = rand(N)
y = rand(N)

scatter(x, y) # matplot::scatter(x, y)
matplot_show() # matplot::show()


text:
"""
    FEMMDeforNonlinearExplModule

Base module for operations on interiors of domains to construct
system vectors for explicit-dynamics nonlinear deformation models.
"""
module FEMMDeforNonlinearExplModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.FENodeSetModule: FENodeSet
import FinEtools.FESetModule: AbstractFESet, gradN!, nodesperelem, manifdim
import FinEtools.IntegDomainModule: IntegDomain, integrationdata, Jacobianvolume
import FinEtools.FieldModule: ndofs, gatherdofnums!, gatherfixedvalues_asvec!, gathervalues_asvec!, gathervalues_asmat!
import FinEtools.NodalFieldModule: NodalField, nnodes
import FinEtools.AssemblyModule: AbstractSysvecAssembler, AbstractSysmatAssembler, SysmatAssemblerSparseSymm, startassembly!, assemble!, makematrix!, makevector!, SysvecAssembler
import ..FEMMDeforNonlinearBaseModule: AbstractFEMMDeforNonlinear
import FinEtools.CSysModule: CSys, updatecsmat!
import FinEtoolsDeforLinear.DeforModelRedModule: AbstractDeforModelRed, DeforModelRed2DAxisymm, nstressstrain, nthermstrain, Blmat!
import FinEtoolsDeforLinear.MatDeforModule: rotstressvec!
import FinEtools.MatModule: massdensity
import ..MatDeforNonlinearModule: AbstractMatDeforNonlinear, tangentmoduli!, update!, newstate
import FinEtools.SurfaceNormalModule: SurfaceNormal, updatenormal!
import FinEtools.MatrixUtilityModule: add_btdb_ut_only!, complete_lt!, locjac!, add_nnt_ut_only!, add_gkgt_ut_only!, add_btsigma!, jac!
import ..FEMMDeforNonlinearBaseModule: restoringforce
import FinEtools.MatrixUtilityModule: mulCAB!, mulCAtB!, mulCABt!, detC

rotF!(Fr, F, RmTF, Rm) = begin
mulCAtB!(Val(3), RmTF, Rm, F); mulCAB!(Val(3), Fr, RmTF, Rm)
end

struct _Buffers
	initialized::Bool
	dofnums::FIntVec
	loc::FFltMat
	J::FFltMat
	csmatTJ::FFltMat
	gradXN::FFltMat
	gradxmN::FFltMat
	gradxN::FFltMat
	D::FFltMat
	B::FFltMat
	DB::FFltMat
	elmat::FFltMat
	elvec::FFltVec
	pu::FFltVec
	X::FFltMat
	xn::FFltMat
	xn1::FFltMat
	Un::FFltMat
	Un1::FFltMat
	Fn::FFltMat
	Fn1::FFltMat
	Fnm::FFltMat
	Fn1m::FFltMat
	RmTF::FFltMat
	cauchy::FFltVec
	output::FFltVec
	sigma::FFltMat
	idx::Array{StepRange{Int64,Int64},1}
	c1::FFltMat
	sg::FFltMat
end

struct _Bfuns
	npts::FInt
	Ns::Array{Array{Float64,2},2}
	gradNparams::Array{Array{Float64,2},2}
	w::FFltMat
	pc::FFltMat
end

function _makebuffers(initialized, ndn, nne, sdim, mdim, nstrs)
    elmatdim = ndn*nne;             # dimension of the element matrix
    # Prepare _buffers
    dofnums = zeros(FInt, elmatdim); # degree of freedom array -- buffer
    loc = fill(zero(FFlt), 1, sdim); # quadrature point location -- buffer
    J = fill(zero(FFlt), sdim, mdim); # Jacobian matrix -- buffer
    csmatTJ = fill(zero(FFlt), mdim, mdim); # intermediate result -- buffer
    gradXN = fill(zero(FFlt), nne, mdim); # intermediate result -- buffer
    gradxmN = fill(zero(FFlt), nne, mdim); # intermediate result -- buffer
    gradxN = fill(zero(FFlt), nne, mdim); # intermediate result -- buffer
    # Prepare _buffers
    D = fill(zero(FFlt), nstrs, nstrs); # material stiffness matrix -- buffer
    B = fill(zero(FFlt), nstrs, elmatdim); # strain-displacement matrix -- buffer
    DB = fill(zero(FFlt), nstrs, elmatdim); # strain-displacement matrix -- buffer
    elmat = fill(zero(FFlt), elmatdim, elmatdim);      # element matrix -- buffer
    # Prepare _buffers
    X = fill(zero(FFlt), nne, sdim)
    xn = fill(zero(FFlt), nne, sdim)
    xn1 = fill(zero(FFlt), nne, sdim)
    Un = fill(zero(FFlt), nne, sdim)
    Un1 = fill(zero(FFlt), nne, sdim)
    Fn = fill(zero(FFlt), sdim, mdim);
    Fn1 = fill(zero(FFlt), sdim, mdim);
    Fnm = fill(zero(FFlt), sdim, mdim);
    Fn1m = fill(zero(FFlt), sdim, mdim);
    RmTF = fill(zero(FFlt), sdim, mdim);
    elvec = fill(zero(FFlt), elmatdim);      # element matrix -- buffer
    pu = fill(zero(FFlt), elmatdim);      # element matrix -- buffer
    # Prepare _buffers
    cauchy = fill(zero(FFlt), nstrs);      # element matrix -- buffer
    output = fill(zero(FFlt), nstrs);      # element matrix -- buffer
    sigma = fill(zero(FFlt), sdim, sdim);      # element matrix -- buffer
    # Prepare _buffers
    idx = [(1:ndn:(nne-1)*ndn+1).+(i-1) for i in 1:sdim]; # indexing vector -- buffer
    # to compute c1 = gradxmN*sigma*gradxmN';
    c1 = fill(zero(FFlt), nne, nne); # strain-displacement matrix -- buffer
    sg = fill(zero(FFlt), sdim, nne); # strain-displacement matrix -- buffer
    return _buffers = _Buffers(initialized,
    	dofnums, loc, J, csmatTJ, gradXN, gradxN, gradxmN,
    	D, B, DB, elmat, elvec, pu,
    	X, xn, xn1, Un, Un1, Fn, Fn1, Fnm, Fn1m, RmTF,
    	cauchy, output, sigma,
    	idx, c1, sg
    	)
end

"""
    FEMMDeforNonlinearExpl <: AbstractFEMMDeforNonlinear

Type of FEMM for explicit-dynamics simulations of nonlinear deformation.

!!! note
Julia threads do not cooperate with BLAS threads. All the BLAS calls are
eliminated from the implementation of this FEMM in order to obtain good parallel
efficiency.
"""
mutable struct FEMMDeforNonlinearExpl{MR<:AbstractDeforModelRed,  S<:AbstractFESet, F<:Function, M<:AbstractMatDeforNonlinear} <: AbstractFEMMDeforNonlinear
	mr::Type{MR} # model reduction type
	integdomain::IntegDomain{S, F} # integration domain data
	mcsys::CSys # updater of the material orientation matrix
	material::M # material object
	statev::Vector{Vector{FFltVec}} # vector of state vectors, one for each element, and then for each integration point
	_buffers::_Buffers # workspace: pre-allocate so that the memory is handled efficiently
	_bfuns::_Bfuns # workspace: pre-allocate so that the memory is handled efficiently
end

"""
    FEMMDeforNonlinearExpl(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:AbstractDeforModelRed, S<:AbstractFESet, F<:Function, M<:AbstractMatDeforNonlinear}

Constructor of nonlinear deformation finite element modeling machine.

The material coordinate system is the default.
"""
function FEMMDeforNonlinearExpl(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:AbstractDeforModelRed, S<:AbstractFESet, F<:Function, M<:AbstractMatDeforNonlinear}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (integdomain.axisymmetric) || (mr != DeforModelRed2DAxisymm) "Axially symmetric requires axisymmetric to be true"
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(integdomain);
    statev = Vector{FFltVec}[]
    for i in 1:count(integdomain.fes)
        push!(statev, [newstate(material) for j in 1:npts])
    end
    npts, Ns, gradNparams, w, pc = integrationdata(integdomain);
    return FEMMDeforNonlinearExpl(mr, integdomain, CSys(manifdim(integdomain.fes)), material, statev, _makebuffers(false, 1, 1, 1, 1, 1), _Bfuns(npts, Ns, gradNparams, w, pc))
end

"""
    FEMMDeforNonlinearExpl(mr::Type{MR}, integdomain::IntegDomain{S, F}, material::M) where {MR<:AbstractDeforModelRed, S<:AbstractFESet, F<:Function, M<:AbstractMatDeforNonlinear}

Constructor of nonlinear deformation finite element modeling machine.
"""
function FEMMDeforNonlinearExpl(mr::Type{MR}, integdomain::IntegDomain{S, F}, mcsys::CSys, material::M) where {MR<:AbstractDeforModelRed, S<:AbstractFESet, F<:Function, M<:AbstractMatDeforNonlinear}
    @assert mr == material.mr "Model reduction is mismatched"
    @assert (integdomain.axisymmetric) || (mr != DeforModelRed2DAxisymm) "Axially symmetric requires axisymmetric to be true"
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(integdomain);
    statev = Vector{FFltVec}[]
    for i in 1:count(integdomain.fes)
        push!(statev, [newstate(material) for j in 1:npts])
    end
    npts, Ns, gradNparams, w, pc = integrationdata(integdomain);
    return FEMMDeforNonlinearExpl(mr, integdomain, mcsys, material, statev, _makebuffers(false, 1, 1, 1, 1, 1), _Bfuns(npts, Ns, gradNparams, w, pc))
end

function _buff1(self::FEMMDeforNonlinearExpl, geom::NodalField, u::NodalField)
    return self._buffers.dofnums, self._buffers.loc, self._buffers.J, self._buffers.csmatTJ, self._buffers.gradXN, self._buffers.gradxN, self._buffers.gradxmN
end

function _buff2(self::FEMMDeforNonlinearExpl, geom::NodalField, u::NodalField)
    return self._buffers.D, self._buffers.B, self._buffers.DB, self._buffers.elmat
end

function _buff3(self::FEMMDeforNonlinearExpl, geom::NodalField, u::NodalField)
    return self._buffers.X, self._buffers.xn, self._buffers.xn1, self._buffers.Un, self._buffers.Un1, self._buffers.Fn, self._buffers.Fn1, self._buffers.Fnm, self._buffers.Fn1m, self._buffers.RmTF
end

function _buff4(self::FEMMDeforNonlinearExpl, geom::NodalField, u::NodalField)
    return self._buffers.B, self._buffers.elvec, self._buffers.pu
end

function _buff5(self::FEMMDeforNonlinearExpl, geom::NodalField, u::NodalField)
    return self._buffers.cauchy, self._buffers.output, self._buffers.sigma
end

function _buff6(self::FEMMDeforNonlinearExpl, geom::NodalField, u::NodalField)
    return self._buffers.idx, self._buffers.elmat, self._buffers.c1, self._buffers.sg
end

function _makebuffers!(self::FEMMDeforNonlinearExpl,  geom::NodalField{FFlt}, u::NodalField)
	if !self._buffers.initialized
		fes = self.integdomain.fes
		ndn = ndofs(u); # number of degrees of freedom per node
		nne = nodesperelem(fes); # number of nodes for element
		sdim = ndofs(geom);            # number of space dimensions
		mdim = manifdim(fes); # manifold dimension of the element
		nstrs = nstressstrain(self.mr);  # number of stresses
		self._buffers = _makebuffers(true, ndn, nne, sdim, mdim, nstrs)
	end
    return self # default is no-op
end

"""
    restoringforce(self::AbstractFEMMDeforNonlinear, assembler::A, geom::NodalField{FFlt}, un1::NodalField{T}, un::NodalField{T}, tn::FFlt, dtn::FFlt, savestate = false) where {A<:AbstractSysvecAssembler, T<:Number}

Compute the restoring force vector.

Note: This method *UPDATES* the state of the FEMM object.  In
particular, the material state gets updated (refer to `savestate`).

# Arguments
- `assembler` = vector assembler,
- `geom` = geometry field: reference coordinates of the nodes,
- `un1` = displacement field at time `tn1 = tn + dtn`,
- `un` = displacement field at time `tn`,
- `tn` = time in step `n`
- `dtn` = increment of time
- `savestate` = bool flag: should we modify the material states (`savestate = true`)? Otherwise work with a copy of the material state.
"""
function restoringforce(self::FEMMDeforNonlinearExpl, assembler::A, geom::NodalField{FFlt}, un1::NodalField{T}, un::NodalField{T}, tn::FFlt, dtn::FFlt, savestate = false) where {A<:AbstractSysvecAssembler, T<:Number}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = self._bfuns.npts,  self._bfuns.Ns,  self._bfuns.gradNparams,  self._bfuns.w,  self._bfuns.pc;
    _makebuffers!(self, geom, un1)
    dofnums, loc, J, csmatTJ, gradXN, gradxN, gradxmN = _buff1(self, geom, un1)
    B, elvec, pu = _buff4(self, geom, un1)
    X, xn, xn1, Un, Un1, Fn, Fn1, Fnm, Fn1m, RmTF = _buff3(self, geom, un1)
    cauchy, output, sigma = _buff5(self, geom, un1)
    statev = self.statev # Initially we are assuming that the state can be overridden
    if !savestate  # unless we are instructed to work with a copy
        statev = deepcopy(self.statev) # work with copies of material state
    end
    startassembly!(assembler, un1.nfreedofs);
    for i = 1:count(fes) # Loop over elements
        gathervalues_asmat!(geom, X, fes.conn[i]);
        gathervalues_asmat!(un, Un, fes.conn[i]);
        gathervalues_asmat!(un1, Un1, fes.conn[i]);
        @. xn = X + Un; # previous known coordinates
        @. xn1 = X + Un1; # current coordinates
        fill!(elvec,  0.0); # Initialize element matrix
        for j = 1:npts # Loop over quadrature points
            locjac!(loc, J, X, Ns[j], gradNparams[j])
            Jac = Jacobianvolume(self.integdomain, J, loc, fes.conn[i], Ns[j]);
            updatecsmat!(self.mcsys, loc, J, fes.label[i]); Rm = self.mcsys.csmat
            gradN!(fes, gradXN, gradNparams[j], J); # wrt global c.s.
            jac!(Fn, xn, gradXN) # Previous deformation gradient
            jac!(Fn1, xn1, gradXN)
            rotF!(Fnm, Fn, RmTF, Rm) # wrt local c.s.
            rotF!(Fn1m, Fn1, RmTF, Rm) # wrt local c.s.
            update!(self.material, statev[i][j], cauchy, output, Fn1m, Fnm, tn, dtn, loc, fes.label[i])
            gradN!(fes, gradxN, gradXN, Fn1); # wrt current c.s.
            mulCAB!(gradxmN, gradxN, Rm) # wrt material c.s.
            Blmat!(self.mr, B, Ns[j], gradxmN, loc, Rm); # local strain-global disp
            add_btsigma!(elvec, B, - Jac * w[j] * detC(Val(3), Fn1), cauchy)
        end # Loop over quadrature points
        gatherdofnums!(un1, dofnums, fes.conn[i]); # retrieve degrees of freedom
        assemble!(assembler, elvec, dofnums); # assemble symmetric matrix
    end # Loop over elements
    return makevector!(assembler);
end

function restoringforce(self::FEMMDeforNonlinearExpl, geom::NodalField{FFlt}, un1::NodalField{T}, un::NodalField{T}, tn::FFlt, dtn::FFlt, savestate = false) where {T<:Number}
    assembler = SysvecAssembler();
    return restoringforce(self, assembler, geom, un1, un, tn, dtn, savestate);
end

end


text:
###########################################################
#
#   (non-weighted) sampling
#
###########################################################

using Random: Sampler, Random.GLOBAL_RNG

### Algorithms for sampling with replacement

function direct_sample!(rng::AbstractRNG, a::UnitRange, x::AbstractArray)
    s = Sampler(rng, 1:length(a))
    b = a[1] - 1
    if b == 0
        for i = 1:length(x)
            @inbounds x[i] = rand(rng, s)
        end
    else
        for i = 1:length(x)
            @inbounds x[i] = b + rand(rng, s)
        end
    end
    return x
end
direct_sample!(a::UnitRange, x::AbstractArray) = direct_sample!(Random.GLOBAL_RNG, a, x)

"""
    direct_sample!([rng], a::AbstractArray, x::AbstractArray)

Direct sampling: for each `j` in `1:k`, randomly pick `i` from `1:n`,
and set `x[j] = a[i]`, with `n=length(a)` and `k=length(x)`.

This algorithm consumes `k` random numbers.
"""
function direct_sample!(rng::AbstractRNG, a::AbstractArray, x::AbstractArray)
    s = Sampler(rng, 1:length(a))
    for i = 1:length(x)
        @inbounds x[i] = a[rand(rng, s)]
    end
    return x
end
direct_sample!(a::AbstractArray, x::AbstractArray) = direct_sample!(Random.GLOBAL_RNG, a, x)

### draw a pair of distinct integers in [1:n]

"""
    samplepair([rng], n)

Draw a pair of distinct integers between 1 and `n` without replacement.

Optionally specify a random number generator `rng` as the first argument
(defaults to `Random.GLOBAL_RNG`).
"""
function samplepair(rng::AbstractRNG, n::Int)
    i1 = rand(rng, 1:n)
    i2 = rand(rng, 1:n-1)
    return (i1, ifelse(i2 == i1, n, i2))
end
samplepair(n::Int) = samplepair(Random.GLOBAL_RNG, n)

"""
    samplepair([rng], a)

Draw a pair of distinct elements from the array `a` without replacement.

Optionally specify a random number generator `rng` as the first argument
(defaults to `Random.GLOBAL_RNG`).
"""
function samplepair(rng::AbstractRNG, a::AbstractArray)
    i1, i2 = samplepair(rng, length(a))
    return a[i1], a[i2]
end
samplepair(a::AbstractArray) = samplepair(Random.GLOBAL_RNG, a)

### Algorithm for sampling without replacement

"""
    knuths_sample!([rng], a, x)

*Knuth's Algorithm S* for random sampling without replacement.

Reference: D. Knuth. *The Art of Computer Programming*. Vol 2, 3.4.2, p.142.

This algorithm consumes `length(a)` random numbers. It requires no additional
memory space. Suitable for the case where memory is tight.
"""
function knuths_sample!(rng::AbstractRNG, a::AbstractArray, x::AbstractArray;
                        initshuffle::Bool=true)
    n = length(a)
    k = length(x)
    k <= n || error("length(x) should not exceed length(a)")

    # initialize
    for i = 1:k
        @inbounds x[i] = a[i]
    end
    if initshuffle
        @inbounds for j = 1:k
            l = rand(rng, j:k)
            if l != j
                t = x[j]
                x[j] = x[l]
                x[l] = t
            end
        end
    end

    # scan remaining
    s = Sampler(rng, 1:k)
    for i = k+1:n
        if rand(rng) * i < k  # keep it with probability k / i
            @inbounds x[rand(rng, s)] = a[i]
        end
    end
    return x
end
knuths_sample!(a::AbstractArray, x::AbstractArray; initshuffle::Bool=true) =
    knuths_sample!(Random.GLOBAL_RNG, a, x; initshuffle=initshuffle)

"""
    fisher_yates_sample!([rng], a::AbstractArray, x::AbstractArray)

Fisher-Yates shuffling (with early termination).

Pseudo-code:
```
n = length(a)
k = length(x)

# Create an array of the indices
inds = collect(1:n)

for i = 1:k
    # swap element `i` with another random element in inds[i:n]
    # set element `i` in `x`
end
```

This algorithm consumes `k=length(x)` random numbers. It uses an integer array of
length `n=length(a)` internally to maintain the shuffled indices. It is considerably
faster than Knuth's algorithm especially when `n` is greater than `k`.
It is ``O(n)`` for initialization, plus ``O(k)`` for random shuffling
"""
function fisher_yates_sample!(rng::AbstractRNG, a::AbstractArray, x::AbstractArray)
    n = length(a)
    k = length(x)
    k <= n || error("length(x) should not exceed length(a)")

    inds = Vector{Int}(undef, n)
    for i = 1:n
        @inbounds inds[i] = i
    end

    @inbounds for i = 1:k
        j = rand(rng, i:n)
        t = inds[j]
        inds[j] = inds[i]
        inds[i] = t
        x[i] = a[t]
    end
    return x
end
fisher_yates_sample!(a::AbstractArray, x::AbstractArray) =
    fisher_yates_sample!(Random.GLOBAL_RNG, a, x)

"""
    self_avoid_sample!([rng], a::AbstractArray, x::AbstractArray)

Self-avoid sampling: use a set to maintain the index that has been sampled.
Each time draw a new index, if the index has already been sampled,
redraw until it draws an unsampled one.

This algorithm consumes about (or slightly more than) `k=length(x)` random numbers,
and requires ``O(k)`` memory to store the set of sampled indices.
Very fast when ``n >> k``, with `n=length(a)`.

However, if `k` is large and approaches ``n``, the rejection rate would increase
drastically, resulting in poorer performance.
"""
function self_avoid_sample!(rng::AbstractRNG, a::AbstractArray, x::AbstractArray)
    n = length(a)
    k = length(x)
    k <= n || error("length(x) should not exceed length(a)")

    s = Set{Int}()
    sizehint!(s, k)
    rgen = Sampler(rng, 1:n)

    # first one
    idx = rand(rng, rgen)
    x[1] = a[idx]
    push!(s, idx)

    # remaining
    for i = 2:k
        idx = rand(rng, rgen)
        while idx in s
            idx = rand(rng, rgen)
        end
        x[i] = a[idx]
        push!(s, idx)
    end
    return x
end
self_avoid_sample!(a::AbstractArray, x::AbstractArray) =
    self_avoid_sample!(Random.GLOBAL_RNG, a, x)

"""
    seqsample_a!([rng], a::AbstractArray, x::AbstractArray)

Random subsequence sampling using algorithm A described in the following paper (page 714):
Jeffrey Scott Vitter. "Faster Methods for Random Sampling". Communications of the ACM,
27 (7), July 1984.

This algorithm consumes ``O(n)`` random numbers, with `n=length(a)`.
The outputs are ordered.
"""
function seqsample_a!(rng::AbstractRNG, a::AbstractArray, x::AbstractArray)
    n = length(a)
    k = length(x)
    k <= n || error("length(x) should not exceed length(a)")

    i = 0
    j = 0
    while k > 1
        u = rand(rng)
        q = (n - k) / n
        while q > u  # skip
            i += 1
            n -= 1
            q *= (n - k) / n
        end
        @inbounds x[j+=1] = a[i+=1]
        n -= 1
        k -= 1
    end

    if k > 0  # checking k > 0 is necessary: x can be empty
        s = trunc(Int, n * rand(rng))
        x[j+1] = a[i+(s+1)]
    end
    return x
end
seqsample_a!(a::AbstractArray, x::AbstractArray) = seqsample_a!(Random.GLOBAL_RNG, a, x)

"""
    seqsample_c!([rng], a::AbstractArray, x::AbstractArray)

Random subsequence sampling using algorithm C described in the following paper (page 715):
Jeffrey Scott Vitter. "Faster Methods for Random Sampling". Communications of the ACM,
27 (7), July 1984.

This algorithm consumes ``O(k^2)`` random numbers, with `k=length(x)`.
The outputs are ordered.
"""
function seqsample_c!(rng::AbstractRNG, a::AbstractArray, x::AbstractArray)
    n = length(a)
    k = length(x)
    k <= n || error("length(x) should not exceed length(a)")

    i = 0
    j = 0
    while k > 1
        l = n - k + 1
        minv = l
        u = n
        while u >= l
            v = u * rand(rng)
            if v < minv
                minv = v
            end
            u -= 1
        end
        s = trunc(Int, minv) + 1
        x[j+=1] = a[i+=s]
        n -= s
        k -= 1
    end

    if k > 0
        s = trunc(Int, n * rand(rng))
        x[j+1] = a[i+(s+1)]
    end
    return x
end
seqsample_c!(a::AbstractArray, x::AbstractArray) = seqsample_c!(Random.GLOBAL_RNG, a, x)

"""
    seqsample_d!([rng], a::AbstractArray, x::AbstractArray)

Random subsequence sampling using algorithm D described in the following paper (page 716-17):
Jeffrey Scott Vitter. "Faster Methods for Random Sampling". Communications of the ACM,
27 (7), July 1984.

This algorithm consumes ``O(k)`` random numbers, with `k=length(x)`.
The outputs are ordered.
"""
function seqsample_d!(rng::AbstractRNG, a::AbstractArray, x::AbstractArray)
    N = length(a)
    n = length(x)
    n <= N || error("length(x) should not exceed length(a)")

    i = 0
    j = 0

    vprime = exp(-randexp(rng)/n)
    q1 = N - n + 1
    q2 = q1 / N
    alpha = 1 / 13 # choose alpha value
    threshold = alpha * n

    while n > 1 && threshold < N
        while true
           local X
            while true
                X = N * (1 - vprime)
                s = trunc(Int, X)
                if s < q1
                    break
                end
                vprime = exp(-randexp(rng)/n)
            end

            y = rand(rng) / q2
            lhs = exp(log(y) / (n - 1))
            rhs = ((q1 - s) / q1) * (N / (N - X))

            if lhs <= rhs
                vprime = lhs / rhs
                break
            end

            if n - 1 > s
                bottom = N - n
                limit = N - s
            else
                bottom = N - s - 1
                limit = q1
            end

            top = N - 1

            while top >= limit
                y = y * top / bottom
                bottom -= 1
                top -= 1
            end

            if log(y) < (n - 1)*(log(N) - log(N - X))
                vprime = exp(-randexp(rng) / (n-1))
                break
            end
            vprime = exp(-randexp(rng)/n)
        end

        j += 1
        i += s+1
        @inbounds x[j] = a[i]
        N = N - s - 1
        n -= 1
        q1 -= s
        q2 = q1 / N
        threshold -= alpha
    end

    if n > 1
        seqsample_a!(rng, a[i+1:end], @view x[j+1:end])
    else
        s = trunc(Int, N * vprime)
        @inbounds x[j+=1] = a[i+=s+1]
    end
end

seqsample_d!(a::AbstractArray, x::AbstractArray) = seqsample_d!(Random.GLOBAL_RNG, a, x)

### Interface functions (poly-algorithms)
"""
    sample([rng], a, [wv::AbstractWeights])

Select a single random element of `a`. Sampling probabilities are proportional to
the weights given in `wv`, if provided.

Optionally specify a random number generator `rng` as the first argument
(defaults to `Random.GLOBAL_RNG`).
"""
sample(rng::AbstractRNG, a::AbstractArray) = a[rand(rng, 1:length(a))]
sample(a::AbstractArray) = sample(Random.GLOBAL_RNG, a)

"""
    sample!([rng], a, [wv::AbstractWeights], x; replace=true, ordered=false)

Draw a random sample of `length(x)` elements from an array `a`
and store the result in `x`. A polyalgorithm is used for sampling.
Sampling probabilities are proportional to the weights given in `wv`,
if provided. `replace` dictates whether sampling is performed with
replacement and `order` dictates whether an ordered sample, also called
a sequential sample, should be taken.

Optionally specify a random number generator `rng` as the first argument
(defaults to `Random.GLOBAL_RNG`).
"""
function sample!(rng::AbstractRNG, a::AbstractArray, x::AbstractArray;
                 replace::Bool=true, ordered::Bool=false)
    n = length(a)
    k = length(x)
    k == 0 && return x

    if replace  # with replacement
        if ordered
            sort!(direct_sample!(rng, a, x))
        else
            direct_sample!(rng, a, x)
        end

    else  # without replacement
        k <= n || error("Cannot draw more samples without replacement.")

        if ordered
            if n > 10 * k * k
                seqsample_c!(rng, a, x)
            else
                seqsample_a!(rng, a, x)
            end
        else
            if k == 1
                @inbounds x[1] = sample(rng, a)
            elseif k == 2
                @inbounds (x[1], x[2]) = samplepair(rng, a)
            elseif n < k * 24
                fisher_yates_sample!(rng, a, x)
            else
                self_avoid_sample!(rng, a, x)
            end
        end
    end
    return x
end
sample!(a::AbstractArray, x::AbstractArray; replace::Bool=true, ordered::Bool=false) =
    sample!(Random.GLOBAL_RNG, a, x; replace=replace, ordered=ordered)

"""
    sample([rng], a, [wv::AbstractWeights], n::Integer; replace=true, ordered=false)

Select a random, optionally weighted sample of size `n` from an array `a`
using a polyalgorithm. Sampling probabilities are proportional to the weights
given in `wv`, if provided. `replace` dictates whether sampling is performed
with replacement and `order` dictates whether an ordered sample, also called
a sequential sample, should be taken.

Optionally specify a random number generator `rng` as the first argument
(defaults to `Random.GLOBAL_RNG`).
"""
function sample(rng::AbstractRNG, a::AbstractArray{T}, n::Integer;
                replace::Bool=true, ordered::Bool=false) where T
    sample!(rng, a, Vector{T}(undef, n); replace=replace, ordered=ordered)
end
sample(a::AbstractArray, n::Integer; replace::Bool=true, ordered::Bool=false) =
    sample(Random.GLOBAL_RNG, a, n; replace=replace, ordered=ordered)

"""
    sample([rng], a, [wv::AbstractWeights], dims::Dims; replace=true, ordered=false)

Select a random, optionally weighted sample from an array `a` specifying
the dimensions `dims` of the output array. Sampling probabilities are
proportional to the weights given in `wv`, if provided. `replace` dictates
whether sampling is performed with replacement and `order` dictates whether
an ordered sample, also called a sequential sample, should be taken.

Optionally specify a random number generator `rng` as the first argument
(defaults to `Random.GLOBAL_RNG`).
"""
function sample(rng::AbstractRNG, a::AbstractArray{T}, dims::Dims;
                replace::Bool=true, ordered::Bool=false) where T
    sample!(rng, a, Array{T}(undef, dims); replace=replace, ordered=ordered)
end
sample(a::AbstractArray, dims::Dims; replace::Bool=true, ordered::Bool=false) =
    sample(Random.GLOBAL_RNG, a, dims; replace=replace, ordered=ordered)

################################################################
#
#  Weighted sampling
#
################################################################

"""
    sample([rng], wv::AbstractWeights)

Select a single random integer in `1:length(wv)` with probabilities
proportional to the weights given in `wv`.

Optionally specify a random number generator `rng` as the first argument
(defaults to `Random.GLOBAL_RNG`).
"""
function sample(rng::AbstractRNG, wv::AbstractWeights)
    t = rand(rng) * sum(wv)
    n = length(wv)
    i = 1
    cw = wv[1]
    while cw < t && i < n
        i += 1
        @inbounds cw += wv[i]
    end
    return i
end
sample(wv::AbstractWeights) = sample(Random.GLOBAL_RNG, wv)

sample(rng::AbstractRNG, a::AbstractArray, wv::AbstractWeights) = a[sample(rng, wv)]
sample(a::AbstractArray, wv::AbstractWeights) = sample(Random.GLOBAL_RNG, a, wv)

"""
    direct_sample!([rng], a::AbstractArray, wv::AbstractWeights, x::AbstractArray)

Direct sampling.

Draw each sample by scanning the weight vector.

Noting `k=length(x)` and `n=length(a)`, this algorithm:
* consumes `k` random numbers
* has time complexity ``O(n k)``, as scanning the weight vector each time takes ``O(n)``
* requires no additional memory space.
"""
function direct_sample!(rng::AbstractRNG, a::AbstractArray,
                        wv::AbstractWeights, x::AbstractArray)
    n = length(a)
    length(wv) == n || throw(DimensionMismatch("Inconsistent lengths."))
    for i = 1:length(x)
        x[i] = a[sample(rng, wv)]
    end
    return x
end
direct_sample!(a::AbstractArray, wv::AbstractWeights, x::AbstractArray) =
    direct_sample!(Random.GLOBAL_RNG, a, wv, x)

function make_alias_table!(w::AbstractVector{Float64}, wsum::Float64,
                           a::AbstractVector{Float64},
                           alias::AbstractVector{Int})
    # Arguments:
    #
    #   w [in]:         input weights
    #   wsum [in]:      pre-computed sum(w)
    #
    #   a [out]:        acceptance probabilities
    #   alias [out]:    alias table
    #
    # Note: a and w can be the same way, then that away will be
    #       overriden inplace by acceptance probabilities
    #
    # Returns nothing
    #

    n = length(w)
    length(a) == length(alias) == n ||
        throw(DimensionMismatch("Inconsistent array lengths."))

    ac = n / wsum
    for i = 1:n
        @inbounds a[i] = w[i] * ac
    end

    larges = Vector{Int}(undef, n)
    smalls = Vector{Int}(undef, n)
    kl = 0  # actual number of larges
    ks = 0  # actual number of smalls

    for i = 1:n
        @inbounds ai = a[i]
        if ai > 1.0
            larges[kl+=1] = i  # push to larges
        elseif ai < 1.0
            smalls[ks+=1] = i  # push to smalls
        end
    end

    while kl > 0 && ks > 0
        s = smalls[ks]; ks -= 1  # pop from smalls
        l = larges[kl]; kl -= 1  # pop from larges
        @inbounds alias[s] = l
        @inbounds al = a[l] = (a[l] - 1.0) + a[s]
        if al > 1.0
            larges[kl+=1] = l  # push to larges
        else
            smalls[ks+=1] = l  # push to smalls
        end
    end

    # this loop should be redundant, except for rounding
    for i = 1:ks
        @inbounds a[smalls[i]] = 1.0
    end
    nothing
end

"""
    alias_sample!([rng], a::AbstractArray, wv::AbstractWeights, x::AbstractArray)

Alias method.

Build an alias table, and sample therefrom.

Reference: Walker, A. J. "An Efficient Method for Generating Discrete Random Variables
with General Distributions." *ACM Transactions on Mathematical Software* 3 (3): 253, 1977.

Noting `k=length(x)` and `n=length(a)`, this algorithm takes ``O(n \\log n)`` time
for building the alias table, and then ``O(1)`` to draw each sample. It consumes ``2 k`` random numbers.
"""
function alias_sample!(rng::AbstractRNG, a::AbstractArray, wv::AbstractWeights, x::AbstractArray)
    n = length(a)
    length(wv) == n || throw(DimensionMismatch("Inconsistent lengths."))

    # create alias table
    ap = Vector{Float64}(undef, n)
    alias = Vector{Int}(undef, n)
    make_alias_table!(wv, sum(wv), ap, alias)

    # sampling
    s = Sampler(rng, 1:n)
    for i = 1:length(x)
        j = rand(rng, s)
        x[i] = rand(rng) < ap[j] ? a[j] : a[alias[j]]
    end
    return x
end
alias_sample!(a::AbstractArray, wv::AbstractWeights, x::AbstractArray) =
    alias_sample!(Random.GLOBAL_RNG, a, wv, x)

"""
    naive_wsample_norep!([rng], a::AbstractArray, wv::AbstractWeights, x::AbstractArray)

Naive implementation of weighted sampling without replacement.

It makes a copy of the weight vector at initialization, and sets the weight to zero
when the corresponding sample is picked.

Noting `k=length(x)` and `n=length(a)`, this algorithm consumes ``O(k)`` random numbers,
and has overall time complexity ``O(n k)``.
"""
function naive_wsample_norep!(rng::AbstractRNG, a::AbstractArray,
                              wv::AbstractWeights, x::AbstractArray)
    n = length(a)
    length(wv) == n || throw(DimensionMismatch("Inconsistent lengths."))
    k = length(x)

    w = Vector{Float64}(undef, n)
    copyto!(w, wv)
    wsum = sum(wv)

    for i = 1:k
        u = rand(rng) * wsum
        j = 1
        c = w[1]
        while c < u && j < n
            @inbounds c += w[j+=1]
        end
        @inbounds x[i] = a[j]

        @inbounds wsum -= w[j]
        @inbounds w[j] = 0.0
    end
    return x
end
naive_wsample_norep!(a::AbstractArray, wv::AbstractWeights, x::AbstractArray) =
    naive_wsample_norep!(Random.GLOBAL_RNG, a, wv, x)

# Weighted sampling without replacement
# Instead of keys u^(1/w) where u = random(0,1) keys w/v where v = randexp(1) are used.
"""
    efraimidis_a_wsample_norep!([rng], a::AbstractArray, wv::AbstractWeights, x::AbstractArray)

Weighted sampling without replacement using Efraimidis-Spirakis A algorithm.

Reference: Efraimidis, P. S., Spirakis, P. G. "Weighted random sampling with a reservoir."
*Information Processing Letters*, 97 (5), 181-185, 2006. doi:10.1016/j.ipl.2005.11.003.

Noting `k=length(x)` and `n=length(a)`, this algorithm takes ``O(n + k \\log k)``
processing time to draw ``k`` elements. It consumes ``n`` random numbers.
"""
function efraimidis_a_wsample_norep!(rng::AbstractRNG, a::AbstractArray,
                                     wv::AbstractWeights, x::AbstractArray)
    n = length(a)
    length(wv) == n || throw(DimensionMismatch("a and wv must be of same length (got $n and $(length(wv)))."))
    k = length(x)

    # calculate keys for all items
    keys = randexp(rng, n)
    for i in 1:n
        @inbounds keys[i] = wv.values[i]/keys[i]
    end

    # return items with largest keys
    index = sortperm(keys; alg = PartialQuickSort(k), rev = true)
    for i in 1:k
        @inbounds x[i] = a[index[i]]
    end
    return x
end
efraimidis_a_wsample_norep!(a::AbstractArray, wv::AbstractWeights, x::AbstractArray) =
    efraimidis_a_wsample_norep!(Random.GLOBAL_RNG, a, wv, x)

# Weighted sampling without replacement
# Instead of keys u^(1/w) where u = random(0,1) keys w/v where v = randexp(1) are used.
"""
    efraimidis_ares_wsample_norep!([rng], a::AbstractArray, wv::AbstractWeights, x::AbstractArray)

Implementation of weighted sampling without replacement using Efraimidis-Spirakis A-Res algorithm.

Reference: Efraimidis, P. S., Spirakis, P. G. "Weighted random sampling with a reservoir."
*Information Processing Letters*, 97 (5), 181-185, 2006. doi:10.1016/j.ipl.2005.11.003.

Noting `k=length(x)` and `n=length(a)`, this algorithm takes ``O(k \\log(k) \\log(n / k))``
processing time to draw ``k`` elements. It consumes ``n`` random numbers.
"""
function efraimidis_ares_wsample_norep!(rng::AbstractRNG, a::AbstractArray,
                                        wv::AbstractWeights, x::AbstractArray)
    n = length(a)
    length(wv) == n || throw(DimensionMismatch("a and wv must be of same length (got $n and $(length(wv)))."))
    k = length(x)
    k > 0 || return x

    # initialize priority queue
    pq = Vector{Pair{Float64,Int}}(undef, k)
    i = 0
    s = 0
    @inbounds for _s in 1:n
        s = _s
        w = wv.values[s]
        w < 0 && error("Negative weight found in weight vector at index $s")
        if w > 0
            i += 1
            pq[i] = (w/randexp(rng) => s)
        end
        i >= k && break
    end
    i < k && throw(DimensionMismatch("wv must have at least $k strictly positive entries (got $i)"))
    heapify!(pq)

    # set threshold
    @inbounds threshold = pq[1].first

    @inbounds for i in s+1:n
        w = wv.values[i]
        w < 0 && error("Negative weight found in weight vector at index $i")
        w > 0 || continue
        key = w/randexp(rng)

        # if key is larger than the threshold
        if key > threshold
            # update priority queue
            pq[1] = (key => i)
            percolate_down!(pq, 1)

            # update threshold
            threshold = pq[1].first
        end
    end

    # fill output array with items in descending order
    @inbounds for i in k:-1:1
        x[i] = a[heappop!(pq).second]
    end
    return x
end
efraimidis_ares_wsample_norep!(a::AbstractArray, wv::AbstractWeights, x::AbstractArray) =
    efraimidis_ares_wsample_norep!(Random.GLOBAL_RNG, a, wv, x)

# Weighted sampling without replacement
# Instead of keys u^(1/w) where u = random(0,1) keys w/v where v = randexp(1) are used.
"""
    efraimidis_aexpj_wsample_norep!([rng], a::AbstractArray, wv::AbstractWeights, x::AbstractArray)

Implementation of weighted sampling without replacement using Efraimidis-Spirakis A-ExpJ algorithm.

Reference: Efraimidis, P. S., Spirakis, P. G. "Weighted random sampling with a reservoir."
*Information Processing Letters*, 97 (5), 181-185, 2006. doi:10.1016/j.ipl.2005.11.003.

Noting `k=length(x)` and `n=length(a)`, this algorithm takes ``O(k \\log(k) \\log(n / k))``
processing time to draw ``k`` elements. It consumes ``O(k \\log(n / k))`` random numbers.
"""
function efraimidis_aexpj_wsample_norep!(rng::AbstractRNG, a::AbstractArray,
                                         wv::AbstractWeights, x::AbstractArray)
    n = length(a)
    length(wv) == n || throw(DimensionMismatch("a and wv must be of same length (got $n and $(length(wv)))."))
    k = length(x)
    k > 0 || return x

    # initialize priority queue
    pq = Vector{Pair{Float64,Int}}(undef, k)
    i = 0
    s = 0
    @inbounds for _s in 1:n
        s = _s
        w = wv.values[s]
        w < 0 && error("Negative weight found in weight vector at index $s")
        if w > 0
            i += 1
            pq[i] = (w/randexp(rng) => s)
        end
        i >= k && break
    end
    i < k && throw(DimensionMismatch("wv must have at least $k strictly positive entries (got $i)"))
    heapify!(pq)

    # set threshold
    @inbounds threshold = pq[1].first
    X = threshold*randexp(rng)

    @inbounds for i in s+1:n
        w = wv.values[i]
        w < 0 && error("Negative weight found in weight vector at index $i")
        w > 0 || continue
        X -= w
        X <= 0 || continue

        # update priority queue
        t = exp(-w/threshold)
        pq[1] = (-w/log(t+rand(rng)*(1-t)) => i)
        percolate_down!(pq, 1)

        # update threshold
        threshold = pq[1].first
        X = threshold * randexp(rng)
    end

    # fill output array with items in descending order
    @inbounds for i in k:-1:1
        x[i] = a[heappop!(pq).second]
    end
    return x
end
efraimidis_aexpj_wsample_norep!(a::AbstractArray, wv::AbstractWeights, x::AbstractArray) =
    efraimidis_aexpj_wsample_norep!(Random.GLOBAL_RNG, a, wv, x)

function sample!(rng::AbstractRNG, a::AbstractArray, wv::AbstractWeights, x::AbstractArray;
                 replace::Bool=true, ordered::Bool=false)
    n = length(a)
    k = length(x)

    if replace
        if ordered
            sort!(direct_sample!(rng, a, wv, x))
        else
            if n < 40
                direct_sample!(rng, a, wv, x)
            else
                t = ifelse(n < 500, 64, 32)
                if k < t
                    direct_sample!(rng, a, wv, x)
                else
                    alias_sample!(rng, a, wv, x)
                end
            end
        end
    else
        k <= n || error("Cannot draw $n samples from $k samples without replacement.")

        efraimidis_aexpj_wsample_norep!(rng, a, wv, x)
        if ordered
            sort!(x)
        end
    end
    return x
end
sample!(a::AbstractArray, wv::AbstractWeights, x::AbstractArray) =
    sample!(Random.GLOBAL_RNG, a, wv, x)

sample(rng::AbstractRNG, a::AbstractArray{T}, wv::AbstractWeights, n::Integer;
       replace::Bool=true, ordered::Bool=false) where {T} =
    sample!(rng, a, wv, Vector{T}(undef, n); replace=replace, ordered=ordered)
sample(a::AbstractArray, wv::AbstractWeights, n::Integer;
       replace::Bool=true, ordered::Bool=false) =
    sample(Random.GLOBAL_RNG, a, wv, n; replace=replace, ordered=ordered)

sample(rng::AbstractRNG, a::AbstractArray{T}, wv::AbstractWeights, dims::Dims;
       replace::Bool=true, ordered::Bool=false) where {T} =
    sample!(rng, a, wv, Array{T}(undef, dims); replace=replace, ordered=ordered)
sample(a::AbstractArray, wv::AbstractWeights, dims::Dims;
       replace::Bool=true, ordered::Bool=false) =
    sample(Random.GLOBAL_RNG, a, wv, dims; replace=replace, ordered=ordered)

# wsample interface

"""
    wsample!([rng], a, w, x; replace=true, ordered=false)

Select a weighted sample from an array `a` and store the result in `x`. Sampling
probabilities are proportional to the weights given in `w`. `replace` dictates
whether sampling is performed with replacement and `order` dictates whether an
ordered sample, also called a sequential sample, should be taken.

Optionally specify a random number generator `rng` as the first argument
(defaults to `Random.GLOBAL_RNG`).
"""
wsample!(rng::AbstractRNG, a::AbstractArray, w::RealVector, x::AbstractArray;
         replace::Bool=true, ordered::Bool=false) =
    sample!(rng, a, weights(w), x; replace=replace, ordered=ordered)
wsample!(a::AbstractArray, w::RealVector, x::AbstractArray;
         replace::Bool=true, ordered::Bool=false) =
    sample!(Random.GLOBAL_RNG, a, weights(w), x; replace=replace, ordered=ordered)

"""
    wsample([rng], [a], w)

Select a weighted random sample of size 1 from `a` with probabilities proportional
to the weights given in `w`. If `a` is not present, select a random weight from `w`.

Optionally specify a random number generator `rng` as the first argument
(defaults to `Random.GLOBAL_RNG`).
"""
wsample(rng::AbstractRNG, w::RealVector) = sample(rng, weights(w))
wsample(w::RealVector) = wsample(Random.GLOBAL_RNG, w)
wsample(rng::AbstractRNG, a::AbstractArray, w::RealVector) = sample(rng, a, weights(w))
wsample(a::AbstractArray, w::RealVector) = wsample(Random.GLOBAL_RNG, a, w)

"""
    wsample([rng], [a], w, n::Integer; replace=true, ordered=false)

Select a weighted random sample of size `n` from `a` with probabilities proportional
to the weights given in `w` if `a` is present, otherwise select a random sample of size
`n` of the weights given in `w`. `replace` dictates whether sampling is performed with
replacement and `order` dictates whether an ordered sample, also called a sequential
sample, should be taken.

Optionally specify a random number generator `rng` as the first argument
(defaults to `Random.GLOBAL_RNG`).
"""
wsample(rng::AbstractRNG, a::AbstractArray{T}, w::RealVector, n::Integer;
        replace::Bool=true, ordered::Bool=false) where {T} =
    wsample!(rng, a, w, Vector{T}(undef, n); replace=replace, ordered=ordered)
wsample(a::AbstractArray, w::RealVector, n::Integer;
        replace::Bool=true, ordered::Bool=false) =
    wsample(Random.GLOBAL_RNG, a, w, n; replace=replace, ordered=ordered)

"""
    wsample([rng], [a], w, dims::Dims; replace=true, ordered=false)

Select a weighted random sample from `a` with probabilities proportional to the
weights given in `w` if `a` is present, otherwise select a random sample of size
`n` of the weights given in `w`. The dimensions of the output are given by `dims`.

Optionally specify a random number generator `rng` as the first argument
(defaults to `Random.GLOBAL_RNG`).
"""
wsample(rng::AbstractRNG, a::AbstractArray{T}, w::RealVector, dims::Dims;
        replace::Bool=true, ordered::Bool=false) where {T} =
    wsample!(rng, a, w, Array{T}(undef, dims); replace=replace, ordered=ordered)
wsample(a::AbstractArray, w::RealVector, dims::Dims;
        replace::Bool=true, ordered::Bool=false) =
    wsample(Random.GLOBAL_RNG, a, w, dims; replace=replace, ordered=ordered)


text:
module LFPSQP

using ForwardDiff
using ForwardDiff: GradientConfig, JacobianConfig, Chunk, Dual, Tag
using ReverseDiff
using LinearAlgebra
using LinearAlgebra.BLAS: gemv!, axpy!, ger!
using Parameters
using Printf
using Random: randn!
using LinearMaps

# parameters
import Base.show

# redifine tags to avoid errors
# NOTE that this may break nested differentation, but I don't believe this can be the case
# based on the way that I use these dual types in newton_direction.jl
@inline function ForwardDiff.:≺(::Type{Tag{F1,V1}}, Nothing) where {F1,V1}
    return false
end

@inline function ForwardDiff.:≺(Nothing, ::Type{Tag{F1,V1}}) where {F1,V1}
    return true
end

@enum DisplayOption begin
    off
    iter
end

@enum LinesearchOption begin
    armijo
    exact
end

@enum TerminationCondition begin
	f_tol
	x_tol
	kkt_tol
	max_iter
	armijo_error
end

struct TerminationInfo
	condition::TerminationCondition
	f_diff::Float64
	step_diff::Float64
	kkt_diff::Float64
	iter::Int64
end

Base.show(io::IO, ti::TerminationInfo) = print(io, "TerminationInfo:\ncondition = $(ti.condition)\n" *
	"       Δf = $(ti.f_diff)\n   ||Δx|| = $(ti.step_diff)\n||P(∇f)|| = $(ti.kkt_diff)\n    iters = $(ti.iter)")

@with_kw struct LFPSQPParams
	α::Float64 = 1.0
	β::Float64 = 0.0
	t_β::Int64 = 0
	s::Float64 = 0.5
	σ::Float64 = 1e-4
	ϵ_c::Float64 = 1e-6
	ϵ_f::Float64 = 1e-6
	ϵ_x::Float64 = 0.0
	ϵ_kkt::Float64 = 1e-6
	ϵ_rank::Float64 = 1e-10
	maxiter::Int64 = 10000
	maxiter_retract::Int64 = 100
	maxiter_pcg::Int64 = 100
	μ0::Float64 = 1e-2
	disable_linesearch::Bool = false
	do_project_retract::Bool = true
	disp::DisplayOption = iter
	callback::Union{Nothing,Function} = nothing
	callback_period::Int64 = 100
	linesearch::LinesearchOption = armijo
	do_newton::Bool = true
	tn_maxiter::Int64 = 10000
	tn_κ::Float64 = 0.5
end

# include files
include("projcg.jl")
include("autodiff_generators.jl")
include("la_helper.jl")
include("inequality_helper.jl")

include("linesearch.jl")
include("retractions.jl")

include("optimize.jl")

export LFPSQPParams
export optimize

end # module


text:
## integration.jl

"""
riemann: compute Riemann sum approximations to a definite integral. As well, implement trapezoid and Simpson's rule.

Example:
```
f(x) = exp(x^2)
riemann(f, 0, 1, 1000)   # default right-Riemann sums
riemann(f, 0, 1, 1000, method="left")       # left sums
riemann(f, 0, 1, 1000, method="trapezoid")  # use trapezoid rule
riemann(f, 0, 1, 1000, method="simpsons")   # use Simpson's rule
```

"""
function riemann(f::Function, a::Real, b::Real, n::Int; method="right")
    if method == "right"
        meth = f -> (lr -> begin l,r = lr; f(r) * (r-l) end)
    elseif method == "left"
        meth = f -> (lr -> begin l,r = lr; f(l) * (r-l) end)
    elseif method == "trapezoid"
        meth = f -> (lr -> begin l,r = lr; (1/2) * (f(l) + f(r)) * (r-l) end)
    elseif method == "simpsons"
        meth = f -> (lr -> begin l,r=lr; (1/6) * (f(l) + 4*(f((l+r)/2)) + f(r)) * (r-l) end)
    end

    xs = a .+ (0:n) * (b-a)/n

    sum(meth(f), zip(xs[1:end-1], xs[2:end]))
end

#######
## simplified multivariable integrals

# limits of integration
endpoints(ys,x) = ((f,x) -> isa(f, Function) ? f(x...) : f).(ys, Ref(x))
# avoid specialization in quadgk
struct FWrapper
    f
end
(F::FWrapper)(x) = F.f(x)

"""
fubini(f, dy, dx)
fubini(f, dz, dy, dx)

Computes numeric integral of `f` over region specified by `dz`, `dy`, `dx`. These are a tuple of values of numbers or univariate functions depending on the value of the term on the right (`dy` can depend on `dx`s value).

*Much* slower than `hcubature` from the `HCubature` package, as it refines flat areas too many times, allocates too much, etc. But does allow a more flexible specification of the region to integrate over, as `hcubature` requires box-like regions.

```
f(x,y,z) = x * y^2 * z^3
fubini(f, (0,1), (0,2), (0,3))  # int_0^3 int_0^2 int_0^1 f(x,y,z) dz dy dx
g(v) = f(v...)
hcubature(g, (0,0,0), (3,2,1))  # same. Not order switched

# triangular like region
fubini(f, (0, y->y), (0, x->x), (0,3))
```
"""
fubini(f, dx)     = quadgk(FWrapper(f), dx...)[1]
fubini(f, ys, xs) = fubini(x -> fubini(y -> f(x,y), endpoints(ys, x)), xs)
fubini(f, zs, ys, xs) = fubini(x ->
    fubini(y ->
        fubini(z -> f(x,y,z),
            endpoints(zs, (x,y))),
        endpoints(ys,x)),
    xs)


text:
# include("Walk.jl")

abstract Simulation

type MetropolisSimulation <: Simulation
    walk::Walk

    accepted_steps::Int
    rejected_steps::Int
    # automatically_rejected_steps::Int # FIXME: add this
    trivial_steps::Int
    total_steps::Int

    # beeps::Int # FIXME: add this to keep track of numerical error rate

    measurements::Dict{ASCIIString, Measurement}

    function MetropolisSimulation(walk::Walk)
        return new(walk, 0, 0, 0, 0, Dict{ASCIIString, Measurement}())
    end
end

function add_accepted_step!(sim::MetropolisSimulation)
    sim.accepted_steps += 1
    sim.total_steps += 1
end

function add_rejected_step!(sim::MetropolisSimulation)
    sim.rejected_steps += 1
    sim.total_steps += 1
end

# FIXME: add this
# function add_automatically_rejected_step(sim::MetropolisSimulation)
#     sim.rejected_steps += 1
#     sim.automatically_rejected_steps +=1
#     sim.total_steps += 1
# end

function add_trivial_step!(sim::MetropolisSimulation)
    sim.trivial_steps += 1
    sim.total_steps += 1
end

function consistent_MCinfo(sim::MetropolisSimulation)
    return sim.accepted_steps + sim.rejected_steps + sim.trivial_steps == sim.total_steps
end

# FIXME: instead, just report one acceptance rate, with separate functions giving trival and auto reject %s
function acceptance_rate_1(sim::MetropolisSimulation)
    @assert consistent_MCinfo(sim)
    return (sim.accepted_steps + sim.trivial_steps) / sim.total_steps
end

function acceptance_rate_2(sim::MetropolisSimulation)
    @assert consistent_MCinfo(sim)
    return sim.accepted_steps / (sim.accepted_steps + sim.rejected_steps)
end

function acceptance_rate_3(sim::MetropolisSimulation)
    @assert consistent_MCinfo(sim)
    return sim.accepted_steps / sim.total_steps
end

# by default this will do nothing, but it can be overwritten by a given wavefunction implementation
function print_wf_info(wf)
end

function iterate!(sim::MetropolisSimulation, steps_to_iterate::Int, steps_per_check::Int = 100, print_rate::Float64 = 0.05)

    printlnf("Beginning to iterate a Metropolis simulation $(steps_to_iterate) steps!")
    for k in 1:steps_to_iterate

        if mod(k, print_rate*steps_to_iterate) == 0
        # if true
            printlnf("Starting step #$(k)/$(steps_to_iterate) --> $(100k/steps_to_iterate)%")
            print_wf_info(sim.walk.wf)
        end

        # check for error periodically
        if mod(sim.total_steps, steps_per_check) == 0
            # printlnf("Checking for numerical error!")
            check_for_numerical_error!(sim.walk)
        end

        propose_and_assess_random_transition!(sim.walk)
        if !is_nontrivial_transition(sim.walk)
            add_trivial_step!(sim)
            # printlnf("Step is trivial!")
            # printlnf()
            # FIXME: for trivial transitions, we shouldn't have to recompute the measurements!
            for k in keys(sim.measurements)
                possibly_measure(sim.walk.wf, sim.measurements[k].data, sim.measurements[k].measure_function)
            end
            continue
        else
            probrat = get_probability_ratio(sim.walk)
        end
        # printlnf("probrat = $(probrat)")

        # Metropolis!
        if probrat >= 1.0 || rand() < probrat # accept transition
            accept_transition!(sim.walk)
            add_accepted_step!(sim)
            # printlnf("Step accepted!")
        else # reject transition
            reject_transition!(sim.walk)
            add_rejected_step!(sim)
            # printlnf("Step rejected!")
        end
        # printlnf()

        for k in keys(sim.measurements)
            possibly_measure(sim.walk.wf, sim.measurements[k].data, sim.measurements[k].measure_function)
        end

    end
end


text:
# unit.jl

using Unitful
using MIRTjim: jim
using Test: @test, @testset, @inferred

include("isplot.jl")

@testset "unit" begin
    x = 2u"mm" * (1:7)
    y = 2u"s" * (1:5)
    f = x * y'

#   @test @inferred jim(x, y, f) # Plots.Plot{Plots.GRBackend} vs Plot

    r = f / oneunit(f[1]) / 100 # real
    p1 = jim(1:7, 1:5, r)
    @isplot p1
    p2 = jim(1:7, 1:5, f)
    @isplot p2
    p3 = jim(x, y, r)
    @isplot p3
    p4 = jim(x, y, f)
    @isplot p4
    p5 = jim(1:7, 1:5, f; xlabel="x", ylabel="y")
    @isplot p5
    p6 = jim(x, y, f; xlabel="x", ylabel="y", clim=(50,100).*u"mm*s")
    @isplot p6
    jim(p1, p2, p3, p4, p5, p6)

    z = 2 .+ rand(5,5) + 4im * ones(5,5)
    @test z isa AbstractMatrix{<:Complex}
    zu = z * Unitful.m # complex array with units
    @test !(zu isa AbstractMatrix{<:Complex}) # beware!
    @test zu isa AbstractMatrix{<:Number}
    @isplot jim(real(zu), "units") # colorbar roughly 2-3
    @isplot jim(zu, "units") # colorbar roughly 4-5

    z3u = (randn(3,4,5) + ones(3,4,5)) * Unitful.m # 3D, complex, units
    @isplot jim(z3u)
end


text:
module Helper
export @implicit,
       unzip, remove,
       linesof, exprmodify, clean_code,
       Looping, next!, current, restart!

#-----------------------------------------------------------------------------------
# Looping
#-----------------------------------------------------------------------------------

mutable struct Looping
  iter
  state

  Looping(iter) = new(iter, start(iter))
end

current(l::Looping) =
  next(l.iter, l.state)[1]

restart!(l::Looping) =
  l.state = start(l.iter)

function next!(l::Looping)
  l.state = next(l.iter, l.state)[2]

  done(l.iter, l.state) && restart!(l)
end

#-----------------------------------------------------------------------------------
# Iterable utilities
#-----------------------------------------------------------------------------------

function unzip(A)
  F = []; S = []
  for (first, second) in A
    push!(F, first)
    push!(S, second)
  end
  return F,S
end

function remove(f, A)
  B = []
  for x in A
    !f(x) && push!(B, x)
  end
  return B
end

#-----------------------------------------------------------------------------------
# Metaprogramming utilities
#-----------------------------------------------------------------------------------

is_line_number(ex::Expr) = ex.head == :line
is_line_number(ex::LineNumberNode) = true
is_line_number(ex) = false

function linesof(ex::Expr)
  remove(is_line_number, ex.args)
end

function clean_code(ex::Expr)
  nex = Expr(ex.head, map(clean_code, ex.args)...)

  nex.head == :block && return Expr(:block, linesof(nex)...)
  nex.head == :(->)  && return Expr(:(->),  nex.args[1], nex.args[2].args[1])
  nex
end
clean_code(x) = x

end


text:
import DoubleFloats: Double64, ComplexDF64

@testset "LinearAlgebra" begin
    @testset "Constructor + basics" begin
        A = rand(ComplexF32, 12, 12)
        WS = HC.MatrixWorkspace(A)
        @test WS isa HC.MatrixWorkspace{Float64}
        @test_throws ArgumentError HC.update!(WS, rand(12, 13))

        WS.factorized[] = true
        B = randn(ComplexF64, 12, 12)
        HC.update!(WS, B)
        @test WS.A == B
        @test WS.A !== B
        @test WS.factorized[] == false

        HC.factorization!(WS, HC.LU_FACT)
        @test WS.fact[] == HC.LU_FACT
        WS.factorized[] = true
        HC.factorization!(WS, HC.QR_FACT)
        @test WS.fact[] == HC.QR_FACT
        @test WS.factorized[] == false

        HC.factorization!(WS, HC.LU_FACT)
        @test (@allocated HC.factorize!(WS)) == 0
        HC.factorization!(WS, HC.QR_FACT)
        @test WS.factorized[] == false
        HC.factorize!(WS) # first factorize resizes the work buffers
        WS.factorized[] = false
        @test (@allocated HC.factorize!(WS)) == 0

# test indexing etc
        @test WS[2, 1] == B[2, 1]
        @test WS[2] == B[2]
        WS[2, 3] = 5.0
        @test WS[2, 3] == 5.0
        WS[5] = 2.32
        @test WS[5] == 2.32
        @test size(WS) == (12, 12)
        C = randn(ComplexF64, 12, 12)
        WS .= C
        @test WS.A == C
    end

    @testset "ldiv" begin
        A = randn(ComplexF64, 13, 13)
        b = randn(ComplexF64, 13)
        x = zeros(ComplexF64, 13)
        WS = HC.MatrixWorkspace(A)
        HC.update!(WS, A)
        ldiv!(x, WS, b)
        WS.factorized[] = false
        @test (@allocated ldiv!(x, WS, b)) == 0
        @test (lu(A) \ b) ≈ x rtol = 1e-12

        HC.factorization!(WS, HC.QR_FACT)
        x .= 0
        ldiv!(x, WS, b)
        @test (qr(A, Val(true)) \ b) ≈ x rtol = 1e-12
        WS.factorized[] = false
        @test (@allocated ldiv!(x, WS, b)) == 0

# overdetermined
        A = randn(ComplexF64, 13, 8)
        x = randn(ComplexF64, 8)
        b = A * x
        WS = HC.MatrixWorkspace(A)
        HC.update!(WS, A)
        @test_throws ArgumentError HC.factorization!(WS, HC.LU_FACT)
        ldiv!(x, WS, b)
        WS.factorized[] = false
        x .= 0
        @test (@allocated ldiv!(x, WS, b)) == 0
        @test (qr(A, Val(true)) \ b) ≈ x rtol = 1e-12
    end

    @testset "(Mixed precision) residual" begin
        A = randn(ComplexF64, 13, 13)
        x = randn(ComplexF64, 13)
        b = A * x
        û = convert.(ComplexF64, ComplexDF64.(A) * ComplexDF64.(x) - ComplexDF64.(b))

        u = zeros(ComplexF64, 13)
# fixed precision residual should be of order unit roundoff
        HC.residual!(u, A, x, b)
        @test maximum(abs, ComplexDF64.(u) - ComplexDF64.(û)) < 1e-14
        @test maximum(abs, ComplexDF64.(u) - ComplexDF64.(û)) > 1e-18

        HC.residual!(u, A, x, b, Complex{Double64})
        @test maximum(abs, ComplexDF64.(u) - ComplexDF64.(û)) < 1e-26
    end

    @testset "(Mixed precision) iterative refinement" begin
        A = randn(ComplexF64, 4, 4)
        b = randn(ComplexF64, 4)
        x̂ = zeros(ComplexF64, 4)
        x = copy(x̂)
        WS = HC.MatrixWorkspace(A)
        U, S, VT = svd(A)
        S[end] *= 1e-4
        A2 = U * diagm(0 => S) * VT'
        HC.update!(WS, A2)
        ldiv!(x̂, WS, b)
        δx1 = HC.iterative_refinement_step!(x, WS, x̂, b, HC.InfNorm(), ComplexDF64)
        @test δx1 / norm(x, Inf) > 1e-14
        δx2 = HC.iterative_refinement_step!(WS, x, b, HC.InfNorm(), ComplexDF64)
        @test δx2 / norm(x, Inf) < eps()

        δx = HC.iterative_refinement_step!(x, WS, x̂, b)
        @test δx1 / norm(x, Inf) > 1e-14

# overdetermined
        A = randn(ComplexF64, 4, 3)
        x̂ = zeros(ComplexF64, 3)
        x = copy(x̂)
        WS = HC.MatrixWorkspace(A)
        U, S, VT = svd(A)
        S[end] *= 1e-4
        A2 = U * diagm(0 => S) * VT'
        HC.update!(WS, A2)
        b = A2 * randn(ComplexF64, 3)
        ldiv!(x̂, WS, b)
        δx1 = HC.iterative_refinement_step!(x, WS, x̂, b, HC.InfNorm(), ComplexDF64)
        @test δx1 / norm(x, Inf) > 1e-14
        δx2 = HC.iterative_refinement_step!(WS, x, b, HC.InfNorm(), ComplexDF64)
        @test δx2 / norm(x, Inf) < eps()

        δx = HC.iterative_refinement_step!(x, WS, x̂, b)
        @test δx1 / norm(x, Inf) > 1e-14
    end

    @testset "rcond" begin
        A = randn(ComplexF64, 13, 13)
        WS = HC.MatrixWorkspace(A)
        U, S, VT = svd(A)
        S[end] *= 1e-5
        A2 = U * diagm(0 => S) * VT'
        HC.update!(WS, A2)
        HC.factorization!(WS, HC.LU_FACT)
        HC.factorize!(WS)
        @test cond(A2, Inf) ≈ inv(HC.rcond(WS))
        HC.factorization!(WS, HC.QR_FACT)
        @test 0.5 ≤ cond(A2, Inf) / inv(HC.rcond(WS)) ≤ 1.5
# @test (@allocated rcond(WS)) == 0

        A = randn(ComplexF64, 13, 10)
        WS = HC.MatrixWorkspace(A)
        U, S, VT = svd(A)
        S[end] *= 1e-8
        A2 = U * diagm(0 => S) * VT'
        HC.update!(WS, A2)
        HC.factorize!(WS)
        @test 1 / HC.rcond(WS) ≈ cond(UpperTriangular(WS.qr.R), Inf) rtol = 1e-14
# @test (@allocated rcond(WS)) == 0

# meddle aroudn with the internals to check dead branch
        WS.fact[] = HC.LU_FACT
        @test isnan(HC.rcond(WS))
    end

    @testset "inf norm est" begin
        A = randn(ComplexF64, 5, 5)
        U, s, Vt = svd(A)
        s[end] *= 1e-6
        B = U * diagm(0 => s) * Vt'

        d = rand(5) * 1000
        D = diagm(0 => d)
        g = rand(5) .* 1000
        G = diagm(0 => g)

        @test 0.1 ≤ opnorm(inv(B), Inf) / HC.inf_norm_est(lu(B)) ≤ 10
        @test 0.1 ≤ opnorm(inv(D) * inv(B), Inf) / HC.inf_norm_est(lu(B), nothing, d) ≤ 10
        @test 0.1 ≤ opnorm(inv(B) * G, Inf) / HC.inf_norm_est(lu(B), g, nothing) ≤ 10
        @test 0.1 ≤ opnorm(inv(D) * inv(B) * G, Inf) / HC.inf_norm_est(lu(B), g, d) ≤ 10

        WS = HC.MatrixWorkspace(B)
        @test 0.1 ≤ opnorm(inv(B), Inf) / HC.inf_norm_est(WS) ≤ 10
        @test 0.1 ≤ opnorm(inv(D) * inv(B), Inf) / HC.inf_norm_est(WS, nothing, d) ≤ 10
        @test 0.1 ≤ opnorm(inv(B) * G, Inf) / HC.inf_norm_est(WS, g, nothing) ≤ 10
        @test 0.1 ≤ opnorm(inv(D) * inv(B) * G, Inf) / HC.inf_norm_est(WS, g, d) ≤ 10
    end

    @testset "JacobianMonitor" begin
        A = randn(ComplexF64, 6, 6)
        b = randn(ComplexF64, 6)
        x̂ = similar(b)
        x = similar(b)
        JM = HC.JacobianMonitor(zeros(ComplexF64, 6, 6))
        @test JM isa HC.JacobianMonitor
        HC.jacobian(JM) .= A
        HC.updated!(JM)
        ldiv!(x̂, JM, b)
        ldiv!(x, jacobian(JM), b)
        @test x == x̂
    end
    @testset "Forward error" begin
        A = [
            1.81451 + 0.46473im -0.473651 + 0.813117im 0.130822 - 0.254704im -0.203277 +
                                                                             0.397106im
            -2.04468 - 1.59228im 0.441157 - 1.39994im -0.0984804 + 1.42619im -0.818292 -
                                                                             0.487898im
            0.35501 + 0.0235574im 1.32618 + 0.577264im 1.12561 + 0.43383im 0.579533 +
                                                                           1.49484im
            0.144381 + 0.104175im 0.665384 - 1.43389im 0.22097 + 0.258016im 0.374103 -
                                                                            0.635066im
        ]
        b = [
            -0.14581 - 0.868847im,
            0.0162454 + 0.91156im,
            -0.133769 + 0.510834im,
            -0.867022 + 0.752196im,
        ]
        x̂ = similar(b)
        x = similar(b)
        JM = HC.JacobianMonitor(A)
        HC.updated!(JM)
        ldiv!(x̂, JM, b)

        ferr = HC.forward_err!(x, JM, x̂, b, HC.InfNorm())
        @test HC.forward_err(JM) == ferr
        @test ferr < eps() * cond(A)
        @test ferr > 1e-12
        ferrD64 = HC.forward_err!(x, JM, x̂, b, HC.InfNorm(), ComplexDF64)
        @test 0.5 ≤ ferrD64 / ferr ≤ 2

        ferr = HC.forward_err!(JM, x̂, b, HC.InfNorm(), ComplexF64)
        @test ferr < eps() * cond(A)
        @test ferr > 1e-12

        D = diagm(0 => [10.0^(2i) for i = 1:4])
        JM = HC.JacobianMonitor(D * A)
        HC.updated!(JM)
        ldiv!(x̂, JM, D * b)
        conderr = HC.cond(D * A) * eps()
        sferr = HC.strong_forward_err!(x, JM, x̂, D * b, HC.InfNorm())
        @test HC.forward_err(JM) == sferr
        ferr = HC.forward_err!(x, JM, x̂, D * b, HC.InfNorm())
        @test ferr < sferr < conderr * 1e-4 # check that error estimate is actually substantially less

        sferr2 = HC.strong_forward_err!(JM, x̂, D * b, HC.InfNorm())
        @test sferr == sferr2
        @test x̂ ≈ x

        DR = diagm(0 => [10.0^(2i - 4) for i = 1:4])
        JM = HC.JacobianMonitor(A * DR)
        ldiv!(x̂, JM, b)

        WN = HC.WeightedNorm(HC.InfNorm(), x̂)
        WN .= abs.(x̂)

        wferr = HC.forward_err!(x, JM, x̂, b, WN)
        ferr = HC.forward_err!(x, JM, x̂, b, HC.InfNorm())

        strong_wferr = HC.strong_forward_err!(x, JM, x̂, b, WN)
        strong_ferr = HC.strong_forward_err!(x, JM, x̂, b, HC.InfNorm())
        @test ferr ≤ strong_ferr ≤ eps() * inv(HC.rcond(jacobian(JM)))
        @test wferr ≤ strong_wferr ≤ eps() * inv(HC.rcond(jacobian(JM)))

# overdetermined
        A = randn(ComplexF64, 3, 2)
        x = randn(ComplexF64, 2)
        b = A * x
        x̂ = similar(x)

        JM = HC.JacobianMonitor(A)
        ldiv!(x̂, JM, b)

        ferr = HC.forward_err!(x, JM, x̂, b, HC.InfNorm())
        strong_ferr = HC.strong_forward_err!(x, JM, x̂, b, HC.InfNorm())
        @test ferr < 1e-13
        @test ferr == strong_ferr # test fallback to ferr for qr
    end

    @testset "cond" begin
        A = randn(ComplexF64, 13, 13)
        D = diagm(0 => [2.0^(i - 5) for i = 1:13])
        JM = HC.JacobianMonitor(D * A)
        comp_cond = HC.cond!(JM)
        @test comp_cond === cond(JM)
        B = diagm(0 => inv.(JM.J.row_scaling)) * D * A
        @test abs(cond(JM) - cond(B, Inf)) / cond(JM) < 2

        A = randn(ComplexF64, 13, 12)
        JM = HC.JacobianMonitor(D * A)
        comp_cond = HC.cond!(JM)
        @test comp_cond === cond(JM)
    end

    @testset "JacobianMonitor ldiv updates" begin
        A = randn(ComplexF64, 6, 6)
        b = randn(ComplexF64, 6)
        x̂ = similar(b)
        x = similar(b)
        JM = HC.JacobianMonitor(zeros(ComplexF64, 6, 6))
        @test JM isa HC.JacobianMonitor
        ldiv!(x, JM, b, HC.InfNorm(), HC.JAC_MONITOR_UPDATE_FERR)
        @test HC.forward_err(JM) != 0.0
        ldiv!(x, JM, b, HC.InfNorm(), HC.JAC_MONITOR_UPDATE_COND)
        @test cond(JM) != 1.0
        JM = HC.JacobianMonitor(zeros(ComplexF64, 6, 6))
        ldiv!(x, JM, b, HC.InfNorm(), HC.JAC_MONITOR_UPDATE_ALL)
        @test cond(JM) != 1.0
        @test HC.forward_err(JM) != 0.0
    end

    @testset "Hermite Normal Form" begin
        A = [0 3 1 0; -2 2 -1 2; 1 -1 2 3; -3 3 3 2]
        H, U = HC.hnf(A)
        @test A * U == H

        D = zeros(Int, 1, 1)
        D[1, 1] = -2
        H, U = HC.hnf(D)
        @test H[1, 1] == 2
        @test U[1, 1] == -1

# overflow A
        A = [
            2 4 4 2 -4 -3 -4 3 0 3
            -4 -2 -5 -5 1 0 -3 0 -1 -2
            -3 5 -3 1 4 3 -1 -4 -1 0
            1 3 3 4 -2 -3 -2 -5 -4 -3
            2 2 0 2 -4 4 3 -4 -2 -4
            -3 -3 1 1 -4 4 0 4 4 -4
            -3 -1 -1 -1 2 -4 -3 -4 4 -4
            0 1 -2 4 5 4 3 1 -5 2
            -5 -4 -5 3 1 5 0 -3 -3 -1
            2 -5 -3 -1 -1 5 -4 5 -3 1
        ]
        @test_throws OverflowError HC.hnf(A)
        @test_throws OverflowError HC.hnf(A, Int128)
        H, U = HC.hnf(A, BigInt)
        @test A * U == H
    end
end