text:
"""
    fihess(A, E; fast = true, finite_infinite = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (At, Et, Q, Z, ν, blkdims)

Reduce the regular matrix pencil `A - λE` to an equivalent form `At - λEt = Q'*(A - λE)*Z` using
orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `At` and `Et` are in one of the
following block upper-triangular forms:

(1) if `finite_infinite = false`, then

                   | Ai-λEi |   *     |
        At - λEt = |--------|---------|,
                   |    O   | Af-λEf  |

where the `ni x ni` subpencil `Ai-λEi` contains the infinite elementary divisors and
the `nf x nf` subpencil `Af-λEf` contains the finite eigenvalues of the pencil `A-λE`.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vector `ν` contains the dimensions of the diagonal blocks
of the staircase form  `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[i]-ν[i+1]` for `i = 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i`
(with `ν[nb+1] = 0`).

The subpencil `Af-λEf` is with `Af` in an upper Hessenberg form and `Ef` nonsingular and upper triangular.

The dimensions of the diagonal blocks are returned in `blkdims = (ni, nf)`.

(2) if `finite_infinite = true`, then

                   | Af-λEf |   *    |
        At - λEt = |--------|--------|,
                   |   O    | Ai-λEi |

where the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular,
contains the finite eigenvalues of the pencil `A-λE` and the `ni x ni` subpencil `Ai-λEi`
contains the infinite elementary divisors.

The subpencil `Af-λEf` is with `Af` in an upper Hessenberg form and `Ef` nonsingular and upper triangular.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular  and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vectors `ν` contains the dimensions of the diagonal blocks
of the staircase form `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[nb-j+1]-ν[nb-j]` for `j = 1, 2, ..., nb` is the number of infinite elementary
divisors of degree `j` (with `ν[0] = 0`).

The dimensions of the diagonal blocks are returned in `blkdims = (nf, ni)`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` and `E`.

The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`.
Otherwise, `Q` is set to `nothing`.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` if `withZ = true`.
Otherwise, `Z` is set to `nothing`.
"""
function fihess(A::AbstractMatrix, E::AbstractMatrix;
                   fast::Bool = true, finite_infinite::Bool = false,
                   atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))),
                   rtol::Real = (size(A,1)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)),
                   withQ::Bool = true, withZ::Bool = true)

    n = LinearAlgebra.checksquare(A)
    n == LinearAlgebra.checksquare(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
    T = promote_type(eltype(A), eltype(E))
    T <: BlasFloat || (T = promote_type(Float64,T))

    A1 = copy_oftype(A,T)
    E1 = copy_oftype(E,T)

    withQ ? (Q = Matrix{T}(I,n,n)) : (Q = nothing)
    withZ ? (Z = Matrix{T}(I,n,n)) : (Z = nothing)

    ν, blkdims = fisplit!(A1, E1, Q, Z, missing, missing;
                         fast  = fast, finite_infinite = finite_infinite,
                         atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = withQ, withZ = withZ)
    (ilo, ihi) = finite_infinite ? (1, blkdims[1]) : (blkdims[1]+1, n)
    ilo > ihi || gghrd!(withQ ? 'V' : 'N', withZ ? 'V' : 'N', ilo, ihi, A1, E1, Q, Z)

    return A1, E1, Q, Z, ν, blkdims
end
"""
    fischur(A, E; fast = true, finite_infinite = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (At, Et, Q, Z, ν, blkdims)

Reduce the regular matrix pencil `A - λE` to an equivalent form `At - λEt = Q'*(A - λE)*Z` using
orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `At` and `Et` are in one of the
following block upper-triangular forms:

(1) if `finite_infinite = false`, then

                   | Ai-λEi |   *     |
        At - λEt = |--------|---------|,
                   |    O   | Af-λEf  |

where the `ni x ni` subpencil `Ai-λEi` contains the infinite elementary divisors and
the `nf x nf` subpencil `Af-λEf` contains the finite eigenvalues of the pencil `A-λE`.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vector `ν` contains the dimensions of the diagonal blocks
of the staircase form  `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[i]-ν[i+1]` for `i = 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i`
(with `ν[nb+1] = 0`).

The subpencil `Af-λEf` is with the pair `(Af,Ef)` in a generalized Schur form, with `Af` (quasi) upper triangular and
`Ef` nonsingular and upper triangular.

The dimensions of the diagonal blocks are returned in `blkdims = (ni, nf)`.

(2) if `finite_infinite = true`, then

                   | Af-λEf |   *    |
        At - λEt = |--------|--------|,
                   |   O    | Ai-λEi |

where the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular,
contains the finite eigenvalues of the pencil `A-λE` and the `ni x ni` subpencil `Ai-λEi`
contains the infinite elementary divisors.

The subpencil `Af-λEf` is with the pair `(Af,Ef)` in a generalized Schur form, with `Af` (quasi) upper triangular and
`Ef` nonsingular and upper triangular.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular
and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vectors `ν` contains the dimensions of the diagonal blocks
of the staircase form `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[nb-j+1]-ν[nb-j]` for `j = 1, 2, ..., nb` is the number of infinite elementary
divisors of degree `j` (with `ν[0] = 0`).

The dimensions of the diagonal blocks are returned in `blkdims = (nf, ni)`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` and `E`.

The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`.
Otherwise, `Q` is set to `nothing`.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` if `withZ = true`.
Otherwise, `Z` is set to `nothing`.
"""
function fischur(A::AbstractMatrix, E::AbstractMatrix;
                      fast::Bool = true, finite_infinite::Bool = false,
                      atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))),
                      rtol::Real = (size(A,1)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)),
                      withQ::Bool = true, withZ::Bool = true)

    n = LinearAlgebra.checksquare(A)
    n == LinearAlgebra.checksquare(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
    T = promote_type(eltype(A), eltype(E))
    T <: BlasFloat || (T = promote_type(Float64,T))

    A1 = copy_oftype(A,T)
    E1 = copy_oftype(E,T)

    withQ ? (Q = Matrix{T}(I,n,n)) : (Q = nothing)
    withZ ? (Z = Matrix{T}(I,n,n)) : (Z = nothing)

    ν, blkdims = fisplit!(A1, E1, Q, Z, missing, missing;
                         fast  = fast, finite_infinite = finite_infinite,
                         atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = withQ, withZ = withZ)
    (ilo, ihi) = finite_infinite ? (1, blkdims[1]) : (blkdims[1]+1, n)
    if ilo < ihi
        compq = withQ ? 'V' : 'N'
        compz = withZ ? 'V' : 'N'
        gghrd!(compq, compz, ilo, ihi, A1, E1, Q, Z)
        hgeqz!(compq, compz, ilo, ihi, A1, E1, Q, Z)
    end

    return A1, E1, Q, Z, ν, blkdims
end
"""
    fischursep(A, E; smarg, disc = false, fast = true, finite_infinite = false, stable_unstable = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (At, Et, Q, Z, ν, blkdims)

Reduce the regular matrix pencil `A - λE` to an equivalent block upper triangular form `At - λEt = Q'*(A - λE)*Z`
using orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `At` and `Et`
have separated infinite, stable and unstable eigenvalues with respect to a stability domain `Cs` defined by the stability margin parameter `smarg` and
the stability type parameter `disc`. If `disc = false`, `Cs` is the set of complex numbers with real parts less than `smarg`,
while if `disc = true`, `Cs` is the set of complex numbers with moduli less than `smarg` (i.e., the interior of a disc
of radius `smarg` centered in the origin). If `smarg = missing`, the default value used is `smarg = 0`, if  `disc = false`,
and `smarg = 1`, if `disc = true`.

The pencil `At - λEt` results in one of the following block upper-triangular forms:

(1) if `finite_infinite = false`, then

                   | Ai-λEi   *      *    |
        At - λEt = |    O   A1-λE1   *    |
                   |    0     0    A2-λE2 |

where the `ni x ni` subpencil `Ai-λEi` contains the infinite elementary divisors,
the `n1 x n1` subpencil `A1-λE1` is with the pair `(A1,E1)` in a generalized Schur form, and the
`n2 x n2` subpencil `A2-λE2` is with the pair `(A2,E2)` in a generalized Schur form.
The pencil `A1-λE1` has unstable finite eigenvalues and `A2-λE2` has stable finite eigenvalues if `stable_unstable = false`,
while `A1-λE1` has stable finite eigenvalues and `A2-λE2` has unstable finite eigenvalues if `stable_unstable = true`.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vector `ν` contains the dimensions of the diagonal blocks
of the staircase form  `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[i]-ν[i+1]` for `i = 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i`
(with `ν[nb+1] = 0`).

The dimensions of the diagonal blocks are returned in `blkdims = (ni, n1, n2)`.

(2) if `finite_infinite = true`, then

                   | A1-λE1   *      *    |
        At - λEt = |    O   A2-λE2   *    |
                   |    0     0    Ai-λEi |

where the `ni x ni` subpencil `Ai-λEi` contains the infinite elementary divisors,
the `n1 x n1` subpencil `A1-λE1` is with the pair `(A1,E1)` in a generalized Schur form, and the
`n2 x n2` subpencil `A2-λE2` is with the pair `(A2,E2)` in a generalized Schur form.
The pencil `A1-λE1` has unstable finite eigenvalues and `A2-λE2` has stable finite eigenvalues if `stable_unstable = false`,
while `A1-λE1` has stable finite eigenvalues and `A2-λE2` has unstable finite eigenvalues if `stable_unstable = true`.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular
and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vectors `ν` contains the dimensions of the diagonal blocks
of the staircase form `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[nb-j+1]-ν[nb-j]` for `j = 1, 2, ..., nb` is the number of infinite elementary
divisors of degree `j` (with `ν[0] = 0`).

The dimensions of the diagonal blocks are returned in `blkdims = (n1, n2, ni)`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` and `E`.

The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`.
Otherwise, `Q` is set to `nothing`.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` if `withZ = true`.
Otherwise, `Z` is set to `nothing`.
"""
function fischursep(A::AbstractMatrix, E::AbstractMatrix;
                    smarg::Union{Real,Missing} = missing, disc::Bool = false,
                    fast::Bool = true, finite_infinite::Bool = false, stable_unstable::Bool = false,
                    atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))),
                    rtol::Real = (size(A,1)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)),
                    withQ::Bool = true, withZ::Bool = true)

    n = LinearAlgebra.checksquare(A)
    n == LinearAlgebra.checksquare(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
    disc && (!ismissing(smarg) && smarg < 0 && error("sdeg must be non-negative if disc = true"))

    T = promote_type(eltype(A), eltype(E))
    T <: BlasFloat || (T = promote_type(Float64,T))

    A1 = copy_oftype(A,T)
    E1 = copy_oftype(E,T)

    withQ ? (Q = Matrix{T}(I,n,n)) : (Q = nothing)
    withZ ? (Z = Matrix{T}(I,n,n)) : (Z = nothing)

    ν, blkdims = fisplit!(A1, E1, Q, Z, missing, missing;
                         fast  = fast, finite_infinite = finite_infinite,
                         atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = withQ, withZ = withZ)
    (ilo, ihi) = finite_infinite ? (1, blkdims[1]) : (blkdims[1]+1, n)
    if ilo <= ihi
        compq = withQ ? 'V' : 'N'
        compz = withZ ? 'V' : 'N'
        gghrd!(compq, compz, ilo, ihi, A1, E1, Q, Z)
        _, _, α, β, _, _ = hgeqz!(compq, compz, ilo, ihi, A1, E1, Q, Z)
        i2 = ilo:ihi
        ismissing(smarg) && (smarg = disc ? one(real(T)) : zero(real(T)))
        select2 = disc ? abs.(α[i2]) .< smarg*abs.(β[i2]) : real.(α[i2] ./ β[i2]) .< smarg
        stable_unstable || (select2 = .!select2)

        if finite_infinite
           n3 = blkdims[2]
           n1 = length(select2[select2 .== true])
           n2 = n-n3-n1
           select = [Int.(select2); zeros(Int,n3)]
        else
           n1 = blkdims[1]
           n2 = length(select2[select2 .== true])
           n3 = n-n2-n1
           select = [ones(Int,n1);Int.(select2)]
        end
        tgsen!(withQ, withZ, select, A1, E1, Q, Z)
    else
        (n1,n2,n3) = finite_infinite ? (0,0,blkdims[2]) : (blkdims[1],0,0)
    end

    return A1, E1, Q, Z, ν, (n1,n2,n3)
end
"""
    sfischursep(A, E; smarg, disc = false, fast = true, finite_infinite = false, stable_unstable = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (At, Et, Q, Z, ν, blkdims)

Reduce the regular matrix pencil `A - λE` to an equivalent block upper triangular form `At - λEt = Q'*(A - λE)*Z`
using orthogonal or unitary transformation matrices `Q` and `Z` such that the transformed matrices `At` and `Et`
have separated infinite, stable and unstable eigenvalues with respect to a stability domain `Cs` defined by the stability margin parameter `smarg` and
the stability type parameter `disc`. If `disc = false`, `Cs` is the set of complex numbers with real parts less than `smarg`,
while if `disc = true`, `Cs` is the set of complex numbers with moduli less than `smarg` (i.e., the interior of a disc
of radius `smarg` centered in the origin). If `smarg = missing`, the default value used is `smarg = 0`, if  `disc = false`,
and `smarg = 1`, if `disc = true`.

The pencil `At - λEt` results in one of the following block upper-triangular forms:

(1) if `finite_infinite = true`, then

                   | Ai1   *      *       *     |
        At - λEt = |  0  A1-λE1   *       *     |
                   |  0    0    A2-λE2    *     |
                   |  0    0      0    Ai2-λEi2 |

where the `ni1 x ni1` matrix `Ai1` and the  `ni2 x ni2` subpencil `Ai2-λEi2` contain the infinite elementary divisors,
the `n1 x n1` subpencil `A1-λE1` is with the pair `(A1,E1)` in a generalized Schur form, and the
`n2 x n2` subpencil `A2-λE2` is with the pair `(A2,E2)` in a generalized Schur form.
The pencil `A1-λE1` has unstable finite eigenvalues and `A2-λE2` has stable finite eigenvalues if `stable_unstable = false`,
while `A1-λE1` has stable finite eigenvalues and `A2-λE2` has unstable finite eigenvalues if `stable_unstable = true`.

The subpencil `Ai2-λEi2` is in a staircase form, with `Ai2` nonsingular and upper triangular and `Ei2` nilpotent and upper triangular.
The `nb`-dimensional vector `ν` contains the dimensions of the diagonal blocks
of the staircase form  `Ai2-λEi2` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[nb-j+2]-ν[nb-j+1]` for `j = 1, 2, ..., nb+1` is the number of infinite elementary
divisors of degree `j` (with `ν[0] := 0` and `ν[nb+1] := ni1`).

The dimensions of the diagonal blocks are returned in `blkdims = (ni1, n1, n2, ni2)`.

(2) if `finite_infinite = false`, then

                   | Ai1-λEi1   *      *     *  |
                   |    0     A1-λE1   *     *  |
        At - λEt = |    0       0    A2-λE2  *  |
                   |    0       0      0    Ai2 |

where the `ni1 x ni1` subpencil `Ai1-λEi1` and the `ni2 x ni2` upper triangular and nonsingular matrix `Ai2` contain
the infinite elementary divisors of `A-λE`,
the `n1 x n1` subpencil `A1-λE1` is with the pair `(A1,E1)` in a generalized Schur form, and the
`n2 x n2` subpencil `A2-λE2` is with the pair `(A2,E2)` in a generalized Schur form.
The pencil `A1-λE1` has unstable finite eigenvalues and `A2-λE2` has stable finite eigenvalues if `stable_unstable = false`,
while `A1-λE1` has stable finite eigenvalues and `A2-λE2` has unstable finite eigenvalues if `stable_unstable = true`.

The subpencil `Ai1-λEi1` is in a staircase form, with `Ai1` nonsingular and upper triangular
and `Ei1` nilpotent and upper triangular.
The `nb`-dimensional vectors `ν` contains the dimensions of the diagonal blocks
of the staircase form `Ai1-λEi1` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[i]-ν[i+1]` for `i = 0, 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i`
(with `ν[nb+1] = 0` and `ν[0] = ni2`).

The dimensions of the diagonal blocks are returned in `blkdims = (ni1, n1, n2, ni2)`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` and `E`.

The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

The performed left orthogonal or unitary transformations are accumulated in the matrix `Q` if `withQ = true`.
Otherwise, `Q` is set to `nothing`.
The performed right orthogonal or unitary transformations are accumulated in the matrix `Z` if `withZ = true`.
Otherwise, `Z` is set to `nothing`.
"""
function sfischursep(A::AbstractMatrix, E::AbstractMatrix;
                    smarg::Union{Real,Missing} = missing, disc::Bool = false,
                    fast::Bool = true, finite_infinite::Bool = false, stable_unstable::Bool = false,
                    atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))),
                    rtol::Real = (size(A,1)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)),
                    withQ::Bool = true, withZ::Bool = true)

    n = LinearAlgebra.checksquare(A)
    n == LinearAlgebra.checksquare(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
    disc && (!ismissing(smarg) && smarg < 0 && error("sdeg must be non-negative if disc = true"))

    T = promote_type(eltype(A), eltype(E))
    T <: BlasFloat || (T = promote_type(Float64,T))

    A1 = copy_oftype(A,T)
    E1 = copy_oftype(E,T)

    withQ ? (Q = Matrix{T}(I,n,n)) : (Q = nothing)
    withZ ? (Z = Matrix{T}(I,n,n)) : (Z = nothing)

    ν, blkdims1 = sfisplit!(A1, E1, Q, Z, missing, missing;
                         fast  = fast, finite_infinite = finite_infinite,
                         atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = withQ, withZ = withZ)
    ni1 = blkdims1[1]
    nf  = blkdims1[2]
    ni2 = blkdims1[3]
    (ilo, ihi) = (ni1+1, ni1+nf)
    if ilo <= ihi
        compq = withQ ? 'V' : 'N'
        compz = withZ ? 'V' : 'N'
        gghrd!(compq, compz, ilo, ihi, A1, E1, Q, Z)
        _, _, α, β, _, _ = hgeqz!(compq, compz, ilo, ihi, A1, E1, Q, Z)
        i2 = ilo:ihi
        ismissing(smarg) && (smarg = disc ? one(real(T)) : zero(real(T)))
        select2 = disc ? abs.(α[i2]) .< smarg*abs.(β[i2]) : real.(α[i2] ./ β[i2]) .< smarg
        stable_unstable || (select2 = .!select2)

        select = [ones(Int,ni1); Int.(select2); zeros(Int,ni2) ]
        n1 = length(select2[select2 .== true])
        n2 = nf-n1
        blkdims = (ni1, n1, n2, ni2)
        tgsen!(withQ, withZ, select, A1, E1, Q, Z)
    else
        blkdims = (ni1, 0, 0, ni2)
    end

    return A1, E1, Q, Z, ν, blkdims
end
"""
    fiblkdiag(A, E, B, C; fast = true, finite_infinite = false, trinv = false, atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (At, Et, Bt, Ct, Q, Z, ν, blkdims, sep)

Reduce the regular matrix pencil `A - λE` to an equivalent form `At - λEt = Q*(A - λE)*Z` using
the transformation matrices `Q` and `Z` such that the transformed matrices `At` and `Et` are in one of the
following block diagonal forms:

(1) if `finite_infinite = false`, then

                   | Ai-λEi |   0     |
        At - λEt = |--------|---------|,
                   |    O   | Af-λEf  |

where the `ni x ni` subpencil `Ai-λEi` contains the infinite elementary
divisors and the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular, contains the finite eigenvalues of the pencil `A-λE`.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vector `ν` contains the dimensions of the diagonal blocks
of the staircase form  `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[i]-ν[i+1]` for `i = 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i`
(with `ν[nb+1] = 0`).

The dimensions of the diagonal blocks are returned in `blkdims = (ni, nf)`.

(2) if `finite_infinite = true`, then

                   | Af-λEf |   0    |
        At - λEt = |--------|--------|,
                   |   O    | Ai-λEi |

where the `nf x nf` subpencil `Af-λEf`, with `Ef` nonsingular and upper triangular,
contains the finite eigenvalues of the pencil `A-λE` and the `ni x ni` subpencil `Ai-λEi`
contains the infinite elementary divisors.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular  and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vectors `ν` contains the dimensions of the diagonal blocks
of the staircase form `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[nb-j+1]-ν[nb-j]` for `j = 1, 2, ..., nb` is the number of infinite elementary
divisors of degree `j` (with `ν[0] = 0`).

The dimensions of the diagonal blocks are returned in `blkdims = (nf, ni)`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` and `E`.

The reduction is performed using rank decisions based on rank revealing QR-decomdiagpositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

If `withQ = true`, `Q` contains the left transformation matrix, if `trinv = false`, or its inverse, if `trinv = true`.
If `withQ = false`, `Q` is set to `nothing`.
If `withZ = true`, `Z` contains the right transformation matrix, if `trinv = false`, or its inverse, if `trinv = true`.
If `withZ = false`, `Z` is set to `nothing`.

`Bt = Q*B`, unless `B = missing`, in which case `Bt = missing` is returned, and `Ct = C*Z`,
unless `C = missing`, in which case `Ct = missing` is returned .

An estimation of the separation of the spectra of `Ai-λEi` and `Af-λEf` is returned in `sep`, where  `0 < sep ≤ 1`.
"""
function fiblkdiag(A::AbstractMatrix, E::AbstractMatrix, B::Union{AbstractMatrix,Missing}, C::Union{AbstractMatrix,Missing};
                   fast::Bool = true, finite_infinite::Bool = false, trinv::Bool = false,
                   atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))),
                   rtol::Real = (size(A,1)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)),
                   withQ::Bool = true, withZ::Bool = true)

    compq = withQ || !ismissing(B)
    compz = withZ || !ismissing(C)
    A1, E1, Q, Z, ν, blkdims = fischur(A, E; fast = fast, finite_infinite = finite_infinite, atol1 = atol1, atol2 = atol2,
                                      rtol = rtol, withQ = compq, withZ = compz)
    T = eltype(A1)
    ONER = one(real(T))
    ismissing(B) ? B1 = missing : B1 = Q'*copy_oftype(B,T)
    ismissing(C) ? C1 = missing : C1 = copy_oftype(C,T)*Z

    minimum(blkdims) == 0 && (return A1, E1, B1, C1, withQ ? (trinv ? Q : adjis!(Q)) : nothing,
                              withZ ? (trinv ? adjis!(Z) : Z) : nothing, ν, blkdims, ONER)

    trinv ? (withZ && adjis!(Z)) : (withQ && adjis!(Q))

    n1 = blkdims[1]
    i1 = 1:n1
    i2 = n1+1:size(A,1)
    X = view(E1,i1,i2)
    Y = view(A1,i1,i2)
    _, _, scale = tgsyl!(view(A1,i1,i1), view(A1,i2,i2), Y, view(E1,i1,i1), view(E1,i2,i2), X)

    ONE = one(T)
    ZERO = zero(T)
    scale != 0 && (scale = ONE/scale)

    ismissing(B) || mul!(view(B1,i1,:),X,view(B1,i2,:),scale,ONE)
    ismissing(C) || mul!(view(C1,:,i2),view(C1,:,i1),Y,-scale,ONE)

    if trinv
       withQ && mul!(view(Q,:,i2),view(Q,:,i1),X,-scale,ONE)
       withZ && mul!(view(Z,i1,:),Y,view(Z,i2,:),scale,ONE)
    else
       withQ && mul!(view(Q,i1,:),X,view(Q,i2,:),scale,ONE)
       withZ && mul!(view(Z,:,i2),view(Z,:,i1),Y,-scale,ONE)
    end

    scale == 0 ? sep = ONER : sep = ONER/max(1+norm(X)*abs(scale),1+norm(Y)*abs(scale))

    fill!(X,ZERO)
    fill!(Y,ZERO)

    return A1, E1, B1, C1, Q, Z, ν, blkdims, sep

end
"""
    gsblkdiag(A, E, B, C; smarg, disc = false, fast = true, finite_infinite = false, stable_unstable = false, trinv = false,
              atol1 = 0, atol2 = 0, rtol, withQ = true, withZ = true) -> (At, Et, Bt, Ct, Q, Z, ν, blkdims, sep)

Reduce the regular matrix pencil `A - λE` to an equivalent block diagonal triangular form `At - λEt = Q*(A - λE)*Z`
using the transformation matrices `Q` and `Z` such that the transformed matrices `At` and `Et`
have separated infinite, stable and unstable eigenvalues with respect to a stability domain `Cs`
defined by the stability margin parameter `smarg` and the stability type parameter `disc`.
If `disc = false`, `Cs` is the set of complex numbers with real parts less than `smarg`,
while if `disc = true`, `Cs` is the set of complex numbers with moduli less than `smarg` (i.e., the interior of a disc
of radius `smarg` centered in the origin). If `smarg = missing`, the default value used is `smarg = 0`, if  `disc = false`,
and `smarg = 1`, if `disc = true`.

The pencil `At - λEt` results in one of the following block upper-triangular forms:

(1) if `finite_infinite = false`, then

                   | Ai-λEi   *      0    |
        At - λEt = |    O   A1-λE1   0    |
                   |    0     0    A2-λE2 |

where the `ni x ni` subpencil `Ai-λEi` contains the infinite elementary divisors,
the `n1 x n1` subpencil `A1-λE1` is with the pair `(A1,E1)` in a generalized Schur form, and the
`n2 x n2` subpencil `A2-λE2` is with the pair `(A2,E2)` in a generalized Schur form.
The pencil `A1-λE1` has unstable finite eigenvalues and `A2-λE2` has stable finite eigenvalues if `stable_unstable = false`,
while `A1-λE1` has stable finite eigenvalues and `A2-λE2` has unstable finite eigenvalues if `stable_unstable = true`.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vector `ν` contains the dimensions of the diagonal blocks
of the staircase form  `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[i]-ν[i+1]` for `i = 1, 2, ..., nb` is the number of infinite elementary divisors of degree `i`
(with `ν[nb+1] = 0`).

The dimensions of the diagonal blocks are returned in `blkdims = (ni, n1, n2)`.

(2) if `finite_infinite = true`, then

                   | A1-λE1   0      0    |
        At - λEt = |    O   A2-λE2   *    |
                   |    0     0    Ai-λEi |

where the `ni x ni` subpencil `Ai-λEi` contains the infinite elementary divisors,
the `n1 x n1` subpencil `A1-λE1` is with the pair `(A1,E1)` in a generalized Schur form, and the
`n2 x n2` subpencil `A2-λE2` is with the pair `(A2,E2)` in a generalized Schur form.
The pencil `A1-λE1` has unstable finite eigenvalues and `A2-λE2` has stable finite eigenvalues if `stable_unstable = false`,
while `A1-λE1` has stable finite eigenvalues and `A2-λE2` has unstable finite eigenvalues if `stable_unstable = true`.

The subpencil `Ai-λEi` is in a staircase form, with `Ai` nonsingular and upper triangular
and `Ei` nilpotent and upper triangular.
The `nb`-dimensional vectors `ν` contains the dimensions of the diagonal blocks
of the staircase form `Ai-λEi` such that `i`-th block has dimensions `ν[i] x ν[i]`.
The difference `ν[nb-j+1]-ν[nb-j]` for `j = 1, 2, ..., nb` is the number of infinite elementary
divisors of degree `j` (with `ν[0] = 0`).

The dimensions of the diagonal blocks are returned in `blkdims = (n1, n2, ni)`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` and `E`.

The reduction is performed using rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.

If `withQ = true`, `Q` contains the left transformation matrix, if `trinv = false`, or its inverse, if `trinv = true`.
If `withQ = false`, `Q` is set to `nothing`.
If `withZ = true`, `Z` contains the left transformation matrix, if `trinv = false`, or its inverse, if `trinv = true`.
If `withZ = false`, `Z` is set to `nothing`.

`Bt = Q*B`, unless `B = missing`, in which case `Bt = missing` is returned, and `Ct = C*Z`,
unless `C = missing`, in which case `Ct = missing` is returned .

An estimation of the separation of the spectra of the two underlying diagonal blocks is returned in `sep`,
where  `0 ≤ sep ≤ 1`. A value `sep ≈ 0` indicates that `A1-λE1` and `A2-λE2` have some almost equal eigenvalues.
"""
function gsblkdiag(A::AbstractMatrix, E::AbstractMatrix, B::Union{AbstractMatrix,Missing}, C::Union{AbstractMatrix,Missing};
                   smarg::Union{Real,Missing} = missing, disc::Bool = false, trinv::Bool = false,
                   fast::Bool = true, finite_infinite::Bool = false, stable_unstable::Bool = false,
                   atol1::Real = zero(real(eltype(A))), atol2::Real = zero(real(eltype(E))),
                   rtol::Real = (size(A,1)*eps(real(float(one(eltype(A))))))*iszero(min(atol1,atol2)),
                   withQ::Bool = true, withZ::Bool = true)

    compq = withQ || !ismissing(B)
    compz = withZ || !ismissing(C)
    A1, E1, Q, Z, ν, blkdims = fischursep(A, E; smarg = smarg, disc = disc,
                                      fast = fast, finite_infinite = finite_infinite, stable_unstable = stable_unstable,
                                      atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = compq, withZ = compz)
    T = eltype(A1)
    ONER = one(real(T))
    ismissing(B) ? B1 = missing : B1 = Q'*copy_oftype(B,T)
    ismissing(C) ? C1 = missing : C1 = copy_oftype(C,T)*Z

    finite_infinite ? n1 = blkdims[1] : n1 = blkdims[1] + blkdims[2]
    (n1 == 0 || size(A1,1) == n1) && (return A1, E1, B1, C1, withQ ? (trinv ? Q : adjis!(Q)) : nothing,
               withZ ? (trinv ? adjis!(Z) : Z) : nothing, ν, blkdims,ONER)

    trinv ? (withZ && adjis!(Z)) : (withQ && adjis!(Q))

    i1 = 1:n1
    i2 = n1+1:size(A,1)
    X = view(E1,i1,i2)
    Y = view(A1,i1,i2)
    _, _, scale = tgsyl!(view(A1,i1,i1), view(A1,i2,i2), Y, view(E1,i1,i1), view(E1,i2,i2), X)

    ONE = one(T)
    ZERO = zero(T)
    scale != 0 && (scale = ONE/scale)

    ismissing(B) || mul!(view(B1,i1,:),X,view(B1,i2,:),scale,ONE)
    ismissing(C) || mul!(view(C1,:,i2),view(C1,:,i1),Y,-scale,ONE)

    if trinv
       withQ && mul!(view(Q,:,i2),view(Q,:,i1),X,-scale,ONE)
       withZ && mul!(view(Z,i1,:),Y,view(Z,i2,:),scale,ONE)
    else
       withQ && mul!(view(Q,i1,:),X,view(Q,i2,:),scale,ONE)
       withZ && mul!(view(Z,:,i2),view(Z,:,i1),Y,-scale,ONE)
    end

    scale == 0 ? sep = ONER : sep = ONER/max(1+norm(X)*abs(scale),1+norm(Y)*abs(scale))

    fill!(X,ZERO)
    fill!(Y,ZERO)

    return A1, E1, B1, C1, Q, Z, ν, blkdims, sep

end
@inline function adjis!(A::AbstractMatrix)
    # in-situ adjoint of a square matrix
    n = LinearAlgebra.checksquare(A)
    for i = 1:n
        for j = i:n
            t = A[i,j]'
            A[i,j] = A[j,i]'
            A[j,i] = t
        end
    end
    return A
end
"""
    ssblkdiag(A, B, C; smarg, disc = false, stable_unstable = false, withQ = true, withZ = true) -> (At, Bt, Ct, Q, Z, blkdims, sep)

Reduce the regular matrix pencil `A - λI` to an equivalent block diagonal triangular form `At - λI = Q*(A - λI)*Z`
using the transformation matrices `Q` and `Z`, where `Q = inv(Z)`, such that the transformed matrix `At` have
separated stable and unstable eigenvalues with respect to a stability domain `Cs` defined by the
stability margin parameter `smarg` and the stability type parameter `disc`.
If `disc = false`, `Cs` is the set of complex numbers with real parts less than `smarg`,
while if `disc = true`, `Cs` is the set of complex numbers with moduli less than `smarg` (i.e., the interior of a disc
of radius `smarg` centered in the origin). If `smarg = missing`, the default value used is `smarg = 0`, if  `disc = false`,
and `smarg = 1`, if `disc = true`.
The matrix `At` results in the following block diagonal form

        At = | A1  0  |
             | 0   A2 |

where the `n1 x n1` matrix `A1` and the `n2 x n2` matrix `A2` are in Schur form.
The matrix `A1` has unstable eigenvalues and `A2` has stable eigenvalues if `stable_unstable = false`,
while `A1` has stable eigenvalues and `A2` has unstable eigenvalues if `stable_unstable = true`.
The dimensions of the diagonal blocks are returned in `blkdims = (n1, n2)`.
If `withQ = true`, `Q` contains the left transformation matrix. If `withQ = false`, `Q` is set to `nothing`.
If `withZ = true`, `Z` contains the right transformation matrix. If `withZ = false`, `Z` is set to `nothing`.
`Bt = Q*B`, unless `B = missing`, in which case `Bt = missing` is returned, and `Ct = C*Z`,
unless `C = missing`, in which case `Ct = missing` is returned .
An estimation of the separation of the spectra of the two underlying diagonal blocks is returned in `sep`,
where  `0 ≤ sep ≤ 1`. A value `sep ≈ 0` indicates that `A1` and `A2` have some almost equal eigenvalues.
"""
function ssblkdiag(A::AbstractMatrix{T}, B::Union{AbstractMatrix{T},Missing}, C::Union{AbstractMatrix{T},Missing};
                   smarg::Union{Real,Missing} = missing, disc::Bool = false, stable_unstable::Bool = false,
                   withQ::Bool = true, withZ::Bool = true ) where T <: BlasFloat

    n = LinearAlgebra.checksquare(A)
    S = schur(A)
    ONER = real(float(one(T)))
    if disc
       ismissing(smarg) && (smarg = ONER)
       stable_unstable ? select = abs.(S.values) .< smarg : select = abs.(S.values) .> smarg
    else
       ismissing(smarg) && (smarg = real(float(zero(T))))
       stable_unstable ? select = real(S.values) .< smarg : select = real(S.values) .> smarg
    end
    n1 = length(filter(y-> y == true,select))

    (n1 == 0 || n1 == n) && (return A, B, C, withQ ? copy(S.Z') : nothing, withZ ? S.Z : nothing, (n1,n-n1), ONER)

    ordschur!(S, select)

    i1 = 1:n1
    i2 = n1+1:size(A,1)
    Y = view(S.T,i1,i2)  # Y will contain the negative of solution -Y
    _, scale = LAPACK.trsyl!('N', 'N', view(S.T,i1,i1), view(S.T,i2,i2), Y, -1)

    ONE = one(T)
    ZERO = zero(T)
    scale != 0 && (scale = ONE/scale)

    ismissing(B) ? B1 = missing : (B1 = S.Z'*B; mul!(view(B1,i1,:),Y,view(B1,i2,:),scale,ONE))
    ismissing(C) ? C1 = missing : (C1 = C*S.Z; mul!(view(C1,:,i2),view(C1,:,i1),Y,-scale,ONE))

    withQ ? Q = copy(S.Z') : Q = nothing
    withQ && mul!(view(Q,i1,:),Y,view(Q,i2,:),scale,ONE)
    withZ && mul!(view(S.Z,:,i2),view(S.Z,:,i1),Y,-scale,ONE)

    scale == 0 ? sep = zero(real(T)) : sep = ONER/(1+norm(Y)*abs(scale))

    fill!(Y,ZERO)
    return S.T, B1, C1, Q, withZ ? S.Z : nothing, (n1, n-n1), sep
end
# fallback for versions prior 1.3
VERSION > v"1.3.0" || (mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}, α::T, β::T) where {T<:BlasFloat} =
                           BLAS.gemm!('N', 'N', α, A, B, β, C))


text:
using CSPSolver
using Random
using Test

@testset "sudoku easy 1" begin
    board = [
         3 -1 -1 -1 -1 -1 -1  6  8;
         7  6  2  1 -1 -1  3  4 -1;
        -1  4 -1 -1 -1 -1 -1 -1 -1;
        -1 -1  9 -1 -1  3 -1  5  7;
        -1 -1  5  2  6 -1 -1  8 -1;
        -1  8  7 -1  1 -1 -1 -1 -1;
         8 -1 -1  6  7  1 -1 -1 -1;
        -1  2 -1 -1  8  4  5  7  3;
         9 -1 -1  3  5  2  8  1 -1
    ]

    s = Sudoku(board)
    res = backtrack(s, 1)
    @test !isnothing(res)
end

@testset "sudoku hard 1" begin
    board = [
        -1 -1  5 -1 -1  2 -1  4 -1;
        -1  7 -1 -1  4 -1 -1  9  8;
        -1 -1  6 -1  3 -1 -1  7 -1;
         4 -1 -1  3 -1 -1 -1 -1  1;
        -1  9 -1 -1 -1 -1 -1  2 -1;
        -1 -1 -1 -1 -1 -1  6 -1 -1;
        -1 -1  3  9  7 -1 -1 -1 -1;
         5 -1 -1 -1 -1 -1  2 -1 -1;
        -1 -1 -1  1 -1 -1 -1 -1 -1
    ]

    s = Sudoku(board)
    res = backtrack(s, 1)
    @test !isnothing(res)
end

@testset "sudoku generate" begin
    Random.seed!(1)
    board, solution = generate_random_sudoku(1)
    s = Sudoku(board)
    solutions = backtrack(s, 2)
    # even though we asked for 2 solutions, there should only be 1
    @test length(solutions) == 1
    # the one solution found must be the solution provided
    @test solutions[1].variables == solution
end

ext: jl

max_stars_repo_name: gindij/CSPSolver

max_issues_repo_name: gindij/CSPSolver

max_forks_repo_name: gindij/CSPSolver


text:
"""
# Rolling Functions

```julia
rollapply(fun::Function, ts::TS, column::Any, windowsize::Int)
```

Apply a function to a column of `ts` for each continuous set of rows
of size `windowsize`. `column` could be any of the `DataFrame` column
selectors.

The output is a TS object with `(nrow(ts) - windowsize + 1)` rows
indexed with the last index value of each window.

This method uses `RollingFunctions` package to implement this
functionality.

# Examples

```jldoctest; setup = :(using TSx, DataFrames, Dates, Random, Statistics)
julia> ts = TS(1:12, Date("2022-02-01"):Month(1):Date("2022-02-01")+Month(11))

julia> show(ts)
(12 x 1) TS with Dates.Date Index

 Index       x1
 Date        Int64
───────────────────
 2022-02-01      1
 2022-03-01      2
 2022-04-01      3
 2022-05-01      4
 2022-06-01      5
 2022-07-01      6
 2022-08-01      7
 2022-09-01      8
 2022-10-01      9
 2022-11-01     10
 2022-12-01     11
 2023-01-01     12

julia> rollapply(sum, ts, :x1, 10)
(3 x 1) TS with Dates.Date Index

 Index       x1_rolling_sum
 Date        Float64
────────────────────────────
 2022-11-01            55.0
 2022-12-01            65.0
 2023-01-01            75.0

julia> rollapply(Statistics.mean, ts, 1, 5)
(8 x 1) TS with Dates.Date Index

 Index       x1_rolling_mean
 Date        Float64
─────────────────────────────
 2022-06-01              3.0
 2022-07-01              4.0
 2022-08-01              5.0
 2022-09-01              6.0
 2022-10-01              7.0
 2022-11-01              8.0
 2022-12-01              9.0
 2023-01-01             10.0

```
"""
function rollapply(fun::Function, ts::TS, column::Any, windowsize::Int) # TODO: multiple columns
    if windowsize < 1
        error("windowsize must be greater than or equal to 1")
    end
    col = Int(1)
    if typeof(column) <: Int
        col = copy(column)
        col = col+1             # index is always 1
    else
        col = column
    end
    res = RollingFunctions.rolling(fun, ts.coredata[!, col], windowsize)
    idx = TSx.index(ts)[windowsize:end]
    colname = names(ts.coredata[!, [col]])[1]
    res_df = DataFrame([idx, res], ["Index", "$(colname)_rolling_$(fun)"])
    return TS(res_df)
end


text:
"""
    water_density(temp)

Calculate water density as a function of temperature (`temp`).

# Parameters

`temp`: Temperature in Kelvin.

# Returns

`rho_w`: Water density (kg m^-3).

# References

[WP02] Wagner, W. and Pruß, A. (2002). The IAPWS Formulation 1995 for the
    Thermodynamic Properties of Ordinary Water Substance for General and
    Scientific Use. *J. Phys. Chem. Ref. Data*, 31, 387--535.
    <https://doi.org/10.1063/1.1461829>.

# Examples

```jldoctest
julia> @printf "%g" water_density(298.15)
996.999
```
"""
function water_density(temp)
    # critical temperature and density
    T_crit = 647.096
    rho_crit = 322.

    theta = 1. - temp / T_crit
    rho_ratio = 1. + 1.99274064 * theta ^ (1. / 3.) +
        1.09965342 * theta ^ (2. / 3.) - 0.510839303 * theta ^ (5. / 3.) -
        1.7549349 * theta ^ (16. / 3.) - 45.5170352 * theta ^ (43. / 3.) -
        6.74694450e5 * theta ^ (110. / 3.)

    return rho_crit * rho_ratio
end

"""
    water_dissoc(temp)

Calculate water dissociation constant (pK_w) as a function of temperature.

# Parameters

`temp`: Temperature in Kelvin.

# Returns

`pK_w`:  Water dissociation constant.

# References

[BL06] Bandura, A. V. and Lvov, S. N. (2006). The Ionization Constant of Water
    over Wide Ranges of Temperature and Density. *J. Phys. Chem. Ref. Data*,
    35, 15--30. <https://doi.org/10.1063/1.1928231>.

# Examples

```jldoctest
julia> @printf "%g" water_dissoc(298.15)
13.9949
```
"""
function water_dissoc(temp)
    n = 6.
    alpha_0 = -0.864671
    alpha_1 = 8659.19
    alpha_2 = -22786.2
    beta_0 = 0.642044
    beta_1 = -56.8534
    beta_2 = -0.375754

    rho_w = water_density(temp) * 1e-3  # convert to g cm^-3 here
    temp_sq = temp * temp

    Z = rho_w * exp(alpha_0 + alpha_1 / temp +
                    alpha_2 / temp_sq * rho_w ^ (2. / 3.))
    pK_w_G = 0.61415 + 48251.33 / temp - 67707.93 / temp_sq +
        10102100. / (temp * temp_sq)

    return -2. * n * (log10(1. + Z) - Z / (Z + 1.) * rho_w *
                      (beta_0 + beta_1 / temp + beta_2 * rho_w)) +
        pK_w_G + 2. * log10(Constants.M_w)
end


text:
"""
        atkinson94(setting, iters, crit)

Runs the Atkinson94 algorithm to find out outliers using LMS method.

# Arguments
- `setting::RegressionSetting`: A regression setting object.
- `iters::Int`: Number of random samples.
- `crit::Float64`: Critical value for residuals

# Description
The algorithm randomly selects initial basic subsets and performs a very robust method, e.g `lms`
to enlarge the basic subset. In each iteration of forward search, the best objective value and parameter
estimates are stored. These values are also used in Atkinson's Stalactite Plot for a visual investigation of
outliers. See `atkinsonstalactiteplot`.

# Output
- `["optimum_index"]`: The iteration number in which the minimum objective is obtained
- `["residuals_matrix"]`: Matrix of residuals obtained in each iteration
- `["outliers"]`: Array of indices of detected outliers
- `["objective"]`: Minimum objective value
- `["coef"]`: Estimated regression coefficients
- `["crit"]`: Critical value given by the user.

# Examples
```julia-repl
julia> reg = createRegressionSetting(@formula(stackloss ~ airflow + watertemp + acidcond), stackloss)
julia> atkinson94(reg)
Dict{Any,Any} with 6 entries:
  "optimum_index"    => 10
  "residuals_matrix" => [0.0286208 0.0620609 … 0.0796249 0.0; 0.0397778 0.120547 … 0.118437 0.0397778; … ; 1.21133 1.80846 … 0.690327 4.14366; 1.61977 0.971592 … 0.616204 3.58098]
  "outliers"         => [1, 3, 4, 21]
  "objective"        => 0.799134
  "coef"             => [-38.3133, 0.745659, 0.432794, 0.0104587]
  "crit"             => 3.0

```
# References
Atkinson, Anthony C. "Fast very robust methods for the detection of multiple outliers."
Journal of the American Statistical Association 89.428 (1994): 1329-1339.
"""
function atkinson94(setting::RegressionSetting; iters=nothing, crit=3.0)
    X, y = @extractRegressionSetting setting
    return atkinson94(X, y, iters=iters, crit=crit)
end

function atkinson94(X::Array{Float64,2}, y::Array{Float64,1}; iters=nothing, crit=3.0)
    n, p = size(X)

    # the median index
    h = Int(floor((n + p + 1) / 2.0))

    if iters === nothing
        iters = minimum([500 * p, 3000])
    end

    bestobjective = Inf
    bestparameters = zeros(Float64, p)
    bestres = zeros(Float64, n, n)
    bestindex = 0
    indices = collect(1:n)

    # store all the sigma values across all iterations
    sigmas = zeros(Float64, iters, n)

    for iter in 1:iters
        m_subset_indices = sample(indices, p, replace=false)

        # stores the n - p sigma values of a forward run
        sigma_tilde = zeros(Float64, n)

        # stores the (n - p) * n residuals during the forward run
        studentized_residuals = zeros(Float64, n, n)
        copy_parameters = false

        @inbounds for m = p:n
            olsreg = ols(X[m_subset_indices,:], y[m_subset_indices])
            betas = coef(olsreg)
            e = (y .- X * betas)
            r = e.^2

            # sigma is the median of the residuals
            sigma_tilde[m] = sqrt(sort(r)[h])
            XXinv = pinv(X[m_subset_indices, :]'X[m_subset_indices, :])

            # scale the residual according to whether it belongs to m_subset_indices or not
            @inbounds for index in 1:n
                if index in m_subset_indices
                    h_i = X[index, :]' * XXinv * X[index, :]
                    studentized_residuals[m, index] = abs(e[index]) / (sigma_tilde[m] * sqrt(abs(1 - h_i)))
                else
                    d_i = X[index, :]' * XXinv * X[index, :]
                    studentized_residuals[m, index] = abs(e[index]) / (sigma_tilde[m] * sqrt(1 + d_i))
                    r[index] = (e[index]^2) / (1 + d_i)
                end
            end

            # increase the size of subset by 1
            ordered_residue_indices = sortperm(r)
            if m != n
                m_subset_indices = ordered_residue_indices[1:m + 1]
            end

            # if sigma was lowest, record the parameters
            if sigma_tilde[m] < bestobjective
                bestobjective = sigma_tilde[m]
                bestparameters = betas
                bestindex = m
                copy_parameters = true
        end
        end
#         sigmas[iter, :] = sigma_tilde

        # copy the entire n^2 residuals for the forward search to generate stalactite plot
        if copy_parameters
            bestres = copy(studentized_residuals)
        end
    end
#     sigma_bar = mean(sigmas, dims=1)
#     for m = p:n
#         bestres[m, :] = bestres[m, :] .* sigma_bar[m]
#     end
    d = Dict()
#     d["sigma_bar"] = sigma_bar
#     d["sigmas"] = sigmas[:, p:end]
    d["coef"] = bestparameters
    d["objective"] = bestobjective
    d["optimum_index"] = bestindex
    d["residuals_matrix"] = bestres[p:end, :]
    d["crit"] = crit
    d["outliers"] = filter(i -> bestres[bestindex, i] > crit, 1:n)
    return d
end

max_stars_repo_name: jbytecode/LinRegOutliers

max_issues_repo_name: jbytecode/LinRegOutliers

max_forks_repo_name: jbytecode/LinRegOutliers


text:
#using advanced_timer
using Main.advanced_timer
using Test
using LinearAlgebra

# write your own tests here
begin
  timer =  Main.advanced_timer.class_advanced_timer()
  Main.advanced_timer.start_advanced_timer(timer);
  Main.advanced_timer.start_advanced_timer(timer,"test");
  Main.advanced_timer.pause_advanced_timer(timer,"test");
  Main.advanced_timer.start_advanced_timer(timer,"takes time");
  LinearAlgebra.dot(rand(1000000),rand(1000000))
  Main.advanced_timer.pause_advanced_timer(timer,"takes time");
  Main.advanced_timer.start_advanced_timer(timer,"abc");
  LinearAlgebra.dot(randn(100),randn(100))
  Main.advanced_timer.pause_advanced_timer(timer,"abc");
  Main.advanced_timer.pause_advanced_timer(timer);
  Main.advanced_timer.print_timer_stats(timer)

  timer2 = Main.advanced_timer.class_advanced_timer()
  Main.advanced_timer.start_advanced_timer(timer2);
  Main.advanced_timer.start_advanced_timer(timer2,"test");
  Main.advanced_timer.pause_advanced_timer(timer2,"test");
  Main.advanced_timer.start_advanced_timer(timer2,"takes time");
  LinearAlgebra.dot(rand(1000000),rand(1000000))
  Main.advanced_timer.pause_advanced_timer(timer2,"takes time");
  Main.advanced_timer.pause_advanced_timer(timer2);
  Main.advanced_timer.print_timer_stats(timer2)

  merged_timers = Main.advanced_timer.merge_timers(timer,timer2)
  Main.advanced_timer.print_timer_stats(merged_timers)
end


text:
@testset "Iteration" begin
    t = SegmentMap{Int,Int}()
    @test isempty([x for x in t])

    @testset "from" begin
        n = 100
        global maxend = 1000000
        segments = [randsegment(1, maxend) for _ in 1:n]
        startpos = 50000
        expected_count = 0
        t = SegmentMap{Int,Int}()
        for (i, segment) in enumerate(segments)
            t[segment] = i
            if segment[2] >= startpos
                expected_count += 1
            end
        end

        @test length(collect(from(t, startpos))) === expected_count
        @test length(collect(from(t, 0))) === n
        @test length(collect(from(t, maxend + 1))) === 0
        @test length(collect(from(SegmentTree{Int,Int}(), 0))) === 0
    end
end


text:
using NARX
using Test

@testset "NARX.jl" begin
    # Write your tests here.
end

max_stars_repo_name: wmkouw/NARX

max_issues_repo_name: wmkouw/NARX

max_forks_repo_name: wmkouw/NARX


text:
# log2(x) =  vgetexp(x) + log2(vgetmant8(x))
@inline function vlog_v1(x1::VectorizationBase.AbstractSIMD{W,Float64}) where {W} # Testing if an assorted mix of operations
  x14 = vgetmant(x1)
  x8 = vgetexp(1.3333333333333333 * x1)
  # x2 = reinterpret(UInt64, x1)
  # x3 = x2 >>> 0x0000000000000020
  # greater_than_zero = x1 > zero(x1)
  # alternative = VectorizationBase.ifelse(x1 == zero(x1), -Inf, NaN)
  # isinf = x1 == Inf
  # x5 = x3 + 0x0000000000095f62
  # x6 = x5 >>> 0x0000000000000014
  # x7 = x6 - 0x00000000000003ff
  # x8 = convert(Float64, x7 % Int)
  # @show x3 x5 x6 x7 x8
  # x9 = x5 << 0x0000000000000020
  # x10 = x9 & 0x000fffff00000000
  # x11 = x10 + 0x3fe6a09e00000000
  # x12 = x2 & 0x00000000ffffffff
  # x13 = x11 | x12
  # x14 = reinterpret(Float64, x13)
  x15 = x14 - 1.0
  x18 = x14 + 1.0
  x16 = x15 * x15
  x17 = 0.5 * x16
  x19 = x15 / x18
  x20 = x19 * x19
  x21 = x20 * x20
  x22 = vfmadd(x21, 0.15313837699209373, 0.22222198432149784)
  x23 = vfmadd(x21, x22, 0.3999999999940942)
  x24 = x23 * x21
  x25 = vfmadd(x21, 0.14798198605116586, 0.1818357216161805)
  x29 = x24 + x17
  x26 = vfmadd(x21, x25, 0.2857142874366239)
  x27 = vfmadd(x21, x26, 0.6666666666666735)
  x28 = x27 * x20
  x30 = x29 + x28
  x31 = x8 * 1.9082149292705877e-10
  x32 = vfmadd(x19, x30, x31)
  x33 = x15 - x17
  x34 = x33 + x32
  x35 = vfmadd(x8, 0.6931471803691238, x34)
  return x35
  # @show x35
  # x36 = VectorizationBase.ifelse(greater_than_zero, x35, alternative)
  # VectorizationBase.ifelse(isinf, Inf, x36)
end

const LOGTABLE = [
  0.0,
  0.007782140442054949,
  0.015504186535965254,
  0.02316705928153438,
  0.030771658666753687,
  0.0383188643021366,
  0.0458095360312942,
  0.053244514518812285,
  0.06062462181643484,
  0.06795066190850775,
  0.07522342123758753,
  0.08244366921107459,
  0.08961215868968714,
  0.09672962645855111,
  0.10379679368164356,
  0.11081436634029011,
  0.11778303565638346,
  0.12470347850095724,
  0.13157635778871926,
  0.13840232285911913,
  0.1451820098444979,
  0.15191604202584197,
  0.15860503017663857,
  0.16524957289530717,
  0.17185025692665923,
  0.1784076574728183,
  0.184922338494012,
  0.19139485299962947,
  0.19782574332991987,
  0.2042155414286909,
  0.21056476910734964,
  0.21687393830061436,
  0.22314355131420976,
  0.22937410106484582,
  0.2355660713127669,
  0.24171993688714516,
  0.24783616390458127,
  0.25391520998096345,
  0.25995752443692605,
  0.26596354849713794,
  0.27193371548364176,
  0.2778684510034563,
  0.2837681731306446,
  0.28963329258304266,
  0.2954642128938359,
  0.3012613305781618,
  0.3070250352949119,
  0.3127557100038969,
  0.3184537311185346,
  0.324119468654212,
  0.329753286372468,
  0.3353555419211378,
  0.3409265869705932,
  0.34646676734620857,
  0.3519764231571782,
  0.3574558889218038,
  0.3629054936893685,
  0.3683255611587076,
  0.37371640979358406,
  0.37907835293496944,
  0.38441169891033206,
  0.3897167511400252,
  0.394993808240869,
  0.4002431641270127,
  0.4054651081081644,
  0.4106599249852684,
  0.415827895143711,
  0.42096929464412963,
  0.4260843953109001,
  0.4311734648183713,
  0.43623676677491807,
  0.4412745608048752,
  0.44628710262841953,
  0.45127464413945856,
  0.4562374334815876,
  0.46117571512217015,
  0.46608972992459924,
  0.470979715218791,
  0.4758459048699639,
  0.4806885293457519,
  0.4855078157817008,
  0.4903039880451938,
  0.4950772667978515,
  0.4998278695564493,
  0.5045560107523953,
  0.5092619017898079,
  0.5139457511022343,
  0.5186077642080457,
  0.5232481437645479,
  0.5278670896208424,
  0.5324647988694718,
  0.5370414658968836,
  0.5415972824327444,
  0.5461324375981357,
  0.5506471179526623,
  0.5551415075405016,
  0.5596157879354227,
  0.564070138284803,
  0.5685047353526688,
  0.5729197535617855,
  0.5773153650348236,
  0.5816917396346225,
  0.5860490450035782,
  0.5903874466021763,
  0.5947071077466928,
  0.5990081896460834,
  0.6032908514380843,
  0.6075552502245418,
  0.6118015411059929,
  0.616029877215514,
  0.6202404097518576,
  0.6244332880118935,
  0.6286086594223741,
  0.6327666695710378,
  0.6369074622370692,
  0.6410311794209312,
  0.6451379613735847,
  0.6492279466251099,
  0.6533012720127457,
  0.65735807270836,
  0.661398482245365,
  0.6654226325450905,
  0.6694306539426292,
  0.6734226752121667,
  0.6773988235918061,
  0.6813592248079031,
  0.6853040030989194,
  0.689233281238809,
  0.6931471805599453,
];

logb(::Type{Float32}, ::Val{2}) = 1.4426950408889634
logb(::Type{Float32}, ::Val{:ℯ}) = One()
logb(::Type{Float32}, ::Val{10}) = 0.4342944819032518
logbU(::Type{Float64}, ::Val{2}) = 1.4426950408889634
logbL(::Type{Float64}, ::Val{2}) = 2.0355273740931033e-17
logbU(::Type{Float64}, ::Val{:ℯ}) = One()
logbL(::Type{Float64}, ::Val{:ℯ}) = Zero()
logbU(::Type{Float64}, ::Val{10}) = 0.4342944819032518
logbL(::Type{Float64}, ::Val{10}) = 1.098319650216765e-17

@inline function vlog_base(v::AbstractSIMD{W,Float64}) where {W}
  y = vgetmant12(v)
  mf = vgetexp(v)

  y128 = 128.0 * y
  f128 = vsreduce(y128, Val(0))
  F128 = y128 - f128# - 128.0
  jp = convert(UInt, F128) - 0x0000000000000080 # - 128
  # jp = convert(UInt,F128 - 128.0) # - 128
  # jp = convert(UInt,F128) - 0x000000000000007f # - 127

  hi = vload(zstridedpointer(LOGTABLE), (jp,))
  # l_hi = muladd(0.6931471805601177, mf, hi)
  l_hi = muladd(0.6931471805599453, mf, hi)

  u = (2.0 * f128) / (y128 + F128)
  # u = (2.0*f128)*vinv_fast(y128+F128)
  v = u * u
  q = u * v * muladd(0.012500053168098584, v, 0.08333333333303913)

  ## Step 4
  m_hi = logbU(Float64, Val(:ℯ))
  m_lo = logbL(Float64, Val(:ℯ))
  return fma(m_hi, l_hi, fma(m_hi, (u + q), m_lo * l_hi))
end

const LOG_TABLE_1 = Vec((
  Core.VecElement(-0.6931471805599453),
  Core.VecElement(-0.6325225587435105),
  Core.VecElement(-0.5753641449035618),
  Core.VecElement(-0.5212969236332861),
  Core.VecElement(-0.4700036292457356),
  Core.VecElement(-0.42121346507630353),
  Core.VecElement(-0.3746934494414107),
  Core.VecElement(-0.33024168687057687),
))
const LOG_TABLE_2 = Vec((
  Core.VecElement(-0.2876820724517809),
  Core.VecElement(-0.24686007793152578),
  Core.VecElement(-0.2076393647782445),
  Core.VecElement(-0.16989903679539747),
  Core.VecElement(-0.13353139262452263),
  Core.VecElement(-0.09844007281325252),
  Core.VecElement(-0.06453852113757118),
  Core.VecElement(-0.0317486983145803),
))
const LOG2_TABLE_1 = Vec((
  Core.VecElement(-1.0),
  Core.VecElement(-0.9125371587496606),
  Core.VecElement(-0.8300749985576876),
  Core.VecElement(-0.7520724865564145),
  Core.VecElement(-0.6780719051126377),
  Core.VecElement(-0.6076825772212398),
  Core.VecElement(-0.5405683813627028),
  Core.VecElement(-0.4764380439429871),
))
const LOG2_TABLE_2 = Vec((
  Core.VecElement(-0.4150374992788438),
  Core.VecElement(-0.3561438102252753),
  Core.VecElement(-0.2995602818589078),
  Core.VecElement(-0.24511249783653147),
  Core.VecElement(-0.19264507794239588),
  Core.VecElement(-0.14201900487242788),
  Core.VecElement(-0.09310940439148147),
  Core.VecElement(-0.04580368961312479),
))

@inline function logkern_5(x)
  # c5 = 1.442695040910480151153610022148007414933621426294604120324193438507895065110221
  # c4 = -0.72134752044717115185862113438227555323895953631804827391698847566536447158971
  # c3 = 0.4808981720904872241142144122992626169360496403735493070672998980156544684558863
  # c2 = -0.360673694648970133332498335009073996964643657971752173300330580172269152150781
  # c1 = 0.288894281140228357047827707854224408808934973997456175085011693761970856737952
  # c0 = -0.2406712024835816804281715317202731627491220621210104128295077305793044245650939
  # c5 = 1.000000000004163234554825918991486005405540498226964649982446185635197444806964
  # c4 = -0.5000000000218648309257667907050168790946684146132007063399540958947352513227305
  # c3 = 0.3333332566080562266977844691232384812586343066695585749657855454332169568535108
  # c2 = -0.2499998582263749067900056202432438929369945607183640079498376650719435347943278
  # c1 = 0.2002094490921391961924703669397489236496809044839667749701044680362284866088264
  # c0 = -0.1669110816794161712236614512036245351950242871767130610770379611369082994542161
  c5 = 1.442695040894969685206229605569847427762291596137745860941223179497930808392396
  c4 = -0.7213475204760259868264417322489103259481854442059435485772582433076339753304261
  c3 = 0.4808982362718110098786870255779326384270826960362097549483653590161925169305216
  c2 = -0.3606735556861350010148197635190486357722791582949450655689225850005157925747294
  c1 = 0.2888411793443405953870653225945083997249994212173033378047806528684724370915621
  c0 = -0.2408017898083064242160316054616972180231905999993934707245482866675760802361285
  x * vmuladd_fast(
    vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(c0, x, c1), x, c2), x, c3), x, c4),
    x,
    c5,
  )
end
@inline function logkern_8(x)
  c8 = 1.442695040888963407816086035203962841602383176205260263758551036257573600888366
  c7 = -0.7213475204444770815507673261267005370242876005369959065788590988093119140411675
  c6 = 0.4808983469629686528945376971132664305807537133089662498292066034260696561965772
  c5 = -0.360673760285340441741375667112629490399269780822657391947176916709523570522447
  c4 = 0.2885390083116887609918806358733183651489107200432558827730599551390737454696275
  c3 = -0.2404489409194524614303773871505793252012695382714584724417731408820873961880021
  c2 = 0.20609895474146851962537542011143008396458423224465328706920087552924788983302
  c1 = -0.180654273148821532237028827440264938703627661154871403363931834846318112573294
  c0 = 0.1606521926313221723948643034905571707676419356400609690156705532267686492918868
  x * vmuladd_fast(
    vmuladd_fast(
      vmuladd_fast(
        vmuladd_fast(
          vmuladd_fast(
            vmuladd_fast(vmuladd_fast(vmuladd_fast(c0, x, c1), x, c2), x, c3),
            x,
            c4,
          ),
          x,
          c5,
        ),
        x,
        c6,
      ),
      x,
      c7,
    ),
    x,
    c8,
  )
end
@inline function logkern_9(x)
  c9 = 1.442695040888963407587958876357516873223860332904434711906078224014586198339343
  c8 = -0.7213475204444817057711203479236040464424888160067823280588699959470059217833771
  c7 = 0.480898346962976128170017519528240603877475137860740135494448418612334521964808
  c6 = -0.3606737602222044528926485178925468152698544632914287402177742224438034341355266
  c5 = 0.2885390082734153359755295484533201676037233143490110170185338945645707433116462
  c4 = -0.2404491736638805517094434611266213341107664043116792359649609367170386448753448
  c3 = 0.2060990174486639706727491347594731948792605354335332681880496880277396108033051
  c2 = -0.1803364995021978052782243890647317433194251060644842283875753942662516598663569
  c1 = 0.1606200865414345528500280443621751155379142697854250695344975262087060612426921
  c0 = -0.1446222903636375678552838798467050326611974502493221959482776774515654211467748
  x * vmuladd_fast(
    vmuladd_fast(
      vmuladd_fast(
        vmuladd_fast(
          vmuladd_fast(
            vmuladd_fast(
              vmuladd_fast(vmuladd_fast(vmuladd_fast(c0, x, c1), x, c2), x, c3),
              x,
              c4,
            ),
            x,
            c5,
          ),
          x,
          c6,
        ),
        x,
        c7,
      ),
      x,
      c8,
    ),
    x,
    c9,
  )
end
@inline function logkern_7(x)
  c7 = 1.442695040888962266311320415441181965256745101864924017704730474290825179713386
  c6 = -0.7213475204444734849561916107798427186725025515630934214552341603017790847666743
  c5 = 0.4808983470003734483568608991659360631657097156445816019788134536213757250523677
  c4 = -0.3606737603148037624174654830835163329225507561910479104865837838596482628691627
  c3 = 0.2885388167992364548202154388340846682584745440281798969580530860342799719518947
  c2 = -0.2404488805785497044335968015014331102915764909299868970339917618519049206072763
  c1 = 0.2064127286335890836287774329056506792116396261746810041599847620260037837288407
  c0 = -0.1806895812056507095939337124225687685035687284931465242170306355644855623118632
  # c7 = 1.442695040888962266311320415440673226691509089399247319200433080789551133947874
  # c6 = -0.7213475204444734849561916106463536879222835723990934589483454213207143172699565
  # c5 = 0.4808983470003734483568609116686822486155905301580259810441289744019151890859293
  # c4 = -0.3606737603148037624174668433501054879718352818773442452949795967039102009618056
  # c3 = 0.2885388167992364548201727628199776711909314676794911979761485405529591339257216
  # c2 = -0.2404488805785497044302647238681858349839732924674030562755638231437855429968676
  # c1 = 0.2064127286335890836637376093489490181845837682333985654448331246682934490703086
  # c0 = -0.1806895812056507118633103913062383463992924441916921974716424429045368328209497
  x * vmuladd_fast(
    vmuladd_fast(
      vmuladd_fast(
        vmuladd_fast(
          vmuladd_fast(vmuladd_fast(vmuladd_fast(c0, x, c1), x, c2), x, c3),
          x,
          c4,
        ),
        x,
        c5,
      ),
      x,
      c6,
    ),
    x,
    c7,
  )
end
@inline function logkern_6(x)
  c6 = 1.44269504088896112481574913926580271258562668309654298701097420877374680031478
  c5 = -0.7213475204628791729494574154147857065200754915843724028524643433219030089552769
  c4 = 0.4808983470214135184345854322580473442742138878734777669129377480116915710051349
  c3 = -0.3606736095355376179878548883602333628571188158521650924319649055538192403728833
  c2 = 0.2885387593457860540619732799086814986482101737149964973666904711488634555687257
  c1 = -0.2407576763559816924463087687841796294903125418141104155734727426854999960321288
  c0 = 0.2064519502023243513328580164119739425968179159770070826936661162470222233739379
  x * vmuladd_fast(
    vmuladd_fast(
      vmuladd_fast(
        vmuladd_fast(vmuladd_fast(vmuladd_fast(c0, x, c1), x, c2), x, c3),
        x,
        c4,
      ),
      x,
      c5,
    ),
    x,
    c6,
  )
end

# Probably the best
@inline function vlog2_fast(a::AbstractSIMD{8,Float64})
  # log2(vgetmant(a, Val(8))) + vgetexp(a) == log2(a)
  m = vgetmant(a, Val(8)) # m ∈ [1,2)
  e = vgetexp(a)
  # m ∈ [1,2); We reduce further...
  y = inv_approx(m) # y ∈ (0.5,1]
  # return y
  y₀ = vroundscale(y, Val(80)) # 80 = 16*(1+4)
  # y₀ is a multiple of 32; 0.5:(1/32):1.0
  # log2(m) = y * x = log2(y) + log2(x)
  # r = x - 1
  # r + 1 = m*y;
  # log(r+1) = log(m) + log(y₀)
  # log(m) = log(1+r) - log(y₀)
  r = vfmsub(m, y₀, 1.0) # m * y - 1
  log1pr = logkern_5(r)
  # log1pr = logkern_6(r)
  # log1pr = logkern_7(r)
  # log1pr = logkern_9(r)
  # @show r log1pr m
  # log1pr = logkern_8(r)
  y₀notone = y₀ ≠ 1.0
  inds = (reinterpret(UInt, y₀) >>> 48) & 0x0f
  # @show y₀ inds r
  # logy₀ = vpermi2pd(inds, LOG2_TABLE_1, LOG2_TABLE_2)
  logy₀ = vpermi2pd(inds, LOG2_TABLE_1, LOG2_TABLE_2)
  # emlogy₀ = e - logy₀
  emlogy₀ = ifelse(y₀notone, e - logy₀, e)
  log1pr + emlogy₀
  # logm = ifelse(y₀notone, log1pr - logy₀, log1pr)
  # @show r y₀ e
  # return r, m, y₀, logy₀, e
  # logkern_5(r) - logy₀ + e
  # logm = ifelse(y₀isone, log1pr, log1pr - logy₀)
  # log1pr += e

  # logm = ifelse(y₀notone, log1pr - logy₀, log1pr)
  # @show logm log1pr logy₀ e
  # logm + e
end
@inline vlog_fast(x::T) where {T} = vlog2_fast(x) * convert(T, 0.6931471805599453)
@inline vlog10_fast(x::T) where {T} = vlog2_fast(x) * convert(T, 0.3010299956639812)
# @inline Base.FastMath.log_fast(v::AbstractSIMD) = vlog_fast(float(v))
# @inline Base.FastMath.log2_fast(v::AbstractSIMD) = vlog2_fast(float(v))
# @inline Base.FastMath.log10_fast(v::AbstractSIMD) = vlog10_fast(float(v))

@inline function log2_kern_5_256(x)
  c5 = 1.442695040888963430242018860033667730529246266516984575206349078995811595149434
  c4 = -0.7213475204444818238116304451613321925305725926194819917643836452178866375548279
  c3 = 0.480898346935995048804761640936874918715608085587068841980211938177173076051027
  c2 = -0.3606737601723789903264415189400867418476075783130968675245051405508304609911771
  c1 = 0.2885437254815682617764304037319501313717679206671517967342968352103864725323379
  c0 = -0.2404546770228689226321452773936263107841762894627528631584290286000941540854746
  x * vmuladd_fast(
    vmuladd_fast(vmuladd_fast(vmuladd_fast(vmuladd_fast(c0, x, c1), x, c2), x, c3), x, c4),
    x,
    c5,
  )
end

const LOG2_TABLE_128 =
  Float64[log2(x) for x ∈ range(big"0.5", step = 1 / 256, length = 128)]

@inline function vlog2_v2(a::AbstractSIMD{W,Float64}) where {W}
  # log2(vgetmant(a, Val(8))) + vgetexp(a) == log2(a)
  m = vgetmant(a, Val(8)) # m ∈ [1,2)
  e = vgetexp(a)
  # m ∈ [1,2); We reduce further...
  y = inv_approx(m) # y ∈ (0.5,1]
  # return y
  y₀ = vroundscale(y, Val(128)) # 80 = 16*(1+7)

  # y₀ is a multiple of 32; 0.5:
  # r + 1 = m*y;
  # log(r+1) = log(m) + log(y₀)
  # log(m) = log(1+r) - log(y₀)
  r = vfmsub(m, y₀, 1.0)
  # log1pr = logkern_5(r)
  log1pr = log2_kern_5_256(r)
  y₀notone = y₀ ≠ 1.0
  inds = (reinterpret(Int, y₀) >> 45) & 0xff
  # inds = (reinterpret(UInt, y₀) >>> 48) & 0x0f
  # @show y₀ inds r
  # logy₀ = vpermi2pd(inds, LOG2_TABLE_1, LOG2_TABLE_2)
  logy₀ = vload(zstridedpointer(LOG2_TABLE_128), (inds,))

  emlogy₀ = ifelse(y₀notone, e - logy₀, e)
  log1pr + emlogy₀
end

# @inline function Base.log(x1::AbstractSIMD{W,Float64}) where {W}
# @inline function vlog(x1::Float64)
# @inline Base.log(v::AbstractSIMD) = log(float(v))

# @inline function vlog_fast(x1::AbstractSIMD{W,Float32}) where {W}
#     notzero = x1 != zero(x1)
#     greater_than_zero = x1 > zero(x1)
#     x3 = x1 < 1.1754944f-38#3.4332275f-5
#     # x6 = true if x3 entirely false
#     x7 = x1 * 8.388608f6#14.0f0
#     x8 = ifelse(x3, x7, x1)
#     x10 = reinterpret(UInt32, x8)
#     x11 = x10 + 0x004afb0d
#     x12 = x11 >>> 0x00000017
#     x13 = ifelse(x3, 0xffffff6a, 0xffffff81)
#     x15 = x12 + x13
#     x16 = x11 & 0x007fffff
#     x17 = x16 + 0x3f3504f3
#     x18 = reinterpret(Float32, x17)
#     x19 = x18 - 1f0
#     x20 = x18 + 1f0
#     x21 = x19 / x20
#     x22 = x21 * x21
#     x23 = x22 * x22
#     x24 = vfmadd(x23, 0.24279079f0, 0.40000972f0)
#     x25 = x24 * x23
#     x26 = vfmadd(x23, 0.6666666f0, 0.6666666f0)
#     x27 = x26 * x22
#     x28 = x19 * x19
#     x29 = x28 * 5f-1
#     x30 = convert(Float32, x15 % Int32)
#     x31 = x27 + x29
#     x32 = x31 + x25
#     x33 = x30 * -0.00021219444f0
#     x34 = vfmadd(x21, x32, x33)
#     x35 = x19 - x29
#     x36 = x35 + x34
#     x37 = vfmadd(x30, 0.6933594f0, x36)
#     # x37 = vfmadd(x30, 0.6931472f0, x36)
#     x39 = ifelse(x1 == Inf32, Inf32, x37)
#     x40 = ifelse(notzero, x39, -Inf32)
#     ifelse(x1 < zero(x1), NaN32, x40)
# end
# @inline Base.FastMath.log_fast(v::AbstractSIMD) = vlog_fast(float(v))


text:
using Surrogates

#1D
n = 10
lb = 0.0
ub = 5.0
x = sample(n,lb,ub,SobolSample())
f = x-> x^2
y1 = f.(x)
der = x->2*x
y2 = der.(x)
y = vcat(y1,y2)

my_gek = GEK(x,y,lb,ub)
val = my_gek(2.0)
std_err = std_error_at_point(my_gek,1.0)
add_point!(my_gek,2.5,2.5^2)

#ND
n = 10
d = 2
lb = [0.0,0.0]
ub = [5.0,5.0]
x = sample(n,lb,ub,SobolSample())
f = x -> x[1]^2 + x[2]^2
y1 = f.(x)
grad1 = x -> 2*x[1]
grad2 = x -> 2*x[2]
function create_grads(n,d,grad1,grad2,y)
    c = 0
    y2 = zeros(eltype(y[1]),n*d)
    for i = 1:n
        y2[i+c] = grad1(x[i])
        y2[i+c+1] = grad2(x[i])
        c = c+1
    end
    return y2
end
y2 = create_grads(n,d,grad1,grad2,y)
y = vcat(y1,y2)

my_gek_ND = GEK(x,y,lb,ub)
val = my_gek_ND((1.0,1.0))
std_err = std_error_at_point(my_gek_ND,(1.0,1.0))
add_point!(my_gek_ND,(2.0,2.0),8.0)


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((false, 2, 3, 4, 0, 0, 0, 0, 0),),
    ((false, 2, 3, 4, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((true, 5, 3, 4, 0, 0, 0, 0, 0),),
    ((true, 5, 3, 4, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ]
insts["fast"] = [
    ((false, 5, 3, 4, 0, 0, 0, 0, 0),),
    ((false, 5, 3, 4, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((false, 5, 3, 4, 0, 0.1, 0.1, 0, 0),),
    ((false, 3, 4, 5, 0, 0, 0, 0, 0),),
    ((false, 3, 4, 5, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((false, 3, 4, 5, 0, 0.1, 0.1, 0, 0),),
    ((true, 5, 3, 4, 0, 0, 0, 0, 0),),
    ((true, 5, 3, 4, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((true, 5, 3, 4, 0, 0.1, 0.1, 0, 0),),
    ((true, 3, 4, 5, 0, 0, 0, 0, 0),),
    ((true, 3, 4, 5, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((true, 3, 4, 5, 0, 0.1, 0.1, 0, 0),),
    ((false, 15, 10, 20, 0, 0, 0, 0, 0),),
    ((false, 15, 10, 20, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((false, 15, 10, 20, 0, 0.1, 0.1, 0, 0),),
    ((true, 15, 10, 20, 0, 0, 0, 0, 0),),
    ((true, 15, 10, 20, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((true, 15, 10, 20, 0, 0.1, 0.1, 0, 0),),
    ((false, 100, 8, 12, 0, 0, 0, 0, 0),),
    ((false, 100, 8, 12, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((false, 100, 8, 12, 0, 0.1, 0.1, 0, 0),),
    ((true, 100, 8, 12, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((true, 100, 8, 12, 0, 0.1, 0.1, 0, 0),),
    ]
return (MatrixRegressionNative, insts)


text:
# functions to take a .net file and return a ParsedCircuit

function parse_2nodecomponent!(pc::ParsedCircuit{N}, line) where {N<:Number}
    regex1 = regex = r"^([RLCVI]\S+)\s+(\S+)\s+(\S+)\s+(.*)"i
    m = match(regex, line) # two nodes
    m === nothing && return false
    name = Symbol(m.captures[1])
    node_strings = (m.captures[2],m.captures[3])
    restofline = m.captures[4]
    regex2 = r"[a-z0-9_)}\'\"]()\s+[a-z_]"i
    m = match(regex2,restofline)
    if m === nothing
        value_string = restofline
        parameters_string = ""
    else
        value_string = restofline[1:m.offsets[1]]
        parameters_string = restofline[nextind(line,m.offsets[1]):end]
    end
    value = parse_spiceexpression(N, value_string)
    nodes = update_nodedict!(pc, node_strings)
    if line[1] == 'R'
        newcomponent = resistor(name, nodes, value, parameters_string)
    elseif line[1] == 'C'
        newcomponent = capacitor(name, nodes, value, parameters_string)
    elseif line[1] =='L'
        newcomponent = inductor(name, nodes, value, parameters_string)
    elseif line[1] =='V'
        newcomponent = voltageSource(name, nodes, value, parameters_string)
    elseif line[1] =='I'
        newcomponent = currentSource(name, nodes, value, parameters_string)
    end
    pc.max_element +=1
    push!(pc.netlist, newcomponent)
    return true
end

function update_nodedict!(pc::ParsedCircuit, node_strings)
    nodearray = Vector{Int}()
        for node_string in node_strings
            node = Symbol(node_string)
            if ~(node ∈ keys(pc.nodedict))
                pc.max_node +=1
                pc.nodedict[node]=pc.max_node
            end
        push!(nodearray, pc.nodedict[node])
        end
    Tuple(nodearray)
end

function parse_netlist(filename::AbstractString)
    io = open(filename)
    try
        pc = parse_netlist(io)
    finally
        close(io)
    end
    return pc
end
function parse_netlist(io::IO, N::Type=Float64)
    #seekstart(io)
    pc = ParsedCircuit{N}()
    pc.titleline = readline(io)
    for card in readcards(io)
        parse_2nodecomponent!(pc,card) && continue
        @warn "Not Processed: " card
    end
    return pc
end

function parse_spiceexpression(N,s)
    done, value = parse_spicevalue(N,s); done && return value
    done, expression = parse_spiceexpression(s); done && return expression
    throw(ErrorException("Could Not Process: $s"))
end
function parse_spicevalue(N,s)
    m = match(r"^\s*((?:[0-9]+(?:[.][0-9]*)?|[.][0-9]+)(?:[e][-+]?[0-9]+)?)
                (k|meg|mil|g|t|m|u|μ|n|p|f){0,1}[a-z0-9_.]*\s*$"ix,s)
    m === nothing && return (false, nothing)
    value = parse(N,m.captures[1])
    unitstring = m.captures[2]
    unitstring === nothing || (value *= N(units[lowercase(unitstring)]))
    return (true, value)
end
function parse_spiceexpression(s)
    s = fix_spiceunits(s)
    s = fix_netvoltages(s)
    s = fix_branchcurrents(s)
    (true,Base.Meta.parse(s))
end
function fix_spiceunits(s)
    regex = r"(?<![a-z_])(?<value>(?:[0-9]+(?:[.][0-9]*)?|[.][0-9]+)(?:[e][-+]?[0-9]+)?)
              (?<unit>k|meg|mil|g|t|m|u|μ|n|p|f)(?<junk>[a-z0-9_]*)()"ix
    m = match(regex,s)
    m === nothing && return s
    value,unit,junk,null = m.captures
    value_offset, unit_offset, junk_offset, next_offset = m.offsets
    s[1:prevind(s,value_offset)] * value * stringunits[lowercase(unit)] * fix_spiceunits(s[next_offset:end])
end
function fix_netvoltages(s)
    s = replace(s,r"v\(\s*(?<name>[^ ,]+)\s*\)"i => s"netvoltage(Val(Symbol(\"\g<name>\")))")
    regex = r"v\(\s*(?<name1>\S+)\s*,\s*(?<name2>\S+)\s*\)"i
    subex = s"netvoltage(Val(Symbol(\"\g<name1>\")))-netvoltage(Val(Symbol(\"\g<name2>\")))"
    replace(s,regex => subex)
end
function fix_branchcurrents(s)
    s = replace(s,r"i\(\s*(?<name>[^ ,]+)\s*\)"i => s"branchcurrent(Val(Symbol(\"\g<name>\")))")
    regex = r"i\(\s*(?<name1>\S+)\s*,\s*(?<name2>\S+)\s*\)"i
    subex = s"branchcurrent(Val(Symbol(\"\g<name1>\")))-branchcurrent(Val(Symbol(\"\g<name2>\")))"
    replace(s,regex => subex)
end


text:
"""
Currencies

This package provides the `Currency` singleton type, based on the ISO 4217 standard together with five methods:

- `symbol`: The symbol of the currency.
- `currency`: The singleton type for a particular currency symbol.
- `name`: The ISO 4217 name of the currency.
- `code`: The ISO 4217 code for the currency.
- `unit`: The ISO 4217 minor unit, i.e. number of decimal places, for the currency.

See README.md for the full documentation

Copyright 2019-2020, Eric Forgy, Scott P. Jones and other contributors

Licensed under MIT License, see LICENSE.md
"""
module Currencies

export Currency, symbol, currency, unit, code, name

"""
This is a singleton type, intended to be used as a label for dispatch purposes
"""
struct Currency{S}
    function Currency{S}() where {S}
        haskey(_currency_data, S) || error("Currency $S is not defined.")
        new{S}()
    end
end
Currency(S::Symbol) = Currency{S}()

include(joinpath(@__DIR__, "..", "deps", "currency-data.jl"))

"Returns the symbol associated with a currency"
function symbol end

symbol(::Type{Currency{S}}) where {S} = S

"Returns a singleton type Currency{S}"
function currency end

currency(S::Symbol) = _currency_data[S][1]

currency(::Type{C}) where {C<:Currency} = C

"Returns the ISO 4217 minor unit associated with a currency"
function unit end

"Returns the ISO 4217 code associated with a currency"
function code end

"Returns the ISO 4217 name associated with a currency"
function name end

ms = [:unit, :code, :name]
for (i,m) in enumerate(ms)
    @eval $m(S::Symbol) = _currency_data[S][$(i+1)]
    @eval $m(::Type{Currency{S}}) where {S} = $m(S)
end

"Returns all currency symbols"
allsymbols()  = keys(_currency_data)

"Returns all currency data as a pairs Symbol => (Currency,Unit,Code,Name)"
allpairs() = pairs(_currency_data)

end # module Currencies

max_stars_repo_name: TotalVerb/Currencies

max_issues_repo_name: TotalVerb/Currencies

max_forks_repo_name: TotalVerb/Currencies


text:
"""
    pred_vgp(x_pred, gp_params)
Obtain predictions from variational Gaussian process at locations x_pred.
Returns tuple with predicted means and predicted standard deviations
"""
function pred_vgp(x_pred, gp_params::SVGP_params)
    opt_xi    = gp_params.inducing_locs
    opt_m     = gp_params.inducing_mean
    opt_cov   = Hermitian(gp_params.inducing_L * transpose(gp_params.inducing_L))

    ni, nd    = size(opt_xi)

    rho        = exp.(gp_params.log_rho)
    kappa      = exp(gp_params.log_kappa[1])
    sigsq      = exp(gp_params.log_sigma[1])^2

    nugg_scale = 1.e-6 * kappa
    if nugg_scale < 1.e-8
        nugg_scale = 1.e-8
    end

    cov_pred   = covfn(x_pred, x_pred, gp_params) + I * nugg_scale
    cov_mat_i  = covfn(gp_params.inducing_locs, gp_params.inducing_locs, gp_params) + I * nugg_scale
    cross_mat  = covfn(gp_params.inducing_locs, x_pred, gp_params)
    cross_inv  = cov_mat_i \ cross_mat

    p_mean = cross_inv' * (opt_m .- gp_params.const_mean[1]) .+ gp_params.const_mean[1]
    p_cov  = diag( cov_pred + cross_inv' * ( opt_cov - cov_mat_i) * cross_inv);

    return p_mean, sqrt.(p_cov)
end

"""
    pred_vgp(x_pred, gp_obj)
Obtain predictions from variational Gaussian process at locations x_pred.
Returns tuple with predicted means and predicted standard deviations.
Included for backward compatibility, but only part of gp_obj used is the params
"""

function pred_vgp(x_pred, gp_obj::SVGP_obj)
    gp_params = gp_obj.params
    opt_xi    = gp_params.inducing_locs
    opt_m     = gp_params.inducing_mean
    opt_cov   = Hermitian(gp_params.inducing_L * transpose(gp_params.inducing_L))

    ni, nd    = size(opt_xi)

    rho        = exp.(gp_params.log_rho)
    kappa      = exp(gp_params.log_kappa[1])
    sigsq      = exp(gp_params.log_sigma[1])^2

    nugg_scale = 1.e-6 * kappa
    if nugg_scale < 1.e-8
        nugg_scale = 1.e-8
    end

    cov_pred   = covfn(x_pred, x_pred, gp_params) + I * nugg_scale
    cov_mat_i  = covfn(gp_params.inducing_locs, gp_params.inducing_locs, gp_params) + I * nugg_scale
    cross_mat  = covfn(gp_params.inducing_locs, x_pred, gp_params)
    cross_inv  = cov_mat_i \ cross_mat

    p_mean = cross_inv' * (opt_m .- gp_params.const_mean[1]) .+ gp_params.const_mean[1]
    p_cov  = diag( cov_pred + cross_inv' * ( opt_cov - cov_mat_i) * cross_inv);

    return p_mean, sqrt.(p_cov)
end

max_stars_repo_name: ayanbiswas/PRISM

max_issues_repo_name: ayanbiswas/PRISM

max_forks_repo_name: ayanbiswas/PRISM


text:
#===============================================================================

    McInput.jl

    defines struct of data and all analysis parameters
    to be passed to mc_analysis() in melting_curve.jl

    the constructor is intended as the only interface
    to the melting curve analysis and the only place
    where argument defaults are applied

    Author: Tom Price
    Date:   July 2019

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

import DataStructures.OrderedDict
import StaticArrays: SVector
import DataFrames.DataFrame

#===============================================================================
    defaults >>
===============================================================================#

## defaults for mc_analysis()
const DEFAULT_MC_AUTO_SPAN_SMOOTH               = false
const DEFAULT_MC_SPAN_SMOOTH_DEFAULT            = Float_T(0.015)
const DEFAULT_MC_SPAN_SMOOTH_FACTOR             = Float_T(7.2)
const DEFAULT_MC_MAX_TEMPERATURE                = Int_T(1000)

## defaults for mc_peak_analysis()
const DEFAULT_MC_TEMP_TOO_CLOSE_FRAC            = Float_T(0.05)
const DEFAULT_MC_TEMPERATURE_BANDWIDTH          = Float_T(1.0)
const DEFAULT_MC_DENSER_FACTOR                  = Int_T(10)    ## Integer
const DEFAULT_MC_SMOOTH_FLUO_SPLINE             = false
const DEFAULT_MC_PEAK_SPAN_TEMPERATURE          = Float_T(2.0)
const DEFAULT_MC_NEGDERIV_RANGE_LOW_QUANTILE    = Float_T(0.21)
const DEFAULT_MC_MAX_NUM_CROSS_POINTS           = Int_T(10)    ## Integer
const DEFAULT_MC_NOISE_FACTOR                   = Float_T(0.2)
const DEFAULT_MC_NORM_NEGDERIV_QUANTILE         = Float_T(0.64)
const DEFAULT_MC_MAX_NORM_NEGDERIV              = Float_T(0.8)
# const DEFAULT_MC_TOP1_FROM_MAX_UB               = 1.0
const DEFAULT_MC_MAX_NUM_PEAKS                  = Int_T(4)     ## Integer
const DEFAULT_MC_MIN_NORMALIZED_AREA            = Float_T(0.1)

## default for mutate_dups()
const DEFAULT_MC_JITTER_CONSTANT                = Float_T(0.01)

## defaults for mc_peak_analysis() overridden by mc_analysis()
# const OVERRIDDEN_DEFAULT_AUTO_SPAN_SMOOTH       = true
# const OVERRIDDEN_DEFAULT_MC_SPAN_SMOOTH_DEFAULT = 0.05
# const OVERRIDDEN_DEFAULT_MC_SPAN_SMOOTH_FACTOR  = 7.2

#===============================================================================
    Field definitions >>
===============================================================================#

## name, DataType, default value
const MC_FIELD_DEFS = [
    Field(:raw_df,                      DataFrame),
    Field(:num_wells,                   Int_T),
    Field(:num_channels,                Int_T),
    Field(:wells,                       SVector{W,Symbol} where {W}),
    Field(:channels,                    SVector{C,Int_T} where {C}),
    Field(:calibration_data,            CalibrationData{<: NumberOfChannels, <: Union{Int_T,Float_T}}),
    Field(:calibration_args,            CalibrationParameters, DEFAULT_CAL_ARGS),
    Field(:max_temperature,             Union{Int_T,Float_T},   DEFAULT_MC_MAX_TEMPERATURE),
    Field(:temp_too_close_frac,         Float_T,        	 	DEFAULT_MC_TEMP_TOO_CLOSE_FRAC),
    Field(:temperature_bandwidth,       Float_T,         		DEFAULT_MC_TEMPERATURE_BANDWIDTH),
    Field(:auto_span_smooth,            Bool,                	DEFAULT_MC_AUTO_SPAN_SMOOTH),
    Field(:span_smooth_default,         Float_T,         		DEFAULT_MC_SPAN_SMOOTH_DEFAULT),
    Field(:span_smooth_factor,          Float_T,         		DEFAULT_MC_SPAN_SMOOTH_FACTOR),
    Field(:denser_factor,               Int_T,               	DEFAULT_MC_DENSER_FACTOR),
    Field(:smooth_fluo_spline,          Bool,                	DEFAULT_MC_SMOOTH_FLUO_SPLINE),
    Field(:peak_span_temperature,       Float_T,         		DEFAULT_MC_PEAK_SPAN_TEMPERATURE),
    Field(:negderiv_range_low_quantile, Float_T,         		DEFAULT_MC_NEGDERIV_RANGE_LOW_QUANTILE),
    Field(:max_num_cross_points,        Int_T,               	DEFAULT_MC_MAX_NUM_CROSS_POINTS),
    Field(:noise_factor,                Float_T,         		DEFAULT_MC_NOISE_FACTOR),
    Field(:norm_negderiv_quantile,      Float_T,         		DEFAULT_MC_NORM_NEGDERIV_QUANTILE, "qt_prob"),
    Field(:max_norm_negderiv,           Float_T,         		DEFAULT_MC_MAX_NORM_NEGDERIV, "max_normd_qtv"),
    Field(:max_num_peaks,               Int_T,               	DEFAULT_MC_MAX_NUM_PEAKS, "top_N"),
    Field(:min_normalized_area,         Float_T,         		DEFAULT_MC_MIN_NORMALIZED_AREA),
    Field(:jitter_constant,             Float_T,         		DEFAULT_MC_JITTER_CONSTANT),
    Field(:out_format,                  OutputFormat,        	pre_json_output),
    Field(:reporting,                   Function,            	roundoff(JSON_DIGITS))]

#===============================================================================
    macro calls >>
===============================================================================#

## generate struct and constructor
SCHEMA = MC_FIELD_DEFS
@make_struct_from_SCHEMA McInput Input
@make_constructor_from_SCHEMA McInput

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
mutable struct VBGaussianMixtureZBer <: NaiveVariationalRule{GaussianMixture} end
outboundType(::Type{VBGaussianMixtureZBer}) = Message{Bernoulli}
function isApplicable(::Type{VBGaussianMixtureZBer}, input_types::Vector{<:Type})
    (length(input_types) == 6) || return false
    for (i, input_type) in enumerate(input_types)
        if (i == 2)
            (input_type == Nothing) || return false
        else
            matches(input_type, ProbabilityDistribution) || return false
        end
    end
    return true
end

mutable struct VBGaussianMixtureZCat <: NaiveVariationalRule{GaussianMixture} end
outboundType(::Type{VBGaussianMixtureZCat}) = Message{Categorical}
function isApplicable(::Type{VBGaussianMixtureZCat}, input_types::Vector{<:Type})
    (length(input_types) > 6) || return false
    iseven(length(input_types)) || return false
    for (i, input_type) in enumerate(input_types)
        if (i == 2)
            (input_type == Nothing) || return false
        else
            matches(input_type, ProbabilityDistribution) || return false
        end
    end
    return true
end

function matchPVInputs(input_types::Vector{<:Type})
    Nothing_positions = []
    p_positions = []
    for (i, input_type) in enumerate(input_types)
        if matches(input_type, ProbabilityDistribution)
            push!(p_positions, i)
        elseif (input_type == Nothing)
            push!(Nothing_positions, i)
        end
    end

    return (Nothing_positions, p_positions)
end

mutable struct VBGaussianMixtureM <: NaiveVariationalRule{GaussianMixture} end
outboundType(::Type{VBGaussianMixtureM}) = Message{GaussianMeanPrecision}
function isApplicable(::Type{VBGaussianMixtureM}, input_types::Vector{<:Type})
    n_inputs = length(input_types)
    iseven(n_inputs) || return false

    (Nothing_positions, p_positions) = matchPVInputs(input_types)
    n_Nothings = length(Nothing_positions)
    n_ps = length(p_positions)

    (n_Nothings == 1) || return false
    (n_Nothings + n_ps == n_inputs) || return false
    (1 in p_positions) || return false
    (2 in p_positions) || return false
    isodd(Nothing_positions[1]) || return false

    return true
end

mutable struct VBGaussianMixtureW <: NaiveVariationalRule{GaussianMixture} end
outboundType(::Type{VBGaussianMixtureW}) = Message{Union{Gamma, Wishart}}
function isApplicable(::Type{VBGaussianMixtureW}, input_types::Vector{<:Type})
    n_inputs = length(input_types)
    iseven(n_inputs) || return false

    (Nothing_positions, p_positions) = matchPVInputs(input_types)
    n_Nothings = length(Nothing_positions)
    n_ps = length(p_positions)

    (n_Nothings == 1) || return false
    (n_Nothings + n_ps == n_inputs) || return false
    (1 in p_positions) || return false
    (2 in p_positions) || return false
    iseven(Nothing_positions[1]) || return false

    return true
end

mutable struct VBGaussianMixtureOut <: NaiveVariationalRule{GaussianMixture} end
outboundType(::Type{VBGaussianMixtureOut}) = Message{GaussianWeightedMeanPrecision}
function isApplicable(::Type{VBGaussianMixtureOut}, input_types::Vector{<:Type})
    iseven(length(input_types)) || return false
    for (i, input_type) in enumerate(input_types)
        if (i == 1)
            (input_type == Nothing) || return false
        else
            matches(input_type, ProbabilityDistribution) || return false
        end
    end
    return true
end


text:
#= Requires enough llvm

while i <= L makes function >5x slower the while i < (L+1) #31416
ndinsmore opened this issue on Mar 20, 2019 · 27 comments
https://github.com/JuliaLang/julia/issues/31416
=#
#=vchuravy commented on Mar 21, 2019 •

While I do not recommend this, "there is always a way."
=#

    function unsafe_inc(i::Int64)
        Base.llvmcall("""
        %i = add nsw i64 %0, 1
        ret i64 %i
        """, Int64, Tuple{Int64}, i)
    end

function slow_entitlement_test_func(f,iter::Union{Base.KeySet{<:Any, <:Dict}, Base.ValueIterator{<:Dict}},N=1)
    ret=0
    @inbounds for n in 1:N
        ret=0
        d=iter.dict
        slots=d.slots
        L=length(slots)
        vals=d.vals
        i=d.idxfloor
        while i <= L
         #for i = d.idxfloor:L
            @inbounds if slots[i] === 0x1
                @inbounds v=vals[i]
                ret+=f(v)
            end
            i = unsafe_inc(i)
        end
    end
    return ret
end
#=
Since this is not a bug in Julia, but due to the defined semantics I am going to close this issue.
Please continue discussion on discourse.julialang.org
=#

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
module ITensorNetworkAD

using Reexport

include("ITensorChainRules/ITensorChainRules.jl")
include("ITensorAutoHOOT/ITensorAutoHOOT.jl")
include("ITensorNetworks/ITensorNetworks.jl")
include("Optimizations/Optimizations.jl")

end


text:
struct ArgIterator <: AbstractArgIter
    dict::Dict
    stable_arg::Vector{String}
    arg_list::Vector{String}
    done::Bool
    make_args
    ArgIterator(dict, stable_arg; arg_list=nothing, make_args=nothing) = new(dict, stable_arg, arg_list==nothing ? collect(keys(dict)) : arg_list, false, make_args)
end

function make_arguments(iter::ArgIterator, state)
    arg_list = Vector{String}()

    if iter.make_args == nothing
        new_ret_list = Vector{String}()
        for (arg_idx, arg) in enumerate(iter.arg_list)
            push!(new_ret_list, arg)
            push!(new_ret_list, string(iter.dict[arg][state[2][arg_idx]]))
        end
        arg_list = [new_ret_list; iter.stable_arg]
    else
        d = Dict{String, Union{String, Tuple}}()
        for (arg_idx, arg) in enumerate(iter.arg_list)
            if typeof(iter.dict[arg][state[2][arg_idx]]) <: Tuple
                d[arg] = string.(iter.dict[arg][state[2][arg_idx]])
            else
                d[arg] = string(iter.dict[arg][state[2][arg_idx]])
            end
        end
        arg_list = [iter.make_args(d); iter.stable_arg]
    end
end

function Base.iterate(iter::ArgIterator)
    state = (1, ones(Int64, length(iter.arg_list)))
    arg_list = make_arguments(iter, state)
    return (state[1], arg_list), next_state(iter, state)
end

function next_state(iter::ArgIterator, _state)
    state = _state[2]
    n_state = _state[1]

    state[end] += 1

    for (arg_idx, arg) in Iterators.reverse(enumerate(iter.arg_list))
        if arg_idx == 1
            return (n_state+1, state)
        end

        if state[arg_idx] > length(iter.dict[arg])
            state[arg_idx] = 1
            state[arg_idx - 1] += 1
        end

    end
end

function Base.iterate(iter::ArgIterator, state)
    if state[2][1] > length(iter.dict[iter.arg_list[1]])
        return nothing
    end
    arg_list = make_arguments(iter, state)
    return (state[1], arg_list), next_state(iter, state)
end

function Base.length(iter::ArgIterator)
    return *([length(iter.dict[key]) for key in iter.arg_list]...)
end