text:
"""

    SGS(E, b, g)

The stochastically generated skewed (SGS) distribution with parameters E, b and g, derived by Sardeshmukh and Sura (2009) and following the form of Sardeshmukh et al. (2015). The SGS probability density function (pdf) is

```math
p(x) = \\frac{1}{\\mathcal{N}} \\left[ \\left( Ex + g \\right)^{2} + b^{2}  \\right]^{-\\left(1 + \\left(1/E^{2}\\right)\\right)} \\mathrm{exp} \\! \\left[ \\frac{2g}{E^{2}b} \\mathrm{arctan} \\! \\left(\\frac{Ex + g}{b} \\right) \\right]
```

where the normalization constant, ``\\mathcal{N}``, is given as

```math
\\mathcal{N} = \\frac{2 \\pi \\nu^{1/2} \\left( 2b \\right)^{-\\left(2 \\nu + 1 \\right)} \\Gamma \\! \\left(2\\nu + 1\\right)}{\\Gamma \\! \\left( \\nu + 1 - iq/2 \\right) \\Gamma \\! \\left( \\nu + 1 + iq/2 \\right)}
```

See also:

>Sardeshmukh, P. D., and P. Sura, 2009: Reconciling Non-Gaussian Climate Statistics with Linear Dynamics. *Journal of Climate*, **22**, 1193–1207, [https://doi.org/10.1175/2008JCLI2358.1](https://doi.org/10.1175/2008JCLI2358.1).
>
>Sardeshmukh, P. D., G. P. Compo, and C. Penland, 2015: Need for Caution in Interpreting Extreme Weather Statistics. *Journal of Climate*, **28**, 9166–9187, [https://doi.org/10.1175/JCLI-D-15-0020.1](https://doi.org/10.1175/JCLI-D-15-0020.1).

# Examples
```jldoctest
julia> d = SGS(0.6236095644623235, 1.1709106481844573, 0.48997894350611143)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

```

"""
struct SGS{T<:Real} <: ContinuousUnivariateDistribution
    E::T
    b::T
    g::T

    function SGS{T}(E::T, b::T, g::T) where {T}
        # This constraint check should be required if estimating the SGS
        # distribution using the method of moments.
        # @check_args(SGS, E ≥ zero(E) && b > zero(b))
        new(E, b, g)
    end
end

# Outer constructors.
SGS(E::T, b::T, g::T) where {T<:Real} = SGS{T}(E, b, g)
# NOTE: dots after promote()... expand the tuple.
SGS(E::Real, b::Real, g::Real) = SGS(promote(E, b, g)...)
SGS(E::Integer, b::Integer, g::Integer) = SGS(Float64(E), Float64(b), Float64(g))
SGS() = SGS(1.0, 1.0, 1.0)

# Type conversions.
# TODO: More type checking here.

# Still can't seem to find much on the @distr_support macro.
@distr_support SGS -Inf Inf

# Parameter retrieval
# See https://juliastats.github.io/Distributions.jl/stable/univariate/#Parameter-Retrieval-1

"""

    params(d::SGS)

Return a tuple of SGS parameters `E`, `b` and `g`. Let `d` be a distribution of type `D`, then `D(params(d)...)` will construct exactly the same distribution as ``d``.

# Examples
```jldoctest
julia> d = SGS(0.6236095644623235, 1.1709106481844573, 0.48997894350611143)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> params(d)
(0.6236095644623235, 1.1709106481844573, 0.48997894350611143)

```

"""
params(d::SGS) = (d.E, d.b, d.g)

# Required methods for a UnivariateDistribution type
# See https://juliastats.github.io/Distributions.jl/latest/extends.html#Univariate-Distribution-1
# rand(d::Power)
# sampler(d::Power)
# pdf(d::Power, x::Real)
# logpdf(d::Power, x::Real)
# cdf(d::Power, x::Real)
# quantile(d::Power, q::Real)
# minimum(d::Power)
# maximum(d::Power)
# insupport(d::Power, x::Real)

# Sampling

"""

    invitp(d::SGS)

Use a linear interpolation to produce an inversely transformed sample of the
SGS distribution.

# External links

* [Inverse transform sampling on Wikipedia](https://en.wikipedia.org/wiki/Inverse_transform_sampling)

"""
function invitp(d::SGS)
    xarr = collect(-10:0.01:10)
    A = cdf.(d, xarr)
    itp = LinearInterpolation(A, xarr, extrapolation_bc=Line())
    return itp
end

"""

    rand(d::SGS)

Generate a scalar sample from the SGS distribution `d`.

"""
function rand(d::SGS)
    return rand(GLOBAL_RNG, d)
end

"""

    rand(d::SGS, n::Int64)

Generate an array of `n` samples from the SGS distribution `d`.

"""
function rand(d::SGS, n::Int64)
    return rand(GLOBAL_RNG, d, n)
end

"""

    rand(d::SGS)

Generate a scalar sample from the SGS distribution `d`.

"""
function rand(rng::AbstractRNG, d::SGS)
    itp = invitp(d)
    return itp(rand(rng))
end

"""

    rand(rng::AbstractRNG, d::SGS, n::Int64)

Generate an array of `n` samples from the SGS distribution `d`.

"""
function rand(rng::AbstractRNG, d::SGS, n::Int64)
    itp = invitp(d)
    return itp(rand(rng, n))
end

# Statistics

"""

    mean(d::SGS)

Compute the mean of the SGS distribution.

# Examples
```jldoctest
julia> d = SGS(0.6236095644623235, 1.1709106481844573, 0.48997894350611143)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> mean(d)
0.0

```

"""
function mean(d::SGS)
    return 0.0
end

"""

    var(d::SGS)

Compute the variance of the SGS distribution.

# Examples
```jldoctest
julia> d = SGS(0.6236095644623235, 1.1709106481844573, 0.48997894350611143)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> var(d)
1.0

```

"""
function var(d::SGS)
    (E, b, g) = params(d)
    return (g^2 + b^2) / (2 - E^2)
end

"""

    std(d::SGS)

Compute the standard deviation of the SGS distribution.

# Examples
```jldoctest
julia> d = SGS(0.6236095644623235, 1.1709106481844573, 0.48997894350611143)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> std(d)
1.0

```

"""
function std(d::SGS)
    return sqrt(var(d))
end

"""

    skewness(d::SGS)

Compute the skewness of the SGS distribution.

# Examples
```jldoctest
julia> d = SGS(0.6236095644623235, 1.1709106481844573, 0.48997894350611143)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> skewness(d)
1.0

```

"""
function skewness(d::SGS)
    (E, b, g) = params(d)
    return moment_constraint(3, E) ? ((2.0 * E) / (1.0 - E^2)) * (g / sqrt(var(d))) : NaN
end

"""

    kurtosis(d::SGS; correction::Bool=true)

Compute the excess kurtosis of the SGS distribution. Excess kurtosis is returned by default with `correction=true`.

# Examples
```jldoctest
julia> d = SGS(0.6236095644623235, 1.1709106481844573, 0.48997894350611143)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> kurtosis(d, true)
4.999999999999998

```

"""
function kurtosis(d::SGS; excess::Bool=true)
    (E, b, g) = params(d)
    if moment_constraint(4, E)
        term1 = 1.5 * ((1.0 - E^2) / (1.0 - 1.5 * E^2)) * skewness(d)^2
        term2 = (3.0 * E^2) / (1.0 - 1.5*E^2)
        kurt = term1 + term2
        if excess
            return kurt         # excess kurtosis
        else
            return kurt + 3     # kurtosis
        end
    else
        return NaN  # kurtosis does not exist
    end
    # Skewness-kurtosis constraint.
    # if excess_kurtosis ≤ 1.5 * skewness(d)^2
    #     return NaN
    # else
    #     return excess_kurtosis
end

# Evaluation

"""

    pdf(d::SGS, x::Real; norm::Bool=true)

Compute the pdf of the SGS distribution. The SGS pdf is normalized by default,
i.e., `norm=true`. If `norm` is set to `false`, then the SGS normalization
constant will not be applied.

Here, the SGS probability density function (pdf) is given in the form of [Sardeshmukh et al. (2015)](https://doi.org/10.1175/JCLI-D-15-0020.1) as

```math
p(x) = \\frac{1}{\\mathcal{N}} \\left[ \\left( Ex + g \\right)^{2} + b^{2}  \\right]^{-\\left(1 + \\left(1/E^{2}\\right)\\right)} \\mathrm{exp} \\! \\left[ \\frac{2g}{E^{2}b} \\mathrm{arctan} \\! \\left(\\frac{Ex + g}{b} \\right) \\right]
```

where the normalization constant, ``\\mathcal{N}``, is

```math
\\mathcal{N} = \\frac{2 \\pi \\nu^{1/2} \\left( 2b \\right)^{-\\left(2 \\nu + 1 \\right)} \\Gamma \\! \\left(2\\nu + 1\\right)}{\\Gamma \\! \\left( \\nu + 1 - iq/2 \\right) \\Gamma \\! \\left( \\nu + 1 + iq/2 \\right)}
```

See also:

>Sardeshmukh, P. D., G. P. Compo, and C. Penland, 2015: Need for Caution in Interpreting Extreme Weather Statistics. *Journal of Climate*, **28**, 9166–9187, [https://doi.org/10.1175/JCLI-D-15-0020.1](https://doi.org/10.1175/JCLI-D-15-0020.1).

# Examples
```jldoctest
julia> d = SGS(0.6236095644623235, 1.1709106481844573, 0.48997894350611143)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> pdf(d, 0.0)
0.46210984589674214

```

"""
function pdf(d::SGS, x::Real; norm::Bool=true)
    return _pdf(d, x, norm)
end

"""

    pdf(d::SGS, x::AbstractArray; norm::Bool=true)

Compute the pdf of the SGS distribution. The SGS pdf is normalized by default,
i.e., `norm=true`. If `norm` is set to `false`, then the SGS normalization
constant will not be applied.

Here, the SGS probability density function (pdf) is given in the form of [Sardeshmukh et al. (2015)](https://doi.org/10.1175/JCLI-D-15-0020.1) as

```math
p(x) = \\frac{1}{\\mathcal{N}} \\left[ \\left( Ex + g \\right)^{2} + b^{2}  \\right]^{-\\left(1 + \\left(1/E^{2}\\right)\\right)} \\mathrm{exp} \\! \\left[ \\frac{2g}{E^{2}b} \\mathrm{arctan} \\! \\left(\\frac{Ex + g}{b} \\right) \\right]
```

where the normalization constant, ``\\mathcal{N}``, is

```math
\\mathcal{N} = \\frac{2 \\pi \\nu^{1/2} \\left( 2b \\right)^{-\\left(2 \\nu + 1 \\right)} \\Gamma \\! \\left(2\\nu + 1\\right)}{\\Gamma \\! \\left( \\nu + 1 - iq/2 \\right) \\Gamma \\! \\left( \\nu + 1 + iq/2 \\right)}
```

See also:

>Sardeshmukh, P. D., G. P. Compo, and C. Penland, 2015: Need for Caution in Interpreting Extreme Weather Statistics. *Journal of Climate*, **28**, 9166–9187, [https://doi.org/10.1175/JCLI-D-15-0020.1](https://doi.org/10.1175/JCLI-D-15-0020.1).

# Examples
```jldoctest
julia> d = SGS(0.6236095644623235, 1.1709106481844573, 0.48997894350611143)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> x = [-0.2, 0.4, 1.5, 5.0]
4-element Array{Float64,1}:
 -0.2
  0.4
  1.5
  5.0

julia> pdf(d, x)
4-element Array{Float64,1}:
 0.4821612875311956
 0.3552310712921366
 0.091209594116695
 0.0011832665173870727

```

"""
function pdf(d::SGS, x::AbstractArray; norm::Bool=true)
    pdfs = Array{Float64}(undef, length(x))
    for (indx, num) in enumerate(x)
        pdfs[indx] = _pdf(d, num, norm)
    end
    return pdfs
end

function pdf(E::Real, b::Real, g::Real, x::Real)
    return _pdf(SGS(E, b, g), x, false)
end

function _pdf(d::SGS, x::Real, norm::Bool)
    (E, b, g) = params(d)

    @debug "E = $E"
    @debug "b = $b"
    @debug "g = $g"

    # if E < 0.05
    #     return _pdf_E_lim(d, x)
    # end

    # Calculate the log of the normalization constant.
    if norm
        log_N = logN(E, b, g)
    else
        log_N = 0.0
    end

    # Calculate the log pdf.
    q = (2.0 * g) / (E^2 * b)
    @debug "p(x) q = $q"
    atan_term = atan((E * x + g) / b)
    @debug "p(x) atan_term = $atan_term"
    log_px = -log_N + q * atan_term - (1.0 + 1.0/E^2) * log((E * x + g)^2 + b^2)
    @debug "log(p(x)) = $log_px"

    # Calculate the pdf.
    px = exp(log_px)
    @debug "p(x) = $px"

    return px

end

function _pdf(E::Real, b::Real, g::Real, x::Real)

    @debug "E = $E"
    @debug "b = $b"
    @debug "g = $g"

    # Calculate the log of the normalization constant.
    # if norm
    #     log_N = logN(E, b, g)
    # else
    #     log_N = 0.0
    # end

    log_N = 0.0

    # Calculate the log pdf.
    q = (2.0 * g) / (E^2 * b)
    @debug "p(x) q = $q"
    atan_term = atan((E * x + g) / b)
    @debug "p(x) atan_term = $atan_term"
    log_px = -log_N + q * atan_term - (1.0 + 1.0/E^2) * log((E * x + g)^2 + b^2)
    @debug "log(p(x)) = $log_px"

    # Calculate the pdf.
    px = exp(log_px)
    @debug "p(x) = $px"

    return px

end

function _pdf_E_lim(d::SGS, x::Real)
    (E, b, g) = params(d)
    px = 1.0 / sqrt(2*π) * exp((-x^2 / 2.0) + ((E * g * x) / 3.0) * (x^2 - 3.0))
    return px
end

"""

    pdf_xmax(d::SGS)

Compute the location of the unique pdf maximum.

# Examples
```jldoctest
julia> d = SGS(0.6236095644623235, 1.1709106481844573, 0.48997894350611143)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> pdf_xmax(d)
-0.22000000000000003

```

"""
function pdf_xmax(d::SGS)
    (E, b, g) = params(d)
    return -1.0 * E * g / (1.0 + E^2)
end

"""

    logN(E::Real, b::Real, g::Real)

Compute the log of the normalization constant `N` of the SGS distribution.

# Examples
```jldoctest
julia> d = SGS(0.6236095644623235, 1.1709106481844573, 0.48997894350611143)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> logN(d.E, d.b, d.g)
-0.07844022339848955

```

"""
function logN(E::Real, b::Real, g::Real)

    # Calculate parameters for the normalization constant.
    ν = 1.0 / E^2
    q = (2.0 * g * ν) / b
    @debug "ν = $ν"
    @debug "q = $q"

    # Calculate first two terms in numerator, taking logs for numerical stability.
    term1 = log(2.0 * π * sqrt(ν))
    term2 = log((2.0 * b)^(-(2.0 * ν + 1.0)))
    @debug "log(2.0 * π * sqrt(ν)) = $term1"
    @debug "log((2.0 * b)^(1.0 - 2.0 * ν)) = $term2"

    # Calculate parameters for the gamma function.
    logΓ0 = loggamma(2.0 * ν + 1.0)
    logΓ1 = loggamma(ν + 1.0 - (1im * q / 2.0))
    logΓ2 = loggamma(ν + 1.0 + (1im * q / 2.0))
    @debug "log(Γ0) = $logΓ0"
    @debug "log(Γ1) = $logΓ1"
    @debug "log(Γ2) = $logΓ2"

    # Find the log of the normalization constant.
    logN = real(term1 + term2 + logΓ0 - logΓ1 - logΓ2)
    @debug "log(N) = $logN"
    # N = exp(logN)
    # @debug "N = $N"
    #
    # @debug "Returning logN..."
    return logN

end

"""

    logpdf(d::SGS, x::AbstractArray)

Calculate the log of the probability density function of the SGS distribution.

# Examples
```jldoctest
julia> d = fit_mm(SGS, 1, 1, 5)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> x = [-10.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0, 5.0]
8-element Array{Float64,1}:
 -10.0
  -1.0
  -0.5
   0.0
   0.5
   1.0
   2.0
   5.0

julia> logpdf(d, x)
8-element Array{Float64,1}:
 -15.504184293152903
  -1.3393373142985014
  -0.8053427806087643
  -0.7719526544799672
  -1.1296521694192523
  -1.7130802377432286
  -3.099847429101194
  -6.739476429937023

```

"""
function logpdf(d::SGS, x::AbstractArray)
    lpdfs = Array{Float64}(undef, length(x))
    for (indx, num) in enumerate(x)
        lpdfs[indx] = log(pdf(d, num))
    end
    return lpdfs
end

"""

    logpdf(d::SGS, x::Real)

Calculate the log of the cumulative distribution function of the SGS
distribution.

# Examples
```jldoctest
julia> d = fit_mm(SGS, 1, 1, 5)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> x = 0.5
0.5

julia> logpdf(d, x)
-1.1296521694192523

```

"""
function logpdf(d::SGS, x::Real)
    return log(pdf(d, x))
end

"""

    cdf(d::SGS, x::AbstractArray)

Calculate the cumulative distribution function of the SGS distribution.

# Examples
```jldoctest
julia> d = fit_mm(SGS, 1, 1, 5)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> x = [-10.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0, 5.0]
8-element Array{Float64,1}:
 -10.0
  -1.0
  -0.5
   0.0
   0.5
   1.0
   2.0
   5.0

julia> cdf(d, x)
8-element Array{Float64,1}:
 2.7059681332492457e-7
 0.1268532553530455
 0.3074639317904149
 0.5437968749444915
 0.7430905390656029
 0.86706141301943
 0.9654409769867174
 0.9986564347567426

```

"""
function cdf(d::SGS, x::AbstractArray)
    cdfs = Array{Float64}(undef, length(x))
    for (indx, num) in enumerate(x)
        cdfs[indx] = cdf(d, num)
    end
    return cdfs
end

"""

    cdf(d::SGS, x::Real)

Calculate the cumulative distribution function of the SGS distribution.

# Examples
```jldoctest
julia> d = fit_mm(SGS, 1, 1, 5)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> x = 0.5
0.5

julia> cdf(d, x)
0.7430905390656029

```

"""
function cdf(d::SGS, x::Real)
    integral, error = quadgk(x -> pdf(d, x), -Inf, x)
    return integral
end

# function cdf(d::SGS, x::Real)
#
#     dx = 0.01
#     σval = collect(-10:dx:x)
#
#     cdf = 0.0
#     for i in σval
#         cdf += pdf(d, i) * dx
#     end
#
#     return cdf
#
# end
#
# function cdf(d::SGS, x::AbstractArray)
#
#     dx = 0.01
#     cdf = zeros(length(x))
#     for i=1:length(x)
#         for xval in -10:dx:x[i]
#             cdf[i] += pdf(d, xval) * dx
#         end
#     end
#     return cdf
#
# end

"""

    logcdf(d::SGS, x::AbstractArray)

Calculate the log of the cumulative distribution function of the SGS
distribution.

# Examples
```jldoctest
julia> d = fit_mm(SGS, 1, 1, 5)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> x = [-10.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0, 5.0]
8-element Array{Float64,1}:
 -10.0
  -1.0
  -0.5
   0.0
   0.5
   1.0
   2.0
   5.0

julia> logcdf(d, x)
8-element Array{Float64,1}:
 -15.12263589760972
  -2.064724330254346
  -1.1793974936056804
  -0.6091794935003683
  -0.2969373856106992
  -0.14264547077775455
  -0.035170311051895865
  -0.0013444686363078492

```

"""
function logcdf(d::SGS, x::AbstractArray)
    lcdfs = Array{Float64}(undef, length(x))
    for (indx, num) in enumerate(x)
        lcdfs[indx] = log(cdf(d, num))
    end
    return lcdfs
end

"""

    logcdf(d::SGS, x::Real)

Calculate the log of the cumulative distribution function of the SGS
distribution.

# Examples
```jldoctest
julia> d = fit_mm(SGS, 1, 1, 5)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> x = 0.5
0.5

julia> logcdf(d, x)
-0.2969373856106992

```

"""
function logcdf(d::SGS, x::Real)
    integral, error = quadgk(x -> pdf(d, x), -Inf, x)
    return log(integral)
end

# function logcdf(d::SGS, x::Real)
#     return log(cdf(d, x))
# end

"""

    ccdf(d::SGS, x::AbstractArray)

Calculate the complementary cumulative distribution function of the SGS
distribution.

# Examples
```jldoctest
julia> d = fit_mm(SGS, 1, 1, 5)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> x = [-10.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0, 5.0]
8-element Array{Float64,1}:
 -10.0
  -1.0
  -0.5
   0.0
   0.5
   1.0
   2.0
   5.0

julia> ccdf(d, x)
8-element Array{Float64,1}:
 0.9999997294031867
 0.8731467446469545
 0.6925360682095851
 0.45620312505550853
 0.25690946093439715
 0.13293858698057004
 0.03455902301328262
 0.0013435652432574052

```

"""
function ccdf(d::SGS, x::AbstractArray)
    ccdfs = Array{Float64}(undef, length(x))
    for (indx, num) in enumerate(x)
        ccdfs[indx] = 1.0 - cdf(d, num)
    end
    return ccdfs
end

# function ccdf(d::SGS, x::AbstractArray)
#     ccdfs = Array{Float64}(undef, length(x))
#     for (indx, num) in enumerate(x)
#         ccdfs[indx] = ccdf(d, num)
#     end
#     return ccdfs
# end

"""

    ccdf(d::SGS, x::Real)

Calculate the complementary cumulative distribution function of the SGS
distribution.

# Examples
```jldoctest
julia> d = fit_mm(SGS, 1, 1, 5)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> x = 0.5
0.5

julia> ccdf(d, x)
0.25690946093439715

```

"""
function ccdf(d::SGS, x::Real)
    integral, error = quadgk(x -> pdf(d, x), -Inf, x)
    return 1.0 - integral
end

# function ccdf(d::SGS, x::Real)
#     return 1.0 - cdf(d, x)
# end

"""

    logccdf(d::SGS, x::AbstractArray)

Calculate the log of the complementary cumulative distribution function of the
SGS distribution.

# Examples
```jldoctest
julia> d = fit_mm(SGS, 1, 1, 5)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> x = [-10.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0, 5.0]
8-element Array{Float64,1}:
 -10.0
  -1.0
  -0.5
   0.0
   0.5
   1.0
   2.0
   5.0

julia> logccdf(d, x)
8-element Array{Float64,1}:
 -2.70596849901216e-7
 -0.1356516448893756
 -0.3673949582198083
 -0.7848171189678753
 -1.3590315482404094
 -2.017868009426468
 -3.365086604737457
 -6.612428568931222

```

"""
function logccdf(d::SGS, x::AbstractArray)
    ccdfs = Array{Float64}(undef, length(x))
    for (indx, num) in enumerate(x)
        ccdfs[indx] = log(1.0 - cdf(d, num))
    end
    return ccdfs
end

# function logccdf(d::SGS, x::AbstractArray)
#     logccdfs = Array{Float64}(undef, length(x))
#     for (indx, num) in enumerate(x)
#         logccdfs[indx] = log(ccdf(d, num))
#     end
#     return logccdfs
# end

"""

    logccdf(d::SGS, x::Real)

Calculate the log of the complementary cumulative distribution function of the
SGS distribution.

# Examples
```jldoctest
julia> d = fit_mm(SGS, 1, 1, 5)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> x = 0.5
0.5

julia> logccdf(d, x)
-1.3590315482404094

```

"""
function logccdf(d::SGS, x::Real)
    integral, error = quadgk(x -> pdf(d, x), -Inf, x)
    return log(1.0 - integral)
end

# function logccdf(d::SGS, x::Real)
#     return log(ccdf(d, x))
# end


text:
Base.firstindex(::DailyDatesRange) = 1
Base.isempty(::DailyDatesRange) = false # has at least one item
Base.length(d::DailyDatesRange) = daycount(d.daycountconvention, d.startdate, d.enddate) + 1
Base.issorted(d::DailyDatesRange) = isascending(d) || length(d) == 1
Base.minimum(dd::DailyDatesRange) = dd.startdate
Base.maximum(dd::DailyDatesRange) = dd.enddate

function Base.in(date::Date, dd::DailyDatesRange{A, T}) :: Bool where {A,T<:Union{Actual360, Actual365}}
    return dd.startdate <= date <= dd.enddate
end

function Base.in(date::Date, dd::DailyDatesRange{A, BDays252}) :: Bool where {A}
    return (dd.startdate <= date <= dd.enddate) && BusinessDays.isbday(dd.daycountconvention.hc, date)
end

yearfraction(d::DailyDatesRange) = yearfraction(d.daycountconvention, 1)
yearfractionvalue(d::DailyDatesRange) = yearfractionvalue(yearfraction(d))

function Base.getindex(d::DailyDatesRange{true}, i::Integer) :: Dates.Date
    result = advancedays(d.daycountconvention, d.startdate, i - 1)
    @assert result <= d.enddate "Index $i out of bounds."
    return result
end

function Base.getindex(d::DailyDatesRange{false}, i::Integer) :: Dates.Date
    result = advancedays(d.daycountconvention, d.enddate, -(i - 1))
    @assert result >= d.startdate
    return result
end

function isascending(::DailyDatesRange{A}) :: Bool where {A}
    A
end

Base.reverse(dd::DailyDatesRange) = DailyDatesRange(dd.startdate, dd.enddate, dd.daycountconvention, !isascending(dd))
Base.show(io::IO, dd::DailyDatesRange) = print(io, "DailyDatesRange($(dd.startdate), $(dd.enddate), $(dd.daycountconvention))")

Base.iterate(dd::DailyDatesRange{true}) = (dd.startdate, dd.startdate)
Base.iterate(dd::DailyDatesRange{false}) = (dd.enddate, dd.enddate)

function Base.iterate(dd::DailyDatesRange{true}, current_date::Date)
    next_date = advancedays(dd.daycountconvention, current_date, 1)
    if next_date > dd.enddate
        return nothing
    else
        return (next_date, next_date)
    end
end

function Base.iterate(dd::DailyDatesRange{false}, current_date::Date)
    next_date = advancedays(dd.daycountconvention, current_date, -1)
    if next_date < dd.startdate
        return nothing
    else
        return (next_date, next_date)
    end
end

ext: jl


text:
################################################################################
#
#  Small groups
#
################################################################################

using Pkg.Artifacts

#export number_of_small_groups, small_groups_limit
export small_group
#CF: as of now, nothing works here and the exports cause errors in Oscar

mutable struct SmallGroupDB
  path::String
  max_order::Int
  db::Vector{Vector{NamedTuple{(:id, :name, :gens, :rels, :nontrivrels,
                                :orderdis, :ordersubdis, :isabelian, :iscyclic,
                                :issolvable, :isnilpotent, :autorder,
                                :aut_gens, :nchars, :dims, :schur, :galrep,
                                :fields, :mod),
                              Tuple{Tuple{Int, Int}, String, Vector{Perm{Int}},
                                    Vector{Vector{Int}}, Vector{Vector{Int}},
                                    Vector{Tuple{Int,Int}},
                                    Vector{Tuple{Int,Int}},
                                    Bool, Bool, Bool, Bool, fmpz,
                                    Vector{Vector{Vector{Int}}}, Int,
                                    Vector{Int},Vector{Int}, Vector{Int},
                                    Vector{Vector{fmpq}},
                                    Vector{Vector{Vector{Vector{fmpq}}}}}}}}

end

function SmallGroupDB(path::String)
  db = Vector{NamedTuple{(:id, :name, :gens, :rels, :nontrivrels, :orderdis,
                          :ordersubdis, :isabelian, :iscyclic, :issolvable,
                          :isnilpotent, :autorder, :aut_gens, :nchars, :dims,
                          :schur, :galrep, :fields, :mod),
                         Tuple{Tuple{Int, Int}, String, Vector{Perm{Int}},
                               Vector{Vector{Int}}, Vector{Vector{Int}},
                               Vector{Tuple{Int,Int}},
                               Vector{Tuple{Int,Int}}, Bool, Bool, Bool,
                               Bool, fmpz, Vector{Vector{Vector{Int}}}, Int,
                               Vector{Int},Vector{Int}, Vector{Int},
                               Vector{Vector{fmpq}},
                               Vector{Vector{Vector{Vector{fmpq}}}}}}}[]
  z = eltype(eltype(db))[]
  open(path) do io
    while !eof(io)
    e = _parse_row(io)
    if isempty(z) || e.id[1] == z[1].id[1]
      push!(z, e)
    else
      push!(db, z)
      z = eltype(eltype(db))[]
      push!(z, e)
    end
  end
  push!(db, z)
  end
  max_order = length(db)
  return SmallGroupDB(path, max_order, db)
end

function SmallGroupDBLegacy(path::String)

  f = function(z)
    if z isa BigInt
      return fmpz(z)
    elseif z isa Rational{BigInt}
      return fmpq()
    elseif z isa Vector{Vector{Rational{BigInt}}}
      return Vector{fmpq}[ fmpq.(v) for v in z]
    elseif z isa Vector{Vector{Vector{Vector{Rational{BigInt}}}}}
      zz = Vector{Vector{Vector{Vector{fmpq}}}}(undef, length(z))
      for i in 1:length(z)
        zz[i] = Vector{Vector{Vector{fmpq}}}(undef, length(z[i]))
        for j in 1:length(z[i])
          zz[i][j] = Vector{Vector{fmpq}}(undef, length(z[i][j]))
          for k in 1:length(z[i][j])
            zz[i][j][k] = fmpq.(z[i][j][k])
          end
        end
      end
      return zz
    else
      return z
    end
  end

  dbnew = Vector{NamedTuple{(:id, :name, :gens, :rels, :nontrivrels, :orderdis,
                          :ordersubdis, :isabelian, :iscyclic, :issolvable,
                          :isnilpotent, :autorder, :aut_gens, :nchars, :dims,
                          :schur, :galrep, :fields, :mod),
                         Tuple{Tuple{Int, Int}, String, Vector{Perm{Int}},
                               Vector{Vector{Int}}, Vector{Vector{Int}},
                               Vector{Tuple{Int,Int}},
                               Vector{Tuple{Int,Int}}, Bool, Bool, Bool,
                               Bool, fmpz, Vector{Vector{Vector{Int}}}, Int,
                               Vector{Int},Vector{Int}, Vector{Int},
                               Vector{Vector{fmpq}},
                               Vector{Vector{Vector{Vector{fmpq}}}}}}}[]

  db = Hecke.eval(Meta.parse(Base.read(path, String)))
  for i in 1:length(db)
    push!(dbnew, eltype(dbnew)(undef, length(db[i])))
    for j in 1:length(db[i])
      v = db[i][j]
      vv = map(f, v)
      dbnew[i][j] = (id = (i, j), vv...)
    end
  end

  max_order = length(dbnew)
  return SmallGroupDB(path, max_order, dbnew)
end

function show(io::IO, L::SmallGroupDB)
  print(io, "Database of small groups (order limit = ", L.max_order, ")")
end

const legacy_default_small_group_db = [joinpath(pkgdir, "data/small_groups_extended"), joinpath(pkgdir, "data/small_groups_default")]

function small_group_database()
  st = artifact"SmallGroupDB"
  return SmallGroupDB(joinpath(st, "SmallGroupDB", "data"))
end

const _DefaultSmallGroupDB = Ref{Any}(nothing)

function DefaultSmallGroupDB()
  _DB = _DefaultSmallGroupDB[]
  if _DefaultSmallGroupDB[] === nothing
    DB = small_group_database()
    _DefaultSmallGroupDB[] = DB
    return DB::SmallGroupDB
  else
    return _DB::SmallGroupDB
  end
end

isfrom_db(G::GrpGen) = G.isfromdb

function small_group(i, j; DB = DefaultSmallGroupDB())
  data = DB.db[i][j]
  #i < 1 || i > 63 && error("Group order ($i) must be between 1 and $small_groups_limit")
  #j < 1 || j > number_of_small_groups(i) && error("Index ($j) must be between 1 and $(number_of_small_groups(i))")
  P = SymmetricGroup(i)
  Gens = [one(P)]
  for p in data.gens
    push!(Gens, P(p))
  end
  G, PtoG, GtoP = generic_group(closure(Gens, *, one(P)), *)

  # The small groups have generators, for which we know
  # generators for the relations.
  # This is used in the homomorphism computation.
  # To make this work, the generators must much exactly
  # with the generators in the database.
  G.gens = Vector{Int}(undef, length(data.gens))
  for k in 1:length(data.gens)
    p = data.gens[k]
    G.gens[k] = PtoG[P(p)][]
  end

  G.issolvable = data.issolvable
  G.isnilpotent = data.isnilpotent
  G.isfromdb = true
  G.small_group_id = (i, j)
  return G
end

################################################################################
#
#  Better parsing
#
################################################################################

function _parse_row(io::IO)
  b = Base.read(io, UInt8)
  # id
  b, id = _parse(Tuple{Int, Int}, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, name = _parse(String, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, gens = _parse(Vector{Perm{Int}}, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, rels = _parse(Vector{Vector{Int}}, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, nontrivrels = _parse(Vector{Vector{Int}}, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, orderdis = _parse(Vector{Tuple{Int, Int}}, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, ordersubdis = _parse(Vector{Tuple{Int, Int}}, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, isabelian = _parse(Bool, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, iscyclic = _parse(Bool, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, issolvable = _parse(Bool, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, isnilpotent = _parse(Bool, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, autorder = _parse(fmpz, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, aut_gens = _parse(Vector{Vector{Vector{Int}}}, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, nchars = _parse(Int, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, dims = _parse(Vector{Int}, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, schur = _parse(Vector{Int}, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, galrep = _parse(Vector{Int}, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, fields = _parse(Vector{Vector{fmpq}}, io, b)
  @assert b == UInt(',')
  b = Base.read(io, UInt8)
  b, mod = _parse(Vector{Vector{Vector{Vector{fmpq}}}}, io, b)
  return (id = id, name = name, gens = gens, rels = rels, nontrivrels = nontrivrels, orderdis = orderdis,
          ordersubdis = ordersubdis, isabelian = isabelian, iscyclic = iscyclic, issolvable = issolvable,
          isnilpotent = isnilpotent, autorder = autorder, aut_gens = aut_gens, nchars = nchars, dims = dims,
          schur = schur, galrep = galrep, fields = fields, mod = mod)
end


text:
#=
test_src:
- Author: Berk
- Date: 2020-06-16
This tests everything to do with the core of the OCTHaGOn code,
without *many* machine learning components
=#

##########################
# SOURCE TESTS #
##########################

""" Tests expression parsing. """
function test_expressions()
    model, x, y, z, a = test_model()
    expr = :((x, y, z) -> sum(x[i] for i=1:4) - y[1] * y[2] + z)
    simp_expr = :(x -> sum(5 .* x))
    f = functionify(expr)

    # Testing function evaluation
    res = Base.invokelatest(f, (ones(4), ones(2),5)...)
    @test res isa Float64
    res = Base.invokelatest(f, (x, y, z)...)
    @test res isa JuMP.GenericQuadExpr
    @test res == sum(x[1:4]) - y[1] * y[2] + z

    # Testing variable input parsing from expression
    @test vars_from_expr(expr, model) == [x,y,z]
    @test vars_from_expr(simp_expr, model) == [x]
    @test vars_from_expr(:((x, a) -> x[2] + a[2,2]),model) == [x, a]

    # Testing "flattening of expressions" for nonlinearization
    expr_vars = vars_from_expr(expr, model)
    @test OCTHaGOn.get_var_ranges(expr_vars) == [(1:5),(6:8),9]
    @test OCTHaGOn.zeroarray([(1:5),(6:8),9]) == [zeros(5), zeros(3), 0]
    flat_expr = :((x...) -> $(expr)([x[i] for i in $(OCTHaGOn.get_var_ranges(expr_vars))]...))
    fn = functionify(flat_expr)
    @test Base.invokelatest(fn, [1,2,3,4,1,5,-6,-7,7]...) == Base.invokelatest(f, ([1,2,3,4,1], [5,-6,-7], 7)...)
    @test Base.invokelatest(fn, flat(expr_vars)...) == res

    # Testing proper mapping for expressions
    flatvars = flat([y[2], z, x[1:4]])
    vars = vars_from_expr(expr, model)
    @test OCTHaGOn.get_varmap(vars, flatvars) == [(2,2), (3,0), (1,1), (1, 2), (1,3), (1,4)]
    @test OCTHaGOn.get_datamap(vars, flatvars) == [7, 9, 1, 2, 3, 4]

    # Testing gradientify
    grad = gradientify(expr, vars_from_expr(expr, model))
    @test all(grad(ones(9)) .≈ [1, 1, 1, 1, 0, -1, -1, 0, 1])
    other_grad = gradientify(simp_expr, vars_from_expr(simp_expr, model))
    @test all(other_grad(ones(5)) .≈ 5 .* ones(5))
    @test all(grad(ones(9)) .≈ [1, 1, 1, 1, 0, -1, -1, 0, 1]) # in case of world age problems
    con = @constraint(model, Base.invokelatest(f, (x,y,z)...) >= 1) # also on JuMP constraints
    gradfn = gradientify(con, expr_vars)
    @test all(gradfn(ones(9)) .≈ grad(ones(9)))
    @test_throws OCTHaGOnException gradientify(@constraint(model, [x[1] x[2];
                                                              x[3] x[4]] in PSDCone()), x[1:4])
    # Testing vars_from_constraint as well
    @test all([var in [x[1], x[2], x[3], x[4], y[1], y[2], z] for var in vars_from_constraint(con)])
    m = pool1(false)
    l_constrs, nl_constrs = classify_constraints(m)
    l_vars = vars_from_constraint.(l_constrs)
    @test length.(l_vars) == [4,4,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
    nl_vars = vars_from_constraint.(nl_constrs)
    @test length.(nl_vars) == [3,3,3,2]
    g_fns = gradientify.(nl_constrs, nl_vars)
    @test g_fns[1](ones(length(nl_vars[1]))) == [-1, 2, -8]
end

function test_variables()
    model, x, y, z, a  = test_model()

    # Testing getting variables
    varkeys = ["x[1]", x[1], :z, :x];
    vars = [x[1], x[1], z, x[:]];
    @test all(vars .==  fetch_variable(model, varkeys))
end

function test_bounds()
    model, x, y, z, a = test_model()

    # Bounds and fixing variables
    bounds = get_bounds(model);
    @test bounds[z] == [-30., Inf]
    bounds = Dict(:x => [-10,1], z => [-10, 10])
    bound!(model, bounds)
    @test get_bounds(model)[z] == [-10, 10]
    @test get_bounds(model)[x[3]] == [-5, 1]

    # Check infeasible bounds
    new_bound = x[4] => [-10,-6]
    @test_throws OCTHaGOnException OCTHaGOn.check_infeasible_bound(new_bound)
    @test_throws OCTHaGOnException bound!(model, new_bound)

    # Check unbounds
    @test all(collect(values(get_unbounds(model)))[i] == [-2., Inf] for i = 1:4)

    # Check that JuMP does not set Inf bounds
    m = JuMP.Model()
    @variable(m, x[1:2])
    @variable(m, y >= 5)
    bound!(m, Dict(var => [-Inf, Inf] for var in all_variables(m)))
    @test !JuMP.has_lower_bound(x[1]) && JuMP.has_lower_bound(y) && !JuMP.has_upper_bound(x[2])
    bound!(m, Dict(var => [-10, 10] for var in all_variables(m)))
    @test isnothing(get_unbounds(m))

    # Test bounded auxiliary variables
    model, x, y, z, a = test_model()
    # For classifiers
    aux_vars, cons = bounded_aux(x, @variable(model, binary=true))
    @test length(cons) == 2*length(x)
    # For regressors
    aux_vars, cons = bounded_aux(x[1:end-1], x[end], @variable(model, binary=true))
    @test length(cons) == 2*length(x)
end

function test_sets()
    sets = [MOI.GreaterThan(2), MOI.EqualTo(0), MOI.SecondOrderCone(3), MOI.GeometricMeanCone(2), MOI.SOS1([1,2,3])]
    @test get_constant.(sets) == [2, 0, nothing, nothing, nothing]

    model, x, y, z, a = test_model()
    restrict_to_set(x[3], [1,2,3])
    @test length(all_variables(model)) == 18
end

function test_linearize()
    model, x, y, z, a = test_model()
    # Linearization of objective
    linearize_objective!(model);
    @objective(model, Min, x[3]^2)
    linearize_objective!(model);
    @test JuMP.objective_function(model) isa JuMP.VariableRef

    @constraint(model, sum(x[4]^2 + x[5]^2) <= z)
    @constraint(model, sum(y[:]) >= -2)
    # Separation of model constraints
    l_constrs, nl_constrs = classify_constraints(model)
    @test length(l_constrs) == 27 && length(nl_constrs) == 1
end

function test_nonlinearize(gm::GlobalModel = minlp(true))
    nonlinearize!(gm)
    set_optimizer(gm, SOLVER_SILENT)
    @test_throws ErrorException("The solver does not support nonlinear problems (i.e., NLobjective and NLconstraint).") optimize!(gm)
    @test true
end

function test_bbc()
    model, x, y, z, a = test_model()
    nl_constr = @constraint(model, sum(x[4]^2 + x[5]^2) <= z)
    expr = :((x, y, z) -> sum(x[i] for i=1:4) - y[1] * y[2] + z)

    # "Sanitizing" data
    inp = Dict(x[1] => 1., x[2] => 2, x[3] => 3, x[4] => 4, x[5] => 1, y[1] => 5, y[2] => -6, y[3] => -7, z => 7)
    inp_dict = Dict(string(key) => val for (key, val) in inp)
    inp_df = DataFrame(inp_dict)
    @test data_to_DataFrame(inp) == data_to_DataFrame(inp_dict) == data_to_DataFrame(inp_df) == inp_df
    @test data_to_Dict(inp_df, model) == data_to_Dict(inp, model) == data_to_Dict(inp_dict, model) == inp

    # Test bbl creation
    @test isnothing(functionify(nl_constr))
    @test functionify(expr) isa Function
    bbls = [BlackBoxClassifier(constraint = nl_constr, vars = [x[4], x[5], z], expr_vars = [x[4], x[5], z]),
        BlackBoxClassifier(constraint = expr, vars = flat([x[1:4], y[1:2], z]),
                         expr_vars = [x,y,z])]

    # Evaluation (scalar)
    # Quadratic (JuMP compatible) constraint
    @test evaluate(bbls[1], inp) == evaluate(bbls[1], inp_dict) == evaluate(bbls[1], inp_df) == -10.
    # Nonlinear expression
    @test evaluate(bbls[2], inp) == evaluate(bbls[2], inp_dict) == evaluate(bbls[2], inp_df)

    # Evaluation (vector)
    inp_df = DataFrame(-5 .+ 10 .*rand(3, size(inp_df,2)), string.(keys(inp)))
    inp_dict = data_to_Dict(inp_df, model)
    @test evaluate(bbls[1], inp_dict) == evaluate(bbls[1], inp_df) == inp_df[!, "z"] -
                                                    inp_df[!, "x[4]"].^2 - inp_df[!, "x[5]"].^2
    @test evaluate(bbls[2], inp_dict) == evaluate(bbls[2], inp_df)

    # bbl CHECKS
    bbl = bbls[1]

    # Check unbounded sampling
    @test_throws OCTHaGOnException X = lh_sample(bbl, n_samples=100);

    # Check knn_sampling without previous samples
    @test_throws OCTHaGOnException knn_sample(bbl)

    # Check evaluation of samples
    samples = DataFrame(randn(10, length(bbl.vars)),string.(bbl.vars))
    vals = bbl(samples);
    @test vals ≈ -1*samples[!, "x[4]"].^2 - samples[!, "x[5]"].^2 + samples[!, "z"]

    # Checks different kinds of sampling
    bound!(model, Dict(z => [-Inf, 10]))
    X_bound = boundary_sample(bbl);
    @test size(X_bound, 1) == 2^(length(bbl.vars)+1)
    @test_throws OCTHaGOnException knn_sample(bbl)
    X_lh = lh_sample(bbl, lh_iterations=3);

    # Check sample_and_eval
    uniform_sample_and_eval!(bbl);

    # Sampling, learning and showing...
    set_param(bbl, :ignore_feasibility, true)
    learn_constraint!(bbl);

    # Check feasibility and accuracy
    @test 0 <= feasibility(bbl) <= 1
    @test 0 <= evaluate_accuracy(bbl) <= 1

    # Training a model
    jc = add_feas_constraints!(model, bbl.vars, bbl.learners[1]);
    @test true
end

""" Testing some IAI kwarging. """
function test_kwargs()
    # Classification kwargs first...
    sample_kwargs = Dict(:localsearch => false,
                       :invalid_kwarg => :hello,
                       :ls_num_tree_restarts => 20)

    dict_fit = fit_classifier_kwargs(; sample_kwargs...)
    dict_fit2 = fit_classifier_kwargs(localsearch = false, invalid_kwarg = :hello,
                           ls_num_tree_restarts = 20)
    @test dict_fit == dict_fit2 == Dict(:sample_weight => :autobalance)

    dict_lnr = classifier_kwargs(; sample_kwargs...)
    dict_lnr2 = classifier_kwargs(localsearch = false, invalid_kwarg = :hello, ls_num_tree_restarts = 20)
    @test dict_lnr == dict_lnr2 == Dict(:localsearch => false, :ls_num_tree_restarts => 20)

    # Regression kwargs next...
    dict_fit = fit_regressor_kwargs(; sample_kwargs)
    @test dict_fit == Dict()

    dict_lnr = regressor_kwargs(; sample_kwargs...)
    dict_lnr2 = regressor_kwargs(localsearch = false, invalid_kwarg = :hello, ls_num_tree_restarts = 20)
    @test dict_lnr == dict_lnr2 == Dict(:localsearch => false, :ls_num_tree_restarts => 20)
end

""" Tests different kinds of regression. """
function test_regress()
    X = DataFrame(:x => 3*rand(100) .- 1, :y => 3*rand(100) .- 1);
    Y = Array(X[!,:x].^3 .* sin.(X[!,:y]));
    α0, α = u_regress(X, Y, SOLVER_SILENT)
    β0, β = l_regress(X, Y, SOLVER_SILENT)
    γ0, γ = ridge_regress(X, Y, SOLVER_SILENT)
    lowers = β0 .+ Matrix(X) * β;
    uppers = α0 .+ Matrix(X) * α;
    best_fit = γ0 .+ Matrix(X) * γ;
    @test all(lowers .<= Y) && all(uppers .>= Y) && all(uppers .>= lowers)
    errors = [sum((lowers-Y).^2), sum((uppers-Y).^2), sum((best_fit-Y).^2)]
    @test errors[3] <= errors[1] && errors[3] <= errors[2]

    X = DataFrame(:x => rand(100), :y => rand(100))
    Y = X[!,:y] - X[!,:x] .+ 0.1
    β0, β = svm(Matrix(X), Y, SOLVER_SILENT)
    predictions = Matrix(X) * β .+ β0
    @test sum((predictions-Y).^2) <= 1e-10
end

""" Tests various ways to train a regressor. """
function test_bbr()
    gm = minlp(true)
    set_optimizer(gm, SOLVER_SILENT)
    uniform_sample_and_eval!(gm)
    bbr = gm.bbls[3]

    # Checking that sampling adds bounds to BBRs.
    #TODO: COMPLETE.
    # check_unbounds([bbr.])

    # Make sure to train and add bbcs, and forget for a while...
    bbcs = [bbl for bbl in gm.bbls if bbl isa BlackBoxClassifier]
    learn_constraint!(bbcs)
    for bbc in bbcs
        add_tree_constraints!(gm, bbc)
    end
    @test all(evaluate_accuracy.(bbcs) .>= 0.95)

    # Threshold training
    learn_constraint!(bbr, "upper" => 20.)
    lnr = bbr.learners[end]
    @test lnr isa IAI.OptimalTreeClassifier
    all_leaves = find_leaves(lnr)
    # Add a binary variable for each leaf
    feas_leaves =
        [i for i in all_leaves if Bool(IAI.get_classification_label(lnr, i))];
    @test sort(-1 .*collect(keys(bbr.ul_data[end]))) == sort(feas_leaves) # upper bounding leaves have negative idxs
    @test bbr.thresholds[end] == ("upper" => 20.)

    # Check adding of upper bounding constraint to empty model
    init_constraints = count_constraints(gm)
    init_variables = count_variables(gm)
    add_tree_constraints!(gm, bbr)
    @test sort(union(feas_leaves, [1])) == sort(abs.(collect(keys(bbr.mi_constraints))))
    @test sort(abs.(collect(keys(bbr.leaf_variables)))) == sort(feas_leaves)

    # Checking that correct numbers of constraints and variables are added
    final_constraints = count_constraints(gm.model)
    final_variables = count_variables(gm.model)
    @test final_constraints == init_constraints + length(all_mi_constraints(bbr)) + length(bbr.leaf_variables)
    @test final_variables == init_variables + length(bbr.leaf_variables)*(length(bbr.vars) + 2)

    # Check clearing of all constraints and variables
    clear_tree_constraints!(gm, bbr)
    @test init_constraints == count_constraints(gm.model)
    @test init_variables == count_variables(gm.model)
    # Since all_variables(gm) doesn't count auxiliary variables...
    @test length(all_variables(gm.model)) == length(all_variables(gm)) +
        length(bbcs[1].leaf_variables)*(length(bbcs[1].vars) + 1) +
        length(bbcs[2].leaf_variables)*(length(bbcs[2].vars) + 1)

    # Flat prediction training
    learn_constraint!(bbr, regression_sparsity = 0, max_depth = 2)
    lnr = bbr.learners[end]
    all_leaves = find_leaves(lnr)
    @test lnr isa IAI.OptimalTreeRegressor
    @test length(bbr.ul_data[end]) == length(all_leaves) * 2 # double since contains lower and upper data
    @test bbr.thresholds[end] == Pair("reg", nothing)

    # Full regression training
    learn_constraint!(bbr, max_depth = 0)
    lnr = bbr.learners[end]
    @test lnr isa IAI.OptimalTreeRegressor
    all_leaves = find_leaves(lnr)
    @test length(bbr.ul_data[end]) == length(all_leaves) * 2
    @test bbr.thresholds[end] == Pair("reg", nothing)

    # Lower regression training
    learn_constraint!(bbr, "lower" => 5.)
    @test all(sign.(collect(keys(bbr.ul_data[end]))) .== 1)
    @test bbr.thresholds[end] == ("lower" => 5.)

    # upperlower regression training (sets an upper bound, but lower bounds in leaves)
    learn_constraint!(bbr, "upperlower" => 10.)
    lnr = bbr.learners[end]
    @test sum(collect(keys(bbr.ul_data[end]))) == 0

    # Checking all possible update scenarios
    # Note: [1] => upper bounds, [2,3] => regressors, [4] => lower bounds
    update_tree_constraints!(gm, bbr, 1) # Updating nothing with single upper bound
    @test bbr.active_trees == Dict(1 => bbr.thresholds[1]) &&
        all(sign.(collect(keys(bbr.mi_constraints))) .== -1)
    @test isnothing(active_lower_tree(bbr))
    @test active_upper_tree(bbr) == 1

    clear_tree_constraints!(gm, bbr)
    update_tree_constraints!(gm, bbr, 2) # Updating nothing with single regressor
    @test all(sign.(collect(keys(bbr.leaf_variables))) .== 1) &&
            length(bbr.leaf_variables) + 1 == length(bbr.mi_constraints) &&
                bbr.active_trees == Dict(2 => bbr.thresholds[2])
    @test active_lower_tree(bbr) == 2 && active_upper_tree(bbr) == 2

    update_tree_constraints!(gm, bbr, 3) # Replacing regressor with a regressor
    # No upper and lower bounding, just regressor here
    @test all(sign.(collect(keys(bbr.leaf_variables))) .== 1) &&
            bbr.active_trees == Dict(3 => bbr.thresholds[3])
    @test active_lower_tree(bbr) == 3 && active_upper_tree(bbr) == 3

    update_tree_constraints!(gm, bbr, 4) # Replacing regressor with lower bound
    @test all(sign.(collect(keys(bbr.leaf_variables))) .== 1) && # since lower bounding
        length(bbr.leaf_variables) + 1 == length(bbr.mi_constraints) &&
            bbr.active_trees == Dict(4 => bbr.thresholds[4])
    @test active_lower_tree(bbr) == 4
    @test isnothing(active_upper_tree(bbr))

    update_tree_constraints!(gm, bbr, 1) # Adding upper bound to lower bound
    @test length(bbr.leaf_variables) + 2 == length(bbr.mi_constraints) && # checking leaf variables
        length(bbr.active_trees) == 2
    optimize!(gm)
    @test active_lower_tree(bbr) == 4 && active_upper_tree(bbr) == 1

    clear_tree_constraints!(gm, bbr)
    update_tree_constraints!(gm, bbr, 1)
    update_tree_constraints!(gm, bbr, 4) # Adding lower bound to upper bound
    @test length(bbr.leaf_variables) + 2 == length(bbr.mi_constraints) && # checking leaf variables
        length(bbr.active_trees) == 2
    optimize!(gm)
    @test active_lower_tree(bbr) == 4 && active_upper_tree(bbr) == 1

    update_tree_constraints!(gm, bbr, 5) # Replace separate u/l bounds with upperlower
    @test length(bbr.leaf_variables)*2 + 1 == length(bbr.mi_constraints) &&
        length(bbr.active_trees) == 1
    @test active_lower_tree(bbr) == 5 && active_upper_tree(bbr) == 5

    # Make sure that all solutions are the same.
    update_tree_constraints!(gm, bbr, 1)
    update_tree_constraints!(gm, bbr, 4)
    optimize!(gm)
    update_tree_constraints!(gm, bbr, 1)
    optimize!(gm)

    # Testing leaf sampling
    df = last_leaf_sample(bbr)
    active_tree_idxs = collect(keys(bbr.active_trees))
    @test sort(active_tree_idxs) == [1, 4]
    @test size(df, 1) == get_param(bbr, :n_samples)
    for bbc in bbcs
        df = last_leaf_sample(bbc)
        @test size(df, 1) == get_param(bbc, :n_samples)
    end
    update_tree_constraints!(gm, bbcs[1])
    @test_throws OCTHaGOnException last_leaf_sample(bbcs[1])

    @test all(Array(gm.solution_history[:,"obj"]) .≈ gm.solution_history[1, "obj"])

    # Checking proper storage
    @test all(length.([bbr.ul_data, bbr.thresholds, bbr.learners, bbr.learner_kwargs]) .== 5)
    clear_tree_data!(bbr)
    @test all(length.([bbr.ul_data, bbr.thresholds, bbr.learners, bbr.learner_kwargs]) .== 0)
    @test !isempty(bbr.mi_constraints) && !isempty(bbr.leaf_variables)
    clear_tree_constraints!(gm, bbr)
    @test isempty(bbr.mi_constraints) && isempty(bbr.leaf_variables)
end

""" Tests basic functionalities in GMs. """
function test_basic_gm()
    gm = sagemark_to_GlobalModel(3, false)
    set_optimizer(gm, SOLVER_SILENT)

    # Actually trying to optimize...
    find_bounds!(gm, all_bounds=true)
    uniform_sample_and_eval!(gm)
    learn_constraint!(gm)

    # Solving of model
    set_param(gm, :ignore_accuracy, true)
    add_tree_constraints!(gm)
    println("Approximation accuracies: ", evaluate_accuracy(gm))
    optimize!(gm)
    vals = solution(gm);
    @test all(gm.bbls[i].active_leaves[1] in keys(gm.bbls[i].leaf_variables) for i=1:length(gm.bbls))
    println("X values: ", vals)
    println("Optimal X: ", vcat(exp.([5.01063529, 3.40119660, -0.48450710]), [-147-2/3]))

    # Testing constraint addition and removal
    clear_tree_constraints!(gm) # Clears all bbl constraints
    @test !any(is_valid(gm.model, constraint) for constraint in all_mi_constraints(gm.bbls[2]))
    add_tree_constraints!(gm, gm.bbls[2])
    @test all(is_valid(gm.model, constraint) for constraint in all_mi_constraints(gm.bbls[2]))
    clear_tree_constraints!(gm);
    add_tree_constraints!(gm);
    clear_tree_constraints!(gm, gm.bbls[1])
    @test !any(is_valid(gm.model, constraint) for constraint in all_mi_constraints(gm.bbls[1]))
    clear_tree_constraints!(gm) # Finds and clears the one remaining bbl constraint.
    @test all([!is_valid(gm.model, constraint) for constraint in gm.bbls[1].mi_constraints])
    @test all([!is_valid(gm.model, var) for var in values(gm.bbls[1].leaf_variables)])

    # Saving fit for test_load_fits()
    # save_fit(gm)

    # Testing clearing all data
    clear_data!(gm)
    @test all([size(bbl.X, 1) == 0 for bbl in gm.bbls])
    @test all([length(bbl.learners) == 0 for bbl in gm.bbls])

    # Testing surveysolve, and add_infeasibility_cuts for
    @test_throws OCTHaGOnException surveysolve(gm) # no data error
    uniform_sample_and_eval!(gm)
    surveysolve(gm)
    update_leaf_vexity(gm.bbls[1])
    @test true

    # Testing relaxations
    add_relaxation_variables!(gm)
    # Adding new variables for each nonlinear constraint
    @test sum(length(all_constraints(gm.model, type[1], type[2]))
            for type in JuMP.list_of_constraint_types(gm.model)) == 10
    relax_objective!(gm)
    surveysolve(gm)
    @test [getvalue(bbl.relax_var) for bbl in gm.bbls] == zeros(length(gm.bbls))
    clear_relaxation_variables!(gm)
    optimize!(gm)
    @test gm.cost[end] == gm.cost[end-1]
    @test length(JuMP.objective_function(gm.model).terms) == 1
    tight_objective!(gm)
    @test JuMP.objective_function(gm.model) isa JuMP.VariableRef
end

function test_convex_objective()
    gm = test_gqp()
    bbl = gm.bbls[1]
    set_param(bbl, :n_samples, 100)
    uniform_sample_and_eval!(gm)
    update_gradients(bbl, [1,2])
    @test !any(ismissing.(flat(Matrix(bbl.X[1:2,:]))))
    update_gradients(bbl)
    hand_calcs = DataFrame(hcat([[6*x[1] + 2*x[2] + 1, 2*x[2] + 2*x[1] + 6] for x in eachrow(Matrix(bbl.X))]...)', string.(bbl.vars))
    @test all(all(Array(bbl.gradients[i,:]) .≈ Array(hand_calcs[i,:])) for i = 1:100)
    update_vexity(bbl)
    @test bbl.local_convexity == 1.0
    @test bbl.convex == true

    # "Learning" convex functions should not result in trees.
    learn_constraint!(bbl)
    @test isempty(bbl.learners)
    add_tree_constraints!(gm, bbl)
    @test length(bbl.mi_constraints[1]) >= 1
    @test isempty(bbl.leaf_variables)

    # Testing adding gradient cuts
    optimize!(gm)
    @test gm.bbls[1].active_leaves == []
    add_infeasibility_cuts!(gm)
    optimize!(gm)
    while abs(gm.cost[end] - gm.cost[end-1]) >  get_param(gm, :abstol)
        add_infeasibility_cuts!(gm)
        optimize!(gm)
    end
    @test all([gm.cost[i] <= gm.cost[i+1] for i= 1:length(gm.cost)-1])
    update_leaf_vexity(gm.bbls[1])
    @test gm.bbls[1].vexity[1][2] == 1.0
end

function test_data_driven()
    model, x, y, z, a = test_model()
    rand_data = DataFrame("d" => [1,2,3])
    @test_throws UndefVarError bound_to_data!(model, rand_data)
    add_variables_from_data!(model, rand_data)
    @test model[:d] isa JuMP.VariableRef
    bound_to_data!(model, rand_data)
    @test get_bounds(model)[model[:d]] == [1,3]

    # Testing more complex constraints from afpm_model
    gm = afpm_model()
    @test length(gm.vars) == 15
    @test isnothing(get_unbounds(gm))
end

""" Test RandomForest learners for BlackBoxRegression.
NOTE: Although rfreg is tested, IT IS NOT A SUPPORTED FEATURE SINCE IT
IS NOT USEFUL FOR MOST PROBLEMS. """
function test_rfs()
    gm = minlp(true)
    set_optimizer(gm, SOLVER_SILENT)
    uniform_sample_and_eval!(gm)
    init_constraints = sum(length(all_constraints(gm.model, type[1], type[2]))
        for type in JuMP.list_of_constraint_types(gm.model))
    bbr = gm.bbls[3]
    for bbl in gm.bbls
        if bbl isa BlackBoxClassifier
            learn_constraint!(bbl)
            add_tree_constraints!(gm, bbl)
        elseif bbl isa BlackBoxRegressor
            learn_constraint!(bbl, "rfreg" => nothing)
            add_tree_constraints!(gm, bbl)
        end
    end
    optimize!(gm)
    descend!(gm)
    clear_tree_constraints!(gm)
    clear_relaxation_variables!(gm)
    # TODO: for some reason, relaxation variables are not removed from just clear_tree_constraints! Figure out.
    # TODO: fix this test. On the backburner though, since RFs are not fully supported.
    @test init_constraints == sum(length(all_constraints(gm.model, type[1], type[2]))
                                    for type in JuMP.list_of_constraint_types(gm.model))
end

""" Tests convexity checking. """
function test_convexcheck()
    # An example from : https://arxiv.org/pdf/0903.1287v1.pdf
    m = Model()
    @variable(m, -1 <= x[1:3] <= 1)
    @variable(m, -10 <= y <= 10)
    gm = GlobalModel(model = m)
    add_nonlinear_constraint(gm, :((x,y) -> y-(32*x[1]^8 + 118*x[1]^6*x[2]^2 + 40*x[1]^6*x[3]^2 + 25*x[1]^4*x[2]^4 -
                                43*x[1]^4*x[2]^2*x[3]^2 - 35*x[1]^4*x[3]^4 + 3*x[1]^2*x[2]^4*x[3]^2 -
                                16*x[1]^2*x[2]^2*x[3]^4 + 24*x[1]^2*x[3]^6 + 16*x[2]^8 +
                                44*x[2]^6*x[3]^2 + 70*x[2]^4*x[3]^4 + 60*x[2]^2*x[3]^6 + 30*x[3]^8)))
    uniform_sample_and_eval!(gm.bbls[end])
    update_vexity(gm.bbls[end])
    @test gm.bbls[end].convex

    # Easy examples to check method
    add_nonlinear_constraint(gm, :(x -> x[1] - 0.25x[2] + 0.01))
    uniform_sample_and_eval!(gm.bbls[end])
    update_vexity(gm.bbls[end])
    @test gm.bbls[end].convex

    # LSE
    add_nonlinear_constraint(gm, :(x -> 1 - 0.5*exp(-10*x[1]+x[2])))
    uniform_sample_and_eval!(gm.bbls[end])
    update_vexity(gm.bbls[end])
    @test gm.bbls[end].convex

    # Quadratic
    add_nonlinear_constraint(gm, :(x -> 1 - x[1]^2 - 5x[2]^2))
    uniform_sample_and_eval!(gm.bbls[end])
    update_vexity(gm.bbls[end])
    @test gm.bbls[end].convex

    # A simple nonconvex example
    add_nonlinear_constraint(gm, :(x -> minimum([1 - x[1]^2 - x[2]^2, x[1]^2 + x[2]^2 - 0.5])))
    uniform_sample_and_eval!(gm.bbls[end])
    update_vexity(gm.bbls[end])
    @test !gm.bbls[end].convex
end

function test_mergelinked()
    gm = hanoi_network(true)
    bins = detect_linked_constraints!(gm)
    @test length(bins) == 3
    merge_linked_constraints!(gm, bins)
    @test length(gm.bbls) == 3
    @test length([ll for bbl in gm.bbls for ll in bbl.lls]) == 99
end

function test_linking()
    gm = foxes_and_rabbits()
    init_constraints = count_constraints(gm)
    init_variables = count_variables(gm)
    set_param(gm, :max_iterations, 20)
    set_param(gm, :tighttol, 1e-3)
    set_param(gm, :equality_penalty, 1e5)
    set_param(gm, :step_penalty, 1e4)
    add_relaxation_variables!(gm)
    relax_objective!(gm)
    uniform_sample_and_eval!(gm)
    learn_constraint!(gm)
    add_tree_constraints!(gm)
    optimize!(gm)

    # # Plotting temporal population data (for visual debugging)
    # using Plots
    # plot(gm.soldict[:x], label = "Prey")
    # plot!(gm.soldict[:y], label = "Predators", xlabel = "Time", ylabel = "Normalized population")
    # # OR simultaneously in the population dimensions
    # plot!(gm.soldict[:x], gm.soldict[:y],
    #     xlabel = "Prey", ylabel = "Predators",
    #     label = 1:20, legend = false)

    # Checking constraint generation
    bbl_mic_count = length(all_mi_constraints(gm.bbls[1]))
    @test all(bbl_mic_count .== length(all_mi_constraints(ll)) for ll in gm.bbls[1].lls)
    clear_tree_constraints!(gm)
    clear_relaxation_variables!(gm)
    @test init_constraints == count_constraints(gm)
    @test init_variables == count_variables(gm)

    # Trying the same with big-M constraints (BBC with equalities)
    set_bigM(gm, true)
    add_relaxation_variables!(gm)
    relax_objective!(gm)
    add_tree_constraints!(gm)
    optimize!(gm)

    @test sum(isapprox.(values(gm.solution_history[1,:]),
                        values(gm.solution_history[2,:]), rtol = 1e-1))/size(gm.solution_history, 2) >= 0.9

    # Checking constraints generation for big-M
    clear_tree_constraints!(gm)
    clear_relaxation_variables!(gm)
    @test init_constraints == count_constraints(gm)
    @test init_variables == count_variables(gm)
end

function test_oos()
    op = oos_params()
    gm = oos_gm!()
    print_details(gm)
    init_variables = count_variables(gm)
    init_constraints = count_constraints(gm)
    uniform_sample_and_eval!(gm)
    learn_constraint!(gm, max_depth=8)
    for bbl in gm.bbls
        add_tree_constraints!(gm, bbl)
        clear_tree_constraints!(gm, bbl)
        init_constraints ==  count_constraints(gm) ||
            println("$(bbl.name)")
    end
    add_tree_constraints!(gm)
    optimize!(gm)

    # # Testing big-M formulation
    clear_tree_constraints!(gm)
    @test count_constraints(gm) == init_constraints
    @test count_variables(gm) == init_variables
    set_bigM(gm, true)
    add_tree_constraints!(gm)
    optimize!(gm)
    bin_vals = round.(JuMP.getvalue.(gm.model[:xx]))
    clear_tree_constraints!(gm)
    @test all(isapprox.(Array(gm.solution_history[1,:]),
                             Array(gm.solution_history[2,:]), rtol = 1e-2))
    @test count_constraints(gm) == init_constraints

    # NOTE: PGD may fail due to CPLEX issues.
    fix.(gm.model[:xx], bin_vals, force = true)
    descend!(gm, max_iterations = 2, tighttol = 1e-5, step_penalty = 1e8, equality_penalty = 1e6)

    # Post-processing
    plot_r = round.((getvalue.(gm.soldict[:r_orbit]) .- op.rE)./1e3, sigdigits = 5)
    ords = Int.(round.(getvalue.(gm.soldict[:sat_order])))
    tas = round.(getvalue.(gm.soldict[:ta]), sigdigits=3)
    times = round.(getvalue.(gm.soldict[:t_maneuver])./(24*3600), sigdigits=3) # days
    revs = round.(getvalue.(gm.soldict[:N_orbit]), sigdigits=5)
    refuels = round.(getvalue.(gm.soldict[:fuel_needed]), sigdigits=5)
    transfuels = round.(getvalue.(gm.soldict[:masses][:,1] .- gm.soldict[:masses][:,5]), sigdigits=5)
    n_trans = op.n_sats - 1

    # # Printing results
    println("Orbit altitudes (km) : $(plot_r)")
    println("Satellite order: $(ords)")
    println("True anomalies (radians): $(tas)")
    println("Orbital revolutions: $(revs)")
    println("Time for maneuvers (days): $(times)")
    println("Maneuver fuel (kg): $(transfuels)")
    println("Total mission time (years): $(sum(times)/365)")
    print_feas_gaps(gm)
    @test true

    # using Plots
    # colors = palette(:darktest, n_trans)

    # # # Orbit plotting
    # # plts = []
    # # using Plots
    # # colors = palette(:darktest, n_trans)
    # # for i=1:n_trans
    # #     plt = plot(x -> (op.r_sat-op.rE)/1e3,  0, 2pi, proj = :polar, label = false)
    # #     plt = plot!(x -> plot_r[i], 0, 2pi, proj = :polar, lims = (750, 870),
    # #                 color = colors[i], title = "Transfer $(i)", titlefontsize = 9, label = false,
    # #                 xtickfontsize=2, ylabel="altitude", yguidefontsize=5)
    # #     push!(plts, plt)
    # # end
    # # orb_plots = plot([plt for plt in plts]...)
    # # @show orb_plots

    # # Bar plotting
    # using Plots
    # colors = palette(:darktest, n_trans)
    # bar1 = bar(1:op.n_sats, refuels,  xlabel = "Satellite order", ylabel = "Satellite refuel (kg)", legend = false)
    # bar2 = bar(1:n_trans, times, xlabel = "Transfer index", ylabel = "Maneuver time (days)", legend = false)
    # bar3 = bar(1:n_trans, transfuels, xlabel = "Transfer index", ylabel = "Transfer fuel (kg)", legend = false)
    # bar_plots = plot(bar1, bar2, bar3, layout = (1,3))
    # @show bar_plots

    # Trying Ipopt's hand (doesn't work)
    # using Ipopt
    # m = oos_gm!()
    # JuMP.unset_binary.(m.model[:xx])
    # @constraint(m.model, m.model[:sat_order] .==  [4, 3, 2, 1, 7, 6, 5])
    # nonlinearize!(m)
    # uniform_sample_and_eval!(m)
    # set_optimizer(m, Ipopt.Optimizer)
    # optimize!(m)
end

test_expressions()

test_variables()

test_bounds()

test_sets()

test_linearize()

test_nonlinearize()

test_bbc()

test_kwargs()

test_regress()

test_bbr()

test_basic_gm()

test_convex_objective()

test_data_driven()

test_rfs()

test_convexcheck()

test_mergelinked()

test_linking()

test_oos()


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(
        full_bridge_optimizer,
        (MOI.AbstractOptimizer, Type{Float64}),
    )
    Base.precompile(
        unbridged_function,
        (LazyBridgeOptimizer, MOI.AbstractScalarFunction),
    )   # time: 0.1381832
    return
end

_precompile_()


text:
"""
    eom_henon2(dx, x, p, n)

Equations of motion for a 2D Henon system consisting of two identical
Henon maps with unidirectional forcing from `X` to`Y` [1].

# References

1. Krakovská, A., Jakubík, J., Chvosteková, M., Coufal, D.,
Jajcay, N., & Paluš, M. (2018). Comparison of six methods for the detection of
causality in a bivariate time series. Physical Review E, 97(4), 042207.
"""
function eom_henon2(x, p, n)
    c = p[1]
    x₁, x₂, y₁, y₂ = (x...,)
    dx₁ = 1.4 - x₁^2 + 0.3*x₂
    dx₂ = x₁
    dy₁ = 1.4 - (c * x₁ * y₁  +  (1 - c)*y₁^2) + 0.3*y₂
    dy₂ = y₁
    return SVector{4}(dx₁, dx₂, dy₁, dy₂)
end

"""
    henon2(u₀, c) -> DiscreteDynamicalSystem

Initialize a 2D Henon system consisting of two identical Henon maps with
unidirectional forcing from `X` to`Y` [1].

The difference equations are:

```math
\\begin{aligned}
x_1(t+1) &= 1.4 - x_1^2(t) + 0.3x_2(t) \\
x_2(t+1) &= x_1(t) \\
y_1(t+1) &= 1.4 - [c x_1(t) y_1(t) + (1-c) y_1^2(t)] + 0.3 y_2(t) \\
y_2(t+1) &= y_1(t)
\\end{aligned}
```

This system was investigated by Krakovská to study the performance of different
causality detection algorithms.

# References
Krakovská, A., Jakubík, J., Chvosteková, M., Coufal, D., Jajcay, N., & Paluš, M. (2018). Comparison of six methods for the detection of causality in a bivariate time series. Physical Review E, 97(4), 042207.
"""
function henon2(u₀, c)
    p = [c]
    logistic_system = DiscreteDynamicalSystem(eom_henon2, u₀, p)
    return logistic_system
end

"""
    henon2(;u₀ = rand(4), c = 2.0) -> DiscreteDynamicalSystem

Initialize a 2D Henon system consisting of two identical Henon maps with
unidirectional forcing from `X` to`Y` [1].

The difference equations are:

```math
\begin{aligned}
x_1(t+1) &= 1.4 - x_1^2(t) + 0.3x_2(t) \\
x_2(t+1) &= x_1(t) \\
y_1(t+1) &= 1.4 - [c*x_1(t)*y_1(t) + (1-c)*y_1^2(t)] + 0.3*y_2(t) \\
y_2(t+1) &= y_1(t)
\end{aligned}
```

This system was investigated by Krakovská to study the performance of different
causality detection algorithms.

# References
Krakovská, A., Jakubík, J., Chvosteková, M., Coufal, D., Jajcay, N., & Paluš, M. (2018). Comparison of six methods for the detection of causality in a bivariate time series. Physical Review E, 97(4), 042207.
"""
henon2(;u₀ = rand(4), c = 2.0) = henon2(u₀, c)


text:
using CSV, DataFrames, Gadfly, StatsBase

df = CSV.read("distances_tiff.csv")

function lincor(x, y)
    n = length(x)
    prod = sum(x .* y)
    ux, sx = mean_and_std(x)
    uy, sy = mean_and_std(y)
    cor = prod / n - ux * uy
    cor = cor / (sx * sy)
    return cor
end
function getKind(df::DataFrame, k::String)
    tmp = filter(e->e[:kind]==k, df)
    deletecols!(tmp, :kind)
    rename!(tmp, Dict(:value => Symbol(k)))
    sort!(tmp, Symbol(k))
    return tmp
end

symbolsStr = ["c2c_full_eucli", "c2c_proj_eucli", "c2c_maha",
              "s2s_full_eucli", "s2s_proj_eucli", "s2s_maha"]
println("\nDistance-Hist_Div correlations:")
for str in symbolsStr
    sym = Symbol(str)
    div_df = join(getKind(df, str), getKind(df, "hist_divergence"), on = [:index1, :index2])
    div_df[:id] = 1:(size(div_df)[1])
    s_cor = corspearman(div_df[sym], div_df[:hist_divergence])
    p_cor = lincor(div_df[sym], div_df[:hist_divergence])
    println("[$(str)]\tspearman: $(s_cor)\tp_cor: $(p_cor) .")
end

df_divs = join(getKind(df, "divergence"), getKind(df, "mle_divergence"),
                   on = [:index1, :index2])
println(size(df_divs))
df_divs = join(df_divs, getKind(df, "hist_divergence"), on = [:index1, :index2])
println(size(df_divs))

s_cor = corspearman(df_divs[:divergence], df_divs[:hist_divergence])
p_cor = lincor(df_divs[:divergence], df_divs[:hist_divergence])
println("[Normal-Hist divergences correlation]\tSpearman: $(s_cor)\tPearson: $(p_cor) .")

s_cor = corspearman(df_divs[:mle_divergence], df_divs[:hist_divergence])
p_cor = lincor(df_divs[:mle_divergence], df_divs[:hist_divergence])
println("[MLE-Hist divergences correlation]\tSpearman: $(s_cor)\tPearson: $(p_cor) .")

s_cor = corspearman(df_divs[:mle_divergence], df_divs[:divergence])
p_cor = lincor(df_divs[:mle_divergence], df_divs[:divergence])
println("[MLE-Normal divergences correlation]\tSpearman: $(s_cor)\tPearson: $(p_cor) .")

max_stars_repo_name: Jiliac/hemipt

max_issues_repo_name: Jiliac/hemipt

max_forks_repo_name: Jiliac/hemipt


text:
function test_with_reference(X, y, ref)
  for i in 1:nrow(ref)
    loss = eval(symbol(ref[i, :loss]))()
    regtype = eval(symbol(ref[i, :reg]))
    solver = eval(symbol(ref[i, :solver]))()
    tol = ref[i, :tol]
    C = ref[i, :C]
    iter = ref[i, :iter]
    fval = ref[i, :f_val]
    tnsv = ref[i, :nsv]
    acc = ref[i, :acc]
    msg2(ref[i, :solver], ", ", ref[i, :loss], ", C = "); @printf("%0.0e : ", C)
    retdual = svm(X, y, C = C, ftol = tol, iterations = 50_000, dual = true, solver = solver, regtype = regtype, loss = loss)
    retprimal = svm(X, y, C = C, ftol = tol, iterations = 50_000, dual = false, solver = solver, regtype = regtype, loss = loss)
    println("niters = $(iterations(retdual)) (ref: $iter)")
    if C > 1.
      # Dual prediction can take a while if there are a lot of SV
      @test acc - tol <= _accuracy(classify(retdual, X), y) <= acc + tol
    end
    @test acc - tol <= _accuracy(classify(retprimal, full(X)), y) <= acc + tol
    @test tnsv == nsv(retdual)
    @test tnsv - 2 <= nsv(retprimal) <= tnsv + 2
    @test fval - 0.00001 <= minimum(retdual) <= fval + 0.00001
  end
end

for csvname in getFilenames(_referenceDir(), "csv")
  dataname = splitext(csvname)[1]
  loadIfExists(test_with_reference, dataname, "Dense: Compare accuracy to Reference", sparse = false)
end

for dataname in ["a1a", "a9a", "data1.git", "data2.git"]
  loadIfExists(test_with_reference, dataname, "Sparse: Compare accuracy to Reference", sparse = true)
end


text:
"""
    cor_bounds(dA::UnivariateDistribution, dB::UnivariateDistribution, C::Type{<:Correlation}=Pearson; n_samples::Real=100_000)

Compute the stochastic lower and upper correlation bounds between two marginal
distributions.

This method relies on sampling from each distribution and then estimating the
specified correlation between the sorted samples. Because the samples are random,
there will be some variation in the answer for each call to `cor_bounds`. Increasing
the number of samples will increase the accuracy of the estimate, but will also
take longer to sort. Therefore ≈100,000 samples (the default) are recommended so
that it runs fast while still returning a good estimate.

The possible correlation types are:
  * [`Pearson`](@ref)
  * [`Spearman`](@ref)
  * [`Kendall`](@ref)

# Examples
```julia-repl
julia> using Distributions

julia> A = Normal(78, 10); B = LogNormal(3, 1);

julia> cor_bounds(A, B)
(lower = -0.7646512417819491, upper = 0.7649206069306482)

julia> cor_bounds(A, B, n_samples=Int(1e9))
(lower = -0.7629776825238167, upper = 0.7629762333824238)

julia> cor_bounds(A, B, n_samples=Int(1e4))
(lower = -0.7507010142250724, upper = 0.7551879647701095)

julia> cor_bounds(A, B, Spearman)
(lower = -1.0, upper = 1.0)
```
"""
function cor_bounds(dA::UD, dB::UD, C::Type{<:Correlation}=Pearson; n_samples::Real=100_000)
    a = rand(dA, Int(n_samples))
    b = rand(dB, Int(n_samples))

    sort!(a)
    sort!(b)
    upper = cor(a, b, C)

    reverse!(b)
    lower = cor(a, b, C)

    (lower = lower, upper = upper)
end

"""
    cor_bounds(margins::Vector{<:UD}, C::Type{<:Correlation}=Pearson; n_samples::Real=100_000)

Compute the stochastic pairwise lower and upper correlation bounds between a set
of marginal distributions.

The possible correlation types are:
  * [`Pearson`](@ref)
  * [`Spearman`](@ref)
  * [`Kendall`](@ref)
"""
function cor_bounds(margins::Vector{<:UD}, C::Type{<:Correlation}=Pearson; n_samples::Real=100_000)
    d = length(margins)
    lower, upper = zeros(Float64, d, d), zeros(Float64, d, d)
    n_samples = Int(n_samples)

    @threads for i in collect(subsets(1:d, Val{2}()))
        l, u = cor_bounds(margins[i[1]], margins[i[2]], C, n_samples=n_samples)
        lower[i...] = l
        upper[i...] = u
    end

    lower .= Symmetric(lower)
    cor_constrain!(lower)

    upper .= Symmetric(upper)
    cor_constrain!(upper)

    (lower = lower, upper = upper)
end


text:
using Test
using StaticArrays: SVector
using QuantumLattices.Essentials.SpinPackage
using QuantumLattices.Essentials.Spatials: PID,Point,Bond
using QuantumLattices.Essentials.DegreesOfFreedom: Table,OID,isHermitian,IDFConfig,Operators,oidtype,otype,script,iid,latexformat
using QuantumLattices.Essentials.Terms: Couplings,@subscript,statistics,abbr
using QuantumLattices.Interfaces: dims,inds,expand,matrix,permute,rank
using QuantumLattices.Prerequisites: Float
using QuantumLattices.Mathematics.Combinatorics: Permutations
using QuantumLattices.Mathematics.AlgebraOverFields: ID
using QuantumLattices.Mathematics.VectorSpaces: IsMultiIndexable,MultiIndexOrderStyle

@testset "SID" begin
    @test SID|>fieldnames==(:orbital,:spin,:tag)
    @test SID(orbital=1,spin=0.5,tag='z')'==SID(orbital=1,spin=0.5,tag='z')
    @test SID(orbital=1,spin=0.5,tag='x')'==SID(orbital=1,spin=0.5,tag='x')
    @test SID(orbital=1,spin=0.5,tag='y')'==SID(orbital=1,spin=0.5,tag='y')
    @test SID(orbital=1,spin=0.5,tag='+')'==SID(orbital=1,spin=0.5,tag='-')
    @test SID(orbital=1,spin=0.5,tag='-')'==SID(orbital=1,spin=0.5,tag='+')
end

@testset "matrix" begin
    @test isapprox(matrix(SID(1,0.5,'z')),[[-0.5,0.0] [0.0,0.5]])
    @test isapprox(matrix(SID(1,0.5,'x')),[[0.0,0.5] [0.5,0.0]])
    @test isapprox(matrix(SID(1,0.5,'y')),[[0.0,-0.5im] [0.5im,0.0]])
    @test isapprox(matrix(SID(1,0.5,'+')),[[0.0,1.0] [0.0,0.0]])
    @test isapprox(matrix(SID(1,0.5,'-')),[[0.0,0.0] [1.0,0.0]])

    @test isapprox(matrix(SID(1,1.0,'z')),[[-1.0,0.0,0.0] [0.0,0.0,0.0] [0.0,0.0,1.0]])
    @test isapprox(matrix(SID(1,1.0,'x')),[[0.0,√2/2,0.0] [√2/2,0.0,√2/2] [0.0,√2/2,0.0]])
    @test isapprox(matrix(SID(1,1.0,'y')),[[0.0,-√2im/2,0.0] [√2im/2,0.0,-√2im/2] [0.0,√2im/2,0.0]])
    @test isapprox(matrix(SID(1,1.0,'+')),[[0.0,√2,0.0] [0.0,0.0,√2] [0.0,0.0,0.0]])
    @test isapprox(matrix(SID(1,1.0,'-')),[[0.0,0.0,0.0] [√2,0.0,0.0] [0.0,√2,0.0]])
end

@testset "Spin" begin
    spin=Spin(atom=1,norbital=2,spin=1.0)
    @test IsMultiIndexable(Spin)==IsMultiIndexable(true)
    @test MultiIndexOrderStyle(Spin)==MultiIndexOrderStyle('C')
    @test dims(spin)==(2,5)
    @test inds(SID(1,1.0,'z'),spin)==(1,3)
    @test SID((1,1),spin)==SID(1,1.0,'x')
    @test spin|>collect==[  SID(1,1.0,'x'),SID(1,1.0,'y'),SID(1,1.0,'z'),SID(1,1.0,'+'),SID(1,1.0,'-'),
                            SID(2,1.0,'x'),SID(2,1.0,'y'),SID(2,1.0,'z'),SID(2,1.0,'+'),SID(2,1.0,'-')
                            ]
end

@testset "SIndex" begin
    @test SIndex|>fieldnames==(:scope,:site,:orbital,:spin,:tag)
    @test union(PID{Char},SID)==SIndex{Char}
end

@testset "script" begin
    @test script(OID(SIndex('S',1,2,0.5,'z')),Val(:site))==1
    @test script(OID(SIndex('S',1,2,0.5,'z')),Val(:orbital))==2
    @test script(OID(SIndex('S',1,2,0.5,'z')),Val(:spin))==0.5
    @test script(OID(SIndex('S',1,2,0.5,'z')),Val(:tag))=='z'
end

@testset "oidtype" begin
    @test oidtype(SID,Point{2,PID{Int}},Nothing,Val(true))==OID{SIndex{Int},SVector{2,Float},SVector{2,Float},Nothing}
    @test oidtype(SID,Point{2,PID{Int}},Table,Val(true))==OID{SIndex{Int},SVector{2,Float},SVector{2,Float},Int}
    @test oidtype(SID,Point{2,PID{Int}},Nothing,Val(false))==OID{SIndex{Int},Nothing,Nothing,Nothing}
    @test oidtype(SID,Point{2,PID{Int}},Table,Val(false))==OID{SIndex{Int},Nothing,Nothing,Int}
end

@testset "SOperator" begin
    @test latexformat(SOperator)==soptdefaultlatex
    opt=SOperator(1.0,(SIndex('a',1,1,0.5,'+'),SIndex('a',1,1,0.5,'-')))
    @test opt|>statistics==opt|>typeof|>statistics=='B'
    @test opt'==SOperator(1.0,(SIndex('a',1,1,0.5,'+'),SIndex('a',1,1,0.5,'-')))
    @test isHermitian(opt)
    @test repr(opt)=="S^{+}_{1,1}S^{-}_{1,1}"
end

@testset "permute" begin
    soptrep(opt::SOperator)=opt.value*prod([matrix(opt.id[i].index|>iid) for i=1:rank(opt)])
    for S in (0.5,1.0,1.5)
        oids=[OID(SIndex('S',1,2,S,tag),seq=i) for (i,tag) in enumerate(('x','y','z','+','-'))]
        table=Dict(oid.index=>oid.seq for oid in oids)
        for (id1,id2) in Permutations{2}(oids)
            left=soptrep(SOperator(1,ID(id1,id2)))
            right=sum([soptrep(opt) for opt in permute(SOperator,id1,id2,table)])
            @test isapprox(left,right)
        end
    end
    id1=OID(SIndex('S',1,2,0.5,'z'))
    id2=OID(SIndex('S',2,2,0.5,'z'))
    @test permute(SOperator,id1,id2,nothing)==(SOperator(1,ID(id2,id1)),)
end

@testset "SCID" begin
    @test SCID(center=1,atom=1,orbital=1,tag='z')|>typeof|>fieldnames==(:center,:atom,:orbital,:tag,:subscript)
end

@testset "SpinCoupling" begin
    @test SpinCoupling{2}(1.0,tags=('+','-'))|>string=="SpinCoupling(value=1.0,tags=(+,-))"
    @test SpinCoupling{2}(1.0,atoms=(1,1),tags=('z','z'))|>string=="SpinCoupling(value=1.0,atoms=(1,1),tags=(z,z))"
    @test SpinCoupling{2}(1.0,atoms=(1,1),orbitals=(1,2),tags=('-','+'))|>string=="SpinCoupling(value=1.0,atoms=(1,1),orbitals=(1,2),tags=(-,+))"
    @test SpinCoupling{2}(2.0,tags=('x','y'))|>repr=="2.0 SxSy"

    sc1=SpinCoupling{2}(1.5,tags=('+','-'),atoms=(1,2),orbitals=(@subscript (x,y) with x>y),centers=(1,2))
    sc2=SpinCoupling{2}(2.0,tags=('+','-'),atoms=(1,2),orbitals=(@subscript (x,y) with x<y),centers=(1,2))
    @test sc1|>repr=="1.5 S+S- sl(1,2)⊗ob(x,y)@(1,2) with $(sc1.id[1].subscript) && $(sc1.id[2].subscript)"
    @test sc2|>repr=="2.0 S+S- sl(1,2)⊗ob(x,y)@(1,2) with $(sc2.id[1].subscript) && $(sc2.id[2].subscript)"

    sc=sc1*sc2
    @test sc|>repr=="3.0 S+S-S+S- sl(1,2,1,2)⊗ob(x,y,x,y)@(1,2,1,2) with $(sc1.id[1].subscript) && $(sc1.id[2].subscript) && $(sc2.id[1].subscript) && $(sc2.id[2].subscript)"

    ex=expand(SpinCoupling{2}(2.0,tags=('+','-'),atoms=(1,1)),PID(1,1),Spin(atom=2,norbital=2,spin=1.0))
    @test collect(ex)==[]

    ex=expand(SpinCoupling{2}(2.0,tags=('+','-'),atoms=(1,2),orbitals=(1,2)),(PID(1,1),PID(1,2)),(Spin(atom=1,norbital=2,spin=1.0),Spin(atom=2,norbital=2,spin=1.0)))
    @test IsMultiIndexable(typeof(ex))==IsMultiIndexable(true)
    @test MultiIndexOrderStyle(typeof(ex))==MultiIndexOrderStyle('C')
    @test dims(ex)==(1,)
    @test collect(ex)==[(2.0,(SIndex(1,1,1,1.0,'+'),SIndex(1,2,2,1.0,'-')))]

    ex=expand(SpinCoupling{4}(2.0,tags=('+','-','+','-'),orbitals=(@subscript (α,α,β,β) with α<β)),PID(1,1),Spin(norbital=3,spin=1.0))
    @test dims(ex)==(3,)
    @test collect(ex)==[(2.0,(SIndex(1,1,1,1.0,'+'),SIndex(1,1,1,1.0,'-'),SIndex(1,1,2,1.0,'+'),SIndex(1,1,2,1.0,'-'))),
                        (2.0,(SIndex(1,1,1,1.0,'+'),SIndex(1,1,1,1.0,'-'),SIndex(1,1,3,1.0,'+'),SIndex(1,1,3,1.0,'-'))),
                        (2.0,(SIndex(1,1,2,1.0,'+'),SIndex(1,1,2,1.0,'-'),SIndex(1,1,3,1.0,'+'),SIndex(1,1,3,1.0,'-')))
    ]
end

@testset "Heisenberg" begin
    @test Heisenberg(orbitals=(1,2))==Couplings(SpinCoupling{2}(1//1,tags=('z','z'),orbitals=(1,2)),
                                                SpinCoupling{2}(1//2,tags=('+','-'),orbitals=(1,2)),
                                                SpinCoupling{2}(1//2,tags=('-','+'),orbitals=(1,2))
                                                )
    @test Heisenberg("xyz")==Couplings( SpinCoupling{2}(1,tags=('x','x')),
                                        SpinCoupling{2}(1,tags=('y','y')),
                                        SpinCoupling{2}(1,tags=('z','z'))
                                        )
end

@testset "Ising" begin
    @test Ising('x',atoms=(1,2))==Couplings(SpinCoupling{2}(1,tags=('x','x'),atoms=(1,2)))
    @test Ising('y',atoms=(1,2))==Couplings(SpinCoupling{2}(1,tags=('y','y'),atoms=(1,2)))
    @test Ising('z',atoms=(1,2))==Couplings(SpinCoupling{2}(1,tags=('z','z'),atoms=(1,2)))
end

@testset "Gamma" begin
    @test Gamma('x',orbitals=(1,1))==SpinCoupling{2}(1,tags=('y','z'),orbitals=(1,1))+SpinCoupling{2}(1,tags=('z','y'),orbitals=(1,1))
    @test Gamma('y',atoms=(1,2))==SpinCoupling{2}(1,tags=('z','x'),atoms=(1,2))+SpinCoupling{2}(1,tags=('x','z'),atoms=(1,2))
    @test Gamma('z')==SpinCoupling{2}(1,tags=('x','y'))+SpinCoupling{2}(1,tags=('y','x'))
end

@testset "DM" begin
    @test DM('x',orbitals=(1,1))==SpinCoupling{2}(1,tags=('y','z'),orbitals=(1,1))-SpinCoupling{2}(1,tags=('z','y'),orbitals=(1,1))
    @test DM('y',atoms=(1,2))==SpinCoupling{2}(1,tags=('z','x'),atoms=(1,2))-SpinCoupling{2}(1,tags=('x','z'),atoms=(1,2))
    @test DM('z')==SpinCoupling{2}(1,tags=('x','y'))-SpinCoupling{2}(1,tags=('y','x'))
end

@testset "Sᵅ" begin
    @test Sˣ(atom=1,orbital=1)==Couplings(SpinCoupling{1}(1,tags=('x',),atoms=(1,),orbitals=(1,)))
    @test Sʸ(atom=1)==Couplings(SpinCoupling{1}(1,tags=('y',),atoms=(1,)))
    @test Sᶻ(orbital=1)==Couplings(SpinCoupling{1}(1,tags=('z',),orbitals=(1,)))
end

@testset "spincoupling" begin
    sc=sc"1.0 S+S- sl(1,1)⊗ob(α,β)@(1,2) with α<β"
    ob=sc.subscripts[1].identifier
    @test repr(sc)=="1.0 S+S- sl(1,1)⊗ob(α,β)@(1,2) with $ob && $ob"

    sc=sc"1.0 S+S- sl(1,1)⊗ob(α,β)@(1,2)"
    ob=sc.subscripts[1].identifier
    @test repr(sc)=="1.0 S+S- sl(1,1)⊗ob(α,β)@(1,2) with $ob && $ob"

    sc=sc"1.0 S+S- sl(1,1)⊗ob(α,β)"
    ob=sc.subscripts[1].identifier
    @test repr(sc)=="1.0 S+S- sl(1,1)⊗ob(α,β) with $ob && $ob"

    sc=sc"1.0 S+S- sl(1,1)⊗ob(1,2)@(1,2)"
    @test repr(sc)=="1.0 S+S- sl(1,1)⊗ob(1,2)@(1,2)"

    sc=sc"1.0 S+S- sl(1,1)@(1,2)"
    @test repr(sc)=="1.0 S+S- sl(1,1)@(1,2)"

    sc=sc"1.0 S+S- @(1,2)"
    @test repr(sc)=="1.0 S+S- @(1,2)"

    sc=sc"1.0 S+S-"
    @test repr(sc)=="1.0 S+S-"
end

@testset "spincouplings" begin
    @test heisenberg"sl(1,1)⊗ob(1,3)@(1,2)"==Heisenberg(centers=(1,2),atoms=(1,1),orbitals=(1,3))
    @test heisenberg"sl(1,1)⊗ob(1,3)"==Heisenberg(atoms=(1,1),orbitals=(1,3))
    @test heisenberg"@(1,2)"==Heisenberg(centers=(1,2))
    @test heisenberg"ob(1,3)"==Heisenberg(orbitals=(1,3))
    @test heisenberg"sl(1,1)"==Heisenberg(atoms=(1,1))
    @test heisenberg""==heisenberg"+-z"==Heisenberg()
    @test heisenberg"xyz"==Heisenberg("xyz")

    @test ising"x"==Ising('x') && ising"x sl(1,1)⊗ob(1,3)@(1,2)"==Ising('x',centers=(1,2),atoms=(1,1),orbitals=(1,3))
    @test ising"y"==Ising('y') && ising"y sl(1,1)⊗ob(1,3)@(1,2)"==Ising('y',centers=(1,2),atoms=(1,1),orbitals=(1,3))
    @test ising"z"==Ising('z') && ising"z sl(1,1)⊗ob(1,3)@(1,2)"==Ising('z',centers=(1,2),atoms=(1,1),orbitals=(1,3))

    @test gamma"x"==Gamma('x') && gamma"x sl(1,1)⊗ob(1,3)@(1,2)"==Gamma('x',centers=(1,2),atoms=(1,1),orbitals=(1,3))
    @test gamma"y"==Gamma('y') && gamma"y sl(1,1)⊗ob(1,3)@(1,2)"==Gamma('y',centers=(1,2),atoms=(1,1),orbitals=(1,3))
    @test gamma"z"==Gamma('z') && gamma"z sl(1,1)⊗ob(1,3)@(1,2)"==Gamma('z',centers=(1,2),atoms=(1,1),orbitals=(1,3))

    @test dm"x"==DM('x') && dm"x sl(1,1)⊗ob(1,3)@(1,2)"==DM('x',centers=(1,2),atoms=(1,1),orbitals=(1,3))
    @test dm"y"==DM('y') && dm"y sl(1,1)⊗ob(1,3)@(1,2)"==DM('y',centers=(1,2),atoms=(1,1),orbitals=(1,3))
    @test dm"z"==DM('z') && dm"z sl(1,1)⊗ob(1,3)@(1,2)"==DM('z',centers=(1,2),atoms=(1,1),orbitals=(1,3))

    @test sˣ""==Sˣ() && sˣ"sl(1)⊗ob(2)"==Sˣ(atom=1,orbital=2)
    @test sʸ""==Sʸ() && sʸ"sl(1)⊗ob(2)"==Sʸ(atom=1,orbital=2)
    @test sᶻ""==Sᶻ() && sᶻ"sl(1)⊗ob(2)"==Sᶻ(atom=1,orbital=2)
end

@testset "SpinTerm" begin
    term=SpinTerm{1}(:h,1.5,0,couplings=Sᶻ())
    @test term|>abbr==:sp
    @test otype(term|>typeof,OID{SIndex{Int},SVector{2,Float},SVector{2,Float},Nothing})==SOperator{Float,ID{OID{SIndex{Int},SVector{2,Float},SVector{2,Float},Nothing},1}}
    @test otype(term|>typeof,OID{SIndex{Int},SVector{2,Float},SVector{2,Float},Int})==SOperator{Float,ID{OID{SIndex{Int},SVector{2,Float},SVector{2,Float},Int},1}}

    point=Point(PID('a',1),(0.5,0.5),(0.0,0.0))
    config=IDFConfig{Spin}(pid->Spin(atom=pid.site%2,norbital=2,spin=0.5),[point.pid])
    table=Table(config,usualspinindextotuple)
    term=SpinTerm{1}(:h,1.5,0,couplings=Sᶻ())
    operators=Operators(SOperator(1.5,ID(OID(SIndex('a',1,1,0.5,'z'),[0.5,0.5],[0.0,0.0],1))),
                        SOperator(1.5,ID(OID(SIndex('a',1,2,0.5,'z'),[0.5,0.5],[0.0,0.0],2)))
    )
    @test expand(term,point,config,table)==operators

    bond=Bond(1,Point(PID('a',1),(0.0,0.0),(0.0,0.0)),Point(PID('b',1),(0.5,0.5),(0.0,0.0)))
    config=IDFConfig{Spin}(pid->Spin(atom=pid.site%2,norbital=2,spin=0.5),[bond.spoint.pid,bond.epoint.pid])
    table=Table(config,usualspinindextotuple)
    term=SpinTerm{2}(:J,1.5,1,couplings=Heisenberg())
    operators=Operators(SOperator(1.50,ID(OID(SIndex('b',1,2,0.5,'z'),[0.5,0.5],[0.0,0.0],4),OID(SIndex('a',1,2,0.5,'z'),[0.0,0.0],[0.0,0.0],2))),
                        SOperator(0.75,ID(OID(SIndex('b',1,2,0.5,'-'),[0.5,0.5],[0.0,0.0],4),OID(SIndex('a',1,2,0.5,'+'),[0.0,0.0],[0.0,0.0],2))),
                        SOperator(0.75,ID(OID(SIndex('b',1,1,0.5,'-'),[0.5,0.5],[0.0,0.0],3),OID(SIndex('a',1,1,0.5,'+'),[0.0,0.0],[0.0,0.0],1))),
                        SOperator(0.75,ID(OID(SIndex('b',1,1,0.5,'+'),[0.5,0.5],[0.0,0.0],3),OID(SIndex('a',1,1,0.5,'-'),[0.0,0.0],[0.0,0.0],1))),
                        SOperator(1.50,ID(OID(SIndex('b',1,1,0.5,'z'),[0.5,0.5],[0.0,0.0],3),OID(SIndex('a',1,1,0.5,'z'),[0.0,0.0],[0.0,0.0],1))),
                        SOperator(0.75,ID(OID(SIndex('b',1,2,0.5,'+'),[0.5,0.5],[0.0,0.0],4),OID(SIndex('a',1,2,0.5,'-'),[0.0,0.0],[0.0,0.0],2)))
    )
    @test expand(term,bond,config,table)==operators
end


text:
using spc
using Documenter

makedocs(;
    modules=[spc],
    authors="Daan Van Hauwermeiren",
    repo="https://github.com/DaanVanHauwermeiren/spc.jl/blob/{commit}{path}#L{line}",
    sitename="spc.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://DaanVanHauwermeiren.github.io/spc.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/DaanVanHauwermeiren/spc.jl",
)

max_stars_repo_name: DaanVanHauwermeiren/spc

max_issues_repo_name: DaanVanHauwermeiren/spc

max_forks_repo_name: DaanVanHauwermeiren/spc


text:
using SurrogateModelOptim
using PlotlyJS
using Statistics
using Distributed

# Optimize the test function Rosenbrock
function rosenbrock_2D(x)
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end
search_range=[(-5.0,5.0),(-5.0,5.0)]

result = smoptimize(rosenbrock_2D, search_range;
                    options=SurrogateModelOptim.Options(
                    iterations=25,
                    num_interpolants=20, #Preferably even number of added processes
                    num_start_samples=5,
                        ));

function plot_fun_2D(fun,sr,title)
    N = 51
    x = range(sr[1][1], stop = sr[1][2], length = N)
    y = range(sr[2][1], stop = sr[2][2], length = N)
    grid(x,y) = [x,y]
    z = @. fun(grid(x,y'))

    trace = surface(x=x,y=y,z=z, colorscale="Viridis")
    layout = Layout(title=title)
    p = plot(trace,layout)
end

display(plot_fun_2D(rosenbrock_2D,search_range,"Original function"))
display(plot_fun_2D(x->median(result.sm_interpolant(x)),search_range,"Estimated function"))


text:
using JuMP, GLPK

function takumi_LP(; verbose = true)
    model = Model(with_optimizer(GLPK.Optimizer))

    @variable(model, 20 <= x <= 150, Int)
    @variable(model, 10 <= y <= 60, Int)

    @objective(model, Max, 8x + 17y)
    @constraint(model, 8x + 18y <= 1000)

    if verbose
        print(model)
    end

    JuMP.optimize!(model)

    obj_value = JuMP.objective_value(model)
    x_value = JuMP.value(x)
    y_value = JuMP.value(y)

    if verbose
        println("Objective value: ", obj_value)
        println("x: ", x_value)
        println("y: ", y_value)
    end

end

takumi_LP(verbose = true)


text:
using ManifoldsBase
using ManifoldsBase: ℝ
using Test

struct ErrorTestManifold <: AbstractManifold{ℝ} end

function ManifoldsBase.check_size(::ErrorTestManifold, p)
    size(p) != (2,) && return DomainError(size(p), " size $p not (2,)")
    return nothing
end
function ManifoldsBase.check_size(::ErrorTestManifold, p, X)
    size(X) != (2,) && return DomainError(size(X), " size $X not (2,)")
    return nothing
end
function ManifoldsBase.check_point(::ErrorTestManifold, x)
    if any(u -> u < 0, x)
        return DomainError(x, "<0")
    end
    return nothing
end
function ManifoldsBase.check_vector(M::ErrorTestManifold, x, v)
    mpe = ManifoldsBase.check_point(M, x)
    mpe === nothing || return mpe
    if any(u -> u < 0, v)
        return DomainError(v, "<0")
    end
    return nothing
end

@testset "Domain errors" begin
    M = ErrorTestManifold()
    @test isa(ManifoldsBase.check_point(M, [-1, 1]), DomainError)
    @test isa(ManifoldsBase.check_size(M, [-1, 1, 1]), DomainError)
    @test isa(ManifoldsBase.check_size(M, [-1, 1], [1, 1, 1]), DomainError)
    @test ManifoldsBase.check_point(M, [1, 1]) === nothing
    @test !is_point(M, [-1, 1])
    @test !is_point(M, [1, 1, 1]) # checksize fails
    @test_throws DomainError is_point(M, [-1, 1, 1], true) # checksize errors
    @test is_point(M, [1, 1])
    @test_throws DomainError is_point(M, [-1, 1], true)

    @test isa(ManifoldsBase.check_vector(M, [1, 1], [-1, 1]), DomainError)
    @test ManifoldsBase.check_vector(M, [1, 1], [1, 1]) === nothing
    @test !is_vector(M, [1, 1], [-1, 1])
    @test !is_vector(M, [1, 1], [1, 1, 1])
    @test_throws DomainError is_vector(M, [1, 1], [-1, 1, 1], true)
    @test !is_vector(M, [1, 1, 1], [1, 1, 1], false, true)
    @test_throws DomainError is_vector(M, [1, 1, 1], [1, 1], true, true)
    @test is_vector(M, [1, 1], [1, 1])
    @test_throws DomainError is_vector(M, [1, 1], [-1, 1], true)
end


text:
using DataFrames, StructuralEquationModels, Symbolics,
    LinearAlgebra, SparseArrays, Optim, LineSearches,
    BenchmarkTools, CSV, Statistics

date = string(ARGS...)

cd("cfa")

include("functions.jl")

config = DataFrame(CSV.File("config.csv"))
config2 = copy(config)
config.backend .= "Optim.jl"
config2.backend .= "NLopt.jl"

config = [config; config2]

config = filter(row -> (row.backend == "Optim.jl") & (row.meanstructure == 0), config)

data_vec = read_files("data", get_data_paths(config))
par_vec = read_files("parest", get_data_paths(config))
# start_vec = read_files("start", get_data_paths(config))

##############################################

models = gen_models(config, data_vec)

fits = get_fits(models)

correct = compare_estimates(fits, par_vec, config)

config.correct = correct

##############################################
# using MKL

benchmarks = benchmark_models(models)

results = select(config, :Estimator, :n_factors, :n_items, :meanstructure, :backend, :correct)

results.median_time = median.(getfield.(benchmarks, :times))
results.mean_time = median.(getfield.(benchmarks, :times))
results.sd_time = std.(getfield.(benchmarks, :times))
results.n_repetitions = vec(getfield.(getfield.(benchmarks, :params), :samples))

# results.n_par = 2*(results.n_factors.*results.n_items) + results.n_factors.*(results.n_factors.-1)/2,

CSV.write("results/benchmarks_julia_"*date*".csv", results, delim = ";")


text:
if VERSION < v"1.0.0-rc1"
    import Base.expand
end
IMPLEMENT_ONE_ARG_FUNC(:expand, :expand)

if get_symbol(:basic_cse) != C_NULL
    function cse(exprs...)
        vec = convert(CVecBasic, exprs...)
        replacement_syms = CVecBasic()
        replacement_exprs = CVecBasic()
        reduced_exprs = CVecBasic()
        ccall((:basic_cse, libsymengine), Nothing, (Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid}),
            replacement_syms.ptr, replacement_exprs.ptr, reduced_exprs.ptr,vec.ptr)
        return replacement_syms, replacement_exprs, reduced_exprs
    end
else
    function cse(exprs...)
        error("libsymengine is too old")
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineCaptureResult <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    properties::Any # spec type: Union{ Nothing, VirtualMachineCaptureResultProperties } # spec name: properties

    function VirtualMachineCaptureResult(;id=nothing, properties=nothing)
        o = new()
        validate_property(VirtualMachineCaptureResult, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(VirtualMachineCaptureResult, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type VirtualMachineCaptureResult

const _property_map_VirtualMachineCaptureResult = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("properties")=>Symbol("properties"))
const _property_types_VirtualMachineCaptureResult = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("properties")=>"VirtualMachineCaptureResultProperties")
Base.propertynames(::Type{ VirtualMachineCaptureResult }) = collect(keys(_property_map_VirtualMachineCaptureResult))
Swagger.property_type(::Type{ VirtualMachineCaptureResult }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineCaptureResult[name]))}
Swagger.field_name(::Type{ VirtualMachineCaptureResult }, property_name::Symbol) =  _property_map_VirtualMachineCaptureResult[property_name]

function check_required(o::VirtualMachineCaptureResult)
    true
end

function validate_property(::Type{ VirtualMachineCaptureResult }, name::Symbol, val)
end


text:
using DelimitedFiles
using Plots
using Random

#using DataFrames
#using CSV

# read files and
# create two arrays to store the sequences and expression values
#arti_rep1 = CSV.read("/home/rljahn/Dropbox/Julia repo/mutational distance/arti_SDR_rep1.csv")

path = "/home/rljahn/FitnessLandscape_SD/data/Read Count Table/UnionTable/"
function theoretical_value(n::Int64)
    result::Float64 = 0.0
    for i::Int64 = 1:(n+1)
        result += (i-1) * binomial(n,i-1) * 3^(n-i+1)
        println(result)
    end
    result = result/(4^n)
    return result
end

function mutational_distance(file, shuffule::Bool = false)
    arti_rep1 = readdlm(path * file, ',')
    arti_seq = arti_rep1[2:end,1]
    arti_seq = convert(Array{String,1}, arti_seq)

    arti_exp = arti_rep1[2:end, end]
    for i in 1:length(arti_exp)
        if arti_exp[i] == "NA"
            arti_exp[i] = NaN
        end
    end
    arti_exp = convert(Array{Float64,1}, arti_exp)
    rng = MersenneTwister(1234);
    if shuffule == true
        println("Shuffled!!!")
        shuffle!(rng, arti_exp)
    end

    #gpmapping = Dict(zip(arti_seq,arti_exp))

    # a function to compute hamming distance
    # e.g. hamming_distance(arti_seq[1],arti_seq[end])
    function hamming_similarity(seq1::String,seq2::String, len::Int64 = 9)::Float64
        return sum([1 for m = 1:len if seq1[m] == seq2[m]])
        #if length(seq1) == length(seq1)
        #    return length(seq1) - sum(collect(seq1) .== collect(seq2))
        #end
    end

    # a function to classify sequences into bins/buckets
    # based on theirs expression values
    function bucket_classifier(array::Array{Float64,1}, mini::Float64 = 0.0, maxi::Float64 = 3.2, number_bins::Int64 = 20)::Array{Union{Nothing, Int64},1}
        function comparator(input_value::Float64, bins = mini:(maxi-mini)/number_bins:maxi )
            for i = 1:length(bins)
                if input_value > bins[i] && input_value <= bins[i+1]
                    return i
                end
            end
        end
        values = comparator.(array)
        return values
    end

    buckets = [[] for i = 1:20]
    bucket_indeces = bucket_classifier(arti_exp)
    for (index, i) in enumerate(bucket_indeces)
        if i != nothing
            push!(buckets[i],arti_seq[index])
        end
    end

    # check the variants in each bins
    println([length(buckets[i]) for i in 1:20])
    println("The first bin has: ", buckets[1])
    println("The last bin has: ", buckets[end])
    #dictionary = Dict(zip(arti_seq,bucket_classifier(arti_exp)))

    n::Int64 = 20
    distance_matrix = zeros(n,n)

    for i::Int64 = 1:n
        for j::Int64 = 1:i
            println("Current progress is ($i,$j)!!!")
            t_start = time()
            dist::Float64 = 0.0
            num::Float64 = 0.0
            if i == j
                for k::Int64 = 1:length(buckets[i])
                    for l::Int64 = 1:(k-1)
                        dist += (hamming_similarity(buckets[i][k],buckets[j][l])*2)
                        num += (1.0*2)
                    end
                    dist += hamming_similarity(buckets[i][k],buckets[j][k])
                    num += 1.0
                end
            else
                for seq1::String in buckets[i]
                    for seq2::String in buckets[j]
                        dist += hamming_similarity(seq1,seq2)
                        num += 1.0
                    end
                end
            end

            #for seq1::String in buckets[i]
            #    for seq2::String in buckets[j]
            #        dist += hamming_similarity(seq1,seq2)
            #        num += 1.0
            #    end
            #end

            println("similarity score = $dist")
            println("number = $num")
            if num != 0
                distance_matrix[i,j] = dist/num
            end
            println("result = $(distance_matrix[i,j])")
            tmp = time()-t_start
            println("Time used: $tmp seconds")
        end
    end

    distance_matrix = distance_matrix + transpose(distance_matrix)
    for i = 1: size(distance_matrix)[1]
        distance_matrix[i,i] /= 2
    end

    heatmap(distance_matrix, c=:RdBu, aspect_ratio =1)
    writedlm("/Users/imac/Desktop/Julia_repo/" * "matrix_"* file * "_20190311.csv", distance_matrix, ',')
    #savefig("/Users/imac/Desktop/Julia_repo/" * "heatmap_" * file * ".svg")
    #println("Maximum time used in each loop: $max_time")
end

println(">>>>>> arti !!!")
@time mutational_distance("arti_SDR_union_count25.csv",true)

println(">>>>>> dmsC !!!")
@time mutational_distance("dmsC_SDR_union_count25.csv",true)

println(">>>>>> fepB !!!")
@time mutational_distance("fepB_SDR_union_count25.csv",true)

function plot_part(file)
    table1 = readdlm("/Users/imac/Desktop/Julia_repo/" * "matrix_"* file * "_20190311.csv",  ',')
    heatmap(table1, c=cgrad([:blue, :green, :yellow, :orange, :red]), clim = (1.5,4), aspect_ratio =1)
    savefig("/Users/imac/Desktop/Julia_repo/" * "heatmap_part" * file * "_20190311.svg")
end

plot_part("arti_SDR_union_count25.csv")
plot_part("fepB_SDR_union_count25.csv")
plot_part("dmsC_SDR_union_count25.csv")