text:
using LinearAlgebra: checksquare
using LazySets.Arrays: isinvertible

# ==========================
# Exponentiation functions
# ==========================

"""
    AbstractExpAlg

Abstract supertype for all exponentiation algorithms.
"""
abstract type AbstractExpAlg end

# FIXME use _alias_exp to distinguish from setops
# no-op
_alias(alg::AbstractExpAlg) = alg

# alias defined on value types
_alias(alg::Symbol) = _alias(Val(alg))

"""
    BaseExpAlg <: AbstractExpAlg

Matrix exponential using the scaling and squaring algorithm implemented
in Julia Base.

### Notes

If the array is static, the method implemented in `StaticArrays.jl` is applied.
This algorithm admits the alias `:base`.
"""
struct BaseExpAlg <: AbstractExpAlg end

const BaseExp = BaseExpAlg()
_alias(alg::Val{:base}) = BaseExp

"""
    LazyExpAlg <: AbstractExpAlg

Lazy wrapper for the matrix exponential, operations defined in `LazySets.jl`.

### Fields

- `m`   -- (optional, default: `30`) size of the Krylov subspace
- `tol` -- (optional, default: `1e-10`) tolerance

### Notes

This algorithm admits the alias `:lazy` and also `:krylov`(using default options).
"""
struct LazyExpAlg <: AbstractExpAlg
    m::Int        # size of the Krylov subspace
    tol::Float64  # tolerance
end

# default constructor
LazyExpAlg() = LazyExpAlg(30, 1e-10)

const LazyExp = LazyExpAlg()
_alias(alg::Union{Val{:lazy}, Val{:krylov}}) = LazyExp

"""
    IntervalExpAlg <: AbstractExpAlg

Matrix exponential using an interval enclosure of the Taylor series remainder.

### Fields

- `order` -- order of the Taylor series expansion of the matrix exponential

### Notes

This method allows to overapproximate ``exp(Aδ)`` with an interval matrix.
It also accepts an interval matrix ``A`` given as input.
This method admits the alias `:interval` and `:taylor`.
"""
struct IntervalExpAlg <: AbstractExpAlg
    order::Int
end

# default constructor
IntervalExpAlg() = IntervalExpAlg(10)

const IntervalExp = IntervalExpAlg()
_alias(alg::Val{:interval}) = IntervalExp
_alias(alg::Val{:taylor}) = IntervalExp

"""
    PadeExpAlg <: AbstractExpAlg

Matrix exponential for sparse matrices using Pade approximants.

### Notes

Requires `Expokit.jl`. This algorithm admits the alias `:pade`.
"""
struct PadeExpAlg <: AbstractExpAlg end

const PadeExp = PadeExpAlg()
_alias(alg::Val{:pade}) = PadeExp

# general case: convert to Matrix
@inline _exp(A::AbstractMatrix, ::BaseExpAlg) = exp(Matrix(A))
@inline _exp(A::Matrix, ::BaseExpAlg) = exp(A)

# static arrays have their own exp method
@inline _exp(A::StaticArray, ::BaseExpAlg) = exp(A)

# exponential of an identity multiple, defined in MathematicalSystems.jl
@inline _exp(A::IdentityMultiple, ::BaseExpAlg) = exp(A)

# lazy wrapper of the matrix exponential
@inline _exp(A::AbstractMatrix, ::LazyExpAlg) = SparseMatrixExp(A)

# pade approximants (requires Expokit.jl)
@inline _exp(A::SparseMatrixCSC, ::PadeExpAlg) = padm(A)

function _exp(A::AbstractMatrix, alg::IntervalExpAlg)
    return exp_overapproximation(A, one(eltype(A)), alg.order)
end

function _exp(A::AbstractMatrix, δ, alg::IntervalExpAlg)
    return exp_overapproximation(A, δ, alg.order)
end

"""
    _exp(A::AbstractMatrix, δ, [alg]::AbstractExpAlg=BaseExp)

Compute the matrix exponential ``e^{Aδ}``.

### Input

- `A`    -- matrix
- `δ`    -- step size
- `alg`  -- (optional, default: `BaseExp`) the algorithm used to take the matrix
            exponential of `Aδ`, possible options are `BaseExp`, `LazyExp`, `PadeExp` and `IntervalExp`
            see details in the *Algorithm* section below

### Output

A matrix or a lazy wrapper of the matrix exponential, depending on `alg`.

### Algorithm

- `BaseExp` -- (alias: `:base`) use Higham's scaling and squaring method implemented
               in Julia standard library; see `?exp` for details; if `A` is a static array,
               uses the implementation in `StaticArrays.jl`

- `LazyExp` -- (alias: `:lazy`) return a lazy wrapper type around the matrix exponential
               using the implementation `LazySets.SparseMatrixExp`

- `PadeExp` -- (alias: `pade`) apply Pade approximant method to compute the matrix
               exponential of a sparse matrix (requires `Expokit.jl`)

- `IntervalExp` -- (alias: `interval`, `taylor`) apply the Taylor series expansion of the matrix
                    exponential with an interval remainder; works if `A` is an interval matrix

### Notes

If the algorithm `LazyExp` is used, evaluations of the action of the matrix
exponential are done with the `expmv` implementation from `Expokit`
(but see `LazySets#1312` for the planned generalization to other backends).
"""
function _exp(A::AbstractMatrix, δ, alg::AbstractExpAlg=BaseExp)
    n = checksquare(A)
    return _exp(A * δ, alg)
end

function _exp(A::AbstractMatrix, δ, alg::PadeExpAlg)
    n = checksquare(A)
    @requires Expokit
    return _exp(A * δ, alg)
end

@inline function _P_2n(A::AbstractMatrix{N}, δ, n) where {N}
    return [Matrix(A*δ) Matrix(δ*I, n, n);
            zeros(n, 2*n)]::Matrix{N}
end

@inline function _P_2n(A::SparseMatrixCSC{N, M}, δ, n) where {N, M}
    return [sparse(A*δ) sparse(δ*I, n, n);
            spzeros(n, 2*n)]::SparseMatrixCSC{N, M}
end

@inline function _P_3n(A::AbstractMatrix{N}, δ, n) where {N}
    return [Matrix(A*δ)     Matrix(δ*I, n, n)  zeros(n, n)    ;
            zeros(n, 2*n          )  Matrix(δ*I, n, n)        ;
            zeros(n, 3*n          )                           ]::Matrix{N}
end

@inline function _P_3n(A::SparseMatrixCSC{N, M}, δ, n) where {N, M}
    return [sparse(A*δ)     sparse(δ*I, n, n)  spzeros(n, n) ;
            spzeros(n, 2*n          )  sparse(δ*I, n, n)     ;
            spzeros(n, 3*n          )                        ]::SparseMatrixCSC{N, M}
end

"""
    Φ₁(A::AbstractMatrix, δ, [alg]::AbstractExpAlg=BaseExp, [isinv]::Bool=false, [Φ]=nothing)

Evaluate the series

```math
Φ₁(A, δ) = ∑_{i=0}^∞ \\dfrac{δ^{i+1}}{(i+1)!}A^i,
```
where ``A`` is a square matrix of order ``n`` and ``δ > 0``.

### Input

- `A`      -- coefficients matrix
- `δ`      -- step size
- `alg`    -- (optional, default: `BaseExp`) the method used to take the matrix
               exponential of the coefficient matrix; see the documentation of
               [`_exp`](@ref) for available options
- `isinv`  -- (optional, default: `false`) if `true`, assume that the coefficients
              matrix is invertible and compute ``A^{-1}``
- `Φ`      -- (optional, default: `nothing`) optionally pass the matrix exponential ``e^{Aδ}``

### Output

A matrix.

### Algorithm

If ``A`` is invertible, ``Φ₁`` can be computed as

```math
Φ₁(A, δ) = A^{-1}(e^{Aδ} - I_n),
```
where ``I_n`` is the identity matrix of order ``n``.

In the general case, implemented in this function, it can be
computed as submatrices of the block matrix

```math
P_{2n} = \\exp \\begin{pmatrix}
Aδ && δI_n \\\\
0 && 0
\\end{pmatrix}.
```
It can be shown that
```math
\\exp(P_{2n}) = \\begin{pmatrix}
Φ(A, δ) && Φ₁(A, δ) \\\\
0 &&  δI_n
\\end{pmatrix}.
```
where ``Φ(A, δ) = e^{Aδ}``. In particular, `Φ₁(A, δ) = P[1:n, (n+1):2*n]`.
This method can be found in [[FRE11]](@ref).
"""
function Φ₁(A::AbstractMatrix, δ, alg::AbstractExpAlg=BaseExp, isinv=false, Φ=nothing)
    return _Φ₁(A, δ, alg, Val(isinv), Φ)
end

# dispatch for discretization methods
_Φ₁(A, δ, alg, isinv::Val{:false}, Φ) = _Φ₁_blk(A, δ, alg)
_Φ₁(A, δ, alg, isinv::Val{:true}, Φ) = _Φ₁_inv(A, δ, alg, Φ)

# evaluate the series Φ₁(A, δ) = ∑_{i=0}^∞ \\dfrac{δ^{i+1}}{(i+1)!}A^i
# without assuming invertibility of A, by taking the exponential of a 2n x 2n matrix
function _Φ₁_blk(A, δ, alg)
    n = checksquare(A)
    P_2n = _P_2n(A, δ, n)
    Q = _exp(P_2n, alg)
    return _P₁_blk(Q, n, alg)
end

@inline _P₁_blk(P, n::Int, ::AbstractExpAlg) = P[1:n, (n+1):2*n]
@inline _P₁_blk(P, n::Int, ::LazyExpAlg) = sparse(get_columns(P, (n+1):2*n)[1:n, :])

# compute the matrix Φ₁(A, δ) = A^{-1}(exp(Aδ) - I), assuming that A is invertible
# and explicitly computing inv(A); this function optionally receives the matrix Φ = exp(Aδ)
function _Φ₁_inv(A::AbstractMatrix, δ, alg, Φ=nothing)
    if isnothing(Φ)
        Φ = _exp(A, δ, alg)
    end
    n = size(A, 1)
    N = eltype(A)
    In = Matrix(one(N)*I, n, n)
    Ainv = inv(A)
    return Ainv * (Φ - In)
end

# special case when A is a multiple of the identity
function _Φ₁_inv(A::IdentityMultiple, δ, alg, Φ=nothing)
    λ = A.M.λ
    @assert !iszero(λ) "the given identity multiple is not invertible"
    α = (1/λ) * (exp(δ*λ) - 1)
    return IdentityMultiple(α, size(A, 1))
end

# method to compute Φ₁ * u, where u is a vector
function _Φ₁_u(A, δ, alg, isinv::Val{:true}, u::AbstractVector, Φ=nothing)
    return _Φ₁_u(A, δ, alg, Φ, u)
end

function _Φ₁_u(A, δ, alg, isinv::Val{:false}, u::AbstractVector, Φ=nothing)
    M = _Φ₁_blk(A, δ, alg)
    return M * u
end

# compute the vector Φ₁(A, δ)u = A^{-1}(exp(Aδ) - I) u
# assuming that A is invertible and solving the associated linear system
# (the inverse of A is not explicitly computed)
function _Φ₁_u(A, δ, alg, u::AbstractVector, Φ=nothing)
    if isnothing(Φ)
        Φ = _exp(A, δ, alg)
    end
    x = Φ * u - u
    return A \ x
end

function load_Φ₁_krylov()
return quote

# compute Φ₁(A, δ)u = A^{-1}(exp(Aδ) - I) u without explicitly computing exp(Aδ)
# and assuming that A is invertible
function _Φ₁_u(A, δ, alg::LazyExpAlg, u::AbstractVector, ::Nothing)
    w = expv(1.0 * δ, A, u, m=alg.m, tol=alg.tol)
    x = w - u
    return A \ x
end

end end  # quote / load_Φ₁_krylov()

"""
    Φ₂(A::AbstractMatrix, δ, [alg]::AbstractExpAlg=BaseExp, [isinv]::Bool=false, [Φ]=nothing)

Evaluate the series

```math
Φ₂(A, δ) = ∑_{i=0}^∞ \\dfrac{δ^{i+2}}{(i+2)!}A^i,
```
where ``A`` is a square matrix of order ``n`` and ``δ > 0``.

### Input

- `A`      -- coefficients matrix
- `δ`      -- step size
- `alg`    -- (optional, default: `BaseExp`) the method used to take the matrix
               exponential of the coefficient matrix; see the documentation of
               [`_exp`](@ref) for available options
- `isinv`  -- (optional, default: `false`) if `true`, assume that the coefficients
              matrix is invertible and compute ``A^{-1}``
- `Φ`      -- (optional, default: `nothing`) optionally pass the matrix exponential ``e^{Aδ}``

### Output

A matrix.

### Algorithm

If ``A`` is invertible, ``Φ₂`` can be computed as

```math
Φ₂(A, δ) = A^{-2}(e^{δA} - I_n - δA).
```

In the general case, implemented in this function, it can be computed as
submatrices of the block matrix

```math
P_{3n} = \\exp \\begin{pmatrix}
Aδ && δI_n && 0 \\\\
0 && 0 && δI_n \\\\
0 && 0 && 0
\\end{pmatrix}.
```
It can be shown that
```math
\\exp(P_{3n}) = \\begin{pmatrix}
Φ(A, δ) && Φ₁(A, δ) && Φ₂(A, δ) \\\\
0 && I_n     && δI_n \\\\
0 && 0     && I_n
\\end{pmatrix}.
```
where ``Φ(A, δ) = e^{Aδ}``. In particular, `Φ₂ = P_{3n}[1:n, (2*n+1):3*n]`.
This method can be found in [[FRE11]](@ref).
"""
function Φ₂(A::AbstractMatrix, δ, alg::AbstractExpAlg=BaseExp, isinv=false, Φ=nothing)
    return _Φ₂(A, δ, alg, Val(isinv), Φ)
end

# dispatch for discretization methods
_Φ₂(A, δ, alg, isinv::Val{:false}, Φ) = _Φ₂_blk(A, δ, alg)
_Φ₂(A, δ, alg, isinv::Val{:true}, Φ) = _Φ₂_inv(A, δ, alg, Φ)

@inline _P₂_blk(P, n, ::AbstractExpAlg) = P[1:n, (2*n+1):3*n]
@inline _P₂_blk(P, n, ::LazyExpAlg) = sparse(get_columns(P, (2*n+1):3*n)[1:n, :])

# evaluate the series Φ₂(A, δ) = ∑_{i=0}^∞ \\dfrac{δ^{i+2}}{(i+2)!}A^i
# without assuming invertibility of A, by taking the exponential of a 3n x 3n matrix
function _Φ₂_blk(A, δ, alg)
    n = checksquare(A)
    B = _P_3n(A, δ, n)
    P = _exp(B, alg)
    return _P₂_blk(P, n, alg)
end

# compute the matrix Φ₂ = A^{-2} (exp(A*δ) - I - A*δ) assuming that A is invertible
# and explicitly computing inv(A); this function optionally receives Φ = exp(Aδ)
# FIXME don't pass algorithm since it is ignore if Φ is given
#@time Φ2 = ReachabilityAnalysis._Φ₂_inv(abs.(A), δ, ReachabilityAnalysis.BaseExp, Φ);
function _Φ₂_inv(A::AbstractMatrix, δ, alg, Φ=nothing)
    Aδ = A * δ
    if isnothing(Φ)
        Φ = _exp(Aδ, alg)
    end
    n = size(A, 1)
    N = eltype(A)
    In = Matrix(one(N)*I, n, n)
    B = Φ - In - Aδ
    Ainv = inv(Matrix(A))
    Ainvsqr = Ainv^2
    return Ainvsqr * B
end

# special case when A is a multiple of the identity
function _Φ₂_inv(A::IdentityMultiple, δ, alg, Φ=nothing)
    λ = A.M.λ
    @assert !iszero(λ) "the given identity multiple is not invertible"
    δλ = δ * λ
    α = (1/λ)^2 * (exp(δλ) - 1 - δλ)
    return IdentityMultiple(α, size(A, 1))
end

# ================
# Absolute values
# ================

@inline _elementwise_abs(A::AbstractMatrix) = abs.(A)
@inline _elementwise_abs(A::SparseMatrixCSC) = SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, abs.(nonzeros(A)))
@inline _elementwise_abs(A::IdentityMultiple) = IdentityMultiple(abs(A.M.λ), size(A, 1))

# ====================================
# Exponentiation of interval matrices
# ====================================

# TODO: outsource to IntervalMatrices.jl

# compute Iδ + 1/2 * δ^2 * A + 1/6 * δ^3 * A² + ... + 1/(η+1)! * δ^(η+1) * A^η + E(δ) * δ
function _Cδ(A, δ, order)
    n = size(A, 1)
    A² = A * A
    if isa(A, IntervalMatrix)
        Iδ = IntervalMatrix(Diagonal(fill(IA.Interval(δ), n)))
    else
        Iδ = Matrix(δ * I, n, n)
    end

    IδW = Iδ + 1/2 * δ^2 * A + 1/6 * δ^3 * A²
    M = IδW

    if order > 2
        # i = 2
        αᵢ₊₁ = 6 # factorial of (i+1)
        Aⁱ = A²
        δⁱ⁺¹ = δ^3
        @inbounds for i in 3:order
            αᵢ₊₁ *= i+1
            δⁱ⁺¹ *= δ
            Aⁱ *= A
            M += (δⁱ⁺¹/αᵢ₊₁) * Aⁱ
        end
    end
    E = IntervalMatrices._exp_remainder(A, δ, order; n=n)
    return M + E * δ
end


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

Utility functions for ImagePhaseCongruency

Copyright (c) 2015 Peter Kovesi
peterkovesi.com

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, subject to the following conditions:

The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.

The Software is provided "as is", without warranty of any kind.

PK August 2015
   October 2018 Julia 0.7/1.0

---------------------------------------------------------------------=#
import Images

export replacenan, fillnan
export hysthresh, imgnormalize

# Note the following two functions are in ImageProjective Geometry but are
# duplicated here to minimise dependencies.
export imgnormalise, histtruncate

#----------------------------------------------------------------------
# fillnan
"""
Fill NaN values in an image with closest non NaN value.

This can be used as a crude (but quick) 'inpainting' function to allow a FFT to
be computed on an image containing NaN values.  While the 'inpainting' is crude
it is typically good enough to remove most of the edge effects one might get at
the boundaries of the NaN regions.  The NaN regions should then be remasked out
of the final processed image.
```
Usage:  (newim, mask) = fillnan(img)

  Argument:  img   - Image to be 'filled'.
  Returns:   newim - Filled image.
             mask  - Binary image indicating the valid, non-NaN, regions in
                     the original image.
```
See also: [`replacenan`](@ref)
"""
function fillnan(img::AbstractArray{T,N}) where T where N

    mask = .!isnan.(img) # Non-NaN image regions

    if all(mask)
        newimg = copy(img)
        @warn("All elements are NaN, no filling possible");
        return newimg, mask
    end

    # Generate feature transform from non NaN regions of the image.
    # F will contain cartesian indices of closest non NaN points in the image
    F = Images.feature_transform(mask)

    # Fill NaN locations with value of closest non NaN pixel
    newimg = copy(img)

    for i in eachindex(F)
        newimg[i] = img[F[i]]
    end

    return newimg, mask
end
#----------------------------------------------------------------------
# replacenan
"""
Replace NaNs in an array with a specified value.
```
Usage: (newimg, mask) = replacenan(img, defaultval=0)

Arguments:
   img        - The Array containing NaN values.
   defaultval - The default value to replace NaNs.

Returns:
        newim - Filled image,
        mask  - Boolean image indicating non-NaN regions in the original
                image.
```
See also: [`fillnan`](@ref)
"""
function replacenan(img::AbstractArray{T,N}, defaultval::Real = 0) where T <: AbstractFloat where N
    mask = .!(isnan.(img))
    newimg = copy(img)
    newimg[isnan.(img)] .= defaultval
    return newimg, mask
end

#----------------------------------------------------------------------
# hysthresh
"""
Hysteresis thresholding of an image.
```
Usage: bw = hysthresh(img, T1, T2)

Arguments:
           img  - Image to be thresholded
        T1, T2  - Upper and lower threshold values.  T1 and T2
                  can be entered in any order, the larger of the
                  two values is used as the upper threshold.
Returns:
            bw  - The binary thresholded image as a BitArray
```
All pixels with values above threshold T1 are marked as edges. All pixels that
are connected to points that have been marked as edges and with values above
threshold T2 are also marked as edges. Eight connectivity is used.

"""
function hysthresh(img::AbstractArray{T0,2}, T1::Real, T2::Real) where T0 <: Real

    bw = falses(size(img))

    if T1 < T2    # Swap T1 and T2
        T1,T2 = T2,T1
    end

    # Form 8-connected components of pixels with a value above the
    # lower threshold and get the indices of pixels in each component.
    label = Images.label_components(img .>= T2, trues(3,3))
    pix = Images.component_indices(label)

    # For each list of pixels in pix test to see if there are any
    # image values above T1.  If so, set these pixels in the output
    # image.  Note we ignore pix[1] as these are the background pixels.
    for n = 2:length(pix)
        for i in eachindex(pix[n])
            if img[pix[n][i]] >= T1
                bw[pix[n]] .= true
                break
            end
        end
    end

    return bw
end

#----------------------------------------------------------------------
"""
imgnormalise/imgnormalize - Normalises image values to 0-1, or to desired mean and variance

```
Usage 1:      nimg = imgnormalise(img)
```
Offsets and rescales image so that the minimum value is 0
and the maximum value is 1.

```
Usage 2:      nimg = imgnormalise(img, reqmean, reqvar)

Arguments:  img     - A grey-level input image.
            reqmean - The required mean value of the image.
            reqvar  - The required variance of the image.
```
Offsets and rescales image so that nimg has mean reqmean and variance
reqvar.
"""
function imgnormalise(img::Array) # Normalise 0 - 1
    n = img .- minimum(img)
    return n = n/maximum(n)
end

# Normalise to desired mean and variance
function imgnormalise(img::Array, reqmean::Real, reqvar::Real)
    n = img .- mean(img)
    n /= std(img)      # Zero mean, unit std dev
    return n = reqmean .+ n*sqrt(reqvar)
end

# For those who spell normalise with a 'z'
"""
imgnormalize - Normalizes image values to 0-1, or to desired mean and variance
```
Usage 1:      nimg = imgnormalize(img)
```
Offsets and rescales image so that the minimum value is 0
and the maximum value is 1.
```
Usage 2:      nimg = imgnormalize(img, reqmean, reqvar)

Arguments:  img     - A grey-level input image.
            reqmean - The required mean value of the image.
            reqvar  - The required variance of the image.
```
Offsets and rescales image so that nimg has mean reqmean and variance
reqvar.
"""
function imgnormalize(img::Array)
    return imgnormalise(img)
end

function imgnormalize(img::Array, reqmean::Real, reqvar::Real)
    return imgnormalise(img, reqmean, reqvar)
end

#----------------------------------------------------------------------
"""
histtruncate - Truncates ends of an image histogram.

Function truncates a specified percentage of the lower and
upper ends of an image histogram.

This operation allows grey levels to be distributed across
the primary part of the histogram.  This solves the problem
when one has, say, a few very bright values in the image which
have the overall effect of darkening the rest of the image after
rescaling.

```
Usage:
1)   newimg = histtruncate(img, lHistCut, uHistCut)
2)   newimg = histtruncate(img, HistCut)

Arguments:
 Usage 1)
   img         -  Image to be processed.
   lHistCut    -  Percentage of the lower end of the histogram
                  to saturate.
   uHistCut    -  Percentage of the upper end of the histogram
                  to saturate.  If omitted or empty defaults to the value
                  for lHistCut.
 Usage 2)
   HistCut     -  Percentage of upper and lower ends of the histogram to cut.

Returns:
   newimg      -  Image with values clipped at the specified histogram
                  fraction values.  If the input image was colour the
                  lightness values are clipped and stretched to the range
                  0-1.  If the input image is greyscale no stretching is
                  applied. You may want to use imgnormalise() to achieve this.
```
See also: imgnormalise()
"""
function  histtruncate(img::Array, lHistCut::Real, uHistCut::Real)

    if lHistCut < 0 || lHistCut > 100 || uHistCut < 0 || uHistCut > 100
	error("Histogram truncation values must be between 0 and 100")
    end

    if ndims(img) > 2
	error("histtruncate only defined for grey scale images")
    end

    newimg = copy(img)
    sortv = sort(newimg[:])   # Generate a sorted array of pixel values.

    # Any NaN values will end up at the end of the sorted list. We
    # need to ignore these.
#    N = sum(.!isnan.(sortv))  # Number of non NaN values. v0.6
    N = sum(broadcast(!,isnan.(sortv)))  # compatibity for v0.5 and v0.6

    # Compute indicies corresponding to specified upper and lower fractions
    # of the histogram.
    lind = floor(Int, 1 + N*lHistCut/100)
    hind =  ceil(Int, N - N*uHistCut/100)

    low_val  = sortv[lind]
    high_val = sortv[hind]

    # Adjust image
    newimg[newimg .< low_val] .= low_val
    newimg[newimg .> high_val] .= high_val

    return newimg
end

function  histtruncate(img::Array, HistCut::Real)
    return histtruncate(img, HistCut, HistCut)
end

ext: jl


text:
const UNREGISTERED_PACKAGES = MirrorUpdater.Types.SrcDestPair[
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/bcbi/PredictMD.jl",
        destination_repo_name =
            "PredictMD.jl-3e7d7328-36f8-4388-bd01-4613c92c7370",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/bcbi/PredictMDExtra.jl",
        destination_repo_name =
            "PredictMDExtra.jl-d14d998a-9e6b-11e8-16d3-6f2879ea456d",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url = "https://github.com/bcbi/PredictMDFull.jl",
        destination_repo_name =
            "PredictMDFull.jl-5c0c5c38-9dd5-11e8-3ab7-453bd9ce6c97",
        ),
    MirrorUpdater.Types.SrcDestPair(
        ;
        source_url =
            "https://github.com/UnofficialJuliaMirror/MirrorUpdater.jl",
        destination_repo_name =
            "MirrorUpdater.jl-7f0fc4bf-4ffe-4e21-ab91-7bf1358d5ab3",
        ),
    ]


text:
module Compat

import LinearAlgebra
using LinearAlgebra: Adjoint, Diagonal, Transpose, UniformScaling, RealHermSymComplexHerm, BLAS

include("compatmacro.jl")

# https://github.com/JuliaLang/julia/pull/29440
if VERSION < v"1.1.0-DEV.389"
    Base.:(:)(I::CartesianIndex{N}, J::CartesianIndex{N}) where N =
        CartesianIndices(map((i,j) -> i:j, Tuple(I), Tuple(J)))
end

# https://github.com/JuliaLang/julia/pull/29442
if VERSION < v"1.1.0-DEV.403"
    Base.oneunit(::CartesianIndex{N}) where {N} = oneunit(CartesianIndex{N})
    Base.oneunit(::Type{CartesianIndex{N}}) where {N} = CartesianIndex(ntuple(x -> 1, Val(N)))
end

# https://github.com/JuliaLang/julia/pull/30268
if VERSION < v"1.1.0-DEV.811"
    Base.get(A::AbstractArray, I::CartesianIndex, default) = get(A, I.I, default)
end

# https://github.com/JuliaLang/julia/pull/29679
if VERSION < v"1.1.0-DEV.472"
    export isnothing
    isnothing(::Any) = false
    isnothing(::Nothing) = true
end

# https://github.com/JuliaLang/julia/pull/29749
if VERSION < v"1.1.0-DEV.792"
    export eachrow, eachcol, eachslice
    eachrow(A::AbstractVecOrMat) = (view(A, i, :) for i in axes(A, 1))
    eachcol(A::AbstractVecOrMat) = (view(A, :, i) for i in axes(A, 2))
    @inline function eachslice(A::AbstractArray; dims)
        length(dims) == 1 || throw(ArgumentError("only single dimensions are supported"))
        dim = first(dims)
        dim <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
        idx1, idx2 = ntuple(d->(:), dim-1), ntuple(d->(:), ndims(A)-dim)
        return (view(A, idx1..., i, idx2...) for i in axes(A, dim))
    end
end

function rangeargcheck(;step=nothing, length=nothing, kwargs...)
    if step===nothing && length===nothing
        throw(ArgumentError("At least one of `length` or `step` must be specified"))
    end
end

if VERSION < v"1.1.0-DEV.506"
    function Base.range(start, stop; kwargs...)
        rangeargcheck(;kwargs...)
        range(start; stop=stop, kwargs...)
    end
end

# https://github.com/JuliaLang/julia/pull/30496
if VERSION < v"1.2.0-DEV.272"
    Base.@pure hasfield(::Type{T}, name::Symbol) where T =
        Base.fieldindex(T, name, false) > 0
    export hasfield
    hasproperty(x, s::Symbol) = s in propertynames(x)
    export hasproperty
end

if VERSION < v"1.3.0-DEV.349"
    Base.findfirst(ch::AbstractChar, string::AbstractString) = findfirst(==(ch), string)
    Base.findnext(ch::AbstractChar, string::AbstractString, ind::Integer) =
        findnext(==(ch), string, ind)
    Base.findlast(ch::AbstractChar, string::AbstractString) = findlast(==(ch), string)
    Base.findprev(ch::AbstractChar, string::AbstractString, ind::Integer) =
        findprev(==(ch), string, ind)
end

# https://github.com/JuliaLang/julia/pull/29259
if VERSION < v"1.1.0-DEV.594"
    Base.merge(a::NamedTuple, b::NamedTuple, cs::NamedTuple...) = merge(merge(a, b), cs...)
    Base.merge(a::NamedTuple) = a
end

# https://github.com/JuliaLang/julia/pull/33129
if VERSION < v"1.4.0-DEV.142"
    export only

    Base.@propagate_inbounds function only(x)
        i = iterate(x)
        @boundscheck if i === nothing
            throw(ArgumentError("Collection is empty, must contain exactly 1 element"))
        end
        (ret, state) = i
        @boundscheck if iterate(x, state) !== nothing
            throw(ArgumentError("Collection has multiple elements, must contain exactly 1 element"))
        end
        return ret
    end

    # Collections of known size
    only(x::Ref) = x[]
    only(x::Number) = x
    only(x::Char) = x
    only(x::Tuple{Any}) = x[1]
    only(x::Tuple) = throw(
        ArgumentError("Tuple contains $(length(x)) elements, must contain exactly 1 element")
    )
    only(a::AbstractArray{<:Any, 0}) = @inbounds return a[]
    only(x::NamedTuple{<:Any, <:Tuple{Any}}) = first(x)
    only(x::NamedTuple) = throw(
        ArgumentError("NamedTuple contains $(length(x)) elements, must contain exactly 1 element")
    )
end

# https://github.com/JuliaLang/julia/pull/32628
if VERSION < v"1.3.0-alpha.8"
    Base.mod(i::Integer, r::Base.OneTo) = mod1(i, last(r))
    Base.mod(i::Integer, r::AbstractUnitRange{<:Integer}) = mod(i-first(r), length(r)) + first(r)
end

# https://github.com/JuliaLang/julia/pull/32739
# This omits special methods for more exotic matrix types, Triangular and worse.
if VERSION < v"1.4.0-DEV.92" # 2425ae760fb5151c5c7dd0554e87c5fc9e24de73

    # stdlib/LinearAlgebra/src/generic.jl
    LinearAlgebra.dot(x, A, y) = LinearAlgebra.dot(x, A*y) # generic fallback

    function LinearAlgebra.dot(x::AbstractVector, A::AbstractMatrix, y::AbstractVector)
        (axes(x)..., axes(y)...) == axes(A) || throw(DimensionMismatch())
        T = typeof(LinearAlgebra.dot(first(x), first(A), first(y)))
        s = zero(T)
        i₁ = first(eachindex(x))
        x₁ = first(x)
        @inbounds for j in eachindex(y)
            yj = y[j]
            if !iszero(yj)
                temp = zero(adjoint(A[i₁,j]) * x₁)
                @simd for i in eachindex(x)
                    temp += adjoint(A[i,j]) * x[i]
                end
                s += LinearAlgebra.dot(temp, yj)
            end
        end
        return s
    end
    LinearAlgebra.dot(x::AbstractVector, adjA::Adjoint, y::AbstractVector) =
         adjoint(LinearAlgebra.dot(y, adjA.parent, x))
    LinearAlgebra.dot(x::AbstractVector, transA::Transpose{<:Real}, y::AbstractVector) =
        adjoint(LinearAlgebra.dot(y, transA.parent, x))

    # stdlib/LinearAlgebra/src/diagonal.jl
    function LinearAlgebra.dot(x::AbstractVector, D::Diagonal, y::AbstractVector)
        mapreduce(t -> LinearAlgebra.dot(t[1], t[2], t[3]), +, zip(x, D.diag, y))
    end

    # stdlib/LinearAlgebra/src/symmetric.jl
    function LinearAlgebra.dot(x::AbstractVector, A::RealHermSymComplexHerm, y::AbstractVector)
        require_one_based_indexing(x, y)
        (length(x) == length(y) == size(A, 1)) || throw(DimensionMismatch())
        data = A.data
        r = zero(eltype(x)) * zero(eltype(A)) * zero(eltype(y))
        if A.uplo == 'U'
            @inbounds for j = 1:length(y)
                r += LinearAlgebra.dot(x[j], real(data[j,j]), y[j])
                @simd for i = 1:j-1
                    Aij = data[i,j]
                    r += LinearAlgebra.dot(x[i], Aij, y[j]) +
                        LinearAlgebra.dot(x[j], adjoint(Aij), y[i])
                end
            end
        else # A.uplo == 'L'
            @inbounds for j = 1:length(y)
                r += LinearAlgebra.dot(x[j], real(data[j,j]), y[j])
                @simd for i = j+1:length(y)
                    Aij = data[i,j]
                    r += LinearAlgebra.dot(x[i], Aij, y[j]) +
                        LinearAlgebra.dot(x[j], adjoint(Aij), y[i])
                end
            end
        end
        return r
    end

    # stdlib/LinearAlgebra/src/uniformscaling.jl
    LinearAlgebra.dot(x::AbstractVector, J::UniformScaling, y::AbstractVector) =
        LinearAlgebra.dot(x, J.λ, y)
    LinearAlgebra.dot(x::AbstractVector, a::Number, y::AbstractVector) =
        sum(t -> LinearAlgebra.dot(t[1], a, t[2]), zip(x, y))
    LinearAlgebra.dot(x::AbstractVector, a::Union{Real,Complex}, y::AbstractVector) =
        a*LinearAlgebra.dot(x, y)
end

# https://github.com/JuliaLang/julia/pull/30630
if VERSION < v"1.2.0-DEV.125" # 1da48c2e4028c1514ed45688be727efbef1db884
    require_one_based_indexing(A...) = !Base.has_offset_axes(A...) || throw(ArgumentError(
        "offset arrays are not supported but got an array with index other than 1"))
# At present this is only used in Compat inside the above dot(x,A,y) functions, #32739
elseif VERSION < v"1.4.0-DEV.92"
    using Base: require_one_based_indexing
end

# https://github.com/JuliaLang/julia/pull/33568
if VERSION < v"1.4.0-DEV.329"
    Base.:∘(f, g, h...) = ∘(f ∘ g, h...)
end
# https://github.com/JuliaLang/julia/pull/34251
if VERSION < v"1.5.0-DEV.56"
    Base.:∘(f) = f
end

# https://github.com/JuliaLang/julia/pull/33128
if VERSION < v"1.4.0-DEV.397"
    export pkgdir
    function pkgdir(m::Module)
        rootmodule = Base.moduleroot(m)
        path = pathof(rootmodule)
        path === nothing && return nothing
        return dirname(dirname(path))
    end
end

# https://github.com/JuliaLang/julia/pull/33736/
if VERSION < v"1.4.0-DEV.493"
    Base.Order.ReverseOrdering() = Base.Order.ReverseOrdering(Base.Order.Forward)
end

# https://github.com/JuliaLang/julia/pull/32968
if VERSION < v"1.4.0-DEV.551"
    Base.filter(f, xs::Tuple) = Base.afoldl((ys, x) -> f(x) ? (ys..., x) : ys, (), xs...)
    Base.filter(f, t::Base.Any16) = Tuple(filter(f, collect(t)))
end

# https://github.com/JuliaLang/julia/pull/34652
if VERSION < v"1.5.0-DEV.247"
    export ismutable
    ismutable(@nospecialize(x)) = (Base.@_pure_meta; typeof(x).mutable)
end

# https://github.com/JuliaLang/julia/pull/28761
export uuid5
if VERSION < v"1.1.0-DEV.326"
    import SHA
    import UUIDs: UUID
    function uuid5(ns::UUID, name::String)
        nsbytes = zeros(UInt8, 16)
        nsv = ns.value
        for idx in Base.OneTo(16)
            nsbytes[idx] = nsv >> 120
            nsv = nsv << 8
        end
        hash_result = SHA.sha1(append!(nsbytes, convert(Vector{UInt8}, codeunits(unescape_string(name)))))
        # set version number to 5
        hash_result[7] = (hash_result[7] & 0x0F) | (0x50)
        hash_result[9] = (hash_result[9] & 0x3F) | (0x80)
        v = zero(UInt128)
        #use only the first 16 bytes of the SHA1 hash
        for idx in Base.OneTo(16)
            v = (v << 0x08) | hash_result[idx]
        end
        return UUID(v)
    end
else
    using UUIDs: uuid5
end

# https://github.com/JuliaLang/julia/pull/34773
if VERSION < v"1.5.0-DEV.301"
    Base.zero(::AbstractIrrational) = false
    Base.zero(::Type{<:AbstractIrrational}) = false

    Base.one(::AbstractIrrational) = true
    Base.one(::Type{<:AbstractIrrational}) = true
end

# https://github.com/JuliaLang/julia/pull/32753
if VERSION < v"1.4.0-DEV.513"
    function evalpoly(x, p::Tuple)
        if @generated
            N = length(p.parameters)
            ex = :(p[end])
            for i in N-1:-1:1
                ex = :(muladd(x, $ex, p[$i]))
            end
            ex
        else
            _evalpoly(x, p)
        end
    end

    evalpoly(x, p::AbstractVector) = _evalpoly(x, p)

    function _evalpoly(x, p)
        N = length(p)
        ex = p[end]
        for i in N-1:-1:1
            ex = muladd(x, ex, p[i])
        end
        ex
    end

    function evalpoly(z::Complex, p::Tuple)
        if @generated
            N = length(p.parameters)
            a = :(p[end])
            b = :(p[end-1])
            as = []
            for i in N-2:-1:1
                ai = Symbol("a", i)
                push!(as, :($ai = $a))
                a = :(muladd(r, $ai, $b))
                b = :(muladd(-s, $ai, p[$i]))
            end
            ai = :a0
            push!(as, :($ai = $a))
            C = Expr(:block,
            :(x = real(z)),
            :(y = imag(z)),
            :(r = x + x),
            :(s = muladd(x, x, y*y)),
            as...,
            :(muladd($ai, z, $b)))
        else
            _evalpoly(z, p)
        end
    end
    evalpoly(z::Complex, p::Tuple{<:Any}) = p[1]

    evalpoly(z::Complex, p::AbstractVector) = _evalpoly(z, p)

    function _evalpoly(z::Complex, p)
        length(p) == 1 && return p[1]
        N = length(p)
        a = p[end]
        b = p[end-1]

        x = real(z)
        y = imag(z)
        r = 2x
        s = muladd(x, x, y*y)
        for i in N-2:-1:1
            ai = a
            a = muladd(r, ai, b)
            b = muladd(-s, ai, p[i])
        end
        ai = a
        muladd(ai, z, b)
    end
    export evalpoly
end

# https://github.com/JuliaLang/julia/pull/35304
if VERSION < v"1.5.0-DEV.574"
    Base.similar(A::PermutedDimsArray, T::Type, dims::Base.Dims) = similar(parent(A), T, dims)
end

# https://github.com/JuliaLang/julia/pull/34548
if VERSION < v"1.5.0-DEV.314"
    macro NamedTuple(ex)
        Meta.isexpr(ex, :braces) || Meta.isexpr(ex, :block) ||
            throw(ArgumentError("@NamedTuple expects {...} or begin...end"))
        decls = filter(e -> !(e isa LineNumberNode), ex.args)
        all(e -> e isa Symbol || Meta.isexpr(e, :(::)), decls) ||
            throw(ArgumentError("@NamedTuple must contain a sequence of name or name::type expressions"))
        vars = [QuoteNode(e isa Symbol ? e : e.args[1]) for e in decls]
        types = [esc(e isa Symbol ? :Any : e.args[2]) for e in decls]
        return :(NamedTuple{($(vars...),), Tuple{$(types...)}})
    end

    export @NamedTuple
end

# https://github.com/JuliaLang/julia/pull/34296
if VERSION < v"1.5.0-DEV.182"
    export mergewith, mergewith!
    _asfunction(f::Function) = f
    _asfunction(f) = (args...) -> f(args...)
    mergewith(f, dicts...) = merge(_asfunction(f), dicts...)
    mergewith!(f, dicts...) = merge!(_asfunction(f), dicts...)
    mergewith(f) = (dicts...) -> mergewith(f, dicts...)
    mergewith!(f) = (dicts...) -> mergewith!(f, dicts...)
end

# https://github.com/JuliaLang/julia/pull/32003
if VERSION < v"1.4.0-DEV.29"
    hasfastin(::Type) = false
    hasfastin(::Union{Type{<:AbstractSet},Type{<:AbstractDict},Type{<:AbstractRange}}) = true
    hasfastin(x) = hasfastin(typeof(x))
else
    const hasfastin = Base.hasfastin
end

# https://github.com/JuliaLang/julia/pull/34427
if VERSION < v"1.5.0-DEV.124"
    const FASTIN_SET_THRESHOLD = 70

    function isdisjoint(l, r)
        function _isdisjoint(l, r)
            hasfastin(r) && return !any(in(r), l)
            hasfastin(l) && return !any(in(l), r)
            Base.haslength(r) && length(r) < FASTIN_SET_THRESHOLD &&
                return !any(in(r), l)
            return !any(in(Set(r)), l)
        end
        if Base.haslength(l) && Base.haslength(r) && length(r) < length(l)
            return _isdisjoint(r, l)
        end
        _isdisjoint(l, r)
    end

    export isdisjoint
end

# https://github.com/JuliaLang/julia/pull/35577
if VERSION < v"1.5.0-DEV.681"
    Base.union(r::Base.OneTo, s::Base.OneTo) = Base.OneTo(max(r.stop,s.stop))
end

# https://github.com/JuliaLang/julia/pull/35929
# and also https://github.com/JuliaLang/julia/pull/29135 -> Julia 1.5
if VERSION < v"1.5.0-rc1.13" || v"1.6.0-" < VERSION < v"1.6.0-DEV.323"

    # Compat.stride not Base.stride, so as not to overwrite the method, and not to create ambiguities:
    function stride(A::AbstractArray, k::Integer)
        st = strides(A)
        k ≤ ndims(A) && return st[k]
        return sum(st .* size(A))
    end
    stride(A,k) = Base.stride(A,k) # Fall-through for other methods.

    # These were first defined for Adjoint{...,StridedVector} etc in #29135
    Base.strides(A::Adjoint{<:Real, <:AbstractVector}) = (stride(A.parent, 2), stride(A.parent, 1))
    Base.strides(A::Transpose{<:Any, <:AbstractVector}) = (stride(A.parent, 2), stride(A.parent, 1))
    Base.strides(A::Adjoint{<:Real, <:AbstractMatrix}) = reverse(strides(A.parent))
    Base.strides(A::Transpose{<:Any, <:AbstractMatrix}) = reverse(strides(A.parent))
    Base.unsafe_convert(::Type{Ptr{T}}, A::Adjoint{<:Real, <:AbstractVecOrMat}) where {T} = Base.unsafe_convert(Ptr{T}, A.parent)
    Base.unsafe_convert(::Type{Ptr{T}}, A::Transpose{<:Any, <:AbstractVecOrMat}) where {T} = Base.unsafe_convert(Ptr{T}, A.parent)

    Base.elsize(::Type{<:Adjoint{<:Real, P}}) where {P<:AbstractVecOrMat} = Base.elsize(P)
    Base.elsize(::Type{<:Transpose{<:Any, P}}) where {P<:AbstractVecOrMat} = Base.elsize(P)

end

# https://github.com/JuliaLang/julia/pull/27516
if VERSION < v"1.2.0-DEV.77"
    import Test: @inferred
    using Core.Compiler: typesubtract

    macro inferred(allow, ex)
        _inferred(ex, __module__, allow)
    end

    function _inferred(ex, mod, allow = :(Union{}))
        if Meta.isexpr(ex, :ref)
            ex = Expr(:call, :getindex, ex.args...)
        end
        Meta.isexpr(ex, :call)|| error("@inferred requires a call expression")
        farg = ex.args[1]
        if isa(farg, Symbol) && first(string(farg)) == '.'
            farg = Symbol(string(farg)[2:end])
            ex = Expr(:call, GlobalRef(Base.Test, :_materialize_broadcasted),
                farg, ex.args[2:end]...)
        end
        Base.remove_linenums!(quote
            let
                allow = $(esc(allow))
                allow isa Type || throw(ArgumentError("@inferred requires a type as second argument"))
                $(if any(a->(Meta.isexpr(a, :kw) || Meta.isexpr(a, :parameters)), ex.args)
                    # Has keywords
                    args = gensym()
                    kwargs = gensym()
                    quote
                        $(esc(args)), $(esc(kwargs)), result = $(esc(Expr(:call, _args_and_call, ex.args[2:end]..., ex.args[1])))
                        inftypes = $(gen_call_with_extracted_types(mod, Base.return_types, :($(ex.args[1])($(args)...; $(kwargs)...))))
                    end
                else
                    # No keywords
                    quote
                        args = ($([esc(ex.args[i]) for i = 2:length(ex.args)]...),)
                        result = $(esc(ex.args[1]))(args...)
                        inftypes = Base.return_types($(esc(ex.args[1])), Base.typesof(args...))
                    end
                end)
                @assert length(inftypes) == 1
                rettype = result isa Type ? Type{result} : typeof(result)
                rettype <: allow || rettype == typesubtract(inftypes[1], allow) || error("return type $rettype does not match inferred return type $(inftypes[1])")
                result
            end
        end)
    end
    #export @inferred
end

# https://github.com/JuliaLang/julia/pull/36360
if VERSION < v"1.6.0-DEV.322" # b8110f8d1ec6349bee77efb5022621fdf50bd4a5

    function guess_vendor()
        # like determine_vendor, but guesses blas in some cases
        # where determine_vendor returns :unknown
        ret = BLAS.vendor()
        if Base.Sys.isapple() && (ret == :unknown)
            ret = :osxblas
        end
        ret
    end

    """
        Compat.set_num_threads(n)

    Set the number of threads the BLAS library should use.

    Also accepts `nothing`, in which case julia tries to guess the default number of threads.
    Passing `nothing` is discouraged and mainly exists because,
    on exotic variants of BLAS, `nothing` may be returned by `get_num_threads()`.
    Thus the following pattern may fail to set the number of threads, but will not error:
    ```julia
    old = get_num_threads()
    set_num_threads(1)
    @threads for i in 1:10
        # single-threaded BLAS calls
    end
    set_num_threads(old)
    ```
    """
    set_num_threads(n)::Nothing = _set_num_threads(n)

    function _set_num_threads(n::Integer; _blas = guess_vendor())
        if _blas === :openblas || _blas == :openblas64
            return ccall((BLAS.@blasfunc(openblas_set_num_threads), BLAS.libblas), Cvoid, (Cint,), n)
        elseif _blas === :mkl
           # MKL may let us set the number of threads in several ways
            return ccall((:MKL_Set_Num_Threads, BLAS.libblas), Cvoid, (Cint,), n)
            elseif _blas === :osxblas
            # OSX BLAS looks at an environment variable
            ENV["VECLIB_MAXIMUM_THREADS"] = n
        else
            @assert _blas === :unknown
            @warn "Failed to set number of BLAS threads." maxlog=1
        end
        return nothing
    end
    _tryparse_env_int(key) = tryparse(Int, get(ENV, key, ""))

    function _set_num_threads(::Nothing; _blas = guess_vendor())
        n = something(
            _tryparse_env_int("OPENBLAS_NUM_THREADS"),
            _tryparse_env_int("OMP_NUM_THREADS"),
            max(1, Base.Sys.CPU_THREADS ÷ 2),
        )
        _set_num_threads(n; _blas = _blas)
    end

    """
        Compat.get_num_threads()

    Get the number of threads the BLAS library is using.

    On exotic variants of `BLAS` this function can fail,
    which is indicated by returning `nothing`.

    In Julia 1.6 this is `LinearAlgebra.BLAS.get_num_threads()`
    """
    get_num_threads(;_blas=guess_vendor())::Union{Int, Nothing} = _get_num_threads()

    function _get_num_threads(; _blas = guess_vendor())::Union{Int, Nothing}
        if _blas === :openblas || _blas === :openblas64
            return Int(ccall((BLAS.@blasfunc(openblas_get_num_threads), BLAS.libblas), Cint, ()))
        elseif _blas === :mkl
            return Int(ccall((:mkl_get_max_threads, BLAS.libblas), Cint, ()))
        elseif _blas === :osxblas
            key = "VECLIB_MAXIMUM_THREADS"
            nt = _tryparse_env_int(key)
            if nt === nothing
                @warn "Failed to read environment variable $key" maxlog=1
            else
                return nt
            end
        else
            @assert _blas === :unknown
        end
        @warn "Could not get number of BLAS threads. Returning `nothing` instead." maxlog=1
        return nothing
    end

else
    # Ensure that these can still be accessed as Compat.get_num_threads() etc:
    import LinearAlgebra.BLAS: set_num_threads, get_num_threads
end

# https://github.com/JuliaLang/julia/pull/30915
if VERSION < v"1.2.0-DEV.257" # e7e726b3df1991e1306ef0c566d363c0a83b2dea
    Base.:(!=)(x) = Base.Fix2(!=, x)
    Base.:(>=)(x) = Base.Fix2(>=, x)
    Base.:(<=)(x) = Base.Fix2(<=, x)
    Base.:(>)(x) = Base.Fix2(>, x)
    Base.:(<)(x) = Base.Fix2(<, x)
end

# https://github.com/JuliaLang/julia/pull/35132
if VERSION < v"1.5.0-DEV.639" # cc6e121386758dff6ba7911770e48dfd59520199
    export contains
    contains(haystack::AbstractString, needle) = occursin(needle, haystack)
    contains(needle) = Base.Fix2(contains, needle)
end

# https://github.com/JuliaLang/julia/pull/35052
if VERSION < v"1.5.0-DEV.438" # 0a43c0f1d21ce9c647c49111d93927369cd20f85
    Base.endswith(s) = Base.Fix2(endswith, s)
    Base.startswith(s) = Base.Fix2(startswith, s)
end

# https://github.com/JuliaLang/julia/pull/37517
if VERSION < v"1.6.0-DEV.1037"
    export ComposedFunction
    # https://github.com/JuliaLang/julia/pull/35980
    if VERSION < v"1.6.0-DEV.85"
        const ComposedFunction = let h = identity ∘ convert
            Base.typename(typeof(h)).wrapper
        end
        @eval ComposedFunction{F,G}(f, g) where {F,G} =
            $(Expr(:new, :(ComposedFunction{F,G}), :f, :g))
        ComposedFunction(f, g) = ComposedFunction{Core.Typeof(f),Core.Typeof(g)}(f, g)
    else
        using Base: ComposedFunction
    end
    function Base.getproperty(c::ComposedFunction, p::Symbol)
        if p === :f
            return getfield(c, :f)
        elseif p === :g
            return getfield(c, :g)
        elseif p === :outer
            return getfield(c, :f)
        elseif p === :inner
            return getfield(c, :g)
        end
        error("type ComposedFunction has no property ", p)
    end
    Base.propertynames(c::ComposedFunction) = (:f, :g, :outer, :inner)
else
    using Base: ComposedFunction
end

# https://github.com/JuliaLang/julia/pull/37244
if VERSION < v"1.6.0-DEV.873" # 18198b1bf85125de6cec266eac404d31ccc2e65c
    export addenv
    function addenv(cmd::Cmd, env::Dict)
        new_env = Dict{String,String}()
        if cmd.env !== nothing
            for (k, v) in split.(cmd.env, "=")
                new_env[string(k)::String] = string(v)::String
            end
        end
        for (k, v) in env
            new_env[string(k)::String] = string(v)::String
        end
        return setenv(cmd, new_env)
    end

    function addenv(cmd::Cmd, pairs::Pair{<:AbstractString}...)
        return addenv(cmd, Dict(k => v for (k, v) in pairs))
    end

    function addenv(cmd::Cmd, env::Vector{<:AbstractString})
        return addenv(cmd, Dict(k => v for (k, v) in split.(env, "=")))
    end
end

# https://github.com/JuliaLang/julia/pull/37559
if VERSION < v"1.6.0-DEV.1083"
    """
        reinterpret(reshape, T, A::AbstractArray{S}) -> B

    Change the type-interpretation of `A` while consuming or adding a "channel dimension."

    If `sizeof(T) = n*sizeof(S)` for `n>1`, `A`'s first dimension must be
    of size `n` and `B` lacks `A`'s first dimension. Conversely, if `sizeof(S) = n*sizeof(T)` for `n>1`,
    `B` gets a new first dimension of size `n`. The dimensionality is unchanged if `sizeof(T) == sizeof(S)`.

    # Examples

    ```
    julia> A = [1 2; 3 4]
    2×2 Matrix{$Int}:
     1  2
     3  4

    julia> reinterpret(reshape, Complex{Int}, A)    # the result is a vector
    2-element reinterpret(reshape, Complex{$Int}, ::Matrix{$Int}):
     1 + 3im
     2 + 4im

    julia> a = [(1,2,3), (4,5,6)]
    2-element Vector{Tuple{$Int, $Int, $Int}}:
     (1, 2, 3)
     (4, 5, 6)

    julia> reinterpret(reshape, Int, a)             # the result is a matrix
    3×2 reinterpret(reshape, $Int, ::Vector{Tuple{$Int, $Int, $Int}}):
     1  4
     2  5
     3  6
    ```
    """
    function Base.reinterpret(::typeof(reshape), ::Type{T}, a::A) where {T,S,A<:AbstractArray{S}}
        isbitstype(T) || throwbits(S, T, T)
        isbitstype(S) || throwbits(S, T, S)
        if sizeof(S) == sizeof(T)
            N = ndims(a)
        elseif sizeof(S) > sizeof(T)
            rem(sizeof(S), sizeof(T)) == 0 || throwintmult(S, T)
            N = ndims(a) + 1
        else
            rem(sizeof(T), sizeof(S)) == 0 || throwintmult(S, T)
            N = ndims(a) - 1
            N > -1 || throwsize0(S, T, "larger")
            axes(a, 1) == Base.OneTo(sizeof(T) ÷ sizeof(S)) || throwsize1(a, T)
        end
        paxs = axes(a)
        new_axes = if sizeof(S) > sizeof(T)
            (Base.OneTo(div(sizeof(S), sizeof(T))), paxs...)
        elseif sizeof(S) < sizeof(T)
            Base.tail(paxs)
        else
            paxs
        end
        reshape(reinterpret(T, vec(a)), new_axes)
    end

    @noinline function throwintmult(S::Type, T::Type)
        throw(ArgumentError("`reinterpret(reshape, T, a)` requires that one of `sizeof(T)` (got $(sizeof(T))) and `sizeof(eltype(a))` (got $(sizeof(S))) be an integer multiple of the other"))
    end
    @noinline function throwsize1(a::AbstractArray, T::Type)
        throw(ArgumentError("`reinterpret(reshape, $T, a)` where `eltype(a)` is $(eltype(a)) requires that `axes(a, 1)` (got $(axes(a, 1))) be equal to 1:$(sizeof(T) ÷ sizeof(eltype(a))) (from the ratio of element sizes)"))
    end
    @noinline function throwbits(S::Type, T::Type, U::Type)
        throw(ArgumentError("cannot reinterpret `$(S)` as `$(T)`, type `$(U)` is not a bits type"))
    end
    @noinline function throwsize0(S::Type, T::Type, msg)
        throw(ArgumentError("cannot reinterpret a zero-dimensional `$(S)` array to `$(T)` which is of a $msg size"))
    end
end

if VERSION < v"1.3.0-alpha.115"
    # https://github.com/JuliaLang/julia/pull/29634
    # Note this is much less performant than real 5-arg mul!, but is provided so old versions of julia don't error at least

    function _mul!(C, A, B, alpha, beta)
        Y = similar(C)
        LinearAlgebra.mul!(Y, A, B)
        C .= Y .* alpha .+ C .* beta
        return C
    end

    # all combination of Number and AbstractArray for A and B except both being Number
    function LinearAlgebra.mul!(C::AbstractArray, A::Number, B::AbstractArray, alpha::Number, beta::Number)
        return _mul!(C, A, B, alpha, beta)
    end
    function LinearAlgebra.mul!(C::AbstractArray, A::AbstractArray, B::Number, alpha::Number, beta::Number)
        return _mul!(C, A, B, alpha, beta)
    end
    function LinearAlgebra.mul!(C::AbstractArray, A::AbstractArray, B::AbstractArray, alpha::Number, beta::Number)
        return _mul!(C, A, B, alpha, beta)
    end
end

# https://github.com/JuliaLang/julia/pull/35243
if VERSION < v"1.6.0-DEV.15"
    _replace_filename(@nospecialize(x), filename, line_offset=0) = x
    function _replace_filename(x::LineNumberNode, filename, line_offset=0)
        return LineNumberNode(x.line + line_offset, filename)
    end
    function _replace_filename(ex::Expr, filename, line_offset=0)
        return Expr(
            ex.head,
            Any[_replace_filename(i, filename, line_offset) for i in ex.args]...,
        )
    end

    function parseatom(text::AbstractString, pos::Integer; filename="none")
        ex, i = Meta.parse(text, pos, greedy=false)
        return _replace_filename(ex, Symbol(filename)), i
    end

    function _skip_newlines(text, line, i)
        while i <= lastindex(text) && isspace(text[i])
            line += text[i] == '\n'
            i = nextind(text, i)
        end
        return line, i
    end

    function parseall(text::AbstractString; filename="none")
        filename = Symbol(filename)
        ex = Expr(:toplevel)
        line, prev_i = _skip_newlines(text, 1, firstindex(text))
        ex_n, i = Meta.parse(text, prev_i)
        while ex_n !== nothing
            push!(ex.args, LineNumberNode(line, filename))
            push!(ex.args, _replace_filename(ex_n, filename, line-1))
            line += count(==('\n'), SubString(text, prev_i:prevind(text, i)))
            line, prev_i = _skip_newlines(text, line, i)
            ex_n, i = Meta.parse(text, prev_i)
        end
        return ex
    end
else
    using .Meta: parseatom, parseall
end

include("iterators.jl")
include("deprecated.jl")

end # module Compat


text:
#=
# start of assignment 4 code
"""
    normalizeDNA(sequence)

Normalizes DNA sequences or bases to uppercase `String`s.
Throws an error if invalid bases are encountered.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> normalizeDNA("aaatg")
    "AAATG"

    julia> normalizeDNA("XTG")
    ERROR: Invalid base X encountered

    julia> normalizeDNA(5)
    ERROR: Invalid base 5 encountered

    julia> normalizeDNA('G')
    "G"

    julia> b = normalizeDNA('G');

    julia> typeof(b)
    String
"""
function normalizeDNA(sequence)
    # 1. make sequence into a String

    sequence = string(sequence)
    # 2. make sequence uppercase
    sequence = uppercase(sequence)

    for base in sequence
        if !occursin(base, "ACGT")
            error("Invalid base $base encountered")
        end
    end

    return sequence

    # return the normalized sequence
end

"""
    basecomposition(sequence)

Counts the number of each type of base
in a DNA sequence and returns a tuple those counts
in the order A, C, G, T

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> basecomposition("AATCGGG")
    (2, 1, 3, 1)

    julia> basecomposition('C')
    (0, 1, 0, 0)

    julia> A,C,G,T = basecomposition("accgggtttt")
    (1, 2, 3, 4)

    julia> A
    1

    julia> T
    4

    julia> basecomposition("BACCGGGTTTT")
    ERROR: Invalid base B encountered
"""
function basecomposition(sequence)
    sequence = normalizeDNA(sequence) # make uppercase string, check invalid bases
    a = c = g = t = 0 # sets all 4 variables to `0`

    for base in sequence
        # add 1 to each base as it occurs
        #println(base)
        if base == 'A'
            a = a+1
        elseif base == 'C'
            c = c+1
        elseif base == 'G'
            g = g+1
        elseif base == 'T'
            t = t+1
        end

    end

    return a,c,g,t
end

"""
   kmercount(sequence, k)

Finds all kmers in a sequence,
returning a dictionary of those kmers
and the number of times they appear in the sequence.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> kmercount("ggg", 3)
    Dict{Any,Any} with 1 entry:
    "GGG" => 1

    julia> kmercount("ATATATATA", 4)
    Dict{Any,Any} with 2 entries:
    "TATA" => 3
    "ATAT" => 3

    julia> kmercount("ATATATATAx", 4)
    ERROR: Invalid base X encountered

    julia> kmercount("A", 2)
    ERROR: k must be a positive integer less than the length of the sequence
"""
function  kmercount(sequence, k)
    1 <= k <= length(sequence) || error("k must be a positive integer less than the length of the sequence")
    kmers = Dict() # initialize dictionary

    # We're going to loop through the string with numerical index,
    # each time grabbing the bases at position i through i+k-1.
    # What is the last index that we should search?
    stopindex = length(sequence)-k+1

    for i in 1:stopindex
        kmer = sequence[i:i+k-1] # Change to index the sequence from i to i+k-1
        #println(kmer)
        kmer = normalizeDNA(kmer)
        #   if this kmer is a key the dictionary
        if haskey(kmers, kmer)
        #       add 1 to the value referenced by that kmer
            kmers[kmer] = kmers[kmer] + 1
        #   otherwise
        else
        #       make a new entry in the dictionary with a value of 1
            kmers[kmer] = 1
            #setindex!(kmers, 1, kmer)
        end
        #println(kmers)

    end
    return kmers
end
# end of assignment 4 code
=#


text:
const SString = Union{AbstractString, Symbol}
const Maybe{T} = Union{Missing, Nothing, T}

function middle(range::AbstractRange)
    array = collect(range)
    N = length(array)
    return array[ceil(Int, N / 2)]
end

function convert(::Type{GdkRGBA}, c::T) where T <: ColorAlpha
    color = convert(RGBA, c)
    return GdkRGBA(color.r, color.g, color.b, color.alpha)
end

convert(::Type{GdkRGBA}, color::T) where T <: Colorant   = convert(GdkRGBA, convert(RGBA, color))
convert(::Type{T}, color::GdkRGBA) where T <: ColorAlpha = convert(T, RGBA(color.r, color.g, color.b, color.a))
convert(::Type{T}, color::GdkRGBA) where T <: Color      = convert(T, RGB(color.r, color.g, color.b))

set!(w::GtkColorButton, k::SString, color::Colorant) = set!(w, k, convert(Gtk.GdkRGBA, color))

function set!(w::GtkWidget, k::SString, props::Vector{Pair})
    for (key, value) in props
        set!(w, key, value)
    end
end
"""
        @margin values...

Useful macro for setting the `margins` of a widget.

| Example             | Behavior                                                      |
|:--------------------|:--------------------------------------------------------------|
|`@margin(s)`         | Sets all margins to `s`                                       |
|`@margin(v, h)`      | Sets `top-bottom` and `left-right` margins respectively       |
|`@margin(t, b, h)`   | Sets `top`, `bottom` and `left-right` margins respectively    |
|`@margin(t, r, b, l)`| Sets `top`, `right`, `bottom` and `left` margins respectively |

# Examples

```julia
julia> Button(:Foo, @margin(10)) # Sets all margins to 10
GtkButton...

julia> Button(:Foo, @margin(10, 20)) # Sets top-bottom margins to 10 and left-right margins to 20
GtkButton...
```
"""
macro margin(t, r, b, l)
    return esc(Expr(:kw, :margin, :((
        top = $t,
        left = $l,
        right = $r,
        bottom = $b,
    ))))
end

macro margin(t, b, lr)
    esc(:( @margin($t, $lr, $b, $lr) ))
end

macro margin(tb, lr)
    esc(:( @margin($tb, $lr, $tb, $lr) ))
end

macro margin(s)
    esc(:( @margin($s, $s, $s, $s) ))
end

function align(value::Symbol)
    galign = Symbol(uppercase(string(value)))
    @assert isdefined(GtkAlign, galign) "Alignment type $value unknow."
    return getfield(GtkAlign, galign)
end

macro valign(ex)
    Expr(:kw, :valign, :( align($(esc(ex))) ))
end

macro halign(ex)
    Expr(:kw, :halign, :( align($(esc(ex))) ))
end

"""
        @align v h

Sets the vertical and horizontal alignment of Gtk widget.

# Examples

```julia
julia> Button("Foo", @align(:center)) # Sets both alignments to center
GtkButton...

julia> Button("Foo", @align(:center, :end))
GtkButton...
```
"""
macro align(v, h)
    esc(Expr(:kw, :align, :(Pair[
        :valign => align($v),
        :halign => align($h)
    ])))
end

macro align(value)
    esc(:(@align($value, $value)))
end

"""
        @spacing r c

Sets the `row` and `column` spacing of a `GtkGrid` widget.

# Examples

```julia
julia> Grid(@spacing(10)) # sets both row and column spacing.
GtkGrid...

julia> Grid(@spacing(10, 20))
GtkGrid...
```
"""
macro spacing(r, c)
    esc(Expr(:kw, :____, :(Pair[
        :row_spacing => $r,
        :column_spacing => $c
    ])))
end

macro spacing(value)
    esc(:(@spacing($value, $value)))
end

macro expand()
    Expr(:kw, :expand, true)
end

macro expand(ex)
    esc(Expr(:kw, :expand, ex))
end

macro hexpand()
    Expr(:kw, :hexpand, true)
end

macro hexpand(ex)
    esc(Expr(:kw, :hexpand, ex))
end

macro vexpand()
    Expr(:kw, :vexpand, true)
end

macro vexpand(ex)
    esc(Expr(:kw, :vexpand, ex))
end

macro homogeneous(r, c)
    esc(Expr(:kw, :homogeneous, :(Pair[
        :column_homogeneous => $(r),
        :row_homogeneous => $(c),
    ])))
end

macro homogeneous(v)
    esc(:(@homogeneous($v, $v)))
end

macro homogeneous()
   :(@homogeneous(true, true))
end

macro on(ex...)
    if length(ex) == 1
        throw(error("The @on macro must have at least 2 arguments (event name and callback function)"))
    end
    eventname = ("on" * join(map(string, ex[1:end-1]), "_")) |> Symbol
    esc(Expr(:kw, eventname, last(ex)))
end

"""
        keyboard(key)

Returns a number that represents a `key` of the keyboard.
"""
function keyboard(key::Union{AbstractString, Symbol})
    try
        gdk_key = Symbol("GDK_KEY_$key")
        getfield(Gtk.GConstants, gdk_key)
    catch e
        throw(error("Keyboard '$key' key not found."))
    end
end


text:
function make_idmap(table; col_user = :userid, col_item = :itemid)
    table, user2uidx = reindex_id_column(table, col_user)
    table, item2iidx = reindex_id_column(table, col_item)
    iidx2item = Dict(iidx => itemid for (itemid, iidx) in item2iidx)
    return table, user2uidx, item2iidx, iidx2item
end

struct StopTrain <: Exception end

abstract type AbstractCallback end

function initialize!(
    cb::AbstractCallback,
    model::AbstractRecommender;
    col_user = :userid,
    col_item = :itemid,
)

    throw("initialize! not implemented.")
end

struct LogTrainLoss <: AbstractCallback end

initialize!(
    cb::LogTrainLoss,
    model::AbstractRecommender;
    col_user = :userid,
    col_item = :itemid,
) = nothing

function (cb::LogTrainLoss)(model::AbstractRecommender, train_loss, epoch, verbose = -1)
    if verbose >= 1 && (epoch % verbose == 0)
        @info "epoch=$epoch: train_loss=$train_loss"
    end
end

"""
    EvaluateValidData(valid_metric::MeanMetric, valid_table, early_stopping_rounds, name = "val_metric")

Callback to monitor the validation metrics during training, and raise `StopTrain` exception if early stopping is requred.

# Constructor arguments
- `valid_metric`: monotring metric. See [Evaluation metrics](@ref) for the available ones.
- `valid_table`: any `Tables.jl`-compatible object for validation dataset.
- `early_stopping_rounds`: If the validation metric does not improve more than this epochs, the early stopping is invoked. If set to be less than 1, no early stopping is applied.
- `name`: metrics name to show on logger.

# Example
Use in the matrix factorizaion training.
```julia
ndcg10 = MeanNDCG(10)
cb = EvaluateValidData(ndcg10, test_table, 1, "val_NDCG")

model = ImplicitMF(16, true, 0.01)
fit!(
    model,
    train_table,
    10,
    callbacks = [cb],
    col_item = :movieid,
    n_epochs = 20,
    n_negatives = 1,
    learning_rate = 0.01,
    verbose = 1,
)
```
"""
mutable struct EvaluateValidData <: AbstractCallback
    valid_metric::MeanMetric
    name::Union{AbstractString,Symbol}
    valid_table::Any
    early_stopping_rounds::Int
    best_epoch::Int
    best_val_metric::Float64
    topk::Int
    val_xs::Any
    val_ys::Any

    EvaluateValidData(
        valid_metric::MeanMetric,
        valid_table,
        early_stopping_rounds,
        name = "val_metric",
    ) = new(
        valid_metric,
        name,
        valid_table,
        early_stopping_rounds,
        1,
        0.0,
        -1,
        nothing,
        nothing,
    )
end

function initialize!(
    cb::EvaluateValidData,
    model::AbstractRecommender;
    col_user = :userid,
    col_item = :itemid,
)
    # TODO: loss as metric
    cb.topk = cb.valid_metric.base_metric.k
    cb.val_xs, cb.val_ys =
        make_u2i_dataset(cb.valid_table, col_user = col_user, col_item = col_item)
    preds = predict_u2i(model, cb.val_xs, cb.topk, drop_history = true)
    cb.best_val_metric = cb.valid_metric(preds, cb.val_ys)
end

function (cb::EvaluateValidData)(
    model::AbstractRecommender,
    train_loss,
    epoch,
    verbose = -1,
)
    preds = predict_u2i(model, cb.val_xs, cb.topk, drop_history = true)
    current_metric = cb.valid_metric(preds, cb.val_ys)

    if current_metric > cb.best_val_metric
        cb.best_epoch = epoch
        cb.best_val_metric = current_metric
    end

    if verbose >= 1 && (epoch % verbose == 0)
        @info "epoch=$epoch: $(cb.name)=$current_metric, best_$(cb.name)=$(cb.best_val_metric), best_epoch=$(cb.best_epoch)"
    end

    if cb.early_stopping_rounds >= 1 &&
       ((epoch - cb.best_epoch) >= cb.early_stopping_rounds)
        throw(StopTrain)
    end
end


text:
module NaD

end

max_stars_repo_name: Nazri-Ahmad/NaD

max_issues_repo_name: Nazri-Ahmad/NaD

max_forks_repo_name: Nazri-Ahmad/NaD


text:
# test the structures
@info "Testing FT and NaN for the structures...";
@testset "StomtaModels --- structures" begin
    for FT in [Float32, Float64]
        can_3 = CanopyLayer{FT}(n_leaf=2);
        esm_1 = ESMBallBerry{FT}();
        esm_2 = ESMGentine{FT}();
        esm_3 = ESMLeuning{FT}();
        esm_4 = ESMMedlyn{FT}();
        osm_1 = OSMEller{FT}();
        osm_2 = OSMSperry{FT}();
        osm_3 = OSMWang{FT}();
        osm_4 = OSMWAP{FT}();
        osm_5 = OSMWAPMod{FT}();

        # test the structures
        for result in [ can_3, esm_1, esm_2, esm_3, esm_4, osm_4, osm_5 ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end
    end
end

# test the structures
println();
@info "Testing gas exchange functions...";
@testset "StomtaModels --- gas exchange functions" begin
    for FT in [Float32, Float64]
        envir  = AirLayer{FT}();
        leaf_3 = Leaf{FT}();
        leaf_4 = Leaf{FT}();
        mod_3  = C3CLM(FT);
        mod_4  = C4CLM(FT);
        rand_T = rand(FT) + 298;
        can_3  = CanopyLayer{FT}(n_leaf=2);
        can_4  = CanopyLayer{FT}(n_leaf=2);
        esm_1  = ESMBallBerry{FT}();
        esm_2  = ESMGentine{FT}();
        esm_3  = ESMLeuning{FT}();
        esm_4  = ESMMedlyn{FT}();
        osm_1  = OSMEller{FT}();
        osm_2  = OSMSperry{FT}();
        osm_3  = OSMWang{FT}();
        osm_4  = OSMWAP{FT}();
        osm_5  = OSMWAPMod{FT}();
        hs     = LeafHydraulics{FT}();

        # test the refresh functions
        update_leaf_TP!(mod_3, can_3, hs, envir);
        update_leaf_TP!(mod_4, can_4, hs, envir);
        @test NaN_test(can_3);
        @test NaN_test(can_4);

        update_leaf_AK!(mod_3, can_3, hs, envir);
        update_leaf_AK!(mod_4, can_4, hs, envir);
        @test NaN_test(can_3);
        @test NaN_test(can_4);

        gas_exchange!(mod_3, can_3, envir, GlcDrive(), 1, FT(0.1));
        gas_exchange!(mod_4, can_4, envir, GlcDrive(), 1, FT(0.1));
        @test NaN_test(can_3);
        @test NaN_test(can_4);

        gas_exchange!(mod_3, can_3, envir, GswDrive(), 1, FT(0.05));
        gas_exchange!(mod_4, can_4, envir, GswDrive(), 1, FT(0.05));
        @test NaN_test(can_3);
        @test NaN_test(can_4);

        can_3.g_sw[2] = 0;
        can_4.g_sw[2] = 0;
        gsw_control!(mod_3, can_3, envir, 2);
        gsw_control!(mod_4, can_4, envir, 2);
        @test NaN_test(can_3);
        @test NaN_test(can_4);
    end
end

# test the structures
println();
@info "Testing empirical formulations...";
@testset "StomtaModels --- empirical formulations" begin
    for FT in [Float32, Float64]
        envir  = AirLayer{FT}();
        leaf_3 = Leaf{FT}();
        leaf_4 = Leaf{FT}();
        mod_3  = C3CLM(FT);
        mod_4  = C4CLM(FT);
        rand_T = rand(FT) + 298;
        can_3  = CanopyLayer{FT}(n_leaf=2);
        can_4  = CanopyLayer{FT}(n_leaf=2);
        esm_1  = ESMBallBerry{FT}();
        esm_2  = ESMGentine{FT}();
        esm_3  = ESMLeuning{FT}();
        esm_4  = ESMMedlyn{FT}();
        osm_1  = OSMEller{FT}();
        osm_2  = OSMSperry{FT}();
        osm_3  = OSMWang{FT}();
        osm_4  = OSMWAP{FT}();
        osm_5  = OSMWAPMod{FT}();
        hs     = LeafHydraulics{FT}();

        # test the empirical model formulations
        for result in [ stomatal_conductance(esm_1, leaf_3, envir, FT(1)),
                        stomatal_conductance(esm_2, leaf_3, envir, FT(1)),
                        stomatal_conductance(esm_3, leaf_3, envir, FT(1)),
                        stomatal_conductance(esm_4, leaf_3, envir, FT(1)),
                        stomatal_conductance(esm_1, can_3, envir, FT(1)),
                        stomatal_conductance(esm_2, can_3, envir, FT(1)),
                        stomatal_conductance(esm_3, can_3, envir, FT(1)),
                        stomatal_conductance(esm_4, can_3, envir, FT(1)),
                        stomatal_conductance(esm_1, can_3, envir, FT(1), 1),
                        stomatal_conductance(esm_2, can_3, envir, FT(1), 1),
                        stomatal_conductance(esm_3, can_3, envir, FT(1), 1),
                        stomatal_conductance(esm_4, can_3, envir, FT(1), 1) ]
            @test FT_test(result, FT);
            @test NaN_test(result);
        end
    end
end

# test the structures
println();
@info "Testing stomatal models...";
@testset "StomtaModels --- stomatal models" begin
    for FT in [Float32, Float64]
        envir  = AirLayer{FT}();
        leaf_3 = Leaf{FT}();
        leaf_4 = Leaf{FT}();
        mod_3  = C3CLM(FT);
        mod_4  = C4CLM(FT);
        rand_T = rand(FT) + 298;
        can_3  = CanopyLayer{FT}(n_leaf=2);
        can_4  = CanopyLayer{FT}(n_leaf=2);
        esm_1  = ESMBallBerry{FT}();
        esm_2  = ESMGentine{FT}();
        esm_3  = ESMLeuning{FT}();
        esm_4  = ESMMedlyn{FT}();
        osm_1  = OSMEller{FT}();
        osm_2  = OSMSperry{FT}();
        osm_3  = OSMWang{FT}();
        osm_4  = OSMWAP{FT}();
        osm_5  = OSMWAPMod{FT}();
        hs     = LeafHydraulics{FT}();

        # test the solution functions
        for (mod,can) in zip([mod_3, mod_4], [can_3, can_4])
            for sm in [esm_1, esm_2, esm_3, esm_4]
                for result in [ solution_diff!(FT(0.1), mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaGLinearPleaf{FT}(), GlcDrive(), 1),
                                solution_diff!(FT(0.1), mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaGLinearPsoil{FT}(), GlcDrive(), 1),
                                solution_diff!(FT(0.1), mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaGLinearSWC{FT}(), GlcDrive(), 1),
                                solution_diff!(FT(0.1), mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaVLinearPleaf{FT}(), GlcDrive(), 1),
                                solution_diff!(FT(0.1), mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaVLinearPsoil{FT}(), GlcDrive(), 1),
                                solution_diff!(FT(0.1), mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaVLinearSWC{FT}(), GlcDrive(), 1),
                                ]
                    @test FT_test(result, FT);
                    @test NaN_test(result);
                end
            end
            for sm in [osm_1, osm_2, osm_3, osm_4, osm_5]
                result = solution_diff!(FT(0.1), mod, can, hs, envir, sm, GlcDrive(), 1);
                @test FT_test(result, FT);
                @test NaN_test(result);
            end
        end

        # test the stomata solutions
        for (mod,can) in zip([mod_3, mod_4], [can_3, can_4])
            for sm in [esm_1, esm_2, esm_3, esm_4]
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaGLinearKleaf{FT}(), 1);
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaGLinearPleaf{FT}(), 1);
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaGLinearPsoil{FT}(), 1);
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaGLinearSWC{FT}(), 1);
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaVLinearKleaf{FT}(), 1);
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaVLinearPleaf{FT}(), 1);
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaVLinearPsoil{FT}(), 1);
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaVLinearSWC{FT}(), 1);
                @test NaN_test(can);
            end
            for sm in [osm_1, osm_2, osm_3, osm_4, osm_5]
                gas_exchange!(mod, can, hs, envir, sm, 1);
                @test NaN_test(can);
            end
            gas_exchange!(mod, can, TreeSimple{FT}(), envir, osm_3);
            @test NaN_test(can);
        end

        for (mod,can) in zip([mod_3, mod_4], [can_3, can_4])
            for sm in [esm_1, esm_2, esm_3, esm_4]
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaGLinearKleaf{FT}());
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaGLinearPleaf{FT}());
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaGLinearPsoil{FT}());
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaGLinearSWC{FT}());
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaVLinearKleaf{FT}());
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaVLinearPleaf{FT}());
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaVLinearPsoil{FT}());
                @test NaN_test(can);
                gas_exchange!(mod, can, hs, FT(-1), FT(0.4), envir, sm, BetaVLinearSWC{FT}());
                @test NaN_test(can);
            end
            for sm in [osm_1, osm_2, osm_3, osm_4, osm_5]
                gas_exchange!(mod, can, hs, envir, sm);
                @test NaN_test(can);
            end
        end

        # test the nighttime stomatal conductance solution for Wang model
        can_3.APAR[1] = 0;
        gas_exchange!(mod_3, can_3, hs, envir, osm_3);

        # test the prognostic g_sw functions
        prognostic_gsw!(can_3, envir, esm_1, FT(1), FT(120));
        prognostic_gsw!(mod_3, can_3, hs, envir, osm_3, FT(120));
        @test true;
    end
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using OMEinsumContractionOrders
using Test

@testset "kahypar" begin
    include("kahypar.jl")
end


text:
using Benchmark
################################################################################
# Benchmark to see a better way of implementing crop
#
# | Row | Function     | Average   | Relative | Replications |
# |-----|--------------|-----------|----------|--------------|
# | 1   | "crop_loop"  | 0.0260486 | 1.0      | 50           |
# | 2   | "crop_whole" | 0.109872  | 4.21796  | 50           |
################################################################################

function crop_loop(input, crop_wh, output)
  w_off = div(size(input,1) - crop_wh[1], 2)
  h_off = div(size(input,2) - crop_wh[2], 2)

  for n = 1:size(input,4)
    for c = 1:size(input,3)
      for h = 1:crop_wh[2]
        @simd for w = 1:crop_wh[1]
          @inbounds output[w,h,c,n] = input[w+w_off,h+h_off,c,n]
        end
      end
    end
  end
end

function crop_whole(input, crop_wh, output)
  w_off = div(size(input,1) - crop_wh[1], 2)
  h_off = div(size(input,2) - crop_wh[2], 2)

  output[:] = input[w_off+1:w_off+crop_wh[1], h_off+1:h_off+crop_wh[2],:,:]
end

input = rand(256, 256, 3, 128)
crop_wh = [227, 227]
out1 = zeros(crop_wh..., size(input,3), size(input,4))
out2 = similar(out1)
crop_loop(input, crop_wh, out1)
crop_whole(input, crop_wh, out2)
@assert all(abs.(out1-out2) .< 1e-10)

crop_loop() = crop_loop(input, crop_wh, out1)
crop_whole() = crop_whole(input, crop_wh, out2)

df = compare([crop_loop, crop_whole], 50)
println("$df")


text:
module HessenbergBench

using IterativeSolvers
using BenchmarkTools

function backslash_versus_givens(; n = 1_000, ms = 10 : 10 : 100)
    results = BenchmarkGroup()
    results["backslash"] = BenchmarkGroup()
    results["givens_qr"] = BenchmarkGroup()

    # Example matrix
    A = Tridiagonal(fill(-1.5, n - 1), fill(3.0, n), fill(-1.0, n - 1))

    # Different sizes of the Hessenberg matrix
    for m = ms
        println(m)

        H = zeros(m + 1, m)

        # Create an orthonormal basis for the Krylov subspace
        V = rand(n, m + 1)
        V[:, 1] /= norm(V[:, 1])

        for i = 1 : m
            # Expand
            V[:, i + 1] = A * V[:, i]

            # Orthogonalize
            H[1 : i, i] = view(V, :, 1 : i)' * V[:, i + 1]
            V[:, i + 1] -= view(V, :, 1 : i) * H[1 : i, i]

            # Re-orthogonalize
            update = view(V, :, 1 : i)' * V[:, i + 1]
            H[1 : i, i] += update
            V[:, i + 1] -= view(V, :, 1 : i) * update

            # Normalize
            H[i + 1, i] = norm(V[:, i + 1])
            V[:, i + 1] /= H[i + 1, i]
        end

        # Fist standard basis vector as rhs
        rhs = [i == 1 ? 1.0 : 0.0 for i = 1 : size(H, 1)]

        # Run the benchmark
        results["givens_qr"][m] = @benchmark A_ldiv_B!(IterativeSolvers.Hessenberg(myH), myRhs) setup = (myH = copy($H); myRhs = copy($rhs))
        results["backslash"][m] = @benchmark $H \ $rhs
    end

    return results
end

end


text:
using EchelleCCFs
using Test

@testset "Check CCF accuracy" begin  # TODO repeat for multiple mask shapes
    @testset "Tophat mask w/ 1 line" begin
        Δv = 1000
        λ1 = 5025.0
        d1 = 0.5
        σ_v = 5000.0
        λlo = 5000
        λhi = 5050
        npix = 2000
        @test_nowarn BasicLineList([λ1], [d1])
        ll = BasicLineList([λ1], [d1])
        ll = BasicLineList([λ1], [d1])
        @test_nowarn BasicCCFPlan(line_list=ll, mask_shape=TopHatCCFMask(Δv))
        ccfpl = BasicCCFPlan(line_list=ll, mask_shape=TopHatCCFMask(Δv))
        @test_nowarn calc_ccf_v_grid(ccfpl)
        v_grid = calc_ccf_v_grid(ccfpl)
        @test length(v_grid) == calc_length_ccf_v_grid(ccfpl)
        λs = exp.(range(log(λlo),stop=log(λhi),length=npix))
        flux = ones(size(λs))
        flux .*= 1 .- d1.*exp.(-0.5.*((λs.-λ1)./λ1.*(EchelleCCFs.speed_of_light_mps./σ_v)).^2)
        @test_nowarn ccf_1D(λs, flux, ccfpl)
        ccf = ccf_1D(λs, flux, ccfpl)
        @test sum(ccf) > 0
        @test_nowarn  MeasureRvFromCCFGaussian()
        mrfcg = MeasureRvFromCCFGaussian(frac_of_width_to_fit=0.75)
        @test_nowarn measure_rv_from_ccf(v_grid,ccf, alg=mrfcg )
        vfit = measure_rv_from_ccf(v_grid,ccf, alg=mrfcg )
        @test vfit.rv ≈ 0  atol = 0.1   # TODO: tune tolerance better
        var = flux./10^5
        @test_nowarn ccf_1D(λs, flux, var, ccfpl)
        (ccf, ccf_var) = ccf_1D(λs, flux, var, ccfpl)
        @test sum(ccf) > 0
        @test_nowarn measure_rv_from_ccf(v_grid,ccf, ccf_var, alg=mrfcg )
        vfit = measure_rv_from_ccf(v_grid,ccf,ccf_var, alg=mrfcg )  # TODO: Check Why is this doing so much worse?
        #=
        using Plots
        global plt = plot()
        plot!(plt,v_grid,ccf,label=:none)
        scatter!(plt,v_grid,ccf,color=1,label=:none)
        plot!(plt,v_grid,sqrt.(ccf_var),label=:none)
        =#
        @test vfit.rv ≈ 0  atol = 0.5   # TODO: tune tolerance better
    end

    @testset "Tophat mask w/ 2 lines" begin
        Δv = 1000
        λ1 = 5025.0
        λ2 = 5075.0
        d1 = 0.5
        d2 = 0.25
        σ_v = 5000.0
        λlo = 5000
        λhi = 5100
        npix = 8000
        @test_nowarn BasicLineList([λ1,λ2], [d1,d2])
        ll = BasicLineList([λ1,λ2], [d1,d2])
        ll = BasicLineList([λ1,λ2], [d1,d2])
        @test_nowarn BasicCCFPlan(line_list=ll, mask_shape=TopHatCCFMask(Δv))
        ccfpl = BasicCCFPlan(line_list=ll, mask_shape=TopHatCCFMask(Δv))
        @test_nowarn calc_ccf_v_grid(ccfpl)
        v_grid = calc_ccf_v_grid(ccfpl)
        @test length(v_grid) == calc_length_ccf_v_grid(ccfpl)
        λs = exp.(range(log(λlo),stop=log(λhi),length=npix))
        flux = ones(size(λs))
        flux .*= 1 .- d1.*exp.(-0.5.*((λs.-λ1)./λ1.*(EchelleCCFs.speed_of_light_mps./σ_v)).^2)
        flux .*= 1 .- d2.*exp.(-0.5.*((λs.-λ2)./λ2.*(EchelleCCFs.speed_of_light_mps./σ_v)).^2)
        @test_nowarn ccf_1D(λs, flux, ccfpl)
        var = flux./10^5
        @test_nowarn ccf_1D(λs, flux, var, ccfpl)
        (ccf, ccf_var) = ccf_1D(λs, flux, var, ccfpl)

        @test_nowarn  MeasureRvFromCCFGaussian()
        mrfcg = MeasureRvFromCCFGaussian()
        @test_nowarn measure_rv_from_ccf(v_grid,ccf, alg=mrfcg )
        @test_nowarn measure_rv_from_ccf(v_grid,ccf, ccf_var, alg=mrfcg )
        vfit = measure_rv_from_ccf(v_grid,ccf, alg=mrfcg )
        @test vfit.rv ≈ 0  atol = 0.3   # TODO: tune tolerance better
        vfit = measure_rv_from_ccf(v_grid,ccf,ccf_var, alg=mrfcg )  # TODO: Check Why is this doing so much worse?
        #=
        using Plots
        global plt = plot()
        plot!(plt,v_grid,ccf,label=:none)
        scatter!(plt,v_grid,ccf,color=1,label=:none)
        plot!(plt,v_grid,sqrt.(ccf_var),label=:none)
        =#
        @test vfit.rv ≈ 0  atol = 0.3   # TODO: tune tolerance better

    end
    @testset "Gaussian mask" begin
        Δv = 7e3
        λ1 = 5025.0
        λ2 = 5075.0
        d1 = 0.5
        d2 = 0.25
        σ_v = 5000.0
        λlo = 5000
        λhi = 5100
        npix = 4000
        @test_nowarn BasicLineList([λ1,λ2], [d1,d2])
        ll = BasicLineList([λ1,λ2], [d1,d2])
        @test_nowarn BasicCCFPlan(line_list=ll, mask_shape=GaussianCCFMask(σ_v))
        ccfpl = BasicCCFPlan(line_list=ll, mask_shape=GaussianCCFMask(σ_v))
        @test_nowarn calc_ccf_v_grid(ccfpl)
        v_grid = calc_ccf_v_grid(ccfpl)
        @test length(v_grid) == calc_length_ccf_v_grid(ccfpl)
        λs = exp.(range(log(λlo),stop=log(λhi),length=npix))
        flux = ones(size(λs))
        flux .*= 1 .- d1.*exp.(-0.5.*((λs.-λ1)./λ1.*(EchelleCCFs.speed_of_light_mps./σ_v)).^2)
        flux .*= 1 .- d2.*exp.(-0.5.*((λs.-λ2)./λ2.*(EchelleCCFs.speed_of_light_mps./σ_v)).^2)
        @test_nowarn ccf_1D(λs, flux, ccfpl)
        var = flux./10^5
        @test_nowarn ccf_1D(λs, flux, var, ccfpl)
        (ccf, ccf_var) = ccf_1D(λs, flux, var, ccfpl)

        @test_nowarn measure_rv_from_ccf(v_grid,ccf, alg=MeasureRvFromCCFQuadratic() )
        @test_nowarn measure_rv_from_ccf(v_grid,ccf,ccf_var, alg=MeasureRvFromCCFQuadratic() )
        vfit = measure_rv_from_ccf(v_grid,ccf, alg=MeasureRvFromCCFQuadratic() )    # TODO: Check Why is this doing so much worse?
        @test vfit.rv ≈ 0  atol = 0.3   # TODO: tune tolerance better
        vfit = measure_rv_from_ccf(v_grid,ccf,ccf_var, alg=MeasureRvFromCCFQuadratic() )
        @test vfit.rv ≈ 0  atol = 0.3   # TODO: tune tolerance better

        @test_nowarn measure_rv_from_ccf(v_grid,ccf, alg=MeasureRvFromMinCCF() )
        vfit = measure_rv_from_ccf(v_grid,ccf, alg=MeasureRvFromMinCCF())
        @test vfit.rv ≈ 0  atol = 0.1   # TODO: tune tolerance better
        vfit = measure_rv_from_ccf(v_grid,ccf,ccf_var, alg=MeasureRvFromMinCCF())
        @test vfit.rv ≈ 0  atol = 0.5   # TODO: tune tolerance better

        #=
        @test_nowarn measure_rv_from_ccf(v_grid,ccf, alg=MeasureRvFromCCFCentroid() )
        vfit = measure_rv_from_ccf(v_grid,ccf)
        @test vfit.rv ≈ 0  atol = 500   # TODO: tune tolerance better
        =#
    end
end


text:
"""
    init_model_breathing!(
        J :: Vector{Vector{Float64}},
        l :: Lattice
        ) :: Nothing

Init Heisenberg model on a breathing pyrochlore or kagome lattice by overwriting the respective bonds.
Here, J[n] is the coupling to the n-th nearest neighbor (Euclidean norm). J[1] has to be an array of
length 2, specifying the breathing, anisotropic nearest neighbor couplings.
If there are m symmetry inequivalent n-th nearest neighbors (n > 1), these are
* uniformly initialized if J[n] is a single value
* initialized in ascending bond distance from the origin, if J[n] is an array of length m
"""
function init_model_breathing!(
    J :: Vector{Vector{Float64}},
    l :: Lattice
    ) :: Nothing

    @assert l.name in ["kagome", "pyrochlore"] "Breathing model requires Pyrochlore or Kagome lattice."
    @assert length(J[1]) == 2 "Breathing model needs two nearest neighbor couplings."

    # iterate over sites and add Heisenberg couplings to lattice bonds
    for i in eachindex(l.sites)
        for n in eachindex(J)
            # find n-th nearest neighbors
            nbs = get_nbs(n, l.sites[i], l.sites)

            # treat nearest neighbors according to breathing
            if n == 1
                for j in nbs
                    # if bond is within unit cell, it belongs to first kind of breathing coupling
                    if (l.sites[j].int - l.sites[i].int)[1 : 3] == [0, 0, 0]
                        add_bond!(J[1][1], l.bonds[i, j], 1, 1)
                        add_bond!(J[1][1], l.bonds[i, j], 2, 2)
                        add_bond!(J[1][1], l.bonds[i, j], 3, 3)
                    # if it crosses unit cell boundary, initalizes with second kind
                    else
                        add_bond!(J[1][2], l.bonds[i, j], 1, 1)
                        add_bond!(J[1][2], l.bonds[i, j], 2, 2)
                        add_bond!(J[1][2], l.bonds[i, j], 3, 3)
                    end
                end
            # uniform initialization for n-th nearest neighbor, if no further couplings provided
            elseif length(J[n]) == 1
                for j in nbs
                    add_bond!(J[n][1], l.bonds[i, j], 1, 1)
                    add_bond!(J[n][1], l.bonds[i, j], 2, 2)
                    add_bond!(J[n][1], l.bonds[i, j], 3, 3)
                end
            # initialize symmetry non-equivalent bonds interactions in ascending bond-length order
            else
                # get bond distances of neighbors
                dist = Int64[get_metric(l.sites[j], l.sites[i], l.uc) for j in nbs]

                # filter out classes of bond distances
                nbkinds = sort(unique(dist))

                # sanity check
                @assert length(J[n]) == length(nbkinds) "$(l.name) has $(length(nbkinds)) inequivalent $(n)-th nearest neighbors, but $(length(J[n])) couplings were supplied."

                for nk in eachindex(nbkinds)
                    # filter out neighbors with dist == nbkinds[nk]
                    nknbs = nbs[findall(x -> x == nbkinds[nk], dist)]

                    for j in nknbs
                        add_bond!(J[n][nk], l.bonds[i, j], 1, 1)
                        add_bond!(J[n][nk], l.bonds[i, j], 2, 2)
                        add_bond!(J[n][nk], l.bonds[i, j], 3, 3)
                    end
                end
            end
        end
    end

    return nothing
end

"""
    init_model_pyrochlore_breathing_c3!(
        J :: Vector{Vector{Float64}},
        l :: Lattice
        ) :: Nothing

Init Heisenberg model on a breathing pyrochlore lattice with broken C3 symmetry by overwriting the respective bonds.
Here, J[1] = [J1, J2, δ] specifies the breathing nearest-neighbor couplings J1 (for up tetrahedra),
J2 (for down tetrahedra) and δ, which quantifies the C3 breaking perturbation.
"""
function init_model_pyrochlore_breathing_c3!(
    J :: Vector{Vector{Float64}},
    l :: Lattice
    ) :: Nothing

    @assert l.name in ["pyrochlore"] "Model requires Pyrochlore lattice."
    @assert length(J[1]) == 3 "Model requires two nearest neighbor couplings and C3 breaking perturbation."
    @assert length(J) == 1 "Model supports only nearest neighbor couplings."

    # iterate over sites and add Heisenberg couplings to lattice bonds
    for i in eachindex(l.sites)
        # find nearest neighbors
        nbs = get_nbs(1, l.sites[i], l.sites)

        # treat nearest neighbors according to breathing anisotropy
        for j in nbs
            # get basis indices
            idxs = (l.sites[j].int[4], l.sites[i].int[4])

            # set coupling to J1 (± δ) for up tetrahedra
            if (l.sites[j].int - l.sites[i].int)[1 : 3] == [0, 0, 0]
                # add δ on bonds connecting basis sites 1 and 2
                if idxs == (1, 2) || idxs == (2, 1)
                    add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 1, 1)
                    add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 2, 2)
                    add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 3, 3)
                # add δ on bonds connecting basis sites 3 and 4
                elseif idxs == (3, 4) || idxs == (4, 3)
                    add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 1, 1)
                    add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 2, 2)
                    add_bond!(J[1][1] + J[1][3], l.bonds[i, j], 3, 3)
                # subtract δ for all other bonds
                else
                    add_bond!(J[1][1] - J[1][3], l.bonds[i, j], 1, 1)
                    add_bond!(J[1][1] - J[1][3], l.bonds[i, j], 2, 2)
                    add_bond!(J[1][1] - J[1][3], l.bonds[i, j], 3, 3)
                end
            # set coupling to J2 (± δ) for down tetrahedra
            else
                # add δ on bonds connecting basis sites 1 and 2
                if idxs == (1, 2) || idxs == (2, 1)
                    add_bond!(J[1][2] + J[1][3], l.bonds[i, j], 1, 1)
                    add_bond!(J[1][2] + J[1][3], l.bonds[i, j], 2, 2)
                    add_bond!(J[1][2] + J[1][3], l.bonds[i, j], 3, 3)
                # add δ on bonds connecting basis sites 3 and 4
                elseif idxs == (3, 4) || idxs == (4, 3)
                    add_bond!(J[1][2] + J[1][3], l.bonds[i, j], 1, 1)
                    add_bond!(J[1][2] + J[1][3], l.bonds[i, j], 2, 2)
                    add_bond!(J[1][2] + J[1][3], l.bonds[i, j], 3, 3)
                # ignore δ for all other bonds
                else
                    add_bond!(J[1][2] - J[1][3], l.bonds[i, j], 1, 1)
                    add_bond!(J[1][2] - J[1][3], l.bonds[i, j], 2, 2)
                    add_bond!(J[1][2] - J[1][3], l.bonds[i, j], 3, 3)
                end
            end
        end
    end

    return nothing
end