text:
@doc raw"""
    AbstractMetric

Abstract type for the pseudo-Riemannian metric tensor ``g``, a family of smoothly
varying inner products on the tangent space. See [`inner`](@ref).

# Functor

    (metric::Metric)(M::Manifold)

Generate the `MetricManifold` that wraps the manifold `M` with given `metric`.
This works for both a variable containing the metric as well as a subtype `T<:AbstractMetric`,
where a zero parameter constructor `T()` is availabe.
"""
abstract type AbstractMetric end

# piping syntax for decoration
(metric::AbstractMetric)(M::AbstractManifold) = MetricManifold(M, metric)
(::Type{T})(M::AbstractManifold) where {T<:AbstractMetric} = MetricManifold(M, T())

"""
    MetricManifold{𝔽,M<:AbstractManifold{𝔽},G<:AbstractMetric} <: AbstractDecoratorManifold{𝔽}

Equip a [`AbstractManifold`](@ref) explicitly with a [`AbstractMetric`](@ref) `G`.

For a Metric AbstractManifold, by default, assumes, that you implement the linear form
from [`local_metric`](@ref) in order to evaluate the exponential map.

If the corresponding [`AbstractMetric`](@ref) `G` yields closed form formulae for e.g.
the exponential map and this is implemented directly (without solving the ode),
you can of course still implement that directly.

# Constructor

    MetricManifold(M, G)

Generate the [`AbstractManifold`](@ref) `M` as a manifold with the [`AbstractMetric`](@ref) `G`.
"""
struct MetricManifold{𝔽,M<:AbstractManifold{𝔽},G<:AbstractMetric} <:
       AbstractConnectionManifold{𝔽}
    manifold::M
    metric::G
end

@doc raw"""
    RiemannianMetric <: AbstractMetric

Abstract type for Riemannian metrics, a family of positive definite inner
products. The positive definite property means that for ``X  ∈ T_p \mathcal M``, the
inner product ``g(X, X) > 0`` whenever ``X`` is not the zero vector.
"""
abstract type RiemannianMetric <: AbstractMetric end

@doc raw"""
    christoffel_symbols_first(
        M::MetricManifold,
        p,
        B::AbstractBasis;
        backend::AbstractDiffBackend = diff_backend(),
    )

Compute the Christoffel symbols of the first kind in local coordinates of basis `B`.
The Christoffel symbols are (in Einstein summation convention)

```math
Γ_{ijk} = \frac{1}{2} \Bigl[g_{kj,i} + g_{ik,j} - g_{ij,k}\Bigr],
```

where ``g_{ij,k}=\frac{∂}{∂ p^k} g_{ij}`` is the coordinate
derivative of the local representation of the metric tensor. The dimensions of
the resulting multi-dimensional array are ordered ``(i,j,k)``.
"""
christoffel_symbols_first(::AbstractManifold, ::Any, B::AbstractBasis)
function christoffel_symbols_first(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    ∂g = local_metric_jacobian(M, p, B; backend=backend)
    n = size(∂g, 1)
    Γ = allocate(∂g, Size(n, n, n))
    @einsum Γ[i, j, k] = 1 / 2 * (∂g[k, j, i] + ∂g[i, k, j] - ∂g[i, j, k])
    return Γ
end
@decorator_transparent_signature christoffel_symbols_first(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

function christoffel_symbols_second(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    Ginv = inverse_local_metric(M, p, B)
    Γ₁ = christoffel_symbols_first(M, p, B; backend=backend)
    Γ₂ = allocate(Γ₁)
    @einsum Γ₂[l, i, j] = Ginv[k, l] * Γ₁[i, j, k]
    return Γ₂
end

"""
    connection(::MetricManifold)

Return the [`LeviCivitaConnection`](@ref) for a metric manifold.
"""
connection(::MetricManifold) = LeviCivitaConnection()

@doc raw"""
    det_local_metric(M::AbstractManifold, p, B::AbstractBasis)

Return the determinant of local matrix representation of the metric tensor ``g``, i.e. of the
matrix ``G(p)`` representing the metric in the tangent space at ``p`` with as a matrix.

See also [`local_metric`](@ref)
"""
det_local_metric(::AbstractManifold, p, ::AbstractBasis)
function det_local_metric(M::AbstractManifold, p, B::AbstractBasis)
    return det(local_metric(M, p, B))
end
@decorator_transparent_signature det_local_metric(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis,
)
"""
    einstein_tensor(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend())

Compute the Einstein tensor of the manifold `M` at the point `p`, see [https://en.wikipedia.org/wiki/Einstein_tensor](https://en.wikipedia.org/wiki/Einstein_tensor)
"""
einstein_tensor(::AbstractManifold, ::Any, ::AbstractBasis)
function einstein_tensor(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    Ric = ricci_tensor(M, p, B; backend=backend)
    g = local_metric(M, p, B)
    Ginv = inverse_local_metric(M, p, B)
    S = sum(Ginv .* Ric)
    G = Ric - g .* S / 2
    return G
end
@decorator_transparent_signature einstein_tensor(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

@doc raw"""
    flat(N::MetricManifold{M,G}, p, X::FVector{TangentSpaceType})

Compute the musical isomorphism to transform the tangent vector `X` from the
[`AbstractManifold`](@ref) `M` equipped with [`AbstractMetric`](@ref) `G` to a cotangent by
computing

````math
X^♭= G_p X,
````
where ``G_p`` is the local matrix representation of `G`, see [`local_metric`](@ref)
"""
flat(::MetricManifold, ::Any...)

@decorator_transparent_fallback function flat!(
    M::MetricManifold,
    ξ::CoTFVector,
    p,
    X::TFVector,
)
    g = local_metric(M, p, ξ.basis)
    copyto!(ξ.data, g * X.data)
    return ξ
end

@doc raw"""
    inverse_local_metric(M::AbstractcManifold, p, B::AbstractBasis)

Return the local matrix representation of the inverse metric (cometric) tensor, usually
written ``g^{ij}``.

See also [`local_metric`](@ref)
"""
inverse_local_metric(::AbstractManifold, ::Any, ::AbstractBasis)
function inverse_local_metric(M::AbstractManifold, p, B::AbstractBasis)
    return inv(local_metric(M, p, B))
end
@decorator_transparent_signature inverse_local_metric(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis,
)

default_decorator_dispatch(M::MetricManifold) = default_metric_dispatch(M)

"""
    is_default_metric(M, G)

Indicate whether the [`AbstractMetric`](@ref) `G` is the default metric for
the [`AbstractManifold`](@ref) `M`. This means that any occurence of
[`MetricManifold`](@ref)(M,G) where `typeof(is_default_metric(M,G)) = true`
falls back to just be called with `M` such that the [`AbstractManifold`](@ref) `M`
implicitly has this metric, for example if this was the first one implemented
or is the one most commonly assumed to be used.
"""
function is_default_metric(M::AbstractManifold, G::AbstractMetric)
    return _extract_val(default_metric_dispatch(M, G))
end

default_metric_dispatch(::AbstractManifold, ::AbstractMetric) = Val(false)
function default_metric_dispatch(M::MetricManifold)
    return default_metric_dispatch(base_manifold(M), metric(M))
end

"""
    is_default_metric(MM::MetricManifold)

Indicate whether the [`AbstractMetric`](@ref) `MM.G` is the default metric for
the [`AbstractManifold`](@ref) `MM.manifold,` within the [`MetricManifold`](@ref) `MM`.
This means that any occurence of
[`MetricManifold`](@ref)`(MM.manifold, MM.G)` where `is_default_metric(MM.manifold, MM.G)) = true`
falls back to just be called with `MM.manifold,` such that the [`AbstractManifold`](@ref) `MM.manifold`
implicitly has the metric `MM.G`, for example if this was the first one
implemented or is the one most commonly assumed to be used.
"""
function is_default_metric(M::MetricManifold)
    return _extract_val(default_metric_dispatch(M))
end

function Base.convert(::Type{MetricManifold{𝔽,MT,GT}}, M::MT) where {𝔽,MT,GT}
    return _convert_with_default(M, GT, default_metric_dispatch(M, GT()))
end

function _convert_with_default(
    M::MT,
    T::Type{<:AbstractMetric},
    ::Val{true},
) where {MT<:AbstractManifold}
    return MetricManifold(M, T())
end
function _convert_with_default(
    M::MT,
    T::Type{<:AbstractMetric},
    ::Val{false},
) where {MT<:AbstractManifold}
    return error(
        "Can not convert $(M) to a MetricManifold{$(MT),$(T)}, since $(T) is not the default metric.",
    )
end

@doc raw"""
    inner(N::MetricManifold{M,G}, p, X, Y)

Compute the inner product of `X` and `Y` from the tangent space at `p` on the
[`AbstractManifold`](@ref) `M` using the [`AbstractMetric`](@ref) `G`. If `G` is the default
metric (see [`is_default_metric`](@ref)) this is done using `inner(M, p, X, Y)`,
otherwise the [`local_metric`](@ref)`(M, p)` is employed as

````math
g_p(X, Y) = ⟨X, G_p Y⟩,
````
where ``G_p`` is the loal matrix representation of the [`AbstractMetric`](@ref) `G`.
"""
inner(::MetricManifold, ::Any, ::Any, ::Any)

@decorator_transparent_fallback :intransparent function inner(
    M::MetricManifold,
    p,
    X::TFVector,
    Y::TFVector,
)
    X.basis === Y.basis ||
        error("calculating inner product of vectors from different bases is not supported")
    return dot(X.data, local_metric(M, p, X.basis) * Y.data)
end

@doc raw"""
    local_metric(M::AbstractManifold, p, B::AbstractBasis)

Return the local matrix representation at the point `p` of the metric tensor ``g`` with
respect to the [`AbstractBasis`](@ref) `B` on the [`AbstractManifold`](@ref) `M`, usually written ``g_{ij}``.
The matrix has the property that ``g(X, Y)=X^\mathrm{T} [g_{ij}] Y = g_{ij} X^i Y^j``,
where the latter expression uses Einstein summation convention.
The metric tensor is such that the formula works for the given [`AbstractBasis`](@ref) `B`.
"""
local_metric(::AbstractManifold, ::Any, ::AbstractBasis)
@decorator_transparent_signature local_metric(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

@doc raw"""
    local_metric_jacobian(
        M::AbstractManifold,
        p,
        B::AbstractBasis;
        backend::AbstractDiffBackend = diff_backend(),
    )

Get partial derivatives of the local metric of `M` at `p` in basis `B` with respect to the
coordinates of `p`, ``\frac{∂}{∂ p^k} g_{ij} = g_{ij,k}``. The
dimensions of the resulting multi-dimensional array are ordered ``(i,j,k)``.
"""
local_metric_jacobian(::AbstractManifold, ::Any, B::AbstractBasis)
function local_metric_jacobian(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    n = size(p, 1)
    ∂g = reshape(_jacobian(q -> local_metric(M, q, B), p, backend), n, n, n)
    return ∂g
end
@decorator_transparent_signature local_metric_jacobian(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

@doc raw"""
    log(N::MetricManifold{M,G}, p, q)

Copute the logarithmic map on the [`AbstractManifold`](@ref) `M` equipped with the [`AbstractMetric`](@ref) `G`.

If the metric was declared the default metric using [`is_default_metric`](@ref), this method
falls back to `log(M,p,q)`. Otherwise, you have to provide an implementation for the non-default
[`AbstractMetric`](@ref) `G` metric within its [`MetricManifold`](@ref)`{M,G}`.
"""
log(::MetricManifold, ::Any...)

@doc raw"""
    log_local_metric_density(M::AbstractManifold, p, B::AbstractBasis)

Return the natural logarithm of the metric density ``ρ`` of `M` at `p`, which
is given by ``ρ = \log \sqrt{|\det [g_{ij}]|}`` for the metric tensor expressed in basis `B`.
"""
log_local_metric_density(::AbstractManifold, ::Any, ::AbstractBasis)
function log_local_metric_density(M::AbstractManifold, p, B::AbstractBasis)
    return log(abs(det_local_metric(M, p, B))) / 2
end
@decorator_transparent_signature log_local_metric_density(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis,
)

@doc raw"""
    metric(M::MetricManifold)

Get the metric ``g`` of the manifold `M`.
"""
metric(::MetricManifold)

function metric(M::MetricManifold)
    return M.metric
end

@doc raw"""
    ricci_curvature(M::AbstractManifold, p, B::AbstractBasis; backend::AbstractDiffBackend = diff_backend())

Compute the Ricci scalar curvature of the manifold `M` at the point `p` using basis `B`.
The curvature is computed as the trace of the Ricci curvature tensor with respect to
the metric, that is ``R=g^{ij}R_{ij}`` where ``R`` is the scalar Ricci curvature at `p`,
``g^{ij}`` is the inverse local metric (see [`inverse_local_metric`](@ref)) at `p` and
``R_{ij}`` is the Riccie curvature tensor, see [`ricci_tensor`](@ref). Both the tensor and
inverse local metric are expressed in local coordinates defined by `B`, and the formula
uses the Einstein summation convention.
"""
ricci_curvature(::AbstractManifold, ::Any, ::AbstractBasis)
function ricci_curvature(
    M::AbstractManifold,
    p,
    B::AbstractBasis;
    backend::AbstractDiffBackend=diff_backend(),
)
    Ginv = inverse_local_metric(M, p, B)
    Ric = ricci_tensor(M, p, B; backend=backend)
    S = sum(Ginv .* Ric)
    return S
end
@decorator_transparent_signature ricci_curvature(
    M::AbstractDecoratorManifold,
    p,
    B::AbstractBasis;
    kwargs...,
)

@doc raw"""
    sharp(N::MetricManifold{M,G}, p, ξ::FVector{CotangentSpaceType})

Compute the musical isomorphism to transform the cotangent vector `ξ` from the
[`AbstractManifold`](@ref) `M` equipped with [`AbstractMetric`](@ref) `G` to a tangent by
computing

````math
ξ^♯ = G_p^{-1} ξ,
````
where ``G_p`` is the local matrix representation of `G`, i.e. one employs
[`inverse_local_metric`](@ref) here to obtain ``G_p^{-1}``.
"""
sharp(::MetricManifold, ::Any, ::CoTFVector)

function sharp!(M::N, X::TFVector, p, ξ::CoTFVector) where {N<:MetricManifold}
    Ginv = inverse_local_metric(M, p, X.basis)
    copyto!(X.data, Ginv * ξ.data)
    return X
end

function Base.show(io::IO, M::MetricManifold)
    return print(io, "MetricManifold($(M.manifold), $(M.metric))")
end

#
# Introduce transparency
# (a) new functions & other parents
for f in [
    christoffel_symbols_first,
    det_local_metric,
    einstein_tensor,
    inverse_local_metric,
    local_metric,
    local_metric_jacobian,
    log_local_metric_density,
    ricci_curvature,
]
    eval(
        quote
            function decorator_transparent_dispatch(
                ::typeof($f),
                ::AbstractConnectionManifold,
                args...,
            )
                return Val(:parent)
            end
        end,
    )
end

function decorator_transparent_dispatch(
    ::typeof(christoffel_symbols_second),
    ::MetricManifold,
    args...,
)
    return Val(:parent)
end


text:
# common utilities

## convenient type alias
#
#  These types signficantly reduces the need of using
#  type parameters in functions (which are often just
#  for the purpose of restricting the arrays to real)
#
# These could be removed when the Base supports
# covariant type notation, i.e. AbstractVector{<:Real}
#

typealias RealArray{T<:Real,N} AbstractArray{T,N}
typealias RealVector{T<:Real} AbstractArray{T,1}
typealias RealMatrix{T<:Real} AbstractArray{T,2}

typealias IntegerArray{T<:Integer,N} AbstractArray{T,N}
typealias IntegerVector{T<:Integer} AbstractArray{T,1}
typealias IntegerMatrix{T<:Integer} AbstractArray{T,2}

@compat typealias RealFP Union{Float32, Float64}

## conversion from real to fp types

@compat fptype{T<:Union{Float32,Bool,Int8,UInt8,Int16,UInt16}}(::Type{T}) = Float32
@compat fptype{T<:Union{Float64,Int32,UInt32,Int64,UInt64,Int128,UInt128}}(::Type{T}) = Float64
fptype(::Type{Complex64}) = Complex64
fptype(::Type{Complex128}) = Complex128

ext: jl


text:
"""
    setSelector(method::Function, varArgs...; genj::GeneticJuliaStructure = GenJ,
                nSelected::Real = 1.0, needsComparison::Bool = true,
                samplingWithRep::Bool = true)

Sets the selector of parents for `genj`, receiving a function that selects a
single individual among the population and the number of individuals that are
going to be selected. This function constructs a structure of type `Selector`
and adds it to the main structure.

# Arguments
- `method::Function`: function that receives the individuals of the population as
    first argument and `"varArgs"`, which is a list of arguments, as own arguments,
    and returns the index of the selected individual as an `Integer`.
- `varArgs`: list of arguments for the selection method. See `"Examples"`.

# Keyword Arguments
- `genj::GeneticJuliaStructure = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.
- `nSelected::Real = 1.0`: if `Integer`, it represents the exact number
    of parents, whereas if `Float`, it represents a ratio of the population.
- `needsComparison::Bool = true`: tells wether the user wants the system to compare
    between the individuals in order to select the ones that will replace the old
    population.
- `samplingWithRep::Bool = true`: specifies if the selection is done with or without
    replacement.

If the user specifies that comparison between individuals is necessary, the
selection method must receive the `compareFunction` function as it is
shown in `"Examples"`.

# Examples
```jldoctest
julia> function SelectAnIndividual(population::Array{Individual})
    #···Code···
    return selectedIndex
end
SelectAnIndividual (generic function with 1 method)

julia> setSelector(SelectAnIndividual, nSelected = 1.0, needsComparison = false)
SelectionOperator(SelectAnIndividual, 1.0, false, true, Any[])
```

```jldoctest
# Fitness comparison was not set here
julia> function AlsoSelectAnIndividual(population::Array{Individual}, arg1, arg2)
    #···Code···
    return selectedIndex
end
AlsoSelectAnIndividual (generic function with 1 method)

julia> setSelector(AlsoSelectAnIndividual, arg1, arg2, nSelected = 2, needsComparison = false)
SelectionOperator(AlsoSelectAnIndividual, 1.0, false, true, Any[arg1, arg2])
```

```jldoctest
julia> function YetAnotherSelectAnIndividual(population::Array{Individual}, compareFitness::Function)
    #···Code···
    return selectedIndex
end
YetAnotherSelectAnIndividual (generic function with 1 method)

julia> setSelector(YetAnotherSelectAnIndividual, nSelected = 1.0, needsComparison = true)
SelectionOperator(YetAnotherSelectAnIndividual, 1.0, true, true, Any[])
```

```jldoctest
# Fitness comparison was set here
julia> function SelectAnIndividualWithoutRep(population::Array{Individual}, compareFitness::Function)
    #···Code···
    return selectedIndex
end
SelectAnIndividualWithoutRep (generic function with 1 method)

julia> setSelector(SelectAnIndividualWithoutRep, nSelected = 1.0, needsComparison = true, samplingWithRep = false)
SelectionOperator(SelectAnIndividualWithoutRep, 1.0, true, false, Any[])
```

See also: [`Selector`](@ref)
"""
function setSelector(method::Function, varArgs...; genj::GeneticJuliaStructure = GenJ,
                     nSelected::Real = 1.0,
                     needsComparison::Bool = true,
                     samplingWithRep::Bool = true)

    if typeof(nSelected) <: Integer
        if nSelected < 1
            error("Number of selected individuals must be greater than 0")
        end
    else
        if nSelected < 0.0
            error("Number of selected individuals must be a ratio greater than 0")
        end
    end

    varArgs = [i for i in varArgs]
    genj._selector = SelectionOperator(method, nSelected, needsComparison, samplingWithRep, varArgs)
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using NDTensors
using Test

@testset "EmptyStorage" begin
  T = Tensor(EmptyStorage(NDTensors.EmptyNumber), (2, 2))
  @test size(T) == (2, 2)
  @test eltype(T) == NDTensors.EmptyNumber
  @test T[1, 1] == NDTensors.EmptyNumber()
  @test T[1, 2] == NDTensors.EmptyNumber()
  # TODO: This should fail with an out of bounds error!
  #@test T[1, 3] == NDTensors.EmptyNumber()

  Tc = complex(T)
  @test size(Tc) == (2, 2)
  @test eltype(Tc) == Complex{NDTensors.EmptyNumber}
  @test Tc[1, 1] == Complex(NDTensors.EmptyNumber(), NDTensors.EmptyNumber())
  @test Tc[1, 2] == Complex(NDTensors.EmptyNumber(), NDTensors.EmptyNumber())
end


text:
export mult_gauss, mult_gauss!, mult_gauss_chn, mult_gauss_chn!

@inline function comb_mult_gauss(x, n)
    return x * n
end

mult_gauss!(X, σ=0.1, μ=1; clip=false) = apply_noise!(comb_mult_gauss, f_awg(σ, μ, complex_copy(σ), complex_copy(μ)), X, clip)
mult_gauss(X, σ=0.1, μ=1; clip=false) = mult_gauss!(copy(X), σ, μ, clip=clip)

"""
    mult_gauss(X; clip=false[, σ=0.1, μ=1])

Returns the array `X` with the array value multiplied with a gauss distribution (standard deviation `σ` and mean `μ`) .
`σ` and `μ` are optional arguments representing standard deviation and mean of gauss.
If keyword argument `clip` is provided the values are clipped to be in [0, 1].
If `X` is a RGB{Normed} or Gray{Normed} image, then the values will be automatically clipped and the keyword
`clip` is meaningless.

If `X<:Complex`, `μ` and `σ` are applied to the imaginary in the same way as for the real part.
If you want to have different behaviour for real and imaginary part, simply
choose `μ` or `σ` complex.
"""
mult_gauss

mult_gauss_chn(X, σ=0.1, μ=1; clip=false) = mult_gauss_chn!(copy(X), σ, μ, clip=clip)
mult_gauss_chn!(X, σ=0.1, μ=1; clip=false) = apply_noise_chn!(comb_mult_gauss, f_awg(σ, μ, complex_copy(σ), complex_copy(μ)), X, clip)

"""
    mult_gauss_chn(X; clip=false[, σ=0.1, μ=0.0])

Returns the RGB image `X` with the values of the pixel multiplied with a
gauss distribution (standard deviation `σ` and mean `μ`) pixelwise.
However, every channel of one pixel receives the same amount of noise.
The noise therefore acts roughly as intensity - but not color - changing noise.
If keyword argument `clip` is provided the values are clipped to be in [0, 1].
`σ` and `μ` are optional arguments representing standard deviation and mean of gauss.

If `X<:Complex`, `μ` and `σ` are applied to the imaginary in the same way as for the real part.
If you want to have different behaviour for real and imaginary part, simply
choose `μ` or `σ` complex.
"""
mult_gauss_chn

max_issues_repo_name: roflmaostc/Noise

max_forks_repo_name: roflmaostc/Noise


text:
struct Quaternion{T <: Real} <: FieldVector{4,T}
    s::T
    v1::T
    v2::T
    v3::T
end

# Constructors
Quaternion(s::Real,v1::Real,v2::Real,v3::Real) = Quaternion(promote(s, v1, v2, v3)...)
Quaternion(s::Real) = Quaternion(s, 0, 0, 0)
# Quaternion(v::Vector) = Quaternion(0, v[1], v[2], v[3])
Quaternion(v::SVector{3,T}) where T = Quaternion(0, v[1], v[2], v[3])
Quaternion(s::T,v::SVector{3,T}) where T = Quaternion(s, v[1], v[2], v[3])
Quaternion(R::Rotation) = Quaternion(Quat(R).w, Quat(R).x, Quat(R).y, Quat(R).z)
Quaternion{T}() where T = Quaternion{T}(1, 0, 0, 0)

# Basic quaternion operations
LinearAlgebra.real(q::Quaternion) = q[1]
LinearAlgebra.imag(q::Quaternion) = q[SUnitRange(2, 4)]

LinearAlgebra.conj(q::Quaternion) = Quaternion(q.s, -q.v1, -q.v2, -q.v3)
Base.abs(q::Quaternion) = sqrt(q.s * q.s + q.v1 * q.v1 + q.v2 * q.v2 + q.v3 * q.v3)
Base.abs2(q::Quaternion) = q.s * q.s + q.v1 * q.v1 + q.v2 * q.v2 + q.v3 * q.v3
Base.inv(q::Quaternion) = conj(q)

Base.:*(q1::Quaternion, q2::Quaternion) = Quaternion(  q1.s * q2.s - q1.v1 * q2.v1 - q1.v2 * q2.v2 - q1.v3 * q2.v3,
                                                       q1.s * q2.v1 + q1.v1 * q2.s + q1.v2 * q2.v3 - q1.v3 * q2.v2,
                                                       q1.s * q2.v2 - q1.v1 * q2.v3 + q1.v2 * q2.s + q1.v3 * q2.v1,
                                                       q1.s * q2.v3 + q1.v1 * q2.v2 - q1.v2 * q2.v1 + q1.v3 * q2.s)
Base.:*(q::Quaternion, x::Number) = Quaternion(q.s * x, q.v1 * x, q.v2 * x, q.v3 * x)
Base.:*(x::Number, q::Quaternion) = q * x

Base.:/(q1::Quaternion, q2::Quaternion) = q1 * inv(q2)
Base.:\(q1::Quaternion, q2::Quaternion) = inv(q1) * q2

# Base.:-(q::Quaternion) = Quaternion(-q.s,-q.v1,-q.v2,-q.v3)

angleaxis(q::Quaternion) = angle(q), axis(q)
angle(q::Quaternion) = 2 * atan(sqrt(q.v1^2 + q.v2^2 + q.v3^2), q.s)
function axis(q::Quaternion{T}) where T
    qv = SVector(q.v1, q.v2, q.v3)
    if norm(qv) == 0
        return SVector{3,T}(0, 0, 0)
    else
        return qv / norm(qv)
    end
end

qrotate(x::Quaternion,q::Quaternion) = q * x / q
vrotate(x::SVector,q::Quaternion) = imag(qrotate(Quaternion(x), q))

# Matrix equivalences
# 𝟙(::Type{T}) where T = Quaternion(one(T))
# 𝟙() = 𝟙(Float64)
Vmat(::Type{T}) where T = SMatrix{3,4,T,12}(0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1)
Vmat() = Vmat(Float64)
Vmat(q::SVector) = q[SUnitRange(2, 4)]
Vmat(A::SMatrix) = A[SUnitRange(2, 4),:]
Vmat(q::Quaternion) = imag(q)
Vᵀmat(::Type{T}) where T = SMatrix{4,3,T,12}(0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1)
Vᵀmat() = Vᵀmat(Float64)
Tmat(::Type{T}) where T = SMatrix{4,4,T,16}(1, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1, 0, 0, 0, 0, -1)
Tmat() = Tmat(Float64)

Lmat(q::Quaternion{T}) where T = SMatrix{4,4,T,16}(q.s, q.v1, q.v2, q.v3, -q.v1, q.s, q.v3, -q.v2, -q.v2, -q.v3, q.s, q.v1, -q.v3, q.v2, -q.v1, q.s)
Lᵀmat(q::Quaternion) = Lmat(q)'
Rmat(q::Quaternion{T})  where T = SMatrix{4,4,T,16}(q.s, q.v1, q.v2, q.v3, -q.v1, q.s, -q.v3, q.v2, -q.v2, q.v3, q.s, -q.v1, -q.v3, -q.v2, q.v1, q.s)
Rᵀmat(q::Quaternion) = Rmat(q)'

VLmat(q::Quaternion{T}) where T = SMatrix{3,4,T,12}(q.v1, q.v2, q.v3, q.s, q.v3, -q.v2, -q.v3, q.s, q.v1, q.v2, -q.v1, q.s)
VLᵀmat(q::Quaternion{T}) where T = SMatrix{3,4,T,12}(-q.v1, -q.v2, -q.v3, q.s, -q.v3, q.v2, q.v3, q.s, -q.v1, -q.v2, q.v1, q.s)
VRmat(q::Quaternion{T})  where T = SMatrix{3,4,T,12}(q.v1, q.v2, q.v3, q.s, -q.v3, q.v2, q.v3, q.s, -q.v1, -q.v2, q.v1, q.s)
VRᵀmat(q::Quaternion{T})  where T = SMatrix{3,4,T,12}(-q.v1, -q.v2, -q.v3, q.s, q.v3, -q.v2, -q.v3, q.s, q.v1, q.v2, -q.v1, q.s)

LVᵀmat(q::Quaternion{T}) where T = SMatrix{4,3,T,12}(-q.v1, q.s, q.v3, -q.v2, -q.v2, -q.v3, q.s, q.v1, -q.v3, q.v2, -q.v1, q.s)
LᵀVᵀmat(q::Quaternion{T}) where T = SMatrix{4,3,T,12}(q.v1, q.s, -q.v3, q.v2, q.v2, q.v3, q.s, -q.v1, q.v3, -q.v2, q.v1, q.s)
RVᵀmat(q::Quaternion{T}) where T = SMatrix{4,3,T,12}(-q.v1, q.s, -q.v3, q.v2, -q.v2, q.v3, q.s, -q.v1, -q.v3, -q.v2, q.v1, q.s)
RᵀVᵀmat(q::Quaternion{T}) where T = SMatrix{4,3,T,12}(q.v1, q.s, q.v3, -q.v2, q.v2, -q.v3, q.s, q.v1, q.v3, q.v2, -q.v1, q.s)

skewplusdiag(v::AbstractVector{T},w::T) where T = SMatrix{3,3,T,9}(w, v[3], -v[2], -v[3], w, v[1], v[2], -v[1], w)


text:
mutable struct FastSpikingNeuron <: SpecialDiffSys
    C::Float64
    G_Na::Float64
    G_K::Float64
    G_L::Float64
    E_Na::Float64
    E_K::Float64
    E_L::Float64
    I_bias::Float64
    noise::AbstractArray
end

name(::FastSpikingNeuron) = "Fast spiking neuron equation"
dimension(::FastSpikingNeuron) = 4

FastSpikingNeuron(; C=1.0, G_Na=112.0, G_K=224.0, G_L=0.1, E_Na=55.0, E_K=-97.0, E_L=-70.0, I_bias=6.0, noise=zeros(4)) = FastSpikingNeuron(C, G_Na, G_K, G_L, E_Na, E_K, E_L, I_bias, noise)

function vec_field(fsn::FastSpikingNeuron, x::AbstractArray)
    V, m, h, n = x
    C = fsn.C
    G_Na, G_K, G_L = fsn.G_Na, fsn.G_K, fsn.G_L
    E_Na, E_K, E_L = fsn.E_Na, fsn.E_K, fsn.E_L
    I_bias = fsn.I_bias
    _dV = G_Na*m^3*h*(E_Na-V) + G_K*n^2*(E_K-V) + G_L*(E_L-V) + I_bias
    dV = _dV/C
    αₘ = 40(V-75)/(1-exp((75-V)/13.5))
    βₘ = 1.2262exp(-V/42.248)
    dm = αₘ*(1-m) - βₘ*m
    αₕ = 0.0035exp(-V/24.186)
    βₕ = 0.017(-51.25-V)/(exp((-51.25-V)/5.2)-1)
    dh = αₕ*(1-h) - βₕ*h
    αₙ = (V-95)/(1-exp((95-V)/11.8))
    βₙ = 0.025exp(-V/22.222)
    dn = αₙ*(1-n) - βₙ*n
    return [dV, dm, dh, dn]
end

function jacobian(fsn::FastSpikingNeuron, x::AbstractArray)
    V, m, h, n = x
    C = fsn.C
    G_Na, G_K, G_L = fsn.G_Na, fsn.G_K, fsn.G_L
    E_Na, E_K = fsn.E_Na, fsn.E_K
    _dVdV = -G_Na*m^3*h - G_K*n^2 - G_L
    _dVdm = 3m^2*G_Na*h*(E_Na-V)
    _dVdh = G_Na*m^3*(E_Na-V)
    _dVdn = 2G_K*n*(E_K-V)
    dVs = [_dVdV _dVdm _dVdh _dVdn]/C
    αₘ = 40(V-75)/(1-exp((75-V)/13.5))
    βₘ = 1.2262exp(-V/42.248)
    dαₘ = begin
        y = (V-75)/13.5
        40*(1-exp(-y)-y*exp(-y))/(1-exp(-y))^2
    end
    dβₘ = -(1.2262/42.248)*exp(-V/42.248)
    dms = [dαₘ*(1-m)-dβₘ*m -αₘ-βₘ 0 0]
    αₕ = 0.0035exp(-V/24.186)
    βₕ = 0.017(-51.25-V)/(exp((-51.25-V)/5.2)-1)
    dαₕ = -(0.0035/24.186)*exp(-V/24.186)
    dβₕ = begin
        y = (V+51.25)/5.2
        0.017*(1-exp(-y)-y*exp(-y))/(1-exp(-y))^2
    end
    dhs = [dαₕ*(1-h)-dβₕ*h 0 -αₕ-βₕ 0]
    αₙ = (V-95)/(1-exp((95-V)/11.8))
    βₙ = 0.025exp(-V/22.222)
    dαₙ = begin
        y = (V-95)/11.8
        (1-exp(-y)-y*exp(-y))/(1-exp(-y))^2
    end
    dβₙ = -(0.025/22.222)*exp(-V/22.222)
    dns = [dαₙ*(1-n)-dβₙ*n 0 0 -αₙ-βₙ]
    return vcat(dVs, dms, dhs, dns)
end


text:
function Initi_model_struct(H::T,Hc::T,N::T,Wf::T,B::T,αf::T,β₁::T,β₂::T,
    Nₛ::T,rₘ::T,a_Jmax::T,b_Jmax::T,i::T,Kₓₗ₀::T) where {T<:Float64}
    cons = Constants()
    env = EnvironmentStruct()
    kinetic = PhotoKineticRates()
    photo = PhotoPar(kinetic,env.Tₐ)
    hydPar = HydraulicsPar(;Kₓₗ₀=Kₓₗ₀,i=i)
    treepar = TreePar(;αf=αf,β₁=β₁,β₂=β₂,rₘ=rₘ,Nₛ=Nₛ,a_Jmax=a_Jmax,b_Jmax=b_Jmax)

    Hs = H-Hc
    As = Wf/treepar.αf
    Ww = treepar.ρw*As*(treepar.β₁*H+treepar.β₂*Hs)
    treesize = TreeSize(Wf,Ww,H,Hs,As,B,N)

    model = CCPHStruct(cons,env,treepar,treesize,photo,hydPar)

    return model,kinetic
end

function Initi_model_struct(j::Integer,data::RoData,
    αf::T,β₁::T,β₂::T,Nₛ::T,rₘ::T,a_Jmax::T,b_Jmax::T,i::T,Kₓₗ₀::T) where {T<:Float64}
    data_ind = Find_data_ind(j)
    H = data.H[data_ind]
    N = data.N[data_ind]
    Wf = data.Wf[data_ind]
    B = data.B[data_ind]
    Hc = data.Hc
    model,kinetic = Initi_model_struct(H,Hc,N,Wf,B,αf,β₁,β₂,Nₛ,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

    return model,kinetic
end

function Run_RO_C_F_CCPH(par::Array{Float64,1},weatherts_F::WeatherTS,
    weatherts_C::WeatherTS,data_F::RoData,data_C::RoData;sim_steps::Integer=83)

    αf,β₁,β₂,b_Jmax = 460.0,1.27,-0.27,0.0
    Nₛ,rₘ,a_Jmax,Kₓₗ₀,Nₛ_C,i = par[1:6]

    #Fertilized stand
    GPP_model_F = zeros(sim_steps+1)
    EC_model_F = zeros(sim_steps+1)
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_F,αf,β₁,β₂,Nₛ,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

        GPP_model_F[j],EC_model_F[j], modeloutput, gₛ_opt, Nₘ_f_opt = CalcModelOutput!(j,model,weatherts_F,kinetic)
    end

    #Control stand
    GPP_model_C = zeros(sim_steps+1)
    EC_model_C = zeros(sim_steps+1)
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_C,αf,β₁,β₂,Nₛ_C,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

        GPP_model_C[j],EC_model_C[j], modeloutput, gₛ_opt, Nₘ_f_opt = CalcModelOutput!(j,model,weatherts_C,kinetic)
    end

    return GPP_model_F,GPP_model_C,EC_model_F,EC_model_C
end

function Get_Data_RO_C_F_CCPH(par::Array{Float64,1},weatherts_F::WeatherTS,
    weatherts_C::WeatherTS,data_F::RoData,data_C::RoData;sim_steps::Integer=83)

    αf,β₁,β₂,b_Jmax = 460.0,1.27,-0.27,0.0
    Nₛ,rₘ,a_Jmax,Kₓₗ₀,Nₛ_C,i = par[1:6]

    #Fertilized stand
    GPP_model_F = zeros(sim_steps+1)
    EC_model_F = zeros(sim_steps+1)
    gₛ_model_F = zeros(sim_steps+1)
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_F,αf,β₁,β₂,Nₛ,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

        GPP_model_F[j],EC_model_F[j], modeloutput, gₛ_model_F[j], Nₘ_f_opt = CalcModelOutput!(j,model,weatherts_F,kinetic)
    end

    #Control stand
    GPP_model_C = zeros(sim_steps+1)
    EC_model_C = zeros(sim_steps+1)
    gₛ_model_C = zeros(sim_steps+1)
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_C,αf,β₁,β₂,Nₛ_C,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

        GPP_model_C[j],EC_model_C[j], modeloutput, gₛ_model_C[j], Nₘ_f_opt = CalcModelOutput!(j,model,weatherts_C,kinetic)
    end

    return GPP_model_F,GPP_model_C,EC_model_F,EC_model_C,gₛ_model_F,gₛ_model_C
end

function calc_error_term(ydata::T,ymodel::T) where {T<:Float64}
    e = (ydata-ymodel)^2

    return e
end

function calc_logP_term(ydata::T,ymodel::T,a::T,b::T) where {T<:Float64}
    e = abs(ydata-ymodel)
    σ = a+b*ymodel
    σ>0.0||error("Negative variance")
    return e/σ+log(σ)
end

function obj_logP(x::Array{T,1},GPP_model_F::Array{T,1},Ec_model_F::Array{T,1},data_F::RoData,
    GPP_model_C::Array{T,1},Ec_model_C::Array{T,1},data_C::RoData) where {T<:Float64}
    a_GPP,b_GPP,a_Ec,b_Ec = x

    logP = 0.0
    for i = 1:length(data_F.GPP)
        logP += calc_logP_term(data_F.GPP[i],GPP_model_F[i],a_GPP,b_GPP)
        logP += calc_logP_term(data_C.GPP[i],GPP_model_C[i],a_GPP,b_GPP)

        logP += calc_logP_term(data_F.E_C[i],Ec_model_F[i],a_Ec,b_Ec)
        logP += calc_logP_term(data_C.E_C[i],Ec_model_C[i],a_Ec,b_Ec)
    end

    return logP
end

function calc_error(
    GPP_model_F::Array{T,1},Ec_model_F::Array{T,1},data_F::RoData,
    GPP_model_C::Array{T,1},Ec_model_C::Array{T,1},data_C::RoData) where {T<:Float64}

    e_tot = 0.0

    #=
    for i = 1:length(data_F.GPP)
        e_tot += calc_error_term(data_F.GPP[i],GPP_model_F[i])
        e_tot += calc_error_term(data_C.GPP[i],GPP_model_C[i])

        e_tot += calc_error_term(data_F.E_C[i],Ec_model_F[i])
        e_tot += calc_error_term(data_C.E_C[i],Ec_model_C[i])
    end
    =#
    e_tot += calcR²(data_F.GPP,GPP_model_F)
    e_tot += calcR²(data_C.GPP,GPP_model_C)

    e_tot += calcR²(data_F.E_C,Ec_model_F)
    e_tot += calcR²(data_C.E_C,Ec_model_C)

    return -e_tot
end

function obj_RO_C_F_CCPH(par::Array{Float64,1},RO_data::RO_raw_data)
    try
        X0,τ,Smax,τ_C = par[7:10]
        weatherts_F,GPP_data_F = create_weather_struct_RO(RO_data;stand_type="Fertilized",X0=X0,τ=τ,Smax=Smax)
        weatherts_C,GPP_data_C = create_weather_struct_RO(RO_data;stand_type="Control",X0=X0,τ=τ_C,Smax=Smax)

        data_F,data_C = Create_RoData_C_F(GPP_data_F,GPP_data_C,weatherts_F,weatherts_C)

        GPP_model_F,GPP_model_C,EC_model_F,EC_model_C = Run_RO_C_F_CCPH(par,weatherts_F,weatherts_C,data_F,data_C)

        e_tot = calc_error(
        GPP_model_F,EC_model_F,data_F,
        GPP_model_C,EC_model_C,data_C)

        return e_tot
    catch err
        println("Parameters Error: ", err)
        return Inf
    end
end

function run_opt()
    #file_name = "opt_no_error_RO_C_GPP_Ec_20220223_2"
    file_name = "opt_no_error_RO_C_GPP_Ec_20220224"

    RO_data = Load_RO_data()

    #Nₛ,rₘ,a_Jmax,Kₓₗ₀,Nₛ_C,i,X0,τ,Smax,τ_C,a_GPP,b_GPP,a_Ec,b_Ec
    ranges = [(0.001,0.1),(10.0,40.0),(0.01,1.0),(0.001,0.1),(0.001,0.1),
    (0.1,6.0),(-18.0,-0.9),(1.0,15.0),(10.0,25.0),(1.0,15.0)]

    res = BlackBoxOptim.bboptimize(x::Array{Float64,1}->obj_RO_C_F_CCPH(x::Array{Float64,1},RO_data)
    ; SearchRange = ranges,PopulationSize = 200,MaxSteps=30000)

    @show xopt = BlackBoxOptim.best_candidate(res)

    save("./output/"*file_name*".jld","xopt",xopt)
end
function run_opt_par()
    file_name = "opt_no_error_RO_C_GPP_Ec_20220224"

    par = load("./output/"*file_name*".jld","xopt")

    #=
    [0.012854764847976128, 14.865444332753684, 0.26557235210247054,
    0.03784509575809752, 0.013224669959225405, 0.36827028167982595,
    -14.259069117859587, 8.754639146359056, 14.402015724819176, 1.2640782289025818]
    =#

    RO_data = Load_RO_data()

    X0,τ,Smax,τ_C = par[7:10]
    weatherts_F,GPP_data_F = create_weather_struct_RO(RO_data;stand_type="Fertilized",X0=X0,τ=τ,Smax=Smax)
    weatherts_C,GPP_data_C = create_weather_struct_RO(RO_data;stand_type="Control",X0=X0,τ=τ_C,Smax=Smax)

    data_F,data_C = Create_RoData_C_F(GPP_data_F,GPP_data_C,weatherts_F,weatherts_C)

    GPP_model_F,GPP_model_C,EC_model_F,EC_model_C,gₛ_model_F,gₛ_model_C = Get_Data_RO_C_F_CCPH(par,weatherts_F,weatherts_C,data_F,data_C)

    @show GPP_R²_F =  calcR²(data_F.GPP,GPP_model_F)
    @show GPP_cor_F = cor(data_F.GPP,GPP_model_F)
    @show gₛ_R²_F =  calcR²(data_F.E_C,EC_model_F)
    @show gₛ_cor_F = cor(data_F.E_C,EC_model_F)

    @show GPP_R²_C =  calcR²(data_C.GPP,GPP_model_C)
    @show GPP_cor_C = cor(data_C.GPP,GPP_model_C)
    @show gₛ_R²_C =  calcR²(data_C.E_C,EC_model_C)
    @show gₛ_cor_C = cor(data_C.E_C,EC_model_C)

    println("Fertilized")
    clac_GPP_R²_annual(data_C.GPP,GPP_model_F,weatherts_F.date)
    println("Control")
    clac_GPP_R²_annual(data_C.GPP,GPP_model_C,weatherts_C.date)

    plot(weatherts_F.date,data_F.GPP,label="Data",ylabel="GPP")
    pl1 = plot!(weatherts_F.date,GPP_model_F,label="Model")
    plot(weatherts_C.date,data_C.GPP,label="Data",ylabel="GPP")
    pl2 = plot!(weatherts_C.date,GPP_model_C,label="Model")

    plot(weatherts_F.date,data_F.E_C,label="Data",ylabel="E_C")
    pl3 = plot!(weatherts_F.date,EC_model_F,label="Model")
    plot(weatherts_C.date,data_C.E_C,label="Data",ylabel="E_C")
    pl4 = plot!(weatherts_C.date,EC_model_C,label="Model")

    plot(pl1,pl3,pl2,pl4,layout=(2,2),legends=false)
    savefig("./plots/"*file_name*"_result.svg")

    plot(weatherts_F.date,data_F.gₛ,label="Data",ylabel="gₛ")
    pl1 = plot!(weatherts_F.date,gₛ_model_F,label="Model",ylabel="gₛ")
    plot(weatherts_C.date,data_C.gₛ,label="Data",ylabel="gₛ")
    pl2 = plot!(weatherts_C.date,gₛ_model_C,label="Model",ylabel="gₛ")

    plot(pl1,pl2,layout=(1,2),legends=false)
    savefig("./plots/"*file_name*"_result_gs.svg")

    #Test optimizaiton
    inner_optimizer = Optim.BFGS(linesearch=Optim.LineSearches.BackTracking())

    lower = [0.0001,0.0001,0.0001,0.0001]
    upper = [5.0,3.0,5.0,3.0]
    x0 = [0.237,0.0845,0.0688,0.146]

    opt_trait = Optim.optimize(x->obj_logP(x,GPP_model_F,EC_model_F,data_F,
    GPP_model_C,EC_model_C,data_C),
    lower,upper,x0,Optim.Fminbox(inner_optimizer))

    opt_trait.minimizer
end

#run_opt()
run_opt_par()


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

## Binary arithmetic operators ##

for f in (:+, :-)
    @eval function ($f)(A::AbstractArray, B::AbstractArray)
        promote_shape(A, B) # check size compatibility
        broadcast_preserving_zero_d($f, A, B)
    end
end

function +(A::Array, Bs::Array...)
    for B in Bs
        promote_shape(A, B) # check size compatibility
    end
    broadcast_preserving_zero_d(+, A, Bs...)
end

for f in (:/, :\, :*)
    if f !== :/
        @eval ($f)(A::Number, B::AbstractArray) = broadcast_preserving_zero_d($f, A, B)
    end
    if f !== :\
        @eval ($f)(A::AbstractArray, B::Number) = broadcast_preserving_zero_d($f, A, B)
    end
end

## data movement ##

"""
    reverse(A; dims=:)

Reverse `A` along dimension `dims`, which can be an integer (a
single dimension), a tuple of integers (a tuple of dimensions)
or `:` (reverse along all the dimensions, the default).  See
also [`reverse!`](@ref) for in-place reversal.

# Examples
```jldoctest
julia> b = Int64[1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> reverse(b, dims=2)
2×2 Matrix{Int64}:
 2  1
 4  3

julia> reverse(b)
2×2 Matrix{Int64}:
 4  3
 2  1
```

!!! compat "Julia 1.6"
    Prior to Julia 1.6, only single-integer `dims` are supported in `reverse`.
"""
reverse(A::AbstractArray; dims=:) = _reverse(A, dims)
_reverse(A, dims) = reverse!(copymutable(A); dims)

"""
    reverse!(A; dims=:)

Like [`reverse`](@ref), but operates in-place in `A`.

!!! compat "Julia 1.6"
    Multidimensional `reverse!` requires Julia 1.6.
"""
reverse!(A::AbstractArray; dims=:) = _reverse!(A, dims)
_reverse!(A::AbstractArray{<:Any,N}, ::Colon) where {N} = _reverse!(A, ntuple(identity, Val{N}()))
_reverse!(A, dim::Integer) = _reverse!(A, (Int(dim),))
_reverse!(A, dims::NTuple{M,Integer}) where {M} = _reverse!(A, Int.(dims))
function _reverse!(A::AbstractArray{<:Any,N}, dims::NTuple{M,Int}) where {N,M}
    dimrev = ntuple(k -> k in dims, Val{N}()) # boolean tuple indicating reversed dims

    if N < M || M != sum(dimrev)
        throw(ArgumentError("invalid dimensions $dims in reverse!"))
    end
    M == 0 && return A # nothing to reverse

    # swapping loop only needs to traverse ≈half of the array
    halfsz = ntuple(k -> k == dims[1] ? size(A,k) ÷ 2 : size(A,k), Val{N}())

    last1 = ntuple(k -> lastindex(A,k)+firstindex(A,k), Val{N}()) # offset for reversed index
    for i in CartesianIndices(ntuple(k -> firstindex(A,k):firstindex(A,k)-1+@inbounds(halfsz[k]), Val{N}()))
        iₜ = Tuple(i)
        iᵣ = CartesianIndex(ifelse.(dimrev, last1 .- iₜ, iₜ))
        @inbounds A[iᵣ], A[i] = A[i], A[iᵣ]
    end
    if M > 1 && isodd(size(A, dims[1]))
        # middle slice for odd dimensions must be recursively flipped
        mid = firstindex(A, dims[1]) + (size(A, dims[1]) ÷ 2)
        midslice = CartesianIndices(ntuple(k -> k == dims[1] ? (mid:mid) : (firstindex(A,k):lastindex(A,k)), Val{N}()))
        _reverse!(view(A, midslice), dims[2:end])
    end
    return A
end
# fix ambiguity with array.jl:
_reverse!(A::AbstractVector, dim::Tuple{Int}) = _reverse!(A, first(dim))

"""
    rotl90(A)

Rotate matrix `A` left 90 degrees.

# Examples
```jldoctest
julia> a = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> rotl90(a)
2×2 Matrix{Int64}:
 2  4
 1  3
```
"""
function rotl90(A::AbstractMatrix)
    ind1, ind2 = axes(A)
    B = similar(A, (ind2,ind1))
    n = first(ind2)+last(ind2)
    for i=axes(A,1), j=ind2
        B[n-j,i] = A[i,j]
    end
    return B
end

"""
    rotr90(A)

Rotate matrix `A` right 90 degrees.

# Examples
```jldoctest
julia> a = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> rotr90(a)
2×2 Matrix{Int64}:
 3  1
 4  2
```
"""
function rotr90(A::AbstractMatrix)
    ind1, ind2 = axes(A)
    B = similar(A, (ind2,ind1))
    m = first(ind1)+last(ind1)
    for i=ind1, j=axes(A,2)
        B[j,m-i] = A[i,j]
    end
    return B
end
"""
    rot180(A)

Rotate matrix `A` 180 degrees.

# Examples
```jldoctest
julia> a = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> rot180(a)
2×2 Matrix{Int64}:
 4  3
 2  1
```
"""
function rot180(A::AbstractMatrix)
    B = similar(A)
    ind1, ind2 = axes(A,1), axes(A,2)
    m, n = first(ind1)+last(ind1), first(ind2)+last(ind2)
    for j=ind2, i=ind1
        B[m-i,n-j] = A[i,j]
    end
    return B
end
"""
    rotl90(A, k)

Left-rotate matrix `A` 90 degrees counterclockwise an integer `k` number of times.
If `k` is a multiple of four (including zero), this is equivalent to a `copy`.

# Examples
```jldoctest
julia> a = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> rotl90(a,1)
2×2 Matrix{Int64}:
 2  4
 1  3

julia> rotl90(a,2)
2×2 Matrix{Int64}:
 4  3
 2  1

julia> rotl90(a,3)
2×2 Matrix{Int64}:
 3  1
 4  2

julia> rotl90(a,4)
2×2 Matrix{Int64}:
 1  2
 3  4
```
"""
function rotl90(A::AbstractMatrix, k::Integer)
    k = mod(k, 4)
    k == 1 ? rotl90(A) :
    k == 2 ? rot180(A) :
    k == 3 ? rotr90(A) : copy(A)
end
"""
    rotr90(A, k)

Right-rotate matrix `A` 90 degrees clockwise an integer `k` number of times.
If `k` is a multiple of four (including zero), this is equivalent to a `copy`.

# Examples
```jldoctest
julia> a = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> rotr90(a,1)
2×2 Matrix{Int64}:
 3  1
 4  2

julia> rotr90(a,2)
2×2 Matrix{Int64}:
 4  3
 2  1

julia> rotr90(a,3)
2×2 Matrix{Int64}:
 2  4
 1  3

julia> rotr90(a,4)
2×2 Matrix{Int64}:
 1  2
 3  4
```
"""
rotr90(A::AbstractMatrix, k::Integer) = rotl90(A,-k)
"""
    rot180(A, k)

Rotate matrix `A` 180 degrees an integer `k` number of times.
If `k` is even, this is equivalent to a `copy`.

# Examples
```jldoctest
julia> a = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> rot180(a,1)
2×2 Matrix{Int64}:
 4  3
 2  1

julia> rot180(a,2)
2×2 Matrix{Int64}:
 1  2
 3  4
```
"""
rot180(A::AbstractMatrix, k::Integer) = mod(k, 2) == 1 ? rot180(A) : copy(A)

max_stars_repo_name: wentasah/julia

max_issues_repo_name: wentasah/julia

max_forks_repo_name: wentasah/julia


text:
module EMLHeatExchangerSimplified_cost
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...heat_exchangers.EMLHeatex
	import EMLtypes.length
	include("HeatExchangerSimplified_cost/HeatExchanger_LMTD_cost.jl")
	include("HeatExchangerSimplified_cost/Shell_and_tubes_LMTD_cost.jl")
	include("HeatExchangerSimplified_cost/HeatExchanger_NTU_cost.jl")
	include("HeatExchangerSimplified_cost/Shell_and_tubes_NTU_cost.jl")
end


text:
module TestConsistency
    verbose = false

    using Base.Test

    import IntertemporalChoiceHeuristics:
        load,
        fit,
        models,
        simulate,
        make_closures,
        stderrs

    import Optim: minimizer

    import ForwardDiff: hessian

    import Distributions: ccdf, Chisq

    inputs, weights = load()
    n_trials = 1_000

    for model in models
        res = fit(model, inputs, weights)
        parameters = minimizer(res)
        m = model(parameters)

        for sim in 1:10
            sim_inputs, sim_weights = simulate(m, inputs, n_trials)

            sim_res = fit(model, sim_inputs, sim_weights)
            sim_parameters = minimizer(sim_res)
            sim_m = model(sim_parameters)

            p_value = NaN
            try
                se = stderrs(model, sim_inputs, sim_weights)
                z = (sim_parameters .- parameters) ./ se
                χ = sum(z.^2)
                p_value = ccdf(Chisq(length(parameters)), χ)
            end

            if verbose
                @printf(
                    "%s\t%.2f\t%.4f\n",
                    model,
                    norm(
                        100 * (parameters .- sim_parameters) ./ parameters,
                        Inf,
                    ),
                    p_value,
                )
            end
        end
    end
end


text:
module INOAssistant

using LibSerialPort
using ArgParse

include("read_port.jl")

end


text:
module CeterisParibus

    include("f().jl")

    include("g().jl")

    include("avg().jl")

    include("Wavg.jl")

    export f, g, avg, Wavg
end

#using Pkg
#Pkg.activate("CeterisParibus")


text:
export BERL_NEAT, process

function BERL_NEAT(cfg::Dict, fitness::Function)::Cambrian.Evolution
    evaluate!::Function = e::Evolution -> NEAT_evaluate!(e, fitness)
    selection::Function = i::Array{NEATIndiv} -> i[1]

    if cfg["selection_type"] == "tournament"
        selection =
            i::Array{NEATIndiv} -> NEAT_tournament(i, cfg["tournament_size"])
    elseif cfg["selection_type"] == "random_top"
        selection =
            i::Array{NEATIndiv} -> NEAT_random_top(i, cfg["survival_threshold"])
    else
        throw(ArgumentError("Wrong selection type: " + cfg["selection_type"]))
    end

    populate!::Function = e::Evolution -> NEAT_populate!(e, selection)
    Evolution(NEATIndividual, cfg; id=cfg["id"], evaluate = evaluate!, populate = populate!)
end

function process(i::NEATIndiv, X)
    NeuroEvolution.process(i, X)
end

algorithms["NEAT"]=BERL_NEAT


text:
Tables.istable(::Type{<:File}) = true
Tables.schema(f::File) = Tables.Schema(f.names, f.types)

Tables.columnaccess(::Type{File{t, c, I, P, KW}}) where {t, c, I, P, KW} = c

struct Columns
    names
    columns
end
Base.propertynames(c::Columns) = getfield(c, 1)
function Base.getproperty(c::Columns, nm::Symbol)
    getfield(c, 2)[findfirst(x->x==nm, getfield(c, 1))]
end

function Tables.columns(@nospecialize(f::File{transpose, true})) where {transpose}
    len = length(f)
    types = f.types
    columns = Any[Tables.allocatecolumn(T, len) for T in types]
    for (i, row) in enumerate(f)
        for col = 1:length(types)
            columns[col][i] = getproperty(row, types[col], col, :_)
        end
    end
    return Columns(f.names, columns)
end

# row interfaces
abstract type Row end

struct UntypedRow{F} <: Row
    file::F
    row::Int
end
Base.propertynames(r::UntypedRow) = getfield(r, 1).names

struct RowIterator{names, types, F}
    file::F
end
Tables.schema(r::RowIterator) = Tables.Schema(r.file.names, r.file.types)

struct TypedRow{F} <: Row
    rowiterator::F
    row::Int
end
Base.propertynames(r::TypedRow{F}) where {F <: RowIterator{names}} where {names} = names

function Base.show(io::IO, r::Row)
    println(io, "CSV.Row($(getfield(r, 2))) of:")
    show(io, getfield(r, 1))
end

getranspose(r::UntypedRow{F}) where {F <: File{transpose}} where {transpose} = transpose
getranspose(r::RowIterator{N, T, F}) where {N, T, F <: File{transpose}} where {transpose} = transpose
getranspose(r::TypedRow) = getranspose(getfield(r, 1))

# File iterates UntypedRows
Base.eltype(f::F) where {F <: File} = UntypedRow{F}
Base.length(f::File{transpose}) where {transpose} = transpose ? f.lastparsedcol[] : length(f.positions)
Base.size(f::File) = (length(f), length(f.names))

@inline function Base.iterate(f::File{transpose}, st=1) where {transpose}
    st > length(f) && return nothing
    # println("row=$st")
    if transpose
        st === 1 && (f.positions .= f.originalpositions)
    else
        @inbounds Parsers.fastseek!(f.io, f.positions[st])
        f.currentrow[] = st
        f.lastparsedcol[] = 0
        f.lastparsedcode[] = Parsers.SUCCESS
    end
    return UntypedRow(f, st), st + 1
end

# Tables.rows(f) -> RowIterator
Tables.rowaccess(::Type{<:File}) = true
Tables.rows(f::F) where {F <: File} = RowIterator{Tuple(f.names), Tuple{f.types...}, F}(f)

Base.eltype(f::R) where {R <: RowIterator} = TypedRow{R}
Base.length(f::RowIterator) = getranspose(f) ? f.file.lastparsedcol[] : length(f.file.positions)
Base.size(f::RowIterator) = (length(f.file), length(f.file.names))

# RowIterator iterates TypedRows
@inline function Base.iterate(r::RowIterator, st=1)
    f = r.file
    st > length(f) && return nothing
    transpose = getranspose(r)
    # println("row=$st")
    if transpose
        st === 1 && (f.positions .= f.originalpositions)
    else
        @inbounds Parsers.fastseek!(f.io, f.positions[st])
        f.currentrow[] = st
        f.lastparsedcol[] = 0
        f.lastparsedcode[] = Parsers.SUCCESS
    end
    return TypedRow(r, st), st + 1
end

parsingtype(::Type{Missing}) = Missing
parsingtype(::Type{Union{Missing, T}}) where {T} = T
parsingtype(::Type{T}) where {T} = T

@inline function parsefield(f, ::Type{CatStr}, row, col, strict)
    r = Parsers.parse(f.parsinglayers, f.io, Tuple{Ptr{UInt8}, Int})
    f.lastparsedcode[] = r.code
    if r.result isa Missing
        return missing
    else
        @inbounds pool = f.pools[col]
        val = r.result::Tuple{Ptr{UInt8}, Int}
        i = get(pool, val, nothing)
        if i === nothing
            i = get!(pool, unsafe_string(val[1], val[2]))
            issorted(levels(pool)) || levels!(pool, sort(levels(pool)))
        end
        return CatStr(i, pool)
    end
end

@inline function parsefield(f, T, row, col, strict)
    r = Parsers.parse(f.parsinglayers, f.io, T; f.kwargs...)
    f.lastparsedcode[] = r.code
    if !Parsers.ok(r.code)
        strict ? throw(Error(Parsers.Error(f.io, r), row, col)) :
            println("warning: failed parsing $T on row=$row, col=$col, error=$(Parsers.codes(r.code))")
    end
    return r.result
end

@noinline function skipcells(f, n)
    r = Parsers.Result(Tuple{Ptr{UInt8}, Int})
    for _ = 1:n
        r.code = Parsers.SUCCESS
        Parsers.parse!(f.parsinglayers, f.io, r)
        if newline(r.code)
            f.lastparsedcode[] = r.code
            return false
        end
    end
    f.lastparsedcode[] = r.code
    return true
end

@inline function Base.getproperty(csvrow::UntypedRow, name::Symbol)
    f = getfield(csvrow, 1)
    i = findfirst(x->x===name, f.names)
    return getproperty(csvrow, f.types[i], i, name)
end
@inline Base.getproperty(csvrow::TypedRow{F}, name::Symbol) where {F <: RowIterator{names, types}} where {names, types} =
    getproperty(csvrow, Tables.columntype(names, types, name), Tables.columnindex(names, name), name)

getfile(f::File) = f
getfile(r::RowIterator) = r.file

function Base.getproperty(csvrow::Row, ::Type{T}, col::Int, name::Symbol) where {T}
    col === 0 && return missing
    transpose = getranspose(csvrow)
    f = getfile(getfield(csvrow, 1))
    row = getfield(csvrow, 2)
    if transpose
        @inbounds Parsers.fastseek!(f.io, f.positions[col])
    else
        if f.currentrow[] != row
            @inbounds Parsers.fastseek!(f.io, f.positions[row])
            f.lastparsedcol[] = 0
        end
        lastparsed = f.lastparsedcol[]
        if col === lastparsed + 1
            if newline(f.lastparsedcode[])
                f.lastparsedcol[] = col
                return missing
            end
        elseif col > lastparsed + 1
            # skipping cells
            if newline(f.lastparsedcode[]) || !skipcells(f, col - (lastparsed + 1))
                f.lastparsedcol[] = col
                return missing
            end
        else
            @inbounds Parsers.fastseek!(f.io, f.positions[row])
            # randomly seeking within row
            if !skipcells(f, col - 1)
                f.lastparsedcol[] = col
                return missing
            end
        end
    end
    r = parsefield(f, parsingtype(T), row, col, f.strict)
    if transpose
        @inbounds f.positions[col] = position(f.io)
    else
        f.lastparsedcol[] = col
    end
    return r
end


text:
# "Abstract type that represents a time-domain, either `Discrete` or `Continuous`"
# abstract type TimeDomain end
# "Continuous time domain"
# struct Continuous <: TimeDomain end
# "Discrete (sampled) time domain"
# struct Discrete <: TimeDomain end

import ControlSystems: Continuous, Discrete, TimeEvolution
ControlSystems.Discrete() = Discrete(1)
Base.Broadcast.broadcastable(p::TimeEvolution) = Ref(p)

abstract type AbstractAssignmentMethod end

"""
    SortAssignement{F} <: AbstractAssignmentMethod

Contains a single Function field that determines what to sort roots by.
"""
struct SortAssignement{F} <: AbstractAssignmentMethod
    by::F
end
(a::SortAssignement)(e1,e2) = sortperm(e1, by=a.by), sortperm(e2, by=a.by)

"""
    HungarianAssignement <: AbstractAssignmentMethod

Sort roots using Hungarian method
"""
struct HungarianAssignement <: AbstractAssignmentMethod
end

(a::HungarianAssignement)(e1,e2) = 1:length(e1), hungarian(distmat_euclidean(e1,e2))[1]

abstract type AbstractRoots{T} <: AbstractVector{T} end

"""
    DiscreteRoots{T, V <: AbstractVector{T}} <: AbstractRoots{T}

Represent roots in discrete time
"""
struct DiscreteRoots{T, V <: AbstractVector{T}} <: AbstractRoots{T}
    r::V
    DiscreteRoots(r) = new{eltype(r), typeof(r)}(anglesort(reflectd.(r)))
end
StaticArrays.similar_type(r::DiscreteRoots{T,V}) where {T,V} = DiscreteRoots

"""
    ContinuousRoots{T, V <: AbstractVector{T}} <: AbstractRoots{T}

Represents roots in continuous time
"""
struct ContinuousRoots{T, V <: AbstractVector{T}} <: AbstractRoots{T}
    r::V
    ContinuousRoots(r) = new{eltype(r), typeof(r)}(eigsort(reflectc.(r)))
end
StaticArrays.similar_type(r::ContinuousRoots{T,V}) where {T,V} = ContinuousRoots

"""
    DiscreteRoots(r::ContinuousRoots) = begin

Represents roots of a polynomial in discrete time
"""
DiscreteRoots(r::ContinuousRoots) = domain_transform(Discrete(), r)
DiscreteRoots(r::DiscreteRoots) = r
"""
    ContinuousRoots(r::DiscreteRoots) = begin

Represents roots of a polynomial in continuous time
"""
ContinuousRoots(r::DiscreteRoots) = domain_transform(Continuous(), r)
ContinuousRoots(r::ContinuousRoots) = r

Base.convert(::Type{DiscreteRoots}, r::Vector{<: Complex}) = DiscreteRoots(r)
Base.convert(::Type{ContinuousRoots}, r::Vector{<: Complex}) = ContinuousRoots(r)

Base.:(≈)(r1::T, r2::T) where T <: AbstractRoots = sort(r1, by=LinearAlgebra.eigsortby) ≈ sort(r2, by=LinearAlgebra.eigsortby)

"""
    domain_transform(d::Domain, e::AbstractRoots)

Change the domain of the roots
"""
domain_transform(d::Discrete,e::ContinuousRoots) = exp(e)
domain_transform(d::Continuous,e::DiscreteRoots) = log(e)
domain_transform(d::Continuous,e::ContinuousRoots) = e
domain_transform(d::Discrete,e::DiscreteRoots) = e
domain(::DiscreteRoots) = Discrete()
domain(::ContinuousRoots) = Continuous()

"""
    change_precision(F, m::AbstractModel)

Changes the precision of all fields in `m` to `F`, e.g., `F=Float64`. This can be useful since the default precision for many operations in this package is `Double64`. This ensures that roots are calculated with high accuracy, but the high precision might not be required to evaluate distances etc.
"""
change_precision(F, r::ContinuousRoots) = ContinuousRoots(Complex{F}.(r.r))
change_precision(F, r::DiscreteRoots) = DiscreteRoots(Complex{F}.(r.r))
change_precision(F) = r->change_precision(F,r)

imageigsortby(λ::Real) = λ
imageigsortby(λ::Complex) = (imag(λ),real(λ))
eigsort(e) = sort(e, by=imageigsortby)
anglesort(e) = sort(e, by=angle)

Base.Vector(r::AbstractRoots) = r.r
function Base.log(r::DiscreteRoots)
    checkroots(r)
    ContinuousRoots(log.(r.r))
end
Base.exp(r::ContinuousRoots) = DiscreteRoots(exp.(r.r))
ControlSystems.c2d(r::DiscreteRoots,h=1) = ContinuousRoots(log(r) ./ h)
d2c(r::ContinuousRoots,h=1) = DiscreteRoots(exp.(h .* r.r))

Lazy.@forward DiscreteRoots.r (Base.length, Base.getindex, Base.setindex!, Base.size, Base.enumerate)
Lazy.@forward ContinuousRoots.r (Base.length, Base.getindex, Base.setindex!, Base.size, Base.enumerate)

for R in (:ContinuousRoots, :DiscreteRoots)
    for ff in [+, -, *, /]
        f = nameof(ff)
        @eval Base.$f(r::$R,a::Number) = $R($f.(r.r, a))
        @eval Base.$f(a::Number,r::$R) = $R($f.(a, r.r))
        @eval Base.$f(r::$R,a::AbstractArray) = $R($f.(r.r, a))
        @eval Base.$f(a::AbstractArray,r::$R) = $R($f.(a, r.r))

        @eval Base.$f(r::$R,a::$R) = $R($f.(r.r, a.r))
    end
end
polar2complex(m,a) = m * cis(a)
function polar(e)
    abs.(e), angle.(e)
end

"""
    polar(e::Number)

magnitude and angle of a complex number
"""
function polar(e::Number)
    abs(e), angle(e)
end

# sqrtmag(magang::Tuple) = sqrt.(magang[1]), magang[2]
# logmag(magang::Tuple) = log.(magang[1]), magang[2]
#
# function logmag(e::Number)
#     m,a = polar(e)
#     m = log(m)
#     polar2complex(m,a)
# end
#
#
# function sqrtreim(e)
#     r,i = real(e), imag(e)
#     complex(sqrt(abs(r))*sign(r), sqrt(abs(i))*sign(i))
# end
#
# function sqrtim(e)
#     r,i = real(e), imag(e)
#     complex(r, sqrt(abs(i))*sign(i))
# end
#
# function sqrtre(e)
#     r,i = real(e), imag(e)
#     complex(sqrt(abs(r))*sign(r), i)
# end
#
# function logreim(e, ϵ=1e-4)
#     r,i = real(e), imag(e)
#     complex(log(abs(r) + ϵ)*sign(r), log(abs(i) + ϵ)*sign(i))
# end
#
# function logim(e, ϵ=1e-4)
#     r,i = real(e), imag(e)
#     complex(r, log(abs(i) + ϵ)*sign(i))
# end
#
# function logre(e, ϵ=1e-4)
#     r,i = real(e), imag(e)
#     complex(log(abs(r) + ϵ)*sign(r), i)
# end
#
# function invre(e)
#     r,i = real(e), imag(e)
#     complex(1/r,i)
# end
#
# function squarere(e)
#     r,i = real(e), imag(e)
#     complex(r*r,i)
# end
#
# function expre(e)
#     r,i = real(e), imag(e)
#     complex(exp(r),i)
# end
#
# function projim(e)
#     complex(0,imag(e))
# end

"""
    residueweight(e::AbstractRoots)

Returns a vector where each entry is roughly corresponding to the amount of energy contributed to the spectrum be each pole. See also [`simplex_residueweight`](@ref) for a normalized version.
"""
function residueweight(e::AbstractRoots)
    # @warn "This method disregards the numerator"
    res = residues(ContinuousRoots(e))
    rw = abs.(π*abs2.(res)./ real.(e))
    isderiving() ? complex.(rw) : rw
end

function residueweight(rw,e::AbstractRoots)
    # @warn "This method disregards the numerator"
    abs2residues!(rw,ContinuousRoots(e))
    @. rw = abs(π*rw/ real(e))
    isderiving() ? complex.(rw) : rw
end

"""
    simplex_residueweight(x)

Returns a vector where each entry is roughly corresponding to the amount of energy contributed to the spectrum be each pole, normalized to sum to 1. See [`residueweight`](@ref) for a non-normalized version.
"""
simplex_residueweight(x) = s1(residueweight(x))
function simplex_residueweight(rw, x)
    residueweight(rw, x)
    rw ./= sum(rw)
    rw
end

"""
    unitweight(e)

A weighting function that returns a vector of uniform weights that sum to 1.
"""
function unitweight(e::AbstractArray{T}) where T
    RT = float(real(T))
    N = length(e)
    isderiving() ? complex.(fill(RT(1/N),size(e))) : fill(RT(1/N),size(e))
end

unitweight(s, e::AbstractArray) = s .= unitweight(e) # support inplace version

toreim(x::AbstractVector{<:Complex}) = (real.(x), imag.(x))
toreim(x::Tuple) = x

reflectc(x::Real) = x < 0 ? x : -x
reflectc(x::Complex) = complex(x.re < 0 ? x.re : -x.re,x.im)
function reflectd(x)
    a = abs(x)
    a < 1 && return oftype(cis(angle(x)),x)
    1/a * cis(angle(x))
end
"""
    reflect(r::AbstractRoots)

Reflects unstable roots to a corresponding stable position (in unit circle for disc. in LHP for cont.)
"""
reflect(r::ContinuousRoots) = ContinuousRoots(reflectc.(r.r))
reflect(r::DiscreteRoots) = DiscreteRoots(reflectd.(r.r))
# function remove_negreal(p)
#     pn = filter(x->x.im == 0 && x.re < 0, p)
#     pn = sort(pn, by=real)
#     i = 1
#     if length(pn) % 2 != 0
#         p[findfirst(==(pn[1]),p)] = 0.01
#         i += 1
#     end
#     while i < length(pn)
#         a = conv([1, -pn[i]], [1, -pn[i+1]])
#         w = sqrt(a[3])
#         @show z = a[2]/(2w)
#         @assert z >= 1
#         z = 0.99
#         @show pi = -z*w + sqrt((z*w)^2-4w^2 + 0*im)
#         p[p .== pn[i]] = pi
#         p[p .== pn[i+1]] = conj(pn)
#         i += 2
#     end
#     p
# end

# scalereal(x) = complex(1x.re, x.im)

"""
    hungariansort(p1, p2)

takes two vectors of numbers and sorts and returns `p2` such that it is in the order of the best Hungarian assignement between `p1` and `p2`. Uses `abs` for comparisons, works on complex numbers.
"""
function hungariansort(p1,p2,p=2)
    length(p1) <= 1 && (return p2)
    distmat = abs.(p1 .- transpose(p2)).^p
    ass,cost = hungarian(distmat)
    p2[ass]
end

function move_real_poles(r::DiscreteRoots, tol=1e-4)
    ro = deepcopy(r)
    reals = findall(abs.(imag.(r)) .< tol)
    @assert length(reals) % 2 == 0 "Uneven number of real poles, I can't solve this for you"
    for i = 1:length(reals)÷2
        r1,r2 = real(r[reals[i]]), real(r[reals[end-i+1]])
        r0 = (r1+r2)/2
        ro[reals[i]] = complex(r0, tol)
        ro[reals[end-i+1]] = complex(r0, -tol)
    end
    DiscreteRoots(ro)
end


text:
abstract type Object end

const INTEGER_OBJ = "INTEGER"
const BOOLEAN_OBJ = "BOOLEAN"
const NULL_OBJ = "NULL"
const RETURN_VALUE = "RETURN_VALUE"
const ERROR_OBJ = "ERROR"
const FUNCTION_OBJ = "FUNCTION"
const STRING_OBJ = "STRING"
const BUILTIN_OBJ = "BUILTIN"
const ARRAY_OBJ = "ARRAY"
const HASH_OBJ = "HASH"
const QUOTE_OBJ = "QUOTE"
const MACRO_OBJ = "MACRO"
const COMPILED_FUNCTION_OBJ = "COMPILED_FUNCTION"
const CLOSURE_OBJ = "CLOSURE"

is_truthy(::Object) = true
Base.show(io::IO, object::Object) = print(io, string(object))

struct IntegerObj <: Object
    value::Int64
end

type_of(::IntegerObj) = INTEGER_OBJ
Base.string(i::IntegerObj) = string(i.value)
Base.:(==)(a::IntegerObj, b::IntegerObj) = a.value == b.value
Base.hash(i::IntegerObj) = hash(i.value)
Object(i::Int) = IntegerObj(i)

struct BooleanObj <: Object
    value::Bool
end

const _TRUE = BooleanObj(true)
const _FALSE = BooleanObj(false)
is_truthy(b::BooleanObj) = b.value
type_of(::BooleanObj) = BOOLEAN_OBJ
Base.string(b::BooleanObj) = string(b.value)
Base.:(==)(a::BooleanObj, b::BooleanObj) = a.value == b.value
Base.hash(b::BooleanObj) = hash(b.value)
Object(b::Bool) = BooleanObj(b)

struct NullObj <: Object end

const _NULL = NullObj()
is_truthy(::NullObj) = false
type_of(::NullObj) = NULL_OBJ
Base.string(::NullObj) = "null"
Object(::Nothing) = NullObj()

struct ReturnValue <: Object
    value::Object
end

type_of(::ReturnValue) = RETURN_VALUE
Base.string(r::ReturnValue) = string(r.value)

struct BreakObj <: Object end

struct ContinueObj <: Object end

struct ErrorObj <: Object
    message::String
end

type_of(::ErrorObj) = ERROR_OBJ
Base.string(e::ErrorObj) = "ERROR: " * e.message

struct Environment
    store::Dict{String,Object}
    outer::Union{Environment,Nothing}
    input::IO
    output::IO

    Environment(; input = stdin, output = stdout) = new(Dict(), nothing, input, output)
    Environment(outer::Environment) = new(Dict(), outer, outer.input, outer.output)
end

get(env::Environment, name::String) = begin
    result = Base.get(env.store, name, nothing)
    if isnothing(result) && !isnothing(env.outer)
        return get(env.outer, name)
    end
    return result
end
set!(env::Environment, name::String, value::Object) = push!(env.store, name => value)
reassign!(env::Environment, name::String, value::Object) = begin
    result = Base.get(env.store, name, nothing)
    if isnothing(result) && !isnothing(env.outer)
        return reassign!(env.outer, name, value)
    end

    # The analyzer ensures that the identifier can always be found.
    set!(env, name, value)
    return value
end

struct FunctionObj <: Object
    parameters::Vector{Identifier}
    body::BlockStatement
    env::Environment
end

type_of(::FunctionObj) = FUNCTION_OBJ
Base.string(f::FunctionObj) =
    "fn(" * join(map(string, f.parameters), ", ") * ") {\n" * string(f.body) * "\n}"

struct StringObj <: Object
    value::String
end

type_of(::StringObj) = STRING_OBJ
Base.string(s::StringObj) = "\"" * string(s.value) * "\""
Base.:(==)(a::StringObj, b::StringObj) = a.value == b.value
Base.hash(s::StringObj) = hash(s.value)
Object(s::String) = StringObj(s)

struct ArrayObj <: Object
    elements::Vector{Object}
end

type_of(::ArrayObj) = ARRAY_OBJ
Base.string(a::ArrayObj) = "[" * join(map(string, a.elements), ", ") * "]"
Base.:(==)(a::ArrayObj, b::ArrayObj) = a.elements == b.elements
Base.hash(a::ArrayObj) = hash(a.elements)
Object(a::Vector) = ArrayObj(map(Object, a))

struct Builtin <: Object
    fn::Function
end

type_of(::Builtin) = BUILTIN_OBJ
Base.string(b::Builtin) = "builtin function"

struct HashObj <: Object
    pairs::Dict{Object,Object}
end

type_of(::HashObj) = HASH_OBJ
Base.string(h::HashObj) =
    "{" * join(map(x -> string(x[1]) * ":" * string(x[2]), collect(h.pairs)), ", ") * "}"
Base.:(==)(a::HashObj, b::HashObj) = a.pairs == b.pairs
Base.hash(h::HashObj) = hash(h.pairs)
Object(h::Dict) = HashObj(Dict(map(x -> Object(x.first) => Object(x.second), collect(h))))

struct QuoteObj <: Object
    node::Node
end

type_of(::QuoteObj) = QUOTE_OBJ
Base.string(q::QuoteObj) = "QUOTE(" * string(q.node) * ")"

struct MacroObj <: Object
    parameters::Vector{Identifier}
    body::BlockStatement
    env::Environment
end

type_of(::MacroObj) = MACRO_OBJ
Base.string(m::MacroObj) =
    "macro(" * join(map(string, m.parameters), ", ") * ") {\n" * string(m.body) * "\n}"

struct CompiledFunctionObj <: Object
    instructions::Instructions
    local_count::Int
    param_count::Int
    within_loop::Bool
end

type_of(::CompiledFunctionObj) = COMPILED_FUNCTION_OBJ
Base.string(c::CompiledFunctionObj) = "compiled function"

struct ClosureObj <: Object
    fn::CompiledFunctionObj
    free::Vector{Object}
end

type_of(::ClosureObj) = CLOSURE_OBJ
Base.string(c::ClosureObj) = "closure"


text:
"""
    transformer_encoder(planes, depth, heads, headplanes, mlppanes; dropout = 0.)

Transformer as used in the base ViT architecture.
([reference](https://arxiv.org/abs/2010.11929)).

# Arguments
- `planes`: number of input channels
- `depth`: number of attention blocks
- `heads`: number of attention heads
- `headplanes`: number of hidden channels per head
- `mlppanes`: number of hidden channels in the MLP block
- `dropout`: dropout rate
"""
function transformer_encoder(planes, depth, heads, headplanes, mlpplanes; dropout = 0.)
  layers = [Chain(SkipConnection(prenorm(planes, MHAttention(planes, headplanes, heads; dropout)), +),
                  SkipConnection(prenorm(planes, mlp_block(planes, mlpplanes; dropout)), +))
            for _ in 1:depth]

  Chain(layers...)
end

"""
    vit(imsize::NTuple{2} = (256, 256); inchannels = 3, patch_size = (16, 16), planes = 1024,
        depth = 6, heads = 16, mlppanes = 2048, headplanes = 64, dropout = 0.1, emb_dropout = 0.1,
        pool = :class, nclasses = 1000)

Creates a Vision Transformer (ViT) model.
([reference](https://arxiv.org/abs/2010.11929)).

# Arguments
- `imsize`: image size
- `inchannels`: number of input channels
- `patch_size`: size of the patches
- `embedplanes`: the number of channels after the patch embedding
- `depth`: number of blocks in the transformer
- `heads`: number of attention heads in the transformer
- `mlpplanes`: number of hidden channels in the MLP block in the transformer
- `headplanes`: number of hidden channels per head in the transformer
- `dropout`: dropout rate
- `emb_dropout`: dropout rate for the positional embedding layer
- `pool`: pooling type, either :class or :mean
- `nclasses`: number of classes in the output
"""
function vit(imsize::NTuple{2} = (256, 256); inchannels = 3, patch_size = (16, 16),
             embedplanes = 768, depth = 6, heads = 16, mlpplanes = 2048, headplanes = 64,
             dropout = 0.1, emb_dropout = 0.1, pool = :class, nclasses = 1000)

  @assert pool in [:class, :mean]
  "Pool type must be either :class (class token) or :mean (mean pooling)"

  npatches = prod(imsize .÷ patch_size)

  return Chain(Chain(PatchEmbedding(imsize; inchannels, patch_size, embedplanes),
                     ClassTokens(embedplanes),
                     ViPosEmbedding(embedplanes, npatches + 1),
                     Dropout(emb_dropout),
                     transformer_encoder(embedplanes, depth, heads, headplanes, mlpplanes; dropout),
                     (pool == :class) ? x -> x[:, 1, :] : seconddimmean),
               Chain(LayerNorm(embedplanes), Dense(embedplanes, nclasses)))
end

struct ViT
  layers
end

"""
    ViT(imsize::NTuple{2} = (256, 256); inchannels = 3, patch_size = (16, 16),
        embedplanes = 768, depth = 6, heads = 16, mlpplanes = 2048, headplanes = 64,
        dropout = 0.1, emb_dropout = 0.1, pool = :class, nclasses = 1000)

Creates a Vision Transformer (ViT) model.
([reference](https://arxiv.org/abs/2010.11929)).

# Arguments
- `imsize`: image size
- `inchannels`: number of input channels
- `patch_size`: size of the patches
- `embedplanes`: the number of channels after the patch embedding
- `depth`: number of blocks in the transformer
- `heads`: number of attention heads in the transformer
- `mlpplanes`: number of hidden channels in the MLP block in the transformer
- `headplanes`: number of hidden channels per head in the transformer
- `dropout`: dropout rate
- `emb_dropout`: dropout rate for the positional embedding layer
- `pool`: pooling type, either :class or :mean
- `nclasses`: number of classes in the output

See also [`Metalhead.vit`](#).
"""
function ViT(imsize::NTuple{2} = (256, 256); inchannels = 3, patch_size = (16, 16),
             embedplanes = 768, depth = 12, heads = 16, mlpplanes = 3072, headplanes = 64,
             dropout = 0.1, emb_dropout = 0.1, pool = :class, nclasses = 1000)

  layers = vit(imsize; inchannels, patch_size, embedplanes, depth, heads, mlpplanes, headplanes,
               dropout, emb_dropout, pool, nclasses)

  ViT(layers)
end

(m::ViT)(x) = m.layers(x)

backbone(m::ViT) = m.layers[1]
classifier(m::ViT) = m.layers[2]

@functor ViT


text:
function ilqrallocs(solver)
    allocs = @allocated J_prev = TO.cost(solver, solver.Z̄)

    # Calculate expansions
    # TODO: do this in parallel
    allocs += @allocated Altro.errstate_jacobians!(solver.model, solver.G, solver.Z̄)
    allocs += @allocated Altro.dynamics_expansion!(solver, solver.Z̄)
    allocs += @allocated Altro.cost_expansion!(solver.obj, solver.Efull, solver.Z̄)
    allocs += @allocated Altro.error_expansion!(solver.model, solver.Eerr, solver.Efull, solver.G, solver.Z̄)

    # Get next iterate
    allocs += @allocated Altro.backwardpass!(solver)
    allocs += @allocated Jnew = Altro.forwardpass!(solver, J_prev)

    # Accept the step and update the current trajectory
    # This is kept out of the forward pass function to make it easier to
    # benchmark the forward pass
    allocs += @allocated copyto!(solver.Z, solver.Z̄)

    # Calculate the gradient of the new trajectory
    dJ = J_prev - Jnew
    grad = Altro.gradient!(solver)

    # Record the iteration
    Altro.record_iteration!(solver, Jnew, dJ, grad)

    # Check convergence
    exit = Altro.evaluate_convergence(solver)
    return allocs
end

@testset "iLQR Solver Allocations" begin
solver = Altro.iLQRSolver(Problems.Cartpole()..., use_static=Val(true))
ilqrallocs(solver)
@test ilqrallocs(solver) == 0
solver = Altro.iLQRSolver(Problems.Cartpole()..., use_static=Val(false))
ilqrallocs(solver)
@test ilqrallocs(solver) == 0
end


text:
mutable struct SubSVDD <: SubOCClassifier
    state::ModelState

    # model parameters
    C::Float64
    kernel_fct::Vector{Kernel}
    subspaces::Vector{Vector{Int}}
    weight_update_strategy::Union{WeightUpdateStrategy, Nothing}

    # training data
    data::Array{Float64,2}
    K::Array{Array{Float64,2}}
    adjust_K::Bool
    K_adjusted::Array{Array{Float64,2}}
    pools::Dict{Symbol, Array{Int64,1}}
    v::Vector{Float64}

    # fitted values
    alpha_values::Vector{Vector{Float64}}
    const_term::Vector{Float64}
    R::Vector{Float64}
    C_delta::Vector{Float64}

    function SubSVDD(data, subspaces, pools)
        m = new()
        m.C = 1.0
        m.state = model_created
        m.data = data
        m.adjust_K = false
        m.pools = labelmap(pools)
        m.v = fill(1.0, size(data,2))
        m.subspaces = subspaces
        m.weight_update_strategy = nothing
        return m
    end
end

SubSVDD(data, pools::Vector{Symbol}) = SubSVDD(data, [], pools)

function get_model_params(model::SubSVDD)
    return Dict(:C => model.C,
                :subspaces => model.subspaces,
                :weight_update_strategy => model.weight_update_strategy)
end

is_valid_param_value(model::SubSVDD, x::Type{Val{:C1}}, v) = 0 <= v <= 1
is_valid_param_value(model::SubSVDD, x::Type{Val{:subspaces}}, v) = true
is_valid_param_value(model::SubSVDD, x::Type{Val{:weight_update_strategy}}, v) = true

function invalidate_solution!(model::SubSVDD)
    model.alpha_values = Vector{Vector{Float64}}()
    model.const_term = Vector{Vector{Float64}}()
    model.R = Vector{Vector{Float64}}()
    model.state = model_initialized
    return nothing
end

function set_C!(model::SubSVDD, C::Real)
    model.C = C
    return nothing
end

function fit!(model::SubSVDD, solver)
    debug(LOGGER, "[FIT] Fitting $(typeof(model)).")
    model.state == model_created && throw(ModelStateException(model.state, model_initialized))
    status = solve!(model, solver)
    r_and_const_term = @eachsubspace get_R_and_const_term(model)
    model.R = map(x -> x[:R], r_and_const_term)
    model.const_term = map(x -> x[:const_term], r_and_const_term)
    model.state = model_fitted
    debug(LOGGER, "[FIT] $(typeof(model)) is now in state $(model.state).")
    return status
end

function solve!(model::SubSVDD, solver::JuMP.OptimizerFactory)
    debug(LOGGER, "[SOLVE] Setting up QP for SubSVDD with $(is_K_adjusted(model) ? "adjusted" : "non-adjusted") kernel matrix.")
    QP = Model(solver)
    K = is_K_adjusted(model) ? model.K_adjusted : model.K

    n = size(model.data, 2)

    @variable(QP, α[1:n, 1:length(K)] >= 0.0)
    @objective(QP, Max, sum(sum(α[i, k]*K[k][i,i] for i in 1:n) -
                            sum(α[i, k] * α[j, k] * K[k][i, j] for i in 1:n for j in 1:n)
                        for k in eachindex(model.subspaces)))
    for i in 1:n
        @constraint(QP, sum(α[i, k] for k in eachindex(model.subspaces)) <= model.v[i] * model.C)
    end

    for k in eachindex(model.subspaces)
        @constraint(QP, sum(α[i, k] for i in 1:n) == 1)
    end

    debug(LOGGER, "[SOLVE] Solving QP with $(typeof(solver))...")
    JuMP.optimize!(QP)
    status = JuMP.termination_status(QP)
    debug(LOGGER, "[SOLVE] Finished with status: $(status).")
    model.alpha_values = [JuMP.value.(α)[:, k] for k in eachindex(model.subspaces)]
    return status
end

```
    Calculates the upper limits for support vectors in subspace subspace_idx
```
function calculate_upper_limit(α::Vector{Vector{Float64}}, C::Float64, v::Vector{Float64}, subspace_idx::Int)
    return v.*C .- sum([α[i] for i in 1:length(α) if i != subspace_idx])
end

function find_support_vectors(model::SubSVDD, subspace_idx)::Vector{Int}
    if length(model.subspaces) == 1
        upper_limits = model.v .* model.C
    else
        upper_limits = SVDD.calculate_upper_limit(model.alpha_values, model.C, model.v, subspace_idx)
    end
    return findall((model.alpha_values[subspace_idx] .> SVDD.OPT_PRECISION) .& (model.alpha_values[subspace_idx] .< (upper_limits .- SVDD.OPT_PRECISION)))
end

function find_positive_alpha(model::SubSVDD, subspace_idx)
    return findall(model.alpha_values[subspace_idx] .> SVDD.OPT_PRECISION)
end

function get_R_and_const_term(model::SubSVDD, subspace_idx)
    sv_idx = find_support_vectors(model, subspace_idx)
    α = model.alpha_values[subspace_idx]
    K = model.K[subspace_idx]

    const_term = sum(α[i]*α[j]*K[i,j] for i in eachindex(α) for j in eachindex(α))
    R_squared = 0.0
    tmp = 0.0
    for k in sv_idx
        tmp = K[k,k] - 2 * sum(α[i] * K[i,k] for i in eachindex(α)) + const_term
        R_squared = tmp > R_squared ? tmp : R_squared
    end
    return (R=sqrt(R_squared), const_term=const_term)
end

#TODO speed-up for in-sample predict by directly using kernel matrix
predict(model::SubSVDD, subspace_idx) = predict(model, model.data[model.subspaces[subspace_idx], :], subspace_idx)

# out of sample predict
"""
    target: data projected in subspace
"""
function predict(model::SubSVDD, target::Array{T,2}, subspace_idx) where T <: Real
    model.state == model_fitted || throw(ModelStateException(model.state, model_fitted))
    @assert size(model.data[model.subspaces[subspace_idx], :], 1) == size(target, 1) "Dimension mismatch between model data and target."
    s = model.subspaces[subspace_idx]
    pos_sv_idx = find_positive_alpha(model, subspace_idx)
    function predict_observation(z)
        kernel(model.kernel_fct[subspace_idx], z, z) -
             2 * sum(model.alpha_values[subspace_idx][i] * kernel(model.kernel_fct[subspace_idx], model.data[s, i], z) for i in pos_sv_idx) +
             model.const_term[subspace_idx]
    end
    vec(sqrt.(mapslices(predict_observation, target, dims=1)) .- model.R[subspace_idx])
end

"""
    target: fullspace data
"""
function predict(model::SubSVDD, target::Array{<:Real,2})
    @assert size(model.data, 1) == size(target, 1)
    map(idx -> predict(model, target[model.subspaces[idx], :], idx), eachindex(model.subspaces))
end

function apply_update_strategy!(model::SubSVDD, new_pools::Vector{Symbol}, query_ids::Vector{Int},
    old_idx_remaining::Vector{Int},
    new_idx_remaining::Vector{Int})

    model.weight_update_strategy === nothing && error("Cannot update $(typeof(model)) with update strategy 'nothing'.")
    @assert length(old_idx_remaining) == length(new_idx_remaining)

    old_v = copy(model.v[old_idx_remaining])
    model.v = fill(get_default_v(model.weight_update_strategy), length(new_pools))
    model.v[new_idx_remaining] .= old_v

    for q_id in query_ids
        model.v[q_id] = update_v(model.v[q_id], new_pools[q_id], model.weight_update_strategy)
    end
    return nothing
end


text:
struct TreeField
    width::UInt32
    height::UInt32
    data::Array{Bool, 1}
end

function fromStdin()::TreeField
    height = 0
    width = 0
    data = Bool[]
    for line in eachline()
        height += 1
        width = length(line)
        for i in firstindex(line):lastindex(line)
            push!(data, line[i] != '.')
        end
    end
    TreeField(width, height, data)
end

function getIndex(field::TreeField, x, y)::Bool
    field.data[(y - 1) * field.width + mod1(x, field.width)]
end

function checkSlope(field::TreeField, xStride, yStride)
    x = 1
    y = 1
    trees = 0
    while y <= field.height
        if getIndex(field, x, y)
            trees += 1
        end
        y += yStride
        x += xStride
    end
    trees
end

function main()
    field = fromStdin()
    println(
        checkSlope(field, 1, 1) *
        checkSlope(field, 3, 1) *
        checkSlope(field, 5, 1) *
        checkSlope(field, 7, 1) *
        checkSlope(field, 1, 2))
end

main()


text:
using JuMP, JuMPeR, Gurobi, Random, Distributions, LinearAlgebra, DataFrames, Plots

n = 10 # Number of facilities
m = 50 # Number of customers

# Generating random data (please don't change the seeds.)
facilities = 0.6.*rand(MersenneTwister(5), n,2) .+ 0.2;
customers = rand(MersenneTwister(2), m, 2);
c = [LinearAlgebra.norm(customers[i, :] .- facilities[j, :])[1] for j=1:n, i=1:m];
f = rand(MersenneTwister(3), n)*1 .+ 5;
s = rand(MersenneTwister(4), n)*2 .+ 15;
d = rand(MersenneTwister(5), m)*0.5 .+ 0.75

# P matrix
R_D = 0.25
P = [0.2*exp(-1/R_D .*LinearAlgebra.norm(customers[i, :] .- customers[j, :])[1]) for j=1:m, i=1:m];
P = (P .>= 0.2*exp(-1/R_D .* R_D)) .* P

""" Nominal facility location model. """
function facility_model(c::Matrix, f::Vector)
    n, m = size(c)
    @assert length(f) == n
    model = Model(solver = GurobiSolver())

    # VARIABLES
    @variable(model, x[1:n], Bin)      # Facility locations
    @variable(model, y[1:n, 1:m] >= 0) # Flow variables
    # CONSTRAINTS
    @constraint(model, [j=1:m], sum(y[i, j] for i = 1:n) >= d[j])      # demand constraint
    @constraint(model, [i=1:n], sum(y[i, j] for j = 1:m) <= s[i]*x[i]) # capacity constraint
    # OBJECTIVE
    @objective(model, Min,
            sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * y[i, j] for i=1:n, j=1:m))
    return model, x, y
end

""" Robust facility location model, with JuMPeR. """
function robust_facility_model(c::Matrix, f::Vector, rho::Real, Gamm::Real)
    n, m = size(c)
    @assert length(f) == n
    model = RobustModel(solver = GurobiSolver())
    # model = Model(solver = GurobiSolver())

    # VARIABLES
    @variable(model, x[1:n], Bin)                 # Facility locations
    @variable(model, y[1:n, 1:m] >= 0)            # Flow variables
    @uncertain(model, -rho <= z[1:m] <= rho)      #  Demand uncertainties
    @constraint(model, norm(z, 1) <= Gamm)

    # CONSTRAINTS
    @constraint(model, [j=1:m], sum(y[i, j] for i = 1:n) >= d[j] + (P*z)[j]) # demand constraint
    @constraint(model, [i=1:n], sum(y[i, j] for j = 1:m) <= s[i] * x[i]) # capacity constraint
    # OBJECTIVE
    @objective(model, Min,
            sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * y[i, j] for i=1:n, j=1:m))
    return model, x, y
end

""" Robust facility location model, with JuMP using the robust counterpart. """
function robust_facility_model(c::Matrix, f::Vector, rho::Real, Gamm::Real)
    n, m = size(c)
    @assert length(f) == n
    model = Model(solver = GurobiSolver())

    # VARIABLES
    @variable(model, x[1:n], Bin)                 # Facility locations
    @variable(model, u[1:n, 1:m] >= 0)            # Flow variables

    # CONSTRAINTS
    for j=1:m
        infdummy = @variable(model)
        normdummy = @variable(model, [1:m])
        y = @variable(model, [1:m])
        @constraint(model, [i=1:m], normdummy[i] >= y[i])
        @constraint(model, [i=1:m], normdummy[i] >= -y[i])
        @constraint(model, [i=1:m], infdummy >= -y[i] - P[j,i])
        @constraint(model, [i=1:m], infdummy >= y[i] + P[j,i])
        @constraint(model, sum(u[i, j] for i = 1:n) >= d[j] + rho*sum(normdummy) + Gamm * infdummy)
    end
    @constraint(model, [i=1:n], sum(u[i, j] for j = 1:m) <= s[i] * x[i]) # capacity constraint
    # OBJECTIVE
    @objective(model, Min,
            sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * u[i, j] for i=1:n, j=1:m))
    return model, x, u
end

""" Adaptive facility location model, with JuMPeR and relaxation of binary variables. """
function adaptive_facility_model(c::Matrix, f::Vector, rho::Real, Gamm::Real)
    n, m = size(c)
    @assert length(f) == n
    model = RobustModel(solver = GurobiSolver())

    # VARIABLES
    @variable(model, 0 <= x[1:n] <= 1)                 # Facility locations
    @uncertain(model, -rho <= z[1:m] <= rho)      #  Demand uncertainties
    @constraint(model, norm(z, 1) <= Gamm)
    @adaptive(model, y[i=1:n, j=1:m] >= 0, policy = Affine, depends_on=z) # Flow variables

    # CONSTRAINTS
    @constraint(model, [j=1:m], sum(y[i, j] for i = 1:n) >= d[j] + (P*z)[j]) # demand constraint
    @constraint(model, [i=1:n], sum(y[i, j] for j = 1:m) <= s[i] * x[i]) # capacity constraint
    # OBJECTIVE
    @variable(model, F)
    @constraint(model, F >= sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * y[i, j] for i=1:n, j=1:m))
    @objective(model, Min, F)
    return model, x, y
end

""" Plots the solution of the facility location model.
Blue circles are active facilities with different capacities.
Orange plus signs are other potential facility locations.
Rays describe connections between facilities and demand nodes.
"""
function plot_solution(model, x, y, cost = nothing)
    plt = scatter(facilities[:, 1], facilities[:, 2], markersize = 0.4 .* s .* getvalue(x))
    scatter!(facilities[:, 1], facilities[:, 2], markersize = 0.4 .* s, markershape = :+)
    for i=1:n
        for j=1:m
            if getvalue(y[i,j]) >= 1e-10
                plot!([customers[j, 1], facilities[i,1]], [customers[j,2], facilities[i,2]], linewidth = getvalue(y[i,j]), legend=false)

            end
        end
    end
    if cost == nothing
        scatter!(customers[:, 1], customers[:, 2], markersize = 3*d,
                title = "Total cost: $(round(getobjectivevalue(model), sigdigits=5))")
    else
        scatter!(customers[:, 1], customers[:, 2], markersize = 3*d,
        title = "Total cost: $(round(cost,sigdigits=5))")
    end
    println("Facility cost: $(getvalue(sum(f[j] * x[j] for j = 1:n)))")
    println("Transportation cost: $(getvalue(sum(c[i, j] * y[i, j] for i=1:n, j=1:m)))")
    return plt
end

# Nominal model (3.1)
model, x, y = facility_model(c, f)
solve(model)
plt = plot_solution(model, x, y)
@show plt

# Robust model (3.2)
rho = 1
Gamm = 5
model, x, y = robust_facility_model(c, f, rho, Gamm)
solve(model)
plt = plot_solution(model, x, y)
@show plt

# Adaptive solution with greedy binary heuristic (3.4 and 3.6)
function iterate_adaptive_model(c, f, rho, Gamm)
    ints = []
    xvals = zeros(size(c, 1))
    while !all([val in [0,1] for val in xvals]) || isempty(ints)
        model, x, y = adaptive_facility_model(c, f, rho, Gamm)
        for int in ints
            @constraint(model, x[int] == 1)
        end
        solve(model)
        xvals = getvalue(x)
        nonint = [!(var in [0,1]) for var in xvals]
        if sum(nonint) == 0
            return model, x, y, ints
        end
        max_val, max_idx = findmax(nonint.*xvals)
        push!(ints, max_idx)
    end
    return
end

model, x, y, ints = iterate_adaptive_model(c, f, rho, Gamm)

# Plot generation for adaptive solution (3.4 and 3.6)
cost = getobjectivevalue(model)
nonints = [val for val in 1:n if !(val in ints)]
plot_model, x, y = facility_model(c, f)
@constraint(plot_model, x[ints] .== 1)
@constraint(plot_model, x[nonints] .== 0)
solve(plot_model)
plt = plot_solution(model, x, y, cost)
@show plt

# Different radii of influence
R_Ds = collect(0:0.05:0.75)
R_Ds[1] = 1e-5

# Unscaled P (3.3)
df = DataFrame("R_D" => [], "sumP" => [], "Gamma" => [], "f+c" => [], "f" => [], "c" => [], "nP" => [], "nx" => [])
for k=1:length(R_Ds)
    global P = [0.2*exp(-1/R_Ds[k] * LinearAlgebra.norm(customers[i, :] .- customers[j, :])[1]) for j=1:m, i=1:m];
    global P = (P .>= 0.2*exp(-1/R_Ds[k] .* R_Ds[k])) .* P
    println("Nonzero Ps: $(sum(P .> 0))")
    model, x, y = robust_facility_model(c, f, rho, Gamm)
    solve(model)
    push!(df, Dict("R_D" => R_Ds[k],
                     "sumP" => sum(P),
                     "Gamma" => Gamm,
                     "f+c" => getobjectivevalue(model),
                     "f" => getvalue(sum(f[j] * x[j] for j = 1:n)),
                     "c" => getvalue(sum(c[i, j] * y[i, j] for i=1:n, j=1:m)),
                     "nP" => sum(P .> 0),
                     "nx" => sum(getvalue(x))))
end

# Column-wise scaled P (3.7)
df = DataFrame("R_D" => [], "sumP" => [], "Gamma" => [], "f+c" => [], "f" => [], "c" => [], "nP" => [], "nx" => [])
P_nominal = [0.2*exp(-1/R_D .*LinearAlgebra.norm(customers[i, :] .- customers[j, :])[1]) for j=1:m, i=1:m];
P_nominal = (P_nominal .>= 0.2*exp(-1/R_D .* R_D)) .* P_nominal
for k=1:length(R_Ds)
    global P = [0.2*exp(-1/R_Ds[k] * LinearAlgebra.norm(customers[i, :] .- customers[j, :])[1]) for j=1:m, i=1:m];
    global P = (P .>= 0.2*exp(-1/R_Ds[k] .* R_Ds[k])) .* P
    for i=1:m
        global P[i,:] = P[i,:] .* sum(P_nominal[i,:]) ./ sum(P[i,:])
    end
    println("Nonzero Ps: $(sum(P .> 0))")
    model, x, y = robust_facility_model(c, f, rho, Gamm)
    solve(model)
    push!(df, Dict("R_D" => R_Ds[k],
                     "sumP" => sum(P),
                     "Gamma" => Gamm,
                     "f+c" => getobjectivevalue(model),
                     "f" => getvalue(sum(f[j] * x[j] for j = 1:n)),
                     "c" => getvalue(sum(c[i, j] * y[i, j] for i=1:n, j=1:m)),
                     "nP" => sum(P .> 0),
                     "nx" => sum(getvalue(x))))
end


text:
println("========================")
println("Preparing Run for Evoqus")
println("========================")
#Make sure the package is installed
try
    Pkg.clone("https://github.com/QuantifyingUncertainty/GeneralizedMetropolisHastings.jl")
catch
    warn("Package already installed in previous run")
end

include("/workspace/GMHPhotoReceptor.jl")

import GeneralizedMetropolisHastings
import GMHPhotoReceptor

@everywhere include("/workspace/GMHPhotoReceptor.jl")


text:
### A Pluto.jl notebook ###
# v0.12.8

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Statistics, DataFrames, CSV
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing4.07"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
data = CSV.read(datadir("purchaseData.csv"), copycols=true)

# ╔═╡ bda159e6-2440-11eb-01a6-f5f4d118214e
mean(data.Price)

# ╔═╡ bda18d2e-2440-11eb-101d-f5c5084a72a4
mean(skipmissing(data.Price))

# ╔═╡ bda23ce4-2440-11eb-0b01-f1e78f176026
coalesce.(data.Grade, "QQ")[1:4]

# ╔═╡ bdac24fc-2440-11eb-0d17-71f2e8eeffc9
first(dropmissing(data,:Price), 4)

# ╔═╡ bdaca91a-2440-11eb-3e00-75b670798eea
sum(ismissing.(data.Name))

# ╔═╡ bdb6f90e-2440-11eb-1e8f-1d6fedde42e0
findall(completecases(data))[1:4]

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing4.07"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═bda159e6-2440-11eb-01a6-f5f4d118214e
# ╠═bda18d2e-2440-11eb-101d-f5c5084a72a4
# ╠═bda23ce4-2440-11eb-0b01-f1e78f176026
# ╠═bdac24fc-2440-11eb-0d17-71f2e8eeffc9
# ╠═bdaca91a-2440-11eb-3e00-75b670798eea
# ╠═bdb6f90e-2440-11eb-1e8f-1d6fedde42e0
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
# Resampling schemes for particle filters
# Some references
#  - http://arxiv.org/pdf/1301.4019.pdf
#  - http://people.isy.liu.se/rt/schon/Publications/HolSG2006.pdf
# Code adapted from: http://uk.mathworks.com/matlabcentral/fileexchange/24968-resampling-methods-for-particle-filtering

# default resampling scheme
resample( w::Vector{Float64}, num_particles::Int = length(w) ) = resampleSystematic(w, num_particles)

function resampleMultinomial( w::Vector{Float64}, num_particles::Int )

  s = Distributions.sampler(Categorical(w))
  indx = rand(s, num_particles)

end

function resampleResidual( w::Vector{Float64}, num_particles::Int )

  M = length( w )

  # "Repetition counts" (plus the random part, later on):
  Ns = floor(length(w) .* w)

  # The "remainder" or "residual" count:
  R = Int(sum( Ns ))

  # The number of particles which will be drawn stocastically:
  M_rdn = num_particles-R;

  # The modified weights:
  Ws = (M .* w - floor(M .* w))/M_rdn;

  # Draw the deterministic part:
  indx1 = Array{Int}(R)
  i=1
  for j=1:M
    for k=1:Ns[j]
      indx1[i]=j
      i = i + 1
    end
  end

  # And now draw the stocastic (Multinomial) part:
  s = Distributions.sampler(Categorical(w))
  indx2 = rand(s, M_rdn)
  indx = append!(indx1, indx2)

  return indx

end

function resampleStratified( w::Vector{Float64}, num_particles::Int )

  N = num_particles
  Q = cumsum(w)

  T = Array{Float64}(N+1)
  for i=1:N,
    T[i] = rand()/N + (i-1)/N
  end
  T[N+1] = 1

  i=1
  j=1

  indx = Array{Int}(N)
  while i<=N
    if T[i]<Q[j]
      indx[i]=j
      i=i+1
    else
      j=j+1
    end
  end

  return indx

end

function resampleSystematic( w::Vector{Float64}, num_particles::Int )

  N = num_particles
  Q = cumsum(w)

  T = collect(linspace(0,1-1/N,N)) + rand()/N
  push!(T, 1)

  i=1
  j=1

  indx = Array{Int}(N)
  while i<=N
    if (T[i]<Q[j])
      indx[i]=j
      i=i+1
    else
      j=j+1
    end
  end

  return indx

end