text:
# Note: agreement = output * target
#       Agreement is high when output and target are the same sign and |output| is large.
#       It is an indication that the output represents the correct class in a margin-based model.

# ============================================================

@doc doc"""
    ZeroOneLoss <: MarginLoss

The classical classification loss. It penalizes every misclassified
observation with a loss of `1` while every correctly classified
observation has a loss of `0`.
It is not convex nor continuous and thus seldom used directly.
Instead one usually works with some classification-calibrated
surrogate loss, such as [L1HingeLoss](@ref).

```math
L(a) = \begin{cases} 1 & \quad \text{if } a < 0 \\ 0 & \quad \text{if } a >= 0\\ \end{cases}
```

---
```
              Lossfunction                     Derivative
      ┌────────────┬────────────┐      ┌────────────┬────────────┐
    1 │------------┐            │    1 │                         │
      │            |            │      │                         │
      │            |            │      │                         │
      │            |            │      │_________________________│
      │            |            │      │                         │
      │            |            │      │                         │
      │            |            │      │                         │
    0 │            └------------│   -1 │                         │
      └────────────┴────────────┘      └────────────┴────────────┘
      -2                        2      -2                        2
                y * h(x)                         y * h(x)
```
"""
struct ZeroOneLoss <: MarginLoss end

deriv(loss::ZeroOneLoss, target::Number, output::Number) = zero(output)
deriv2(loss::ZeroOneLoss, target::Number, output::Number) = zero(output)

value(loss::ZeroOneLoss, agreement::T) where {T<:Number} = sign(agreement) < 0 ? one(T) : zero(T)
deriv(loss::ZeroOneLoss, agreement::T) where {T<:Number} = zero(T)
deriv2(loss::ZeroOneLoss, agreement::T) where {T<:Number} = zero(T)
value_deriv(loss::ZeroOneLoss, agreement::T) where {T<:Number} = sign(agreement) < 0 ? (one(T), zero(T)) : (zero(T), zero(T))

isminimizable(::ZeroOneLoss) = true
isdifferentiable(::ZeroOneLoss) = false
isdifferentiable(::ZeroOneLoss, at) = at != 0
istwicedifferentiable(::ZeroOneLoss) = false
istwicedifferentiable(::ZeroOneLoss, at) = at != 0
isnemitski(::ZeroOneLoss) = true
islipschitzcont(::ZeroOneLoss) = true
isconvex(::ZeroOneLoss) = false
isclasscalibrated(loss::ZeroOneLoss) = true
isclipable(::ZeroOneLoss) = true

# ============================================================

@doc doc"""
    PerceptronLoss <: MarginLoss

The perceptron loss linearly penalizes every prediction where the
resulting `agreement <= 0`.
It is Lipschitz continuous and convex, but not strictly convex.

```math
L(a) = \max \{ 0, -a \}
```

---
```
              Lossfunction                     Derivative
      ┌────────────┬────────────┐      ┌────────────┬────────────┐
    2 │\.                       │    0 │            ┌------------│
      │ '..                     │      │            |            │
      │   \.                    │      │            |            │
      │     '.                  │      │            |            │
    L │      '.                 │   L' │            |            │
      │        \.               │      │            |            │
      │         '.              │      │            |            │
    0 │           \.____________│   -1 │------------┘            │
      └────────────┴────────────┘      └────────────┴────────────┘
      -2                        2      -2                        2
                 y ⋅ ŷ                            y ⋅ ŷ
```
"""
struct PerceptronLoss <: MarginLoss end

value(loss::PerceptronLoss, agreement::T) where {T<:Number} = max(zero(T), -agreement)
deriv(loss::PerceptronLoss, agreement::T) where {T<:Number} = agreement >= 0 ? zero(T) : -one(T)
deriv2(loss::PerceptronLoss, agreement::T) where {T<:Number} = zero(T)
value_deriv(loss::PerceptronLoss, agreement::T) where {T<:Number} = agreement >= 0 ? (zero(T), zero(T)) : (-agreement, -one(T))

isdifferentiable(::PerceptronLoss) = false
isdifferentiable(::PerceptronLoss, at) = at != 0
istwicedifferentiable(::PerceptronLoss) = false
istwicedifferentiable(::PerceptronLoss, at) = at != 0
islipschitzcont(::PerceptronLoss) = true
isconvex(::PerceptronLoss) = true
isstrictlyconvex(::PerceptronLoss) = false
isstronglyconvex(::PerceptronLoss) = false
isclipable(::PerceptronLoss) = true

# ============================================================

@doc doc"""
    LogitMarginLoss <: MarginLoss

The margin version of the logistic loss. It is infinitely many
times differentiable, strictly convex, and Lipschitz continuous.

```math
L(a) = \ln (1 + e^{-a})
```

---
```
              Lossfunction                     Derivative
      ┌────────────┬────────────┐      ┌────────────┬────────────┐
    2 │ \.                      │    0 │                  ._--/""│
      │   \.                    │      │               ../'      │
      │     \.                  │      │              ./         │
      │       \..               │      │            ./'          │
    L │         '-_             │   L' │          .,'            │
      │            '-_          │      │         ./              │
      │               '\-._     │      │      .,/'               │
    0 │                    '""*-│   -1 │__.--''                  │
      └────────────┴────────────┘      └────────────┴────────────┘
      -2                        2      -4                        4
                 y ⋅ ŷ                            y ⋅ ŷ
```
"""
struct LogitMarginLoss <: MarginLoss end
value(loss::LogitMarginLoss, agreement::Number) = log1p(exp(-agreement))
deriv(loss::LogitMarginLoss, agreement::Number) = -one(agreement) / (one(agreement) + exp(agreement))
deriv2(loss::LogitMarginLoss, agreement::Number) = (eᵗ = exp(agreement); eᵗ / abs2(one(eᵗ) + eᵗ))
value_deriv(loss::LogitMarginLoss, agreement::Number) = (eᵗ = exp(-agreement); (log1p(eᵗ), -eᵗ / (one(eᵗ) + eᵗ)))

isunivfishercons(::LogitMarginLoss) = true
isdifferentiable(::LogitMarginLoss) = true
isdifferentiable(::LogitMarginLoss, at) = true
istwicedifferentiable(::LogitMarginLoss) = true
istwicedifferentiable(::LogitMarginLoss, at) = true
islipschitzcont(::LogitMarginLoss) = true
isconvex(::LogitMarginLoss) = true
isstrictlyconvex(::LogitMarginLoss) = true
isstronglyconvex(::LogitMarginLoss) = false
isclipable(::LogitMarginLoss) = false

# ============================================================

@doc doc"""
    L1HingeLoss <: MarginLoss

The hinge loss linearly penalizes every predicition where the
resulting `agreement < 1` .
It is Lipschitz continuous and convex, but not strictly convex.

```math
L(a) = \max \{ 0, 1 - a \}
```

---
```
              Lossfunction                     Derivative
      ┌────────────┬────────────┐      ┌────────────┬────────────┐
    3 │'\.                      │    0 │                  ┌------│
      │  ''_                    │      │                  |      │
      │     \.                  │      │                  |      │
      │       '.                │      │                  |      │
    L │         ''_             │   L' │                  |      │
      │            \.           │      │                  |      │
      │              '.         │      │                  |      │
    0 │                ''_______│   -1 │------------------┘      │
      └────────────┴────────────┘      └────────────┴────────────┘
      -2                        2      -2                        2
                 y ⋅ ŷ                            y ⋅ ŷ
```
"""
struct L1HingeLoss <: MarginLoss end
const HingeLoss = L1HingeLoss

value(loss::L1HingeLoss, agreement::T) where {T<:Number} = max(zero(T), one(T) - agreement)
deriv(loss::L1HingeLoss, agreement::T) where {T<:Number} = agreement >= 1 ? zero(T) : -one(T)
deriv2(loss::L1HingeLoss, agreement::T) where {T<:Number} = zero(T)
value_deriv(loss::L1HingeLoss, agreement::T) where {T<:Number} = agreement >= 1 ? (zero(T), zero(T)) : (one(T) - agreement, -one(T))

isfishercons(::L1HingeLoss) = true
isdifferentiable(::L1HingeLoss) = false
isdifferentiable(::L1HingeLoss, at) = at != 1
istwicedifferentiable(::L1HingeLoss) = false
istwicedifferentiable(::L1HingeLoss, at) = at != 1
islipschitzcont(::L1HingeLoss) = true
isconvex(::L1HingeLoss) = true
isstrictlyconvex(::L1HingeLoss) = false
isstronglyconvex(::L1HingeLoss) = false
isclipable(::L1HingeLoss) = true

# ============================================================

@doc doc"""
    L2HingeLoss <: MarginLoss

The truncated least squares loss quadratically penalizes every
predicition where the resulting `agreement < 1`.
It is locally Lipschitz continuous and convex,
but not strictly convex.

```math
L(a) = \max \{ 0, 1 - a \}^2
```

---
```
              Lossfunction                     Derivative
      ┌────────────┬────────────┐      ┌────────────┬────────────┐
    5 │     .                   │    0 │                 ,r------│
      │     '.                  │      │               ,/        │
      │      '\                 │      │             ,/          │
      │        \                │      │           ,/            │
    L │         '.              │   L' │         ./              │
      │          '.             │      │       ./                │
      │            \.           │      │     ./                  │
    0 │              '-.________│   -5 │   ./                    │
      └────────────┴────────────┘      └────────────┴────────────┘
      -2                        2      -2                        2
                 y ⋅ ŷ                            y ⋅ ŷ
```
"""
struct L2HingeLoss <: MarginLoss end

value(loss::L2HingeLoss, agreement::T) where {T<:Number} = agreement >= 1 ? zero(T) : abs2(one(T) - agreement)
deriv(loss::L2HingeLoss, agreement::T) where {T<:Number} = agreement >= 1 ? zero(T) : convert(T,2) * (agreement - one(T))
deriv2(loss::L2HingeLoss, agreement::T) where {T<:Number} = agreement >= 1 ? zero(T) : convert(T,2)
value_deriv(loss::L2HingeLoss, agreement::T) where {T<:Number} = agreement >= 1 ? (zero(T), zero(T)) : (abs2(one(T) - agreement), convert(T,2) * (agreement - one(T)))

isunivfishercons(::L2HingeLoss) = true
isdifferentiable(::L2HingeLoss) = true
isdifferentiable(::L2HingeLoss, at) = true
istwicedifferentiable(::L2HingeLoss) = false
istwicedifferentiable(::L2HingeLoss, at) = at != 1
islocallylipschitzcont(::L2HingeLoss) = true
islipschitzcont(::L2HingeLoss) = false
isconvex(::L2HingeLoss) = true
isstrictlyconvex(::L2HingeLoss) = false
isstronglyconvex(::L2HingeLoss) = false
isclipable(::L2HingeLoss) = true

# ============================================================

@doc doc"""
    SmoothedL1HingeLoss <: MarginLoss

As the name suggests a smoothed version of the L1 hinge loss.
It is Lipschitz continuous and convex, but not strictly convex.

```math
L(a) = \begin{cases} \frac{0.5}{\gamma} \cdot \max \{ 0, 1 - a \} ^2 & \quad \text{if } a \ge 1 - \gamma \\ 1 - \frac{\gamma}{2} - a & \quad \text{otherwise}\\ \end{cases}
```

---
```
              Lossfunction (γ=2)               Derivative
      ┌────────────┬────────────┐      ┌────────────┬────────────┐
    2 │\.                       │    0 │                 ,r------│
      │ '.                      │      │               ./'       │
      │   \.                    │      │              ,/         │
      │     '.                  │      │            ./'          │
    L │      '.                 │   L' │           ,'            │
      │        \.               │      │         ,/              │
      │          ',             │      │       ./'               │
    0 │            '*-._________│   -1 │______./                 │
      └────────────┴────────────┘      └────────────┴────────────┘
      -2                        2      -2                        2
                 y ⋅ ŷ                            y ⋅ ŷ
```
"""
struct SmoothedL1HingeLoss{T<:AbstractFloat} <: MarginLoss
    gamma::T

    function SmoothedL1HingeLoss{T}(γ::T) where T
        γ > 0 || error("γ must be strictly positive")
        new{T}(γ)
    end
end
SmoothedL1HingeLoss(γ::T) where {T<:AbstractFloat} = SmoothedL1HingeLoss{T}(γ)
SmoothedL1HingeLoss(γ) = SmoothedL1HingeLoss(Float64(γ))

function value(loss::SmoothedL1HingeLoss{R}, agreement::T)::promote_type(R,T) where {R,T<:Number}
    if agreement >= 1 - loss.gamma
        R(0.5) / loss.gamma * abs2(max(zero(T), one(T) - agreement))
    else
        one(T) - loss.gamma / R(2) - agreement
    end
end
function deriv(loss::SmoothedL1HingeLoss{R}, agreement::T)::promote_type(R,T) where {R,T<:Number}
    if agreement >= 1 - loss.gamma
        agreement >= 1 ? zero(T) : (agreement - one(T)) / loss.gamma
    else
        -one(T)
    end
end
function deriv2(loss::SmoothedL1HingeLoss{R}, agreement::T)::promote_type(R,T) where {R,T<:Number}
    agreement < 1 - loss.gamma || agreement > 1 ? zero(T) : one(T) / loss.gamma
end

isdifferentiable(::SmoothedL1HingeLoss) = true
isdifferentiable(::SmoothedL1HingeLoss, at) = true
istwicedifferentiable(::SmoothedL1HingeLoss) = false
istwicedifferentiable(loss::SmoothedL1HingeLoss, at) = at != 1 && at != 1 - loss.gamma
islocallylipschitzcont(::SmoothedL1HingeLoss) = true
islipschitzcont(::SmoothedL1HingeLoss) = true
isconvex(::SmoothedL1HingeLoss) = true
isstrictlyconvex(::SmoothedL1HingeLoss) = false
isstronglyconvex(::SmoothedL1HingeLoss) = false
isclipable(::SmoothedL1HingeLoss) = true

# ============================================================

@doc doc"""
    ModifiedHuberLoss <: MarginLoss

A special (4 times scaled) case of the [`SmoothedL1HingeLoss`](@ref)
with `γ=2`. It is Lipschitz continuous and convex,
but not strictly convex.

```math
L(a) = \begin{cases} \max \{ 0, 1 - a \} ^2 & \quad \text{if } a \ge -1 \\ - 4 a & \quad \text{otherwise}\\ \end{cases}
```

---
```
              Lossfunction                     Derivative
      ┌────────────┬────────────┐      ┌────────────┬────────────┐
    5 │    '.                   │    0 │                .+-------│
      │     '.                  │      │              ./'        │
      │      '\                 │      │             ,/          │
      │        \                │      │           ,/            │
    L │         '.              │   L' │         ./              │
      │          '.             │      │       ./'               │
      │            \.           │      │______/'                 │
    0 │              '-.________│   -5 │                         │
      └────────────┴────────────┘      └────────────┴────────────┘
      -2                        2      -2                        2
                 y ⋅ ŷ                            y ⋅ ŷ
```
"""
struct ModifiedHuberLoss <: MarginLoss end

function value(loss::ModifiedHuberLoss, agreement::T) where T<:Number
    agreement >= -1 ? abs2(max(zero(T), one(agreement) - agreement)) : -convert(T,4) * agreement
end
function deriv(loss::ModifiedHuberLoss, agreement::T) where T<:Number
    if agreement >= -1
        agreement > 1 ? zero(T) : convert(T,2)*agreement - convert(T,2)
    else
        -convert(T,4)
    end
end
function deriv2(loss::ModifiedHuberLoss, agreement::T) where T<:Number
    agreement < -1 || agreement > 1 ? zero(T) : convert(T,2)
end

isdifferentiable(::ModifiedHuberLoss) = true
isdifferentiable(::ModifiedHuberLoss, at) = true
istwicedifferentiable(::ModifiedHuberLoss) = false
istwicedifferentiable(loss::ModifiedHuberLoss, at) = at != 1 && at != -1
islocallylipschitzcont(::ModifiedHuberLoss) = true
islipschitzcont(::ModifiedHuberLoss) = true
isconvex(::ModifiedHuberLoss) = true
isstrictlyconvex(::ModifiedHuberLoss) = false
isstronglyconvex(::ModifiedHuberLoss) = false
isclipable(::ModifiedHuberLoss) = true

# ============================================================

@doc doc"""
    L2MarginLoss <: MarginLoss

The margin-based least-squares loss for classification,
which penalizes every prediction where `agreement != 1` quadratically.
It is locally Lipschitz continuous and strongly convex.

```math
L(a) = {\left( 1 - a \right)}^2
```

---
```
              Lossfunction                     Derivative
      ┌────────────┬────────────┐      ┌────────────┬────────────┐
    5 │     .                   │    2 │                       ,r│
      │     '.                  │      │                     ,/  │
      │      '\                 │      │                   ,/    │
      │        \                │      ├                 ,/      ┤
    L │         '.              │   L' │               ./        │
      │          '.             │      │             ./          │
      │            \.          .│      │           ./            │
    0 │              '-.____.-' │   -3 │         ./              │
      └────────────┴────────────┘      └────────────┴────────────┘
      -2                        2      -2                        2
                 y ⋅ ŷ                            y ⋅ ŷ
```
"""
struct L2MarginLoss <: MarginLoss end

value(loss::L2MarginLoss, agreement::T) where {T<:Number} = abs2(one(T) - agreement)
deriv(loss::L2MarginLoss, agreement::T) where {T<:Number} = convert(T,2) * (agreement - one(T))
deriv2(loss::L2MarginLoss, agreement::T) where {T<:Number} = convert(T,2)

isunivfishercons(::L2MarginLoss) = true
isdifferentiable(::L2MarginLoss) = true
isdifferentiable(::L2MarginLoss, at) = true
istwicedifferentiable(::L2MarginLoss) = true
istwicedifferentiable(::L2MarginLoss, at) = true
islocallylipschitzcont(::L2MarginLoss) = true
islipschitzcont(::L2MarginLoss) = false
isconvex(::L2MarginLoss) = true
isstrictlyconvex(::L2MarginLoss) = true
isstronglyconvex(::L2MarginLoss) = true
isclipable(::L2MarginLoss) = true

# ============================================================

@doc doc"""
    ExpLoss <: MarginLoss

The margin-based exponential loss for classification, which
penalizes every prediction exponentially. It is infinitely many
times differentiable, locally Lipschitz continuous and strictly
convex, but not clipable.

```math
L(a) = e^{-a}
```

---
```
              Lossfunction                     Derivative
      ┌────────────┬────────────┐      ┌────────────┬────────────┐
    5 │  \.                     │    0 │               _,,---:'""│
      │   l                     │      │           _r/"'         │
      │    l.                   │      │        .r/'             │
      │     ":                  │      │      .r'                │
    L │       \.                │   L' │     ./                  │
      │        "\..             │      │    .'                   │
      │           '":,_         │      │   ,'                    │
    0 │                ""---:.__│   -5 │  ./                     │
      └────────────┴────────────┘      └────────────┴────────────┘
      -2                        2      -2                        2
                 y ⋅ ŷ                            y ⋅ ŷ
```
"""
struct ExpLoss <: MarginLoss end

value(loss::ExpLoss, agreement::Number) = exp(-agreement)
deriv(loss::ExpLoss, agreement::Number) = -exp(-agreement)
deriv2(loss::ExpLoss, agreement::Number) = exp(-agreement)
value_deriv(loss::ExpLoss, agreement::Number) = (eᵗ = exp(-agreement); (eᵗ, -eᵗ))

isunivfishercons(::ExpLoss) = true
isdifferentiable(::ExpLoss) = true
isdifferentiable(::ExpLoss, at) = true
istwicedifferentiable(::ExpLoss) = true
istwicedifferentiable(::ExpLoss, at) = true
islocallylipschitzcont(::ExpLoss) = true
islipschitzcont(::ExpLoss) = false
isconvex(::ExpLoss) = true
isstrictlyconvex(::ExpLoss) = true
isstronglyconvex(::ExpLoss) = false
isclipable(::ExpLoss) = false

# ============================================================

@doc doc"""
    SigmoidLoss <: MarginLoss

Continuous loss which penalizes every prediction with a loss
within in the range (0,2). It is infinitely many times
differentiable, Lipschitz continuous but nonconvex.

```math
L(a) = 1 - \tanh(a)
```

---
```
              Lossfunction                     Derivative
      ┌────────────┬────────────┐      ┌────────────┬────────────┐
    2 │""'--,.                  │    0 │..                     ..│
      │      '\.                │      │ "\.                 ./" │
      │         '.              │      │    ',             ,'    │
      │           \.            │      │      \           /      │
    L │            "\.          │   L' │       \         /       │
      │              \.         │      │        \.     ./        │
      │                \,       │      │         \.   ./         │
    0 │                  '"-:.__│   -1 │          ',_,'          │
      └────────────┴────────────┘      └────────────┴────────────┘
      -2                        2      -2                        2
                 y ⋅ ŷ                            y ⋅ ŷ
```
"""
struct SigmoidLoss <: MarginLoss end

value(loss::SigmoidLoss, agreement::Number) = one(agreement) - tanh(agreement)
deriv(loss::SigmoidLoss, agreement::Number) = -abs2(sech(agreement))
deriv2(loss::SigmoidLoss, agreement::T) where {T<:Number} = convert(T,2) * tanh(agreement) * abs2(sech(agreement))

isunivfishercons(::SigmoidLoss) = true
isdifferentiable(::SigmoidLoss) = true
isdifferentiable(::SigmoidLoss, at) = true
istwicedifferentiable(::SigmoidLoss) = true
istwicedifferentiable(::SigmoidLoss, at) = true
islocallylipschitzcont(::SigmoidLoss) = true
islipschitzcont(::SigmoidLoss) = true
isclasscalibrated(::SigmoidLoss) = true
isconvex(::SigmoidLoss) = false
isstrictlyconvex(::SigmoidLoss) = false
isstronglyconvex(::SigmoidLoss) = false
isclipable(::SigmoidLoss) = false

# ============================================================

@doc doc"""
    DWDMarginLoss <: MarginLoss

The distance weighted discrimination margin loss. It is a
differentiable generalization of the [L1HingeLoss](@ref) that is
different than the [SmoothedL1HingeLoss](@ref). It is Lipschitz
continuous and convex, but not strictly convex.

```math
L(a) = \begin{cases} 1 - a & \quad \text{if } a \ge \frac{q}{q+1} \\ \frac{1}{a^q} \frac{q^q}{(q+1)^{q+1}} & \quad \text{otherwise}\\ \end{cases}
```

---
```
              Lossfunction (q=1)               Derivative
      ┌────────────┬────────────┐      ┌────────────┬────────────┐
    2 │      ".                 │    0 │                     ._r-│
      │        \.               │      │                   ./    │
      │         ',              │      │                 ./      │
      │           \.            │      │                 /       │
    L │            "\.          │   L' │                .        │
      │              \.         │      │                /        │
      │               ":__      │      │               ;         │
    0 │                   '""---│   -1 │---------------┘         │
      └────────────┴────────────┘      └────────────┴────────────┘
      -2                        2      -2                        2
                 y ⋅ ŷ                            y ⋅ ŷ
```
"""
struct DWDMarginLoss{T<:AbstractFloat} <: MarginLoss
    q::T
    function DWDMarginLoss{T}(q::T) where T
        q > 0 || error("q must be strictly positive")
        new{T}(q)
    end
end
DWDMarginLoss(q::T) where {T<:AbstractFloat} = DWDMarginLoss{T}(q)
DWDMarginLoss(q) = DWDMarginLoss(Float64(q))

function value(loss::DWDMarginLoss{R}, agreement::T)::promote_type(R, T) where {R,T<:Number}
    q = loss.q
    if agreement <= q/(q+1)
        R(1) - agreement
    else
        (q^q/(q+1)^(q+1)) / agreement^q
    end
end

function deriv(loss::DWDMarginLoss{R}, agreement::T)::promote_type(R, T) where {R,T<:Number}
    q = loss.q
    agreement <= q/(q+1) ? -one(T) : -(q/(q+1))^(q+1) / agreement^(q+1)
end

function deriv2(loss::DWDMarginLoss{R}, agreement::T)::promote_type(R, T) where {R,T<:Number}
    q = loss.q
    agreement <= q/(q+1) ? zero(T) : ( (q^(q+1))/((q+1)^q) ) / agreement^(q+2)
end

isdifferentiable(::DWDMarginLoss) = true
isdifferentiable(::DWDMarginLoss, at) = true
istwicedifferentiable(::DWDMarginLoss) = true
istwicedifferentiable(loss::DWDMarginLoss, at) = true
islocallylipschitzcont(::DWDMarginLoss) = true
islipschitzcont(::DWDMarginLoss) = true
isconvex(::DWDMarginLoss) = true
isstrictlyconvex(::DWDMarginLoss) = false
isstronglyconvex(::DWDMarginLoss) = false
isfishercons(::DWDMarginLoss) = true
isunivfishercons(::DWDMarginLoss) = true
isclipable(::DWDMarginLoss) = false


text:
using Polynomials
"""
    _getpolylagrange(k::Int64, order::Int64, origin::Int64)

Function that return the k-th Lagrange Polynomial of a certain order. Coefficients are rational then the return is exact. The polynomial is equal to :
``\\prod_{i=0,\\ i \\neq k}^{order} \\frac{x - i - origin}{k - i}``

# Arguments
- `k::Int64` : number of the Polynomial, `k` must be between `0` and `order` (`0<= k <= order`).
- `order::Int64` : order of the polynomial.
- `origin::Int64` : origin of the first indice.

# Returns
- `Polynomial{Rational{BigInt}}` : the k-th Lagrange polynomial of order `order`

# Throws
- `DommaineError` : when `0 <= k <= order` is `false` or when N ∉ {BInt64, BigInt}
"""
function _getpolylagrange(k::Int, order::Int, origin::Int) where {N<:Integer}
    0 <= k <= order || throw(DomainError("the constant 0 <= k <= order is false"))
    # the computed is made with big rational
    result = Polynomials.Polynomial([big(1 // 1)])
    for l = 0:order
        if l != k
            result *= Polynomials.Polynomial([-(l + origin), 1 // 1] .// (k - l))
        end
    end
    return result
end

"""
    Lagrange{T, edge, order, N} <: AbstractInterpolation{T, edge, order}
    Lagrange(order, T::DataType=Float64; edge::EdgeType=CircEdge)

Type containing Lagrange Polynomials coefficients for Lagrange interpolation

# Type parameters
- `T` : the type of data that is interpolate
- `edge::EdgeType` : type of edge traitment
- `order::Int`: order of lagrange interpolation

# Implementation :
- `tabfct::Vector{Polynomial{T}}` : vector of all lagrange polynomial, per example the k-th Lagrange polynomial for the designed order is tabfct[k+1]

# Arguments :
- `order::Int` : the order of interpolation
- `[T::DataType=Float64]` : The type values to interpolate

# Keywords arguments :
- `edge::EdgeType=CircEdge` : type of edge traitment

"""
struct Lagrange{T,edge,order} <: AbstractInterpolation{T,edge,order}
    tabfct::Vector{Polynomial{T}}
    function Lagrange(order::Int, T::DataType = Float64; edge::EdgeType = CircEdge)
        origin = -div(order, 2)
        tabfct_rat = collect([_getpolylagrange(i, order, origin) for i = 0:order])
        new{T,edge,order}(convert.(Polynomial{T}, tabfct_rat))
    end
end

function _c(k,n)
    p = integrate(_getpolylagrange(k,n,0))
    return p(0)-p(-1)
end
struct ABcoef
    tab::Array{Rational{Int}}
    function ABcoef(ordermax::Int)
        tab = zeros(Rational{Int}, ordermax, ordermax)
        for j=1:ordermax, i=1:j
            tab[i,j] = _c(i-1,j-1)
        end
        return new(tab)
    end
end
c(st::ABcoef, k, n)=st.tab[k,n]

ext: jl


text:
########################################################################################
# Solution
########################################################################################
#! format: off
abstract type AbstractProbODESolution{T,N,S} <: DiffEqBase.AbstractODESolution{T,N,S} end
mutable struct ProbODESolution{
    T, N, uType, puType, uType2, DType, tType, rateType, xType, diffType, llType, P, A, IType, DE
} <: AbstractProbODESolution{T,N,uType}
    u::uType
    pu::puType
    u_analytic::uType2
    errors::DType
    t::tType
    k::rateType
    x_filt::xType
    x_smooth::xType
    diffusions::diffType
    log_likelihood::llType
    prob::P
    alg::A
    interp::IType
    dense::Bool
    tslocation::Int
    destats::DE
    retcode::Symbol
end
ProbODESolution{T,N}(
    u, pu, u_analytic, errors, t, k, x_filt, x_smooth, diffusions, log_likelihood, prob,
    alg, interp, dense, tslocation, destats, retcode,
) where {T,N} = ProbODESolution{
    T, N, typeof(u), typeof(pu), typeof(u_analytic), typeof(errors), typeof(t), typeof(k),
    typeof(x_filt), typeof(diffusions), typeof(log_likelihood), typeof(prob), typeof(alg),
    typeof(interp), typeof(destats),
}(
    u, pu, u_analytic, errors, t, k, x_filt, x_smooth, diffusions, log_likelihood, prob,
    alg, interp, dense, tslocation, destats, retcode,
)

function DiffEqBase.solution_new_retcode(sol::ProbODESolution{T,N}, retcode) where {T,N}
    return ProbODESolution{T,N}(
        sol.u, sol.pu, sol.u_analytic, sol.errors, sol.t, sol.k, sol.x_filt, sol.x_smooth,
        sol.diffusions, sol.log_likelihood, sol.prob, sol.alg, sol.interp, sol.dense,
        sol.tslocation, sol.destats, retcode,
    )
end

# Used to build the initial empty solution in OrdinaryDiffEq.__init
function DiffEqBase.build_solution(
    prob::DiffEqBase.AbstractODEProblem,
    alg::GaussianODEFilter,
    t,
    u;
    k=nothing,
    retcode=:Default,
    destats=nothing,
    dense=true,
    kwargs...,
)
    T = eltype(eltype(u))
    N = length((size(prob.u0)..., length(u)))

    d = length(prob.u0)
    uElType = eltype(prob.u0)
    D = d
    pu_cov =
        alg isa EK0 && !(prob.f isa DynamicalODEFunction) ?
        SRMatrix(zeros(uElType, d, D), Diagonal(zeros(uElType, d, d))) :
        SRMatrix(zeros(uElType, d, D))
    x_cov = SRMatrix(zeros(uElType, d, d))
    pu = StructArray{Gaussian{Vector{uElType},typeof(pu_cov)}}(undef, 0)
    x_filt = StructArray{Gaussian{Vector{uElType},typeof(x_cov)}}(undef, 0)
    x_smooth = copy(x_filt)

    interp = GaussianODEFilterPosterior(alg, prob.u0)

    if DiffEqBase.has_analytic(prob.f)
        u_analytic = Vector{typeof(prob.u0)}()
        errors = Dict{Symbol,real(uElType)}()
    else
        u_analytic = nothing
        errors = nothing
    end

    ll = zero(uElType)
    return ProbODESolution{T,N}(
        u, pu, u_analytic, errors, t, [], x_filt, x_smooth, [], ll, prob, alg, interp,
        dense, 0, destats, retcode,
    )
end

function DiffEqBase.build_solution(
    sol::ProbODESolution{T,N},
    u_analytic,
    errors,
) where {T,N}
    return ProbODESolution{T,N}(
        sol.u, sol.pu, u_analytic, errors, sol.t, sol.k, sol.x_filt, sol.x_smooth,
        sol.diffusions, sol.log_likelihood, sol.prob, sol.alg, sol.interp, sol.dense,
        sol.tslocation, sol.destats, sol.retcode,
    )
end

########################################################################################
# Compat with classic ODE solutions, to enable analysis with DiffEqDevTools.jl
########################################################################################
"""Mean of a probabilistic ODE solution

Used for compatibility with classic ODE solutions, e.g. in `DiffEqDevTools.appxtrue`."""
mutable struct MeanProbODESolution{
    T, N, uType, uType2, DType, tType, rateType, P, A, IType, DE, PSolType
} <: DiffEqBase.AbstractODESolution{T,N,uType}
    u::uType
    u_analytic::uType2
    errors::DType
    t::tType
    k::rateType
    prob::P
    alg::A
    interp::IType
    dense::Bool
    tslocation::Int
    destats::DE
    retcode::Symbol
    probsol::PSolType
end
function mean(sol::ProbODESolution{T,N}) where {T,N}
    return MeanProbODESolution{
        T, N, typeof(sol.u), typeof(sol.u_analytic), typeof(sol.errors), typeof(sol.t),
        typeof(sol.k), typeof(sol.prob), typeof(sol.alg), typeof(sol.interp),
        typeof(sol.destats), typeof(sol),
    }(
        sol.u, sol.u_analytic, sol.errors, sol.t, sol.k, sol.prob, sol.alg, sol.interp,
        sol.dense, sol.tslocation, sol.destats, sol.retcode, sol,
    )
end
#! format: on
(sol::MeanProbODESolution)(t::Real, deriv::Val{N}=Val(0)) where {N} =
    mean(sol.probsol(t, deriv))
(sol::MeanProbODESolution)(t::AbstractVector, deriv=Val(0)) =
    DiffEqArray(mean(sol.probsol(t, deriv).u), t)
DiffEqBase.calculate_solution_errors!(sol::ProbODESolution, args...; kwargs...) =
    DiffEqBase.calculate_solution_errors!(mean(sol), args...; kwargs...)

########################################################################################
# Dense Output
########################################################################################
abstract type AbstractODEFilterPosterior <: DiffEqBase.AbstractDiffEqInterpolation end
struct GaussianODEFilterPosterior{SPType,AType,QType,PType} <: AbstractODEFilterPosterior
    d::Int
    q::Int
    SolProj::SPType
    A::AType
    Q::QType
    P::PType
    PI::PType
    smooth::Bool
end
set_smooth(p::GaussianODEFilterPosterior) =
    GaussianODEFilterPosterior(p.d, p.q, p.SolProj, p.A, p.Q, p.P, p.PI, true)
function GaussianODEFilterPosterior(alg, u0)
    uElType = eltype(u0)
    d = u0 isa ArrayPartition ? length(u0) ÷ 2 : length(u0)
    q = alg.order
    D = d * (q + 1)

    Proj = projection(d, q, uElType)
    SolProj = u0 isa ArrayPartition ? [Proj(1); Proj(0)] : Proj(0)

    A, Q = ibm(d, q, uElType)
    P, PI = init_preconditioner(d, q, uElType)
    return GaussianODEFilterPosterior(d, q, SolProj, A, Q, P, PI, false)
end
DiffEqBase.interp_summary(interp::GaussianODEFilterPosterior) =
    "Gaussian ODE Filter Posterior"

function (posterior::GaussianODEFilterPosterior)(
    tval::Real,
    t,
    x_filt,
    x_smooth,
    diffusions;
    smoothed=posterior.smooth,
)
    @unpack A, Q, d, q = posterior

    if tval < t[1]
        error("Invalid t<t0")
    end
    if tval in t
        idx = sum(t .<= tval)
        @assert t[idx] == tval
        return smoothed ? x_smooth[idx] : x_filt[idx]
    end

    idx = sum(t .<= tval)
    prev_t = t[idx]
    prev_rv = x_filt[idx]
    diffusion = diffusions[minimum((idx, end))]

    # Extrapolate
    h1 = tval - prev_t
    make_preconditioners!(posterior, h1)
    P, PI = posterior.P, posterior.PI
    Qh = apply_diffusion(Q, diffusion)
    goal_pred = predict(P * prev_rv, A, Qh)
    goal_pred = PI * goal_pred

    if !smoothed || tval >= t[end]
        return goal_pred
    end

    @assert length(x_filt) == length(x_smooth)
    next_t = t[idx+1]
    next_smoothed = x_smooth[idx+1]

    # Smooth
    h2 = next_t - tval
    make_preconditioners!(posterior, h2)
    P, PI = posterior.P, posterior.PI
    goal_pred = P * goal_pred
    next_smoothed = P * next_smoothed
    Qh = apply_diffusion(Q, diffusion)

    goal_smoothed, _ = smooth(goal_pred, next_smoothed, A, Qh)

    return PI * goal_smoothed
end
function (sol::ProbODESolution)(t::Real, deriv::Val{N}=Val(0)) where {N}
    @unpack q, d = sol.interp
    return sol.interp.SolProj *
           sol.interp(t, sol.t, sol.x_filt, sol.x_smooth, sol.diffusions)
end
(sol::ProbODESolution)(t::AbstractVector, deriv=Val(0)) =
    DiffEqArray(StructArray(sol.(t, deriv)), t)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: jl:percent
#     text_representation:
#       extension: .jl
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# %%
# Ensure your working directory is: https://github.com/henriquebecker91/phd/tree/master/latex/revised_PPG2KP
include("notebook_setup.jl")

# %%
# Read the data, show nothing.
raw_cdf_path = "./data/lagos2021.csv"
raw_cdf = DataFrame(CSV.File(raw_cdf_path))
#nothing
showtable(raw_cdf)

# %%
# Read the data, show nothing.
raw_idf_path = "./data/dataset_c_info.csv"
raw_idf = DataFrame(CSV.File(raw_idf_path))
showtable(raw_idf)

# %%
# Clean the data a little
idf = let raw_idf = deepcopy(raw_idf)
    raw_idf[!, :instance] = basename.(raw_idf[!, :instance])
    select!(raw_idf, :instance => :instance_name, Not(:instance))
    raw_ne2df_path = "./data/dataset_c_2NE_opt.csv"
    raw_ne2df = DataFrame(CSV.File(raw_ne2df_path))
    raw_idf = innerjoin(raw_idf, raw_ne2df; on = :instance_name)

    raw_idf
end
showtable(idf)

# %%
# Clean the data: replace sentinel values by missing; add a string id combining the parameter combinations;
# create some boolean columns from the data.
cdf = let cdf = deepcopy(raw_cdf)
    # Keep only the instance name (not the path).
    @with(cdf, :instance_name .= basename.(:instance_name))
    @with(cdf, :datafile .= basename.(:datafile))

    cdf[!, :finished] = .!cdf[!, :run_ended_by_exception]

    colnames = names(cdf)
    for name in colnames
        column = cdf[!, name]
        if eltype(column) <: AbstractFloat
            cdf[!, name] = ifelse.(isnan.(column), missing, column)
        elseif eltype(column) <: Integer
            cdf[!, name] = ifelse.(column .< 0, missing, column)
        end
    end

    # Create columns for those attributes that indicate their final state before the final solving.
    for attr in (:cmvars, :pevars, :plates)
        @assert all(cdf[!, "qt_$(attr)_after_preprocess"] .== cdf[!, "qt_$(attr)_after_purge"])
        rename!(cdf, "qt_$(attr)_after_preprocess" => attr)
        select!(cdf, Not("qt_$(attr)_after_purge"))
    end

    # this is brittle, but it works because every run finished with a solution (even if not optimal)
    cdf[!, :g2csp] = (!ismissing).(cdf[!, :solution_value])
    cdf[!, :g2kp] = (!ismissing).(cdf[!, :solution_profit])
    @assert all(cdf[!, :g2csp] .⊻ cdf[!, :g2kp])

    cdf[!, :opt] = @with(cdf, ifelse.(
        :finished, ifelse.(:g2kp, :solution_profit, :solution_value), missing
    ))

    cdf
end
showtable(cdf)
#nothing

# %%
# Remove all gcut instances for now
filter!(:instance_name => s -> isnothing(match(r"^gcut", s)), cdf)

# %%
# Remove enhanced without PSN and original with PSN.
filter!([:faithful, :round2disc] => (f, psn) -> (f && !psn) || (!f  && psn), cdf)

# %%
# Work only with the multi-seed runs there.
let cdf = deepcopy(cdf)
    # Get only the runs with seeds greater than one, and then add the respectives seed one runs.
    extra = filter(:seed => !isequal(1), cdf)
    pk_no_seed = [:instance_name, :round2disc, :faithful, :threads, :g2csp]
    extra = semijoin(cdf, extra; on = pk_no_seed)
    pk = push!(copy(pk_no_seed), :seed)
    @assert all(extra[!, :threads] .== 12)
    @assert all(extra[!, :g2csp])
    summary = sort!(extra, pk) |> df -> select!(df, pk..., :build_and_solve_time) |>
        df -> groupby(df, [:instance_name, :faithful]) |>
        df -> @based_on(df, min_time = minimum(:build_and_solve_time), max_time = maximum(:build_and_solve_time))
    summary[!, :abs_gap] = @with(summary, :max_time .- :min_time)
end

# %%
# Remove all seeds different from one
filter!(:seed => isequal(1), cdf)

# %%
cidf = let cdf = deepcopy(cdf), idf = deepcopy(idf)
    function something_if_possible(itr)
        nmiss = skipmissing(itr)
        return isempty(nmiss) ? missing : first(nmiss)
    end

    select!(idf, :instance_name, :N, :n_, :L, :W, :min_l => :l_, :min_w => :w_,
        :kp_ub => :kp_tub, :csp_lb => :csp_tub, :opt2ne => :csp_opt2ne)
    csp_cdf = @linq filter(r -> r.g2csp, cdf) |> groupby(:instance_name) |>
        based_on(
            csp_opt = something_if_possible(:opt),
            csp_rel = mean(skipmissing(:final_root_relaxation_value))
        )
    kp_cdf = @linq filter(r -> r.g2kp, cdf) |> groupby(:instance_name) |>
        based_on(
            kp_opt = something_if_possible(:opt),
            kp_rel = mean(skipmissing(:final_root_relaxation_value))
        )

    cidf = innerjoin(idf, csp_cdf, kp_cdf; on = :instance_name)
    select!(cidf,
        :instance_name, :N, :n_, :L, :W, :l_, :w_,
        :kp_tub, :kp_rel, :kp_opt, :csp_tub, :csp_rel, :csp_opt, :csp_opt2ne
    )
    sort!(cidf, :instance_name)
end
showtable(cidf)

# %%
cidf[!, :kp_gap_tub_rel] = @with(cidf, :kp_tub ./ :kp_rel) .* 100 .- 100
sort(cidf, :kp_gap_tub_rel; rev = true) |> df -> select(df, :instance_name, :kp_gap_tub_rel) |> showtable
#mean(cidf[!, :kp_gap_tub_rel])
cidf[!, :kp_gap_rel_opt] = @with(cidf, :kp_rel ./ :kp_opt) .* 100 .- 100
#@show mean(cidf[!, :kp_gap_rel_opt])
sort(cidf, :kp_gap_rel_opt; rev = true) |> df -> select(df, :instance_name, :kp_gap_rel_opt) |> showtable

cidf[!, :csp_gap_tub_rel] = @with(cidf, ceil.(:csp_rel) .- ceil.(:csp_tub))
@show sum(iszero.(cidf[!, :csp_gap_tub_rel]))
sort(cidf, :csp_gap_tub_rel) |> df -> select(df, :instance_name, :csp_gap_tub_rel) |> showtable

#=filter(r -> r.instance_name in ("W", "A1", "3"), cidf) |>
    df -> select(df, :instance_name, :csp_gap_tub_rel, :csp_opt) |>
    df -> (df[!, :per_csp_gap_tub_rel] = df[!, :csp_gap_tub_rel] ./ df[!, :csp_opt]; df) |>
    showtable
=#
cidf[!, :csp_gap_rel_opt] = @with(cidf, :csp_opt .- ceil.(:csp_rel))
sort(cidf, :csp_gap_rel_opt) |> df -> select(df, :instance_name, :csp_gap_rel_opt) |> showtable

#cidf[!, :csp_gap_tub_opt] = @with(cidf, :csp_opt .- ceil.(:csp_tub))
#@show sum(iszero.(skipmissing(cidf[!, :csp_gap_tub_opt])))

# %%
latex_insts = let cidf = deepcopy(cidf)
    #:instance_name, :N, :n_, :L, :W, :kp_tub, :kp_rel, :kp_opt, :csp_tub, :csp_rel, :csp_opt
    cidf[!, :kp_tub] = trunc.(Int, cidf[!, :kp_tub]) # the decimal places are not relevant here
    cidf[!, :kp_rel] = trunc.(Int, cidf[!, :kp_rel]) # the decimal places are not relevant here

    for col in names(cidf)
        col == "instance_name" && continue
        cidf[!, col] = number2latex.(cidf[!, col]; enclose = false)
    end

    cidf
end

pretty_table(
    latex_insts; backend = :latex, nosubheader = true,
    alignment = vcat([:l], repeat([:r], ncol(latex_insts) - 1))
)

# %%
# Create a code for identifying the distinct runs for the same instance.

#= Version considering all combinations of PSN = #
function build_variant(threads, faithful, psn, problem)
    e, n = convert.(Int, (!faithful, psn))
    p = problem ? "k" : "c"
    return "$(threads)t_$(e)e_$(n)n_$(p)"
end

cdf[!, :variant] = build_variant.(
    getindex.((cdf,), (!,), [:threads, :faithful, :round2disc, :g2kp])...
)
# =#

#= Version considering only enhanced+PSN e faithful-PSN =#
function build_variant(threads, faithful, psn, problem)
    e, n = convert.(Int, (!faithful, psn))
    @assert e == n
    p = problem ? "k" : "c"
    return "$(threads)t_$(e)e_$(p)"
end

cdf[!, :variant] = build_variant.(
    getindex.((cdf,), (!,), [:threads, :faithful, :round2disc, :g2kp])...
)
# =#

# %%
select(cdf, :instance_name, :faithful, :threads, :g2kp, :g2csp, :total_instance_time, :rnr_solved_by_simplex) |>
    df -> filter(r -> r.rnr_solved_by_simplex && r.threads == 12, df) |>
    df -> select(df, :instance_name, :faithful, :g2kp, :g2csp, :total_instance_time) |>
    df -> sort(df, :instance_name) |>
    df -> groupby(df, :instance_name) |>
    df -> @based_on(df, faithful = sum(:faithful), g2kp = sum(:g2kp), g2csp = sum(:g2csp)) |>
    showtable
    #df -> filter(r -> r.rnr_solved_by_simplex, df)# |>
    #df -> groupby(df, :instance_name) |>
    #df -> based_on(df, )

# %%
@assert all(ifelse.(cdf[!, :had_timeout], .!cdf[!, :finished], true))

# %%
# Check how many did not finish in each group.
@linq cdf |> groupby([:round2disc, :faithful, :threads, :g2kp]) |>
    based_on(; qt_total = length(:instance_name), qt_not_finished = sum(.!:finished)) |>
    sort([:threads, :faithful, :round2disc, :g2kp]) |>
    select!(:threads, :faithful, :round2disc, :g2kp, :qt_total, :qt_not_finished)

# %%
# Check which instances did not finish the "normal" way.
let cdf = deepcopy(cdf)
    atypical = @linq cdf |> where((.!:finished) .| :had_timeout) |>
        select(:instance_name, :threads, :faithful, :round2disc, :g2kp, :build_and_solve_time, :datafile) |>
        sort([:instance_name, :threads, :faithful, :round2disc, :g2kp])
    showtable(atypical)
end

# %%
# Just check if the optimals and the relaxations agree.
let cdf = deepcopy(cdf)
    problem = :g2csp # :g2kp #
    my_field = :solution_value # :solution_profit #
    relaxations = filter(row -> row[:finished] .& row[problem], cdf) |>
        df -> select(df, :instance_name, :variant, :final_root_relaxation_value) |>
        df -> unstack(df, :instance_name, :variant, :final_root_relaxation_value)
    optimals = filter(row -> row[:finished] .& row[problem], cdf) |>
        df -> select(df, :instance_name, :variant, my_field) |>
        df -> unstack(df, :instance_name, :variant, my_field)
    showtable(optimals)
    #showtable(relaxations)
    opt = map(eachrow(optimals)) do row
        instance = row[:instance_name]
        v = filter(!ismissing, collect(values(row[Not(:instance_name)])))
        #println(instance, " -- ", maximum(v) - minimum(v))
        @assert all(iszero(maximum(v) - minimum(v)))
        instance => first(v)
    end

    #=differences = map(eachrow(relaxations)) do row
        instance = row[:instance_name]
        v = filter(!ismissing, collect(values(row[Not(:instance_name)])))
        #println(instance, " -- ", maximum(v) - minimum(v))
        instance => maximum(v) - minimum(v)
    end
    rel_diff_by_opt = sort(first.(opt)) .=> (last.(sort(differences, by = first)) ./
        last.(sort(opt, by = first)))
    sort!(rel_diff_by_opt, by = last, rev = true)=#
end

# %%
latex_times = let cdf = deepcopy(cdf)
    cdf[!, :build_and_solve_time] = @with(cdf,
        ifelse.(:finished, :build_and_solve_time, missing)
    )
    times = filter(row -> row[:finished], cdf) |>
        df -> select(df, :instance_name, :variant, :build_and_solve_time) |>
        df -> unstack(df, :instance_name, :variant, :build_and_solve_time)
    rnr_simplex = cdf |>
        df -> select(df, :instance_name, :variant, :rnr_solved_by_simplex) |>
        df -> unstack(df, :instance_name, :variant, :rnr_solved_by_simplex)
    showtable(rnr_simplex)
    # Put columns in the desired order. rnr_simplex is refereed by column and
    # does not need that.
    col_order = [
        "$(t)t_$(e)e_$(p)" for p in ("k", "c") for e in (0, 1) for t in (1, 12)
    ]
    @show col_order
    select!(times, :instance_name, col_order...)

    for col in names(times)
        col == "instance_name" && continue
        times[!, col] = number2latex.(times[!, col]; enclose = false)
        times[!, col] = ifelse.(rnr_simplex[!, col],
            wrap_in_textit.(times[!, col]),
            times[!, col]
        )
    end
    showtable(times)

    highlight_best_values!(times; columns = 2:5, cleaner = dirt2number)
    highlight_best_values!(times; columns = 6:9, cleaner = dirt2number)

    sort!(times, :instance_name)
end

pretty_table(
    latex_times; backend = :latex, nosubheader = true, alignment = vcat([:l], repeat([:r], ncol(latex_times) - 1))
)


text:
# Autogenerated wrapper script for boostpython_jll for x86_64-apple-darwin
export libboost_python

using boost_jll
using Python_jll
JLLWrappers.@generate_wrapper_header("boostpython")
JLLWrappers.@declare_library_product(libboost_python, "@rpath/libboost_python38.dylib")
function __init__()
    JLLWrappers.@generate_init_header(boost_jll, Python_jll)
    JLLWrappers.@init_library_product(
        libboost_python,
        "lib/libboost_python38.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test
using Piwall

top_left = Piwall.Screen(1920,1080)
top_right = Piwall.Screen(1920,1080)
bottom_left = Piwall.Screen(1920,1080)
bottom_right = Piwall.Screen(1920,1080)

wall_2x2 = Piwall.Wall([top_left, top_right, bottom_left, bottom_right])

println("Piwall tests completed.:wall_2x2")

mypi = 3.1415926535

@test mypi == 3.1415926535

@test mypi ≈ 3.141592654

@test isapprox(mypi, 3.14159265)


text:
import Statistics: mean, cov

export  InflationType, IdentityInflation, AdditiveInflation,
        MultiplicativeInflation, MultiAddInflation,
        exactn

"""
        exactn(N)

A function to create a 1D sample with exactly mean 0 and covariance 1
(The samples are no longer i.i.d but this can be usueful when the initialization of the problem is challenging.)
"""
function exactn(N)
    a = deepcopy(randn(N))
    return (a .- mean(a))./std(a)
end

"""
#     InflationType
#
# An abstract type for Inflation.
"""
abstract type InflationType end

"""
    IdentityInflation <: InflationType

An type to store identity inflation :

Define additive inflation: x <- x
"""
struct IdentityInflation <: InflationType
end

"""
        (A::IdentityInflation)(X)

Apply an `IdentityInflation` `A` on an ensemble matrix `X`, i.e. xⁱ -> xⁱ
"""
function (A::IdentityInflation)(X)
    nothing
end

"""
        AdditiveInflation <: InflationType

An type to store additive inflation :

Define additive inflation: x <- x + ϵ with ϵ a random vector
drawn from the distribution α

## Fields:
$(TYPEDFIELDS)

## Constructors
- `AdditiveInflation(Nx::Int64, α::ContinuousMultivariateDistribution)`
- `AdditiveInflation(Nx::Int64)`
- `AdditiveInflation(Nx::Int64, m::Array{Float64,1}, Σ::Union{Array{Float64,2}, Diagonal{Float64}})`
- `AdditiveInflation(Nx::Int64, m::Array{Float64,1}, σ::Array{Float64,1})`
- `AdditiveInflation(Nx::Int64, m::Array{Float64,1}, σ::Float64)`
- `AdditiveInflation(m::Array{Float64,1}, σ::Float64)`
"""

struct AdditiveInflation <: InflationType
    "Dimension of the state vector"
    Nx::Int64

    "Mean of the additive inflation"
    m::Array{Float64,1}

    "Covariance of the additive inflation"
    Σ::Union{Array{Float64,2}, Diagonal{Float64}}

    "Square-root of the covariance matrix"
    σ::Union{Array{Float64,2}, Diagonal{Float64}}
end

# Some convenient constructors for multivariate Gaussian distributions
# By default, the distribution of the additive inflation α is a multivariate
 # normal distribution with zero mean and identity as the covariance matrix
@inline AdditiveInflation(Nx::Int64) = AdditiveInflation(Nx, zeros(Nx),  Diagonal(ones(Nx)), Diagonal(ones(Nx)))

function AdditiveInflation(Nx::Int64, m::Array{Float64,1}, Σ::Union{Array{Float64,2}, Diagonal{Float64}})
@assert Nx==size(m,1) "Error dimension of the mean"
@assert Nx==size(Σ,1)==size(Σ,2) "Error dimension of the covariance matrix"

return AdditiveInflation(Nx, m, Σ, sqrt(Σ))

end

function AdditiveInflation(Nx::Int64, m::Array{Float64,1}, σ::Array{Float64,1})
@assert Nx==size(m,1) "Error dimension of the mean"
@assert Nx==size(σ,1) "Error dimension of the std vector"

return AdditiveInflation(Nx, m, Diagonal(σ .^2), Diagonal(σ))

end

function AdditiveInflation(Nx::Int64, m::Array{Float64,1}, σ::Float64)
@assert Nx==size(m,1) "Error dimension of the mean"

return AdditiveInflation(Nx, m, Diagonal(σ^2*ones(Nx)), Diagonal(σ*ones(Nx)))

end

function AdditiveInflation(m::Array{Float64,1}, σ::Float64)
    Nx = size(m,1)
    return AdditiveInflation(Nx, m, Diagonal(σ^2*ones(Nx)), Diagonal(σ*ones(Nx)))
end

Base.size(A::AdditiveInflation)  = A.Nx
mean(A::AdditiveInflation) = A.m
cov(A::AdditiveInflation) = A.Σ

"""
        (A::AdditiveInflation)(X, start::Int64, final::Int64)

Apply the additive inflation `A` to the lines `start` to `final` of an ensemble matrix `X`,
i.e. xⁱ -> xⁱ + ϵⁱ with ϵⁱ ∼ `A.α`.
"""
function (A::AdditiveInflation)(X, start::Int64, final::Int64; laplace::Bool=false)
    Ne = size(X,2)
    @assert A.Nx == final - start + 1 "final-start + 1 doesn't match the length of the additive noise"
    # @show X[start:final, 1]
    if laplace == false
        @inbounds for i=1:Ne
            col = view(X, start:final, i)
            col .+= A.m + A.σ*randn(A.Nx)
        end
    else
        @inbounds for i=1:Ne
            col = view(X, start:final, i)
            col .+= A.m + sqrt(2.0)*A.σ*randn(Laplace(), A.Nx)
        end
    end
    # @show X[start:final, 1]
end

"""
        (A::AdditiveInflation)(X)

Apply the additive inflation `A` to an ensemble matrix `X`,
i.e. xⁱ -> xⁱ + ϵⁱ with ϵⁱ ∼ `A.α`.
"""
(A::AdditiveInflation)(X; laplace::Bool=false) = A(X, 1, size(X,1); laplace = laplace)

"""
        (A::AdditiveInflation)(x::Array{Float64,1})

Apply the additive inflation `A` to the vector `x`,
i.e. x -> x + ϵ with ϵ ∼ `A.α`.
"""
function (A::AdditiveInflation)(x::Array{Float64,1})
    x .+= A.m + A.σ*randn(A.Nx)
    return x
end

"""
    MultiplicativeInflation <: InflationType

An type to store multiplicative inflation :

xⁱ -> xⁱ + β*(xⁱ - x̄) with β a scalar

# Fields:
- 'β' : multiplicative inflation factor
"""
struct MultiplicativeInflation <: InflationType
    "Multiplicative inflation factor β"
    β::Real
end

"""
        (A::MultiplicativeInflation)(X, start::Int64, final::Int64)

Apply the multiplicative inflation `A` to the lines `start` to `final` of an ensemble matrix `X`,
i.e. xⁱ -> x̄ + β*(xⁱ - x̄)  with β scalar, usually ∼ 1.0.
"""
function (A::MultiplicativeInflation)(X, start::Int64, final::Int64)
    Ne = size(X,2)
    X̂ = copy(mean(view(X, start:final,:), dims = 2)[:,1])
    @inbounds for i=1:Ne
        col = view(X, start:final, i)
        col .= (1.0-A.β)*X̂ + A.β*col
    end
end

"""
        (A::MultiplicativeInflation)(X)

Apply the multiplicative inflation `A` to an ensemble matrix `X`,
i.e. xⁱ -> x̄ + β*(xⁱ - x̄)  with β scalar, usually ∼ 1.0.
"""
(A::MultiplicativeInflation)(X) = A(X, 1, size(X, 1))

"""
    MultiAddInflation <: InflationType

An type to store multiplico-additive inflation :

Define multiplico-additive inflation: xⁱ -> x̄ + β*(xⁱ - x̄)  + ϵⁱ with ϵⁱ ∼ α and β a scalar

## Fields:
- `Nx` : dimension of the vector
- 'β' : Multiplicative inflation factor
- 'α' : Distribution of the additive inflation

## Constructors:
- `MultiAddInflation(Nx::Int64, β::Real, α::ContinuousMultivariateDistribution)`
- `MultiAddInflation(Nx::Int64, β::Float64, m::Array{Float64,1}, Σ)`
- `MultiAddInflation(Nx::Int64, β::Float64, m::Array{Float64,1}, σ::Array{Float64,1})`
- `MultiAddInflation(Nx::Int64, β::Float64, m::Array{Float64,1}, σ::Float64)`
"""
struct MultiAddInflation <: InflationType
    "Dimension of the state vector"
    Nx::Int64

    "Multiplicative inflation factor β"
    β::Real

    "Mean of the additive inflation"
    m::Array{Float64,1}

    "Covariance of the additive inflation"
    Σ::Union{Array{Float64,2}, Diagonal{Float64}}

    "Square-root of the covariance matrix"
    σ::Union{Array{Float64,2}, Diagonal{Float64}}
end

# Some convenient constructors for multivariate Gaussian additive distributions
# By default, for a Multiplico-additive inflation, the multiplicative inflation
# factor β is set to 1.0, and  α is a  multivariate
# normal distribution with zero mean and identity as the covariance matrix
function MultiAddInflation(Nx::Int)
    return MultiAddInflation(Nx, 1.0, zeros(Nx), Diagonal(ones(Nx)), Diagonal(ones(Nx)))
end

function MultiAddInflation(Nx::Int64, β::Float64, m::Array{Float64,1}, Σ)
    @assert β>0.0 "The multiplicative inflation must be >0.0"
    @assert Nx==size(m,1) "Error dimension of the mean"
    @assert Nx==size(Σ,1)==size(Σ,2) "Error dimension of the covariance matrix"

    return MultiAddInflation(Nx, β, m, Σ, sqrt(Σ))
end

function MultiAddInflation(Nx::Int64, β::Float64, m::Array{Float64,1}, σ::Array{Float64,1})
    @assert β>0.0 "The multiplicative inflation must be >0.0"
    @assert Nx==size(m,1) "Error dimension of the mean"
    @assert Nx==size(σ,1) "Error dimension of the std vector"

    return MultiAddInflation(Nx, β, m, Diagonal(σ .^2), Diagonal(σ))
end

function MultiAddInflation(Nx::Int64, β::Float64, m::Array{Float64,1}, σ::Float64)
    @assert β>0.0 "The multiplicative inflation must be >0.0"
    @assert Nx==size(m,1) "Error dimension of the mean"

    return MultiAddInflation(Nx, β, m, Diagonal(σ^2*ones(Nx)), Diagonal(σ*ones(Nx)))
end

"""
    size(A::MultiAddInflation)

Return the dimension of the additive inflation of `A`.
"""
Base.size(A::MultiAddInflation) = A.Nx
mean(A::MultiAddInflation) = A.m
cov(A::MultiAddInflation)  = A.Σ

"""
        (A::MultiAddInflation)(X, start::Int64, final::Int64)

Apply the multiplicat inflation `A` to the lines `start` to `final` of an ensemble matrix `X`,
i.e. xⁱ -> x̄ + β*(xⁱ - x̄)  + ϵⁱ with ϵⁱ ∼ `A.α` and β a scalar, usually ∼ 1.0.
"""
function (A::MultiAddInflation)(X, start::Int64, final::Int64)
    @assert A.Nx = final - start +1 "Dimension does not match"
    Ne = size(X,2)
    X̂ = copy(mean(view(X, start:final,:), dims = 2)[:,1])
    @inbounds for i=1:Ne
        col = view(X, start:final, i)
        col .= (1.0-A.β)*X̂ + A.β*col + A.m + A.σ*randn(A.Nx)
    end
end

"""
        (A::MultiAddInflation)(X, start::Int64, final::Int64)

Apply the multiplico-additive inflation `A` to the ensemble matrix `X`,
i.e. xⁱ -> x̄ + β*(xⁱ - x̄)  + ϵⁱ with ϵⁱ ∼ `A.α` and β a scalar, usually ∼ 1.0.
"""
(A::MultiAddInflation)(X) = A(X, 1, size(X,1))


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.1
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/is-there-a-maximum-f-op-itrs-in-julia/63868/3

# %%
using BenchmarkTools

a = randn(1000, 1000)
b = randn(1000, 1000)

f(a, b) = maximum(abs, x - y for (x, y) in zip(a, b))
g(a, b) = maximum(Base.splat(abs∘-), zip(a, b))
h(a, b) = mapreduce(abs∘-, max, a, b)

@show f(a, b) == g(a, b) == h(a, b)
@btime f($a, $b)
@btime g($a, $b)
@btime h($a, $b)

# %%

max_stars_repo_path: 0008/Base.splat example - maximum of abs differences.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
@doc raw"""
    NegBinomialLikelihood(r::Real)

## Arguments
- `r::Real` number of failures until the experiment is stopped

---

[Negative Binomial likelihood](https://en.wikipedia.org/wiki/Negative_binomial_distribution) with number of failures `r`
```math
    p(y|r, f) = {y + r - 1 \choose y} (1 - \sigma(f))^r \sigma(f)^y,
```
if ``r\in \mathbb{N}`` or
```math
    p(y|r, f) = \frac{\Gamma(y + r)}{\Gamma(y + 1)\Gamma(r)} (1 - \sigma(f))^r \sigma(f)^y,
```
if ``r\in\mathbb{R}``.
Where ``\sigma`` is the logistic function
"""
struct NegBinomialLikelihood{L,Tr} <: AbstractLikelihood
    invlink::L
    r::Tr
end

NegBinomialLikelihood(r) = NegBinomialLikelihood(LogisticLink(), r)

(l::NegBinomialLikelihood)(f::Real) = NegativeBinomial(l.r, l.invlink(f))

implemented(::NegBinomialLikelihood, ::Union{<:AnalyticVI,<:GibbsSampling}) = true

function (l::NegBinomialLikelihood)(y::Real, f::Real)
    return pdf(l(f), y)
end

function Distributions.loglikelihood(l::NegBinomialLikelihood, y::Real, f::Real)
    return logpdf(l(f), y)
end

function Base.show(io::IO, l::NegBinomialLikelihood)
    return print(io, "Negative Binomial Likelihood (r = $(l.r))")
end

function compute_proba(
    l::NegBinomialLikelihood, μ::AbstractVector{<:Real}, σ²::AbstractVector{<:Real}
)
    N = length(μ)
    T = eltype(μ)
    pred = zeros(T, N)
    sig_pred = zeros(T, N)
    for i in 1:N
        x = pred_nodes .* sqrt(max(σ²[i], zero(T))) .+ μ[i]
        pred[i] = dot(pred_weights, l.invlink.(x))
        sig_pred[i] = dot(pred_weights, l.invlink.(x) .^ 2) - pred[i]^2
    end
    return pred, sig_pred
end

## Local Updates ##
function init_local_vars(::NegBinomialLikelihood, batchsize::Int, T::DataType=Float64)
    return (; c=rand(T, batchsize), θ=zeros(T, batchsize))
end

function local_updates!(
    local_vars,
    l::NegBinomialLikelihood,
    y::AbstractVector,
    μ::AbstractVector,
    diagΣ::AbstractVector,
)
    map!(sqrt_expec_square, local_vars.c, μ, diagΣ) # √E[f^2]
    map!(local_vars.θ, local_vars.c, y) do c, y
        (l.r + y) * tanh(c / 2) / c
    end # E[ω]
    return local_vars
end

function sample_local!(
    local_vars, l::NegBinomialLikelihood, y::AbstractVector, f::AbstractVector
)
    map!(local_vars.θ, f, y) do f, y
        rand(PolyaGamma(y + Int(l.r), abs(f)))
    end
    return local_vars
end

## Global Updates ##

@inline function ∇E_μ(l::NegBinomialLikelihood, ::AOptimizer, y::AbstractVector, state)
    return ((y .- l.r) / 2,)
end
@inline function ∇E_Σ(::NegBinomialLikelihood, ::AOptimizer, y::AbstractVector, state)
    return (state.θ / 2,)
end

## ELBO Section ##

AugmentedKL(l::NegBinomialLikelihood, state, y) = PolyaGammaKL(l, state, y)

function logabsbinomial(n, k)
    return log(binomial(n, k))
end

function negbin_logconst(y, r::Real)
    return loggamma.(y .+ r) - loggamma.(y .+ 1) .- loggamma(r)
end

function negbin_logconst(y, r::Int)
    return logabsbinomial.(y .+ (r - 1), y)
end

function expec_loglikelihood(
    l::NegBinomialLikelihood,
    ::AnalyticVI,
    y::AbstractVector,
    μ::AbstractVector,
    diag_cov::AbstractVector,
    state,
)
    tot = Zygote.@ignore(sum(negbin_logconst(y, l.r))) - log(2.0) * sum(y .+ l.r)
    tot += dot(μ, (y .- l.r)) / 2 - dot(state.θ, μ) / 2 - dot(state.θ, diag_cov) / 2
    return tot
end

function PolyaGammaKL(l::NegBinomialLikelihood, state, y)
    return PolyaGammaKL(y .+ l.r, state.c, state.θ)
end


text:
module BuilderExample

# Car parts
struct Engine
    value
end

struct Wheels
    value
end

struct Chassis
    value
end

# Easy syntax for Union{T, Nothing}
const Maybe{T} = Union{T,Nothing}

# Car object
mutable struct Car
    engine::Maybe{Engine}
    wheels::Maybe{Wheels}
    chassis::Maybe{Chassis}
end

Car() = Car(nothing, nothing, nothing)

function add(wheels::Wheels)
    return function (c::Car)
        c.wheels = wheels
        return c
    end
end

function add(engine::Engine)
    return function (c::Car)
        c.engine = engine
        return c
    end
end

function add(chassis::Chassis)
    return function (c::Car)
        c.chassis = chassis
        return c
    end
end

function test1()
    car = Car() |>
        add(Engine("4-cylinder 1600cc Engine")) |>
        add(Wheels("4x 20-inch wide wheels")) |>
        add(Chassis("Roadster Chassis"))
    println(car)
end

function test2()
    car = Car()
    car.engine = Engine("4-cylinder 1600cc Engine")
    car.wheels = Wheels("4x 20-inch wide wheels")
    car.chassis = Chassis("Roadster Chassis")
    println(car)
end

end #module

using .BuilderExample
BuilderExample.test1()
BuilderExample.test2()


text:
# # A Dirichlet Poisson problem

#md # ```@meta
#md # CurrentModule = ImmersedLayers
#md # ```

#=
Now, let's put together a lot of the tools we've discussed earlier
to solve a standard problem: the solution of Laplace's equation,
with Dirichlet boundary conditions on a surface. Our objective is to solve this
problem with one Dirichlet value for the region external to the surface,
and another Dirichlet value internal to the surface.

To be specific, our problem is

$$\nabla^2\varphi^+ = 0,\qquad \nabla^2\varphi^- = 0$$

where $+$ denotes the exterior and $-$ the interior of the surface. (We will
consider a circle of radius 1.) The boundary conditions on this surface are

$$\varphi^+ = x, \qquad \varphi^- = 0$$

In other words, we seek to set the value on the exterior to the local $x$ coordinate
of the surface, while the interior should have zero field.

Discretizing this problem by the techniques developed in [^1], we seek to solve

$$\begin{bmatrix} L & R \\ R^T & 0 \end{bmatrix} \begin{pmatrix} f \\ -s \end{pmatrix} = \begin{pmatrix} D_s d \\ \overline{f}_b \end{pmatrix}$$

where $\overline{f}_b = (f^+_b + f^-_b)/2 = x/2$ and $d = f^+_b - f^-_b = x$.

This saddle-point problem can be solved by block-LU decomposition. First solve

$$L f^{*} = D_s d$$

for $f^*$. Then solve

$$R^TL^{-1}R s = \overline{f}_b - R^T f^{*}$$

for $s$, and finally, compute

$$f = f^{*} + L^{-1}R s$$

We will demonstrate these steps here.
=#

using ImmersedLayers
#!jl using Plots
using LinearAlgebra

#=
### Set up the grid, shape, and cache
We do this just as we did in [Immersed layer caches](@ref)
=#
Δx = 0.01
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx)
RadC = 1.0
Δs = 1.4*cellsize(g)
body = Circle(RadC,Δs)
cache = SurfaceScalarCache(body,g,scaling=GridScaling);

#=
Let's set up the right-hand side values, $D_s d$ and $\overline{f}_b$.
We will use `fstar` as a placeholder for the right-hand side of the
upper equations.
=#
pts = points(cache)
x = zeros_surface(cache)
fb = zeros_surface(cache)
fstar = zeros_grid(cache)
x .= pts.u
surface_divergence!(fstar,x,cache)
fb .= 0.5*x;

#=
First step in the algorithm, find $f^*$
=#
inverse_laplacian!(fstar,cache);

#=
Now, generate the Schur complement matrix, $S = -R^TL^{-1}R$, using
[`create_RTLinvR`](@ref),
=#
S = create_RTLinvR(cache);

#=
Now we set up the right-hand side of the Schur complement system, and solve
for the Lagrange multipliers, $s$. We will use `s` as a placeholder for the
right-hand side, first holding $R^T f^*$, then $\overline{f}_b - R^T f^*$,
and finally, the solution of $- S s = \overline{f}_b - R^T f^*$.
=#
s = zeros_surface(cache)
interpolate!(s,fstar,cache)
s .= fb - s
s .= -(S\s);

#=
Finally, correct $f^*$ with $L^{-1}Rs$
=#
f = zeros_grid(cache)
regularize!(f,s,cache)
inverse_laplacian!(f,cache)
f .+= fstar;

#=
Now, let's look at the resulting field
=#
#!jl plot(f,cache)

#=
We can see that the field external to the circle has the expected behavior,
and the solution inside the circle is blank, indicating it is zero.

Let's also plot the solution for $s$ on the surface. We will plot the
'raw' result, which is rather noisy, and also apply the filtering matrix $C = \tilde{R}^T R$
several times to smooth it out.
=#
C = create_surface_filter(cache)
#!jl plot(s,cache,label="Original solution")
s .= C^5*s;
#!jl plot!(s,cache,label="Filtered solution")

#md # [^1]: Eldredge, J. D. (2021) "A method of immersed layers on Cartesian grids, with application to incompressible flows," arXiv:2103.04521.


text:
export K2_5KrylovParams

"""
Type to use the K2.5 formulation with a Krylov method, using the package
[`Krylov.jl`](https://github.com/JuliaSmoothOptimizers/Krylov.jl).
The outer constructor

    K2_5KrylovParams(; uplo = :L, kmethod = :minres, preconditioner = :Identity,
                     rhs_scale = true,
                     atol0 = 1.0e-4, rtol0 = 1.0e-4,
                     atol_min = 1.0e-10, rtol_min = 1.0e-10,
                     ρ0 = sqrt(eps()) * 1e5, δ0 = sqrt(eps()) * 1e5,
                     ρ_min = 1e2 * sqrt(eps()), δ_min = 1e2 * sqrt(eps()),
                     mem = 20)

creates a [`RipQP.SolverParams`](@ref) that should be used to create a [`RipQP.InputConfig`](@ref).
The available methods are:
- `:minres`
- `:minres_qlp`
- `:symmlq`

The list of available preconditioners for this solver is displayed here: [`RipQP.PreconditionerDataK2`](@ref)
"""
mutable struct K2_5KrylovParams <: AugmentedParams
  uplo::Symbol
  kmethod::Symbol
  preconditioner::Symbol
  rhs_scale::Bool
  atol0::Float64
  rtol0::Float64
  atol_min::Float64
  rtol_min::Float64
  ρ0::Float64
  δ0::Float64
  ρ_min::Float64
  δ_min::Float64
  mem::Int
end

function K2_5KrylovParams(;
  uplo::Symbol = :L,
  kmethod::Symbol = :minres,
  preconditioner::Symbol = :Identity,
  rhs_scale::Bool = true,
  atol0::T = 1.0e-4,
  rtol0::T = 1.0e-4,
  atol_min::T = 1.0e-10,
  rtol_min::T = 1.0e-10,
  ρ0::T = sqrt(eps()) * 1e5,
  δ0::T = sqrt(eps()) * 1e5,
  ρ_min::T = 1e2 * sqrt(eps()),
  δ_min::T = 1e3 * sqrt(eps()),
  mem::Int = 20,
) where {T <: Real}
  return K2_5KrylovParams(
    uplo,
    kmethod,
    preconditioner,
    rhs_scale,
    atol0,
    rtol0,
    atol_min,
    rtol_min,
    ρ0,
    δ0,
    ρ_min,
    δ_min,
    mem,
  )
end

mutable struct PreallocatedDataK2_5Krylov{
  T <: Real,
  S,
  L <: LinearOperator,
  Pr <: PreconditionerDataK2,
  Ksol <: KrylovSolver,
} <: PreallocatedDataAugmentedKrylov{T, S}
  pdat::Pr
  D::S                                  # temporary top-left diagonal
  sqrtX1X2::S # vector to scale K2 to K2.5
  tmp1::S # temporary vector for products
  tmp2::S # temporary vector for products
  rhs::S
  rhs_scale::Bool
  regu::Regularization{T}
  δv::Vector{T}
  K::L # augmented matrix
  KS::Ksol
  atol::T
  rtol::T
  atol_min::T
  rtol_min::T
end

function opK2_5prod!(
  res::AbstractVector{T},
  nvar::Int,
  Q::Union{AbstractMatrix{T}, AbstractLinearOperator{T}},
  D::AbstractVector{T},
  A::Union{AbstractMatrix{T}, AbstractLinearOperator{T}},
  sqrtX1X2::AbstractVector{T},
  tmp1::AbstractVector{T},
  tmp2::AbstractVector{T},
  δv::AbstractVector{T},
  v::AbstractVector{T},
  α::T,
  β::T,
  uplo::Symbol,
) where {T}
  tmp2 .= @views sqrtX1X2 .* v[1:nvar]
  mul!(tmp1, Q, tmp2, -α, zero(T))
  tmp1 .= @views sqrtX1X2 .* tmp1 .+ α .* D .* v[1:nvar]
  if β == zero(T)
    res[1:nvar] .= tmp1
  else
    res[1:nvar] .= @views tmp1 .+ β .* res[1:nvar]
  end
  if uplo == :U
    @views mul!(tmp1, A, v[(nvar + 1):end], α, zero(T))
    res[1:nvar] .+= sqrtX1X2 .* tmp1
    @views mul!(res[(nvar + 1):end], A', tmp2, α, β)
  else
    @views mul!(tmp1, A', v[(nvar + 1):end], α, zero(T))
    res[1:nvar] .+= sqrtX1X2 .* tmp1
    @views mul!(res[(nvar + 1):end], A, tmp2, α, β)
  end
  res[(nvar + 1):end] .+= @views (α * δv[1]) .* v[(nvar + 1):end]
end

function PreallocatedData(
  sp::K2_5KrylovParams,
  fd::QM_FloatData{T},
  id::QM_IntData,
  itd::IterData{T},
  pt::Point{T},
  iconf::InputConfig{Tconf},
) where {T <: Real, Tconf <: Real}

  # init Regularization values
  D = similar(fd.c, id.nvar)
  if iconf.mode == :mono
    regu = Regularization(T(sp.ρ0), T(sp.δ0), T(sp.ρ_min), T(sp.δ_min), :classic)
    D .= -T(1.0e0) / 2
  else
    regu =
      Regularization(T(sp.ρ0), T(sp.δ0), T(sqrt(eps(T)) * 1e0), T(sqrt(eps(T)) * 1e0), :classic)
    D .= -T(1.0e-2)
  end
  sqrtX1X2 = fill!(similar(D), one(T))
  tmp1 = similar(D)
  tmp2 = similar(D)
  δv = [regu.δ] # put it in a Vector so that we can modify it without modifying opK2prod!
  K = LinearOperator(
    T,
    id.nvar + id.ncon,
    id.nvar + id.ncon,
    true,
    true,
    (res, v, α, β) ->
      opK2_5prod!(res, id.nvar, fd.Q, D, fd.A, sqrtX1X2, tmp1, tmp2, δv, v, α, β, fd.uplo),
  )

  rhs = similar(fd.c, id.nvar + id.ncon)

  KS = init_Ksolver(K, rhs, sp)

  pdat = eval(sp.preconditioner)(id, fd, regu, D, K)

  return PreallocatedDataK2_5Krylov(
    pdat,
    D,
    sqrtX1X2,
    tmp1,
    tmp2,
    rhs,
    sp.rhs_scale,
    regu,
    δv,
    K, #K
    KS, #K_fact
    T(sp.atol0),
    T(sp.rtol0),
    T(sp.atol_min),
    T(sp.rtol_min),
  )
end

function solver!(
  dd::AbstractVector{T},
  pad::PreallocatedDataK2_5Krylov{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
  step::Symbol,
) where {T <: Real}

  # erase dda.Δxy_aff only for affine predictor step with PC method
  pad.rhs[1:(id.nvar)] .= @views dd[1:(id.nvar)] .* pad.sqrtX1X2
  pad.rhs[(id.nvar + 1):end] .= @views dd[(id.nvar + 1):end]
  if pad.rhs_scale
    rhsNorm = kscale!(pad.rhs)
  end
  pad.K.nprod = 0
  ksolve!(pad.KS, pad.K, pad.rhs, pad.pdat.P, verbose = 0, atol = pad.atol, rtol = pad.rtol)
  update_kresiduals_history!(res, pad.K, pad.KS.x, pad.rhs)
  if pad.rhs_scale
    kunscale!(pad.KS.x, rhsNorm)
  end
  pad.KS.x[1:(id.nvar)] .*= pad.sqrtX1X2

  dd .= pad.KS.x

  return 0
end

function update_pad!(
  pad::PreallocatedDataK2_5Krylov{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
) where {T <: Real}
  if cnts.k != 0
    update_regu!(pad.regu)
  end

  if pad.atol > pad.atol_min
    pad.atol /= 10
  end
  if pad.rtol > pad.rtol_min
    pad.rtol /= 10
  end

  # K2.5
  pad.sqrtX1X2 .= one(T)
  pad.sqrtX1X2[id.ilow] .*= sqrt.(itd.x_m_lvar)
  pad.sqrtX1X2[id.iupp] .*= sqrt.(itd.uvar_m_x)
  pad.D .= zero(T)
  pad.D[id.ilow] .-= pt.s_l
  pad.D[id.iupp] .*= itd.uvar_m_x
  pad.tmp1 .= zero(T)
  pad.tmp1[id.iupp] .-= pt.s_u
  pad.tmp1[id.ilow] .*= itd.x_m_lvar
  pad.D .+= pad.tmp1 .- pad.regu.ρ

  pad.δv[1] = pad.regu.δ

  update_preconditioner!(pad.pdat, pad, itd, pt, id, fd, cnts)

  return 0
end


text:
using Base.Test

macro wrappedallocs(expr)
    argnames = [gensym() for a in expr.args]
    quote
        function g($(argnames...))
            @allocated $(Expr(expr.head, argnames...))
        end
        $(Expr(:call, :g, [esc(a) for a in expr.args]...))
    end
end

@testset "computed field types" begin
    include("computed/runtests.jl")
end

include("simple_interface.jl")
include("interface_with_parent.jl")

module NotebookTest
    using NBInclude
    @nbinclude(joinpath(@__DIR__, "..", "demo.ipynb"))
end


text:
############################################################
#  Devuelve el índice de tantos cromosomas a buscar
############################################################
"""
    IndiceCromosomasProbabilidad(tamaño, porcentaje_seleccion, valores)
Devuelve lista con los índices a mutar

Tamaño es el tamaño de los cromosomas
porcentaje_seleccion flotante entre (0,1]
valores: no se utilizan para nada, es por manterner la cabecera de ese tipo de funciones
"""
function IndiceCromosomasProbabilidad(tamaño, porcentaje_seleccion, valores)
    cantidad_indices = round(Int, porcentaje_seleccion* tamaño)
    permutaciones = randperm(tamaño)
    return permutaciones[1:cantidad_indices]
end

"""
    MejoresCromosomas(tamaño, porcentaje_seleccion, valores)
Devuelve una lista con los índices a mutar
Concretamemente un subconjunto de tamaño tamaño*porcentaje_seleccion con los mejores cromosomas
"""
function MejoresCromosomas(tamaño, porcentaje_seleccion, valores)
    cantidad_indices = round(Int, porcentaje_seleccion* tamaño)
    permutaciones = Array{Int, }(undef, cantidad_indices)
    valores_aux = copy(valores)

    for i in 1:cantidad_indices
        indice_mejor = argmax(valores_aux)
        permutaciones[i] = indice_mejor
        valores_aux[i] = -1 # valor por debajo del válido
    end
    return permutaciones
end


text:
#=
import Pkg; Pkg.activate(joinpath(@__DIR__, "../../../"))
using CytofDensityEstimation
const CDE = CytofDensityEstimation
=#

function scenarios(n::Int)
  # NOTE: Later.
  # loc = [-1, 0., 1]
  # scale = [1, 1, 1]

  loc = [-1, 1, 3.]
  scale = [.7, .7, .7]
  etaC = [.5, .5, 0]
  df = [7, 5, 10.]
  skew = [-5, -3, 0.]
  Ni = 1000

  if n == 0  # TEST CASE!
    return CDE.Model.generate_samples(NC=100, NT=100, gammaC=0.3, gammaT=0.3,
                                      etaC=[1., 0, 0], etaT=[1., 0, 0],
                                      loc=loc, scale=scale, df=df, skew=skew)
  elseif n == 1
    return CDE.Model.generate_samples(NC=Ni, NT=Ni, gammaC=0.3, gammaT=0.2,
                                      etaC=etaC, etaT=[.5, .4, .1],
                                      loc=loc, scale=scale, df=df, skew=skew)
  elseif n == 2
    return CDE.Model.generate_samples(NC=Ni, NT=Ni, gammaC=0.3, gammaT=0.3,
                                      etaC=etaC, etaT=[.5, .4, .1],
                                      loc=loc, scale=scale, df=df, skew=skew)
  elseif n == 3
    return CDE.Model.generate_samples(NC=Ni, NT=Ni, gammaC=0.3, gammaT=0.3,
                                      etaC=etaC, etaT=[.5, .45, .05],
                                      loc=loc, scale=scale, df=df, skew=skew)
  elseif n == 4
    return CDE.Model.generate_samples(NC=Ni, NT=Ni, gammaC=0.3, gammaT=0.3,
                                      etaC=etaC, etaT=[.5, .5, .0],
                                      loc=loc, scale=scale, df=df, skew=skew)
  else
    throw(ArgumentError("n=$(n) was provided. But n must be an integer between 0 and 5."))
  end
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
#!/usr/bin/env julia

export pe10

# need Primes.jl
# https://github.com/JuliaMath/Primes.jl
using Primes

pe10(n::Integer=2*10^6)=primes(n)|>sum

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe10()|>print
end

max_stars_repo_path: 10. Summation of primes/10.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 10. Summation of primes/10.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 10. Summation of primes/10.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
function specification(P, df, args, attr_var_pairs)
	dict = Dict(attr_var_pairs...)

	# 0. Take out some special attributes
	layout_vars = let
		grp_x    = pop!(dict, :layout_x, nothing)
		grp_y    = pop!(dict, :layout_y, nothing)
		grp_wrap = pop!(dict, :layout_wrap, nothing)
		linkxaxes  = pop!(dict, :linkxaxes, true)
		linkyaxes  = pop!(dict, :linkyaxes, true)
		linkzcolor = pop!(dict, :linkzcolor, true)

		(; grp_x, grp_y, grp_wrap, linkxaxes, linkyaxes, linkzcolor)
	end

	# 1. Grouping
	@unpack group_pairs, style_pairs, kws = group_style_other(df, dict)

	group_dict = build_group_dict(df, group_pairs)
	style_dict = build_style_dict(df, style_pairs)

	(; P, df, args, layout_vars, group_dict, style_dict, kws, group_pairs, style_pairs)
end

specification(P, df, args...; attr_var_pairs...) = specification(P, df, args, attr_var_pairs)

function build_group_dict(df, group_pairs)
	group_dict = Dict()
	for (attr, var) in pairs(group_pairs)
		if attr ∉ [:group]
			var_levels = categorical_labels(get(df, var))
			thm = AbstractPlotting.current_default_theme().palette

			if haskey(thm, attr)
				group_dict[attr] = [var_levels[i] => thm[attr][][i] for i in categorical_range(get(df, var))]
			else
				group_dict[attr] = [var_levels[i] => i for i in categorical_range(get(df, var))] # TODO: is is just DataAPI.invrefpool, isn't it?
			end
		end
	end

	group_dict
end

function build_style_dict(df, style_pairs)
	style_dict = Dict()
	for (attr, var) in pairs(style_pairs)
		var_extr = extrema(get(df, var))

		style_dict[attr] = var_extr
	end

	style_dict
end

function group_style_other(df, dict)
	group_ = Dict()
	style_ = Dict()
	kws_   = Dict()

	for (attr, var) in pairs(dict)
		if iscolumn(df, var)
			if is_discrete(get(df, var))
	  			group_[attr] = var
			else
				style_[attr] = var
			end
		else
			pop!(dict, attr)
			kws_[attr] = var
		end
	end

	group_pairs = (; group_...)
	style_pairs = (; style_...)
	kws   = (; kws_...)

	(; group_pairs, style_pairs, kws)
end

get_marker(x, marker_dict) = last.(marker_dict)[categorical_positions(x)]

function unique_or_identity(x)
	x_unique = unique(x)
	if length(x_unique) == 1
		return only(x_unique)
	else
		return x
	end
end

is_discrete(x) = !(categorical_trait(x) isa Continuous)

"This function replaces group indicators (numbers, categories) by attributes that can be plotted (:solid, :red, etc...)"
function lookup_symbols(df, group_pairs, style_pairs, group_dict)
	group_style_pairs = (; group_pairs..., style_pairs...)

	map(collect(pairs(group_style_pairs))) do (attr, var)
		x0 = get(df, var)
		if is_discrete(x0)
			x = get_marker(x0, group_dict[attr]) |> unique_or_identity
		else
			x = x0
		end

		attr => x
	end
end

max_stars_repo_name: greimel/TabularMakie

max_issues_repo_name: greimel/TabularMakie

max_forks_repo_name: greimel/TabularMakie


text:
function Base.show(io::IO, ::MIME"text/plain", x::AbstractExplicitContainerLayer)
    if get(io, :typeinfo, nothing) === nothing  # e.g. top level in REPL
        _big_show(io, x)
    elseif !get(io, :compact, false)  # e.g. printed inside a Vector, but not a Matrix
        _layer_show(io, x)
    else
        show(io, x)
    end
end

function _big_show(io::IO, obj, indent::Int=0, name=nothing)
    pre, post = "(", ")"
    children = _get_children(obj)
    if obj isa Function
        println(io, " "^indent, isnothing(name) ? "" : "$name = ", obj)
    elseif all(_show_leaflike, children)
        _layer_show(io, obj, indent, name)
    else
        println(io, " "^indent, isnothing(name) ? "" : "$name = ", nameof(typeof(obj)), pre)
        if obj isa Chain{<:NamedTuple}
            for k in Base.keys(obj)
                _big_show(io, obj.layers[k], indent + 4, k)
            end
        elseif obj isa Parallel{<:Any, <:NamedTuple}
            if obj.connection !== nothing
                _big_show(io, obj.connection, indent + 4)
            end
            for k in Base.keys(obj)
                _big_show(io, obj.layers[k], indent + 4, k)
            end
        elseif obj isa PairwiseFusion
            _big_show(io, obj.connection, indent + 4)
            for k in Base.keys(obj)
                _big_show(io, obj.layers[k], indent + 4, k)
            end
        elseif obj isa BranchLayer
            for k in Base.keys(obj)
                _big_show(io, obj.layers[k], indent + 4, k)
            end
        else
            if children isa NamedTuple
                for (k, c) in pairs(children)
                    _big_show(io, c, indent + 4, k)
                end
            else
                for c in children
                    _big_show(io, c, indent + 4)
                end
            end
        end
        if indent == 0  # i.e. this is the outermost container
            print(io, rpad(post, 2))
            _big_finale(io, obj)
        else
            println(io, " "^indent, post, ",")
        end
    end
end

_show_leaflike(x) = Functors.isleaf(x)  # mostly follow Functors, except for:
_show_leaflike(x::AbstractExplicitLayer) = false
_show_leaflike(::Tuple{Vararg{<:Number}}) = true         # e.g. stride of Conv
_show_leaflike(::Tuple{Vararg{<:AbstractArray}}) = true  # e.g. parameters of LSTMcell

function _get_children(l::AbstractExplicitContainerLayer{names}) where {names}
    return NamedTuple{names}(getfield.((l,), names))
end
function _get_children(p::Parallel)
    p.connection === nothing ? p.layers : (p.connection, p.layers...)
end
_get_children(s::SkipConnection) = (s.layers, s.connection)
_get_children(s::WeightNorm) = (s.layer,)
_get_children(::Any) = ()

function Base.show(io::IO, ::MIME"text/plain", x::AbstractExplicitLayer)
    if !get(io, :compact, false)
        _layer_show(io, x)
    else
        show(io, x)
    end
end

function _layer_show(io::IO, layer, indent::Int=0, name=nothing)
    _str = isnothing(name) ? "" : "$name = "
    str = _str * sprint(show, layer; context=io)
    print(io, " "^indent, str, indent == 0 ? "" : ",")
    paramlength = parameterlength(layer)
    if paramlength > 0
        print(io, " "^max(2, (indent == 0 ? 20 : 39) - indent - length(str)))
        printstyled(io, "# ", underscorise(paramlength), " parameters";
                    color=:light_black)
        nonparam = statelength(layer)
        if nonparam > 0
            printstyled(io, ", plus ", underscorise(nonparam),
                        indent == 0 ? " non-trainable" : ""; color=:light_black)
        end
    end
    return indent == 0 || println(io)
end

function _big_finale(io::IO, m)
    paramlength = parameterlength(m)
    nonparamlength = statelength(m)
    pars = underscorise(paramlength)
    bytes = Base.format_bytes(Base.summarysize(m))
    nonparam = underscorise(nonparamlength)
    printstyled(io, " "^08, "# Total: "; color=:light_black)
    println(io, pars, " parameters,")
    printstyled(io, " "^10, "#        plus "; color=:light_black)
    print(io, nonparam, " states, ")
    printstyled(io, "summarysize "; color=:light_black)
    print(io, bytes, ".")
    return
end

# utility functions

function underscorise(n::Integer)
    join(reverse(join.(reverse.(Iterators.partition(digits(n), 3)))), '_')
end

function _nan_show(io::IO, x)
    if !isempty(x) && _all(iszero, x)
        printstyled(io, "  (all zero)"; color=:cyan)
    elseif _any(isnan, x)
        printstyled(io, "  (some NaN)"; color=:red)
    elseif _any(isinf, x)
        printstyled(io, "  (some Inf)"; color=:red)
    end
end

_any(f, xs::AbstractArray{<:Number}) = any(f, xs)
_any(f, xs) = any(x -> _any(f, x), xs)
_any(f, x::Number) = f(x)

_all(f, xs) = !_any(!f, xs)


text:
export
ruleSPBetaOutNMM,
ruleSPBetaOutMCNMM,
ruleSPBetaMNM,
ruleSPBetaMMN,
ruleVBBetaOut,
ruleVBBetaA,
ruleVBBetaB

ruleSPBetaOutNMM(msg_out::Nothing, msg_a::Message{PointMass, Univariate}, msg_b::Message{PointMass, Univariate}) = Message(Univariate, Beta, a=deepcopy(msg_a.dist.params[:m]), b=deepcopy(msg_b.dist.params[:m]))

function ruleSPBetaOutMCNMM(msg_out::Nothing, msg_a::Message{F1, Univariate}, msg_b::Message{F2, Univariate}) where {F1<:FactorNode, F2<:FactorNode}
    n_samples = default_n_samples
    samples_a = sample(msg_a.dist,n_samples)
    samples_b = sample(msg_b.dist,n_samples)
    s_list = betainvcdf.(samples_a, samples_b, rand(n_samples))
    w_list = ones(n_samples)/n_samples
    Message(Univariate,SampleList,s=s_list,w=w_list)
end

# 1000 (default_n_samples) samples is chosen arbitrarily. Better to be defined by the user.
function ruleSPBetaMNM(msg_out::Message{F1,Univariate}, msg_a::Nothing, msg_b::Message{F2,Univariate}) where {F1<:FactorNode, F2<:FactorNode}
    n_samples = default_n_samples
    samples_out = sample(msg_out.dist,n_samples)
    samples_b = sample(msg_b.dist,n_samples)
    logp(x,a,b) = (a-1)*log(x) + (b-1)*log(1.0-x) - loggamma(a) - loggamma(b) + loggamma(a+b)
    logp(a) = sum(logp.(samples_out,a,samples_b))/n_samples
    Message(Univariate, Function, log_pdf = logp)
end

function ruleSPBetaMMN(msg_out::Message{F1,Univariate}, msg_a::Message{F2,Univariate}, msg_b::Nothing) where {F1<:FactorNode, F2<:FactorNode}
    n_samples = default_n_samples
    samples_out = sample(msg_out.dist,n_samples)
    samples_a = sample(msg_a.dist,n_samples)
    logp(x,a,b) = (a-1)*log(x) + (b-1)*log(1.0-x) - loggamma(a) - loggamma(b) + loggamma(a+b)
    logp(b) = sum(logp.(samples_out,samples_a,b))/n_samples
    Message(Univariate, Function, log_pdf = logp)
end

ruleVBBetaOut(marg_out::Any, dist_a::ProbabilityDistribution{Univariate}, dist_b::ProbabilityDistribution{Univariate}) = Message(Univariate, Beta, a=unsafeMean(dist_a), b=unsafeMean(dist_b))

function ruleVBBetaA( dist_out::ProbabilityDistribution{Univariate},
    dist_a::Any, dist_b::ProbabilityDistribution{Univariate})

    logGamAsumB(a) = sum(loggamma.(a.+sample(dist_b,default_n_samples)))/default_n_samples
    logp(a) = (a-1)*unsafeLogMean(dist_out) + logGamAsumB(a) - loggamma(a)

    Message(Univariate, Function, log_pdf = logp)
end

function ruleVBBetaB( dist_out::ProbabilityDistribution{Univariate},
    dist_a::ProbabilityDistribution{Univariate}, dist_b::Any)

    logGamAsumB(b) = sum(loggamma.(b.+sample(dist_a,default_n_samples)))/default_n_samples
    logp(b) = (b-1)*unsafeMirroredLogMean(dist_out) + logGamAsumB(b) - loggamma(b)

    Message(Univariate, Function, log_pdf = logp)
end


text:
module Trilinos

# Exports for package modules
export Belos
export Ifpack2
export Kokkos
export MueLu
export Teuchos
export Tpetra
export Thyra

# Exports for extra functionality
export TpetraSolver

using CxxWrap

# Load the deps
const depsfile = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if !isfile(depsfile)
  error("$depsfile not found, package Trilinos did not build properly")
end
include(depsfile)

# Overload for size_t
@static if Sys.WORD_SIZE == 64
  CxxWrap.argument_overloads(t::Type{UInt64}) = [Int64]
end

const CxxUnion{T} = Union{T,CxxWrap.SmartPointer{T}}

include("Teuchos.jl")
include("Kokkos.jl")
include("Tpetra.jl")
include("Belos.jl")
include("Thyra.jl")
include("Ifpack2.jl")
include("MueLu.jl")
include("Benchmark.jl")
include("Testing.jl")

include("Solve.jl")

end # module


text:
module MerkleTrees

export MerkleTree,
       merkle_parent, merkle_parent_level, merkle_root,
       root, populate!

include("hash.jl")

"""
    merkle_parent(hash1::Vector{UInt8}, hash2::Vector{UInt8})
    -> Vector{UInt8}

Takes the binary hashes and calculates the hash256
"""
function merkle_parent(hash1::Vector{UInt8}, hash2::Vector{UInt8})
    hash256(append!(copy(hash1), hash2))
end

"""
    merkle_parent_level(Vector{Vector{UInt8}})
    -> Vector{Vector{UInt8}}

Takes a list of binary hashes and returns a list that"s half
the length
Returns an error if the list has exactly 1 element
"""
function merkle_parent_level(hashes::Vector{Vector{UInt8}})
    h = copy(hashes)
    if length(h) == 1
        error("Cannot take a parent level with only 1 item")
    end
    if length(h) % 2 == 1
        push!(h, h[end])
    end
    parent_level = Vector{UInt8}[]
    for i in 1:2:length(h)
        parent = merkle_parent(h[i], h[i + 1])
        push!(parent_level, parent)
    end
    parent_level
end

"""
    merkle_root(hashes::Vector{Vector{UInt8}})
    -> Vector{UInt8}

Takes a list of binary hashes and returns the merkle root
"""
function merkle_root(hashes::Vector{Vector{UInt8}})
    current_level = hashes
    while length(current_level) > 1
        current_level = merkle_parent_level(current_level)
    end
    current_level[1]
end

mutable struct MerkleTree
    total::Integer
    max_depth::Integer
    nodes
    current_depth::Integer
    current_index::Integer
    MerkleTree(total, max_depth, nodes, current_depth::Integer=1, current_index::Integer=1) = new(total, max_depth, nodes, current_depth, current_index)
end

"""
    MerkleTree(x::Integer) -> MerkleTree

Create a MerkleTree of `x` elements
"""
function MerkleTree(total::Integer)
    max_depth = 1 + ceil(Integer, log2(total))
    nodes = Vector{Union{Nothing, Vector{UInt8}}}[]
    for depth in 1:max_depth
        num_items = ceil(Integer, total / 2.0^(max_depth - depth))
        level_hashes = fill(nothing, num_items)
        push!(nodes, level_hashes)
    end
    MerkleTree(total, max_depth, nodes)
end

function show(io::IO, z::MerkleTree)
    result, depth = "", 1
    for level in z.nodes
        items, index = "", 1
        for h in level
            if h == nothing
                short = "  Nothing   "
            else
                short = " " * bytes2hex(h)[1:8] * "..."
            end
            if depth == z.current_depth && index == z.current_index
                result *= (items * "*" * short[2:end-2] * " *")
            else
                result *= (items * short)
            end
            index +=1
        end
        result *= (items * "\n")
        depth += 1
    end
    print(io, result)
end

"""
Reduce depth by 1 and halve the index
"""
function up!(tree::MerkleTree)
    if tree.current_depth != 1
        tree.current_depth -= 1
        tree.current_index = div(tree.current_index + 1, 2)
    end
end

"""
Increase depth by 1 and double the index
"""
function left!(tree::MerkleTree)
    tree.current_depth += 1
    tree.current_index = tree.current_index * 2 - 1
end

"""
Increase depth by 1 and double the index - 1
"""
function right!(tree::MerkleTree)
    tree.current_depth += 1
    tree.current_index *= 2;
end

"""
    root(tree::MerkleTree)
    -> Union{Nothing, Vector{UInt8}}

Returns the value of the tree root
"""
function root(tree::MerkleTree)
    tree.nodes[1][1]
end

"""
    set_current_node!(tree::MerkleTree, value::Vector{UInt8})
    -> Nothing

Set `value` to current node
"""
function set_current_node!(tree::MerkleTree, value::Vector{UInt8})
    tree.nodes[tree.current_depth][tree.current_index] = value
end

"""
    get_left_node(tree::MerkleTree)
    -> Union{Nothing, Vector{UInt8}}

Returns the value of the current node
"""
function get_current_node(tree::MerkleTree)
    tree.nodes[tree.current_depth][tree.current_index]
end

"""
    get_left_node(tree::MerkleTree)
    -> Union{Nothing, Vector{UInt8}}

Returns the value of the left children
"""
function get_left_node(tree::MerkleTree)
    tree.nodes[tree.current_depth + 1][tree.current_index * 2 - 1]
end

"""
    get_right_node(tree::MerkleTree)
    -> Union{Nothing, Vector{UInt8}}

Returns the value of the right children
"""
function get_right_node(tree::MerkleTree)
    tree.nodes[tree.current_depth + 1][tree.current_index * 2]
end

"""
    is_leaf(tree::MerkleTree) -> Bool

Checks current node is a leaf
"""
function is_leaf(tree::MerkleTree)
    tree.current_depth == tree.max_depth
end

"""
    right_exists(tree::MerkleTree) -> Bool

Check if a right children exists
"""
function right_exists(tree::MerkleTree)
    try
        length(tree.nodes[tree.current_depth + 1]) >= tree.current_index * 2
    catch
        false
    end
end

"""
    populate!(tree::MerkleTree, flag_bits::Vector{Bool},
              hashes::Vector{Vector{UInt8}})
    -> Nothing

Fills a Merkle Tree with a list of hashes
"""
function populate!(tree::MerkleTree, flag_bits::Vector{Bool}, hashes::Vector{Vector{UInt8}})
    while root(tree) == nothing
        if is_leaf(tree)
            popfirst!(flag_bits)
            set_current_node!(tree, popfirst!(hashes))
            up!(tree)
        else
            left_hash = get_left_node(tree)
            if left_hash == nothing
                if !(popfirst!(flag_bits))
                    set_current_node!(tree, popfirst!(hashes))
                    up!(tree)
                else
                    left!(tree)
                end
            elseif right_exists(tree)
                right_hash = get_right_node(tree)
                if right_hash == nothing
                    right!(tree)
                else
                    set_current_node!(tree, merkle_parent(left_hash, right_hash))
                    up!(tree)
                end
            else
                set_current_node!(tree, merkle_parent(left_hash, left_hash))
                up!(tree)
            end
        end
    end
    if length(hashes) != 0
        error("hashes not all consumed ", length(hashes))
    end
    for flag_bit in flag_bits
        if flag_bit
            error("flag bits not all consumed")
        end
    end
end

end # module


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
cli.jl defines the `jlcodestyle.cli` module containing functions for the `jlcodestyle` CLI.

Notes
-----
* CLI functions are defined in a .jl file so that testing and code quality tools can by
  applied to the CLI source code.
"""
module cli

# --- Exports

export parse_args, run

# --- Imports

# External packages
using ArgParse: ArgParse
using JuliaFormatter

# Local modules
using ..jlcodestyle

# --- Functions/Methods

"""
    parse_args(; raw_args::Vector{<:AbstractString}=ARGS)::Dict

Parse and return CLI arguments contained in `raw_args`. By default, `raw_args` is set to
`ARGS`, the command-line arguments provided to the executable that called `parse_args()`.
"""
function parse_args(; raw_args::Vector{<:AbstractString}=ARGS)::Dict

    # Define command-line arguments
    description = "Check source code files against Julia style conventions."
    arg_table = ArgParse.ArgParseSettings(; prog="jlcodestyle", description=description)
    ArgParse.@add_arg_table! arg_table begin
        "--overwrite", "-o"
        help = "overwrite files with reformatted source code"
        action = :store_true

        "--style", "-s"
        help = "Julia style convention to apply. Supported styles: blue, yas, default"
        default = "blue"

        "--verbose", "-v"
        help = "display more output to the console"
        action = :store_true

        "--version", "-V"
        help = "show version and exit"
        action = :store_true

        "paths"
        help = "space-separated list of files or directories to apply style conventions to"
        nargs = '*'
    end

    # Parse command-line arguments
    args::Dict = ArgParse.parse_args(raw_args, arg_table)
    args["paths"] = convert(Vector{String}, args["paths"])

    # Set code style
    style_str = lowercase(args["style"])
    if style_str == "default" || style_str == "defaultstyle"
        args["style"] = DefaultStyle()
    elseif style_str == "yas" || style_str == "yasstyle"
        args["style"] = YASStyle()
    else
        if !(style_str == "blue" || style_str == "bluestyle")
            style = args["style"]
            message = "Invalid style: $(style). Using BlueStyle."
            @warn message
        end

        args["style"] = BlueStyle()
    end

    return args
end

"""
    run(paths::Vector; <keyword arguments>)

Run code style checks for files contained in `paths`.

# Keyword Arguments

* `style::JuliaFormatter.AbstractStyle`: code style to apply.
  Default: `BlueStyle()`

* `overwrite::Bool`: overwrite existing files with style-corrected versions.
  Default: `false`

* `verbose::Bool`: print more output to the console. Default: `false`
"""
function run(
    paths::Vector;
    style::JuliaFormatter.AbstractStyle=BlueStyle(),
    overwrite::Bool=false,
    verbose::Bool=false,
)
    # --- Check arguments

    # Ensure that `paths` contains strings
    paths = convert(Vector{String}, paths)

    # --- Preparations

    # Set paths
    if isempty(paths)
        paths = ["."]
    end

    # --- Emit messages

    if verbose
        style_name = split("$(typeof(style))", '.')[2]
        @info "Style = $(style_name)"

        @info "Overwrite = $(overwrite)"
    end

    # --- Check code style

    check_passed = format(paths; style=style, overwrite=overwrite, verbose=verbose)

    if check_passed
        if verbose
            println()
            println("No style errors found.")
        end
    else
        if verbose
            println(stderr)
        end
        if overwrite
            println(stderr, "Style errors found. Files modified to correct errors.")
        else
            println(stderr, "Style errors found. Files not modified.")
        end
    end

    return nothing
end

end  # End of jlcodestyle.cli module


text:
using JuMP, Base.Test, AmplNLWriter

## Solve test problem with sind and cosd functions
 #
 #  min   (7 - (3*cosd(x1) + 5*cosd(x2)))^2 + (0 - (3*sind(x1) + 5*sind(x2)))^2
 #  s.t.  x1, x2 continuous
 #
 #  The optimal objective value is 0
 ##

# solver = AmplNLSolver(Ipopt.amplexe, ["print_level=0"])

@testset "example: jump_nltrig" begin
    m = Model(solver=solver)
    @variable(m, x[1:2])

    @NLobjective(m, Min, (7 - (3*cosd(x[1]) + 5*cosd(x[2])))^2 + (0 - (3*sind(x[1]) + 5*sind(x[2])))^2)

    setvalue(x[1], 30)
    setvalue(x[2], -50)
    @test solve(m) == :Optimal
    @test isapprox(getvalue(x)[:], [38.21321, -21.78678], atol=1e-5)
    @test isapprox(getobjectivevalue(m), 0.0, atol=1e-5)
    # Now try from the other side
    setvalue(x[1], -30)
    setvalue(x[2], 50)
    @test solve(m) == :Optimal
    @test isapprox(getvalue(x)[:], [-38.21321, 21.78678], atol=1e-5)
    @test isapprox(getobjectivevalue(m), 0.0, atol=1e-5)
end