text:
"""
A `Knockoff` is an `AbstractMatrix`, essentially the matrix [X X̃] of
concatenating X̃ to X. If `A` is a `Knockoff`, `A` behaves like a regular matrix
and can be inputted into any function that supports inputs of type `AbstractMatrix`.
Basic operations like @view(A[:, 1]) are supported.
"""
struct Knockoff{T} <: AbstractMatrix{T}
    X::Matrix{T} # n × p design matrix
    X̃::Matrix{T} # n × p knockoff of X
    s::Vector{T} # p × 1 vector. Diagonal(s) and 2Σ - Diagonal(s) are both psd
    Σ::Matrix{T} # p × p gram matrix X'X
    Σinv::Matrix{T} # p × p inv(X'X)
end

function knockoff(X::AbstractMatrix{T}, X̃::AbstractMatrix{T}) where T
    Knockoff(X, X̃, T[], Matrix{T}(undef, 0, 0), Matrix{T}(undef, 0, 0))
end

Base.size(A::Knockoff) = size(A.X, 1), 2size(A.X, 2)
Base.eltype(A::Knockoff) = eltype(A.X)
function Base.getindex(A::Knockoff, i::Int)
    n, p = size(A.X)
    i ≤ n * p ? getindex(A.X, i) : getindex(A.X̃, i - n * p)
end
function Base.getindex(A::Knockoff, i::Int, j::Int)
    n, p = size(A.X)
    j ≤ p ? getindex(A.X, i, j) : getindex(A.X̃, i, j - p)
end
function LinearAlgebra.mul!(C::AbstractMatrix, A::Knockoff, B::AbstractMatrix)
    p = size(A.X, 2)
    mul!(C, A.X, @view(B[1:p, :]))
    mul!(C, A.X̃, @view(B[p+1:end, :]), 1.0, 1.0)
end
function LinearAlgebra.mul!(c::AbstractVector, A::Knockoff, b::AbstractVector)
    p = size(A.X, 2)
    mul!(c, A.X, @view(b[1:p]))
    mul!(c, A.X̃, @view(b[p+1:end]), 1.0, 1.0)
end

function normalize_col!(X::AbstractMatrix)
    n, p = size(X)
    @inbounds for x in eachcol(X)
        μi = mean(x)
        xnorm = norm(x)
        @simd for i in eachindex(x)
            x[i] = (x[i] - μi) / xnorm
        end
    end
    return X
end

# 1 state of a markov chain
struct GenotypeState
    a::Int # int between 1 and K
    b::Int # int between 1 and K
end

"""
Genotype states are index pairs (ka, kb) where ka, kb is unordered haplotype 1 and 2.
If there are K=5 haplotype motifs, then the 15 possible genotype states and their index are

(1, 1) = 1
(1, 2) = 2     (2, 2) = 6
(1, 3) = 3     (2, 3) = 7     (3, 3) = 10
(1, 4) = 4     (2, 4) = 8     (3, 4) = 11     (4, 4) = 13
(1, 5) = 5     (2, 5) = 9     (3, 5) = 12     (4, 5) = 14     (5, 5) = 15
"""
struct MarkovChainTable
    K::Int # number of states
    index_to_pair::Vector{GenotypeState} # index_to_pair[1] = (1, 1), index_to_pair[2] = (1, 2)...etc
end
function MarkovChainTable(K::Int)
    table = GenotypeState[]
    for a in 1:K, b in a:K
        push!(table, GenotypeState(a, b))
    end
    return MarkovChainTable(K, table)
end

Base.enumerate(mc::MarkovChainTable) = enumerate(mc.index_to_pair)
statespace(mc::MarkovChainTable) = ((mc.K + 1) * mc.K) >> 1
function pair_to_index(mc::MarkovChainTable, a::Int, b::Int)
    statespace(mc) - ((mc.K-a+2)*(mc.K-a+1))>>1 + (b - a + 1)
end
function index_to_pair(mc::MarkovChainTable, i::Int)
    return mc.index_to_pair[i].a, mc.index_to_pair[i].b
end
# mc = MarkovChainTable(5)
# index_to_pair(mc, 10)


text:
module JuliaPortTest

"""
Return a binary "ok" in case of a good input and cause an error otherwise.
"""
function check_function(data::AbstractVector{UInt8})::Vector{UInt8}
    if data == b"good"
        return b"ok"
    else
        throw(DomainError)
    end
end

end

ext: jl

max_stars_repo_name: cursorinsight/julia-port

max_issues_repo_name: Laymer/julia-port

max_forks_repo_name: Laymer/julia-port


text:
"""
    AbstractDimStack

Abstract supertype for dimensional stacks.

These have multiple layers of data, but share dimensions.

Notably, their behaviour lies somewhere between a `DimArray` and a `NamedTuple`:

- indexing with a `Symbol` as in `dimstack[:symbol]` returns a `DimArray` layer.
- iteration amd `map` are apply over array layers, as indexed with a `Symbol`.
- `getindex` and many base methods are applied as for `DimArray` - to avoid the need
    to allways use `map`.

This design gives very succinct code when working with many-layered, mixed-dimension objects.
But it may be jarring initially - the most surprising outcome is that `dimstack[1]` will return
a `NamedTuple` of values for the first index in all layers, while `first(dimstack)` will return
the first value of the iterator - the `DimArray` for the first layer.

See [`DimStack`](@ref) for the concrete implementation.
Most methods are defined on the abstract type.

To extend `AbstractDimStack`, implement [`rebuild`](@ref) and
[`rebuild_from_arrays`](@ref).
"""
abstract type AbstractDimStack{L} end

data(s::AbstractDimStack) = s.data
dims(s::AbstractDimStack) = s.dims
refdims(s::AbstractDimStack) = s.refdims
metadata(s::AbstractDimStack) = s.metadata

layerdims(s::AbstractDimStack) = s.layerdims
layerdims(s::AbstractDimStack, key::Symbol) = dims(s, layerdims(s)[key])
layermetadata(s::AbstractDimStack) = s.layermetadata
layermetadata(s::AbstractDimStack, key::Symbol) = layermetadata(s)[key]

Base.parent(s::AbstractDimStack) = s.data
@inline Base.keys(s::AbstractDimStack) = keys(data(s))
Base.haskey(s::AbstractDimStack, k) = k in keys(s)
Base.values(s::AbstractDimStack) = values(layers(s))
Base.first(s::AbstractDimStack) = s[first(keys(s))]
Base.last(s::AbstractDimStack) = s[last(keys(s))]
# Only compare data and dim - metadata and refdims can be different
Base.:(==)(s1::AbstractDimStack, s2::AbstractDimStack) =
    data(s1) == data(s2) && dims(s1) == dims(s2) && layerdims(s1) == layerdims(s2)
Base.length(s::AbstractDimStack) = length(keys(s))
Base.size(s::AbstractDimStack) = map(length, dims(s))
Base.size(A::AbstractDimStack, dims::DimOrDimType) = size(A, dimnum(A, dims))
Base.size(A::AbstractDimStack, dims::Integer) = size(A)[dims]
Base.axes(s::AbstractDimStack) = map(first ∘ axes, dims(s))
Base.axes(A::AbstractDimStack, dims::DimOrDimType) = axes(A, dimnum(A, dims))
Base.axes(A::AbstractDimStack, dims::Integer) = axes(A)[dims]
Base.iterate(s::AbstractDimStack, args...) = iterate(layers(s), args...)
Base.read(s::AbstractDimStack) = map(read, s)

function rebuild(
    s::AbstractDimStack, data, dims=dims(s), refdims=refdims(s),
    layerdims=layerdims(s), metadata=metadata(s), layermetadata=layermetadata(s)
)
    basetypeof(s)(data, dims, refdims, layerdims, metadata, layermetadata)
end

function rebuildsliced(f::Function, s::AbstractDimStack, layers, I)
    layerdims = map(basedims, layers)
    dims, refdims = slicedims(f, s, I)
    rebuild(s; data=map(parent, layers), dims=dims, refdims=refdims, layerdims=layerdims)
end

"""
    rebuild_from_arrays(s::AbstractDimStack, das::NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractDimArray}}}; kw...)

Rebuild an `AbstractDimStack` from a `NamedTuple` of `AbstractDimArray`
and an existing stack.

# Keywords

Keywords are simply the fields of the stack object:

- `data`
- `dims`
- `refdims`
- `metadata`
- `layerdims`
- `layermetadata`
"""
function rebuild_from_arrays(
    s::AbstractDimStack, das::Tuple{Vararg{<:AbstractDimArray}}; kw...
)
    rebuild_from_arrays(s, NamedTuple{keys(s)}(das); kw...)
end
function rebuild_from_arrays(
    s::AbstractDimStack, das::NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractDimArray}}};
    refdims=refdims(s),
    metadata=DD.metadata(s),
    data=map(parent, das),
    dims=DD.combinedims(das...),
    layerdims=map(DD.basedims, das),
    layermetadata=map(DD.metadata, das),
)
    rebuild(s; data, dims, refdims, layerdims, metadata, layermetadata)
end

function layers(s::AbstractDimStack{<:NamedTuple{Keys}}) where Keys
    NamedTuple{Keys}(map(K -> s[K], Keys))
end

Adapt.adapt_structure(to, s::AbstractDimStack) = map(A -> Adapt.adapt(to, A), s)

# Dipatch on Tuple of Dimension, and map
for func in (:index, :lookup, :metadata, :sampling, :span, :bounds, :locus, :order)
    @eval ($func)(s::AbstractDimStack, args...) = ($func)(dims(s), args...)
end

"""
    DimStack <: AbstractDimStack

    DimStack(data::AbstractDimArray...)
    DimStack(data::Tuple{Vararg{<:AbstractDimArray}})
    DimStack(data::NamedTuple{Keys,Vararg{<:AbstractDimArray}})
    DimStack(data::NamedTuple, dims::DimTuple; metadata=NoMetadata())

DimStack holds multiple objects sharing some dimensions, in a `NamedTuple`.

Notably, their behaviour lies somewhere between a `DimArray` and a `NamedTuple`:

- indexing with a `Symbol` as in `dimstack[:symbol]` returns a `DimArray` layer.
- iteration amd `map` are apply over array layers, as indexed with a `Symbol`.
- `getindex` or `view` with `Int`, `Dimension`s or `Selector`s that resolve to `Int` will
    return a `NamedTuple` of values from each layer in the stack.
    This has very good performace, and avoids the need to always use `map`.
- `getindex` or `view` with a `Vector` or `Colon` will return another `DimStack` where
    all data layers have been sliced.
- `setindex!` must pass a `Tuple` or `NamedTuple` maching the layers.
- many base and `Statistics` methods (`sum`, `mean` etc) will work as for a `DimArray`
    again removing the need to use `map`.

For example, here we take the mean over the time dimension for all layers :

```julia
mean(mydimstack; dims=Ti)
```

And this equivalent to:

```julia
map(A -> mean(A; dims=Ti), mydimstack)
```

This design gives succinct code when working with many-layered, mixed-dimension objects.

But it may be jarring initially - the most surprising outcome is that `dimstack[1]` will return
a `NamedTuple` of values for the first index in all layers, while `first(dimstack)` will return
the first value of the iterator - the `DimArray` for the first layer.

`DimStack` can be constructed from multiple `AbstractDimArray` or a `NamedTuple`
of `AbstractArray` and a matching `dims` tuple.

Most `Base` and `Statistics` methods that apply to `AbstractArray` can be used on
all layers of the stack simulataneously. The result is a `DimStack`, or
a `NamedTuple` if methods like `mean` are used without `dims` arguments, and
return a single non-array value.

## Example

```jldoctest
julia> using DimensionalData

julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0];

julia> dimz = (X([:a, :b]), Y(10.0:10.0:30.0))
X Symbol[a, b],
Y 10.0:10.0:30.0

julia> da1 = DimArray(1A, dimz; name=:one);

julia> da2 = DimArray(2A, dimz; name=:two);

julia> da3 = DimArray(3A, dimz; name=:three);

julia> s = DimStack(da1, da2, da3);

julia> s[At(:b), At(10.0)]
(one = 4.0, two = 8.0, three = 12.0)

julia> s[X(At(:a))] isa DimStack
true
```

"""
struct DimStack{L,D<:Tuple,R<:Tuple,LD<:NamedTuple,M,LM<:NamedTuple} <: AbstractDimStack{L}
    data::L
    dims::D
    refdims::R
    layerdims::LD
    metadata::M
    layermetadata::LM
end
DimStack(das::AbstractDimArray...; kw...) = DimStack(das; kw...)
function DimStack(das::Tuple{Vararg{<:AbstractDimArray}}; kw...)
    DimStack(NamedTuple{uniquekeys(das)}(das); kw...)
end
function DimStack(das::NamedTuple{<:Any,<:Tuple{Vararg{<:AbstractDimArray}}};
    data=map(parent, das), dims=combinedims(das...), layerdims=map(basedims, das),
    refdims=(), metadata=NoMetadata(), layermetadata=map(DD.metadata, das)
)
    DimStack(data, dims, refdims, layerdims, metadata, layermetadata)
end
# Same sized arrays
function DimStack(data::NamedTuple, dims::Tuple;
    refdims=(), metadata=NoMetadata(), layermetadata=map(_ -> NoMetadata(), data)
)
    all(map(d -> axes(d) == axes(first(data)), data)) || _stack_size_mismatch()
    layerdims = map(_ -> basedims(dims), data)
    DimStack(data, format(dims, first(data)), refdims, layerdims, metadata, layermetadata)
end

@noinline _stack_size_mismatch() = throw(ArgumentError("Arrays must have identical axes. For mixed dimensions, use DimArrays`"))


text:
using Statistics
using Oceananigans.Architectures: AbstractGPUArchitecture

#####
##### In place reductions!
#####

for function_name in (:sum, :prod, :maximum, :minimum, :all, :any)

    function_name! = Symbol(function_name, '!')

    # Unwrap ReducedField to a view over interior nodes:
    @eval begin
        Base.$(function_name!)(f::Function, r::AbstractReducedField, a::AbstractArray; kwargs...) = Base.$(function_name!)(f, interior(r), a; kwargs...)
        Base.$(function_name!)(r::AbstractReducedField, a::AbstractArray; kwargs...) = Base.$(function_name!)(identity, interior(r), a; kwargs...)
    end
end

function Statistics.norm(a::AbstractField)
    arch = a.architecture
    grid = a.grid

    r = zeros(arch, grid, 1)

    Base.mapreducedim!(x -> x * x, +, r, a)

    return CUDA.@allowscalar sqrt(r[1])
end


text:
import Base: *, ∈, isempty

export LinearMap,
       an_element,
       constraints_list

"""
    LinearMap{N<:Real, S<:LazySet{N}, NM, MAT<:AbstractMatrix{NM}} <: LazySet{N}

Type that represents a linear transformation ``M⋅S`` of a convex set ``S``.

### Fields

- `M` -- matrix/linear map
- `X` -- convex set

### Notes

This type is parametric in the elements of the linear map, `NM`, which is
independent of the numeric type of the wrapped set (`N`).
Typically `NM = N`, but there may be exceptions, e.g., if `NM` is an interval
that holds numbers of type `N`, where `N` is a floating point number type such
as `Float64`.

### Examples

For the examples we create a ``3×2`` matrix and two unit squares, one of them
being two-dimensional and the other one being one-dimensional.

```jldoctest constructors
julia> A = [1 2; 1 3; 1 4]; X = BallInf([0, 0], 1); Y = BallInf([0], 1);
```

The function ``*`` can be used as an alias to construct a `LinearMap` object.

```jldoctest constructors
julia> lm = LinearMap(A, X)
LinearMap{Int64,BallInf{Int64},Int64,Array{Int64,2}}([1 2; 1 3; 1 4], BallInf{Int64}([0, 0], 1))

julia> lm2 = A * X
LinearMap{Int64,BallInf{Int64},Int64,Array{Int64,2}}([1 2; 1 3; 1 4], BallInf{Int64}([0, 0], 1))

julia> lm == lm2
true
```

For convenience, `A` does not need to be a matrix but we also allow to use
vectors (interpreted as an ``n×1`` matrix) and `UniformScaling`s resp. scalars
(interpreted as a scaling, i.e., a scaled identity matrix).
Scaling by ``1`` is ignored.

```jldoctest constructors
julia> using LinearAlgebra: I

julia> [2, 3] * Y
LinearMap{Int64,BallInf{Int64},Int64,Array{Int64,2}}([2; 3], BallInf{Int64}([0], 1))

julia> lm3 = 2 * X
LinearMap{Int64,BallInf{Int64},Int64,SparseArrays.SparseMatrixCSC{Int64,Int64}}(
  [1, 1]  =  2
  [2, 2]  =  2, BallInf{Int64}([0, 0], 1))

julia> 2I * X == lm3
true

julia> 1I * X == X
true
```

Applying a linear map to a `LinearMap` object combines the two maps into a
single `LinearMap` instance.
Again we can make use of the conversion for convenience.

```jldoctest constructors
julia> B = transpose(A); B * lm
LinearMap{Int64,BallInf{Int64},Int64,Array{Int64,2}}([3 9; 9 29], BallInf{Int64}([0, 0], 1))

julia> B = [3, 4, 5]; B * lm
LinearMap{Int64,BallInf{Int64},Int64,Array{Int64,2}}([12 38], BallInf{Int64}([0, 0], 1))

julia> B = 2; B * lm
LinearMap{Int64,BallInf{Int64},Int64,Array{Int64,2}}([2 4; 2 6; 2 8], BallInf{Int64}([0, 0], 1))
```

The application of a `LinearMap` to a `ZeroSet` or an `EmptySet` is simplified
automatically.

```jldoctest constructors
julia> A * ZeroSet{Int}(2)
ZeroSet{Int64}(3)

julia> A * EmptySet{Int}()
EmptySet{Int64}()
```
"""
struct LinearMap{N<:Real, S<:LazySet{N},
                 NM, MAT<:AbstractMatrix{NM}} <: LazySet{N}
    M::MAT
    X::S

    # default constructor with dimension match check
    function LinearMap(M::MAT, X::S) where {N<:Real, S<:LazySet{N}, NM,
                                            MAT<:AbstractMatrix{NM}}
        @assert dim(X) == size(M, 2) "a linear map of size $(size(M)) cannot " *
            "be applied to a set of dimension $(dim(X))"
        return new{N, S, NM, MAT}(M, X)
    end
end

isoperationtype(::Type{<:LinearMap}) = true

"""
```
    *(map::Union{AbstractMatrix, UniformScaling, AbstractVector, Real}, X::LazySet)
```

Alias to create a `LinearMap` object.

### Input

- `map` -- linear map
- `X`   -- convex set

### Output

A lazy linear map, i.e., a `LinearMap` instance.
"""
function *(map::Union{AbstractMatrix, UniformScaling, AbstractVector, Real}, X::LazySet)
    return LinearMap(map, X)
end

# convenience constructor from a vector
function LinearMap(v::AbstractVector, X::LazySet)
    n = dim(X)
    m = length(v)
    if n == m
        M = reshape(v, 1, length(v))
    else
        M = reshape(v, length(v), 1)
    end
    return LinearMap(M, X)
end

# convenience constructor from a UniformScaling
function LinearMap(M::UniformScaling{N}, X::LazySet) where {N<:Real}
    if M.λ == one(N)
        return X
    end
    return LinearMap(Diagonal(fill(M.λ, dim(X))), X)
end

# convenience constructor from a scalar
function LinearMap(α::Real, X::LazySet)
    n = dim(X)
    return LinearMap(sparse(α * I, n, n), X)
end

# combine two linear maps into a single linear map
function LinearMap(M::AbstractMatrix, lm::LinearMap)
    return LinearMap(M * lm.M, lm.X)
end

# disambiguations
function LinearMap(v::AbstractVector, lm::LinearMap)
    return invoke(LinearMap, Tuple{AbstractVector, LazySet}, v, lm)
end

function LinearMap(α::Real, lm::LinearMap)
    return invoke(LinearMap, Tuple{Real, LazySet}, α, lm)
end

# more efficient version
function LinearMap(M::UniformScaling{N}, lm::LinearMap) where {N<:Real}
    if M.λ == one(N)
        return lm
    end
    return LinearMap(M.λ * lm.M, lm.X)
end

# ZeroSet is "almost absorbing" for LinearMap (only the dimension changes)
function LinearMap(M::AbstractMatrix{N}, Z::ZeroSet{N}) where {N<:Real}
    @assert dim(Z) == size(M, 2) "a linear map of size $(size(M)) cannot " *
            "be applied to a set of dimension $(dim(Z))"
    return ZeroSet{N}(size(M, 1))
end

# EmptySet is absorbing for LinearMap
function LinearMap(M::AbstractMatrix{N}, ∅::EmptySet{N}) where {N<:Real}
    return ∅
end

"""
    dim(lm::LinearMap)::Int

Return the dimension of a linear map.

### Input

- `lm` -- linear map

### Output

The ambient dimension of the linear map.
"""
function dim(lm::LinearMap)::Int
    return size(lm.M, 1)
end

"""
    σ(d::AbstractVector{N}, lm::LinearMap{N}) where {N<:Real}

Return the support vector of the linear map.

### Input

- `d`  -- direction
- `lm` -- linear map

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the wrapped set.

### Notes

If ``L = M⋅S``, where ``M`` is a matrix and ``S`` is a convex set, it follows
that ``σ(d, L) = M⋅σ(M^T d, S)`` for any direction ``d``.
"""
function σ(d::AbstractVector{N}, lm::LinearMap{N}) where {N<:Real}
    return lm.M * σ(_At_mul_B(lm.M, d), lm.X)
end

"""
    ρ(d::AbstractVector{N}, lm::LinearMap{N}; kwargs...) where {N<:Real}

Return the support function of the linear map.

### Input

- `d`      -- direction
- `lm`     -- linear map
- `kwargs` -- additional arguments that are passed to the support function
              algorithm

### Output

The support function in the given direction.
If the direction has norm zero, the result depends on the wrapped set.

### Notes

If ``L = M⋅S``, where ``M`` is a matrix and ``S`` is a convex set, it follows
that ``ρ(d, L) = ρ(M^T d, S)`` for any direction ``d``.
"""
function ρ(d::AbstractVector{N}, lm::LinearMap{N}; kwargs...) where {N<:Real}
    return ρ(_At_mul_B(lm.M, d), lm.X; kwargs...)
end

"""
    isbounded(lm::LinearMap; cond_tol::Number=DEFAULT_COND_TOL)::Bool

Determine whether a linear map is bounded.

### Input

- `lm`       -- linear map
- `cond_tol` -- (optional) tolerance of matrix condition (used to check whether
                the matrix is invertible)

### Output

`true` iff the linear map is bounded.

### Algorithm

We first check if the matrix is zero or the wrapped set is bounded.
If not, we perform a sufficient check whether the matrix is invertible.
If the matrix is invertible, then the map being bounded is equivalent to the
wrapped set being bounded, and hence the map is unbounded.
Otherwise, we check boundedness via [`isbounded_unit_dimensions`](@ref).
"""
function isbounded(lm::LinearMap; cond_tol::Number=DEFAULT_COND_TOL)::Bool
    if iszero(lm.M) || isbounded(lm.X)
        return true
    end
    if isinvertible(lm.M; cond_tol=cond_tol)
        return false
    end
    return isbounded_unit_dimensions(lm)
end

"""
    ∈(x::AbstractVector{N}, lm::LinearMap{N})::Bool where {N<:Real}

Check whether a given point is contained in a linear map of a convex set.

### Input

- `x`  -- point/vector
- `lm` -- linear map of a convex set

### Output

`true` iff ``x ∈ lm``.

### Algorithm

Note that ``x ∈ M⋅S`` iff ``M^{-1}⋅x ∈ S``.
This implementation does not explicitly invert the matrix, which is why it also
works for non-square matrices.

### Examples

```jldoctest
julia> lm = LinearMap([2.0 0.0; 0.0 1.0], BallInf([1., 1.], 1.));

julia> [5.0, 1.0] ∈ lm
false
julia> [3.0, 1.0] ∈ lm
true
```

An example with non-square matrix:
```jldoctest
julia> B = BallInf(zeros(4), 1.);

julia> M = [1. 0 0 0; 0 1 0 0]/2;

julia> [0.5, 0.5] ∈ M*B
true
```
"""
function ∈(x::AbstractVector{N}, lm::LinearMap{N})::Bool where {N<:Real}
    return lm.M \ x ∈ lm.X
end

"""
    an_element(lm::LinearMap{N})::Vector{N} where {N<:Real}

Return some element of a linear map.

### Input

- `lm` -- linear map

### Output

An element in the linear map.
It relies on the `an_element` function of the wrapped set.
"""
function an_element(lm::LinearMap{N})::Vector{N} where {N<:Real}
    return lm.M * an_element(lm.X)
end

"""
    isempty(lm::LinearMap)::Bool

Return if a linear map is empty or not.

### Input

- `lm` -- linear map

### Output

`true` iff the wrapped set is empty.
"""
function isempty(lm::LinearMap)::Bool
    return isempty(lm.X)
end

"""
    vertices_list(lm::LinearMap{N}; prune::Bool=true)::Vector{Vector{N}} where {N<:Real}

Return the list of vertices of a (polyhedral) linear map.

### Input

- `lm` -- linear map
- `prune` -- (optional, default: `true`) if true removes redundant vertices

### Output

A list of vertices.

### Algorithm

We assume that the underlying set `X` is polyhedral.
Then the result is just the linear map applied to the vertices of `X`.
"""
function vertices_list(lm::LinearMap{N}; prune::Bool=true)::Vector{Vector{N}} where {N<:Real}
    # for a zero map, the result is just the list containing the origin
    if iszero(lm.M)
        return [zeros(N, dim(lm))]
    end

    # collect low-dimensional vertices lists
    vlist_X = vertices_list(lm.X)

    # create resulting vertices list
    vlist = Vector{Vector{N}}()
    sizehint!(vlist, length(vlist_X))
    for v in vlist_X
        push!(vlist, lm.M * v)
    end

    return prune ? convex_hull(vlist) : vlist
end

"""
    constraints_list(lm::LinearMap{N}) where {N<:Real}

Return the list of constraints of a (polyhedral) linear map.

### Input

- `lm` -- linear map

### Output

The list of constraints of the linear map.

### Notes

We assume that the underlying set `X` is polyhedral, i.e., offers a method
`constraints_list(X)`.

### Algorithm

We fall back to a concrete set representation and apply `linear_map`.
"""
function constraints_list(lm::LinearMap{N}) where {N<:Real}
    return constraints_list(linear_map(lm.M, lm.X))
end

"""
    linear_map(M::AbstractMatrix{N}, lm::LinearMap{N}) where {N<:Real}

Return the linear map of a lazy linear map.

### Input

- `M`  -- matrix
- `lm` -- linear map

### Output

The polytope representing the linear map of the lazy linear map of a set.
"""
function linear_map(M::AbstractMatrix{N}, lm::LinearMap{N}) where {N<:Real}
     return linear_map(M * lm.M, lm.X)
end


text:
using Test
using GradientBoost.LossFunctions

struct DummyLossFunction <: LossFunction; end

y_examples = [
  [1,1],
  [1,0],
  [0,0],
  [-1,1],
  [-1,-1]
]
y_pred_examples = [
  [1,1],
  [1,1],
  [1,1],
  [1,1],
  [1,1]
]
bern_y_examples = [
  [1,1],
  [1,0],
  [0,1],
  [0,0],
  [0,0]
]

function test_loss(lf::LossFunction,
  y_examples, y_pred_examples, expected)

  for i = 1:size(y_examples, 1)
    @test loss(lf, y_examples[i], y_pred_examples[i]) ≈ expected[i]
  end
end
function test_negative_gradient(lf::LossFunction,
  y_examples, y_pred_examples, expected)

  for i = 1:size(y_examples, 1)
    @test negative_gradient(
      lf, y_examples[i], y_pred_examples[i]
    ) ≈ expected[i]
  end
end
function test_minimizing_scalar(lf::LossFunction,
  y_examples, expected)

  for i = 1:size(y_examples, 1)
    @test minimizing_scalar(lf, y_examples[i]) ≈ expected[i]
  end
end

@testset "Loss functions" begin
  @testset "not implemented functions throw an error" begin
    emp_vec = Array[]
    dlf = DummyLossFunction()

    @test_throws UndefVarError loss(dlf, emp_vec, emp_vec)
    @test_throws UndefVarError negative_gradient(dlf, emp_vec, emp_vec)
    @test_throws UndefVarError minimizing_scalar(dlf, emp_vec)
  end

  @testset "LeastSquares loss works" begin
    lf = LeastSquares()
    expected = [ 0.0, 0.5, 1.0, 2.0, 4.0 ]
    test_loss(lf, y_examples, y_pred_examples, expected)
  end
  @testset "LeastSquares negative_gradient works" begin
    lf = LeastSquares()
    expected = [ [0,0], [0,-1], [-1,-1], [-2,0], [-2,-2] ]
    test_negative_gradient(lf, y_examples, y_pred_examples, expected)
  end
  @testset "LeastSquares minimizing_scalar works" begin
    lf = LeastSquares()
    expected = [ 1.0, 0.5, 0.0, 0.0, -1.0 ]
    test_minimizing_scalar(lf, y_examples, expected)
  end

  @testset "LeastAbsoluteDeviation loss works" begin
    lf = LeastAbsoluteDeviation()
    expected = [ 0.0, 0.5, 1.0, 1.0, 2.0 ]
    test_loss(lf, y_examples, y_pred_examples, expected)
  end
  @testset "LeastAbsoluteDeviation negative_gradient works" begin
    lf = LeastAbsoluteDeviation()
    expected = [ [0,0], [0,-1], [-1,-1], [-1,0], [-1,-1] ]
    test_negative_gradient(lf, y_examples, y_pred_examples, expected)
  end
  @testset "LeastAbsoluteDeviation minimizing_scalar works" begin
    lf = LeastAbsoluteDeviation()
    expected = [ 1.0, 0.5, 0.0, 0.0, -1.0 ]
    test_minimizing_scalar(lf, y_examples, expected)
  end

  @testset "BinomialDeviance loss works" begin
    lf = BinomialDeviance()
    expected = [
      0.6265233750364456, 1.6265233750364456, 1.6265233750364456, 2.6265233750364456, 2.6265233750364456,
    ]
    test_loss(lf, bern_y_examples, y_pred_examples, expected)
  end
  @testset "BinomialDeviance negative_gradient works" begin
    lf = BinomialDeviance()
    expected = [
      [0.2689414213699951, 0.2689414213699951],
      [0.2689414213699951, -0.7310585786300049],
      [-0.7310585786300049, 0.2689414213699951],
      [-0.7310585786300049, -0.7310585786300049],
      [-0.7310585786300049, -0.7310585786300049]
    ]
    test_negative_gradient(lf, bern_y_examples, y_pred_examples, expected)
  end
  @testset "BinomialDeviance minimizing_scalar works" begin
    lf = BinomialDeviance()
    expected = [ Inf, 0.0, 0.0, -Inf, -Inf]
    test_minimizing_scalar(lf, bern_y_examples, expected)
  end
end


text:
JSONRPC.@dict_readable struct ReplRunCodeRequestParams <: JSONRPC.Outbound
    filename::String
    line::Int
    column::Int
    code::String
    mod::String
    showCodeInREPL::Bool
    showResultInREPL::Bool
    softscope::Bool
end

struct Frame
    path::String
    line::Int
end
Frame(st::Base.StackFrame) = Frame(fullpath(string(st.file)), st.line)

JSONRPC.@dict_readable struct ReplRunCodeRequestReturn <: JSONRPC.Outbound
    inline::String
    all::String
    stackframe::Union{Nothing,Vector{Frame}}
end
ReplRunCodeRequestReturn(inline, all) = ReplRunCodeRequestReturn(inline, all, nothing)

JSONRPC.@dict_readable mutable struct ReplWorkspaceItem <: JSONRPC.Outbound
    head::String
    id::Int
    haschildren::Bool
    lazy::Bool
    icon::String
    value::String
    canshow::Bool
    type::String
end

const repl_runcode_request_type = JSONRPC.RequestType("repl/runcode", ReplRunCodeRequestParams, ReplRunCodeRequestReturn)
const repl_interrupt_notification_type = JSONRPC.NotificationType("repl/interrupt", Nothing)
const repl_getvariables_request_type = JSONRPC.RequestType("repl/getvariables", Nothing, Vector{ReplWorkspaceItem})
const repl_getlazy_request_type = JSONRPC.RequestType("repl/getlazy", Int, Vector{ReplWorkspaceItem})
const repl_showingrid_notification_type = JSONRPC.NotificationType("repl/showingrid", String)
const repl_loadedModules_request_type = JSONRPC.RequestType("repl/loadedModules", Nothing, Vector{String})
const repl_isModuleLoaded_request_type = JSONRPC.RequestType("repl/isModuleLoaded", String, Bool)
const repl_startdebugger_notification_type = JSONRPC.NotificationType("repl/startdebugger", String)
const repl_showprofileresult_notification_type = JSONRPC.NotificationType("repl/showprofileresult", String)
const repl_showprofileresult_file_notification_type = JSONRPC.NotificationType("repl/showprofileresult_file", String)
const repl_toggle_plot_pane_notification_type = JSONRPC.NotificationType("repl/togglePlotPane", Bool)
const cd_notification_type = JSONRPC.NotificationType("repl/cd", String)
const activate_project_notification_type = JSONRPC.NotificationType("repl/activateProject", String)
const activate_project_from_dir_request = JSONRPC.RequestType("repl/activateProjectFromDir", String, Union{Nothing,String})


text:
#
# Computing optimal approximate posterior parameters.
#

function optimal_q(f::FiniteGP, y::AV{<:Real}, u::FiniteGP)
    U_y, U = cholesky(Symmetric(f.Σy)).U, cholesky(Symmetric(cov(u))).U
    B_εf, b_y = U' \ (U_y' \ cov(f, u))', U_y' \ (y - mean(f))
    Λ_ε = cholesky(Symmetric(B_εf * B_εf' + I))
    m_ε = Λ_ε \ (B_εf * b_y)
    return m_ε, Λ_ε, U
end
optimal_q(c::Observation, u::FiniteGP) = optimal_q(c.f, c.y, u)

# PPGP = "Pseudo-Point GP". Weird internal thing. Please don't use, or try to access.
# There are a couple of hacks associated with this object. Of particular note is the fact
# that cov(u, z) and cov(u, z, z′) return the precision rather than the cov, while
# cov(u, f, z, x) returns the covariance as expected. This odd behaviour is a hack that the
# routines to compute approximate posterior quantities know how to exploit properly. This
# essentially stems from the lack of a "precision" function, which we need in this case
# sometimes.
struct PPGP{Tm<:AV{<:Real}, TΛ<:Cholesky{<:Real}, Tz<:AV} <: AbstractGP
    m::Tm
    Λ::TΛ
    z::Tz
    n::Int
    gpc::GPC
    function PPGP{Tm, TΛ, Tz}(m::Tm, Λ::TΛ, z::Tz, gpc::GPC) where {Tm, TΛ, Tz}
        gp = new{Tm, TΛ, Tz}(m, Λ, z, next_index(gpc), gpc)
        gpc.n += 1
        return gp
    end
end
function PPGP(m::Tm, Λ::TΛ, z::Tz, gpc::GPC) where {Tm, TΛ, Tz}
    return PPGP{Tm, TΛ, Tz}(m, Λ, z, gpc)
end

mean_vector(u::PPGP, z::AV) = z === u.z ? u.m : throw(ArgumentError("Bad z"))
cov(u::PPGP, z::AV) = z === u.z ? u.Λ : throw(ArgumentError("Bad z"))
cov(u::PPGP, z::AV, z′::AV) = z === z′ ? cov(u, z) : throw(ArgumentError("Bad z"))
function cov(u::PPGP, u′::PPGP, z::AV, z′::AV)
    if u.gpc === u′.gpc && u.n === u′.n
        return cov(u, z, z′)
    else
        throw(ArgumentError("Undefined behaviour"))
    end
end

function cov(u::PPGP, f′::GP, z::AV, x′::AV)
    @assert u.z === z
    return zeros(length(z), length(x′))
end
function cov(f::GP, u′::PPGP, x::AV, z′::AV)
    @assert u′.z === z′
    return zeros(length(x), length(z′))
end

"""
    PseudoObs(u::AbstractGP, z::AV, m::AV{<:Real}, Λ::AM{<:Real}, U::AM)

Construct approximate observations of `u` at `z`, with mean and cov. `m` and `Λ`
"""
struct PseudoObs{Tu<:AbstractGP, Tû<:PPGP, Tz<:AV, TU<:AM, Tα<:AV}
    u::Tu
    û::Tû
    z::Tz
    U::TU
    α::Tα
end

function PseudoObs(u::AbstractGP, z::AV, m::AV, Λ::Cholesky, U::AM)
    return PseudoObs(u, PPGP(m, Λ, z, u.gpc), z, U, U \ m)
end

PseudoObs(c::Observation, u::FiniteGP) = PseudoObs(u.f, u.x, optimal_q(c, u)...)
PseudoObs(c::Observation, us::FiniteGP...) = PseudoObs(c, merge(us))

"""
    |(f::AbstractGP, ỹ::PseudoObs)

Condition `f` on approximate observations `ỹ`.
"""
|(f::AbstractGP, ỹ::PseudoObs) = CompositeGP((|, f, ỹ), f.gpc)

const approx_cond = Tuple{typeof(|), AbstractGP, PseudoObs}

mean_vector((_, f, ỹ)::approx_cond, x::AV) = mean(f(x)) + cov(f(x), ỹ.u(ỹ.z)) * ỹ.α

function cov((_, f, ỹ)::approx_cond, x::AV)
    u, z, U, Λ = ỹ.u, ỹ.z, ỹ.U, ỹ.û.Λ
    Ax = U' \ cov(u, f, z, x)
    return cov(f, x) - Ax' * Ax + Xt_invA_X(Λ, Ax)
end
function cov_diag((_, f, ỹ)::approx_cond, x::AV)
    u, z, U, Λ = ỹ.u, ỹ.z, ỹ.U, ỹ.û.Λ
    Ax = U' \ cov(u, f, z, x)
    return cov_diag(f, x) - diag_At_A(Ax) + diag_Xt_invA_X(Λ, Ax)
end

function cov((_, f, ỹ)::approx_cond, x::AV, x′::AV)
    u, z, U, Λ = ỹ.u, ỹ.z, ỹ.U, ỹ.û.Λ
    Ax = U' \ cov(u, f, z, x)
    Ax′ = U' \ cov(u, f, z, x′)
    return cov(f, x, x′) - Ax' * Ax′ + Xt_invA_Y(Ax, Λ, Ax′)
end
function cov_diag((_, f, ỹ)::approx_cond, x::AV, x′::AV)
    u, z, U, Λ = ỹ.u, ỹ.z, ỹ.U, ỹ.û.Λ
    Ax = U' \ cov(u, f, z, x)
    Ax′ = U' \ cov(u, f, z, x′)
    return cov_diag(f, x) - diag_At_B(Ax, Ax′) + diag_Xt_invA_Y(Ax, Λ, Ax′)
end

function cov((_, f, ỹ)::approx_cond, f′::AbstractGP, x::AV, x′::AV)
    u, z, U, Λ = ỹ.u, ỹ.z, ỹ.U, ỹ.û.Λ
    Ax = U' \ cov(u, f, z, x)
    Ax′ = U' \ cov(u, f′, z, x′)
    Kzx′ = cov(ỹ.û, f′, z, x′)
    if Kzx′ isa Cholesky
        return cov(f, f′, x, x′) - Ax' * Ax′ + Ax' * (Λ \ Matrix(U))
    else
        return cov(f, f′, x, x′) - Ax' * Ax′ + Xt_invA_Y(Ax, Λ, U' \ Kzx′)
    end
end

function cov(f::AbstractGP, (_, f′, ỹ)::approx_cond, x::AV, x′::AV)
    u, z, U, Λ = ỹ.u, ỹ.z, ỹ.U, ỹ.û.Λ
    Ax = U' \ cov(u, f, z, x)
    Ax′ = U' \ cov(u, f′, z, x′)
    Kzx = cov(ỹ.û, f, z, x)
    if Kzx isa Cholesky
        return cov(f, f′, x, x′) - Ax' * Ax′ + Matrix(U)' * (Λ \ Ax′)
    else
        return cov(f, f′, x, x′) - Ax' * Ax′ + Xt_invA_Y(U' \ Kzx, Λ, Ax′)
    end
end

function cov_diag(args::approx_cond, f′::AbstractGP, x::AV, x′::AV)
    return diag(cov(args, f′, x, x′))
end

function cov_diag(f::AbstractGP, args::approx_cond, x::AV, x′::AV)
    return diag(cov(f, args, x, x′))
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

### Parsing utilities

_directives(::Type{SimpleDateFormat{S,T}}) where {S,T} = T.parameters

character_codes(df::Type{SimpleDateFormat{S,T}}) where {S,T} =
    character_codes(_directives(df))
function character_codes(directives::Core.SimpleVector)
    letters = sizehint!(Char[], length(directives))
    for (i, directive) in enumerate(directives)
        if directive <: PeriodPart
            letter = first(directive.parameters)
            push!(letters, letter)
        end
    end
    return letters
end

genvar(t::Type{T}) where {T<:SimpleDate} = Symbol(lowercase(string(nameof(t))))
genvar(t::DataType) = Symbol(lowercase(string(nameof(t))))

"""
    tryparsenext_core(str::AbstractString, pos::Int, len::Int, pf::SimpleDateFormat, raise=false)

Parse the string according to the directives within the `SimpleDateFormat`. Parsing will start at
character index `pos` and will stop when all directives are used or we have parsed up to
the end of the string, `len`. When a directive cannot be parsed the returned value
will be `nothing` if `raise` is false otherwise an exception will be thrown.

If successful, return a 3-element tuple `(values, pos, num_parsed)`:
* `values::Tuple`: A tuple which contains a value
  for each `PeriodPart` within the `SimpleDateFormat` in the order
  in which they occur. If the string ends before we finish parsing all the directives
  the missing values will be filled in with default values.
* `pos::Int`: The character index at which parsing stopped.
* `num_parsed::Int`: The number of values which were parsed and stored within `values`.
  Useful for distinguishing parsed values from default values.
"""
@generated function tryparsenext_core(
    str::AbstractString,
    pos::Int,
    len::Int,
    pf::SimpleDateFormat,
    raise::Bool = false,
)
    directives = _directives(pf)
    letters = character_codes(directives)

    tokens = Type[CONVERSION_SPECIFIERS[letter] for letter in letters]
    value_names = Symbol[genvar(t) for t in tokens]
    value_defaults = Tuple(CONVERSION_DEFAULTS[t] for t in tokens)

    # Pre-assign variables to defaults. Allows us to use `@goto done` without worrying about
    # unassigned variables.
    assign_defaults = Expr[]
    for (name, default) in zip(value_names, value_defaults)
        push!(assign_defaults, quote
            $name = $default
        end)
    end

    vi = 1
    parsers = Expr[]
    for i = 1:length(directives)
        if directives[i] <: PeriodPart
            name = value_names[vi]
            vi += 1
            push!(parsers, quote
                pos > len && @goto done
                let val = tryparsenext(directives[$i], str, pos, len)
                    val === nothing && @goto error
                    $name, pos = val
                end
                num_parsed += 1
                directive_index += 1
            end)
        else
            push!(parsers, quote
                pos > len && @goto done
                let val = tryparsenext(directives[$i], str, pos, len)
                    val === nothing && @goto error
                    delim, pos = val
                end
                directive_index += 1
            end)
        end
    end
    return quote
        directives = pf.tokens

        num_parsed = 0
        directive_index = 1

        $(assign_defaults...)
        $(parsers...)

        pos > len || @goto error

        @label done
        return $(Expr(:tuple, value_names...)), pos, num_parsed

        @label error
        if raise
            if directive_index > length(directives)
                throw(ArgumentError("Found extra characters at the end of period string"))
            else
                d = directives[directive_index]
                throw(
                    ArgumentError(
                        "Unable to parse period. Expected directive $d at char $pos",
                    ),
                )
            end
        end
        return nothing
    end
end

"""
    tryparsenext_internal(::Type{<:TimeType}, str, pos, len, pf::SimpleDateFormat, raise=false)

Parse the string according to the directives within the `SimpleDateFormat`. The specified `TimeType`
type determines the type of and order of tokens returned. If the given `SimpleDateFormat` or string
does not provide a required token a default value will be used. When the string cannot be
parsed the returned value will be `nothing` if `raise` is false otherwise an exception will
be thrown.

If successful, returns a 2-element tuple `(values, pos)`:
* `values::Tuple`: A tuple which contains a value
  for each token as specified by the passed in type.
* `pos::Int`: The character index at which parsing stopped.
"""
@generated function tryparsenext_internal(
    ::Type{T},
    str::AbstractString,
    pos::Int,
    len::Int,
    pf::SimpleDateFormat,
    raise::Bool = false,
) where {T<:SimpleDate}
    letters = character_codes(pf)

    tokens = Type[CONVERSION_SPECIFIERS[letter] for letter in letters]
    value_names = Symbol[genvar(t) for t in tokens]

    output_tokens = CONVERSION_TRANSLATIONS[T]
    output_names = Symbol[genvar(t) for t in output_tokens]
    output_defaults = Tuple(CONVERSION_DEFAULTS[t] for t in output_tokens)

    # Pre-assign output variables to defaults. Ensures that all output variables are
    # assigned as the value tuple returned from `tryparsenext_core` may not include all
    # of the required variables.
    assign_defaults = Expr[
        quote
            $name = $default
        end for (name, default) in zip(output_names, output_defaults)
    ]

    # Unpacks the value tuple returned by `tryparsenext_core` into separate variables.
    value_tuple = Expr(:tuple, value_names...)

    return quote
        val = tryparsenext_core(str, pos, len, pf, raise)
        val === nothing && return nothing
        values, pos, num_parsed = val
        $(assign_defaults...)
        $value_tuple = values
        return $(Expr(:tuple, output_names...)), pos
    end
end

@inline function tryparsenext_base10(
    str::AbstractString,
    i::Int,
    len::Int,
    min_width::Int = 1,
    max_width::Int = 0,
)
    i > len && return nothing
    min_pos = min_width <= 0 ? i : i + min_width - 1
    max_pos = max_width <= 0 ? len : min(i + max_width - 1, len)
    d::Int64 = 0
    @inbounds while i <= max_pos
        c, ii = iterate(str, i)::Tuple{Char,Int}
        if '0' <= c <= '9'
            d = d * 10 + (c - '0')
        else
            break
        end
        i = ii
    end
    if i <= min_pos
        return nothing
    else
        return d, i
    end
end

@inline function tryparsenext_word(str::AbstractString, i, len, maxchars = 0)
    word_start, word_end = i, 0
    max_pos = maxchars <= 0 ? len : min(len, nextind(str, i, maxchars - 1))
    @inbounds while i <= max_pos
        c, ii = iterate(str, i)::Tuple{Char,Int}
        if isletter(c)
            word_end = i
        else
            break
        end
        i = ii
    end
    if word_end == 0
        return nothing
    else
        return SubString(str, word_start, word_end), i
    end
end

function Base.parse(
    ::Type{T},
    str::AbstractString,
    pf::SimpleDateFormat = default_format(T),
) where {T<:SimpleDate}
    pos, len = firstindex(str), lastindex(str)
    val = tryparsenext_internal(T, str, pos, len, pf, true)
    @assert val !== nothing
    values, endpos = val
    return T(values...)::T
end

#=
function Base.tryparse(::Type{T}, str::AbstractString, pf::SimpleDateFormat=default_format(T)) where T<:TimeType
    pos, len = firstindex(str), lastindex(str)
    res = tryparsenext_internal(T, str, pos, len, pf, false)
    res === nothing && return nothing
    values, endpos = res
    if validargs(T, values...) === nothing
        # TODO: validargs gets called twice, since it's called again in the T constructor
        return T(values...)::T
    end
    return nothing
end
=#

"""
    parse_components(str::AbstractString, pf::SimpleDateFormat) -> Array{Any}

Parse the string into its components according to the directives in the `SimpleDateFormat`.
Each component will be a distinct type, typically a subtype of Period. The order of the
components will match the order of the `PeriodPart` directives within the `SimpleDateFormat`. The
number of components may be less than the total number of `PeriodPart`.
"""
@generated function parse_components(str::AbstractString, pf::SimpleDateFormat)
    letters = character_codes(pf)
    tokens = Type[CONVERSION_SPECIFIERS[letter] for letter in letters]

    return quote
        pos, len = firstindex(str), lastindex(str)
        val = tryparsenext_core(str, pos, len, pf, true) #=raise=#
        @assert val !== nothing
        values, pos, num_parsed = val
        types = $(Expr(:tuple, tokens...))
        result = Vector{Any}(undef, num_parsed)
        for (i, typ) in enumerate(types)
            i > num_parsed && break
            result[i] = typ(values[i])  # Constructing types takes most of the time
        end
        return result
    end
end


text:
report(x) = report(STDOUT, x)

function report(io::IO, solver::LESolver; kwargs...)
    report(io, solver.sol; kwargs...)
end

function report(io::IO, sol::LESolution;
                convergence::Bool = true)
    LEs = lyapunov_exponents(sol)

    print_with_color(:blue, io, "Lyapunov Exponents Solution")
    if sol.converged
        print_with_color(:green, io, " (converged)")
    else
        print_with_color(:red, io, " (NOT converged)")
    end
    println(io)

    table = [
        ("#Orth.", sol.num_orth),
        ("#LEs", length(LEs)),
        ("LEs", LEs),
    ]
    for (name, value) in table
        print_with_color(:yellow, io, name)
        print(io, ": ")
        if value isa String
            print(io, value)
        else
            show(IOContext(io, :limit => true), value)
        end
        println(io)
    end

    if convergence
        report(io, sol.convergence)
    end
end

function report(io::IO, convergence::ConvergenceHistory;
                dim_lyap = min(10, length(convergence.errors)))

    if isempty(convergence.orth)
        print_with_color(:red, io, "NO convergence test is done!", bold=true)
        println(io)
        return
    end

    print_with_color(:blue, io, "Convergence")
    print(io, " #Orth.=$(convergence.orth[end])")
    print(io, " #Checks=$(length(convergence.orth))")
    if convergence.kinds[end] == UnstableConvError
        print(io, " [unstable]")
    else
        print(io, " [stable]")
    end
    println(io)

    if length(convergence.kinds) > 1
        print_with_color(:yellow, io, "Stability")
        print(io, ": ")
        for k in convergence.kinds
            if k == UnstableConvError
                print(io, "x")
            else
                print(io, ".")
            end
        end
        println(io)
    end

    print(io, " "^(length("LE$dim_lyap")), "  ",
          "      error", "   ",
          "  threshold")
    if convergence.kinds[end] == UnstableConvError
        print(io,
              "   variance",
              "   tail cov",
              " small tail?")
    end
    println(io)

    for i in 1:dim_lyap
        err = convergence.errors[i][end]
        th = convergence.thresholds[i][end]

        print(io, "LE$i")
        print(io, ": ")
        @printf(io, " %10.5g", err)
        if err < th
            print_with_color(:green, io, " < ")
        else
            print_with_color(:red, io, " > ")
        end
        @printf(io, " %10.5g", th)

        detail = convergence.details[i][end]
        if convergence.kinds[end] == UnstableConvError
            @printf(io, " %10.5g", detail.var)
            @printf(io, " %10.5g", detail.tail_cov)
            print(io, "  ")
            if detail.tail_ok
                print_with_color(:green, io, "yes")
            else
                print_with_color(:red, io, "no")
            end
        else
            print(io, "  ")
            compact_report(io, detail)
        end

        println(io)
    end
end

compact_report(io, _) = nothing
compact_report(io, ::FixedPointConvDetail) =
    print(io, "\"period\" < 1")
compact_report(io, detail::PeriodicConvDetail) =
    print(io, "\"period\": ", detail.period)
compact_report(io, ::NonNegativeAutoCovConvDetail) =
    print(io, "too short (NN)")
compact_report(io, ::NonPeriodicConvDetail) =
    print(io, "too short (NP)")


text:
module SoilMulch

# Packages [deps]
using Distributions
using DataFrames
using Statistics

export rainfall_poisson, soil_leakage, plant_transpiration, soil_evaporation
export cover_fraction, mulch_leakage, mulch_evaporation
export soil_water_balance, soil_mulch_water_balance, df_dt_day
export sol_swb, sol_mswb, sol_swb_crop, sol_mswb_crop

# ------------------------------------------------------------------------------
# Functions
# ------------------------------------------------------------------------------

# Rainfall Poisson -------------------------------------------------------------
"""
`rainfall_poisson(n, α, λ)`

Generates rainfall using Poisson process.

# Arguments
- `n::Integer`: number of events.
- `α::Float64`: mean rain depth.
- `λ::Float64`: rain event interval (1/d).

If rainfall needs to be similated in timescale smaller than daily, multiply `λ` by dt.

# Example
```
n = 30
α = 0.75
λ = 0.45
rain = rainfall_poisson(n, α, λ)
```
"""
function rainfall_poisson(n, α, λ)
    n = Int(n)
    rain = zeros(n)
    for i in 1:n
        r1 = rand(Uniform(0, 1))
        if r1 < λ
            r2 = rand(Uniform(0, 1))
            rain[i] = -α * log((1 - r2))
        else
            rain[i] = 0
        end
    end
    return rain
end

# Leakage ----------------------------------------------------------------------
"""
`soil_leakage(s, ks, β)`

Compute the soil leakage.

# Arguments
- `s::Float64`: soil moisture.
- `ks::Float64`: saturated hydraulic conductivity.
- `β::Float64`: exponent of hydraulic conductivity curve.
"""
function soil_leakage(s, ks, β)
    lk = ks * s^β
    return lk
end

# Transpiration ----------------------------------------------------------------
"""
`plant_transpiration(s, sw, sstar, ct, eto)`

Compute plant transpiration.

# Arguments
- `s::Float64`: soil moisture.
- `sw::Float64`: soil moisture at wilting point.
- `sstar::Float64`: soil moisture below the field capacity.
- `ct::Float64`: baseline transpiration factor.
- `eto::Float64`: reference evapotranspiration.
"""
function plant_transpiration(s, sw, sstar, ct, eto)
    # Compute soil moisture factor
    if s <= sw
        fs = 0
    elseif s > sw && s <= sstar
        fs = (s - sw)/(sstar - sw)
    else
        fs = 1
    end
    # Compute transpiration
    tp = fs * ct * eto
    return tp
end

# Soil evaporation -------------------------------------------------------------
"""
`soil_evaporation(s, sw, sstar, ct, eto)`

Compute soil evaporation.

# Arguments
- `s::Float64`: soil moisture.
- `sh::Float64`: soil moisture at wilting point.
- `ce::Float64`: baseline evaporation coefficient.
- `eto::Float64`: reference evapotranspiration.
"""
function soil_evaporation(s, sh, ce, eto)
    # Compute soil moisture factor
    if s <= sh
        fs = 0
    else
        fs = (s - sh)/(1 - sh)
    end
    ev = fs * ce * eto
    return ev
end

# Canopy interception ----------------------------------------------------------
"""
`canopy_interception(rain, cws, dt)`

Computes the amount of rain intercepted by the canopy.

# Arguments
- `rain::Float64`: rainfall.
- `cws::Float64`: canopy water storage capacity.
- `dt::Float64`: time interval.

Considering that `cws` is expressed in cm/day, but the rainfall might
be in a smaller timescale, `dt` will adjust `cws` to compute the interception
properly.
"""
function canopy_interception(rain, cws, dt)
    if rain <= (cws * dt)
        ci = rain
    else
        ci = (cws * dt)
    end
    return ci
end

# Cover fraction ---------------------------------------------------------------
"""
`cover_fraction(cf, t, tr, cmax, r, γ, tsen)`

Compute crop cover fraction.

# Arguments
- `cf::Float64`: crop cover fraction.
- `t::Float64`: time [days].
- `tr::Float64`: plant transpiration.
- `cmax::Float64`: maximum crop cover fraction.
- `r::Float64`: parameter of cover fraction model.
- `γ::Float64`: slope of increase senescence after tsen.
- `tsen::Float64`: senescence time.

If plant transpiration is not available, use ETo.
"""
function cover_fraction(cf, t, tr, cmax, r, γ, tsen)
    if (t - tsen) < 0
        hvs = 0
    else
        hvs = 1
    end
    dtc = (r * tr * cf * (1 - (cf/cmax))) - ((γ * (t - tsen) * hvs) * cf^2)
    return dtc
end

# Mulch leakage ----------------------------------------------------------------
"""
`mulch_leakage(ϑ, ϑh, km, g)`

Compute leakage from the mulch layer.

# Arguments
- `ϑ::Float64`: crop cover fraction.
- `ϑh::Float64`: threshold that stops mulch leakage and evaporation.
- `km::Float64`: percolation rate.
- `g::Float64`: exponent of mulching percolation.
"""
function mulch_leakage(ϑ, ϑh, km, g)
    if ϑ <= ϑh
        lm = 0
    else
        lm = km * ϑ^g
    end
    return lm
end

# Mulch evaporation ------------------------------------------------------------
"""
`mulch_evaporation(ϑ, ϑh, ce, eto)`

Compute crop cover fraction.

# Arguments
- `ϑ::Float64`: crop cover fraction.
- `ϑh::Float64`: threshold that stops mulch leakage and evaporation.
- `ce::Float64`: baseline evaporation coefficient.
- `eto::Float64`: reference evapotranspiration.
"""
function mulch_evaporation(ϑ, ϑh, ce, eto)
    if ϑ <= ϑh
        fs = 0
    else
        fs = (ϑ - ϑh) / (1 - ϑh)
    end
    em = fs * ce * eto
    return em
end

# Mulch evaporation reduction factor -------------------------------------------
"""
`mulch_res_zm(μ, zm)`

Computes soil evaporation resistance factor due to the mulching layer thickness.

# Arguments
- `μ::Float64`: an empirical constant to compute resistance factor.
- `zm::Float64`: mulching layer thickness [cm].

Reference: Jalota, S., Prihar, S., 1990. Effect of straw mulch on evaporation
reduction in relation to rates of mulching and evaporativity. J. Indian Soc.
Soil Sci. 38, 728–730.

Suggested value for `μ` = 0.11.
"""
function mulch_res_zm(μ, zm)
    etom = exp.(-μ * zm)
    return etom
end

# Soil water balance -----------------------------------------------------------
"""
`soil_water_balance(rain, eto, cf, s, sh, sw, sstar, ks, β, n, zr, rth, ce, ct, dt)`

Computes the soil water balance under no mulching conditions.

# Arguments
- `rain::Float64`: rainfall.
- `eto::Float64`: reference evapotranspiration.
- `cf::Float64`: crop cover fraction.
- `s::Float64`: soil moisture.
- `sh::Float64`: soil moisture at wilting point.
- `sw::Float64`: soil moisture at wilting point.
- `sstar::Float64`: soil moisture below the field capacity.
- `ks::Float64`: saturated hydraulic conductivity.
- `β::Float64`: exponent of hydraulic conductivity curve.
- `n::Float64`: soil porosity.
- `zr::Float64`: root depth.
- `cws::Float64`: canopy water storage capacity.
- `ce::Float64`: baseline evaporation coefficient.
- `ct::Float64`: baseline transpiration factor.
- `dt::Float64`: time interval.
"""
function soil_water_balance(rain, eto, cf, s, sh, sw, sstar, ks, β, n, zr, cws, ce, ct, dt)
    # Soil water storage capacity
    nzr = n * zr

    # Throughfall
    ci = canopy_interception(rain, cws, dt) * cf
    tf = rain - ci

    s = s + tf / nzr

    # Check if there is runoff
    if s > 1
        q = (s - 1) * nzr
        s = 1
    else
        q = 0
    end

    # Compute plant transpiration
    tr = plant_transpiration(s, sw, sstar, ct, eto) * cf * dt
    if s >= sstar
        trns = tr
    else
        trns = 0
    end
    s = s - tr / nzr

    # Compute soil evaporation
    ev = soil_evaporation(s, sh, ce, eto) * (1 - cf) * dt
    s = s - ev / nzr

    # Compute leakage
    lk = soil_leakage(s, ks, β) * dt
    s = s - lk / nzr

    # Write outputs
    y = (s = s, rain = rain, ci = ci, q = q, tr = tr, trns = trns, ev = ev, lk = lk)
    return y
end

# Soil water balance with mulching ---------------------------------------------

"""
`soil_mulch_water_balance(rain, eto, cf, ϑ, s, sh, sw, sstar, ks, β, n, zr, cws,
                          ce, ct, μ, ϕ, zm, ϑh, km, g, dt)`

Computes the mulching and soil water balance.

# Arguments
- `rain::Float64`: rainfall.
- `eto::Float64`: reference evapotranspiration.
- `cf::Float64`: crop cover fraction.
- `ϑ::Float64`: mulching moisture.
- `s::Float64`: soil moisture.
- `sh::Float64`: soil moisture at wilting point.
- `sw::Float64`: soil moisture at wilting point.
- `sstar::Float64`: soil moisture below the field capacity.
- `ks::Float64`: saturated hydraulic conductivity.
- `β::Float64`: exponent of hydraulic conductivity curve.
- `n::Float64`: soil porosity.
- `zr::Float64`: root depth.
- `cws::Float64`: canopy water storage capacity.
- `ce::Float64`: baseline evaporation coefficient.
- `ct::Float64`: baseline transpiration factor.
- `ϕ::Float64`: mulching porosity.
- `ϑh::Float64`: threshold that stops mulch leakage and evaporation.
- `km::Float64`: mulching percolation rate.
- `g::Float64`: exponent of mulching percolation.
- `dt::Float64`: time interval.
"""
function soil_mulch_water_balance(rain, eto, cf, ϑ, s, sh, sw, sstar, ks, β, n,
                                  zr, cws, ce, ct, μ, ϕ, zm, ϑh, km, g, dt)
    # Soil and mulch water storage capacity
    pzm = ϕ * zm
    nzr = n * zr

    # Throughfall
    ci = canopy_interception(rain, cws, dt) * cf
    tf = rain - ci

    ϑ = ϑ + tf / pzm

    if ϑ > 1
        q = (ϑ - 1) * pzm
        ϑ = 1
    else
        q = 0
    end

    # Mulching leakage
    lm = mulch_leakage(ϑ, ϑh, km, g) * dt
    ϑ = ϑ - lm / pzm

    # Mulching evaporation
    em = mulch_evaporation(ϑ, ϑh, ce, eto) * dt
    ϑ = ϑ - em / pzm

    # Soil layer
    s = s + lm / nzr

    if s > 1
        q = q + (s - 1) * nzr
        s = 1
    else
        q = q + 0
    end

    # Compute plant transpiration
    tr = plant_transpiration(s, sw, sstar, ct, eto) * cf * dt
    if s >= sstar
        trns = tr
    else
        trns = 0
    end
    s = s - tr / nzr

    # Compute soil evaporation
    etom = eto * mulch_res_zm(μ, zm)
    evr = soil_evaporation(s, sh, ce, etom) * (1 - cf) * dt - em
    if evr <= 0
        ev = 0
    else
        ev = evr
    end

    s = s - ev / nzr

    # Compute leakage
    lk = soil_leakage(s, ks, β) * dt
    s = s - lk / nzr

    # Write outputs
    y = (s = s, rain = rain, ci = ci, q = q, tr = tr, trns = trns, ev = ev,
         lk = lk, ϑ = ϑ, em = em, lm = lm)
    return y
end

# Convert results to daily timescale -------------------------------------------
"""
`df_dt_day(df)`

Convert model output from `dt` to daily timescale.

# Arguments
-`df::DataFrame`: a data frame.
"""
function df_dt_day(df)
    df1 = combine(groupby(df, :Days), :Rain=>sum, :CI=>sum, :Q=>sum, :s=>mean,
                  :Lk=>sum, :Es=>sum, :Tr=>sum, :Trns=>sum)
    # Check if there is mulching layer
    if "ϑ" in names(df)
        df2 = combine(groupby(df, :Days), :ϑ=>mean, :Lm=>sum, :Em=>sum)
        df1 = rightjoin(df1, df2, on=:Days)
    end

    # Check if there is crop cover fraction
    if "CF" in names(df)
        df3 = combine(groupby(df, :Days), :CF=>last)
        df1 = rightjoin(df1, df3, on=:Days)
    end

    # Get columns name and remove the operation
    cn = names(df1)
    cn2 = []
    for i in 1:length(cn)
        push!(cn2, split(cn[i], "_")[1])
    end

    rename!(df1, Symbol.(cn2))
    return df1
end

# ------------------------------------------------------------------------------
# Solutions
# ------------------------------------------------------------------------------

# Solution crop cover fraction -------------------------------------------------
"""
`cover_fraction_sol(cf, t, tr, cmax, r, γ, tsen)`

Compute numerical solution for the crop cover fraction.

# Arguments
- `c::Float64`: crop cover fraction.
- `t::Float64`: time [days].
- `tr::Float64`: plant transpiration.
- `cmax::Float64`: maximum crop cover fraction.
- `r::Float64`: parameter of cover fraction model.
- `γ::Float64`: slope of increase senescence after tsen.
- `tsen::Float64`: senescence time.

If plant transpiration is not available, use ETo.
"""
function cover_fraction_sol(cf, t, tr, cmax, r, γ, tsen)
    # Solver cover model for the growing season
    cday = zeros(length(eto))
    cday[1] = cf
    for i in 2:length(cday)
        cday[i] = cday[i-1] + cover_fraction(cday[i-1], t[i], tr[i], cmax, r, γ, tsen) * 1
    end
    return cday
end

# Solve soil water balance -----------------------------------------------------
"""
`sol_swb(rain, eto, cf, s, p)`

Solves the soil water balance function.

# Arguments
- `p::NamedTuple`: p is a named tuple with the parameters mostly considered as
constant.

For the parameters descriptions, please check the help of
`soil_water_balance` function.

"""
function sol_swb(rain, eto, cf, s, p)
    # Unzip parameter
    sh = p.sh
    sw = p.sw
    sstar = p.sstar
    ks = p.ks
    β = p.β
    n = p.n
    zr = p.zr
    cws = p.cws
    ce = p.ce
    ct = p.ct
    dt = p.dt

    nr = length(rain) + 1
    sres = zeros(nr)
    cires = zeros(nr)
    qres = zeros(nr)
    eres = zeros(nr)
    tres = zeros(nr)
    trnsres = zeros(nr)
    lres = zeros(nr)

    # Initial condition for soil moisture
    sres[1] = s

    # Create vector for time
    ndays = Int((nr - 1) * dt)
    days = [1:ndays;]
    days = repeat(days, inner = Int(1 / dt))
    days = vcat(0, days)
    daysc = [(1+dt):dt:(ndays+1);] .- dt

    # Adjust ETo
    if length(eto) == 1
        eto = repeat([eto], length(rain))
    else
        eto = eto
    end

    # Adjust cover fraction
    if length(cf) == 1
        cf = repeat([cf], length(rain))
    else
        cf = cf
    end

    for i in 1:length(rain)
        swb = soil_water_balance(rain[i], eto[i], cf[i], sres[i],
                                 sh, sw, sstar, ks, β,
                                 n, zr, cws, ce, ct, dt)
        sres[i+1] = swb.s
        cires[i+1] = swb.ci
        qres[i+1] = swb.q
        eres[i+1] = swb.ev
        tres[i+1] = swb.tr
        trnsres[i+1] = swb.trns
        lres[i+1] = swb.lk
    end
    df = DataFrame(Days = days, DaysC = vcat(0, daysc), Rain = vcat(0, rain),
                   CI = cires, Q = qres, s = sres, Lk = lres, Es = eres,
                   Tr = tres, Trns = trnsres)
    return df

end

# Solve soil water balance coupled with crop -----------------------------------
"""
`sol_swb_crop(rain, eto, cf, s, p)`

Solves the soil water balance coupled with the crop cover fraction.

# Arguments
- `p::NamedTuple`: p is a named tuple with the parameters mostly considered as
constant.

For the parameters descriptions, please check the help of
`soil_water_balance` function.

"""
function sol_swb_crop(rain, eto, cf, s, p)
    # Unzip parameter
    sh = p.sh
    sw = p.sw
    sstar = p.sstar
    ks = p.ks
    β = p.β
    n = p.n
    zr = p.zr
    cws = p.cws
    ce = p.ce
    ct = p.ct
    cmax = p.cmax
    r = p.r
    γ = p.γ
    tsen = p.tsen
    dt = p.dt

    nr = length(rain) + 1
    sres = zeros(nr)
    cires = zeros(nr)
    qres = zeros(nr)
    eres = zeros(nr)
    tres = zeros(nr)
    trnsres = zeros(nr)
    lres = zeros(nr)
    cfres = zeros(nr)

    # Initial condition for soil moisture
    sres[1] = s
    cfres[1] = cf

    # Create vector for time
    ndays = Int((nr - 1) * dt)
    days = [1:ndays;]
    days = repeat(days, inner = Int(1 / dt))
    days = vcat(0, days)
    daysc = [(1+dt):dt:(ndays+1);] .- dt

    # Adjust ETo
    if length(eto) == 1
        eto = repeat([eto], length(rain))
    else
        eto = eto
    end

    for i in 1:length(rain)
        swb = soil_water_balance(rain[i], eto[i], cfres[i], sres[i],
                                 sh, sw, sstar, ks, β,
                                 n, zr, cws, ce, ct, dt)
        sres[i+1] = swb.s
        cires[i+1] = swb.ci
        qres[i+1] = swb.q
        eres[i+1] = swb.ev
        tres[i+1] = swb.tr
        trnsres[i+1] = swb.trns
        lres[i+1] = swb.lk
        cfres[i+1] = cfres[i] + cover_fraction(cfres[i], daysc[i], swb.tr,
                                               cmax, r, γ, tsen)
    end
    df = DataFrame(Days = days, DaysC = vcat(0, daysc), Rain = vcat(0, rain),
                   CI = cires, Q = qres, s = sres, Lk = lres, Es = eres,
                   Tr = tres, Trns = trnsres, CF = cfres)
    return df

end

# Solve mulching soil water balance --------------------------------------------
"""
`sol_mswb(rain, eto, cf, ϑ, s, p)`

Solves the soil water balance function.

# Arguments
- `p::NamedTuple`: p is a named tuple with the parameters mostly considered as
constant.

For the parameters descriptions, please check the help of
`soil_mulch_water_balance` function.

"""
function sol_mswb(rain, eto, cf, ϑ, s, p)

    # Unzip parameter
    sh = p.sh
    sw = p.sw
    sstar = p.sstar
    ks = p.ks
    β = p.β
    n = p.n
    zr = p.zr
    cws = p.cws
    ce = p.ce
    ct = p.ct
    μ = p.μ
    ϕ = p.ϕ
    zm = p.zm
    ϑh = p.ϑh
    km = p.km
    g = p.g
    dt = p.dt

    # Vectors to receive outputs
    nr = length(rain) + 1
    sres = zeros(nr)
    ϑres = zeros(nr)
    cires = zeros(nr)
    qres = zeros(nr)
    eres = zeros(nr)
    tres = zeros(nr)
    trnsres = zeros(nr)
    lres = zeros(nr)
    lmres = zeros(nr)
    emres = zeros(nr)

    # Initial condition for soil moisture
    sres[1] = s
    ϑres[1] = ϑ

    # Create vector for time
    ndays = Int((nr - 1) * dt)
    days = [1:ndays;]
    days = repeat(days, inner = Int(1 / dt))
    days = vcat(0, days)
    daysc = [(1+dt):dt:(ndays+1);] .- dt

    # Adjust ETo
    if length(eto) == 1
        eto = repeat([eto], length(rain))
    else
        eto = eto
    end

    # Adjust cover fraction
    if length(cf) == 1
        cf = repeat([cf], length(rain))
    else
        cf = cf
    end

    for i in 1:length(rain)
        swb = soil_mulch_water_balance(rain[i], eto[i], cf[i], ϑres[i], sres[i],
                                       sh, sw, sstar, ks, β, n, zr, cws, ce, ct,
                                       μ, ϕ, zm, ϑh, km, g, dt)
        sres[i+1] = swb.s
        ϑres[i+1] = swb.ϑ
        cires[i+1] = swb.ci
        qres[i+1] = swb.q
        eres[i+1] = swb.ev
        tres[i+1] = swb.tr
        trnsres[i+1] = swb.trns
        lres[i+1] = swb.lk
        emres[i+1] = swb.em
        lmres[i+1] = swb.lm
    end
    df = DataFrame(Days = days, DaysC = vcat(0, daysc), Rain = vcat(0, rain),
                   CI = cires, Q = qres, s = sres, Lk = lres, Es = eres,
                   Tr = tres, Trns = trnsres, ϑ = ϑres, Lm = lmres, Em = emres)
    return df

end

# Solve mulching soil water balance coupled with crop --------------------------
"""
`sol_mswb_crop(rain, eto, cf, ϑ, s, p)`

Solves the mulchg and soil water balance coupled with the crop cover fraction.

# Arguments
- `p::NamedTuple`: p is a named tuple with the parameters mostly considered as
constant.

For the parameters descriptions, please check the help of
`soil_mulch_water_balance` function.

"""
function sol_mswb_crop(rain, eto, cf, ϑ, s, p)

    # Unzip parameter
    sh = p.sh
    sw = p.sw
    sstar = p.sstar
    ks = p.ks
    β = p.β
    n = p.n
    zr = p.zr
    cws = p.cws
    ce = p.ce
    ct = p.ct
    μ = p.μ
    ϕ = p.ϕ
    zm = p.zm
    ϑh = p.ϑh
    km = p.km
    g = p.g
    cmax = p.cmax
    r = p.r
    γ = p.γ
    tsen = p.tsen
    dt = p.dt

    # Vectors to receive outputs
    nr = length(rain) + 1
    sres = zeros(nr)
    ϑres = zeros(nr)
    cires = zeros(nr)
    qres = zeros(nr)
    eres = zeros(nr)
    tres = zeros(nr)
    trnsres = zeros(nr)
    lres = zeros(nr)
    lmres = zeros(nr)
    emres = zeros(nr)
    cfres = zeros(nr)

    # Initial condition for soil moisture
    sres[1] = s
    ϑres[1] = ϑ
    cfres[1] = cf

    # Create vector for time
    ndays = Int((nr - 1) * dt)
    days = [1:ndays;]
    days = repeat(days, inner = Int(1 / dt))
    days = vcat(0, days)
    daysc = [(1+dt):dt:(ndays+1);] .- dt

    # Adjust ETo
    if length(eto) == 1
        eto = repeat([eto], length(rain))
    else
        eto = eto
    end

    for i in 1:length(rain)
        swb = soil_mulch_water_balance(rain[i], eto[i], cfres[i], ϑres[i], sres[i],
                                       sh, sw, sstar, ks, β, n, zr, cws, ce, ct,
                                       μ, ϕ, zm, ϑh, km, g, dt)
        sres[i+1] = swb.s
        ϑres[i+1] = swb.ϑ
        cires[i+1] = swb.ci
        qres[i+1] = swb.q
        eres[i+1] = swb.ev
        tres[i+1] = swb.tr
        trnsres[i+1] = swb.trns
        lres[i+1] = swb.lk
        emres[i+1] = swb.em
        lmres[i+1] = swb.lm
        cfres[i+1] = cfres[i] + cover_fraction(cfres[i], daysc[i], swb.tr,
                                               cmax, r, γ, tsen)
    end
    df = DataFrame(Days = days, DaysC = vcat(0, daysc), Rain = vcat(0, rain),
                   CI = cires, Q = qres, s = sres, Lk = lres, Es = eres,
                   Tr = tres, Trns = trnsres, ϑ = ϑres, Lm = lmres, Em = emres,
                   CF = cfres)
    return df

end

end


text:
"""
Abstract decomposition of N-mode tensor.
"""
abstract type TensorDecomposition{T<:Number, N} end

"""
Returns relative error between the re-composed and the original tensors.
"""
rel_residue{T,N}(recomposed::StridedArray{T, N}, tensor::StridedArray{T, N}) = vecnorm(recomposed - tensor) / vecnorm(tensor)

"""
Returns relative error between re-composed and the original tensor.
"""
rel_residue{T,N}(decomp::TensorDecomposition{T, N}, tensor::StridedArray{T, N}) = rel_residue(compose(decomp), tensor)

"""
Returns relative tensor decomposition error, NaN if not available.
"""
rel_residue(decomp::TensorDecomposition) = get(decomp.props, :rel_residue, NaN)

_set_rel_residue(decomp::TensorDecomposition, error::Float64) = setindex!(decomp.props, error, :rel_residue)

_set_rel_residue{T,N}(decomp::TensorDecomposition{T,N}, tensor::StridedArray{T,N}) =
    _set_rel_residue(decomp, rel_residue(decomp, tensor))


text:
#
# Evaluate gamma (and if gamp is passed, evaluate gamma' as well) at the n nodes s
# for the forbidden region (regions 1 and 2). s0 is the lower saddle point.
# Used by find_endpoints() and contour(). No other definitions of this contour in the code.
#

function gammaforbidden!(gam::Vector{Complex{Float64}}, gamp::Vector{Complex{Float64}}, s0::Complex{Float64}, s::Vector{Float64}, n::Int)
    res,ims = reim(s0)
    @simd for i=1:n
        d = s[i] - res
        # region 1 (deep forbidden)
        if ims ≤ -M_PI_3
            @inbounds gam[i] = complex(s[i],THIRD*atan(d) - M_PI_3)
            @inbounds gamp[i] = complex(1.0,THIRD/(1+d*d))
        else
            ims = imshack(res,ims)
            d2 = d^2
            ex = exp(-d2)
            g0 = THIRD*(atan(d)-M_PI)
            @inbounds gam[i] = complex(s[i],ims + (g0 - ims)*(1.0-ex))
            @inbounds gamp[i] = complex(1.0,2.0*(g0 - ims)*d*ex + THIRD / (1.0 + d2) * (1.0-ex))
        end
    end
end

# assume if length n not given that vectors are of length 1 (with a slightly different signature)

function gammaforbidden!(gam::Vector{Complex{Float64}}, s0::Complex{Float64}, s::Float64)
    res,ims = reim(s0)
    d = s - res
    # region 1 (deep forbidden)
    if ims ≤ -M_PI_3
        gam[1] = complex(s,THIRD*atan(d) - M_PI_3)
    else
        ims = imshack(res,ims)
        d2 = d^2
        ex = exp(-d2)
        g0 = THIRD*(atan(d)-M_PI)
        gam[1] = complex(s,ims + (g0 - ims)*(1.0-ex))
    end
end

function imshack(res::Float64,ims::Float64)
    # hack to move ctr just below coalescing saddle at high E. Note exp(res) ~ sqrt(E).
    if ims > -0.1 ims -= max(0.0,min(0.7exp(-res),0.1)) end
    return ims
end


text:
module GLTF

using JSON3

export load, save

# constants
const ARRAY_BUFFER = 34962
const ELEMENT_ARRAY_BUFFER = 34963

const BYTE = 5120
const UNSIGNED_BYTE = 5121
const SHORT = 5122
const UNSIGNED_SHORT = 5123
const UNSIGNED_INT = 5125
const FLOAT = 5126

const POINTS = 0
const LINES = 1
const LINE_LOOP = 2
const LINE_STRIP = 3
const TRIANGLES = 4
const TRIANGLE_STRIP = 5
const TRIANGLE_FAN = 6

const NEAREST = 9728
const LINEAR = 9729
const NEAREST_MIPMAP_NEAREST = 9984
const LINEAR_MIPMAP_NEAREST = 9985
const NEAREST_MIPMAP_LINEAR = 9986
const LINEAR_MIPMAP_LINEAR = 9987
const CLAMP_TO_EDGE = 33071
const MIRRORED_REPEAT = 33648
const REPEAT = 10497

include("zvector.jl")
include("accessor.jl")
include("animation.jl")
include("asset.jl")
include("buffer.jl")
include("bufferView.jl")
include("camera.jl")
include("image.jl")
include("material.jl")
include("mesh.jl")
include("node.jl")
include("sampler.jl")
include("scene.jl")
include("skin.jl")
include("texture.jl")
include("object.jl")
include("show.jl")

function load(path::AbstractString)
    open(path) do f
        JSON3.read(read(f, String), Object)
    end
end

function save(path::AbstractString, obj::Object)
    open(path, write=true) do f
        write(f, JSON3.write(obj))
    end
end

end # module