text:
using QuadGK: quadgk
using Roots: find_zero, Order1, ConvergenceFailed, Newton

## Threshold to avoid numerical overflow of the weight function of L1Estimator and ArctanEstimator
DELTA = 1e-8    # chosen because it minimizes the error between (1-ATWDELTA)/DELTA and 1/3 for the linear approximation
L1WDELTA = 1 / (DELTA)
ATWDELTA = atan(DELTA) / DELTA
#DELTA = 1e-8    # chosen because it minimizes the error between (1-ATWDELTA)/DELTA and 1/3 for the linear approximation
#L1WDELTA = 1/(2*sqrt(DELTA))
#ATWDELTA = atan(sqrt(DELTA))*2*L1WDELTA

#=
TODO: change the API to
rlm(form, data, MEstimator(Tukey))
rlm(X, y, TauEstimator(YohaiZamar); method=:cg)
rlm(X, y, MMEstimator(Tukey(2.1), Tukey(5.0)); σ0=:mad)

=#

"The loss function used for the estimation"
function loss end

"The loss function ρ for the M-estimator."
function rho end

"""
The influence function ψ is the derivative of the loss function for the M-estimator,
multiplied by the square of the tuning constant.
"""
function psi end

"The derivative of ψ, used for asymptotic estimates."
function psider end

"The weights for IRLS, the function ψ divided by r."
function weight end

"The integral of exp(-ρ) used for calculating the full-loglikelihood"
function estimator_norm end

"The limit at ∞ of the loss function. Used for scale estimation of bounded loss."
estimator_bound(::Type{<:LossFunction}) = Inf

"The tuning constant of the loss function, can be optimized to get efficient or robust estimates."
tuning_constant(loss::L) where {L<:LossFunction} = loss.c

"Boolean if the estimator or loss function is convex"
isconvex(f::LossFunction) = false
isconvex(f::ConvexLossFunction) = true

"Boolean if the estimator or loss function is bounded"
isbounded(f::LossFunction) = false
isbounded(f::BoundedLossFunction) = true

"The tuning constant associated to the loss that gives an efficient M-estimator."
estimator_high_breakdown_point_constant(::Type{L}) where {L<:LossFunction} = 1

"The tuning constant associated to the loss that gives a robust (high breakdown point) M-estimator."
estimator_high_efficiency_constant(::Type{L}) where {L<:LossFunction} = 1

# Static creators
"The loss initialized with an efficient tuning constant"
efficient_loss(::Type{L}) where {L<:LossFunction} = L(estimator_high_efficiency_constant(L))

"The loss initialized with a robust (high breakdown point) tuning constant"
robust_loss(::Type{L}) where {L<:LossFunction} =
    L(estimator_high_breakdown_point_constant(L))

rho(l::LossFunction, r) = _rho(l, r / tuning_constant(l))
psi(l::LossFunction, r) = tuning_constant(l) * _psi(l, r / tuning_constant(l))
psider(l::LossFunction, r) = _psider(l, r / tuning_constant(l))
weight(l::LossFunction, r) = _weight(l, r / tuning_constant(l))

"Faster version if you need ρ, ψ and w in the same call"
function values(loss::LossFunction, r::Real)
    c = tuning_constant(loss)
    rr = r / c
    return (_rho(loss, rr), c * _psi(loss, rr), _weight(loss, rr))
end

##
# For reminder and not to get lost with the tuning constant
#   ρ(r) ~ r²/(2c²)
#   ψ = c²ρ' ;  ψ ~ r  ;  _psi ~ r/c
#   w ~ 1
##

"""
The M-estimator norm is computed with:
     +∞                    +∞
Z = ∫  exp(-ρ(r))dr = c . ∫  exp(-ρ_1(r))dr    with ρ_1 the function for c=1
    -∞                    -∞
"""
function estimator_norm(loss::L) where {L<:LossFunction}
    2 * quadgk(x -> exp(-rho(loss, x)), 0, Inf)[1]
end

"""
The tuning constant c is computed so the efficiency for Normal distributed
residuals is 0.95. The efficiency of the mean estimate μ is defined by:

eff_μ = (E[ψ'])²/E[ψ²]
"""
function efficiency_tuning_constant(
    ::Type{L};
    eff::Real=0.95,
    c0::Real=1.0,
) where {L<:LossFunction}
    lpsi(x, c) = RobustModels.psi(L(c), x)
    lpsip(x, c) = RobustModels.psider(L(c), x)

    I1(c) = quadgk(x -> (lpsi(x, c))^2 * 2 * exp(-x^2 / 2) / √(2π), 0, Inf)[1]
    I2(c) = quadgk(x -> lpsip(x, c) * 2 * exp(-x^2 / 2) / √(2π), 0, Inf)[1]
    fun_eff(c) = (I2(c))^2 / I1(c)
    copt = find_zero(c -> fun_eff(c) - eff, c0, Order1())
end

"""
    mscale_loss(loss::L, x)

The rho-function that is used for M-scale estimation.

For monotone (convex) functions, χ(r) = r.ψ(r).

For bounded functions, χ(r) = ρ(r)/ρ(∞) so χ(∞) = 1.
"""
mscale_loss(loss::L, x) where {L<:LossFunction} = x * RobustModels.psi(loss, x)
mscale_loss(loss::L, x) where {L<:BoundedLossFunction} =
    RobustModels.rho(loss, x) / estimator_bound(L)

"""
The M-estimate of scale is computed by solving:

```math
\\dfrac{1}{n} \\sum_i \\chi\\left( \\dfrac{r_i}{\\hat{\\sigma}}\\right) = \\delta
```

For monotone (convex) functions, χ(r) = r.ψ(r) and δ is defined as E[χ(r)] = δ for the Normal distribution N(0,1)
For bounded functions, χ(r) = ρ(r)/ρ(∞) with χ(∞) = 1 and δ = E[χ]/χ(∞) with expectation w.r.t. Normal density.

The tuning constant c corresponding to a high breakdown point (0.5)
is such that δ = 1/2, from  1/n Σ χ(r/ŝ) = δ
"""
function breakdown_point_tuning_constant(
    ::Type{L};
    bp::Real=1 / 2,
    c0::Real=1.0,
) where {L<:LossFunction}
    (0 < bp <= 1 / 2) || error("breakdown-point should be between 0 and 1/2")

    I(c) = quadgk(x -> mscale_loss(L(c), x) * 2 * exp(-x^2 / 2) / √(2π), 0, Inf)[1]
    copt = find_zero(c -> I(c) - bp, c0, Order1())
end

########
###     Loss functions
########

"The (convex) L2 loss function is that of the standard least squares problem."
struct L2Loss <: ConvexLossFunction end
L2Loss(c) = L2Loss()
rho(   ::L2Loss, r::Real) = r^2 / 2
psi(   ::L2Loss, r::Real) = r
psider(::L2Loss, r::Real) = oftype(r, 1)
weight(::L2Loss, r::Real) = oftype(r, 1)
values(::L2Loss, r::Real) = (r^2 / 2, r, oftype(r, 1))
estimator_norm(::L2Loss) = √(2 * π)
tuning_constant(::L2Loss) = 1

"""
The standard L1 loss function takes the absolute value of the residual, and is
convex but non-smooth. It is not a real L1 loss but a Huber loss
with very small tuning constant.
"""
struct L1Loss <: ConvexLossFunction end
L1Loss(c) = L1Loss()
rho(   ::L1Loss, r::Real) = abs(r)
psi(   ::L1Loss, r::Real) = sign(r)
psider(::L1Loss, r::Real) = (abs(r) < DELTA) ? oftype(r, 1) : oftype(r, 0)
weight(::L1Loss, r::Real) = (abs(r) < DELTA) ? L1WDELTA     : 1 / abs(r)
function values(est::L1Loss, r::Real)
    rr = abs(r)
    return (rr, sign(r), ((rr < DELTA) ? L1WDELTA : 1 / rr))
end
estimator_norm(::L1Loss) = 2
tuning_constant(::L1Loss) = 1

"""
The convex Huber loss function switches from between quadratic and linear cost/loss
function at a certain cutoff.
"""
struct HuberLoss <: ConvexLossFunction
    c::Float64

    HuberLoss(c::Real) = new(c)
    HuberLoss() = new(estimator_high_efficiency_constant(HuberLoss))
end

_rho(   l::HuberLoss, r::Real) = (abs(r) <= 1) ? r^2 / 2      : (abs(r) - 1 / 2)
_psi(   l::HuberLoss, r::Real) = (abs(r) <= 1) ? r            : sign(r)
_psider(l::HuberLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : oftype(r, 0)
_weight(l::HuberLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : 1 / abs(r)
function values(l::HuberLoss, r::Real)
    rr = abs(r)
    if rr <= l.c
        return ((rr / l.c)^2 / 2, r, oftype(r, 1))
    else
        return (rr / l.c - 1 / 2, l.c * sign(r), l.c / rr)
    end
end
estimator_norm(l::HuberLoss) = l.c * 2.92431
estimator_high_efficiency_constant(::Type{HuberLoss}) = 1.345
estimator_high_breakdown_point_constant(::Type{HuberLoss}) = 0.6745

"""
The convex L1-L2 loss interpolates smoothly between L2 behaviour for small
residuals and L1 for outliers.
"""
struct L1L2Loss <: ConvexLossFunction
    c::Float64

    L1L2Loss(c::Real) = new(c)
    L1L2Loss() = new(estimator_high_efficiency_constant(L1L2Loss))
end
_rho(   l::L1L2Loss, r::Real) = (sqrt(1 + r^2) - 1)
_psi(   l::L1L2Loss, r::Real) = r / sqrt(1 + r^2)
_psider(l::L1L2Loss, r::Real) = 1 / (1 + r^2)^(3 / 2)
_weight(l::L1L2Loss, r::Real) = 1 / sqrt(1 + r^2)
function values(l::L1L2Loss, r::Real)
    sqr = sqrt(1 + (r / l.c)^2)
    return ((sqr - 1), r / sqr, 1 / sqr)
end
estimator_norm(l::L1L2Loss) = l.c * 3.2723
estimator_high_efficiency_constant(::Type{L1L2Loss}) = 1.287
estimator_high_breakdown_point_constant(::Type{L1L2Loss}) = 0.8252

"""
The (convex) "fair" loss switches from between quadratic and linear
cost/loss function at a certain cutoff, and is C3 but non-analytic.
"""
struct FairLoss <: ConvexLossFunction
    c::Float64

    FairLoss(c::Real) = new(c)
    FairLoss() = new(estimator_high_efficiency_constant(FairLoss))
end
_rho(   l::FairLoss, r::Real) = abs(r) - log(1 + abs(r))
_psi(   l::FairLoss, r::Real) = r / (1 + abs(r))
_psider(l::FairLoss, r::Real) = 1 / (1 + abs(r))^2
_weight(l::FairLoss, r::Real) = 1 / (1 + abs(r))
function values(l::FairLoss, r::Real)
    ir = 1 / (1 + abs(r / l.c))
    return (abs(r) / l.c + log(ir), r * ir, ir)
end
estimator_norm(l::FairLoss) = l.c * 4
estimator_high_efficiency_constant(::Type{FairLoss}) = 1.400
estimator_high_breakdown_point_constant(::Type{FairLoss}) = 1.4503

"""
The convex Log-Cosh loss function
log(cosh(r))
"""
struct LogcoshLoss <: ConvexLossFunction
    c::Float64

    LogcoshLoss(c::Real) = new(c)
    LogcoshLoss() = new(estimator_high_efficiency_constant(LogcoshLoss))
end
_rho(   l::LogcoshLoss, r::Real) = log(cosh(r))
_psi(   l::LogcoshLoss, r::Real) = tanh(r)
_psider(l::LogcoshLoss, r::Real) = 1 / (cosh(r))^2
_weight(l::LogcoshLoss, r::Real) = (abs(r) < DELTA) ? (1 - (r)^2 / 3) : tanh(r) / r
function values(l::LogcoshLoss, r::Real)
    tr = l.c * tanh(r / l.c)
    rr = abs(r / l.c)
    return (log(cosh(rr)), tr, ((rr < DELTA) ? (1 - rr^2 / 3) : tr / r))
end
estimator_norm(l::LogcoshLoss) = l.c * π
estimator_high_efficiency_constant(::Type{LogcoshLoss}) = 1.2047
estimator_high_breakdown_point_constant(::Type{LogcoshLoss}) = 0.7479

"""
The convex Arctan loss function
r * arctan(r) - 1/2*log(1 + r^2)
"""
struct ArctanLoss <: ConvexLossFunction
    c::Float64

    ArctanLoss(c::Real) = new(c)
    ArctanLoss() = new(estimator_high_efficiency_constant(ArctanLoss))
end
_rho(   l::ArctanLoss, r::Real) = r * atan(r) - 1 / 2 * log(1 + r^2)
_psi(   l::ArctanLoss, r::Real) = atan(r)
_psider(l::ArctanLoss, r::Real) = 1 / (1 + r^2)
_weight(l::ArctanLoss, r::Real) = (abs(r) < DELTA) ? (1 - r^2 / 3) : atan(r) / r
function values(l::ArctanLoss, r::Real)
    ar = atan(r / l.c)
    rr = abs(r / l.c)
    return (
        r * ar / l.c - 1 / 2 * log(1 + rr^2),
        l.c * ar,
        ((rr < DELTA) ? (1 - rr^2 / 3) : l.c * ar / r),
    )
end
estimator_norm(l::ArctanLoss) = l.c * 2.98151
estimator_high_efficiency_constant(::Type{ArctanLoss}) = 0.919
estimator_high_breakdown_point_constant(::Type{ArctanLoss}) = 0.612

"""
The non-convex Cauchy loss function switches from between quadratic behaviour to
logarithmic tails. This rejects outliers but may result in multiple minima.
For scale estimate, r.ψ(r) is used as a loss, which is the same as for Geman loss.
"""
struct CauchyLoss <: LossFunction
    c::Float64

    CauchyLoss(c::Real) = new(c)
    CauchyLoss() = new(estimator_high_efficiency_constant(CauchyLoss))
end
_rho(   l::CauchyLoss, r::Real) = log(1 + r^2) # * 1/2  # remove factor 1/2 so the loss has a norm
_psi(   l::CauchyLoss, r::Real) = r / (1 + r^2)
_psider(l::CauchyLoss, r::Real) = (1 - r^2) / (1 + r^2)^2
_weight(l::CauchyLoss, r::Real) = 1 / (1 + r^2)
function values(l::CauchyLoss, r::Real)
    ir = 1 / (1 + (r / l.c)^2)
    return (-log(ir), r * ir, ir)
end
estimator_norm(l::CauchyLoss) = l.c * π
isconvex(::CauchyLoss) = false
isbounded(::CauchyLoss) = false

estimator_high_efficiency_constant(::Type{CauchyLoss}) = 2.385
estimator_high_breakdown_point_constant(::Type{CauchyLoss}) = 1.468

"""
The non-convex Geman-McClure for strong supression of outliers and does not guarantee a unique solution.
For the S-Estimator, it is equivalent to the Cauchy loss.
"""
struct GemanLoss <: BoundedLossFunction
    c::Float64

    GemanLoss(c::Real) = new(c)
    GemanLoss() = new(estimator_high_efficiency_constant(GemanLoss))
end
_rho(   l::GemanLoss, r::Real) = 1 / 2 * r^2 / (1 + r^2)
_psi(   l::GemanLoss, r::Real) = r / (1 + r^2)^2
_psider(l::GemanLoss, r::Real) = (1 - 3 * r^2) / (1 + r^2)^3
_weight(l::GemanLoss, r::Real) = 1 / (1 + r^2)^2
function values(l::GemanLoss, r::Real)
    rr2 = (r / l.c)^2
    ir = 1 / (1 + rr2)
    return (1 / 2 * rr2 * ir, r * ir^2, ir^2)
end
estimator_norm(::GemanLoss) = Inf
isconvex(::GemanLoss) = false
isbounded(::GemanLoss) = true

estimator_bound(::Type{GemanLoss}) = 1 / 2
estimator_high_efficiency_constant(::Type{GemanLoss}) = 3.787
estimator_high_breakdown_point_constant(::Type{GemanLoss}) = 0.61200

"""
The non-convex Welsch for strong supression of ourliers and does not guarantee a unique solution
"""
struct WelschLoss <: BoundedLossFunction
    c::Float64

    WelschLoss(c::Real) = new(c)
    WelschLoss() = new(estimator_high_efficiency_constant(WelschLoss))
end
_rho(   l::WelschLoss, r::Real) = -1 / 2 * Base.expm1(-r^2)
_psi(   l::WelschLoss, r::Real) = r * exp(-r^2)
_psider(l::WelschLoss, r::Real) = (1 - 2 * r^2) * exp(-r^2)
_weight(l::WelschLoss, r::Real) = exp(-r^2)
function values(l::WelschLoss, r::Real)
    rr2 = (r / l.c)^2
    er = exp(-rr2)
    return (-1 / 2 * Base.expm1(-rr2), r * er, er)
end
estimator_norm(::WelschLoss) = Inf
isconvex(::WelschLoss) = false
isbounded(::WelschLoss) = true

estimator_bound(::Type{WelschLoss}) = 1 / 2
estimator_high_efficiency_constant(::Type{WelschLoss}) = 2.985
estimator_high_breakdown_point_constant(::Type{WelschLoss}) = 0.8165

"""
The non-convex Tukey biweight estimator which completely suppresses the outliers,
and does not guaranty a unique solution
"""
struct TukeyLoss <: BoundedLossFunction
    c::Float64

    TukeyLoss(c::Real) = new(c)
    TukeyLoss() = new(estimator_high_efficiency_constant(TukeyLoss))
end
_rho(   l::TukeyLoss, r::Real) = (abs(r) <= 1) ? 1 / 6 * (1 - (1 - r^2)^3) : 1 / 6
_psi(   l::TukeyLoss, r::Real) = (abs(r) <= 1) ? r * (1 - r^2)^2           : oftype(r, 0)
_psider(l::TukeyLoss, r::Real) = (abs(r) <= 1) ? 1 - 6 * r^2 + 5 * r^4     : oftype(r, 0)
_weight(l::TukeyLoss, r::Real) = (abs(r) <= 1) ? (1 - r^2)^2               : oftype(r, 0)
function values(l::TukeyLoss, r::Real)
    pr = (abs(r) <= l.c) * (1 - (r / l.c)^2)
    return (1 / 6 * (1 - pr^3), r * pr^2, pr^2)
end
estimator_norm(::TukeyLoss) = Inf
isconvex(::TukeyLoss) = false
isbounded(::TukeyLoss) = true

estimator_bound(::Type{TukeyLoss}) = 1 / 6
estimator_high_efficiency_constant(::Type{TukeyLoss}) = 4.685
estimator_high_breakdown_point_constant(::Type{TukeyLoss}) = 1.5476

"""
The non-convex (and bounded) optimal Yohai-Zamar loss function that
minimizes the estimator bias. It was originally introduced in
Optimal locally robust M-estimates of regression (1997) by Yohai and Zamar
with a slightly different formula.
"""
struct YohaiZamarLoss <: BoundedLossFunction
    c::Float64

    YohaiZamarLoss(c::Real) = new(c)
    YohaiZamarLoss() = new(estimator_high_efficiency_constant(YohaiZamarLoss))
end
function rho(l::YohaiZamarLoss, r::Real)
    z = (r / l.c)^2
    if (z <= 4 / 9)
        1.3846 * z
    elseif (z <= 1)
        min(1.0, 0.5514 - 2.6917 * z + 10.7668 * z^2 - 11.6640 * z^3 + 4.0375 * z^4)
    else
        oftype(r, 1)
    end
end
function psi(l::YohaiZamarLoss, r::Real)
    z = (r / l.c)^2
    if (z <= 4 / 9)
        2.7692 * r
    elseif (z <= 1)
        r * max(0, -5.3834 + 43.0672 * z - 69.984 * z^2 + 32.3 * z^3)
    else
        oftype(r, 0)
    end
end
function psider(l::YohaiZamarLoss, r::Real)
    z = (r / l.c)^2
    if (z <= 4 / 9)
        2.7692
    elseif (z <= 0.997284)  # from the root of ψ expression
        -5.3834 + 129.2016 * z - 349.92 * z^2 + 226.1 * z^3
    else
        oftype(r, 0)
    end
end
function weight(l::YohaiZamarLoss, r::Real)
    z = (r / l.c)^2
    if (z <= 4 / 9)
        2.7692
    elseif (z <= 1)
        max(0, -5.3834 + 43.0672 * z - 69.984 * z^2 + 32.3 * z^3)
    else
        oftype(r, 0)
    end
end
function values(l::YohaiZamarLoss, r::Real)
    z = (r / l.c)^2
    if (z <= 4 / 9)
        return (1.3846 * z, 2.7692 * r, 2.7692)
    elseif (z <= 1)
        ρ = min(1, 0.5514 - 2.6917 * z + 10.7668 * z^2 - 11.6640 * z^3 + 4.0375 * z^4)
        w = max(0, -5.3834 + 43.0672 * z - 69.984 * z^2 + 32.3 * z^3)
        return (ρ, r * w, w)
    else
        return (oftype(r, 1), oftype(r, 0), oftype(r, 0))
    end
end
estimator_norm(::YohaiZamarLoss) = Inf
isconvex(::YohaiZamarLoss) = false
isbounded(::YohaiZamarLoss) = true

estimator_bound(::Type{YohaiZamarLoss}) = 1
estimator_high_efficiency_constant(::Type{YohaiZamarLoss}) = 3.1806
estimator_high_breakdown_point_constant(::Type{YohaiZamarLoss}) = 1.2139

######
###   Convex sum of loss functions
######
struct CompositeLossFunction{L1<:LossFunction,L2<:LossFunction} <: LossFunction
    α1::Float64
    loss1::L1
    α2::Float64
    loss2::L2
end
function CompositeLossFunction(
    loss1::LossFunction,
    loss2::LossFunction,
    α1::Real=1.0,
    α2::Real=1.0,
)
    α1 >= 0 || throw(
        DomainError(α1, "coefficients of CompositeLossFunction are non-negative."),
    )
    α2 >= 0 || throw(
        DomainError(α2, "coefficients of CompositeLossFunction are non-negative."),
    )

    return CompositeLossFunction{typeof(loss1),typeof(loss2)}(
        float(α1),
        loss1,
        float(α2),
        loss2,
    )
end

function show(io::IO, e::CompositeLossFunction)
    mess = "CompositeLossFunction $(round(e.α1; digits=2)) ."*
           " $(e.loss1) + $(round(e.α2; digits=2)) . $(e.loss2)"
    print(io, mess)
end

Base.first(e::CompositeLossFunction) = e.loss1
Base.last(e::CompositeLossFunction) = e.loss2

isbounded(e::CompositeLossFunction) = isbounded(e.loss1) && isbounded(e.loss2)
isconvex(e::CompositeLossFunction) = isconvex(e.loss1) && isconvex(e.loss2)

estimator_norm(e::CompositeLossFunction) =
    e.α1 * estimator_norm(e.loss1) + e.α2 * estimator_norm(e.loss2)
rho(e::CompositeLossFunction, r::Real) =
    e.α1 * rho(e.loss1) * (tuning_constant(E.loss1))^2 +
    e.α2 * rho(e.loss2) * (tuning_constant(E.loss2))^2
psi(e::CompositeLossFunction, r::Real) = e.α1 * psi(e.loss1) + e.α2 * psi(e.loss2)
psider(e::CompositeLossFunction, r::Real) =
    e.α1 * psider(e.loss1) + e.α2 * psider(e.loss2)
weight(e::CompositeLossFunction, r::Real) =
    e.α1 * weight(e.loss1) + e.α2 * weight(e.loss2)
values(e::CompositeLossFunction, r::Real) =
    @.(e.α1 * values(e.loss1) + e.α2 * values(e.loss2))

######
###   Scale estimation
######

"""
    scale_estimate(loss, res; σ0=1.0, wts=[], verbose=false,
                             order=1, approx=false, nmax=30,
                             rtol=1e-4, atol=0.1)

Compute the M-scale estimate from the loss function.
If the loss is bounded, ρ is used as the function χ in the sum,
otherwise r.ψ(r) is used if the loss is not bounded, to coincide with
the Maximum Likelihood Estimator.
Also, for bounded estimator, because f(s) = 1/(nδ) Σ ρ(ri/s) is decreasing
the iteration step is not using the weights but is multiplicative.
"""
function scale_estimate(
    l::L,
    res::AbstractArray{T};
    verbose::Bool=false,
    bound::Union{Nothing,T}=nothing,
    σ0::T=1.0,
    wts::AbstractArray{T}=T[],
    order::Int=1,
    approx::Bool=false,
    nmax::Int=30,
    rtol::Real=1e-4,
    atol::Real=0.1,
) where {L<:LossFunction,T<:AbstractFloat}

    # Compute δ such that E[ρ] = δ for a Normal N(0, 1)
    if isnothing(bound)
        bound = quadgk(x -> mscale_loss(l, x) * 2 * exp(-x^2 / 2) / √(2π), 0, Inf)[1]
    end

    Nz = sum(iszero, res)
    if Nz > length(res) * (1 - bound)
        # The M-scale cannot be estimated because too many residuals are zeros.
        verbose && println(
            "there are too many zero residuals for M-scale estimation: #(r=0) > n*(1-b),"*
            " $(Nz) > $(length(res)*(1-bound))",
        )
        m = "the M-scale cannot be estimated because too many residuals are zeros: $(res)"
        throw(ConvergenceFailed(m))
    end

    # Approximate the solution with `nmax` iterations
    σn = σ0
    converged = false
    verbose && println("Initial M-scale estimate: $(σn)")
    for n in 1:nmax
        ## For non-bounded loss, we suppose that χ = r.ψ,
        ## therefore the weight ( χ/r² ) is ψ/r which is w.
        ww = weight.(l, res ./ σn)
        if length(wts) == length(res)
            ww .*= wts
        end
        σnp1 = sqrt(mean(res ./ σn, weights(ww)) / bound)

        verbose && println("M-scale update: $(σn) ->  $(σnp1)")

        ε = σnp1 / σn
        σn = σnp1
        if abs(ε - 1) < rtol
            verbose && println("M-scale converged after $(n) steps.")
            converged = true
            break
        end
    end

    if !approx
        converged || @warn "the M-scale did not converge, consider increasing the maximum" *
              " number of iterations nmax=$(nmax) or starting with a better" *
              " initial value σ0=$(σ0). Return the current estimate: $(σn)"
    end

    return σn
end

function scale_estimate(
    l::L,
    res::AbstractArray{T};
    verbose::Bool=false,
    bound::Union{Nothing,Real}=0.5,
    σ0::Real=1.0,
    wts::AbstractArray{T}=T[],
    order::Int=1,
    approx::Bool=false,
    nmax::Int=30,
    rtol::Real=1e-4,
    atol::Real=0.1,
) where {L<:BoundedLossFunction,T<:AbstractFloat}

    # Compute δ such that E[ρ] = δ for a Normal N(0, 1)
    if isnothing(bound)
        bound = quadgk(x -> mscale_loss(l, x) * 2 * exp(-x^2 / 2) / √(2π), 0, Inf)[1]
    end

    Nz = sum(iszero, res)
    if Nz > length(res) * (1 - bound)
        # The M-scale cannot be estimated because too many residuals are zeros.
        verbose && println(
            "there are too many zero residuals for M-scale estimation:"*
            " #(r=0) > n*(1-b), $(Nz) > $(length(res)*(1-bound))",
        )
        m = "the M-scale cannot be estimated because too many residuals are zeros: $(res)"
        throw(ConvergenceFailed(m))
    end

    # Approximate the solution with `nmax` iterations
    if !approx && order >= 2
        bb = tuning_constant(l)^2 * estimator_bound(L) * bound
    end
    σn = σ0
    converged = false
    verbose && println("Initial M-scale estimate: $(σn)")
    ε = 0
    for n in 1:nmax
        rr = res / σn
        ε = if length(wts) == length(res)
            mean(mscale_loss.(l, rr), weights(wts)) / bound
        else
            mean(x -> mscale_loss(l, x), rr) / bound
        end
        verbose && println("M-scale 1st order update: $(ε)")

        ## Implemented, but it gives worst results than 1st order...
        ## Uses Newton's method to find the root of
        ## f(σ) = log( 1/(nb) Σ ρ(r/σ) ) = 0
        if !approx && order >= 2
            εp = if length(wts) == length(res)
                mean(rr .* psi.(l, rr), weights(wts)) / bb
            else
                mean(x -> x * psi(l, x), rr) / bb
            end
            # Use the gradient only if the scale is not too small
            if isnan(εp) || !isfinite(εp) || εp <= atol
                verbose && println("M-scale order is set to 1 for this iteration")
            else
                ε = exp(ε * log(ε) / εp)
            end
        end

        σnp1 = σn * ε
        verbose && println("M-scale update: $(ε) : $(σn) ->  $(σnp1)")

        σn = σnp1
        if abs(ε - 1) < rtol
            verbose && println("M-scale converged after $(n) steps.")
            converged = true
            break
        end
    end

    if !approx && !converged
        @warn("the M-scale did not converge ε=$(round(abs(ε-1); digits=5)),"*
              " consider increasing the maximum number of iterations nmax=$(nmax)"*
              " or starting with a better initial value σ0=$(σ0)."*
              " Return the current estimate: $(σn)"
        )
    end
    return σn
end

######
###   M-Estimators
######
"""
    MEstimator{L<:LossFunction} <: AbstractEstimator

M-estimator for a given loss function.

The M-estimator is obtained by minimizing the loss function:

```math
\\hat{\\mathbf{\\beta}} = \\underset{\\mathbf{\\beta}}{\\textrm{argmin}} \\sum_{i=1}^n \\rho\\left(\\dfrac{r_i}{\\hat{\\sigma}}\\right)
```

with the residuals  ``\\mathbf{r} = \\mathbf{y} - \\mathbf{X} \\mathbf{\\beta}`` ,
and a robust scale estimate ``\\hat{\\sigma}``.

# Fields
- `loss`: the [`LossFunction`](@ref) used for the robust estimation.

"""
struct MEstimator{L<:LossFunction} <: AbstractEstimator
    loss::L
end
MEstimator{L}() where {L<:LossFunction} = MEstimator(efficient_loss(L))
MEstimator(::Type{L}) where {L<:LossFunction} = MEstimator(efficient_loss(L))

loss(e::MEstimator) = e.loss

function show(io::IO, obj::MEstimator)
    print(io, "M-Estimator($(obj.loss))")
end

# Forward all methods to the `loss` field
rho(e::MEstimator, r::Real) = rho(e.loss, r)
psi(e::MEstimator, r::Real) = psi(e.loss, r)
psider(e::MEstimator, r::Real) = psider(e.loss, r)
weight(e::MEstimator, r::Real) = weight(e.loss, r)
values(e::MEstimator, r::Real) = values(e.loss, r)
estimator_norm(e::MEstimator, args...) = estimator_norm(e.loss, args...)
estimator_bound(e::MEstimator) = estimator_bound(typeof(e.loss))

isbounded(e::MEstimator) = isbounded(e.loss)
isconvex(e::MEstimator) = isconvex(e.loss)

scale_estimate(est::E, res; kwargs...) where {E<:MEstimator} =
    scale_estimate(est.loss, res; kwargs...)

"`L1Estimator` is a shorthand name for `MEstimator{L1Loss}`. Using exact QuantileRegression should be prefered."
const L1Estimator = MEstimator{L1Loss}

"`L2Estimator` is a shorthand name for `MEstimator{L2Loss}`, the non-robust OLS."
const L2Estimator = MEstimator{L2Loss}

######
###   S-Estimators
######

"""
    SEstimator{L<:BoundedLossFunction} <: AbstractEstimator

S-estimator for a given bounded loss function.

The S-estimator is obtained by minimizing the scale estimate:

```math
\\hat{\\mathbf{\\beta}} = \\underset{\\mathbf{\\beta}}{\\textrm{argmin }} \\hat{\\sigma}^2
```

where the robust scale estimate ``\\hat{\\sigma}}`` is solution of:

```math
\\dfrac{1}{n} \\sum_{i=1}^n \\rho\\left(\\dfrac{r_i}{\\hat{\\sigma}}\\right) = \\delta
```

with the residuals  ``\\mathbf{r} = \\mathbf{y} - \\mathbf{X} \\mathbf{\\beta}`` ,
``\\rho`` is a bounded loss function with  ``\\underset{r \\to \\infty}{\\lim} \\rho(r) = 1`` and
``\\delta`` is the finite breakdown point, usually 0.5.

# Fields
- `loss`: the [`LossFunction`](@ref) used for the robust estimation.

"""
struct SEstimator{L<:BoundedLossFunction} <: AbstractEstimator
    loss::L
end
SEstimator{L}() where {L<:BoundedLossFunction} = SEstimator(robust_loss(L))
SEstimator(::Type{L}) where {L<:BoundedLossFunction} = SEstimator(robust_loss(L))

loss(e::SEstimator) = e.loss

function show(io::IO, obj::SEstimator)
    print(io, "S-Estimator($(obj.loss))")
end

# Forward all methods to the `loss` field
rho(e::SEstimator, r::Real) = rho(e.loss, r)
psi(e::SEstimator, r::Real) = psi(e.loss, r)
psider(e::SEstimator, r::Real) = psider(e.loss, r)
weight(e::SEstimator, r::Real) = weight(e.loss, r)
values(e::SEstimator, r::Real) = values(e.loss, r)
estimator_norm(e::SEstimator, args...) = Inf
estimator_bound(e::SEstimator) = estimator_bound(typeof(e.loss))
isbounded(e::SEstimator) = true
isconvex(e::SEstimator) = false

scale_estimate(est::E, res; kwargs...) where {E<:SEstimator} =
    scale_estimate(est.loss, res; kwargs...)

######
###   MM-Estimators
######

"""
    MMEstimator{L1<:BoundedLossFunction, L2<:LossFunction} <: AbstractEstimator

MM-estimator for the given loss functions.

The MM-estimator is obtained using a two-step process:

1. compute a robust scale estimate with a high breakdown point using a S-estimate and the loss function `L1`.
2. compute an efficient estimate using a M-estimate with the loss function `L2`.

# Fields
- `loss1`: the [`BoundedLossFunction`](@ref) used for the high breakdown point S-estimation.
- `loss2`: the [`LossFunction`](@ref) used for the efficient M-estimation.
- `scaleest`: boolean specifying the if the estimation is in the S-estimation step (`true`)
or the M-estimation step (`false`).

"""
mutable struct MMEstimator{L1<:BoundedLossFunction,L2<:LossFunction} <: AbstractEstimator
    "high breakdown point loss function"
    loss1::L1

    "high efficiency loss function"
    loss2::L2

    "S-Estimator phase indicator (or M-Estimator phase)"
    scaleest::Bool

    MMEstimator{L1,L2}(
        loss1::L1,
        loss2::L2,
        scaleest::Bool=true,
    ) where {L1<:BoundedLossFunction,L2<:LossFunction} = new(loss1, loss2, scaleest)
end
MMEstimator(
    loss1::L1,
    loss2::L2,
    scaleest::Bool,
) where {L1<:BoundedLossFunction,L2<:LossFunction} =
    MMEstimator{L1,L2}(loss1, loss2, scaleest)
MMEstimator(loss1::L1, loss2::L2) where {L1<:BoundedLossFunction,L2<:LossFunction} =
    MMEstimator{L1,L2}(loss1, loss2, true)
MMEstimator(::Type{L1}, ::Type{L2}) where {L1<:BoundedLossFunction,L2<:LossFunction} =
    MMEstimator(robust_loss(L1), efficient_loss(L2))
MMEstimator{L}() where {L<:BoundedLossFunction} =
    MMEstimator(robust_loss(L), efficient_loss(L))
MMEstimator(::Type{L}) where {L<:BoundedLossFunction} = MMEstimator{L}()

loss(e::MMEstimator) = e.scaleest ? e.loss1 : e.loss2

"MEstimator, set to S-Estimation phase"
set_SEstimator(e::MMEstimator) = (e.scaleest = true; e)

"MEstimator, set to M-Estimation phase"
set_MEstimator(e::MMEstimator) = (e.scaleest = false; e)

function show(io::IO, obj::MMEstimator)
    print(io, "MM-Estimator($(obj.loss1), $(obj.loss2))")
end

# Forward all methods to the selected loss
rho(E::MMEstimator, r::Real) = rho(loss(E), r)
psi(E::MMEstimator, r::Real) = psi(loss(E), r)
psider(E::MMEstimator, r::Real) = psider(loss(E), r)
weight(E::MMEstimator, r::Real) = weight(loss(E), r)
values(E::MMEstimator, r::Real) = values(loss(E), r)

# For these methods, only the SEstimator loss is useful,
# not the MEstimator, so E.loss1 is used instead of loss(E)
estimator_bound(E::MMEstimator) = estimator_bound(typeof(E.loss1))
# For these methods, only the MEstimator loss is useful,
# not the SEstimator, so E.loss2 is used instead of loss(E)
estimator_norm(E::MMEstimator, args...) = estimator_norm(E.loss2, args...)
isbounded(E::MMEstimator) = isbounded(E.loss2)
isconvex(E::MMEstimator) = isconvex(E.loss2)

scale_estimate(est::E, res; kwargs...) where {E<:MMEstimator} =
    scale_estimate(est.loss1, res; kwargs...)

######
###   τ-Estimators
######

"""
    TauEstimator{L1<:BoundedLossFunction, L2<:BoundedLossFunction} <: AbstractEstimator

τ-estimator for the given loss functions.

The τ-estimator corresponds to a M-estimation, where the loss function is a weighted
sum of a high breakdown point loss and an efficient loss. The weight is recomputed at
every step of the Iteratively Reweighted Least Square, so the estimate is both robust
(high breakdown point) and efficient.

# Fields
- `loss1`: the high breakdown point [`BoundedLossFunction`](@ref).
- `loss2`: the high efficiency [`BoundedLossFunction`](@ref).
- `w`: the weight in the sum of losses: `w . loss1 + loss2`.

"""
mutable struct TauEstimator{L1<:BoundedLossFunction,L2<:BoundedLossFunction} <:
               AbstractEstimator
    "high breakdown point loss function"
    loss1::L1

    "high efficiency loss function"
    loss2::L2

    "loss weight"
    w::Float64

    TauEstimator{L1,L2}(
        l1::L1,
        l2::L2,
        w::Real=0.0,
    ) where {L1<:BoundedLossFunction,L2<:BoundedLossFunction} = new(l1, l2, float(w))
end
TauEstimator(
    l1::L1,
    l2::L2,
    args...,
) where {L1<:BoundedLossFunction,L2<:BoundedLossFunction} =
    TauEstimator{L1,L2}(l2, l2, args...)
# Warning: The tuning constant of the the efficient loss is NOT optimized for different loss functions
TauEstimator(
    ::Type{L1},
    ::Type{L2},
    args...,
) where {L1<:BoundedLossFunction,L2<:BoundedLossFunction} =
    TauEstimator(robust_loss(L1), efficient_loss(L2), args...)

# With the same loss function, the tuning constant of the the efficient loss is optimized
TauEstimator{L}() where {L<:BoundedLossFunction} =
    TauEstimator(robust_loss(L), L(estimator_tau_efficient_constant(L)))
TauEstimator(::Type{L}) where {L<:BoundedLossFunction} = TauEstimator{L}()

loss(e::TauEstimator) = CompositeLossFunction(e.loss1, e.loss2, e.w, 1)

function show(io::IO, obj::TauEstimator)
    print(io, "τ-Estimator($(obj.loss1), $(obj.loss2))")
end

"""
    tau_efficiency_tuning_constant(::Type{L1}, ::Type{L2}; eff::Real=0.95, c0::Real=1.0)
        where {L1<:BoundedLossFunction,L2<:BoundedLossFunction}

Compute the tuning constant that corresponds to a high breakdown point for the τ-estimator.
"""
function tau_efficiency_tuning_constant(
    ::Type{L1},
    ::Type{L2};
    eff::Real=0.95,
    c0::Real=1.0,
) where {L1<:BoundedLossFunction,L2<:BoundedLossFunction}
    loss1 = L1(estimator_high_breakdown_point_constant(L1))
    w1 = quadgk(x -> x * psi(loss1, x) * 2 * exp(-x^2 / 2) / √(2π), 0, Inf)[1]

    function τest(c)
        loss2 = L2(c)
        t2 = (tuning_constant(loss2))^2
        w2 = quadgk(
            x -> (2 * rho(loss2, x) * t2 - x * psi(loss2, x)) * 2 * exp(-x^2 / 2) / √(2π),
            0,
            Inf,
        )[1]
        TauEstimator{L1,L2}(loss1, loss2, w2 / w1)
    end

    lpsi(x, c) = psi(τest(c), x)
    lpsip(x, c) = psider(τest(c), x)

    I1(c) = quadgk(x -> (lpsi(x, c))^2 * 2 * exp(-x^2 / 2) / √(2π), 0, Inf)[1]
    I2(c) = quadgk(x -> lpsip(x, c) * 2 * exp(-x^2 / 2) / √(2π), 0, Inf)[1]
    fun_eff(c) = (I2(c))^2 / I1(c)
    copt = find_zero(c -> fun_eff(c) - eff, c0, Order1())
end
tau_efficiency_tuning_constant(::Type{L}; kwargs...) where {L<:BoundedLossFunction} =
    tau_efficiency_tuning_constant(L, L; kwargs...)

"The tuning constant associated to the loss that gives a robust τ-estimator."
estimator_tau_efficient_constant(::Type{GemanLoss}) = 5.632
estimator_tau_efficient_constant(::Type{WelschLoss}) = 4.043
estimator_tau_efficient_constant(::Type{TukeyLoss}) = 6.040
estimator_tau_efficient_constant(::Type{YohaiZamarLoss}) = 3.270

"""
    update_weight!(E::TauEstimator, res::AbstractArray{T}; wts::AbstractArray{T}=T[])

Update the weight between the two estimators of a τ-estimator using the scaled residual.
"""
function update_weight!(
    E::TauEstimator,
    res::AbstractArray{T};
    wts::AbstractArray{T}=T[],
) where {T<:AbstractFloat}
    c² = (tuning_constant(E.loss2))^2
    E.w = if length(wts) == length(res)
        w2 = sum(@.(wts * (2 * rho(E.loss2, res) * c² - res * psi(E.loss2, res))))
        w1 = sum(@.(wts * res * psi(E.loss1, res)))
        w2 / w1
    else
        w2 = sum(r -> 2 * rho(E.loss2, r) * c² - r * psi(E.loss2, r), res)
        w1 = sum(r -> r * psi(E.loss1, r), res)
        w2 / w1
    end
    E
end
update_weight!(E::TauEstimator, res::AbstractArray; wts::AbstractArray=[]) =
    update_weight!(E, float(res); wts=float(wts))
update_weight!(E::TauEstimator, w::Real) = (E.w = w; E)

# Forward all methods to the `loss` fields
rho(E::TauEstimator, r::Real) =
    E.w * rho(E.loss1, r) * (tuning_constant(E.loss1))^2 +
    rho(E.loss2, r) * (tuning_constant(E.loss2))^2
psi(E::TauEstimator, r::Real) = E.w * psi(E.loss1, r) + psi(E.loss2, r)
psider(E::TauEstimator, r::Real) = E.w * psider(E.loss1, r) + psider(E.loss2, r)
weight(E::TauEstimator, r::Real) = E.w * weight(E.loss1, r) + weight(E.loss2, r)
function values(E::TauEstimator, r::Real)
    vals1 = values(E.loss1, r)
    vals2 = values(E.loss2, r)
    c12, c22 = (tuning_constant(E.loss1))^2, (tuning_constant(E.loss2))^2
    return (
        E.w * vals1[1] * c12 + vals2[1] * c22,
        E.w * vals1[2] + vals2[3],
        E.w * vals1[3] + vals2[3],
    )
end
estimator_norm(E::TauEstimator, args...) = Inf
estimator_bound(E::TauEstimator) = estimator_bound(typeof(E.loss1))
isbounded(E::TauEstimator) = true
isconvex(E::TauEstimator) = false

scale_estimate(est::E, res; kwargs...) where {E<:TauEstimator} =
    scale_estimate(est.loss1, res; kwargs...)

"""
    tau_scale_estimate!(E::TauEstimator, res::AbstractArray{T}, σ::Real, sqr::Bool=false;
                        wts::AbstractArray{T}=T[], bound::AbstractFloat=0.5) where {T<:AbstractFloat}

The τ-scale estimate, where `σ` is the scale estimate from the robust M-scale.
If `sqr` is true, return the squared value.
"""
function tau_scale_estimate(
    est::TauEstimator,
    res::AbstractArray{T},
    σ::Real,
    sqr::Bool=false;
    wts::AbstractArray=[],
    bound::AbstractFloat=0.5,
) where {T<:AbstractFloat}
    t = if length(wts) == length(res)
        mean(mscale_loss.(Ref(est.loss2), res ./ σ), weights(r.wts)) / bound
    else
        mean(mscale_loss.(Ref(est.loss2), res ./ σ)) / bound
    end

    return sqr ? σ * √t : σ^2 * t
end

######
###   MQuantile Estimators
######
"""
    quantile_weight(τ::Real, r::Real)

Wrapper function to compute quantile-like loss function.
"""
quantile_weight(τ::Real, r::Real) = oftype(r, 2 * ifelse(r > 0, τ, 1 - τ))

"""
    GeneralizedQuantileEstimator{L<:LossFunction} <: AbstractQuantileEstimator

Generalized Quantile Estimator is an M-Estimator with asymmetric loss function.

For [`L1Loss`](@ref), this corresponds to quantile regression (although it is better
to use [`quantreg`](@ref) for quantile regression because it gives the exact solution).

For [`L2Loss`](@ref), this corresponds to Expectile regression (see [`ExpectileEstimator`](@ref)).

# Fields
- `loss`: the [`LossFunction`](@ref).
- `τ`: the quantile value to estimate, between 0 and 1.

# Properties
- `tau`, `q`, `quantile` are aliases for `τ`.
"""
mutable struct GeneralizedQuantileEstimator{L<:LossFunction} <: AbstractQuantileEstimator
    loss::L
    τ::Float64
end
function GeneralizedQuantileEstimator(l::L, τ::Real=0.5) where {L<:LossFunction}
    (0 < τ < 1) ||
        throw(DomainError(τ, "quantile should be a number between 0 and 1 excluded"))
    GeneralizedQuantileEstimator{L}(l, float(τ))
end
GeneralizedQuantileEstimator{L}(τ::Real=0.5) where {L<:LossFunction} =
    GeneralizedQuantileEstimator(L(), float(τ))

function ==(
    e1::GeneralizedQuantileEstimator{L1},
    e2::GeneralizedQuantileEstimator{L2},
) where {L1<:LossFunction,L2<:LossFunction}
    if (L1 !== L2) || (loss(e1) != loss(e2)) || (e1.τ != e2.τ)
        return false
    end
    return true
end
function show(io::IO, obj::GeneralizedQuantileEstimator)
    print(io, "MQuantile($(obj.τ), $(obj.loss))")
end
loss(e::GeneralizedQuantileEstimator) = e.loss

function Base.getproperty(r::GeneralizedQuantileEstimator, s::Symbol)
    if s ∈ (:tau, :q, :quantile)
        r.τ
    else
        getfield(r, s)
    end
end

function Base.setproperty!(r::GeneralizedQuantileEstimator, s::Symbol, v)
    if s ∈ (:tau, :q, :quantile)
        (0 < v < 1) ||
            throw(DomainError(v, "quantile should be a number between 0 and 1 excluded"))
        r.τ = float(v)
    else
        setfield!(r, s, v)
    end
end

Base.propertynames(r::GeneralizedQuantileEstimator, private=false) =
    (:loss, :τ, :tau, :q, :quantile)

# Forward all methods to the `loss` field
rho(e::GeneralizedQuantileEstimator, r::Real) = quantile_weight(e.τ, r) * rho(e.loss, r)
psi(e::GeneralizedQuantileEstimator, r::Real) = quantile_weight(e.τ, r) * psi(e.loss, r)
psider(e::GeneralizedQuantileEstimator, r::Real) =
    quantile_weight(e.τ, r) * psider(e.loss, r)
weight(e::GeneralizedQuantileEstimator, r::Real) =
    quantile_weight(e.τ, r) * weight(e.loss, r)
function values(e::GeneralizedQuantileEstimator, r::Real)
    w = quantile_weight(e.τ, r)
    vals = values(e.loss, r)
    Tuple([x * w for x in vals])
end
estimator_norm(e::GeneralizedQuantileEstimator, args...) = estimator_norm(e.loss, args...)
estimator_bound(e::GeneralizedQuantileEstimator) = estimator_bound(loss(e))
isbounded(e::GeneralizedQuantileEstimator) = isbounded(e.loss)
isconvex(e::GeneralizedQuantileEstimator) = isconvex(e.loss)

function scale_estimate(est::E, res; kwargs...) where {E<:GeneralizedQuantileEstimator}
    error("the M-scale estimate of a generalized quantile estimator is not defined")
end

"""
The expectile estimator is a generalization of the L2 estimator, for other quantile τ ∈ [0,1].

[1] Schnabel, Eilers - Computational Statistics and Data Analysis 53 (2009) 4168–4177 - Optimal expectile smoothing
doi:10.1016/j.csda.2009.05.002
"""
const ExpectileEstimator = GeneralizedQuantileEstimator{L2Loss}

"Non-exact quantile estimator, `GeneralizedQuantileEstimator{L1Loss}`. Prefer using [`QuantileRegression`](@ref)"
const QuantileEstimator = GeneralizedQuantileEstimator{L1Loss}

const UnionL1 = Union{L1Estimator,GeneralizedQuantileEstimator{L1Loss}}

const UnionMEstimator = Union{MEstimator,GeneralizedQuantileEstimator}


text:
## Author:  Pearson
## Date: April 2020
## Julia version: 1.3.1
## Purpose: Functions for EKNR(2016) Full Model - Baseline Equilibrium

#=
    Subroutine 1
    - This function solves for the price of final products using equation (A.30) and (A.31).
    - Note that this function takes changes of factor price of labor and capital as given.
=#
function factor_price_fixpoint!(
    res_fixpoint::AbstractArray,
    guess_fixpoint::AbstractArray, # Product price, dim = NC*NS
    exos_fixpoint::NamedTuple,
    params_fixpoint::NamedTuple,
    )

    # Unpack exogenous variables and parameters
    @unpack ŵ, r̂, π, d̂, T̂ = exos_fixpoint
    @unpack NC, NS, NK, β̃ᴸ, β̃ᴷ, β̃ᴹ, θ = params_fixpoint

    # Pre-allocate memory
    p̂ = similar(guess_fixpoint)
    b̂ = similar(guess_fixpoint)
    p̂′ = zeros(eltype(guess_fixpoint),size(guess_fixpoint))

    # Resolve guess
    p̂[1:NC,1:NS] = guess_fixpoint

    # Generate other variables
    for n in 1:NC
        for l in 1:NS
            b̂[n,l] = ŵ[n]^β̃ᴸ[n,l]
            for k in 1:NK
                b̂[n,l] *= r̂[n,k]^β̃ᴷ[n,l,k]
            end
            for j in 1:NS
                #=if p̂[n,j]<0
                    println(guess_fixpoint)
                end=#
                b̂[n,l] *= p̂[n,j]^β̃ᴹ[n,l,j]
            end
        end
    end

    # Equations
    for n in 1:NC
        for j in 1:NS
            for i = 1:NC
                p̂′[n,j] += π[n,i,j]*(b̂[i,j]*d̂[n,i,j]/T̂[i,j])^-θ
            end
            p̂′[n,j] = p̂′[n,j]^(-1/θ)
            res_fixpoint[n,j] = p̂′[n,j]/p̂[n,j]-1
        end
    end

    return res_fixpoint, p̂, b̂
end

#=
    Subroutine 2
    - This function solves the changes of GDP in sector 3, Semidurable(S) in appendix or
    Nondurable(N) in the paper.
    - Note that this function takes changes of GDP in sector 1 and 2, Construction(C) and
    Durables(D) as given.
=#
function static_problem!(
    res_static::AbstractArray,
    guess_static::AbstractArray, # Product price, dim = NC*NS
    exos_static::NamedTuple,
    params_static::NamedTuple,
    )
    # Unpack exogenous variables and parameters
    @unpack π, Ŷᴷ, Y, Xᶠ, Dᴿ, wL, L̂, rK, K̂, d̂, T̂ = exos_static
    @unpack NC, NS, NK, β̃ᴸ, β̃ᴷ, ψ, θ, β̃ᴹ = params_static

    # Pre-allocate memory, note that Ωᵣ⋆={C,D,S}, Ωₖ={C,D}
    Ŷ = zeros(eltype(guess_static),NC,NS) # changes of sectoral GDP, (𝒩,Ωᵣ⋆)
    ŵ = zeros(eltype(guess_static),NC) # changes of labor wage, (𝒩)
    r̂ = zeros(eltype(guess_static),NC,NK) # changes of capital rental rate, (𝒩,Ωₖ)
    guess_fixpoint = zeros(eltype(guess_static),NC,NS) # goods price guess, (𝒩,Ωᵣ⋆)
    π̂ = zeros(eltype(guess_static),size(π)) # changes of trade share, (𝒩,𝒩,Ωᵣ⋆)
    π′ = similar(π̂) # level of trade share in the following period, (𝒩,𝒩,Ωᵣ⋆)
    Y′ = zeros(eltype(guess_static),NC,NS) # level of sectoral GDP, (𝒩,Ωᵣ⋆)
    Xˢ = zeros(eltype(guess_static),NC) # level of final demand for Semidurable(S), (𝒩)
    RHS = zeros(eltype(guess_static),NC) # Right hand side of Step 7, (𝒩)

    # Step 1
    # Resolve guess
    Ŷ[1:NC,1:2] = Ŷᴷ
    Ŷ[1:NC,3] = guess_static

    # Step 2 and 3
    # Generate other variables
    for n = 1:NC
        for j = 1:NS
            ŵ[n] += β̃ᴸ[n,j]*Y[n,j]*Ŷ[n,j]
        end
        ŵ[n] += β̃ᴸ[n,4]*(Xᶠ[n,4]-Dᴿ[n])
        ŵ[n] /= wL[n]*L̂[n]

        for k = 1:NK
            for j = 1:NS
                r̂[n,k] += β̃ᴷ[n,j,k]*Y[n,j]*Ŷ[n,j]
            end
            r̂[n,k] += β̃ᴷ[n,4,k]*(Xᶠ[n,4]-Dᴿ[n])+(Xᶠ[n,3]+Xᶠ[n,4])*ψ[k]/ψ[3]
            r̂[n,k] /= rK[n,k]*K̂[n,k]
        end
    end

    # Step 4: Solve for the goods price
    # Form the guess
    for n = 1:NC
        for l = 1:NS
            guess_fixpoint[n,l] = ŵ[n]^β̃ᴸ[n,l]
            for k = 1:NK
                guess_fixpoint[n,l] *= r̂[n,k]^β̃ᴷ[n,l,k]
            end
        end
    end

    # Pack exogenous variables and parameters for solve the price
    myexos_fixpoint = @with_kw (
        ŵ = ŵ[1:NC], r̂ = r̂[1:NC,1:NK], π = π[1:NC,1:NC,1:NS], d̂ = d̂[1:NC,1:NC,1:NS], T̂ = T̂[1:NC,1:NS]
    )
    myparams_fixpoint = @with_kw (
        NC = NC, NS = NS, NK = NK, β̃ᴸ = β̃ᴸ, β̃ᴷ = β̃ᴷ, β̃ᴹ = β̃ᴹ, θ = θ
    )
    exos_fixpoint = myexos_fixpoint()
    params_fixpoint = myparams_fixpoint()

    # Solve the fix point problem
    # println("Start to solve the fix point problem.")
    #println("Run time and memory cost:")
    #@time results_fixpoint =
    #    try
            results_fixpoint = nlsolve(
                (res_fixpoint, guess_fixpoint) -> factor_price_fixpoint!(
                    res_fixpoint, guess_fixpoint, exos_fixpoint, params_fixpoint),
                guess_fixpoint,
                ftol=1e-12,
                method=:anderson,
                show_trace=false,
            )
    #    catch err
    #        if isa(err, DomainError)
    #            error("Failed to solve the fix point problem, please try again.")
    #        end
    #    end

    # Check Convergence
    #converged(results_fixpoint) || error("Failed to converge in $(results_fixpoint.iterations) iterations.")
    #println("Successfully solved the fix point problem.\n")

    # Catch Solutions
    res_fixpoint = similar(results_fixpoint.zero)
    res_fixpoint, p̂, b̂ = factor_price_fixpoint!(
        res_fixpoint, results_fixpoint.zero, exos_fixpoint, params_fixpoint)

    # Step 5 and 6
    # Get changes of trade shares for Semidurable(S) in the appendix or Nondurable(N) in the paper
    # Form the trade share matrix at t+1
    for n = 1:NC
        for i = 1:NC
            for j = 2:NS
                π̂[n,i,j] = (b̂[i,j]*d̂[n,i,j]/T̂[i,j]/p̂[n,j])^-θ
                π′[n,i,j] = π[n,i,j]*π̂[n,i,j]
            end
            π̂[n,i,1] = 1.0
            π′[n,i,1] = 0.0
        end
        π′[n,n,1] = 1.0
    end

    # Step 7: back out Ŷ, here goes the equation
    for n = 1:NC
        for j = 1:NS
            Y′[n,j] = Ŷ[n,j]*Y[n,j]
        end
    end
    Xˢ[1:NC] = π′[1:NC,1:NC,3]'\(Y′[1:NC,3])

    nf = 0
    RHS[1:NC] = Xᶠ[1:NC,3]
    for n = 1:NC
        for j = 1:NS
            RHS[n] += β̃ᴹ[n,j,3]*Y′[n,j]
        end
        RHS[n] += β̃ᴹ[n,4,3]*(Xᶠ[n,4]-Dᴿ[n])
        res_static[nf+1] = Xˢ[n]/RHS[n]-1
        nf += 1
    end

    return res_static, Ŷ[1:NC,3], π′, ŵ, r̂, p̂
end

#=
    Subroutine 3
    - This function solves the whole dynamic problem.
    - Note that this function calls subroutine 2 multiple times, which calls subroutine 1
    multiple times.
=#
function dynamic_problem!(
    res_dynamic::AbstractArray,
    guess_dynamic::AbstractArray, # K̂[1:NC,1:NK,1] and Ŷ[1:NC,1:NS,1:T-1]
    init_dynamic::NamedTuple,
    exos_dynamic::NamedTuple,
    params_dynamic::NamedTuple,
    )
    # Unpack exogenous variables and parameters
    @unpack π₁, Y₁, Xᶠ₁, wL₁, rK₁ = init_dynamic
    @unpack Dᴿ, L̂, d̂, T̂ = exos_dynamic
    @unpack T, NC, NS, NK, β̃ᴸ, β̃ᴷ, ψ, θ, β̃ᴹ, ρ, δ, α = params_dynamic

    # Pre-allocate memory
    π = zeros(eltype(guess_dynamic),NC,NC,NS,T)
    Y = zeros(eltype(guess_dynamic),NC,NS,T)
    Xᶠ = zeros(eltype(guess_dynamic),NC,NS+1,T)
    wL = zeros(eltype(guess_dynamic),NC,T)
    rK = zeros(eltype(guess_dynamic),NC,NK,T)
    X = zeros(eltype(guess_dynamic),NC,NS,T)

    K̂ = ones(eltype(guess_dynamic),size(rK))
    Ŷ = ones(eltype(guess_dynamic),size(Y))
    X̂ᶠ = similar(Xᶠ)

    Ŷ_static = similar(Y)
    ŵ = similar(wL)
    r̂ = similar(rK)
    p̂ = similar(Y)

    # Assign initial conditions
    π[1:NC,1:NC,1:NS,1] = π₁
    Y[1:NC,1:NS,1] = Y₁[1:NC,1:NS]
    Xᶠ[1:NC,1:NS+1,1] = Xᶠ₁
    wL[1:NC,1] = wL₁
    rK[1:NC,1:NK,1] = rK₁
    for s in 1:NS
        X[1:NC,s,1] = π[1:NC,1:NC,s,1]'\(Y[1:NC,s,1])
    end

    # Resolve guess
    # guess_dynamic = reshape(guess_dynamic,NC,NK,T)
    K̂[1:NC,1:NK,1] = guess_dynamic[1:NC,1:NK,1]
    Ŷ[1:NC,1:NK,1:T-1] = guess_dynamic[1:NC,1:NK,2:T]

    # Evaluate Euler
    # println("Start to evaluate euler residuals.")
    for t = 1:T-1
        # Step 1
        # Calls subroutine 2
        myexos_static = @with_kw (
            π = π[1:NC,1:NC,1:NS,t],
            Ŷᴷ = Ŷ[1:NC,1:NK,t],
            Y = Y[1:NC,1:NS,t],
            Xᶠ = Xᶠ[1:NC,1:NS+1,t],
            Dᴿ = Dᴿ[1:NC,t+1],
            wL = wL[1:NC,t],
            L̂ = L̂[1:NC,t],
            rK = rK[1:NC,1:NK,t],
            K̂ = K̂[1:NC,1:NK,t],
            d̂ = d̂[1:NC,1:NC,1:NS,t],
            T̂ = T̂[1:NC,1:NS,t],
        )
        myparams_static = @with_kw (
            NC = NC, NS = NS, NK = NK, β̃ᴸ = β̃ᴸ, β̃ᴷ = β̃ᴷ, ψ = ψ, θ = θ, β̃ᴹ = β̃ᴹ
        )
        exos_static = myexos_static()
        params_static = myparams_static()
        guess_static = Ŷ[1:NC,3,t]

        # Solve the static problem
        # println(string("Start to solve the static problem in period ",t,"."))
        # println("Run time and memory cost:")
        #@time results_static =
        #    try
                results_static = nlsolve(
                    (res_static, guess_static) -> static_problem!(
                        res_static, guess_static, exos_static, params_static),
                    guess_static,
                    ftol=1e-12,
                    method=:anderson,
                    show_trace=false,
                )
        #    catch err
        #        if isa(err, DomainError)
        #            error("Failed to solve the static problem, please try again.")
        #        end
        #    end

        # Check Convergence
        # converged(results_static) || error("Failed to converge in $(results_static.iterations) iterations.")
        # println(string("Successfully solved the static problem in period ",t,".\n"))

        # Catch Solutions
        res_static = similar(results_static.zero)
        res_static, Ŷ_static[1:NC,3,t], π[1:NC,1:NC,1:NS,t+1], ŵ[1:NC,t], r̂[1:NC,1:NK,t], p̂[1:NC,1:NS,t] = static_problem!(
            res_static, results_static.zero, exos_static, params_static)

        # Update level variables Part I
        Ŷ[1:NC,3,t] = Ŷ_static[1:NC,3,t]
        for n = 1:NC
            wL[n,t+1] = wL[n,t]*ŵ[n,t]*L̂[n,t]
            for j = 1:NS
                Y[n,j,t+1] = Y[n,j,t]*Ŷ[n,j,t]
            end
            for k = 1:NK
                rK[n,k,t+1] = rK[n,k,t]*r̂[n,k,t]*K̂[n,k,t]
            end
        end
        X[1:NC,3,t+1] = π[1:NC,1:NC,3,t+1]'\(Y[1:NC,3,t+1])

        # Step 2
        # Solve for X̂ᶠ[:,1,t]
        X[1:NC,1,t+1] = Y[1:NC,1,t+1]
        for n = 1:NC
            X̂ᶠ[n,1,t] = Y[n,1,t+1]
            for j = 1:NS
                X̂ᶠ[n,1,t] -= β̃ᴹ[n,j,1]*Y[n,j,t+1]
            end
            X̂ᶠ[n,1,t] -= β̃ᴹ[n,4,1]*(Xᶠ[n,4,t]-Dᴿ[n,t+1])
            X̂ᶠ[n,1,t] /= Xᶠ[n,1,t]
        end

        # Step 3
        # Form Π[:,:,2,t+1] and get X[:,2,t+1]
        X[1:NC,2,t+1] = π[1:NC,1:NC,2,t+1]'\(Y[1:NC,2,t+1])

        # Step 4
        # Solve for X̂ᶠ[:,2,t]
        for n = 1:NC
            X̂ᶠ[n,2,t] = X[n,2,t+1]
            for j = 1:NS
                X̂ᶠ[n,2,t] -= β̃ᴹ[n,j,2]*Y[n,j,t+1]
            end
            X̂ᶠ[n,2,t] -= β̃ᴹ[n,4,2]*(Xᶠ[n,4,t]-Dᴿ[n,t+1])
            X̂ᶠ[n,2,t] /= Xᶠ[n,2,t]
        end

        # Step 5
        # Update level variables
            #= Update(d) levels(I)
                π[:,:,t+1] = π[:,:,t].*π̂[:,:,t]
                Y[:,:,t+1] = Y[:,:,t].*Ŷ[:,:,t]
                wL[:,t+1] = wL[:,t].*ŵ[:,t].*L̂[:,t]
                rK[:,:,t+1] = rK[:,:,t].*r̂[:,:,t].*K̂[:,:,t]
            =#
        # Update level variables(II)
        for n = 1:NC
            for j = 1:NK
                Xᶠ[n,j,t+1] = Xᶠ[n,j,t]*X̂ᶠ[n,j,t]
            end
            Xᶠ[n,3,t+1] = Xᶠ[n,3,t] #   *ϕ̂[n,t]*ψ̂[n,3,t]
            Xᶠ[n,4,t+1] = Xᶠ[n,4,t] #   *ϕ̂[n,t]*ψ̂[n,4,t]
        end

        # Step 6, 7
        # Use results in step 2 & 4 to evaluate Euler equation
        for n = 1:NC
            for k = 1:NK
                # Euler equation
                #=
                if p̂[n,k,t]*K̂[n,k,t]/X̂ᶠ[n,k,t]<0
                    println(n)
                    println(k)
                    println(t)
                    println(p̂[n,k,t])
                    println(K̂[n,k,t])
                    println(X̂ᶠ[n,k,t])
                    debug_t = t+1
                    colnames = [[string("X[1:NC,",i,",",debug_t,"]") for i in 1:3]; [string("Y[1:NC,",i,",",debug_t,"]") for i in 1:3];[string("XF[1:NC,",i,",",debug_t,"]") for i in 1:4]]
                    colnames = [Symbol(names) for names in colnames]
                    π_colnames = [[string("pi[1:NC,",i,",1,",debug_t,"]") for i in 1:21];[string("pi[1:NC,",i,",2,",debug_t,"]") for i in 1:21];[string("pi[1:NC,",i,",3,",debug_t,"]") for i in 1:21]]
                    π_colnames = [Symbol(names) for names in π_colnames]
                    CSV.write(string("debug_t",debug_t,".csv"),Tables.table([X[1:NC,1:NS,debug_t] Y[1:NC,1:NS,debug_t] Xᶠ[1:NC,1:NS+1,debug_t]];header=colnames))
                    CSV.write(string("pi_debug_t",debug_t,".csv"),Tables.table([π[1:NC,1:NC,1,debug_t] π[1:NC,1:NC,2,debug_t] π[1:NC,1:NC,3,debug_t]];header=π_colnames))
                end
                =#
                res_dynamic[n,k,t] = K̂[n,k,t]/(K̂[n,k,t]-(1-δ[k]))/ρ/
                    (α[k]*rK[n,k,t+1]/Xᶠ[n,k,t]+X̂ᶠ[n,k,t]*((1-α[k])+(p̂[n,k,t]*K̂[n,k,t]/X̂ᶠ[n,k,t])^α[k]*((1-δ[k])/(K̂[n,k,t]-(1-δ[k])))/χ̂[n,k,t]))-1.0
                # Update K̂ᵏₜ₊₁
                K̂[n,k,t+1] = χ̂[n,k,t]*(X̂ᶠ[n,k,t]/p̂[n,k,t]/K̂[n,k,t])^α[k]*(K̂[n,k,t]-(1-δ[k]))+(1-δ[k])
            end
        end

        # Step 8: Next iteration, from t -> T-1
    end
    # println("Successfully evaluate euler residuals.")

    # Step 8
    # Evaluate terminal conditions
    # println("Start to evaluate terminal conditions.")
    for n = 1:NC
        for k = 1:NK
            res_dynamic[n,k,T] = K̂[n,k,T]/1-1
        end
    end
    # println("Successfully evaluate terminal conditions.")
    #=
    colnames = [[string("Y_hat[1:NC,",i,",",T-1,"]") for i in 1:3]; [string("K_hat[1:NC,",i,",",T,"]") for i in 1:2]]
    colnames = [Symbol(names) for names in colnames]
    CSV.write(string("dynamic",T,".csv"),Tables.table([Ŷ[1:NC,1:NS,T-1] K̂[1:NC,1:NK,T]];header=colnames))
    =#
    return res_dynamic, K̂, Ŷ
end

ext: jl


text:
include("country.jl")
include("persona.jl")

    city_prefixes = ("Sur", "Norte",)
    city_adjetives = ("Nueva", "Vieja",)
    city_suffixies = ("de la Montaña", "los bajos", "los altos", )
    street_prefixes = ("Ampliación", "Andador", "Avenida", "Boulevard", "Calle", "Callejón","Calzada", "Cerrada", "Circuito", "Circunvalación", "Continuación","Corredor", "Diagonal", "Eje vial", "Pasaje", "Peatonal", "Periférico","Privada", "Prolongación", "Retorno", "Viaducto")
    building_number_formats = ("#####", "####", "###")
    postcode_formats = ("#####", "#####-####")

    ## States and abbrs from Mexico from INEGI
    ## http://www.inegi.org.mx/geo/contenidos/geoestadistica/CatalogoClaves.aspx
    states = (("AGS","Aguascalientes"), ("BC","Baja California"),("BCS","Baja California Sur"), ("CAMP","Campeche"),("COAH", "Coahuila de Zaragoza"), ("COL", "Colima"),("CHIS", "Chiapas"), ("CHIH", "Chihuahua"),("DF","Distrito Federal"), ("DGO", "Durango"),("GTO", "Guanajuato"), ("GRO", "Guerrero"), ("HGO", "Hidalgo"),("JAL", "Jalisco"), ("MEX", "México"),("MICH", "Michoacán de Ocampo"), ("MOR", "Morelos"),("NAY", "Nayarit"), ("NL", "Nuevo León"), ("OAX", "Oaxaca"),("PUE","Puebla"), ("QRO", "Querétaro"),("Q. ROO", "Quintana Roo"), ("SLP", "San Luis Potosí"),("SIN", "Sinaloa"), ("SON", "Sonora"), ("TAB", "Tabasco"),("TAMPS","Tamaulipas"), ("TLAX","Tlaxcala"),("VER", "Veracruz"),("YUC", "Yucatán"), ("ZAC","Zacatecas"),)

        ## The ZipCodes has a begin & final range
        ## Source: Norma Técnica de Domicilios INEGI
    zip_codes = Dict("AGS"=> (20000, 20999), "BC" => (21000, 22999),"BCS"=> (23000, 23999), "CAMP" => (24000, 24999),"COAH"=> (25000, 27999), "COL"=> (28000, 28999),"CHIS"=> (29000, 30999), "CHIH"=> (31000, 33999),"DF"=> (1000, 19999), "DGO"=> (36000, 35999),"GTO"=> (36000, 38999), "GRO"=> (39000, 41999),"HGO"=> (42000, 43999), "JAL"=> (44000, 49999),"MEX"=> (50000, 57999), "MICH"=> (58000, 61999),"MOR"=> (62000, 62999), "NAY"=> (63000, 63999),"NL"=> (64000, 67999), "OAX"=> (68000, 71999),"PUE"=> (72000, 75999), "QRO"=> (76000, 76999),"Q. ROO"=> (77000, 75999), "SLP"=> (78000, 79999),"SIN"=> (80000, 82999), "SON"=> (83000, 85999),"TAB"=> (86000, 86999), "TAMPS"=> (87000, 89999),"TLAX"=> (90000, 90999), "VER"=> (91000, 97999),"YUC"=> (97000, 97999), "ZAC"=> (98000, 99999))

    secondary_address_formats = ("### ###", "### Interior ###","### Edif. ### , Depto. ###")

city_prefix() = city_prefixes[rand(1:length(city_prefixes))]

city_suffix()= city_suffixies[rand(1:length(city_suffixies))]

city_adjetive() = city_adjetives[rand(1:length(city_adjetives))]

street_prefix()= street_prefixes[rand(1:length(street_prefixes))]

function secondary_address()

 Ρ = secondary_address_formats[rand(1:length(secondary_address_formats))]

  for i=Ρ
     Ρ= replace(Ρ,"#","$(rand(0:9))",1)
    end

 return Ρ

end

state() = states[rand(1:length(states))][2]

state_abbr() = states[rand(1:length(states))][1]

zip_code(p::AbstractString)= zip_codes[p][rand(1:2)]

street_name() = street_name_formats[rand(1:length(street_name_formats))]

street_address() = ("$(street_name()) $(secondary_address())")

address() = ("$(street_address()) $(city()), $(state_abbr()) $(postcode())")

function postcode()
     Ρ = postcode_formats[rand(1:length(postcode_formats))]

     for i=Ρ
     Ρ= replace(Ρ,"#","$(rand(0:9))",1)
     end

     return Ρ

end

function building_number()
     Ρ = building_number_formats[rand(1:length(building_number_formats))]

  for i=Ρ
     Ρ= replace(Ρ,"#","$(rand(0:9))",1)
    end

 return Ρ

 end

city() = city_formats[rand(1:length(city_formats))]

    city_formats = (" $(city_adjetive()) $(country())","San $(first_name()) $(city_suffix())")

    street_name_formats = ( "$(street_prefix()) $(last_name())","$(street_prefix()) $(country())","$(street_prefix()) $(state())","$(street_prefix()) $(city_prefix()) $(last_name())")

function geo_coordinate( center="None", radius=0.001)

        if center=="None"
            center=0
        end

        if center == 0
            return Float16(rand(-17.99 : 17.99))
        else
            return Float16(rand( center - radius : center + radius))
        end
end

    latitude() = geo_coordinate() / 2.0

    longitude() = geo_coordinate()/ 1.0


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

""" Typical properties for near eye VR display """
nominal_system_inputs() = (eye_relief = 20mm, fov = (5°,5°),eyebox = (10mm,8mm),display_radius = 125.0mm, pupil_diameter = 3.5mm,pixel_pitch = .9μm, minfnumber = 2.0, mtf = .2, cycles_per_degree = 11, max_display_size = 250μm, )
export nominal_system_inputs

xcoords(a::SMatrix{3,4}) = a[1,:]
ycoords(a::SMatrix{3,4}) = a[2,:]
zcoords(a::SMatrix{3,4}) = a[3,:]

function matnormal(a::SMatrix{3,4})
    side1 = a[:,4] - a[:,1]
    side2 = a[:,2] - a[:,1]
    return normalize(cross(side1,side2))
end

matcentroid(a::SMatrix{3,4}) = Statistics.mean(eachcol(a))

function matrix2rectangle(a::SMatrix{3,4,T}) where{T}
    center = Statistics.mean(eachcol(a))
    xmin,xmax = extrema(xcoords(a))
    halfsizeu = (xmax-xmin)/2.0
    ymin,ymax = extrema(ycoords(a))
    halfsizev = (ymax-ymin)/2.0
    Rectangle(halfsizeu,halfsizev,matnormal(a),center, interface = opaqueinterface())
end

function test_paraxial_lens()
    lens = ParaxialLensRect(10.0,1.0,1.0,[0.0,0.0,1.0],[0.0,0.0,0.0])
    displaypoint = [0.0,0.0,-2.5]
    direction = [0.0,0.0,1.0]
    r1 = Ray(displaypoint,direction)
    intsct1 = OpticSim.surfaceintersection(lens,r1)

    # compute the refracted ray and project this backwards to its intersection with the optical axis. This should be the same position as returned by virtualpoint()
    refrac1,_,_ = OpticSim.processintersection(OpticSim.interface(lens),OpticSim.point(intsct1),OpticSim.normal(lens),OpticalRay(r1,1.0,.55),OpticSim.TEMP_REF,OpticSim.PRESSURE_REF,false)
    r2 = Ray(-displaypoint,-direction)

    intsct2= OpticSim.surfaceintersection(lens,r2)

    # compute the refracted ray and project this backwards to its intersection with the optical axis. This should be the same position as returned by virtualpoint()
    refrac2,_,_ = OpticSim.processintersection(OpticSim.interface(lens),OpticSim.point(intsct2),OpticSim.normal(lens),OpticalRay(r2,1.0,.55),OpticSim.TEMP_REF,OpticSim.PRESSURE_REF,false)

    # isapprox(-refrac2[1],refrac1[1]), refrac2[1], refrac1[1]
    Vis.draw(lens)
    Vis.draw!([Ray(OpticSim.point(intsct1)+ [.5,.5,0.0],refrac1),r1]) #offset first ray so it can be seen. The other ray will be written over it and obscure it.
    Vis.draw!([Ray(OpticSim.point(intsct2),refrac2),r2], color = "green") #visualization should have two rays in opposite directions coming out of lens plane but only see one. Error in ParaxialLens.
end
export test_paraxial_lens

"""Example that shows how to call setup_system with typical values"""
function setup_nominal_system()::LensletSystem
    (;eye_relief,fov,eyebox,display_radius,pupil_diameter,minfnumber,pixel_pitch) = nominal_system_inputs()
    setup_system(eyebox,fov,eye_relief,pupil_diameter,display_radius,minfnumber,pixel_pitch)
end
export setup_nominal_system

function compute_eyebox_rays(sys = setup_nominal_system())
    (;lenses,projected_eyeboxes,eyebox_rectangle) = sys
    centers = [x for x in centroid.(shape.(lenses))] #use the geometric center of the lenses not the optical center_point
    display_centers = matcentroid.(projected_eyeboxes)
    rays = Ray.(display_centers, centers.-display_centers)
    ray_traces = Vector{LensTrace{Float64,3}}(undef,0)

    tracked_rays = [trace(CSGOpticalSystem(LensAssembly(lens),eyebox_rectangle),OpticalRay(ray,1.0,.530),trackrays = ray_traces) for (lens,ray) in zip(lenses,rays)]

    return [x for x in ray_traces]
end
export compute_eyebox_rays

"""Example call of system_properties function"""
function systemproperties_call()
    (;eye_relief,fov,eyebox,pupil_diameter,minfnumber,pixel_pitch) = nominal_system_inputs()
    system_properties(eye_relief,eyebox,fov,pupil_diameter,.2,20,minfnumber = minfnumber,pixelpitch = pixel_pitch)
end
export systemproperties_call

"""Computes the lenslets on the display sphere and draws them."""
function drawspherelenslets()
    (;fov, eye_relief,eyebox,display_radius,pupil_diameter,pixel_pitch,minfnumber,mtf,cycles_per_degree) = nominal_system_inputs()

    computedprops = system_properties(eye_relief,eyebox,fov,pupil_diameter,mtf,cycles_per_degree,minfnumber = minfnumber,maxdisplaysize = 250μm,pixelpitch = pixel_pitch)
    focallength = ustrip(mm,computedprops[:focal_length])
    lenses,_ = spherelenslets(Plane(0.0,0.0,1.0,0.0,0.0,18.0),eye_relief,focallength,[0.0,0.0,1.0],display_radius,fov[1],fov[2],HexBasis1())
    Vis.draw() #clear screen
    Vis.draw!.(lenses)
    return nothing
end
export drawspherelenslets

"""returns the hex polygons on the sphere that will eventually be turned into hexagonal lenslets"""
function testspherepolygons()
    eyebox = Plane(0.0,0.0,1.0,0.0,0.0,12.0)
    dir = [0.0,0.0,1.0]
    fovθ = 10.0°
    fovϕ = 5.0°
    lattice = HexBasis1()
    displayradius = 125.0mm
    eyerelief = 20mm

    return spherepolygons(eyebox,eyerelief,displayradius,dir,fovθ,fovϕ,lattice)
end
export testspherepolygons

function testprojectonplane()
    verts = tilevertices(HexBasis1())
    verts = vcat(verts,[0 for _ in 1:6]')
    projectonbestfitplane(SMatrix{size(verts)...}(verts),[0.0,0.0,1.0])
end
export testprojectonplane

function testproject()
    norml = SVector(0.0, 0, 1)
    hex = HexBasis1()
    verts =vcat(tilevertices((0, 0), hex), [0.0 0 0 0 0 0])
    surf = Plane(norml, SVector(0.0, 0, 10))

    project(verts, norml, surf)
end
export testproject

function testprojection()
    pts = spherepoints(1.0,-.2,-.2,1.0,1.1)
    surf = Plane(0.0,0.0,-1.0,0.0,0.0,0.0)
    dir = [0.0,0.0,-1.0]
    project(pts,dir,surf)
end
export testprojection

"""returns hexagonal lenslets on the display sphere"""
function testspherelenslets()
    (;fov, eye_relief,eyebox,display_radius,pupil_diameter,pixel_pitch,minfnumber,mtf,cycles_per_degree,max_display_size) = nominal_system_inputs()

    fov = (10°,10°)
    computedprops = system_properties(eye_relief,eyebox,fov,pupil_diameter,mtf,cycles_per_degree,minfnumber = minfnumber,maxdisplaysize = max_display_size,pixelpitch = pixel_pitch)
    focallength = ustrip(mm,computedprops[:focal_length])
    spherelenslets(Plane(0.0,0.0,1.0,0.0,0.0,18.0),eye_relief,focallength,[0.0,0.0,1.0],display_radius,fov[1],fov[2],HexBasis1())
end
export testspherelenslets

function draw_projected_corners(system = setup_nominal_system())
    (;lenslet_eye_boxes, subdivisions_of_eyebox,lenses,displayplanes,projected_eyeboxes) = system
    colors = distinguishable_colors(reduce(*,subdivisions_of_eyebox))

    for (lenslet_eyebox,lens,display_plane,projected_eyebox) in zip(lenslet_eye_boxes,lenses,displayplanes,projected_eyeboxes)
        center = opticalcenter(lens)
        for (eyeboxpt,projected_point) in zip(eachcol(lenslet_eyebox),eachcol(projected_eyebox))
            r = Ray(eyeboxpt,center-eyeboxpt)
            intsct = surfaceintersection(display_plane,r)
            closest = closestintersection(intsct)
            display_intersection = point(closest)
            @assert isapprox(display_intersection,projected_point) "error display_intersection $display_intersection projected_point $projected_point"
            lenstrace = LensTrace(OpticalRay(r,1.0,.5),closest)
            Vis.draw!(lenstrace)
        end
    end
end
export draw_projected_corners

"""assigns each lenslet/display subsystem a rectangular sub part of the eyebox"""
function draw_eyebox_assignment(system = setup_nominal_system(),clear_screen = true;draw_eyebox = true)
    (;eyebox_rectangle,
    lenses,
    lenslet_colors,
    lattice_coordinates,
    subdivisions_of_eyebox,
    projected_eyeboxes,
    displayplanes,
    lenslet_eyebox_numbers) = system
    if clear_screen
        Vis.draw() #clear screen
    end

    if draw_eyebox
        Vis.draw!(eyebox_rectangle)
    end

    for (lens,lattice_coordinates,color) in zip(lenses,lattice_coordinates,lenslet_colors)
        Vis.draw!(lens,color = color)
    end

    num_distinct_eyeboxes = reduce(*,subdivisions_of_eyebox)
    colors = distinguishable_colors(num_distinct_eyeboxes)

    #THIS is almost certainly wrong
    # for (eyebox,plane,eyeboxnum,lens) in zip(projected_eyeboxes,displayplanes,lenslet_eyebox_numbers,lenses)
    #     localframe = OpticSim.translation(-OpticSim.normal(lens))*OpticSim.localframe(shape(lens)) #the local frame of the lens is, unfortunately, stored only in the ConvexPoly shape. No other lens type has a local frame which is terrible design. Should be fixed.
    #     transformedpoints = (inv(localframe)*eyebox)[1:2,:]
    #     pts = [SVector{2}(x...) for x in eachcol(transformedpoints)]
    #     # Vis.draw!(plane,color = colors[eyeboxnum])
    #     Vis.draw!(ConvexPolygon(localframe,pts),color = colors[eyeboxnum])
    # end
end
export draw_eyebox_assignment

#compute the ray from each lenslet's eyebox center to the geometric center of the lenslet. Trace this ray through the lens and draw the refracted ray. The refracted ray should line up exactly with -normal(lenslet).
function draw_eyebox_rays(system = setup_nominal_system())
    (;lenses, lenslet_eyebox_centers,) = system

    for (lens,eyebox_center) in zip(lenses,lenslet_eyebox_centers)
        r = Ray(centroid(lens),centroid(lens)-eyebox_center)
        refrac_ray = ray(trace(LensAssembly(lens),OpticalRay(r,1.0,.5)))
        Vis.draw!(refrac_ray,color = "yellow")
    end
end

function draw_eyebox_polygons(system = setup_nominal_system())
    (;projected_eyebox_polygons) = system
    for poly in projected_eyebox_polygons
        Vis.draw!(poly)
    end
end
export draw_eyebox_polygons

function draw_system(system = setup_nominal_system())
    draw_eyebox_assignment(system,draw_eyebox=false)
    draw_subdivided_eyeboxes(system,false)
    # Vis.draw!(compute_eyebox_rays(system))
    draw_eyebox_rays(system)
    draw_projected_corners(system)
    draw_eyebox_polygons(system)
end
export draw_system

"""draw subdivided eyeboxes to make sure they are in the right place"""
function draw_subdivided_eyeboxes(system = setup_nominal_system, clear_screen = true)
    (;subdivided_eyebox_polys) = system

    colors = distinguishable_colors(length(subdivided_eyebox_polys))

    if clear_screen
        Vis.draw()
    end

    for (eyebox,color) in zip(subdivided_eyebox_polys,colors)
        center = matcentroid(eyebox)
        Vis.draw!(matrix2rectangle(eyebox),color = color)
        r = Ray(center,matnormal(eyebox))
        Vis.draw!(r)
        Vis.draw!(leaf(Sphere(.1),Geometry.translation(center)),color = "white")
    end
end
export draw_subdivided_eyeboxes

function test_project_eyebox_to_display_plane()
    boxpoly = SMatrix{3,4}(
        -1.0,1,0.0,
        1.0,1.0,0.0,
        1.0,-1.0,0.0,
        -1.0,-1.0,0.0
        )
    correct_answer = SMatrix{3,4}(
        1.3,1.0, 11.5,
        1.0,1.0, 11.5,
        1.0,1.3, 11.5,
        1.3,1.3, 11.5)
    fl = 1.5
    lenscenter = [1.0,1.0,10.0]
    lens = ParaxialLensRect(fl,.5,.5,[0.0,0.0,1.0],lenscenter)
    displayplane = Plane([0.0,0.0,1.0],[0.0,0.0,lenscenter[3]+fl])

    answer,polygons = project_eyebox_to_display_plane(boxpoly,lens,displayplane)
    @assert isapprox(answer,correct_answer) "computed points $answer"
    @info "eyebox points $answer"
    @info "eyebox polygons $polygons"
end
export test_project_eyebox_to_display_plane

"""For each display computes a ray from the display center to the geometric center of the lenslet. Then it computes the refracted ray and intersects it with the subdivided eyebox associated with the display. Draws all these rays. The rays from the displays should hit the center points of the subdivided eyeboxes"""
function test_ray_cast_from_display_center(clear_display = true)
    system = setup_nominal_system()
    (;eyebox_rectangle,lenses,projected_eyeboxes) = system

    if clear_display
        Vis.draw()
    end

    # draw_system(system)
    draw_subdivided_eyeboxes(system,false)

    for (lens,display_eyebox) in zip(lenses,projected_eyeboxes)
        lens_geometric_center = centroid(lens)

        display_center = matcentroid(display_eyebox)
        diff = lens_geometric_center-display_center

        r = OpticalRay(Ray(display_center,diff),1.0,.53)

        # Vis.draw!(eyebox_rectangle)
        system = CSGOpticalSystem(LensAssembly(lens),eyebox_rectangle)
        Vis.draw!(r,color = "black")
        tr = trace(system,r)

        if tr !== nothing
            Vis.draw!(tr,color = "black")
        end
    end
end
export test_ray_cast_from_display_center

"""refracts the ray from the center of the display passing through the geometric center of the lens and then computes the closes point on this ray to the eyeboxcentroid. If the lens optical center has been calculated correctly the ray should pass through the eyeboxcentroid"""
function test_replace_optical_center(eyeboxcentroid,lensgeometriccenter,displaycenter,lens)
    rdir = lensgeometriccenter-displaycenter
    r = Ray(displaycenter,rdir)

    assy = LensAssembly(lens)
    refracted_trace = trace(assy,OpticalRay(r,1.0,.530))

    r2 = Ray(eyeboxcentroid,lensgeometriccenter-eyeboxcentroid)

    #this doesn't work because of a bug in processintersection for paraxial interfaces. Ray points in the wrong direction.
    reverse_trace = trace(assy,OpticalRay(r2,1.0,.530))
    reverse_point = closestpointonray(ray(reverse_trace),displaycenter)

    return closestpointonray(ray(refracted_trace),eyeboxcentroid)
end
export test_replace_optical_center

function center_data()
    a = HexBasis1()
    verts = 5*tilevertices(a)
    vecofmat = collect(reinterpret(reshape,SVector{2,eltype(verts)},verts))

    lens = ParaxialLensConvexPoly(
        1.0,
        Geometry.translation(1.0,0.0,1.0)*Geometry.rotation(0.0,Float64(π),0.0),
        vecofmat,
        SVector(0.0,0.0)
    )
    display_center = [1.0,0.0,2.0]
    eyebox_center = [0.0,0.0,0.0]
    return lens,eyebox_center,display_center
end
export center_data

"""verify that paraxial lens properly transorms rays into local coordinate frame to compute intersections with convex poly shape"""
function test_paraxial_convex_poly()
    lens,eyebox_center,display_center = center_data()

    @info "normal of lens $(OpticSim.normal(lens))"
    @info "center of lens $(centroid(lens)) center of lens geometry $(centroid(shape(lens)))"

    display_ray = OpticalRay(Ray([1.0,0.0,2.0],[0.0,0.0,-1.0]),1.0,.53)

    ltrace = trace(LensAssembly(lens),display_ray)
    @info "trace with optical center aligned with geometric center $ltrace"
    new_center = compute_optical_center(eyebox_center,display_center,lens)
    @info "offset center $new_center"
    displaced_lens = replace_optical_center(eyebox_center,display_center,lens)
    @info "optical center of displacedlens $(opticalcenter(displaced_lens)) center of geometry $(centroid(shape(displaced_lens)))"

    dtrace = trace(LensAssembly(displaced_lens),display_ray)
    @info "trace with optical center displaced from geometric center  $dtrace"
    refracted_ray = ray(dtrace)
    intsct = surfaceintersection(Plane([0.0,0.0,1.0],eyebox_center),refracted_ray)

    @info "intsct should be [0,0,0] $(point(OpticSim.lower(intsct)))" #ray was computed to refract from the geometric lens center to the point [0,0,0] on the eyebox plane
    return intsct

end
export test_paraxial_convex_poly


text:
module DiffieHellman

using CryptoGroups

struct DH
    wrap::Function
    unwrap::Function
    G::AbstractGroup
    hash::Function
    rngint::Function
end

### It is actually possible for your own socket to define stack and unstack method in the externa libraries. For that however I need a common dependency. Stackers.jl?

function stack(io::IO,msg::Vector{UInt8})
    frontbytes = reinterpret(UInt8,Int16[length(msg)])
    item = UInt8[frontbytes...,msg...]
    write(io,item)
end

function unstack(io::IO)
    sizebytes = [read(io,UInt8),read(io,UInt8)]
    size = reinterpret(Int16,sizebytes)[1]

    msg = UInt8[]
    for i in 1:size
        push!(msg,read(io,UInt8))
    end
    return msg
end

function diffiehellman(io::IO,wrap::Function,unwrap::Function,G::AbstractGroup,hash::Function,a::Integer)
    Avalue = value(G^a)
    envelopeA = wrap(Avalue)
    stack(io,envelopeA)

    envelopeB = unstack(io)
    Bvalue,id = unwrap(envelopeB)

    B = typeof(G)(Bvalue,G) ### I could also have used a type parameter
    @assert B!=G "Trivial group elements are not allowed."
    key = value(B^a)

    cmsgA = hash(envelopeA,envelopeB,key)
    stack(io,cmsgA)

    cmsgB = unstack(io)
    @assert cmsgB==hash(envelopeB,envelopeA,key) "The key exchange failed."

    return key,id
end

diffiehellman(io::IO,dh::DH) = diffiehellman(io,dh.wrap,dh.unwrap,dh.G,dh.hash,dh.rngint())

export diffiehellman, DH

end # module


text:
struct StrainData
  outside_multiplier::Float32
  household_multiplier::Float32
end

const StrainTable = SVector{3, StrainData}

function make_strains(;base_multiplier::Real=1.0, british_multiplier::Real=1.70, delta_multiplier::Real=1.7*1.5)::StrainTable
  SA[
    StrainData(base_multiplier, base_multiplier),
    StrainData(base_multiplier * british_multiplier, base_multiplier * british_multiplier),
    StrainData(base_multiplier* delta_multiplier, base_multiplier * delta_multiplier)
  ]
end

getdata(table::StrainTable, strain::StrainKind) = table[UInt8(strain)]


text:
#
# C interface wrappers for FGb
#
libfgb_handle = Libdl.dlopen(joinpath(@__DIR__, "..", "deps","libfgb.so"))
macro libfgb(symbol)
    ptr = Libdl.dlsym(libfgb_handle, symbol)
    quote
        $ptr
    end
end

unsafe_store!(
    cglobal(@libfgb(FGb_julia_info), Ptr{Void}),
    cfunction(
        info∘unsafe_string,
        Void,
        (Cstring,)
    ),
)
unsafe_store!(
    cglobal(@libfgb(FGb_julia_error), Ptr{Void}),
    cfunction(
        error∘unsafe_string,
        Void,
        (Cstring,)
    ),
)

enter_INT() = ccall(@libfgb(FGb_int_enter_INT), Void, ())
init_urgent(bytes_coeff, bytes_exponent, cargo_cult_drldrl, cargo_cult_100000, cargo_cult_0) = ccall(@libfgb(FGb_int_init_urgent), Void, (UInt32,UInt32,Cstring,UInt32,UInt32), bytes_coeff, bytes_exponent, cargo_cult_drldrl, cargo_cult_100000, cargo_cult_0)
log_output = Ref{UInt}(0)
init(cargo_cult_1, cargo_cult_1_again, cargo_cult_0, logging_file_handle) = ccall(@libfgb(FGb_int_init), Void, (UInt32,UInt32,UInt32,Ref{UInt}), cargo_cult_1, cargo_cult_1_again, cargo_cult_0, logging_file_handle)
reset_coeffs(cargo_cult_1) = ccall(@libfgb(FGb_int_reset_coeffs), Void, (UInt32,), cargo_cult_1)
reset_expos(nvars1, nvars2, varnames) = ccall(@libfgb(FGb_int_reset_expos), Void, (UInt32,UInt32,Ref{Cstring}), nvars1, nvars2, varnames)
internal_threads(n) = ccall(@libfgb(FGb_int_internal_threads), Void, (UInt32,), n)
mod_internal_threads(n) = ccall(@libfgb(FGb_internal_threads), Void, (UInt32,), n)

creat_poly(num_terms) = ccall(@libfgb(FGb_int_creat_poly), Ptr{Void}, (UInt32,), num_terms)
exp = UInt32[0,1,1,0,0]
set_expos2(poly, term, exp, explen) = ccall(@libfgb(FGb_int_set_expos2), Void, (Ptr{Void}, UInt32, Ref{UInt32},UInt32), poly, term, exp, explen)
set_coeff_gmp(poly, term, coeff) = ccall(@libfgb(FGb_int_set_coeff_gmp), Void, (Ptr{Void}, UInt32, Ref{BigInt}), poly,term,coeff)
full_sort_poly2(poly) = ccall(@libfgb(FGb_int_full_sort_poly2), Void, (Ptr{Void},), poly)

mutable struct SFGB_Comp_Desc
    _compute::UInt32   # in C, an enum
    _nb::Int32
    _force_elim::Int32
    _off::UInt32
    _index::UInt32
    _zone::UInt32
    _memory::UInt32

    _nb2::Int32
    _force_elim2::Int32
    _bk2::UInt32
    _aggressive2::Int32
    _dlim::Int32
    _skip::Int32
    SFGB_Comp_Desc() = begin
        env = new()
        env._compute=1 #FGB_COMPUTE_GBASIS
        env._nb=0
        env._force_elim=0 # if force_elim=1 then return only the result of the elimination
                          # (need to define a monomial ordering DRL(k1,k2) with k2>0 )
        env._off=0        # use to select another subset of prime numbers
        env._index=5000000 # This is is the maximal size of the matrices generated by F4
                          # you can increase this value according to your memory
        env._zone=0       # should be 0
        env._memory=0     # should be 0
        env
    end
end

groebner(input_basis, input_basis_length, output_basis, _mini, _elim, cputime, _bk0, _step0, _elim0, _env) = ccall(@libfgb(FGb_int_groebner), UInt32,
       (Ref{Ptr{Void}},UInt32,Ref{Ptr{Void}},UInt32,UInt32,Ref{Float64},UInt32, Int32, UInt32, Ref{SFGB_Comp_Desc}),
       input_basis, input_basis_length, output_basis, _mini, _elim, cputime, _bk0, _step0, _elim0, _env)

nb_terms(p) = ccall(@libfgb(FGb_int_nb_terms), UInt32, (Ptr{Void},), p)
export_poly_INT_gmp2(n_variables, n_terms, coefficients, exponents, p) = ccall(@libfgb(FGb_int_export_poly_INT_gmp2), UInt32, (UInt32, UInt32, Ref{Ptr{BigInt}}, Ref{UInt32}, Ptr{Void}), n_variables, n_terms, coefficients, exponents, p)

reset_memory() = ccall(@libfgb(FGb_int_reset_memory), Void, ())
exit_INT() = ccall(@libfgb(FGb_int_exit_INT), Void, ())


text:
"""
    power_table(sim, alpha = 0.05)

Returns a row table with fields `coefname` and `power`  based on the proportion of
simulated p-values less than alpha, for `sim`, the output of `parametricbootstrap`.
"""
function power_table(sim::MixedModelBootstrap, alpha = 0.05)
# `se` represents the Bernoulli standard error of the number of trials below signifance using the (1 - alpha) for the probability.
# q = 1 - alpha
# se = sqrt(alpha * q / nsim)
    nsim = length(sim.objective)
    dd = Dict{Symbol, Int}()
    for row in sim.coefpvalues
        coef = row.coefname
        # make sure every even zero-power stuff is initialized
        val = get!(dd, coef, 0)
        if row.p < alpha
            dd[coef] = val + 1
        end
    end
    return [ (; :coefname => string(key), :power => val / nsim) for (key, val) in dd ]
end


text:
include("dependencies_for_runtests.jl")

@testset "Oceananigans" begin
    if test_file != :none
        @testset "Single file test" begin
            include(String(test_file))
        end
    end

    # Core Oceananigans
    if group == :unit || group == :all
        @testset "Unit tests" begin
            include("test_grids.jl")
            include("test_operators.jl")
            include("test_boundary_conditions.jl")
            include("test_field.jl")
            include("test_field_reductions.jl")
            include("test_halo_regions.jl")
            include("test_coriolis.jl")
            include("test_buoyancy.jl")
            include("test_stokes_drift.jl")
            include("test_utils.jl")
            include("test_schedules.jl")
        end
    end

    if group == :abstract_operations || group == :all
        @testset "AbstractOperations and broadcasting tests" begin
            include("test_abstract_operations.jl")
            include("test_conditional_reductions.jl")
            include("test_computed_field.jl")
            include("test_broadcasting.jl")
        end
    end

    if group == :poisson_solvers_1 || group == :all
        @testset "Poisson Solvers 1" begin
            include("test_poisson_solvers.jl")
        end
    end

    if group == :poisson_solvers_2 || group == :all
        @testset "Poisson Solvers 2" begin
            include("test_poisson_solvers_vertically_stretched_grid.jl")
        end
    end

    if group == :matrix_poisson_solvers || group == :all
        @testset "Matrix Poisson Solvers" begin
            include("test_matrix_poisson_solver.jl")
        end
    end

    if group == :general_solvers || group == :all
        @testset "General Solvers" begin
            include("test_batched_tridiagonal_solver.jl")
            include("test_preconditioned_conjugate_gradient_solver.jl")
        end
    end

    # Simulations
    if group == :simulation || group == :all
        @testset "Simulation tests" begin
            include("test_simulations.jl")
            include("test_diagnostics.jl")
            include("test_output_writers.jl")
            include("test_output_readers.jl")
            include("test_lagrangian_particle_tracking.jl")
        end
    end

    # Models
    if group == :time_stepping_1 || group == :all
        @testset "Model and time stepping tests (part 1)" begin
            include("test_nonhydrostatic_models.jl")
            include("test_time_stepping.jl")
        end
    end

    if group == :time_stepping_2 || group == :all
        @testset "Model and time stepping tests (part 2)" begin
            include("test_boundary_conditions_integration.jl")
            include("test_forcings.jl")
        end
    end

    if group == :time_stepping_3 || group == :all
        @testset "Model and time stepping tests (part 3)" begin
            include("test_dynamics.jl")
        end
    end

    if group == :turbulence_closures || group == :all
        @testset "Turbulence closures tests" begin
            include("test_turbulence_closures.jl")
        end
    end

    if group == :shallow_water || group == :all
        include("test_shallow_water_models.jl")
    end

    if group == :hydrostatic_free_surface || group == :all
        @testset "HydrostaticFreeSurfaceModel tests" begin
            include("test_hydrostatic_free_surface_models.jl")
            include("test_ensemble_hydrostatic_free_surface_models.jl")
            include("test_hydrostatic_free_surface_immersed_boundaries.jl")
            include("test_vertical_vorticity_field.jl")
            include("test_implicit_free_surface_solver.jl")
            include("test_split_explicit_free_surface_solver.jl")
            include("test_split_explicit_vertical_integrals.jl")
            include("test_hydrostatic_free_surface_immersed_boundaries_implicit_solve.jl")
        end
    end

    # Model enhancements: cubed sphere, distributed, etc
    if group == :multi_region || group == :all
        @testset "Multi Region tests" begin
            include("test_multi_region_unit.jl")
            include("test_multi_region_advection_diffusion.jl")
            include("test_multi_region_implicit_solver.jl")
        end
    end

    if group == :cubed_sphere || group == :all
        @testset "Cubed sphere tests" begin
            include("test_cubed_spheres.jl")
            include("test_cubed_sphere_halo_exchange.jl")
            include("test_cubed_sphere_circulation.jl")
        end
    end

    if group == :distributed || group == :all
        MPI.Initialized() || MPI.Init()
        include("test_distributed_models.jl")
        include("test_distributed_poisson_solvers.jl")
    end

    if group == :nonhydrostatic_regression || group == :all
        include("test_nonhydrostatic_regression.jl")
    end

    if group == :hydrostatic_regression || group == :all
        include("test_hydrostatic_regression.jl")
    end

    if group == :scripts || group == :all
        @testset "Scripts" begin
            include("test_validation.jl")
        end
    end

    if group == :convergence
        include("test_convergence.jl")
    end
end


text:
mutable struct SerialPort <: IO
    ref::Port
    eof::Bool
    open::Bool
    function SerialPort(ref, eof, open)
        sp = new(ref, eof, open)
        finalizer(destroy!, sp)
        return sp
    end
end

"""
`sp = SerialPort(portname::AbstractString)`

Constructor for the `SerialPort` object.
"""
SerialPort(portname::AbstractString) = SerialPort(sp_get_port_by_name(portname), false, false)

"""
`destroy!(sp::SerialPort)`

Destructor for the `SerialPort` object.
"""
function destroy!(sp::SerialPort)
    close(sp)
    sp_free_port(sp.ref)
end

"""
`set_speed(sp::SerialPort,bps::Integer)`

Set connection speed of `sp` in bits per second. The library will return an
error if bps is not a valid/supported value.
"""
function set_speed(sp::SerialPort, bps::Integer)
     sp_set_baudrate(sp.ref, bps)
     return nothing
end

"""
`set_frame(sp::SerialPort [, ndatabits::Integer, pariry::SPParity, nstopbits::Integer])`

Configure packet framing. Defaults to the most common "8N1" scheme. See
https://en.wikipedia.org/wiki/Universal_asynchronous_receiver/transmitter#Data_framing
for more details.

`ndatabits` is the number of data bits which is `8` in the common "8N1" sceme.

The `parity` is set to none in the "8N1" sceme and can take the values:
`SP_PARITY_NONE`, `SP_PARITY_ODD`, `SP_PARITY_EVEN`, `SP_PARITY_MARK` and
`SP_PARITY_SPACE`.

`nstopbits` is the number of stop bits, which is `1` by default.
"""
function set_frame(sp::SerialPort;
    ndatabits::Integer=8,
    parity::SPParity=SP_PARITY_NONE,
    nstopbits::Integer=1)

    sp_set_bits(sp.ref, ndatabits)
    sp_set_parity(sp.ref, parity)
    sp_set_stopbits(sp.ref, nstopbits)
    return nothing
end

"""
`set_flow_control(sp::SerialPort [,rts::SPrts, cts::SPcts, dtr::SPdtr, dst::SPdsr, xonxoff::SPXonXoff])`

Configure flow control settings. Many systems don't support all options.
If an unsupported option is requested, the library will return SP_ERR_SUPP.

`rts` can take the values: `SP_RTS_OFF`, `SP_RTS_ON` and `SP_RTS_FLOW_CONTROL`
and defaults to `SP_RTS_OFF`.

`cts` can take the values: `SP_CTS_IGNORE` and `SP_CTS_FLOW_CONTROL`. Its
default is `SP_CTS_IGNORE`.

`dtr` can take the values: `SP_DTR_OFF`, `SP_DTR_ON`, and `SP_DTR_FLOW_CONTROL`
and defaults to `SP_DTR_OFF`.

`dsr` can take the values: `SP_DSR_IGNORE` and `SP_DSR_FLOW_CONTROL`. Its
default is SP_DSR_IGNORE`.

`xonxoff` can take values: `SP_XONXOFF_DISABLED`, `SP_XONXOFF_IN`,
`SP_XONXOFF_OUT`, and `SP_XONXOFF_INOUT` and defaults to `SP_XONXOFF_DISABLED`.
"""
function set_flow_control(sp::SerialPort;
    rts::SPrts=SP_RTS_OFF,
    cts::SPcts=SP_CTS_IGNORE,
    dtr::SPdtr=SP_DTR_OFF,
    dsr::SPdsr=SP_DSR_IGNORE,
    xonxoff::SPXonXoff=SP_XONXOFF_DISABLED)

    # CTS and RTS must be enabled or disabled as a pair
    sp_set_rts(sp.ref, rts)
    sp_set_cts(sp.ref, cts)

    sp_set_dtr(sp.ref, dtr)
    sp_set_dsr(sp.ref, dsr)

    sp_set_xon_xoff(sp.ref, xonxoff)
    return nothing
end

"""
`listports([nports_guess::Integer])`

Print a list of currently visible ports, along with some basic info.

`nports_guess` provides the number of ports guessed. Its default is `64`.
"""
function list_ports(;nports_guess::Integer=64)
    ports = sp_list_ports()

    for port in unsafe_wrap(Array, ports, nports_guess, own=false)
        port == C_NULL && return

        println(sp_get_port_name(port))
        println("\tDescription:\t",    sp_get_port_description(port))
        println("\tTransport type:\t", sp_get_port_transport(port))
    end

    sp_free_port_list(ports)
    return nothing
end

"""
`print_port_metadata(sp::SerialPort [,show_config::Bool])

Print info found for this port.
Note: port should be open to obtain a valid FD/handle before accessing fields.

`show_config` is `true` by default and prints out the current port settings.
"""
function print_port_metadata(sp::SerialPort; show_config::Bool=true)
    print_port_metadata(sp.ref, show_config=show_config)
    return nothing
end

function print_port_metadata(port::LibSerialPort.Port; show_config::Bool=true)
    println("\nPort name:\t",       sp_get_port_name(port))
    println("Manufacturer:\t",      sp_get_port_usb_manufacturer(port))
    println("Product:\t",           sp_get_port_usb_product(port))
    println("USB serial number:\t", sp_get_port_usb_serial(port))
    println("Bluetooth address:\t", sp_get_port_bluetooth_address(port))
    println("File descriptor:\t",   sp_get_port_handle(port))

    bus, addr = sp_get_port_usb_bus_address(port)
    if bus != -1
        println("USB bus #:\t",   bus)
        println("Address on bus:\t",  addr)
    end

    vid, pid = sp_get_port_usb_vid_pid(port)
    if vid != -1
        println("Vendor ID:\t",   vid)
        println("Product ID:\t",  pid)
    end
    if show_config
        print_port_settings(port)
    end
    return nothing
end

"""
`print_port_settings(sp::SerialPort)`

Print port settings for `sp`.
"""
print_port_settings(sp::SerialPort) = print_port_settings(sp.ref)

function print_port_settings(port::LibSerialPort.Port)
    println("Configuration for ", sp_get_port_name(port), ":")
    config = sp_get_config(port)
    print_port_settings(config)
    sp_free_config(config)
end

function print_port_settings(config::LibSerialPort.Config)
    println("\tbaudrate\t", sp_get_config_baudrate(config))
    println("\tbits\t",     sp_get_config_bits(config))
    println("\tparity\t",   sp_get_config_parity(config))
    println("\tstopbits\t", sp_get_config_stopbits(config))
    println("\tRTS\t",      sp_get_config_rts(config))
    println("\tCTS\t",      sp_get_config_cts(config))
    println("\tDTR\t",      sp_get_config_dtr(config))
    println("\tDSR\t",      sp_get_config_dsr(config))
    println("\tXonXoff\t",  sp_get_config_xon_xoff(config))
    println("")
end

"""
`open(sp::SerialPort [, mode::SPMode])`

Open the serial port `sp`.

`mode` can take the values: `SP_MODE_READ`, `SP_MODE_WRITE`, and
`SP_MODE_READ_WRITE`
"""
function Base.open(sp::SerialPort; mode::SPMode=SP_MODE_READ_WRITE)
    sp_open(sp.ref, mode)
    sp.open = true
    return sp
end

"""
`open(portname::AbstractString,baudrate::Integer [,mode::SPMode,
    ndatabits::Integer,parity::SPParity,nstopbits::Integer])`

construct, configure and open a `SerialPort` object.

For details on posssible settings see `?set_flow_control` and `?set_frame`.
"""
function Base.open(portname::AbstractString,
                   bps::Integer;
                   mode::SPMode=SP_MODE_READ_WRITE,
                   ndatabits::Integer=8,
                   parity::SPParity=SP_PARITY_NONE,
                   nstopbits::Integer=1)
    sp = SerialPort(sp_get_port_by_name(portname), false, true)
    sp_open(sp.ref, mode)
    set_speed(sp, bps)
    set_frame(sp, ndatabits=ndatabits, parity=parity, nstopbits=nstopbits)
    return sp
end

"""
`open_serial_port(port_address::AbstractString, baudrate::Integer)`

Create and configure a SerialPort object with the standard 8N1 settings
and specified `baudrate`. Example: `open_serial_port("/dev/ttyACM0", 115200)`
"""
function open_serial_port(port_address::AbstractString, baudrate::Integer)
    sp = SerialPort(port_address)
    open(sp)
    set_speed(sp, baudrate)
    set_frame(sp, ndatabits=8, parity=SP_PARITY_NONE, nstopbits=1)
    return sp
end

"""
close(sp::SerialPort)

Close the serial port `sp`.
"""
function Base.close(sp::SerialPort)
    if sp.open
        # Flush first, as is done in other close() methods in Base
        sp_flush(sp.ref, SP_BUF_BOTH)

        sp_close(sp.ref)
        sp.open = false
    end
    return sp
end

"""
`flush(sp::SerialPort [, buffer::SPBuffer])`

Flush `buffer` of serial port `sp`.

`buffer` can take the values: `SP_BUF_INPUT`, `SP_BUF_OUTPUT`, and
`SP_BUF_BOTH`.
"""
function Base.flush(sp::SerialPort; buffer::SPBuffer=SP_BUF_BOTH)
    sp_flush(sp.ref, buffer)
end

"""
`write(sp::SerialPort, data::String)`
`write(sp::SerialPort, data::Array{Char})`
`write(sp::SerialPort, data::Array{UInt8})`

Write sequence of Bytes to `sp`.
"""
function Base.write(sp::SerialPort, data::Array{UInt8})
    sp_nonblocking_write(sp.ref, data)
    return sp_drain(sp.ref)
end

function Base.write(sp::SerialPort, data::Array{Char})
    sp_nonblocking_write(sp.ref, data)
    return sp_drain(sp.ref)
end

function Base.write(sp::SerialPort, data::String)
    # sp_nonblocking_write(sp.ref, convert(Array{UInt8},data))
    sp_nonblocking_write(sp.ref, data)
    return sp_drain(sp.ref)
end

"""
`write(sp::SerialPort, data::UInt8)`
`write(sp::SerialPort, data::Char)`

Write single Byte to `sp`
"""
Base.write(sp::SerialPort, data::Char) = write(sp, [data])
Base.write(sp::SerialPort, data::UInt8) = write(sp, [data])

"""
`write(sp::SerialPort, i::Integer)`

Write string representation of `i` to `sp`.
"""
Base.write(sp::SerialPort, i::Integer) = Base.write(sp, "$i")

"""
`write(sp::SerialPort, f::AbstractFloat)`
`write(sp::SerialPort, f::AbstractFloat, format::AbstractString)`

Write formatted string representation of `f` to `sp`. By default the string is
formated using `format="%.3f"`. For details on the format consult the
documentation of the C library function `sprintf`.
"""
Base.write(sp::SerialPort, f::AbstractFloat, format::AbstractString="%.3f") = Base.write(sp, eval(:@sprintf($format, $f)))

"""
`eof(sp::SerialPort)`

Return EOF state (`true` or `false`) of `sp`.
"""
Base.eof(sp::SerialPort) = sp.eof

"""
`seteof(sp::SerialPort, state::Bool)`

Set EOF of `sp` to `state`
"""
function seteof(sp::SerialPort, state::Bool)
    sp.eof = state
    return nothing
end

"""
`seteof(sp::SerialPort, state::Bool)`

Reset EOF of `sp` to `false`
"""
reseteof(sp::SerialPort) = seteof(sp, false)

"""
`read(sp::SerialPort, T::Type{UInt8})`
`read(sp::SerialPort, T::Type{Char})`

Read a single Byte from the specified port and return it represented as `T`.
`T` might be either `Char or `UInt8`. if no Byte is availible in the port
buffer return zero.
"""
function Base.read(sp::SerialPort, readType::Type{Char})
    nbytes_read, bytes = sp_nonblocking_read(sp.ref, 1)
    return (nbytes_read == 1) ? convert(readType,bytes[1]) : readType(0)
end

function Base.read(sp::SerialPort, readType::Type{UInt8})
    nbytes_read, bytes = sp_nonblocking_read(sp.ref, 1)
    return (nbytes_read == 1) ? convert(readType,bytes[1]) : readType(0)
end

"""
`readuntil(sp::SerialPort,delim::Union{Char,AbstractString,Vector{Char}},timeout_ms::Integer)`

Read until the specified delimiting byte (e.g. `'\\n'`) is encountered, or until
timeout_ms has elapsed, whichever comes first.
"""
function Base.readuntil(sp::SerialPort, delim::Char, timeout_ms::Real)
    return readuntil(sp,[delim],timeout_ms)
end

function Base.readuntil(sp::SerialPort, delim::AbstractString, timeout_ms::Real)
    return readuntil(sp,convert(Vector{Char},delim),timeout_ms)
end

function Base.readuntil(sp::SerialPort, delim::Vector{Char}, timeout_ms::Real)
    start_time = time_ns()
    out = IOBuffer()
    lastchars = Char[0 for i=1:length(delim)]
    while !eof(sp)
        if (time_ns() - start_time)/1e6 > timeout_ms
            break
        end
        if bytesavailable(sp) > 0
            c = read(sp, Char)
            write(out, c)
            lastchars = circshift(lastchars,-1)
            lastchars[end] = c
            if lastchars == delim
                break
            end
        end
    end
    return String(take!(out))
end

"""
`Base.bytesavailable(sp::SerialPort)`

Gets the number of bytes waiting in the input buffer.
"""
Base.bytesavailable(sp::SerialPort) = Int(sp_input_waiting(sp.ref))

"""
`readbytes!(sp::SerialPort,nbytes::Integer)`

Read `nbytes` from the specified serial port `sp`, without blocking. Returns
a `UInt8` `Array`.
"""
function Base.readbytes!(sp::SerialPort, nbytes::Integer)
    nbytes_read, bytes = sp_nonblocking_read(sp.ref, nbytes)
    return bytes
end

"""
`read(sp::SerialPort, ::Type{String})`

Read everything from the specified serial ports `sp` input buffer, one byte at
a time, until it is empty. Returns a `String`.
"""
function Base.read(sp::SerialPort, ::Type{String})
    result = Char[]
    while Int(bytesavailable(sp)) > 0
        byte = readbytes!(sp, 1)[1]
        push!(result, byte)
    end
    return String(join(result))
end