text:
"""
    AbstractGradientProblem{T} <: Problem{T}

An abstract type for all functions that provide a (full) gradient, where
`T` is a [`AbstractEvaluationType`](@ref) for the gradient function.
"""
abstract type AbstractGradientProblem{T} <: Problem{T} end

@doc raw"""
    GradientProblem{T} <: AbstractGradientProblem{T}

specify a problem for gradient based algorithms.

# Fields
* `M`        – a manifold ``\mathcal M``
* `cost`     – a function ``F: \mathcal M → ℝ`` to minimize
* `gradient!!` – the gradient ``\operatorname{grad}F:\mathcal M → \mathcal T\mathcal M`` of the cost function ``F``.

Depending on the [`AbstractEvaluationType`](@ref) `T` the gradient has to be provided

* as a function `x -> X` that allocates memory for `X` itself for an [`AllocatingEvaluation`](@ref)
* as a function `(X,x) -> X` that work in place of `X` for an [`MutatingEvaluation`](@ref)

# Constructors
    GradientProblem(M, cost, gradient; evaluation=AllocatingEvaluation())

# See also
[`gradient_descent`](@ref), [`GradientDescentOptions`](@ref)
"""
struct GradientProblem{T,mT<:Manifold,C,G} <: AbstractGradientProblem{T}
    M::mT
    cost::C
    gradient!!::G
end
function GradientProblem(
    M::mT, cost::C, gradient::G; evaluation::AbstractEvaluationType=AllocatingEvaluation()
) where {mT<:Manifold,C,G}
    return GradientProblem{typeof(evaluation),mT,C,G}(M, cost, gradient)
end

"""
    get_gradient(p::AbstractGradientProblem{T},x)
    get_gradient!(p::AbstractGradientProblem{T}, X, x)

evaluate the gradient of a [`AbstractGradientProblem{T}`](@ref)`p` at the point `x`.

The evaluation is done in place of `X` for the `!`-variant.
The `T=`[`AllocatingEvaluation`](@ref) problem might still allocate memory within.
When the non-mutating variant is called with a `T=`[`MutatingEvaluation`](@ref)
memory for the result is allocated.
"""
get_gradient(p::AbstractGradientProblem, x)

function get_gradient(p::AbstractGradientProblem{AllocatingEvaluation}, x)
    return p.gradient!!(p.M, x)
end
function get_gradient(p::AbstractGradientProblem{MutatingEvaluation}, x)
    X = zero_tangent_vector(p.M, x)
    return p.gradient!!(p.M, X, x)
end

function get_gradient!(p::AbstractGradientProblem{AllocatingEvaluation}, X, x)
    return copyto!(X, p.gradient!!(p.M, x))
end

function get_gradient!(p::AbstractGradientProblem{MutatingEvaluation}, X, x)
    return p.gradient!!(p.M, X, x)
end

"""
    DirectionUpdateRule

A general functor, that handles direction update rules. It's field(s) is usually
only a [`StoreOptionsAction`](@ref) by default initialized to the fields required
for the specific coefficient, but can also be replaced by a (common, global)
individual one that provides these values.
"""
abstract type DirectionUpdateRule end

"""
    IdentityUpdateRule <: DirectionUpdateRule

The default gradient direction update is the identity, i.e. it just evaluates the gradient.
"""
struct IdentityUpdateRule <: DirectionUpdateRule end

"""
    AbstractGradientOptions <: Options

A generic [`Options`](@ref) type for gradient based options data.
"""
abstract type AbstractGradientOptions <: Options end

@doc raw"""
    GradientDescentOptions{P,T} <: AbstractGradientOptions

Describes a Gradient based descent algorithm, with

# Fields
a default value is given in brackets if a parameter can be left out in initialization.

* `x0` – an a point (of type `P`) on a manifold as starting point
* `gradient` – the current gradient ``\operatorname{grad}f(x)``
* `stopping_criterion` – ([`StopAfterIteration`](@ref)`(100)`) a [`StoppingCriterion`](@ref)
* `stepsize` – ([`ConstantStepsize`](@ref)`(1.)`)a [`Stepsize`](@ref)
* `direction` - ([`IdentityUpdateRule`](@ref)) a processor to compute the gradient
* `retraction_method` – (`ExponentialRetraction()`) the rectraction to use, defaults to
  the exponential map

# Constructor

    GradientDescentOptions(x, stop, s [, retr=ExponentialRetraction()])

construct a Gradient Descent Option with the fields and defaults as above

# See also
[`gradient_descent`](@ref), [`GradientProblem`](@ref)
"""
mutable struct GradientDescentOptions{
    P,TStop<:StoppingCriterion,TStepsize<:Stepsize,TRTM<:AbstractRetractionMethod
} <: AbstractGradientOptions
    x::P
    direction::DirectionUpdateRule
    stop::TStop
    stepsize::TStepsize
    gradient::P
    retraction_method::TRTM
    function GradientDescentOptions{P}(
        initialX::P,
        s::StoppingCriterion=StopAfterIteration(100),
        stepsize::Stepsize=ConstantStepsize(1.0),
        retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
        direction::DirectionUpdateRule=IdentityUpdateRule(),
    ) where {P}
        o = new{P,typeof(s),typeof(stepsize),typeof(retraction_method)}()
        o.x = initialX
        o.stop = s
        o.retraction_method = retraction_method
        o.stepsize = stepsize
        o.direction = direction
        return o
    end
end
function GradientDescentOptions(
    x::P;
    stopping_criterion::StoppingCriterion=StopAfterIteration(100),
    stepsize::Stepsize=ConstantStepsize(1.0),
    retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
    direction::DirectionUpdateRule=IdentityUpdateRule(),
) where {P}
    return GradientDescentOptions{P}(
        x, stopping_criterion, stepsize, retraction_method, direction
    )
end

function (s::IdentityUpdateRule)(p::GradientProblem, o::GradientDescentOptions, i)
    return get_stepsize(p, o, i), get_gradient!(p, o.gradient, o.x)
end

"""
    MomentumGradient <: DirectionUpdateRule

Append a momentum to a gradient processor, where the last direction and last iterate are
stored and the new is composed as ``η_i = m*η_{i-1}' - s d_i``,
where ``sd_i`` is the current (inner) direction and ``η_{i-1}'`` is the vector transported
last direction multiplied by momentum ``m``.

# Fields
* `gradient` – (`zero_tangent_vector(M,x0)`) the last gradient/direction update added as momentum
* `last_iterate` - remember the last iterate for parallel transporting the last direction
* `momentum` – (`0.2`) factor for momentum
* `direction` – internal [`DirectionUpdateRule`](@ref) to determine directions to
  add the momentum to.
* `vector_transport_method` vector transport method to use

# Constructors
    MomentumGradient(
        p::GradientProlem,
        x0,
        s::DirectionUpdateRule=Gradient();
        gradient=zero_tangent_vector(p.M, o.x), momentum=0.2
       vector_transport_method=ParallelTransport(),
    )

Add momentum to a gradient problem, where by default just a gradient evaluation is used
Equivalently you can also use a `Manifold` `M` instead of the [`GradientProblem`](@ref) `p`.

    MomentumGradient(
        p::StochasticGradientProblem
        x0
        s::DirectionUpdateRule=IdentityUpdateRule();
        gradient=zero_tangent_vector(p.M, x0), momentum=0.2
       vector_transport_method=ParallelTransport(),
    )

Add momentum to a stochastic gradient problem, where by default just a stochastic gradient evaluation is used
"""
mutable struct MomentumGradient{P,T,R<:Real,VTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    gradient::T
    last_iterate::P
    momentum::R
    direction::DirectionUpdateRule
    vector_transport_method::VTM
end
function MomentumGradient(
    p::GradientProblem,
    x0::P,
    s::DirectionUpdateRule=IdentityUpdateRule();
    last_iterate=x0,
    vector_transport_method::VTM=ParallelTransport(),
    gradient=zero_tangent_vector(p.M, x0),
    momentum=0.2,
) where {P,VTM<:AbstractVectorTransportMethod}
    return MomentumGradient{P,typeof(gradient),typeof(momentum),VTM}(
        deepcopy(x0), gradient, momentum, s, vector_transport_method
    )
end
function MomentumGradient(
    M::Manifold,
    x0::P,
    s::DirectionUpdateRule=IdentityUpdateRule();
    gradient=zero_tangent_vector(M, x0),
    last_iterate=x0,
    momentum=0.2,
    vector_transport_method::VTM=ParallelTransport(),
) where {P,VTM<:AbstractVectorTransportMethod}
    return MomentumGradient{P,typeof(gradient),typeof(momentum),VTM}(
        gradient, deepcopy(x0), momentum, s, vector_transport_method
    )
end
function (m::MomentumGradient)(p::Problem, o::AbstractGradientOptions, i)
    s, d = m.direction(p, o, i) #get inner direction and step size
    old_d =
        m.momentum *
        vector_transport_to(p.M, m.last_iterate, m.gradient, o.x, m.vector_transport_method)
    m.gradient = old_d - s .* d
    m.last_iterate = deepcopy(o.x)
    return s, -m.gradient
end

"""
    AverageGradient <: DirectionUpdateRule

Add an average of gradients to a gradient processor. A set of previous directions (from the
inner processor) and the last iterate are stored, average is taken after vector transporting
them to the current iterates tangent space.

# Fields
* `gradients` – (fill(`zero_tangent_vector(M,x0),n)`) the last `n` gradient/direction updates
* `last_iterate` – last iterate (needed to transport the gradients)
* `direction` – internal [`DirectionUpdateRule`](@ref) to determine directions to
  apply the averaging to
* `vector_transport_method` - vector transport method to use

# Constructors
    AverageGradient(
        p::GradientProlem,
        x0,
        n::Int=10
        s::DirectionUpdateRule=IdentityUpdateRule();
        gradients = fill(zero_tangent_vector(p.M, o.x),n),
        last_iterate = deepcopy(x0),
        vector_transport_method = ParallelTransport()
    )

Add average to a gradient problem, `n` determines the size of averaging and `gradients` can be prefilled with some history
Equivalently you can also use a `Manifold` `M` instead of the [`GradientProblem`](@ref) `p`.

    AverageGradient(
        p::StochasticGradientProblem
        x0
        n::Int=10
        s::DirectionUpdateRule=IdentityUpdateRule();
        gradients = fill(zero_tangent_vector(p.M, o.x),n),
        last_iterate = deepcopy(x0),
        vector_transport_method = ParallelTransport()
    )

Add average to a stochastic gradient problem, `n` determines the size of averaging and `gradients` can be prefilled with some history
"""
mutable struct AverageGradient{P,T,VTM<:AbstractVectorTransportMethod} <:
               DirectionUpdateRule
    gradients::AbstractVector{T}
    last_iterate::P
    direction::DirectionUpdateRule
    vector_transport_method::VTM
end
function AverageGradient(
    M::Manifold,
    x0::P,
    n::Int=10,
    s::DirectionUpdateRule=IdentityUpdateRule();
    gradients=fill(zero_tangent_vector(M, x0), n),
    vector_transport_method::VTM=ParallelTransport(),
) where {P,VTM<:AbstractVectorTransportMethod}
    return AverageGradient{P,eltype(gradients),VTM}(
        gradients, x0, s, vector_transport_method
    )
end
function AverageGradient(
    p::GradientProblem,
    x0::P,
    n::Int=10,
    s::DirectionUpdateRule=IdentityUpdateRule();
    gradients=fill(zero_tangent_vector(p.M, x0), n),
    vector_transport_method::VTM=ParallelTransport(),
) where {P,VTM}
    return AverageGradient{P,eltype(gradients),VTM}(
        gradients, deepcopy(x0), s, vector_transport_method
    )
end
function (a::AverageGradient)(p::Problem, o::AbstractGradientOptions, i)
    pop!(a.gradients)
    s, d = a.direction(p, o, i) #get inner gradient and step
    a.gradients = vcat([deepcopy(d)], a.gradients)
    for i in 1:(length(a.gradients) - 1) #transport & shift in place
        vector_transport_to!(
            p.M,
            a.gradients[i],
            a.last_iterate,
            a.gradients[i + 1],
            o.x,
            a.vector_transport_method,
        )
    end
    a.gradients[1] = deepcopy(d)
    a.last_iterate = deepcopy(o.x)
    return s, 1 / length(a.gradients) .* sum(a.gradients)
end

@doc raw"""
    Nesterov <: DirectionUpdateRule

## Fields
* `γ`
* `μ` the strong convexity coefficient
* `v` (=``=v_k``, ``v_0=x_0``) an interims point to compute the next gradient evaluation point `y_k`
* `shrinkage` (`= i -> 0.8`) a function to compute the shrinkage ``β_k`` per iterate.

Let's assume ``f`` is ``L``-Lipschitz and ``μ``-strongly convex.
Given
* a step size ``h_k<\frac{1}{L}`` (from the [`GradientDescentOptions`](@ref)
* a `shrinkage` parameter ``β_k``
* and a current iterate ``x_k``
* as well as the interims values ``γ_k`` and ``v_k`` from the previous iterate.

This compute a Nesterov type update using the following steps, see [^ZhangSra2018]

1. Copute the positive root, i.e. ``α_k∈(0,1)`` of ``α^2 = h_k\bigl((1-α_k)γ_k+α_k μ\bigr)``.
2. Set ``\bar γ_k+1 = (1-α_k)γ_k + α_kμ``
3. ``y_k = \operatorname{retr}_{x_k}\Bigl(\frac{α_kγ_k}{γ_k + α_kμ}\operatorname{retr}^{-1}_{x_k}v_k \Bigr)``
4. ``x_{k+1} = \operatorname{retr}_{y_k}(-h_k \operatorname{grad}f(y_k))``
5. ``v_{k+1} = `\operatorname{retr}_{y_k}\Bigl(\frac{(1-α_k)γ_k}{\barγ_k}\operatorname{retr}_{y_k}^{-1}(v_k) - \frac{α_k}{\bar γ_{k+1}}\operatorname{grad}f(y_k) \Bigr)``
6. ``γ_{k+1} = \frac{1}{1+β_k}\bar γ_{k+1}``

Then the direction from ``x_k`` to ``x_k+1``, i.e. ``d = \operatorname{retr}^{-1}_{x_k}x_{k+1}`` is returned.

# Constructor
    Nesterov(x0::P, γ=0.001, μ=0.9, schrinkage = k -> 0.8;
        inverse_retraction_method=LogarithmicInverseRetraction())

Initialize the Nesterov acceleration, where `x0` initializes `v`.

[^ZhangSra2018]:
    > H. Zhang, S. Sra: _Towards Riemannian Accelerated Gradient Methods_,
    > Preprint, 2018, arXiv: [1806.02812](https://arxiv.org/abs/1806.02812)
"""
mutable struct Nesterov{P,T<:Real} <: DirectionUpdateRule
    γ::T
    μ::T
    v::P
    shrinkage::Function
    inverse_retraction_method::AbstractInverseRetractionMethod
end
function Nesterov(
    x0::P,
    γ::T=0.001,
    μ::T=0.9,
    shrinkage::Function=i -> 0.8;
    inverse_retraction_method::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(),
) where {P,T}
    return Nesterov{P,T}(γ, μ, deepcopy(x0), shrinkage, inverse_retraction_method)
end
function (s::Nesterov)(p::GradientProblem, o::AbstractGradientOptions, i)
    h = get_stepsize(p, o, i)
    α = (h * (s.γ - s.μ) + sqrt(h^2 * (s.γ - s.μ)^2 + 4 * h * s.γ)) / 2
    γbar = (1 - α) * s.γ + α * s.μ
    y = retract(p.M, o.x, (α * s.γ) / (s.γ + α * s.μ) .* inverse_retract(p.M, o.x, s.v))
    gradf_yk = get_gradient(p, y)
    xn = retract(p.M, y, -h * gradf_yk)
    d =
        ((1 - α) * s.γ) / γbar .*
        inverse_retract(p.M, y, s.v, s.inverse_retraction_method) - α / γbar .* gradf_yk
    s.v = retract(p.M, y, d, o.retraction_method)
    s.γ = 1 / (1 + s.shrinkage(i)) * γbar
    return h, -1 / h .* inverse_retract(p.M, o.x, xn) # outer update
end

@doc raw"""
    DebugGradient <: DebugAction

debug for the gradient evaluated at the current iterate

# Constructors
    DebugGradient([long=false,p=print])

display the short (`false`) or long (`true`) default text for the gradient.

    DebugGradient(prefix[, p=print])

display the a `prefix` in front of the gradient.
"""
mutable struct DebugGradient <: DebugAction
    io::IO
    prefix::String
    function DebugGradient(long::Bool=false, io::IO=stdout)
        return new(io, long ? "Gradient: " : "gradF(x):")
    end
    DebugGradient(prefix::String, io::IO=stdout) = new(io, prefix)
end
function (d::DebugGradient)(::GradientProblem, o::GradientDescentOptions, i::Int)
    print(d.io, (i >= 0) ? d.prefix * "" * string(o.gradient) : "")
    return nothing
end

@doc raw"""
    DebugGradientNorm <: DebugAction

debug for gradient evaluated at the current iterate.

# Constructors
    DebugGradientNorm([long=false,p=print])

display the short (`false`) or long (`true`) default text for the gradient norm.

    DebugGradientNorm(prefix[, p=print])

display the a `prefix` in front of the gradientnorm.
"""
mutable struct DebugGradientNorm <: DebugAction
    io::IO
    prefix::String
    function DebugGradientNorm(long::Bool=false, io::IO=stdout)
        return new(io, long ? "Norm of the Gradient: " : "|gradF(x)|:")
    end
    DebugGradientNorm(prefix::String, io::IO=stdout) = new(io, prefix)
end
function (d::DebugGradientNorm)(p::GradientProblem, o::Options, i::Int)
    print(d.io, (i >= 0) ? d.prefix * "$(norm(p.M,o.x,o.gradient))" : "")
    return nothing
end

@doc raw"""
    DebugStepsize <: DebugAction

debug for the current step size.

# Constructors
    DebugStepsize([long=false,p=print])

display the short (`false`) or long (`true`) default text for the step size.

    DebugStepsize(prefix[, p=print])

display the a `prefix` in front of the step size.
"""
mutable struct DebugStepsize <: DebugAction
    io::IO
    prefix::String
    function DebugStepsize(long::Bool=false, io::IO=stdout)
        return new(io, long ? "step size:" : "s:")
    end
    DebugStepsize(prefix::String, io::IO=stdout) = new(io, prefix)
end
function (d::DebugStepsize)(p::GradientProblem, o::GradientDescentOptions, i::Int)
    print(d.io, (i > 0) ? d.prefix * "$(get_last_stepsize(p,o,i))" : "")
    return nothing
end

#
# Records
#
@doc raw"""
    RecordGradient <: RecordAction

record the gradient evaluated at the current iterate

# Constructors
    RecordGradient(ξ)

initialize the [`RecordAction`](@ref) to the corresponding type of the tangent vector.
"""
mutable struct RecordGradient{T} <: RecordAction
    recorded_values::Array{T,1}
    RecordGradient{T}() where {T} = new(Array{T,1}())
end
RecordGradient(ξ::T) where {T} = RecordGradient{T}()
function (r::RecordGradient{T})(
    ::GradientProblem, o::GradientDescentOptions, i::Int
) where {T}
    return record_or_reset!(r, o.gradient, i)
end

@doc raw"""
    RecordGradientNorm <: RecordAction

record the norm of the current gradient
"""
mutable struct RecordGradientNorm <: RecordAction
    recorded_values::Array{Float64,1}
    RecordGradientNorm() = new(Array{Float64,1}())
end
function (r::RecordGradientNorm)(
    p::P, o::O, i::Int
) where {P<:GradientProblem,O<:GradientDescentOptions}
    return record_or_reset!(r, norm(p.M, o.x, o.gradient), i)
end

@doc raw"""
    RecordStepsize <: RecordAction

record the step size
"""
mutable struct RecordStepsize <: RecordAction
    recorded_values::Array{Float64,1}
    RecordStepsize() = new(Array{Float64,1}())
end
function (r::RecordStepsize)(
    p::P, o::O, i::Int
) where {P<:GradientProblem,O<:GradientDescentOptions}
    return record_or_reset!(r, get_last_stepsize(p, o, i), i)
end


text:
"""
This `AeolianSources` module is to create aeolian sources for use in AIBECS-generated models.
"""
module AeolianSources

using SparseArrays          # For sparse matrix
using DataDeps              # For storage location of data
using Downloads
using NCDatasets
using Unitful               # for units
using Reexport
using Interpolations
using StringDistances
@reexport using OceanGrids            # To store the grid

function fallback_download(remotepath, localdir)
    @assert(isdir(localdir))
    filename = basename(remotepath)  # only works for URLs with filename as last part of name
    localpath = joinpath(localdir, filename)
    Downloads.download(remotepath, localpath)
    return localpath
end

# Fuzzy name matching
const DATASET_NAMES = ("Chien", "Kok")
const DATASET_LONGNAMES = Dict(
    "Chien" => "Chien_etal_2016",
    "Kok"   => "Kok_etal_2021"
)
rename(s) = findnearest(s, DATASET_NAMES, Levenshtein())[1]
renamelong(s) = DATASET_LONGNAMES[rename(s)]

# Aeolian source URLs
const URLs = Dict(
    "Chien" => "http://www.geo.cornell.edu/eas/PeoplePlaces/Faculty/mahowald/dust/Chienetal2016/post.aerosols.2x2.seasonal.nc",
    "Kok" => "https://dustcomm.atmos.ucla.edu/data/K21b/DustCOMM_source_region_wetdep_annual_PM20_abs.nc"
)
url(s) = URLs[rename(s)]

# Hashes
const SHAs = Dict(
    "Chien" => "33b0036bcca87019e0636ec2d59e910be989693819a2d57f2e4ab053a564101f", # On  2 Jun 2020
    "Kok"   => "65fec31419f75064a22a2c65247afc48e1f72218753605b463447a9933c378aa"  # On 14 Jul 2021
)
sha(s) = SHAs[rename(s)]

datadepname(s) = string("AIBECS-", renamelong(s))

# Create registry entry
function register_AeolianSource(s="Chien")
    register(
        DataDep(
            datadepname(s),
            """
            References:

            $(citations(s))
            """,
            url(s),
            sha(s),
            fetch_method = fallback_download
        )
    )
    return nothing
end

"""
    load()

Returns the 2D aeolian deposition fields from Chien et al. (2016).

You can specify a different dataset via, e.g.,

```julia
load("Kok")
```

At this stage, only two datasets are available:
- `"Chien"` (default) for different dust types (fires, biofuels, dust, ...)
- `"Kok"` for dust from different regions of origin
"""
load(s="Chien") = rename(s) == "Chien" ? load_Chien() : load_Kok()

#=============================#
#     Chien et al data        #
#=============================#
const Chien_AEROSOLTYPE_NAMES = (
    :fire,    # Fires
    :bioF,    # Biofuels
    :FF,      # FossilFuel
    :dust,    # Dust
    :salt,    # Seasalts
    :bioP,    # PrimaryBiogenicParticles
    :volc,    # Volcanoes
)
function load_Chien()
    s = "Chien"
    register_AeolianSource(s)
    nc_file = @datadep_str joinpath(datadepname(s), "post.aerosols.2x2.seasonal.nc")
    s_A_2D = Dataset(nc_file, "r") do ds
        Dict(
            :lat => ds["lat"][:],
            :lon => ds["lon"][:],
            (t => ds["dep"][:,:,:,i] for (i,t) in enumerate(Chien_AEROSOLTYPE_NAMES))...
        )
    end # ds is closed
    @info """You are about to use the Chien et al. (2016) data for aeolian deposition.
          If you use it for research, please cite:

          $(citations(s))

          You can find the corresponding BibTeX entries in the CITATION.bib file
          at the root of the AIBECS.jl package repository.
          (Look for the $(renamelong(s)) key.)
          """
    return s_A_2D
end

#=============================#
#       Kok et al data        #
#=============================#
const Kok_REGIONS_NAMES = (
    :NWAf,    # western North Africa
    :NEAf,    # eastern North Africa
    :Sahel,   # southern Sahara and Sahel
    :MECA,    # Middle East and central Asia
    :EAsia,   # East Asia
    :NAm,     # North America
    :Aus,     # Australia
    :SAm,     # South America
    :SAf,     # southern Africa
)

function load_Kok()
    s = "Kok"
    register_AeolianSource(s)
    nc_file = @datadep_str joinpath(datadepname(s), "DustCOMM_source_region_wetdep_annual_PM20_abs.nc")
# Units                = Total deposition (kg/m2/year)
    s_A_2D = Dataset(nc_file, "r") do ds
        Dict(
            :lat => ds["lat"][:],
            :lon => ds["lon"][:],
            (r => ds["Mean"][i,:,:] for (i,r) in enumerate(Kok_REGIONS_NAMES))...
        )
    end # ds is closed
    @info """You are about to use the Kok et al. (2021) data for aeolian deposition.
          If you use it for research, please cite:

          $(citations(s))

          You can find the corresponding BibTeX entries in the CITATION.bib file
          at the root of the AIBECS.jl package repository.
          (Look for the $(renamelong(s)) key.)
          """
    return s_A_2D
end

#=============================#
#          Citations          #
#=============================#
const CITATIONS = Dict(
    "Chien" => "- Chien, C.-T., K. R. M. Mackey, S. Dutkiewicz, N. M. Mahowald, J. M. Prospero, and A. Paytan (2016), Effects of African dust deposition on phytoplankton in the western tropical Atlantic Ocean off Barbados, Global Biogeochem. Cycles, 30, doi:10.1002/2015GB005334.",
    "Kok" => """
    - Description: Kok, J. F., Adebiyi, A. A., Albani, S., Balkanski, Y., Checa-Garcia, R., Chin, M., Colarco, P. R., Hamilton, D. S., Huang, Y., Ito, A., Klose, M., Li, L., Mahowald, N. M., Miller, R. L., Obiso, V., Pérez García-Pando, C., Rocha-Lima, A., and Wan, J. S.: Contribution of the world's main dust source regions to the global cycle of desert dust, Atmos. Chem. Phys., 21, 8169–8193, https://doi.org/10.5194/acp-21-8169-2021, 2021.
    - Dataset: Adebiyi, A. A., Kok, J. F., Wang, Y., Ito, A., Ridley, D. A., Nabat, P., and Zhao, C.: Dust Constraints from joint Observational-Modelling-experiMental analysis (DustCOMM): comparison with measurements and model simulations, Atmos. Chem. Phys., 20, 829–863, https://doi.org/10.5194/acp-20-829-2020, 2020.
    - Source apportionment framework: Kok, J. F., Adebiyi, A. A., Albani, S., Balkanski, Y., Checa-Garcia, R., Chin, M., Colarco, P. R., Hamilton, D. S., Huang, Y., Ito, A., Klose, M., Leung, D. M., Li, L., Mahowald, N. M., Miller, R. L., Obiso, V., Pérez García-Pando, C., Rocha-Lima, A., Wan, J. S., and Whicker, C. A.: Improved representation of the global dust cycle using observational constraints on dust properties and abundance, Atmos. Chem. Phys., 21, 8127–8167, https://doi.org/10.5194/acp-21-8127-2021, 2021.
    - Climate model used: Scanza, R. A., Hamilton, D. S., Perez Garcia-Pando, C., Buck, C., Baker, A., and Mahowald, N. M.: Atmospheric processing of iron in mineral and combustion aerosols: development of an intermediate-complexity mechanism suitable for Earth system models, Atmos. Chem. Phys., 18, 14175–14196, https://doi.org/10.5194/acp-18-14175-2018, 2018.
    """
)
citations(s="Chien") = CITATIONS[rename(s)]

end # end module

export AeolianSources

ext: jl


text:
@testset "Testing decoder.jl" begin

    @testset "Testing decoder and objective" begin
        instance = read_dataset("datasets/movie.toi");
        bucket = [0 0;
                  0 0;
                  0 0;
                  0 0;
                  0 0;
                  0 0]
        keys = [0.16666666666666666 0.159;
                0.3333333333333333  0.32;
                0.5                 0.49;
                0.6666666666666666  0.65;
                0.8333333333333334  0.82;
                1.0                 0.84]

        # Test decoder
        AbrkgaObop.decoder!(1,bucket,keys)
        @test bucket[:,1] == [1,2,3,4,5,6]
        AbrkgaObop.decoder!(2,bucket,keys)
        @test bucket[:,2] == [1,2,3,4,5,6]

        # Test objective
        fitness=[0,0]
        AbrkgaObop.objective!(1,bucket,fitness,instance)
        AbrkgaObop.objective!(2,bucket,fitness,instance)
        @test fitness[1] == fitness[2]
    end;

end

max_stars_repo_name: LuizHNLorena/AbrkgaObop

max_issues_repo_name: LuizHNLorena/AbrkgaObop

max_forks_repo_name: LuizHNLorena/AbrkgaObop


text:
using ProjectAssistant, Test

# TODO


text:
using Flux:onecold
using Test

@testset "onecold" begin
  a = [1, 2, 5, 3.]
  A = [1 20 5; 2 7 6; 3 9 10; 2 1 14]
  labels = ['A', 'B', 'C', 'D']

  @test onecold(a) == 3
  @test onecold(A) == [3, 1, 4]
  @test onecold(a, labels) == 'C'
  @test onecold(A, labels) == ['C', 'A', 'D']

  data = [:b, :a, :c]
  labels = [:a, :b, :c]
  hot = Flux.onehotbatch(data, labels)
  cold = onecold(hot, labels)

  @test cold == data
end

@testset "onehotbatch indexing" begin
  y = Flux.onehotbatch(ones(3), 1:10)
  @test y[:,1] isa Flux.OneHotVector
  @test y[:,:] isa Flux.OneHotMatrix
end

@testset "abstractmatrix onehotvector multiplication" begin
  A = [1 3 5; 2 4 6; 3 6 9]
  b1 = Flux.OneHotVector(1, 3)
  b2 = Flux.OneHotVector(3, 5)

  @test A*b1 == A[:,1]
  @test_throws DimensionMismatch A*b2
end

@testset "OneHotArray" begin
  using Flux: OneHotArray, OneHotVector, OneHotMatrix, OneHotLike

  ov = OneHotVector(rand(1:10), 10)
  om = OneHotMatrix(rand(1:10, 5), 10)
  oa = OneHotArray(rand(1:10, 5, 5), 10)

  # sizes
  @testset "Base.size" begin
    @test size(ov) == (10,)
    @test size(om) == (10, 5)
    @test size(oa) == (10, 5, 5)
  end

  @testset "Indexing" begin
    # vector indexing
    @test ov[3] == (ov.indices == 3)
    @test ov[:] == ov

    # matrix indexing
    @test om[3, 3] == (om.indices[3] == 3)
    @test om[:, 3] == OneHotVector(om.indices[3], 10)
    @test om[3, :] == (om.indices .== 3)
    @test om[:, :] == om

    # array indexing
    @test oa[3, 3, 3] == (oa.indices[3, 3] == 3)
    @test oa[:, 3, 3] == OneHotVector(oa.indices[3, 3], 10)
    @test oa[3, :, 3] == (oa.indices[:, 3] .== 3)
    @test oa[3, :, :] == (oa.indices .== 3)
    @test oa[:, 3, :] == OneHotMatrix(oa.indices[3, :], 10)
    @test oa[:, :, :] == oa

    # cartesian indexing
    @test oa[CartesianIndex(3, 3, 3)] == oa[3, 3, 3]
  end

  @testset "Concatenating" begin
    # vector cat
    @test hcat(ov, ov) == OneHotMatrix(vcat(ov.indices, ov.indices), 10)
    @test vcat(ov, ov) == vcat(collect(ov), collect(ov))
    @test cat(ov, ov; dims = 3) == OneHotArray(cat(ov.indices, ov.indices; dims = 2), 10)

    # matrix cat
    @test hcat(om, om) == OneHotMatrix(vcat(om.indices, om.indices), 10)
    @test vcat(om, om) == vcat(collect(om), collect(om))
    @test cat(om, om; dims = 3) == OneHotArray(cat(om.indices, om.indices; dims = 2), 10)

    # array cat
    @test cat(oa, oa; dims = 3) == OneHotArray(cat(oa.indices, oa.indices; dims = 2), 10)
    @test cat(oa, oa; dims = 1) == cat(collect(oa), collect(oa); dims = 1)
  end

  @testset "Base.reshape" begin
    # reshape test
    @test reshape(oa, 10, 25) isa OneHotLike
    @test reshape(oa, 10, :) isa OneHotLike
    @test reshape(oa, :, 25) isa OneHotLike
    @test reshape(oa, 50, :) isa OneHotLike
    @test reshape(oa, 5, 10, 5) isa OneHotLike
    @test reshape(oa, (10, 25)) isa OneHotLike

    @testset "w/ cat" begin
      r = reshape(oa, 10, :)
      @test hcat(r, r) isa OneHotArray
      @test vcat(r, r) isa Array{Bool}
    end

    @testset "w/ argmax" begin
      r = reshape(oa, 10, :)
      @test argmax(r) == argmax(OneHotMatrix(reshape(oa.indices, :), 10))
      @test Flux._fast_argmax(r) == collect(reshape(oa.indices, :))
    end
  end

  @testset "Base.argmax" begin
    # argmax test
    @test argmax(ov) == argmax(convert(Array{Bool}, ov))
    @test argmax(om) == argmax(convert(Array{Bool}, om))
    @test argmax(om; dims = 1) == argmax(convert(Array{Bool}, om); dims = 1)
    @test argmax(om; dims = 2) == argmax(convert(Array{Bool}, om); dims = 2)
    @test argmax(oa; dims = 1) == argmax(convert(Array{Bool}, oa); dims = 1)
    @test argmax(oa; dims = 3) == argmax(convert(Array{Bool}, oa); dims = 3)
  end
end


text:
# integration kernel for loop function
function compute_dΣ_kernel(
    Λ :: Float64,
    w :: Float64,
    v :: Float64,
    r :: Reduced_lattice,
    m :: Mesh,
    a :: Action_su2
    ) :: Float64

    # get buffers for non-local vertex
    b1s = get_buffer_s( v + w, 0.5 * (-v + w), 0.5 * (v - w), m)
    b1t = get_buffer_t(0.0, w, v, m)
    b1u = get_buffer_u(-v + w, 0.5 * ( v + w), 0.5 * (v + w), m)

    # get buffers for local vertex
    b2s = get_buffer_s(v + w, 0.5 * (-v + w), 0.5 * (-v + w), m)
    b2t = get_buffer_t(v - w, 0.5 * ( v + w), 0.5 * ( v + w), m)
    b2u = get_buffer_u(0.0, w, v, m)

    # compute local contributions
    val = 3.0 * get_Γ_comp(1, 1, b2s, b2t, b2u, r, a, apply_flags_su2) + get_Γ_comp(2, 1, b2s, b2t, b2u, r, a, apply_flags_su2)

    # compute contributions for all lattice sites
    for j in eachindex(r.sites)
        val -= 2.0 * r.mult[j] * (2.0 * a.S)  * get_Γ_comp(2, j, b1s, b1t, b1u, r, a, apply_flags_su2)
    end

    # multiply with single scale propagator
    val *= get_S(Λ, v, m, a) / (2.0 * pi)

    return val
end

# first integration kernel for loop function for self energy corrections
function compute_dΣ_kernel_corr1(
    Λ    :: Float64,
    w    :: Float64,
    v    :: Float64,
    r    :: Reduced_lattice,
    m    :: Mesh,
    a    :: Action_su2,
    da_Σ :: Action_su2
    )    :: Float64

    # get buffers for non-local vertex
    b1s = get_buffer_s( v + w, 0.5 * (-v + w), 0.5 * (v - w), m)
    b1t = get_buffer_t(0.0, w, v, m)
    b1u = get_buffer_u(-v + w, 0.5 * ( v + w), 0.5 * (v + w), m)

    # get buffers for local vertex
    b2s = get_buffer_s(v + w, 0.5 * (-v + w), 0.5 * (-v + w), m)
    b2t = get_buffer_t(v - w, 0.5 * ( v + w), 0.5 * ( v + w), m)
    b2u = get_buffer_u(0.0, w, v, m)

    # compute local contributions
    val = 3.0 * get_Γ_comp(1, 1, b2s, b2t, b2u, r, da_Σ, apply_flags_su2, ch_u = false) + get_Γ_comp(2, 1, b2s, b2t, b2u, r, da_Σ, apply_flags_su2, ch_u = false)

    # compute contributions for all lattice sites
    for j in eachindex(r.sites)
        val -= 2.0 * r.mult[j] * (2.0 * a.S) * get_Γ_comp(2, j, b1s, b1t, b1u, r, da_Σ, apply_flags_su2, ch_t = false)
    end

    # multiply with full propagator
    val *= -get_G(Λ, v, m, a) / (2.0 * pi)

    return val
end

# second integration kernel for loop function for self energy corrections
function compute_dΣ_kernel_corr2(
    Λ    :: Float64,
    w    :: Float64,
    v    :: Float64,
    r    :: Reduced_lattice,
    m    :: Mesh,
    a    :: Action_su2,
    da_Σ :: Action_su2
    )    :: Float64

    # get buffers for non-local vertex
    b1s = get_buffer_s( v + w, 0.5 * (-v + w), 0.5 * (v - w), m)
    b1t = get_buffer_t(0.0, w, v, m)
    b1u = get_buffer_u(-v + w, 0.5 * ( v + w), 0.5 * (v + w), m)

    # get buffers for local vertex
    b2s = get_buffer_s(v + w, 0.5 * (-v + w), 0.5 * (-v + w), m)
    b2t = get_buffer_t(v - w, 0.5 * ( v + w), 0.5 * ( v + w), m)
    b2u = get_buffer_u(0.0, w, v, m)

    # compute local contributions
    val = 3.0 * get_Γ_comp(1, 1, b2s, b2t, b2u, r, a, apply_flags_su2) + get_Γ_comp(2, 1, b2s, b2t, b2u, r, a, apply_flags_su2)

    # compute contributions for all lattice sites
    for j in eachindex(r.sites)
        val -= 2.0 * r.mult[j] * (2.0 * a.S) * get_Γ_comp(2, j, b1s, b1t, b1u, r, a, apply_flags_su2)
    end

    # multiply with full propagator
    val *= get_G(Λ, v, m, a)^2 * get_Σ(v, m, da_Σ) / (2.0 * pi)

    return val
end


text:
using Zygote, LinearAlgebra
using PyCall
# Make sure pytorch is installed

torch = pyimport("torch")

# From Mike Innes' JuliaCon 2019 talk
# https://youtu.be/OcUXjk7DFvU?t=1106
Zygote.@adjoint function pycall(f, x...; kw...)
    x = map(py, x)
    y = pycall(f, x...; kw...)
    y.detach().numpy(), function (ȳ)
        y.backward(gradient = py(ȳ))
        (nothing, map(x->x.grad.numpy(), x)...)
    end
end

# Non-symmetric version, which gives complex eigenvalues
# function torch_eig(A)
#     tA = torch.from_numpy(A)
#     out = torch.eig(tA)
#     evs = out[1].numpy()
#     return [ evs[j, 1] + im*evs[j,2] for j = 1:size(evs,1) ]
# end

function torch_symeig(A)
    tA = torch.from_numpy(A)
    out = torch.symeig(tA)
    evs = out[1].numpy()
    return evs
end

function test_eigmin(A)
    evals = torch_symeig(A)
    return minimum(evals)
end

A = rand(4, 4); A = A + A';
@show eigmin(A)
@show test_eigmin(A)
@show Zygote.gradient(test_eigmin, A)

max_stars_repo_path: Methods/PyTorch via PyCall.jl

max_stars_repo_name: ericphanson/Attempt-diff-eigmin

max_issues_repo_path: Methods/PyTorch via PyCall.jl

max_issues_repo_name: ericphanson/Attempt-diff-eigmin

max_forks_repo_path: Methods/PyTorch via PyCall.jl

max_forks_repo_name: ericphanson/Attempt-diff-eigmin


text:
module TestReturns

using Base.Test
using TimeData
using Econometrics
using DataFrames
using Dates

## println("\n Running return function tests\n")

#########################
## disc2log / log2disc ##
#########################

## percentage returns
##-------------------

ret = Timematr(DataFrame(returns = [2, 3]))

## log2disc(disc2log(x)) = x
ret2 = Econometrics.log2disc(Econometrics.disc2log(ret, percent=true),
                             percent=true)

@test_approx_eq TimeData.core(ret) TimeData.core(ret2)

## no percentage returns
##----------------------

ret = Timematr(DataFrame(returns = [0.2, 0.3]))
ret2 = Econometrics.log2disc(Econometrics.disc2log(ret, percent=false),
                             percent=false)

@test_approx_eq TimeData.core(ret) TimeData.core(ret2)

###############
## price2ret ##
###############

prices = testcase(TimeData.Timenum, 2)
logRet = Econometrics.price2ret(prices, log = true)

## manually determine expected result
##-----------------------------------

idxs = prices.idx[2:end]
expTn = TimeData.Timenum(DataFrame(prices1 = @data([NA, 20, 30, 30]),
                                   prices2 = @data([10, NA, 10, NA])),
                         idxs)

## test
@test isequal(logRet, expTn)

## test Timematr
##--------------

prices = testcase(TimeData.Timematr, 1)
logRet = Econometrics.price2ret(prices, log = true)

idxs = prices.idx[2:end]
expTn = TimeData.Timematr(DataFrame(prices1 = [20, -10, 60],
                                    prices2 = [10, -20, 30]),
                          idxs)

## test
@test isequal(logRet, expTn)

## same for discrete returns
##--------------------------

prices = TimeData.testcase(TimeData.Timenum, 2)
discRet = Econometrics.price2ret(prices, log = false)

@test_approx_eq ((140-120) / 120) TimeData.core(discRet[2, 1])
@test_approx_eq ((120-110) / 110) TimeData.core(discRet[1, 2])
@test isna(TimeData.core(discRet[2, 2])[1])

###############
## ret2price ##
###############

##################################
## 1. example: NAs at beginning ##
##################################

prices = @data([NA, NA, 3, 4, 6])

## using log returns
##------------------

logRet = Econometrics.price2ret(prices, log = true)
expRets = @data([NA, NA, 1, 2])
@test isequal(logRet, expRets)

normedPrices = Econometrics.ret2price(logRet, log = true)
expNormed = @data([NA, NA, 0, 1, 3])
@test isequal(normedPrices, expNormed)

## using dicrete returns
##----------------------

discRet = Econometrics.price2ret(prices, log = false)
expRets = @data([NA, NA, 1/3, 0.5])
@test isequal(discRet, expRets)

normedPrices = Econometrics.ret2price(discRet, log = false)
expNormed = @data([NA, NA, 1, 1+1/3, 2])
@test isequal(normedPrices, expNormed)

################################
## 2. example: NAs in between ##
################################

prices = @data([110, 120, NA, 130, NA])

## using log returns
##------------------

logRet = Econometrics.price2ret(prices, log = true)
expRets = @data([10, NA, 10, NA])
@test isequal(logRet, expRets)

normedPrices = Econometrics.ret2price(logRet, log = true)
expNormed = @data([0, 10, NA, 20, NA])
@test isequal(normedPrices, expNormed)

## using dicrete returns
##----------------------

discRet = Econometrics.price2ret(prices, log = false)
expRets = @data([10/110, NA, 10/120, NA])
@test isequal(discRet, expRets)

normedPrices = Econometrics.ret2price(discRet, log = false)
expNormed = @data([1, 1+10/110, NA, 1+20/110, NA])
@test_approx_eq normedPrices[[1,2,4]] expNormed[[1,2,4]]

########################
## 3. example: no NAs ##
########################

prices = @data([100, 120, 110, 170, 160])

## using log returns
##------------------

logRet = Econometrics.price2ret(prices, log = true)
expRets = @data([20, -10, 60, -10])
@test isequal(logRet, expRets)

normedPrices = Econometrics.ret2price(logRet, log = true)
expNormed = @data([0, 20, 10, 70, 60])
@test isequal(normedPrices, expNormed)

## using dicrete returns
##----------------------

discRet = Econometrics.price2ret(prices, log = false)
expRets = @data([0.2, -10/120, 60/110, -10/170])
@test isequal(discRet, expRets)

normedPrices = Econometrics.ret2price(discRet, log = false)
expNormed = @data([1, 1.2, 1.1, 1.7, 1.6])
@test_approx_eq normedPrices[[1,2,4]] expNormed[[1,2,4]]

############################
## same tests for Timenum ##
############################

df = DataFrame()
df[:prices1] = @data([100, 120, 110, 170, 160])
df[:prices2] = @data([110, 120, NA, 130, NA])
df[:prices3] = @data([NA, NA, 3, 4, 6])

dats = [Date(2010,1,1):Date(2010,1,5)]
prices = Timenum(df, dats)

## using log returns
##------------------

logRet = Econometrics.price2ret(prices, log = true)
normedPrices = Econometrics.ret2price(logRet, log = true)
originalPrices = Econometrics.ret2price(logRet, prices, log = true)
@test isequal(prices, originalPrices)

## using dicrete returns
##----------------------

discRet = Econometrics.price2ret(prices, log = false)
normedPrices = Econometrics.ret2price(discRet, log = false)
originalPrices = Econometrics.ret2price(discRet, prices, log = false)
@test isapprox(prices, originalPrices)

########################
## tests for Timematr ##
########################

df = DataFrame()
df[:prices1] = @data([100, 120, 110, 170, 160])
dats = [Date(2010,1,1):Date(2010,1,5)]
prices = Timenum(df, dats)

## using log returns
##------------------

logRet = Econometrics.price2ret(prices, log = true)
normedPrices = Econometrics.ret2price(logRet, log = true)
originalPrices = Econometrics.ret2price(logRet, prices, log = true)
@test isequal(prices, originalPrices)

## using dicrete returns
##----------------------

discRet = Econometrics.price2ret(prices, log = false)
normedPrices = Econometrics.ret2price(discRet, log = false)
originalPrices = Econometrics.ret2price(discRet, prices, log = false)
@test isapprox(prices, originalPrices)

end


text:
function eigenvec()
  const io = open(joinpath("output", "zvo_eigenvec_0_rank_0.dat"))
  const N = read(io, Int)
  res = zeros(Complex128, N)
  @inbounds for i in 1:N
    re = read(io, Float64)
    img = read(io, Float64)
    res[i] = complex(re,img)
  end
  return res
end

function gauge_fix!(vec, ref)
  return vec
  ip = cis(angle(dot(vec, ref)))
  @inbounds for i in 1:length(vec)
    vec[i] *= ip
  end
  vec
end

function normalize!(vec)
  s = sum(abs2(vec))
  s = 1.0/sqrt(s)
  @inbounds for i in 1:length(vec)
    vec[i] *= s
  end
end

max_stars_repo_name: yomichi/HPhi-samples

max_issues_repo_name: yomichi/HPhi-samples

max_forks_repo_name: yomichi/HPhi-samples


text:
"""
Evaluators are used to assess the success or performance of candidate solutions
on some predetermined set of search objectives.

Unlike most other evolutionary computation frameworks, evaluators do not
immediately calculate the fitness values for individuals. Instead, the raw
objective scores for a given set of individuals is first calculated and
stored, before the fitness values for each individual are computed according
to those values according to some fitness scheme.

This difference in treatment makes implementing dynamic problems, co-evolution,
diversity promotion techniques, multi-objective problems, and fuzzy evaluation
much easier!
"""
module evaluator
importall common
using criterion, state, core, utility, fitness, population, _deme_
export Evaluator, EvaluatorDefinition, evaluate!, call

"""
The base type used by all evaluators.
"""
abstract Evaluator

"""
The base type used by all evaluator definitions.
"""
abstract EvaluatorDefinition

"""
Evaluates a given phenome using a provided fitness scheme, according to a
specified evaluation method. This function should be implemented by each
evaluator.
"""
evaluate!(ev::Evaluator, ::FitnessScheme, ::Any) =
  error("Unimplemented evaluate!(::Evaluator, ::FitnessScheme, phenome) for evaluator: $(typeof(ev)).")

"""
Evaluates all unevaluated individuals within a provided deme according to
a provided evaluation method.
"""
function evaluate!(ev::Evaluator, s::State, d::Deme)
  for (id, phenome) in enumerate(d.offspring.stages[ev.stage])
    d.offspring.fitnesses[id] = evaluate!(ev, d.species.fitness, get(phenome))
  end
  s.evaluations += length(d.offspring.stages[ev.stage])
end

"""
Evaluates all unevaluated individuals contained within a given state according
to a provided evaluation method..
"""
evaluate!(ev::Evaluator, s::State) =
  for deme in s.population.demes; evaluate!(ev, s, deme); end

# Load each of the evaluators.
include("evaluator/_simple.jl")
include("evaluator/_tsp.jl")
include("evaluator/_regression.jl")
#include("evaluator/_multiplexer.jl")
#include("evaluator/_ant.jl")
end


text:
import Base.Sort: select

"""
$(SIGNATURES)

Given `cells`, a vector of integers, and `guides`, a vector of barcodes,
performs simulated FACS sorting of the cells into `bins` with the given
cutoffs. `σ` refers to the spread of the phenotype during the FACS screen.
"""
function select(setup::FacsScreen,
                cells::Vector{Int},
                cell_phenotypes::Vector{Float64},
                guides::Vector{Barcode};
                debug = false)

    σ = setup.σ
    bins = setup.bin_info
    @assert size(cells) == size(cell_phenotypes)
    n_cells = length(cells)
    observed = zeros(n_cells)
    @inbounds for i in 1:n_cells
        observed[i] = rand(Normal(cell_phenotypes[i], σ))
        guides[cells[i]].obs_phenotype = observed[i]
    end
    indices = sortperm(observed)
    cells = cells[indices]
    results = OrderedDict{Symbol, Vector{Int}}()

    for (binname, cutoffs) in bins
        left = clamp(round(Int, cutoffs[1]*n_cells), 1, n_cells)
        right = clamp(round(Int, cutoffs[2]*n_cells), 1, n_cells)
        results[binname] = cells[left:right]
    end
    results
end

function grow!(cells::AbstractArray{Int},
               cell_phenotypes::AbstractArray{Float64},
               output_c::AbstractArray{Int},
               output_p::AbstractArray{Float64},
               setup::GrowthScreen)
    num_inserted::Int = 0
    noise_dist = Normal(0, setup.noise)
    @inbounds for i in 1:length(cells)
        ρ::Float64 = cell_phenotypes[i]
        ρ_noisy = ρ + rand(noise_dist)
        decision = abs(ρ_noisy) < rand() ? 2 : 2^trunc(Int, 1 + sign(ρ_noisy))
        rng = num_inserted+1:num_inserted+decision
        output_c[rng] = cells[i]
        output_p[rng] = ρ
        num_inserted+=decision
    end
    num_inserted
end

"""
$(SIGNATURES)

Growth Screen selection
"""
function select(setup::GrowthScreen,
                initial_cells::AbstractArray{Int},
                initial_cell_phenotypes::AbstractArray{Float64},
                guides::Vector{Barcode};
                debug=false)

    bottleneck_representation = setup.bottleneck_representation
    num_bottlenecks = setup.num_bottlenecks
    # all cells at all timepoints
    cellmat = zeros(Int, length(guides)*bottleneck_representation)
    cpmat = zeros(Float64, size(cellmat))
    output_c = Array{Int}(length(initial_cells)*4);
    output_p = Array{Float64}(size(output_c));
    cells = initial_cells # 1st timepoint slice
    picked = Array{Int}(size(cellmat, 1))
    cell_phenotypes = initial_cell_phenotypes

    for k in 1:num_bottlenecks
        num_inserted = grow!(cells, cell_phenotypes, output_c, output_p, setup)
        cells, cell_phenotypes = view(cellmat, :), view(cpmat, :)
        sample!(1:num_inserted, picked, replace=false)
        copy!(view(cellmat, :), view(output_c, picked))
        copy!(view(cpmat, :), view(output_p, picked))
    end
    return OrderedDict(:bin1 => initial_cells, :bin2 => cellmat)
end


text:
#https://groups.google.com/forum/#!topic/julia-users/6p5vRSYcApY
#https://github.com/dpo/ampl.jl
#http://orfe.princeton.edu/~rvdb/ampl/nlmodels/index.html
#https://github.com/mpf/Optimization-Test-Problems (many of the same problems, without the solve command)
#http://netlib.org/ampl/models/

module Lp_bench
  using Lp
  using Jad
  importall Arg

  type DbProblem
    creator::Function
    descr::String
    check_status::Symbol
    check_z#::Array{typestr, 1}
    check_x#::Array{typestr, 1}
  end

  prob_db = Function[]
  prob_last = Array(Any, 0)

  function problem_t1(params::Params)
    return Lp.create_min_Cf0_problem(merge(params, {"maxit" => 5}),
      [1, 1], [2 3;], [10]
      )
  end; push!(prob_db, DbProblem(problem_t1, "t1", :Optimal, :Unset, [0, 0]))

  function problem_LPFE_p88(params::Params)
    return Lp.create_max_Cf0_problem(merge(params, {"maxit" => 5}),
      [4,3], [1 -1; 2 -1; 0 1], [1, 3, 5]
      )
  end; push!(prob_db, DbProblem(problem_LPFE_p88, "LPFE_p88", :Optimal, 31, [4, 5]))

  function problem_LPFE_p11(params::Params)
    return Lp.create_max_Cf0_problem(merge(params, {"maxit" => 5}),
      [5,4,3], [2 3 1; 4 1 2; 3 4 2], [5, 11, 8]
      )
  end; push!(prob_db, DbProblem(problem_LPFE_p11, "LPFE_p11", :Optimal, 13, [2, 0, 1]))

  function problem_LPFE_p27(params::Params)
    return Lp.create_max_Cf0_problem(merge(params, {"maxit" => 5}),
      [10,-57,-9,-24], [0.5 -5.5 -2.5 9; 0.5 -1.5 -0.5 1; 1 0 0 0], [0, 0, 1]
      )
  end; push!(prob_db, DbProblem(problem_LPFE_p27, "LPFE_p27 degen", :Optimal, 1, [1, 0, 1, 0]))

  function problem_LPFE_m27(params::Params)
    return Lp.create_max_Cf0_problem(merge(params, {"maxit" => 5}),
      [10,-57,-9,-24], [0.5 -5.5 -2.5 9; 0.5 -1.5 -0.5 1; 1 0 0 0], [1.e-5, 0, 1]
      )
  end; push!(prob_db, DbProblem(problem_LPFE_m27, "LPFE_m27 pert", :Optimal, 1, [1, 0, 1, 0]))

  function problem_LPFE_p18(params::Params)
    return Lp.create_max_Cf0_problem(merge(params, {"maxit" => 5}),
      [-2,-1], [-1 1; -1 -2; 0 1], [-1, -2, 1]
      )
  end; push!(prob_db, DbProblem(problem_LPFE_p18, "LPFE_p18 phaseI", :Optimal, -3, [4/3, 1/3]))

  function problem_LPFE_p23_9(params::Params)
    return Lp.create_min_Cf0_problem(merge(params, {"maxit" => 10}),
      [2,3,4], [0 2 3; 1 1 2; 1 2 3], [5, 4, 7]
      )
  end; push!(prob_db, DbProblem(problem_LPFE_p23_9, "problem_LPFE_p23_9", :Optimal, 0, [0, 0, 0]))

  function problem_LPFE_p23_8(params::Params)
    return Lp.create_max_Cf0_problem(merge(params, {"maxit" => 10}),
      [3,2], [1 -2; 1 -1; 2 -1; 1 0; 2 1; 1 1; 1 2; 0 1], [1, 2, 6, 5, 16, 12, 21, 10]
      )
  end; push!(prob_db, DbProblem(problem_LPFE_p23_8, "problem_LPFE_p23_8", :Optimal, 28, [4, 8]))

  function problem_LPFE_p22_4(params::Params)
    return Lp.create_max_Cf0_problem(merge(params, {"maxit" => 10}),
      [-1,-3,-1], [2 -5 1; 2 -1 2], [-5, 4]
      )
  end; push!(prob_db, DbProblem(problem_LPFE_p22_4, "problem_LPFE_p22_4", :Optimal, -3, [0, 1, 0]))

  function random_problem_impl(form, seed, scale, dense, n, m, maxit, use_lohi, show, params::Params)
    rng = MersenneTwister(seed)

    c = scale * rand(rng, (n))
    b = scale * rand(rng, (m))
    A = scale * rand(rng, (m, n))
    dtype = typeof(1.0)
    lohi = Array(dtype, 0)
    has_lohi = false

    if (form == string(:Cf2b))
      has_lohi = true
      if (use_lohi)
        lohi = scale * rand(rng, 2*n)
        for i = 1:length(lohi):2
          if (lohi[i] > lohi[i+1])
            lohi[i], lohi[i+1] = lohi[i+1], lohi[i]
          end
        end
      else
        lohi = Array(dtype, 2*n)
        for i = 1:length(lohi):2
          lohi[i] = zero(dtype); lohi[i] = inf(dtype)
        end
      end
    end

    if dense != 100
      #zeros = rand(rng, 1:(m*n), zero_n) # Needs next version of Julia
      zero_n = int(round((m*n) * (1.0 - dense / 100.0)))
      if (zero_n > 0)
        #sprand is useless, the output will not be the same even with fixed seed.
        toindex = x -> clamp(int((m*n)* 0.5*(1.0+x)), 1, m*n)
        xis = [ toindex(x) for x in randn(rng, zero_n) ]
        for i=2:length(xis)
          it = 0
          while (it < 5 && xis[i] in xis[1:i-1])
            xis[i] = toindex(rand(rng))
            it = it + 1
          end
        end
        for i in xis
          A[i] = 0
        end
      end
    end

    if (show)
      println("c", c)
      println("b", b)
      println("A", A)
      if (has_lohi) println("lohi", lohi) end
    end

    nparams = merge(params, {"maxit" => maxit})
    if (form == string(:Cf0))
      return Lp.create_max_Cf0_problem(nparams, c, A, b)
    else
      return Lp.create_max_Cf2b_problem(nparams, c, A, b, lohi)
    end
  end

  function random_problem(params::Params)
    seed = Arg.dict_get(params, "seed", int(time_ns()) % 32768 )
    scale = int(Arg.dict_get(params, "scale", 1.0))
    dense = Arg.dict_get(params, "dense", 100 )
    n = Arg.dict_get(params, "n", 10)
    m = Arg.dict_get(params, "m", n)
    maxit = Arg.dict_get(params, "maxit", 2*(m+n))
    show = Arg.dict_get(params, "show", false)
    use_lohi = Arg.dict_get(params, "lohi", false)
    form = Arg.dict_get(params, "form", string(:Cf0))
    println("/seed:", seed)
    return random_problem_impl(form, seed, scale, dense, n, m, maxit, use_lohi, show, params)
  end; push!(prob_db, DbProblem(random_problem, "random", :Unset, :Unset, :Unset))

  function seed_problem_1(params::Params)
    return random_problem_impl(15578, 1, 60, 20, 20, 2*(20+20), false, params)
  end; push!(prob_db, DbProblem(seed_problem_1, "seed_1", :Unset, :Unset, :Unset))

  function last_problem(params::Params)
    if (length(Lp_bench.prob_last) > 0)
      return Lp_bench.prob_last[end]
    else
      return problem_LPFE_p88(params)
    end
  end; push!(prob_db, DbProblem(last_problem, "last", :Unset, :Unset, :Unset))

  #p184 (general)

  function check_sol(dbprob, lp_prob, sol, params)

    prefix = " "

    if (dbprob.check_status != :Unset)
      if (dbprob.check_status != sol.status)
        print_with_color(:red, "$(prefix)status: '$(sol.status)' should be '$(dbprob.check_status)' \n")
      else (dbprob.check_status != sol.status)
        print_with_color(:green, "$(prefix)status: '$(sol.status)' \n")
      end
    else
      println("$(prefix)$(sol.status)")
    end

    if (dbprob.check_z != :Unset)
      check_z = lp_prob.conv.f(dbprob.check_z)
      if (check_z != sol.z)
        print_with_color(:red, "$(prefix)z: '$(sol.z)' should be '$(check_z)' \n")
      else
        print_with_color(:green, "$(prefix)z: '$(sol.z)' \n")
      end
    else
      if (sol.solved)
        println("$(prefix)z: $(sol.z)")
      end
    end

    if (dbprob.check_x != :Unset)
      check_x = [ lp_prob.conv.f(x) for x in dbprob.check_x ]
      if (sol.solved == false || check_x != sol.x)
        sol_x = sol.solved ? sol.x : ()
        print_with_color(:red, "$(prefix)x: '$(sol_x)' should be '$(check_x)' \n")
      else
        print_with_color(:green, "$(prefix)x: '$(sol.x)' \n")
      end
    else
      if (sol.solved)
        println("$(prefix)x: $(sol.x)")
      else
        if (params["Dcd"] == true)
          println("$(prefix)nbasis:")
          Lp.dcd_nbasis(sol)
        end
      end
    end

    if (sol.iter != 0)
      println("$(prefix)it: $(sol.iter)")
    end

    println();
  end

  function introspect(arg_str::String = "", code_module::Module = lp_defaults)
    def_params = { "type"=>"Float32", "kind"=>"native" }
    params = deepcopy(def_params)
    arg_get(arg_create(arg_str), params)

    function do_introspect(kind, out, code_func, intr_func, intr_arg_types)
      str = ""
      if (out == :Str)
        str = join(code_func(intr_func, intr_arg_types))
      else
        orig_stdout = STDOUT
        try
          (outRead, outWrite) = redirect_stdout()
          code_func(intr_func, intr_arg_types)
          str = readavailable(outRead)
        catch exc
          redirect_stdout(orig_stdout)
          println(exc)
        end
        redirect_stdout(orig_stdout)
      end
      println("--------------: ", kind)
      println(str)
      nlines = length(matchall(r"\n", str))
      println("--------------> ", nlines, " lines")
    end

    code_info_map = {"typed" => (:Str, code_typed), "lowered" => (:Str, code_lowered), "llvm" => (:Stdout, code_llvm), "native" => (:Stdout, code_native) }
    kind = params["kind"]
    intr_func = code_module.solve_problem
    typestr = params["type"]
    solve_arg_types = ( eval(parse("Lp.Cf0_problem{$(typestr)}")), )
    intr_arg_types = solve_arg_types

    # refactor
    intr_func = code_module.solve_dat
    intr_arg_types = ( eval(parse("lp_I_1.Dat{$(typestr)}")), Lp.Solution, )

    if (kind != "all")
      code_info = code_info_map[params["kind"]]
      do_introspect(kind, code_info[1], code_info[2], intr_func, intr_arg_types)
    else
      for k in code_func_map
        code_info = code_info_map[k]
        do_introspect(code_info[1], code_info[2], intr_func, intr_arg_types)
      end
    end
  end

  function solve_problem(code_module::Module, params::Params, lp_prob)
    T = Lp.get_t(lp_prob)
    dat = code_module.construct_dat(T)
    code_module.fill_dat(lp_prob, dat)
    sol = Lp.construct_solution(T, params)
    code_module.solve_dat(dat, sol)
    return sol
  end

  function solve_problem(code_module::Symbol, params::Params, lp_prob::Lp.Cf0_problem)
    code_module = ( code_module == :Unset ? Jad.lp_I_1 : code_module )
    return solve_problem(code_module, params, lp_prob)
  end

  function solve_problem(code_module::Symbol, params::Params, lp_prob::Lp.Cf2b_problem)
    code_module = ( code_module == :Unset ? Jad.lp_III_1 : code_module )
    return solve_problem(code_module, params, lp_prob)
  end

  function solve(arg_str::String = "/prob:p88", code_module = :Unset)
    def_params = { "prob"=>"p88", "type"=>"Float32", "Dcd"=>false }
    params = deepcopy(def_params)
    arg_get(arg_create(arg_str), params)

    prob_key = params["prob"]
    println(prob_key)
    if (length(prob_key) == 0 || prob_key == "all")
      for i = 1:length(prob_db) solve(i, arg_str); end;
      return
    end

    found_prob = false; dbprob = 0;
    for i = 1:length(prob_db)
      iprob = prob_db[i]
      if (contains(iprob.descr, prob_key))
        found_prob = true; dbprob = iprob;
        break
      end
    end

    if (found_prob)
      println()
      print("Problem: '", dbprob.descr, "'")
      println("\n------------")
      lp_prob = dbprob.creator(params)
      #Profile.print()
      if (length(Lp_bench.prob_last) > 0) pop!(Lp_bench.prob_last) end
      push!(Lp_bench.prob_last, deepcopy(lp_prob))
      #Lp_bench.prob_last[1] =
      println("+++++++",
        " n:", Lp.get_n(lp_prob), ", m:", Lp.get_m(lp_prob),
        ", density:", Shared_funcs.format_percent(Lp.comp_density(lp_prob)),
        ", type:", Lp.get_t(lp_prob),
        ", form:", Lp.get_form(lp_prob),
        " +++++++")

      @time raw_sol = solve_problem(code_module, params, lp_prob)
      println("************\n")

      sol = Lp.translate_solution(lp_prob, raw_sol)
      check_sol(dbprob, lp_prob, sol, params)
      return sol
    end
    return 0
  end

end

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
using Base.Cartesian

import Base: show, ==, push!, append!, float
import LinearAlgebra: norm, normalize, normalize!

## Fast getindex function for multiple arrays, returns a tuple of array elements
@inline Base.@propagate_inbounds @generated function _multi_getindex(i::Union{Integer, CartesianIndex}, c::AbstractArray...)
    N = length(c)
    result_expr = Expr(:tuple)
    for j in 1:N
        push!(result_expr.args, :(c[$j][i]))
    end
    result_expr
end

# Need a generated function to promote edge types, because a simple
# promote_type(map(eltype, h.edges)...) isn't type stable (tested
# with Julia v0.5).
@generated function _promote_edge_types(edges::NTuple{N,AbstractVector}) where N
    promote_type(map(eltype, edges.parameters)...)
end

## nice-valued ranges for histograms
function histrange(v::AbstractArray{T}, n::Integer, closed::Symbol=:left) where T
    F = float(T)
    nv = length(v)
    if nv == 0 && n < 0
        throw(ArgumentError("number of bins must be ≥ 0 for an empty array, got $n"))
    elseif nv > 0 && n < 1
        throw(ArgumentError("number of bins must be ≥ 1 for a non-empty array, got $n"))
    elseif nv == 0
        return zero(F):zero(F)
    end

    lo, hi = extrema(v)
    histrange(F(lo), F(hi), n, closed)
end

function histrange(lo::F, hi::F, n::Integer, closed::Symbol=:left) where F
    if hi == lo
        start = F(hi)
        step = one(F)
        divisor = one(F)
        len = one(F)
    else
        bw = (F(hi) - F(lo)) / n
        lbw = log10(bw)
        if lbw >= 0
            step = exp10(floor(lbw))
            r = bw / step
            if r <= 1.1
                nothing
            elseif r <= 2.2
                step *= 2
            elseif r <= 5.5
                step *= 5
            else
                step *= 10
            end
            divisor = one(F)
            start = step*floor(lo/step)
            len = ceil((hi - start)/step)
        else
            divisor = exp10(-floor(lbw))
            r = bw * divisor
            if r <= 1.1
                nothing
            elseif r <= 2.2
                divisor /= 2
            elseif r <= 5.5
                divisor /= 5
            else
                divisor /= 10
            end
            step = one(F)
            start = floor(lo*divisor)
            len = ceil(hi*divisor - start)
        end
    end
    # fix up endpoints
    if closed == :right #(,]
        while lo <= start/divisor
            start -= step
        end
        while (start + (len-1)*step)/divisor < hi
            len += one(F)
        end
    else
        while lo < start/divisor
            start -= step
        end
        while (start + (len-1)*step)/divisor <= hi
            len += one(F)
        end
    end
    Base.floatrange(start,step,len,divisor)
end

histrange(vs::NTuple{N,AbstractVector},nbins::NTuple{N,Integer},closed::Symbol) where {N} =
    map((v,n) -> histrange(v,n,closed),vs,nbins)
histrange(vs::NTuple{N,AbstractVector},nbins::Integer,closed::Symbol) where {N} =
    map(v -> histrange(v,nbins,closed),vs)

## histograms ##
function sturges(n)  # Sturges' formula
    n==0 && return one(n)
    ceil(Integer, log2(n))+1
end

abstract type AbstractHistogram{T<:Real,N,E} end

# N-dimensional histogram object
"""
    Histogram <: AbstractHistogram

The `Histogram` type represents data that has been tabulated into intervals
(known as *bins*) along the real line, or in higher dimensions, over a real space.
Histograms can be fitted to data using the `fit` method.

# Fields
* edges: An iterator that contains the boundaries of the bins in each dimension.
* weights: An array that contains the weight of each bin.
* closed: A symbol with value `:right` or `:left` indicating on which side bins
  (half-open intervals or higher-dimensional analogues thereof) are closed.
  See below for an example.
* isdensity: There are two interpretations of a `Histogram`. If `isdensity=false` the weight of a bin corresponds to the amount of a quantity in the bin.
  If `isdensity=true` then it corresponds to the density (amount / volume) of the quantity in the bin. See below for an example.

# Examples
## Example illustrating `closed`
```jldoctest
julia> using StatsBase

julia> fit(Histogram, [2.],  1:3, closed=:left)
Histogram{Int64,1,Tuple{UnitRange{Int64}}}
edges:
  1:3
weights: [0, 1]
closed: left
isdensity: false

julia> fit(Histogram, [2.],  1:3, closed=:right)
Histogram{Int64,1,Tuple{UnitRange{Int64}}}
edges:
  1:3
weights: [1, 0]
closed: right
isdensity: false
```
## Example illustrating `isdensity`
```julia
julia> using StatsBase, LinearAlgebra

julia> bins = [0,1,7]; # a small and a large bin

julia> obs = [0.5, 1.5, 1.5, 2.5]; # one observation in the small bin and three in the large

julia> h = fit(Histogram, obs, bins)
Histogram{Int64,1,Tuple{Array{Int64,1}}}
edges:
  [0, 1, 7]
weights: [1, 3]
closed: left
isdensity: false

julia> # observe isdensity = false and the weights field records the number of observations in each bin

julia> normalize(h, mode=:density)
Histogram{Float64,1,Tuple{Array{Int64,1}}}
edges:
  [0, 1, 7]
weights: [1.0, 0.5]
closed: left
isdensity: true

julia> # observe isdensity = true and weights tells us the number of observation per binsize in each bin
```
"""
mutable struct Histogram{T<:Real,N,E} <: AbstractHistogram{T,N,E}
    edges::E
    weights::Array{T,N}
    closed::Symbol
    isdensity::Bool
    function Histogram{T,N,E}(edges::NTuple{N,AbstractArray}, weights::Array{T,N},
                              closed::Symbol, isdensity::Bool=false) where {T,N,E}
        closed == :right || closed == :left || error("closed must :left or :right")
        isdensity && !(T <: AbstractFloat) && error("Density histogram must have float-type weights")
        _edges_nbins(edges) == size(weights) || error("Histogram edge vectors must be 1 longer than corresponding weight dimensions")
        new{T,N,E}(edges,weights,closed,isdensity)
    end
end

Histogram(edges::NTuple{N,AbstractVector}, weights::AbstractArray{T,N},
          closed::Symbol=:left, isdensity::Bool=false) where {T,N} =
    Histogram{T,N,typeof(edges)}(edges,weights,closed,isdensity)

Histogram(edges::NTuple{N,AbstractVector}, ::Type{T}, closed::Symbol=:left,
          isdensity::Bool=false) where {T,N} =
    Histogram(edges,zeros(T,_edges_nbins(edges)...),closed,isdensity)

Histogram(edges::NTuple{N,AbstractVector}, closed::Symbol=:left,
          isdensity::Bool=false) where {N} =
    Histogram(edges,Int,closed,isdensity)

function show(io::IO, h::AbstractHistogram)
    println(io, typeof(h))
    println(io,"edges:")
    for e in h.edges
        println(io,"  ",e)
    end
    println(io,"weights: ",h.weights)
    println(io,"closed: ",h.closed)
    print(io,"isdensity: ",h.isdensity)
end

(==)(h1::Histogram,h2::Histogram) = (==)(h1.edges,h2.edges) && (==)(h1.weights,h2.weights) && (==)(h1.closed,h2.closed) && (==)(h1.isdensity,h2.isdensity)

binindex(h::AbstractHistogram{T,1}, x::Real) where {T} = binindex(h, (x,))[1]

binindex(h::Histogram{T,N}, xs::NTuple{N,Real}) where {T,N} =
    map((edge, x) -> _edge_binindex(edge, h.closed, x), h.edges, xs)

@inline function _edge_binindex(edge::AbstractVector, closed::Symbol, x::Real)
    if closed == :right
        searchsortedfirst(edge, x) - 1
    else
        searchsortedlast(edge, x)
    end
end

binvolume(h::AbstractHistogram{T,1}, binidx::Integer) where {T} = binvolume(h, (binidx,))
binvolume(::Type{V}, h::AbstractHistogram{T,1}, binidx::Integer) where {V,T} = binvolume(V, h, (binidx,))

binvolume(h::Histogram{T,N}, binidx::NTuple{N,Integer}) where {T,N} =
    binvolume(_promote_edge_types(h.edges), h, binidx)

binvolume(::Type{V}, h::Histogram{T,N}, binidx::NTuple{N,Integer}) where {V,T,N} =
    prod(map((edge, i) -> _edge_binvolume(V, edge, i), h.edges, binidx))

@inline _edge_binvolume(::Type{V}, edge::AbstractVector, i::Integer) where {V} = V(edge[i+1]) - V(edge[i])
@inline _edge_binvolume(::Type{V}, edge::AbstractRange, i::Integer) where {V} = V(step(edge))
@inline _edge_binvolume(edge::AbstractVector, i::Integer) = _edge_binvolume(eltype(edge), edge, i)

@inline _edges_nbins(edges::NTuple{N,AbstractVector}) where {N} = map(_edge_nbins, edges)

@inline _edge_nbins(edge::AbstractVector) = length(edge) - 1

# 1-dimensional

Histogram(edge::AbstractVector, weights::AbstractVector{T}, closed::Symbol=:left, isdensity::Bool=false) where {T} =
    Histogram((edge,), weights, closed, isdensity)

Histogram(edge::AbstractVector, ::Type{T}, closed::Symbol=:left, isdensity::Bool=false) where {T} =
    Histogram((edge,), T, closed, isdensity)

Histogram(edge::AbstractVector, closed::Symbol=:left, isdensity::Bool=false) =
    Histogram((edge,), closed, isdensity)

push!(h::AbstractHistogram{T,1}, x::Real, w::Real) where {T} = push!(h, (x,), w)
push!(h::AbstractHistogram{T,1}, x::Real) where {T} = push!(h,x,one(T))
append!(h::AbstractHistogram{T,1}, v::AbstractVector) where {T} = append!(h, (v,))
append!(h::AbstractHistogram{T,1}, v::AbstractVector, wv::Union{AbstractVector,AbstractWeights}) where {T} = append!(h, (v,), wv)

fit(::Type{Histogram{T}},v::AbstractVector, edg::AbstractVector; closed::Symbol=:left) where {T} =
    fit(Histogram{T},(v,), (edg,), closed=closed)
fit(::Type{Histogram{T}},v::AbstractVector; closed::Symbol=:left, nbins=sturges(length(v))) where {T} =
    fit(Histogram{T},(v,); closed=closed, nbins=nbins)
fit(::Type{Histogram{T}},v::AbstractVector, wv::AbstractWeights, edg::AbstractVector; closed::Symbol=:left) where {T} =
    fit(Histogram{T},(v,), wv, (edg,), closed=closed)
fit(::Type{Histogram{T}},v::AbstractVector, wv::AbstractWeights; closed::Symbol=:left, nbins=sturges(length(v))) where {T} =
    fit(Histogram{T}, (v,), wv; closed=closed, nbins=nbins)

fit(::Type{Histogram}, v::AbstractVector, wv::AbstractWeights{W}, args...; kwargs...) where {W} = fit(Histogram{W}, v, wv, args...; kwargs...)

# N-dimensional

function push!(h::Histogram{T,N},xs::NTuple{N,Real},w::Real) where {T,N}
    h.isdensity && error("Density histogram must have float-type weights")
    idx = binindex(h, xs)
    if checkbounds(Bool, h.weights, idx...)
        @inbounds h.weights[idx...] += w
    end
    h
end

function push!(h::Histogram{T,N},xs::NTuple{N,Real},w::Real) where {T<:AbstractFloat,N}
    idx = binindex(h, xs)
    if checkbounds(Bool, h.weights, idx...)
        @inbounds h.weights[idx...] += h.isdensity ? w / binvolume(h, idx) : w
    end
    h
end

push!(h::AbstractHistogram{T,N},xs::NTuple{N,Real}) where {T,N} = push!(h,xs,one(T))

function append!(h::AbstractHistogram{T,N}, vs::NTuple{N,AbstractVector}) where {T,N}
    @inbounds for i in eachindex(vs...)
        xs = _multi_getindex(i, vs...)
        push!(h, xs, one(T))
    end
    h
end
function append!(h::AbstractHistogram{T,N}, vs::NTuple{N,AbstractVector}, wv::AbstractVector) where {T,N}
    @inbounds for i in eachindex(wv, vs...)
        xs = _multi_getindex(i, vs...)
        push!(h, xs, wv[i])
    end
    h
end

# Turn kwargs nbins into a type-stable tuple of integers:
function _nbins_tuple(vs::NTuple{N,AbstractVector}, nbins) where N
    template = map(length, vs)
    result = broadcast((t, x) -> typeof(t)(x), template, nbins)
    result::typeof(template)
end

fit(::Type{Histogram{T}}, vs::NTuple{N,AbstractVector}, edges::NTuple{N,AbstractVector}; closed::Symbol=:left) where {T,N} =
    append!(Histogram(edges, T, closed, false), vs)

fit(::Type{Histogram{T}}, vs::NTuple{N,AbstractVector}; closed::Symbol=:left, nbins=sturges(length(vs[1]))) where {T,N} =
    fit(Histogram{T}, vs, histrange(vs,_nbins_tuple(vs, nbins),closed); closed=closed)

fit(::Type{Histogram{T}}, vs::NTuple{N,AbstractVector}, wv::AbstractWeights{W}, edges::NTuple{N,AbstractVector}; closed::Symbol=:left) where {T,N,W} =
    append!(Histogram(edges, T, closed, false), vs, wv)

fit(::Type{Histogram{T}}, vs::NTuple{N,AbstractVector}, wv::AbstractWeights; closed::Symbol=:left, nbins=sturges(length(vs[1]))) where {T,N} =
    fit(Histogram{T}, vs, wv, histrange(vs,_nbins_tuple(vs, nbins),closed); closed=closed)

"""
    fit(Histogram, data[, weight][, edges]; closed=:left, nbins)

Fit a histogram to `data`.

# Arguments

* `data`: either a vector (for a 1-dimensional histogram), or a tuple of
  vectors of equal length (for an *n*-dimensional histogram).

* `weight`: an optional `AbstractWeights` (of the same length as the
  data vectors), denoting the weight each observation contributes to the
  bin. If no weight vector is supplied, each observation has weight 1.

* `edges`: a vector (typically an `AbstractRange` object), or tuple of vectors, that gives
  the edges of the bins along each dimension. If no edges are provided, these
  are determined from the data.

# Keyword arguments

* `closed`: if `:left` (the default), the bin intervals are left-closed [a,b);
  if `:right`, intervals are right-closed (a,b].

* `nbins`: if no `edges` argument is supplied, the approximate number of bins to use
  along each dimension (can be either a single integer, or a tuple of integers).

# Examples

```julia
# Univariate
h = fit(Histogram, rand(100))
h = fit(Histogram, rand(100), 0:0.1:1.0)
h = fit(Histogram, rand(100), nbins=10)
h = fit(Histogram, rand(100), weights(rand(100)), 0:0.1:1.0)
h = fit(Histogram, [20], 0:20:100)
h = fit(Histogram, [20], 0:20:100, closed=:right)

# Multivariate
h = fit(Histogram, (rand(100),rand(100)))
h = fit(Histogram, (rand(100),rand(100)),nbins=10)
```
"""
fit(::Type{Histogram}, args...; kwargs...) = fit(Histogram{Int}, args...; kwargs...)
fit(::Type{Histogram}, vs::NTuple{N,AbstractVector}, wv::AbstractWeights{W}, args...; kwargs...) where {N,W} = fit(Histogram{W}, vs, wv, args...; kwargs...)

# Get a suitable high-precision type for the norm of a histogram.
norm_type(h::Histogram{T,N}) where {T,N} =
    promote_type(T, _promote_edge_types(h.edges))

norm_type(::Type{T}) where {T<:Integer} = promote_type(T, Int64)
norm_type(::Type{T}) where {T<:AbstractFloat} = promote_type(T, Float64)

"""
    norm(h::Histogram)

Calculate the norm of histogram `h` as the absolute value of its integral.
"""
@generated function norm(h::Histogram{T,N}) where {T,N}
    quote
        edges = h.edges
        weights = h.weights
        SumT = norm_type(h)
        v_0 = 1
        s_0 = zero(SumT)
        @inbounds @nloops(
            $N, i, weights,
            d -> begin
                v_{$N-d+1} = v_{$N-d} * _edge_binvolume(SumT, edges[d], i_d)
                s_{$N-d+1} = zero(SumT)
            end,
            d -> begin
                s_{$N-d} += s_{$N-d+1}
            end,
            begin
                $(Symbol("s_$(N)")) += (@nref $N weights i) * $(Symbol("v_$N"))
            end
        )
        s_0
    end
end

float(h::Histogram{T,N}) where {T<:AbstractFloat,N} = h

float(h::Histogram{T,N}) where {T,N} = Histogram(h.edges, float(h.weights), h.closed, h.isdensity)

"""
    normalize!(h::Histogram{T,N}, aux_weights::Array{T,N}...; mode::Symbol=:pdf) where {T<:AbstractFloat,N}

Normalize the histogram `h` and optionally scale one or more auxiliary weight
arrays appropriately. See description of `normalize` for details. Returns `h`.
"""
@generated function normalize!(h::Histogram{T,N}, aux_weights::Array{T,N}...; mode::Symbol=:pdf) where {T<:AbstractFloat,N}
    quote
        edges = h.edges
        weights = h.weights

        for A in aux_weights
            (size(A) != size(weights)) && throw(DimensionMismatch("aux_weights must have same size as histogram weights"))
        end

        if mode == :none
            # nothing to do
        elseif mode == :pdf || mode == :density || mode == :probability
            if h.isdensity
                if mode == :pdf || mode == :probability
                    # histogram already represents a density, just divide weights by norm
                    s = 1/norm(h)
                    weights .*= s
                    for A in aux_weights
                        A .*= s
                    end
                else
                    # :density - histogram already represents a density, nothing to do
                end
            else
                if mode == :pdf || mode == :density
                    # Divide weights by bin volume, for :pdf also divide by sum of weights
                    SumT = norm_type(h)
                    vs_0 = (mode == :pdf) ? sum(SumT(x) for x in weights) : one(SumT)
                    @inbounds @nloops $N i weights d->(vs_{$N-d+1} = vs_{$N-d} * _edge_binvolume(SumT, edges[d], i_d)) begin
                        (@nref $N weights i) /= $(Symbol("vs_$N"))
                        for A in aux_weights
                            (@nref $N A i) /= $(Symbol("vs_$N"))
                        end
                    end
                    h.isdensity = true
                else
                    # :probability - divide weights by sum of weights
                    nf = inv(sum(weights))
                    weights .*= nf
                    for A in aux_weights
                        A .*= nf
                    end
                end
            end
        else
            throw(ArgumentError("Normalization mode must be :pdf, :density, :probability or :none"))
        end
        h
    end
end

"""
    normalize(h::Histogram{T,N}; mode::Symbol=:pdf) where {T,N}

Normalize the histogram `h`.

Valid values for `mode` are:

*  `:pdf`: Normalize by sum of weights and bin sizes. Resulting histogram
   has norm 1 and represents a PDF.
* `:density`: Normalize by bin sizes only. Resulting histogram represents
   count density of input and does not have norm 1. Will not modify the
   histogram if it already represents a density (`h.isdensity == 1`).
* `:probability`: Normalize by sum of weights only. Resulting histogram
   represents the fraction of probability mass for each bin and does not have
   norm 1.
*  `:none`: Leaves histogram unchanged. Useful to simplify code that has to
   conditionally apply different modes of normalization.

Successive application of both `:probability` and `:density` normalization (in
any order) is equivalent to `:pdf` normalization.
"""
normalize(h::Histogram{T,N}; mode::Symbol=:pdf) where {T,N} =
    normalize!(deepcopy(float(h)), mode = mode)

"""
    normalize(h::Histogram{T,N}, aux_weights::Array{T,N}...; mode::Symbol=:pdf) where {T,N}

Normalize the histogram `h` and rescales one or more auxiliary weight arrays
at the same time (`aux_weights` may, e.g., contain estimated statistical
uncertainties). The values of the auxiliary arrays are scaled by the same
factor as the corresponding histogram weight values. Returns a tuple of the
normalized histogram and scaled auxiliary weights.
"""
function normalize(h::Histogram{T,N}, aux_weights::Array{T,N}...; mode::Symbol=:pdf) where {T,N}
    h_fltcp = deepcopy(float(h))
    aux_weights_fltcp = map(x -> deepcopy(float(x)), aux_weights)
    normalize!(h_fltcp, aux_weights_fltcp..., mode = mode)
    (h_fltcp, aux_weights_fltcp...)
end

"""
    zero(h::Histogram)

Create a new histogram with the same binning, type and shape of weights
and the same properties (`closed` and `isdensity`) as `h`, with all weights
set to zero.
"""
Base.zero(h::Histogram{T,N,E}) where {T,N,E} =
    Histogram{T,N,E}(deepcopy(h.edges), zero(h.weights), h.closed, h.isdensity)

"""
    merge!(target::Histogram, others::Histogram...)

Update histogram `target` by merging it with the histograms `others`. See
`merge(histogram::Histogram, others::Histogram...)` for details.
"""
function Base.merge!(target::Histogram, others::Histogram...)
    for h in others
        target.edges != h.edges && throw(ArgumentError("can't merge histograms with different binning"))
        size(target.weights) != size(h.weights) && throw(ArgumentError("can't merge histograms with different dimensions"))
        target.closed != h.closed && throw(ArgumentError("can't merge histograms with different closed left/right settings"))
        target.isdensity != h.isdensity && throw(ArgumentError("can't merge histograms with different isdensity settings"))
    end
    for h in others
        target.weights .+= h.weights
    end
    target
end

"""
    merge(h::Histogram, others::Histogram...)

Construct a new histogram by merging `h` with `others`. All histograms must
have the same binning, shape of weights and properties (`closed` and
`isdensity`). The weights of all histograms are summed up for each bin, the
weights of the resulting histogram will have the same type as those of `h`.
"""
Base.merge(h::Histogram, others::Histogram...) = merge!(zero(h), h, others...)

"""
    StatsBase.midpoints(v)

Calculate the midpoints (pairwise mean of consecutive elements).
"""
midpoints(v::AbstractVector) = [middle(v[i - 1], v[i]) for i in 2:length(v)]

midpoints(r::AbstractRange) = r[1:(end - 1)] .+ (step(r) / 2)


text:
using TableTransforms
using Distributions
using Tables
using TypedTables
using CategoricalArrays
using ScientificTypes: Count, Multiclass
using LinearAlgebra
using Statistics
using Test, Random, Plots
using ReferenceTests, ImageIO

# set default configurations for plots
gr(ms=2, mc=:black, aspectratio=:equal,
   label=false, size=(600,400))

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
isCI = "CI" ∈ keys(ENV)
islinux = Sys.islinux()
visualtests = !isCI || (isCI && islinux)
datadir = joinpath(@__DIR__,"data")

# for functor tests in Functional testset
struct Polynomial{T<:Real}
  coeffs::Vector{T}
end
Polynomial(args::T...) where {T<:Real} = Polynomial(collect(args))
(p::Polynomial)(x) = sum(a * x^(i-1) for (i, a) in enumerate(p.coeffs))

function isequalmissing(a, b)
  length(a) == length(b) || return false
  for (x, y) in zip(a, b)
    x === y || return false
  end
  return true
end

# list of tests
testfiles = [
  "distributions.jl",
  "colspec.jl",
  "transforms.jl"
]

@testset "TableTransforms.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
module PuzzleGame

using Images, ImageView, TestImages
using Gtk, Reactive, GtkReactive
using Random

export main

function rules()
    println("please swap tiles to complete the image")
end

function init_game(img::Array{T,2}; nof_row_blocks::Int=4, nof_col_blocks::Int=4 ) where T  <: Colorant
    block_size = Int.(size(img) ./ (nof_row_blocks, nof_col_blocks))
    row_slices = [ i*block_size[1] + 1 : (i+1) * block_size[1] for i in 0:nof_row_blocks-1]
    col_slices = [ i*block_size[2] + 1 : (i+1) * block_size[2] for i in 0:nof_col_blocks-1]

    blocks = [(x,y) for y in col_slices for x in row_slices]

    shuf = zeros(eltype(img), size(img))
    shuf_blocks = shuffle(blocks)

    for b in 1:length(blocks)
        shuf[shuf_blocks[b][1], shuf_blocks[b][2]] = img[blocks[b][1], blocks[b][2]]
    end
    shuf
end

function main()
    img = testimage("fabio")
    new = init_game(img)
    guidict = imshow(new)
    rules()

    signal_connect(guidict["gui"]["window"], :destroy) do widget
        isAcitve = false
        println("you closded the window")
    end

    coord_sig = map(guidict["gui"]["canvas"].mouse.buttonpress) do btn
        row = Float64(btn.position.x)
        col = Float64(btn.position.y)
        println("you clicked @ [x=$(row), y=$(col)] with $(btn.button)")

        if firstClick
            println("highligth the selected block")
        else
            println("swap the blocks")
        end

        return row,col
    end

    return guidict
end

#Top level code
firstClick = false
isActive = true
if !isinteractive()
    main()
    while isAcitve
        sleep(1)
    end
    println("code done..")
end

end # module PuzzleGame


text:
function TiffFile(uuid::UUID, filepath::String)
    try
        io = open(filepath)
        file = read(io, TiffFile)
        file.uuid = uuid
        file.filepath = filepath
        return file
    catch e # the file probably got renamed
        # if isa(e, SystemError)
        #     dir = dirname(filepath)
        #     filter(x->occursin(".ome.tif", x), readdir(dir))
        # end
        throw(FileIO.LoaderError("OMETIFF", "It looks like this file was renamed, "*
        "but has internal links with the original name. Please rename to $filepath "*
        "to load. See https://github.com/tlnagy/OMETIFF.jl/issues/14 for details.", e))
    end
end

"""
    get_ifds(orig_file, ifd_index, ifd_files) -> Dict, Dict

Run through all the IFDs extracted from the OMEXML and open all the referenced
files to construct a mapping of ZCTP (not guaranteed order) index to IFD object.
This is necessary because there can be multiple files referenced in a single
OMEXML and we need to iterate over the files to identify the actual offsets for
the data since this information isn't found in the OMEXML.

**Output**
- `Dict{Tuple{String, String}, TiffFile}`: a mapping of filepath, UUID to the
   actual TiffFile object
- `OrderedDict{NTuple{4, Int}, IFD}`: a mapping of ZCTP (or other order) index
   to the IFD objects in order that the IFDs are referenced in the OMEXML
"""
function get_ifds(orig_file::TiffFile{O},
                  ifd_index::OrderedDict{Int, NTuple{4, Int}},
                  ifd_files::OrderedDict{Int, Tuple{UUID, String}}) where {O <: Unsigned}

    # open all files referenced by the tiffdatas
    files = OrderedDict{Tuple{UUID, String}, TiffFile{O}}()
    for item in unique(values(ifd_files))
        uuid, filepath = item

        target_uuid = orig_file.uuid

        # if this file is the same as the base file, then don't open it again
        if uuid !== nothing && uuid == target_uuid
            files[item] = orig_file
        # if the UUID is missing or different, attempt to match on filepath
        elseif filepath == orig_file.filepath
            files[item] = orig_file
        # else open a new pointer to this file
        else
            files[item] = TiffFile(uuid, filepath)
        end
    end

    ifds = OrderedDict{NTuple{4, Int}, Tuple{TiffFile{O}, IFD{O}}}()

    for (fileid, file) in files
        # load all IFDs in this file
        rawifds = collect(file)
        foreach(ifd -> load!(file, ifd), rawifds)

        # there is a global IFD index across all files and also a per-file
        # index, e.g. if each IFD is in a separate file than an IFD can have an
        # index of 5 globally yet have an index of 1 in the 5th file.
        ifd_idx_in_file = 0
        for ifd in sort(collect(keys(ifd_index)))
            ifd_file = ifd_files[ifd]

            if ifd_file == fileid
                ifds[ifd_index[ifd]] = (file, rawifds[ifd_idx_in_file+=1])
            end
        end
    end
    files, ifds
end

function loadxml(file::TiffFile)
    ifd = first(file)
    load!(file, ifd)

    xml = nothing
    # grab all IMAGEDESCRIPTION tags in the first ifd, there can be multiple so
    # we need to check all of them
    for descriptions in ifd[Iterable(IMAGEDESCRIPTION)]
        try
            newxml = parsexml(descriptions.data)
            if xml === nothing
                xml = newxml
            else
                throw(ErrorException("Multiple XML entries detected, aborting"))
            end
        catch
        end
    end
    root(xml)
end

make_uuid(uuid::String) = UUID(startswith(uuid, "urn:uuid:") ? uuid[10:end] : uuid)

"""
    load_master_xml(file::TiffFile) -> EzXML.doc

Loads the master OME-XML file from `file` or from a linked file.
"""
function load_master_xml(file::TiffFile)
    omexml = loadxml(file)
    # update the UUID of the file with the value from the OMEXML
    uuid_node = findfirst("/ns:OME/@UUID", omexml, ["ns"=>namespace(omexml)])
    if uuid_node !== nothing
        file.uuid = make_uuid(nodecontent(uuid_node))
    # if the file doesn't have a UUID attempt to match based on name attribute
    # of the Image element. This is often, though not guaranteed, where the
    # file's internal filename is stored. We need this to handle cases where the
    # file was renamed.
    else
        main_filename = nodecontent(findfirst("/ns:OME/ns:Image/@Name", omexml, ["ns"=>namespace(omexml)]))
        main_filename = basename(split(main_filename, ".ome.")[1])

        uuids = findall("/ns:OME//ns:TiffData/ns:UUID", omexml, ["ns"=>namespace(omexml)])

        plane_uuids = unique(nodecontent.(uuids))
        # if there's only one UUID mentioned in the XML then we'll assume it's
        # referring to the current file
        if length(plane_uuids) == 1
            uuid_node = first(uuids)
            if haskey(uuid_node, "FileName")
                plane_filename = basename(split(uuid_node["FileName"], ".ome.")[1])
                if main_filename != plane_filename
                    @warn "TIFF planes refer to a nonexistent file. Proceeding "*
                    "with assumption that they refer to the current file."
                end
            end

            file.uuid = make_uuid(nodecontent(uuid_node))
        end
    end
    try
        # Check if the full OME-XML metadata is stored in another file
        metadata_file = findfirst("/ns:OME/ns:BinaryOnly", omexml, ["ns"=>namespace(omexml)])
        (metadata_file === nothing) && return omexml
        fileuuid, filepath = metadata_file["UUID"], joinpath(dirname(file.filepath), metadata_file["MetadataFile"])

        # we have a companion metadata file
        if endswith(filepath, ".companion.ome")
            xdoc = readxml(filepath)
            omexml = root(xdoc)
        else
            io = open(filepath)
            metadata_file = read(io, TiffFile)
            metadata_file.uuid = make_uuid(fileuuid)
            omexml = loadxml(metadata_file)
            close(metadata_file.io) # clean up
            return omexml
        end
    catch err
        isa(err, BoundsError) && return omexml
        rethrow(err)
    end
end

"""
    dump_omexml(filepath) -> String

Returns the OME-XML embedded inside the OME-TIFF as a prettified string.
"""
function dump_omexml(filepath::String)
    if !endswith(filepath, ".ome.tif")
        error("Passed file is not an OME-TIFF")
    end
    io = IOBuffer()
    open(filepath) do f
        orig_file = read(f, TiffFile)
        omexml = load_master_xml(orig_file)
        prettyprint(io, omexml)
    end
    String(take!(io))
end

"""
    load_comments(file) -> String

Extracts the MicroManager embedded description, if present. Else returns an
empty string.
"""
function load_comments(file)
    seek(file, 24)
    comment_header = Int(read(file, UInt32))
    if comment_header != 99384722
        return ""
    end
    comment_offset = Int(read(file, UInt32))
    seek(file, comment_offset)
    comment_header = read(file, UInt32)
    comment_length = Int(read(file, UInt32))
    metadata_bytes = Array{UInt8}(undef, comment_length)
    read!(file, metadata_bytes)
    metadata = JSON.parse(String(metadata_bytes))
    if !haskey(metadata, "Summary")
        return ""
    end
    metadata["Summary"]
end


text:
using PtFEM

ProjDir = dirname(@__FILE__)
exname = split(@__FILE__)[end][1:end-3]

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 3, 2, 2, 16, :x, Quadrilateral(4, 4), false),
  :properties => [10.92 0.3;],
  :x_coords => collect(range(0.0, stop=0.5, length=3)),
  :y_coords => collect(range(0.0, stop=0.5, length=3)),
  :thickness => 1.0,
  :support => [
    (1, [0 0 0 1]),
    (2, [0 0 1 1]),
    (3, [0 0 1 0]),
    (4, [0 1 0 1]),
    (6, [1 0 1 0]),
    (7, [0 1 0 0]),
    (8, [1 1 0 0]),
    (9, [1 0 0 0])
    ],
  :loaded_nodes => [
    (9, [0.25 0.0 0.0 0.0])
    ]
)

data |> display
println()

@time dis_df, sigma_df = p47(data)
println()

#PtFEM.vtk(data, dis_df, sigma_df, ProjDir, "Ex47.1")
#println()

dis_df |> display
println()
println("\nThe integration point (nip = $(data[:struc_el].nip)) stresses are:")
sigma_df |> display


text:
using Weave
using Base.Test

# Running Changin plotlib in tests segfault, unless they are run
# using separate processes.
#run(`julia --code-coverage=user -e 'include("winston_formats.jl")'`)
#run(`julia --code-coverage=user -e 'include("pyplot_formats.jl")'`)

info("Test: Chunk options")
include("chunk_options.jl")

info("Test: Caching")
include("cache_test.jl")

info("Testing rich output")
include("rich_output.jl")

if VERSION < v"0.5-dev"
  info("Test: Chunk options with Gadfly")
  include("chunk_opts_gadfly.jl")

  #Fails on travis, works locally.
  info("Test: Weaving with Winston")
  include("winston_formats.jl")

  info("Test: Weaving with Gadfly")
  include("gadfly_formats.jl")

  info("Test: Weaving with PyPlot")
  include("pyplot_formats.jl")

  info("Test: Weaving with Plots.jl")
  include("plotsjl_test.jl")
end


text:
function sieve(limit)
    candidates = 2:limit
    out = Int[]
    while length(candidates) > 0
        x = candidates[1]
        push!(out, x)
        candidates = filter(c -> (c .% x) .!= 0, candidates)
    end
    return out
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
function gen_iter(n, d)
    # based on https://twitter.com/evalparse/status/1107964924024635392
    iter = []
    for x in partitions(d + n, n)
        x = x .- 1
        if all(x .<= d)
            ys = Set(multiset_permutations(x, n))
            for y in ys
                push!(iter, Tuple(y))
            end
        end
    end
    iter
end

function construct_iter_all(N::Int, order::Int)

    #Construct an ordered iterator going over all moments
    # sequentially in terms of order

    iters = Tuple{repeat([Int], N)...}[]
    for d in 0:order
        x = Base.sort(gen_iter(N, d), rev=true)
        push!(iters, x...)
    end

    iters

end

# Trim a string of form "(a, b, c, d, ...)" to "abcd..."
trim_key(expr) = filter(x -> !(isspace(x) || x == ')' || x== '(' || x==','), string(expr))

# Expand a symbolic expression (no binomial expansion)
expansion_rule_mod = @acrule ~x * +(~~ys) => sum(map(y-> ~x * y, ~~ys))
expand_mod = Fixpoint(Prewalk(PassThrough(expansion_rule_mod)))
flatten_rule_mod = @rule(~x::isnotflat(+) => flatten_term(+, ~x))
flatten_mod = Fixpoint(PassThrough(flatten_rule_mod))
expand_expr = Fixpoint(PassThrough(Chain([expand_mod, flatten_mod])))

function define_μ(N::Int, order::Int, iter=construct_iter_all(N, order))

    indices = String[]
    for i in iter
        push!(indices, MomentClosure.trim_key(i))
    end

    @parameters t
    t = [t.val]

    μs = OrderedDict()
    for (i, idx) in enumerate(iter)
        if sum(idx) == 0
            μs[idx] = 1
        else
            μs[idx] = Term{Real}(Sym{FnType{Tuple{Any}, Real}}(Symbol("μ$(join(map_subscripts(indices[i]), ""))")), t)
        end
    end

    μs

end

function define_M(N::Int, order::Int, iter=construct_iter_all(N, order))

    indices = String[]
    for i in iter
        push!(indices, MomentClosure.trim_key(i))
    end

    @parameters t
    t = [t.val]

    Ms = OrderedDict()
    for (i, idx) in enumerate(iter)
        if sum(idx) == 0
            Ms[idx] = 1
        elseif sum(idx) == 1
            Ms[idx] = 0
        else
            Ms[idx] = Term{Real}(Sym{FnType{Tuple{Any}, Real}}(Symbol("M$(join(map_subscripts(indices[i]), ""))")), t)
        end
    end

    Ms

end

function extract_variables(eqs::Array{Equation, 1}, N::Int, q_order::Int)

    iters = construct_iter_all(N, q_order)
    iter_μ = filter(x -> sum(x) > 0, iters)
    iter_M = filter(x -> sum(x) > 1, iters)

    μs = values(define_μ(N, q_order, iter_μ))
    Ms = values(define_M(N, q_order, iter_M))
    vars = vcat(μs..., Ms...)
    # extract variables from rhs of each equation
    eq_vars = unique(vcat(get_variables.(eqs)...))
    # need this as get_variables does not extract var from `Differential(t)(var(t))`
    diff_vars = [var_from_nested_derivative(eq.lhs)[1] for eq in eqs]
    # filter out the unique ones
    eq_vars = unique(vcat(eq_vars..., diff_vars...))
    # this should preserve the correct ordering

    vars = intersect!(vars, eq_vars)
    vars

end

## Set of functions to deconstruct polynomial propensities ##

#=
    Consider a polynomial propensity function: a₁ = x^2 + c₁*y²*x,
    where x(t) and y(t) are the molecule numbers variables and c₁ is reaction parameter
    First we need to split a₁ into separate terms, i.e. x^2 and c₁*y²*x. Then we
    determine the independent multiplication factor in each term (1 and c₁). Finally,
    we obtain the power each variable is raised to in each term (x to power 2 in term 1;
    x to power 1 and y to power 2 in term 2). Having this information we can proceed
    in constructing raw moment equations. IF any propensity function is non-polynomial
    then the function `polynomial_propensities` will throw an error.
=#

function isconstant(expr::Union{Symbolic, Number}, vars::Vector, iv::Sym) :: Bool

    # Check that the given expression does NOT depend on the given variables `vars` (expr is constant wrt. vars)
    # A variable here is defined as a function of the independent variable `iv`, e.g,  X(t) is variable, where t
    # is the independent variable

    if istree(expr)
        args = arguments(expr)
        if length(args) == 1
            if isequal(args[1], iv)
                !isvar(expr, vars)
            else
                isconstant(args[1], vars, iv)
            end
        else
            all((isconstant(arg, vars, iv) for arg in args))
        end
    else
        return true
    end

end

isvar(x::Symbolic, vars::Vector) = any((isequal(x, var) for var in vars))

function extract_pwr(expr::Symbolic, smap::Dict, vars::Vector, powers::Vector)
    args = arguments(expr)
    if isvar(args[1], vars) && args[2] isa Int && args[2] >= 0
        idx = smap[args[1]]
        if powers[idx] != 0
            error(args[1], " occurring multiple times in a monomial is unexpected. In ", expr)
        else
            powers[idx] = args[2]
        end
    else
        error("Unexpected term: ", expr)
    end
end

function extract_mul(expr::Symbolic, smap::Dict, vars::Vector, iv::Sym)

    powers = zeros(Int, length(smap))
    factors = []

    for arg in arguments(expr)
        if isconstant(arg, vars, iv)
            push!(factors, arg)
        elseif operation(arg) == ^
            extract_pwr(arg, smap, vars, powers)
        elseif isvar(arg, vars)
            powers[smap[arg]] = 1
        else
            error("Unexpected operation: ", operation(arg), " in ", expr)
        end
    end

    factor = isempty(factors) ? 1 : prod(factors)
    powers, factor

end

function polynomial_propensities(a::Vector, rn::Union{ReactionSystem, ReactionSystemMod}; smap=speciesmap(rn))

    R = length(a)
    N = numspecies(rn)
    vars = [x for (x,y) in Base.sort(collect(smap), by=x->x[2])]
    iv = rn.iv

    all_factors = [[] for i = 1:R]
    all_powers = [[] for i = 1:R]

    for (rind, expr) in enumerate(a)

        expr = expand(expr)

        if isconstant(expr, vars, iv)

            push!(all_factors[rind], expr)
            push!(all_powers[rind], zeros(Int, N))

        elseif isvar(expr, vars)

            push!(all_factors[rind], 1)
            push!(all_powers[rind], map(v -> isequal(expr, v), vars))

        elseif operation(expr) == ^ #Symbolics.Pow

            try
                powers = zeros(Int, N)
                extract_pwr(expr, smap, vars, powers)
                push!(all_factors[rind], 1)
                push!(all_powers[rind], powers)
            catch e
                error("Propensity function ", expr, " is non-polynomial? \n" * string(e))
            end

        elseif operation(expr) == * #Symbolics.Mul

            try
                powers, factor = extract_mul(expr, smap, vars, iv)
                push!(all_powers[rind], powers)
                push!(all_factors[rind], factor)
            catch e
                error("Propensity function ", expr, " is non-polynomial?\n" * string(e))
            end

        elseif operation(expr) == + #Symbolics.Add

            for term in arguments(expr)

                if isconstant(term, vars, iv)
                    push!(all_factors[rind], term)
                    push!(all_powers[rind], zeros(Int, N))
                else
                    try
                        op = operation(term)
                        if op == ^
                            powers = zeros(Int, N)
                            extract_pwr(term, smap, vars, powers)
                            push!(all_factors[rind], 1)
                            push!(all_powers[rind], powers)
                        elseif op == *
                            powers, factor = extract_mul(term, smap, vars, iv)
                            push!(all_factors[rind], factor)
                            push!(all_powers[rind], powers)
                        elseif isvar(term, vars)
                            powers = zeros(Int, N)
                            powers[smap[term]] = 1
                            push!(all_factors[rind], 1)
                            push!(all_powers[rind], powers)
                        else
                            error("Unexpected operation: ", op, " in ", term)
                        end
                    catch e
                        error("Propensity function ", expr, " is non-polynomial?\n" * string(e))
                    end
                end
            end

        else

            error("Expression $expr could not be parsed correctly!")

        end

    end

    max_power = maximum(sum.(vcat(all_powers...)))
    all_factors, all_powers, max_power

end


text:
DIR = @__DIR__

if DIR != nothing
    println("came to if")
    const LIBRDKAFKA = DIR * "/deps" * "/librdkafka.so"
else
    println("came to else")
    const LIBRDKAFKA = Pkg.dir("RDKafka", "deps", "librdkafka.so")
    const LIBRDKAFKA = DIR * "deps/librdkafka.so"
end


text:
"""
    Addition of Two  normal random variables assuming
    zero covariance
"""
+(A::MvNormal, B::MvNormal) = MvNormal(a.μ + b.μ, A.Σ + B.Σ)
-(A::MvNormal, B::MvNormal) = MvNormal(a.μ - b.μ, A.Σ + B.Σ)

chol(A::MvNormal) = A.Σ.chol