text:
mutable struct StandardFieldCache{T1, T2, T3, T4} <: AbstractFieldCache
    Δ::T1
    R::T2
    invRΔ::T2
    IG::T3
    G::T3
    detratio::T4
end

function FieldCache(field, model)
    N = size(field.conf, 1)
    flv = max(nflavors(field), nflavors(model))
    T = interaction_eltype(field)
    GET = greens_eltype(field, model)
    VT = vector_type(GET)
    MT = matrix_type(GET)

    Δ = nflavors(field) == 1 ? zero(T) : vector_type(T)(undef, nflavors(field))

    if flv == 1
        # full greens matrix is flavor symmetric
        R   = zero(GET)
        invRΔ = zero(GET)
        IG  = VT(undef, N)
        G   = VT(undef, N)

    else
        if greens_matrix_type(field, model) <: BlockDiagonal
            # only blocks on diagonals matter
            R   = VT(undef, flv)
            invRΔ = VT(undef, flv)
            IG  = ntuple(_ -> VT(undef, N), flv)
            G   = ntuple(_ -> VT(undef, N), flv)
        else
            # all matter
            R   = MT(undef, flv, flv)
            invRΔ = MT(undef, flv, flv)
            # fast iteration on first index
            IG  = MT(undef, N*flv, flv)
            G   = MT(undef, N*flv, flv)
        end

    end

    detratio = zero(GET)

    return StandardFieldCache(Δ, R, invRΔ, IG, G, detratio)
end

################################################################################
### math for propose_local
################################################################################

# everything should be derivable from Δ and G

# This requires Δ to be set
function calculate_detratio!(cache::StandardFieldCache, G, i)
    calculate_detratio!(cache, cache.Δ, G, i)
end

# TODO Not sure if optimizing this is worth anything
function calculate_detratio!(cache::StandardFieldCache, Δ::Number, G::Union{FMat64, CMat64}, i)
    @inbounds cache.R = 1.0 + Δ * (1.0 - G[i, i])
    cache.detratio = cache.R * cache.R
end

function calculate_detratio!(cache::StandardFieldCache, Δ::FVec64, G::BlockDiagonal{Float64}, i)
    # Unrolled R = I + Δ * (I - G)
    @inbounds for b in eachindex(G.blocks)
        cache.R[b] = 1.0 + Δ[b] * (1.0 - G.blocks[b][i, i])
    end
    # determinant of Diagonal
    cache.detratio = prod(cache.R)
end
# function calculate_detratio!(cache::StandardFieldCache, Δ::FVec64, G::BlockDiagonal{ComplexF64}, i)
#     # Unrolled R = I + Δ * (I - G)
#     @inbounds for b in eachindex(G.blocks)
#         cache.R.re[b] = 1.0 + Δ[b] * (1.0 - G.blocks[b].re[i, i])
#     end
#     @inbounds for b in eachindex(G.blocks)
#         cache.R.im[b] = - Δ[b] * G.blocks.im[b][i, i]
#     end
#     # determinant of Diagonal (don't think this is worth optimizing)
#     cache.detratio = prod(cache.R)
# end
function calculate_detratio!(cache::StandardFieldCache, Δ::CVec64, G::BlockDiagonal{ComplexF64}, i)
    # Unrolled R = I + Δ * (I - G)
    @inbounds for b in eachindex(G.blocks)
        cache.R.re[b] = 1.0 + Δ.re[b] * (1.0 - G.blocks[b].re[i, i])
    end
    @inbounds for b in eachindex(G.blocks)
        cache.R.re[b] += Δ.im[b] * G.blocks[b].im[i, i]
    end
    @inbounds for b in eachindex(G.blocks)
        cache.R.im[b] = Δ.im[b] * (1.0 - G.blocks[b].re[i, i])
    end
    @inbounds for b in eachindex(G.blocks)
        cache.R.im[b] -= Δ.re[b] * G.blocks[b].im[i, i]
    end
    # determinant of Diagonal
    cache.detratio = prod(cache.R)
end

function calculate_detratio!(cache::StandardFieldCache, Δ::FVec64, G::FMat64, i)
    # TODO this is not nice
    N = div(size(G, 1), 2)

    # Unrolled R = I + Δ * (I - G)
    cache.R[1, 1] = 1.0 + Δ[1] * (1.0 - G[i, i])
    cache.R[1, 2] =     - Δ[1] * G[i, i+N]
    cache.R[2, 1] =     - Δ[2] * G[i+N, i]
    cache.R[2, 2] = 1.0 + Δ[2] * (1.0 - G[i+N, i+N])

    # Calculate det of 2x2 Matrix
    # det() vs unrolled: 206ns -> 2.28ns
    cache.detratio = cache.R[1, 1] * cache.R[2, 2] - cache.R[1, 2] * cache.R[2, 1]
end

################################################################################
### Math for accept_local
################################################################################

# This uses cache because it may need to change a float
"""
    vldiv22!(cache, R, Δ)

Computes `cache.invRΔ = R⁻¹ Δ` where `cache.invRΔ`, `R` and `Δ` are a Number,
a 2 component Vector representing a 2x2 Diagonal matrix or a 2x2 Matrix.
"""
function vldiv22!(cache::StandardFieldCache, R::Number, Δ::Number)
    cache.invRΔ = Δ / R
    nothing
end

function vldiv22!(cache::StandardFieldCache, R::FMat64, Δ::FVec64)
    @inbounds begin
        inv_div = 1.0 / cache.detratio
        cache.invRΔ[1, 1] =   R[2, 2] * Δ[1] * inv_div
        cache.invRΔ[1, 2] = - R[1, 2] * Δ[2] * inv_div
        cache.invRΔ[2, 1] = - R[2, 1] * Δ[1] * inv_div
        cache.invRΔ[2, 2] =   R[1, 1] * Δ[2] * inv_div
    end
    nothing
end

function vldiv22!(cache::StandardFieldCache, R::FVec64, Δ::FVec64)
    @inbounds begin
        cache.invRΔ[1] = Δ[1] / R[1]
        cache.invRΔ[2] = Δ[2] / R[2]
    end
    nothing
end

function vldiv22!(cache::StandardFieldCache, R::CVec64, Δ::CVec64)
    @inbounds begin
        # Reminder: 1/c = c* / (cc*) (complex conjugate)
        f1 = 1.0 / (R.re[1] * R.re[1] + R.im[1] * R.im[1])
        f2 = 1.0 / (R.re[2] * R.re[2] + R.im[2] * R.im[2])
        cache.invRΔ.re[1] = f1 * (Δ.re[1] * R.re[1] + Δ.im[1] * R.im[1])
        cache.invRΔ.re[2] = f2 * (Δ.re[2] * R.re[2] + Δ.im[2] * R.im[2])
        cache.invRΔ.im[1] = f1 * (Δ.im[1] * R.re[1] - Δ.re[1] * R.im[1])
        cache.invRΔ.im[2] = f2 * (Δ.im[2] * R.re[2] - Δ.re[2] * R.im[2])
    end
    nothing
end

function update_greens!(cache::StandardFieldCache, G, i, N)
    # calculate Δ R⁻¹
    vldiv22!(cache, cache.R, cache.Δ)

    # calculate (I - G)[:, i:N:end]
    vsub!(cache.IG, I, G, i, N)

    # calculate {Δ R⁻¹} * G[i:N:end, :]
    vmul!(cache.G, cache.invRΔ, G, i, N)

    # update greens function
    # G[m, n] -= {(I - G)[m, i:N:end]} {{Δ R⁻¹} * G[i:N:end, n]}
    vsubkron!(G, cache.IG, cache.G)

    nothing
end

################################################################################
### Field Implementations
################################################################################

maybe_to_float(c::ComplexF64) = abs(imag(c)) < 10eps(real(c)) ? real(c) : c

Base.length(f::AbstractField) = length(conf(f))

function save_field(file, field::AbstractField, entryname="field")
    write(file, entryname * "/VERSION", 1)
    write(file, entryname * "/name", nameof(typeof(field)))
    write(file, entryname * "/conf", conf(field))
end

function load_field(data, ::Val{:Field}, param, model)
    name = data["name"]
    c = data["conf"]
    field = if name == "DensityHirschField"
        DensityHirschField(param, model)
    elseif name == "MagneticHirschField"
        MagneticHirschField(param, model)
    else
        @eval $(Symbol(name))($param, $model)
    end
    copyto!(conf(field), c)

    return field
end

################################################################################
### Hirsch Fields
################################################################################

abstract type AbstractHirschField{T} <: AbstractField end

Base.rand(f::AbstractHirschField) = rand((Int8(-1), Int8(1)), size(f.conf))
Random.rand!(f::AbstractHirschField) = rand!(f.conf, (Int8(-1), Int8(1)))
compress(f::AbstractHirschField) = BitArray(f.conf .== 1)
compressed_conf_type(::AbstractHirschField) = BitArray
decompress(::AbstractHirschField, c) = Int8.(2c .- 1)
decompress!(f::AbstractHirschField, c) = f.conf .= Int8.(2c .- 1)

interaction_eltype(::AbstractHirschField{T}) where {T} = T
interaction_matrix_type(::AbstractHirschField{Float64}, ::Model) = Diagonal{Float64, FVec64}
interaction_matrix_type(::AbstractHirschField{ComplexF64}, ::Model) = Diagonal{ComplexF64, CVec64}

function accept_local!(mc, f::AbstractHirschField, i, slice, args...)
    update_greens!(mc.stack.field_cache, mc.stack.greens, i, size(f.conf, 1))
    @inbounds f.conf[i, slice] *= -1
    nothing
end

########################################
# Density Channel
########################################

"""
    DensityHirschField

Represents the density channel Hirsch transformation of an on-site interaction
term.

exp(Δτ U (n↑ - 0.5)(n↓ - 0.5)) ~ ∑ₓ exp(α x (n↑ + n↓ + 1))
with α = acosh(exp(0.5 Δτ U)) and x ∈ {-1, 1}

This definition uses *positive* U for the attractive case
"""
struct DensityHirschField{T} <: AbstractHirschField{T}
    α::T

    temp_conf::Matrix{Int8}
    conf::Matrix{Int8}
end

function DensityHirschField(param::DQMCParameters, model::Model, U::Number = model.U)
    DensityHirschField(
        maybe_to_float(acosh(exp(0.5 * param.delta_tau * ComplexF64(U)))),
        Matrix{Int8}(undef, length(lattice(model)), param.slices),
        Matrix{Int8}(undef, length(lattice(model)), param.slices)
    )
end

nflavors(::DensityHirschField) = 1

@inline function interaction_matrix_exp!(f::DensityHirschField, result::Diagonal, slice, power)
    # TODO this will index out of bound with 2 flavors
    @inbounds for i in eachindex(result.diag)
        result.diag[i] = exp(power * f.α * f.conf[i, slice])
    end
    nothing
end

function propose_local(mc, f::DensityHirschField, i, slice)
    @inbounds ΔE_boson = -2.0 * f.α * f.conf[i, slice]
    mc.stack.field_cache.Δ = exp(ΔE_boson) - 1
    detratio = calculate_detratio!(mc.stack.field_cache, mc.stack.greens, i)
    return detratio, ΔE_boson, nothing
end

@inline energy_boson(f::DensityHirschField, conf = f.conf) = f.α * sum(conf)

########################################
# Spin Channel
########################################

"""
    MagneticHirschField

Represents the spin channel Hirsch transformation of an on-site interaction
term.

exp(Δτ U (n↑ - 0.5)(n↓ - 0.5)) ~ ∑ₓ exp(α x (n↑ - n↓))
with α = acosh(exp(-0.5 Δτ U)) and x ∈ {-1, 1}
"""
struct MagneticHirschField{T} <: AbstractHirschField{T}
    α::T

    temp_conf::Matrix{Int8}
    conf::Matrix{Int8}
end

function MagneticHirschField(param::DQMCParameters, model::Model, U::Number = model.U)
    MagneticHirschField(
        maybe_to_float(acosh(exp(-0.5 * param.delta_tau * ComplexF64(U)))),
        Matrix{Int8}(undef, length(lattice(model)), param.slices),
        Matrix{Int8}(undef, length(lattice(model)), param.slices)
    )
end

nflavors(::MagneticHirschField) = 2

@inline function interaction_matrix_exp!(f::MagneticHirschField, result::Diagonal, slice, power)
    N = size(f.conf, 1)
    @inbounds for i in axes(f.conf, 1)
        result.diag[i]   = exp( power * f.α * f.conf[i, slice])
    end
    @inbounds for i in axes(f.conf, 1)
        result.diag[i+N] = exp(-power * f.α * f.conf[i, slice])
    end
    nothing
end

function propose_local(mc, f::MagneticHirschField, i, slice)
    ΔE_Boson = -2.0 * f.α * f.conf[i, slice]
    mc.stack.field_cache.Δ[1] = exp(+ΔE_Boson) - 1.0
    mc.stack.field_cache.Δ[2] = exp(-ΔE_Boson) - 1.0
    detratio = calculate_detratio!(mc.stack.field_cache, mc.stack.greens, i)

    # There is no bosonic part (exp(-ΔE_Boson)) to the partition function.
    # Therefore pass 0.0
    return detratio, 0.0, nothing
end

@inline energy_boson(f::MagneticHirschField, conf = f.conf) = 0.0

# DEFAULT
# @bm function propose_global_from_conf(f::MagneticHirschField, mc::DQMC)
# end

################################################################################
### Gauß-Hermite Quadrature
################################################################################

abstract type AbstractGHQField{T} <: AbstractField end

# These represent (-2, -1, +1, 2)
const _GHQVALS = (Int8(1), Int8(2), Int8(3), Int8(4))
Base.rand(f::AbstractGHQField) = rand(_GHQVALS, size(f.conf))
Random.rand!(f::AbstractGHQField) = rand!(f.conf, _GHQVALS)
compressed_conf_type(::AbstractGHQField) = BitArray
function compress(f::AbstractGHQField)
    # converts (1, 2, 3, 4) -> (00, 01, 10, 11)
    BitArray((div(v-1, 2), (v-1) % 2)[step] for v in f.conf for step in (1, 2))
end
function decompress(f::AbstractGHQField, c)
    # converts (00, 01, 10, 11) -> (1, 2, 3, 4)
    conf = similar(f.conf)
    for i in eachindex(conf)
        #  1    +    2 * bit1    +    bit2
        conf[i] = Int8(1) + Int8(2) * c[2i-1] + Int8(c[2i])
    end
    conf
end
function decompress!(f::AbstractGHQField, c)
    for i in eachindex(f.conf)
        #  1    +    2 * bit1    +    bit2
        Int8(1) + Int8(2) * c[i] + Int8(c[i+1])
    end
    f.conf
end

interaction_eltype(::AbstractGHQField{T}) where {T} = T
interaction_matrix_type(::AbstractGHQField{Float64}, ::Model) = Diagonal{Float64, FVec64}
interaction_matrix_type(::AbstractGHQField{ComplexF64}, ::Model) = Diagonal{ComplexF64, CVec64}

"""
    MagneticGHQField

This represents a field generated by using (4 node) Gauss Hermite quadrature to
solve the integral form of the Hubbard Stratonovich transformed interaction. It
is a more general solution than the Hirsch transform, allowing any interaction
that can eb written as V ~ Â² though the current implementation only allows
Hubbard interactions.

In general:
exp(Δτ U Â²) ~ ∑ₓ γ(x) exp(α η(x) Â)
with α = sqrt(Δτ U) and x ∈ {-2, -1, 1, 2}
"""
struct MagneticGHQField{T} <: AbstractGHQField{T}
    α::T
    γ::FVec64
    η::FVec64
    choices::Matrix{Int8}

    temp_conf::Matrix{Int8}
    conf::Matrix{Int8}
end

function MagneticGHQField(param::DQMCParameters, model::Model, U::Number = model.U)
    α = maybe_to_float(sqrt(-0.5 * param.delta_tau * ComplexF64(U)))
    s6 = sqrt(6)
    gammas = Float64[1 - s6/3, 1 + s6/3, 1 + s6/3, 1 - s6/3]
    etas = Float64[-sqrt(6 + 2s6), -sqrt(6 - 2s6), sqrt(6 - 2s6), sqrt(6 + 2s6)]
    choices = Int8[2 3 4; 1 3 4; 1 2 4; 1 2 3]

    MagneticGHQField(
        α, gammas, etas, choices,
        Matrix{Int8}(undef, length(lattice(model)), param.slices),
        Matrix{Int8}(undef, length(lattice(model)), param.slices)
    )
end

nflavors(::MagneticGHQField) = 2
energy_boson(::MagneticGHQField, conf=nothing) = 0.0

# TODO: Maybe worth adding a complex method?
@inline @bm function interaction_matrix_exp!(
        f::MagneticGHQField, result::Diagonal, slice, power
    )
    N = size(f.conf, 1)
    @inbounds for i in 1:N
        result.diag[i]   = exp(+power * f.α * f.η[f.conf[i, slice]])
    end
    @inbounds for i in 1:N
        result.diag[i+N] = exp(-power * f.α * f.η[f.conf[i, slice]])
    end

    return nothing
end

@inline @bm function propose_local(mc, f::MagneticGHQField, i, slice)
    x_old = f.conf[i, slice]
    x_new = @inbounds f.choices[x_old, rand(1:3)]

    exp_ratio = exp(f.α * (f.η[x_new] - f.η[x_old]))
    mc.stack.field_cache.Δ[1] = exp_ratio - 1.0
    mc.stack.field_cache.Δ[2] = 1 / exp_ratio - 1.0
    detratio = calculate_detratio!(mc.stack.field_cache, mc.stack.greens, i)

    return detratio * f.γ[x_new] / f.γ[x_old], 0.0, x_new
end

@inline @bm function accept_local!(
        mc, f::MagneticGHQField, i, slice, detratio, ΔE_boson, x_new
    )
    update_greens!(mc.stack.field_cache, mc.stack.greens, i, size(f.conf, 1))
    @inbounds f.conf[i, slice] = x_new
    return nothing
end

################################################################################
### Util
################################################################################

function ConfigRecorder(::Type{<: AbstractField}, rate = 10)
    ConfigRecorder{BitArray}(rate)
end

function Base.push!(c::ConfigRecorder, field::AbstractField, sweep)
    (sweep % c.rate == 0) && push!(c.configs, compress(field))
    nothing
end

function BufferedConfigRecorder(::Type{<: AbstractField}, filename; rate = 10, chunk_size = 1000)
    BufferedConfigRecorder{BitArray}(filename, rate, chunk_size)
end

function Base.push!(cr::BufferedConfigRecorder, field::AbstractField, sweep)
    (sweep % cr.rate == 0) && _push!(cr, compress(field))
    nothing
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Deprecations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

#                       Deprecations introduced in v0.12
# ==============================================================================

@deprecate(pretty_table(data::Any, header::AbstractVector; kwargs...),
           pretty_table(data; header = header, kwargs...))

@deprecate(pretty_table(io::IO, data::Any, header::AbstractVector; kwargs...),
           pretty_table(io, data; header = header, kwargs...))

@deprecate(pretty_table(::Type{String}, data::Any, header::AbstractVector; kwargs...),
           pretty_table(String, data; header = header, kwargs...))

@deprecate(pretty_table(data::Any, header::AbstractMatrix; kwargs...),
           pretty_table(data;
                        header = tuple([header[i, :] for i = 1:size(header, 1)]...),
                        kwargs...))

@deprecate(pretty_table(io::IO, data::Any, header::AbstractMatrix; kwargs...),
           pretty_table(io, data;
                        header = tuple([header[i, :] for i = 1:size(header, 1)]...),
                        kwargs...))

@deprecate(pretty_table(::Type{String}, data::Any, header::AbstractMatrix; kwargs...),
           pretty_table(String, data;
                        header = tuple([header[i, :] for i = 1:size(header, 1)]...),
                        kwargs...))

# Those definitions are necessary to avoid ambiguities. They **must** be removed
# together with the previous deprecations.
pretty_table(io::IO, data::AbstractVector; kwargs...) =
    _pretty_table(io, data; kwargs...)

pretty_table(io::IO, data::AbstractMatrix; kwargs...) =
    _pretty_table(io, data; kwargs...)

function pretty_table(::Type{String}, data::AbstractVector; kwargs...)
    io = IOBuffer()
    _pretty_table(io, data; kwargs...)
    return String(take!(io))
end

function pretty_table(::Type{String}, data::AbstractMatrix; kwargs...)
    io = IOBuffer()
    _pretty_table(io, data; kwargs...)
    return String(take!(io))
end

ext: jl


text:
using Random
using Test
using SpatialEcology
using DataFrames
using SparseArrays

@testset "Views" begin
    Random.seed!(1337)

    datf = sprand(11,9,0.6)
    datf./=SpatialEcology.colsum(datf)
    com = ComMatrix(datf)

    vcom = view(com, species = [1,2,3,5,7,8], sites = [2,3,6,8,9])
    comc = copy(vcom)

    @test occupancy(vcom) == occupancy(comc)
    @test richness(vcom) == richness(comc)
    @test nsites(vcom) == nsites(comc)
    @test nspecies(vcom) == nspecies(comc)
    @test occurring(vcom) == occurring(comc)
    @test occupied(vcom) == occupied(comc)
    @test noccurring(vcom) == noccurring(comc)
    @test noccupied(vcom) == noccupied(comc)
    @test getspecies(vcom, 4) == getspecies(comc, 4)
    @test nrecords(vcom) == nrecords(comc)
    @test size(vcom) == size(comc)
    @test cooccurring(vcom, 1,3,4) == cooccurring(comc, [1,3,4])
    @test cooccurring(vcom, 3:5) == [false, true, false, true, false]

    @test SpatialEcology.asindices(1:5) == 1:5
    @test SpatialEcology.asindices([1,2,4]) == [1,2,4]
    @test SpatialEcology.asindices(1:3, 2:4) == 1:3
    @test SpatialEcology.asindices(["a", "c"], ["a", "b", "c"]) == [1,3]

end


text:
"""
    with_container(f::Function, ::Type{<:DockerContainer} = DockerContainer; kwargs...)
"""
function with_container(f::F,
                        ::Type{T} = DockerContainer;
                        kwargs...) where {F <: Function, T <: DockerContainer}
    container = T(; kwargs...)
    try
        return f(container)
    finally
        cleanup(container)
    end
end


text:
# Composite Multivariate distribution

## Generic Composite multivariate distribution class

module CompositeDistributions

using Compat

if VERSION >= v"0.7"
  using Statistics
  using Distributed
  import Statistics: mean, median, maximum, minimum, quantile, std, var, cov, cor
else
  using Compat.Statistics
  using Compat.Distributed
  import Base: mean, median, maximum, minimum, quantile, std, var, cov, cor
end

import Base.length, Base.show
using Distributions
import Distributions.params
import Distributions.mean, Distributions.mode, Distributions.var, Distributions.cov
import Distributions.entropy, Distributions.insupport
import Distributions._logpdf, Distributions.logpdf, Distributions._pdf!
import Distributions._rand!
#import Compat.view    # Until v0.5

export AbstractCompositeContinuousDist, ContinuousMultivariateDistribution
export GenericCompositeContinuousDist, CompositeDist
export length, params, set_params!, mean, mode, var, cov, entropy
export insupport, _logpdf, gradlogpdf, _rand!

@compat abstract type AbstractCompositeContinuousDist <: ContinuousMultivariateDistribution end
#@compat abstract type AbstractCompositeDiscreteDist <: DiscreteMultivariateDistribution  end # An idea, but not implemented yet

struct GenericCompositeContinuousDist <: AbstractCompositeContinuousDist
    dist::Vector{ContinuousDistribution}
    indices::Vector{UnitRange{Int64}}

    function GenericCompositeContinuousDist(dist::Vector{ContinuousDistribution})
      length(dist) > 0 || error("number of distributions must be positive")
      #@compat new(dist)
      indices = Array{UnitRange{Int64}}(undef,length(dist))
      idx_start = 1; idx_stop = 0
      for i in 1:length(dist)
         idx_stop += length(dist[i])
         indices[i] = idx_start:idx_stop
         idx_start += length(dist[i])
      end
      new(dist,indices)
    end
end

#typealias CompositeDist  GenericCompositeContinuousDist
const CompositeDist = GenericCompositeContinuousDist

### Parameters
length(d::GenericCompositeContinuousDist, idx::Integer) = length(d.dist[idx])
length(d::GenericCompositeContinuousDist)  = d.indices[end].stop

params(d::GenericCompositeContinuousDist, idx::Integer) = params(d.dist[idx])
function params(d::GenericCompositeContinuousDist)
  n = length(d)
  p = Array{Tuple}(undef,n)
  for i in 1:length(d.dist)
     p[i] = params(d,i)
  end
  return p
end

#= Seems like a good idea, but not clear how to do well
function set_params!(d::GenericCompositeContinuousDist, i::Integer, x::AbstractArray{Tuple,1} ) where T
  param = params(d.dist[i])
  @assert length(param) = length(x)
  for j in 1:length(param)
     param[j] = x[j]
  end
  if length(param) == 1
     d.dist[i] = typeof(d.dist[i])(param[1])
  elseif length(param) == 2
  d.dist[i] = typeof(d.dist[i])(param[1],param[2])
  elseif length(param) == 3
  d.dist[i] = typeof(d.dist[i])(param[1],param[3])
  else
     throw("Too many parameters inside Composite Distribution.")
  end
  return d
end

function set_params!(d::GenericCompositeContinuousDist, x::DenseVector{T}) where T
end
=#

### Basic statistics

mean(d::GenericCompositeContinuousDist, idx::Integer) = mean(d.dist[idx])
function mean(d::GenericCompositeContinuousDist)
  n = length(d)
  mu = Array{Float64}(undef,n)
  for i in 1:length(d.dist)
     mu[d.indices[i]] = mean(d,i)
  end
  return mu
end

mode(d::GenericCompositeContinuousDist, idx::Integer) = mode(d.dist[idx])
function mode(d::GenericCompositeContinuousDist)
  n = length(d)
  mo = Array{Float64}(undef,n)
  for i in 1:length(d.dist)
     mo[d.indices[i]] = mode(d,i)
  end
  return mo
end

var(d::GenericCompositeContinuousDist, idx::Integer) = var(d.dist[idx])
function var(d::GenericCompositeContinuousDist)
  n = length(d)
  v = Array{Float64}(undef,n)
  for i in 1:length(d.dist)
     v[d.indices[i]] = var(d,i)
  end
  return v
end

cov(d::GenericCompositeContinuousDist, idx::Integer) = cov(d.dist[idx])
function cov(d::GenericCompositeContinuousDist)
  n = length(d)
  covar = zeros(Float64,(n,n) )
  for i in 1:length(d.dist)
     if length(d.dist[i]) == 1
        covar[d.indices[i]] = var(d,i)
     else
        covar[d.indices[i],d.indices[i]] = cov(d,i)
     end
  end
  return covar
end

entropy(d::GenericCompositeContinuousDist, idx::Integer) = entropy(d.dist[idx])
function entropy(d::GenericCompositeContinuousDist)
  mapreduce(entropy,+,d.dist)
  #=sum = 0
  for i in 1:length(d.dist)
    sum += entropy(d,i)
  end
  return sum =#
end

### Evaluation

function index(d::GenericCompositeContinuousDist, idx::Integer)
  # note that this function is intentionally not type stable.
  # Needs to provide scalar index for univariate distributions or range of indices for multivariate distributions
  #if d.indices[idx].start == d.indices[idx].stop
  if isa(d.dist[idx], UnivariateDistribution)
     return d.indices[idx].start
  else
     return d.indices[idx]
  end
end

function insupport(d::GenericCompositeContinuousDist, x::AbstractVector{T})  where T<:Real
  if ! (length(d) == length(x)) return false end
  for i in 1:length(d.dist)
    if !insupport(d.dist[i],x[index(d,i)]) return false end
  end
  return true
end

function Distributions._logpdf(d::GenericCompositeContinuousDist, x::AbstractArray{T,2}) where T<:Real
  sum = zeros(T,size(x,2))
  for j in 1:size(x,1)
     for i in 1:length(d.dist)
       sum[j] += logpdf(d.dist[i],x[j,index(d,i)])
     end
  end
  return sum
end

function Distributions._logpdf(d::GenericCompositeContinuousDist, x::AbstractArray{T,1}) where T<:Real
  sum = zero(T)
  for i in 1:length(d.dist)
    #sum += _logpdf(d.dist[i],x[index(d,i)])
    sum += logpdf(d.dist[i],x[index(d,i)])
  end
  return sum
end

function gradlogpdf(d::GenericCompositeContinuousDist, x::DenseVector{T}) where T<:Real
    z = Array{T}(undef,length(d))
    for i = 1:length(d.dist)
        z[index(d,i)] = gradlogpdf(d.dist[i],view(x,index(d,i)))
    end
    return z
end

### Sampling

function Distributions._rand!(d::GenericCompositeContinuousDist, x::DenseVector{T}) where T<:Real
    for i = 1:length(d.dist)
        _rand!(d.dist[i],view(x,index(d,i)))
    end
    return x
end

function _rand!(d::GenericCompositeContinuousDist, x::DenseMatrix{T}) where T<:Real
    for i = 1:length(d.dist)
        _rand!(d.dist[i],view(x,index(d,i),:))   # Check got dimensions right
    end
    return x
end

### Show

distrname(d::GenericCompositeContinuousDist) = "GenericCompositeContinuous"
function Base.show(io::IO, d::GenericCompositeContinuousDist)
  print(io,distrname(d) * ":\n")
  for i in 1:length(d.dist)
     show(io, d.dist[i] )
     print(io,"\n")
  end
end

end # module


text:
################################################################################
#  Copyright 2020, Marta Vanin, Tom Van Acker                                  #
################################################################################
# PowerModelsDistributionStateEstimation.jl                                    #
# An extention package of PowerModelsDistribution.jl for Static Power System   #
# State Estimation.                                                            #
################################################################################

mutable struct ReducedACPUPowerModel <: _PMD.AbstractUnbalancedACPModel _PMD.@pmd_fields end
mutable struct ReducedACRUPowerModel <: _PMD.AbstractUnbalancedACRModel _PMD.@pmd_fields end

AbstractReducedModel = Union{ReducedACRUPowerModel, ReducedACPUPowerModel}

"Power balance constraint for the reduced ACR and ACP formulations.
These formulation are exact for networks like those made available in the ENWL database,
where there are no gound admittance, storage elements and active switches.
Other than this, the function is the same as the constraint_mc_load_power_balance defined in PowerModelsDistribution "
function constraint_mc_power_balance(pm::AbstractReducedModel, i::Int; nw::Int=_IM.nw_id_default)

    bus = _PMD.ref(pm, nw, :bus, i)
    bus_arcs = _PMD.ref(pm, nw, :bus_arcs_conns_branch, i)
    bus_arcs_trans = _PMD.ref(pm, nw, :bus_arcs_conns_transformer, i)
    bus_gens = _PMD.ref(pm, nw, :bus_conns_gen, i)
    bus_loads = _PMD.ref(pm, nw, :bus_conns_load, i)

    p    = get(_PMD.var(pm, nw),    :p, Dict()); _PMD._check_var_keys(p, bus_arcs, "active power", "branch")
    q    = get(_PMD.var(pm, nw),    :q, Dict()); _PMD._check_var_keys(q, bus_arcs, "reactive power", "branch")
    pg   = get(_PMD.var(pm, nw),   :pg, Dict()); _PMD._check_var_keys(pg, bus_gens, "active power", "generator")
    qg   = get(_PMD.var(pm, nw),   :qg, Dict()); _PMD._check_var_keys(qg, bus_gens, "reactive power", "generator")
    pt   = get(_PMD.var(pm, nw),   :pt, Dict()); _PMD._check_var_keys(pt, bus_arcs_trans, "active power", "transformer")
    qt   = get(_PMD.var(pm, nw),   :qt, Dict()); _PMD._check_var_keys(qt, bus_arcs_trans, "reactive power", "transformer")
    pd   = get(_PMD.var(pm, nw),  :pd, Dict()); _PMD._check_var_keys(pd, bus_loads, "active power", "load")
    qd   = get(_PMD.var(pm, nw),  :qd, Dict()); _PMD._check_var_keys(pd, bus_loads, "reactive power", "load")

    terminals = bus["terminals"]
    grounded =  bus["grounded"]

    ungrounded_terminals = [(idx,t) for (idx,t) in enumerate(terminals) if !grounded[idx]]

    for (idx,t) in ungrounded_terminals
        cp = JuMP.@constraint(pm.model,
        sum(  p[a][t] for (a, conns) in bus_arcs if t in conns)
      + sum( pt[a][t] for (a, conns) in bus_arcs_trans if t in conns)
      - sum( pg[g][t] for (g, conns) in bus_gens if t in conns)
      + sum( pd[l][t] for (l, conns) in bus_loads if t in conns)
      == 0.0
      )

        cq = JuMP.@constraint(pm.model,
        sum(  q[a][t] for (a, conns) in bus_arcs if t in conns)
      + sum( qt[a][t] for (a, conns) in bus_arcs_trans if t in conns)
      - sum( qg[g][t] for (g, conns) in bus_gens if t in conns)
      + sum( qd[l][t] for (l, conns) in bus_loads if t in conns)
      == 0.0
      )
   end
end

"If the formulation is not reduced, delegates back to PowerModelsDistribution"
function constraint_mc_power_balance(pm::_PMD.AbstractUnbalancedPowerModel, i::Int; nw::Int=_IM.nw_id_default)
    _PMD.constraint_mc_power_balance(pm, i; nw=nw)
end

function variable_mc_bus_voltage(pm::ReducedACPUPowerModel; bounded = true)
    _PMD.variable_mc_bus_voltage_angle(pm; bounded = bounded)
    _PMD.variable_mc_bus_voltage_magnitude_only(pm; bounded = bounded)
end

function variable_mc_bus_voltage(pm::ReducedACRUPowerModel; bounded = true)
    _PMD.variable_mc_bus_voltage(pm; bounded = bounded)
end

"If the formulation is not reduced, delegates back to PowerModelsDistribution"
function variable_mc_bus_voltage(pm::_PMD.AbstractUnbalancedPowerModel; bounded = true)
    _PMD.variable_mc_bus_voltage(pm; bounded = bounded)
end


text:
# ensure consistency in the command line interface
using ArgParse

include(normpath(joinpath(Base.source_dir(),"..","src", "parsing.jl")))

let s = build_arg_table()

    test_arg(args) = parse_args(args, s)

    @test test_arg(["ls"]) == Dict{AbstractString, Any}("%COMMAND%" => "ls", "ls" => Dict{AbstractString, Any}())
    @test test_arg(["config", "CONFIG", "OUTPUT"]) ==
        Dict{AbstractString,Any}("%COMMAND%"=>"config","config"=>Dict{AbstractString,Any}("config_file"=>"CONFIG","output_dir"=>"OUTPUT", "no-graph" => false))
    @test test_arg(["exp", "-p", "1", "EXP", "OUTPUT"]) ==
        Dict{AbstractString,Any}("exp"=>Dict{AbstractString,Any}("output_file"=>"OUTPUT","analysis_file"=>"EXP","addprocs"=>1,"debug"=>false),"%COMMAND%"=>"exp")
    @test test_arg(["exp", "EXP", "OUTPUT"]) ==
        Dict{AbstractString,Any}("exp"=>Dict{AbstractString,Any}("output_file"=>"OUTPUT","analysis_file"=>"EXP","addprocs"=>0,"debug"=>false),"%COMMAND%"=>"exp")

end


text:
module Heptapus

using TypedTables, CSV, CUDAdrv

export empiricalbandwidth, Roofline

"""
    empiricalbandwidth(nbytes=2*1024^3; devicenumber=0, numtests=10)

Compute the emperical bandwidth in GB/s of a CUDA device using `nbytes` of
memory.  The device to test can be slected with `devicenumber` and the
bandwidth is an average of `ntests`.
"""
function empiricalbandwidth(nbytes=2*1024^3; devicenumber=0, ntests=10)
    dev = CuDevice(devicenumber)
    ctx = CuContext(dev)

    a = Mem.alloc(nbytes÷2)
    b = Mem.alloc(nbytes÷2)

    Mem.transfer!(a, b, nbytes÷2, async=true)
    Mem.transfer!(b, a, nbytes÷2, async=true)

    t = CUDAdrv.@elapsed for n = 1:ntests
        Mem.transfer!(a, b, nbytes÷2, async=true)
        Mem.transfer!(b, a, nbytes÷2, async=true)
    end

    Mem.free(a)
    Mem.free(b)

    bandwidth = 2*nbytes*ntests/(t*1e9)
end

"""
    Roofline(command::Cmd)

Use `nvprof` to profile `command` and compute for each kernel executed:

  - arithmetic intensity;
  - performance (GFLOP/s);
  - kernel max empirical bandwidth (GB/s);
  - max GFLOP/s estimate (GFLOP/s);
  - max empirical bandwidth (GB/s);
  - local loads or stores (`Bool`);
  - floaing point type used; and
  - mixed floating point operations (`Bool`).
"""
struct Roofline
    t::Table

    function Roofline(command::Cmd)
        s = mktemp() do f, _
            metrics = [:dram_read_throughput,
                       :dram_write_throughput,
                       :dram_write_transactions,
                       :dram_read_transactions,
                       :local_load_transactions,
                       :local_store_transactions,
                       :flop_hp_efficiency,
                       :flop_sp_efficiency,
                       :flop_dp_efficiency,
                       :flop_count_hp,
                       :flop_count_sp,
                       :flop_count_dp]
            cmd = `nvprof -u ms --csv --metrics $(join(metrics,",")) --log-file $f $command`
            @info "Getting metrics" cmd
            run(cmd)
            Table(CSV.File(f, comment="=", allowmissing=:none))
        end

        t = mktemp() do f, _
            cmd = `nvprof -u ms --csv --log-file $f $command`
            @info "Getting timings" cmd
            run(cmd)
            Table(CSV.File(f, comment="=", allowmissing=:none, datarow=3))
        end

        kernels = unique(s.Kernel)

        @info "Kernels found in timings" kernels
        @info "Kernels found in measurements" unique(t.Name)

        if !(kernels ⊆ unique(t.Name))
          error("Kernel names do not match")
        end

        getmetric(T, k, m; trim=0) =
            parse(T, s[map(row -> row.Kernel == k &&
                           getproperty(row, Symbol("Metric Name")) == String(m),
                           s)][1].Avg[1:end-trim])

        # get average kernel execution time in seconds
        gettime(k) = t[map(row -> row.Name == k, t)][1].Avg/1e3

        eb = empiricalbandwidth()
        @info "Maximum empirical bandwidth $eb (GB/s)"

        maxempiricalbandwidth = fill(eb, length(kernels))
        arithmeticintensity = zeros(length(kernels))
        performance = similar(arithmeticintensity)
        kernelmaxempiricalbandwidth = similar(arithmeticintensity)
        maxgflopsestimate = similar(arithmeticintensity)
        haslocal = zeros(Bool, length(kernels))
        hasmixedflops = zeros(Bool, length(kernels))
        floptype = Array{Type}(undef, length(kernels))

        @info "Computing results"
        for (i, k) in enumerate(kernels)
            dram_write_throughput = getmetric(Float64, k, :dram_write_throughput, trim=4)
            dram_read_throughput = getmetric(Float64, k, :dram_read_throughput, trim=4)

            dram_write_transactions = getmetric(Int64, k, :dram_write_transactions)
            dram_read_transactions = getmetric(Int64, k, :dram_read_transactions)

            local_load_transactions = getmetric(Int64, k, :local_load_transactions)
            local_store_transactions = getmetric(Int64, k, :local_store_transactions)

            flop_efficiency = (hp=getmetric(Float64, k, :flop_hp_efficiency, trim=1),
                               sp=getmetric(Float64, k, :flop_sp_efficiency, trim=1),
                               dp=getmetric(Float64, k, :flop_dp_efficiency, trim=1))

            flop_count = (hp=getmetric(Int64, k, :flop_count_hp),
                          sp=getmetric(Int64, k, :flop_count_sp),
                          dp=getmetric(Int64, k, :flop_count_dp))

            elapsedtime = gettime(k)

            if local_load_transactions > 0 || local_store_transactions > 0
                haslocal[i] = true
                @warn """Kernel $k has nonzero load or store transactions.

                    This could be evidece of register spilling and the returned
                    roofline numbers will not be accurate for this kernel.

                """ local_load_transactions local_store_transactions
            end

            if sum(flop_count) != maximum(flop_count)
                hasmixedflops[i] = true
                @warn """Kernel $k has floating point operations for multiple types.

                    The performance results will use the floating point type with
                    the maximum flop count.  The may or may not be what is desired.

                    Note: `flop_count` contains half, single, and double precision.
                """ flop_count
            end

            p = argmax(flop_count)
            floptype[i] = (hp=Float16, sp=Float32, dp=Float64)[p]

            # Note that each transaction is 32 bytes
            bytes = 32(dram_write_transactions + dram_read_transactions)
            flops = flop_count[p]

            arithmeticintensity[i] = flops/bytes

            # in GFLOP/s
            performance[i] = (flops/elapsedtime)/1e9
            maxgflopsestimate[i] = 100performance[i]/flop_efficiency[p]

            kernelmaxempiricalbandwidth[i] = empiricalbandwidth(bytes)
        end

        new(Table(kernels = kernels,
                  arithmeticintensity = arithmeticintensity,
                  performance = performance,
                  kernelmaxempiricalbandwidth = kernelmaxempiricalbandwidth,
                  maxgflopsestimate = maxgflopsestimate,
                  maxempiricalbandwidth = maxempiricalbandwidth,
                  haslocal = haslocal,
                  floptype = floptype,
                  hasmixedflops = hasmixedflops))
    end
end

function Base.show(io::IO, r::Roofline)
    print(io, "Roofline containing ")
    show(io, MIME"text/plain"(), r.t)
end

end # module


text:
using ChemicalKinetics

println("    Testing simulation from Fitzhugh Nagumo ODEs...")

odeModel = OdeModel("fitzhughNagumoModel.txt")
# odeModel = OdeModel(joinpath(dirname(@__FILE__), "fitzhughNagumoModel.txt"))

# ckCvode = CkCvode(odeModel, [i for i in 0:0.5:25], 1e-4, [1e-8, 1e-14, 1e-6])
ckCvode = CkCvode(odeModel, [i for i in 0:0.5:25])

# ode_simulation = ckCvode.simulate_odes()
# println("Simulated values of kinetic equations:")
# println(ode_simulation)

# Simulate data from kinetic equations with noise variance = 0.1
ode_simulation, data_simulation = ckCvode.simulate_data(0.1)
println("Simulated values of kinetic equations:")
println(ode_simulation)
println("\nSimulated data:")
println(data_simulation)


text:
# actors.jl

using Rocket

struct StmActor <: Actor{String} end
Rocket.on_next!(actor::StmActor, data::String) = begin
    if data == "READ_STMS"
        stms = AppliAR.read_bank_statements("./bank.csv")
        @show(stms)
        subscribe!(from(stms), ar_actor)
    end
end
Rocket.on_complete!(actor::StmActor) = @info("StmActor completed!")
Rocket.on_error!(actor::StmActor, err) = @info(error(err))

struct SalesActor <: Actor{String} end
Rocket.on_next!(actor::SalesActor, data::String) = begin
    if data == "START"
        #ar_actor = ARActor()
        orders = @fetch AppliSales.process()
        subscribe!(from(orders), ar_actor)
    end
end
Rocket.on_complete!(actor::SalesActor) = @info("SalesActor completed!")
Rocket.on_error!(actor::SalesActor, err) = @info(error(err))

struct ARActor <: Actor{Any}
    ar_pid::Int64
    ARActor(ar_pid) = new(ar_pid)
end
Rocket.on_next!(actor::ARActor, data::AppliSales.Order) = begin
        d = @fetchfrom actor.ar_pid AppliAR.process([data])
        subscribe!(from(d), gl_actor)
end
Rocket.on_next!(actor::ARActor, data::AppliAR.BankStatement) = begin
        unpaid_inv = @fetchfrom actor.ar_pid retrieve_unpaid_invoices()
        entries = @fetchfrom actor.ar_pid AppliAR.process(unpaid_inv, [data])
        subscribe!(from(entries), gl_actor)
end
Rocket.on_complete!(actor::ARActor) = begin
    @info("ARActor Completed!")
end
Rocket.on_error!(actor::ARActor, err) = @info(error(err))

struct GLActor <: Actor{Any}
    gl_pid::Int64
    GLActor(gl_pid) = new(gl_pid)
end
Rocket.on_next!(actor::GLActor, data::Any) = begin
    if data isa AppliGeneralLedger.JournalEntry
        result = @fetchfrom actor.gl_pid AppliGeneralLedger.process([data])
    end
end
Rocket.on_complete!(actor::GLActor) = @info("GLActor completed!")
Rocket.on_error!(actor::GLActor, err) = @info(error(err))


text:
module bases

export Basis, GenericBasis, CompositeBasis, basis,
       tensor, ⊗, ptrace, permutesystems,
       IncompatibleBases,
       samebases, multiplicable,
       check_samebases, check_multiplicable

import Base: ==, ^

using Compat

"""
Abstract base class for all specialized bases.

The Basis class is meant to specify a basis of the Hilbert space of the
studied system. Besides basis specific information all subclasses must
implement a shape variable which indicates the dimension of the used
Hilbert space. For a spin-1/2 Hilbert space this would be the
vector `Int[2]`. A system composed of two spins would then have a
shape vector `Int[2 2]`.

Composite systems can be defined with help of the [`CompositeBasis`](@ref)
class.
"""
abstract type Basis end

==(b1::Basis, b2::Basis) = false

"""
    length(b::Basis)

Total dimension of the Hilbert space.
"""
Base.length(b::Basis) = prod(b.shape)

"""
    basis(a)

Return the basis of an object.

If it's ambiguous, e.g. if an operator has a different left and right basis,
an [`IncompatibleBases`](@ref) error is thrown.
"""
function basis end

"""
    GenericBasis(N)

A general purpose basis of dimension N.

Should only be used rarely since it defeats the purpose of checking that the
bases of state vectors and operators are correct for algebraic operations.
The preferred way is to specify special bases for different systems.
"""
mutable struct GenericBasis <: Basis
    shape::Vector{Int}
end

GenericBasis(N::Int) = GenericBasis(Int[N])

==(b1::GenericBasis, b2::GenericBasis) = equal_shape(b1.shape, b2.shape)

"""
    CompositeBasis(b1, b2...)

Basis for composite Hilbert spaces.

Stores the subbases in a vector and creates the shape vector directly
from the shape vectors of these subbases. Instead of creating a CompositeBasis
directly `tensor(b1, b2...)` or `b1 ⊗ b2 ⊗ …` can be used.
"""
mutable struct CompositeBasis <: Basis
    shape::Vector{Int}
    bases::Vector{Basis}
end
CompositeBasis(bases::Vector{Basis}) = CompositeBasis(Int[prod(b.shape) for b in bases], bases)
CompositeBasis(bases::Basis...) = CompositeBasis(Basis[bases...])

==(b1::CompositeBasis, b2::CompositeBasis) = equal_shape(b1.shape, b2.shape) && equal_bases(b1.bases, b2.bases)

"""
    tensor(x, y, z...)

Tensor product of the given objects. Alternatively, the unicode
symbol ⊗ (\\otimes) can be used.
"""
tensor() = throw(ArgumentError("Tensor function needs at least one argument."))
tensor(b::Basis) = b

"""
    tensor(x::Basis, y::Basis, z::Basis...)

Create a [`CompositeBasis`](@ref) from the given bases.

Any given CompositeBasis is expanded so that the resulting CompositeBasis never
contains another CompositeBasis.
"""
tensor(b1::Basis, b2::Basis) = CompositeBasis(Int[prod(b1.shape); prod(b2.shape)], Basis[b1, b2])
tensor(b1::CompositeBasis, b2::CompositeBasis) = CompositeBasis(Int[b1.shape; b2.shape], Basis[b1.bases; b2.bases])
function tensor(b1::CompositeBasis, b2::Basis)
    N = length(b1.bases)
    shape = Vector{Int}(N+1)
    bases = Vector{Basis}(N+1)
    for i=1:N
        shape[i] = b1.shape[i]
        bases[i] = b1.bases[i]
    end
    shape[end] = prod(b2.shape)
    bases[end] = b2
    CompositeBasis(shape, bases)
end
function tensor(b1::Basis, b2::CompositeBasis)
    N = length(b2.bases)
    shape = Vector{Int}(N+1)
    bases = Vector{Basis}(N+1)
    for i=1:N
        shape[i+1] = b2.shape[i]
        bases[i+1] = b2.bases[i]
    end
    shape[1] = prod(b1.shape)
    bases[1] = b1
    CompositeBasis(shape, bases)
end
tensor(bases::Basis...) = reduce(tensor, bases)
⊗ = tensor

function ^(b::Basis, N::Int)
    if N < 1
        throw(ArgumentError("Power of a basis is only defined for positive integers."))
    end
    tensor([b for i=1:N]...)
end

"""
    equal_shape(a, b)

Check if two shape vectors are the same.
"""
function equal_shape(a::Vector{Int}, b::Vector{Int})
    if a === b
        return true
    end
    if length(a) != length(b)
        return false
    end
    for i=1:length(a)
        if a[i]!=b[i]
            return false
        end
    end
    return true
end

"""
    equal_bases(a, b)

Check if two subbases vectors are identical.
"""
function equal_bases(a::Vector{T}, b::Vector{T}) where T <: Basis
    if a===b
        return true
    end
    for i=1:length(a)
        if a[i]!=b[i]
            return false
        end
    end
    return true
end

"""
Exception that should be raised for an illegal algebraic operation.
"""
mutable struct IncompatibleBases <: Exception end

"""
    samebases(a, b)

Test if two objects have the same bases.
"""
samebases(b1::Basis, b2::Basis) = b1==b2

"""
    check_samebases(a, b)

Throw an [`IncompatibleBases`](@ref) error if the objects don't have
the same bases.
"""
function check_samebases(b1, b2)
    if !samebases(b1, b2)
        throw(IncompatibleBases())
    end
end

"""
    multiplicable(a, b)

Check if two objects are multiplicable.
"""
multiplicable(b1::Basis, b2::Basis) = b1==b2

function multiplicable(b1::CompositeBasis, b2::CompositeBasis)
    if !equal_shape(b1.shape,b2.shape)
        return false
    end
    for i=1:length(b1.shape)
        if !multiplicable(b1.bases[i], b2.bases[i])
            return false
        end
    end
    return true
end

"""
    check_multiplicable(a, b)

Throw an [`IncompatibleBases`](@ref) error if the objects are
not multiplicable.
"""
function check_multiplicable(b1, b2)
    if !multiplicable(b1, b2)
        throw(IncompatibleBases())
    end
end

"""
    ptrace(a, indices)

Partial trace of the given basis, state or operator.

The `indices` argument, which can be a single integer or a vector of integers,
specifies which subsystems are traced out. The number of indices has to be
smaller than the number of subsystems, i.e. it is not allowed to perform a
full trace.
"""
function ptrace(b::CompositeBasis, indices::Vector{Int})
    J = [i for i in 1:length(b.bases) if i ∉ indices]
    if length(J)==0
        throw(ArgumentError("Tracing over all indices is not allowed in ptrace."))
    elseif length(J)==1
        return b.bases[J[1]]
    else
        return CompositeBasis(b.shape[J], b.bases[J])
    end
end
ptrace(a, index::Int) = ptrace(a, Int[index])

"""
    permutesystems(a, perm)

Change the ordering of the subsystems of the given object.

For a permutation vector `[2,1,3]` and a given object with basis `[b1, b2, b3]`
this function results in `[b2, b1, b3]`.
"""
function permutesystems(b::CompositeBasis, perm::Vector{Int})
    @assert length(b.bases) == length(perm)
    @assert isperm(perm)
    CompositeBasis(b.shape[perm], b.bases[perm])
end

end # module


text:
@testset "$TEST $G" begin

g4 = PathDiGraph(5, DG)
d1 =EdgeMap(g4, float([ 0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0]))
d2 =EdgeMap(g4, sparse(float([ 0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])))

y = bellman_ford_shortest_paths(g4, 2, d1)
z = bellman_ford_shortest_paths(g4, 2, d2)
@test y.dists == z.dists == [Inf, 0, 6, 17, 33]
@test enumerate_paths(z)[2] == []
@test enumerate_paths(z)[4] == enumerate_paths(z,4) == [2,3,4]
@test !has_negative_edge_cycle(g4)

z = bellman_ford_shortest_paths(g4, 2)
@test z.dists == [typemax(Int), 0, 1, 2, 3]

end # testset


text:
# TODO: this is not at all tested yet!!!!!!!

# we need this to evaluate the annoying integrand in the displacement field
using JuLIPMaterials: Vec3, Mat3, Ten33, Ten43
using JuLIPMaterials.CLE: onb3D

using Einsum, StaticArrays, LinearAlgebra
using GaussQuadrature: legendre

export Dislocation, IsoEdgeDislocation3D, IsoScrewDislocation3D

abstract type AbstractDislocation{T} end

struct Dislocation3D{T} <: AbstractDislocation{T}
   Nquad::Int
   b::Vec3{T}
   t::Vec3{T}
   cut::Vec3{T}
   C::Ten43{T}
   remove_singularity::Bool
end

"""
`Dislocation3D`
construct a dislocation type
"""
function Dislocation(b::Vec3, t::Vec3, C::Ten43; Nquad = nothing, cut = _autocut_(b,t), remove_singularity = true)
   if Nquad == nothing
      error("still need to implement auto-tuning, please provide Nquad")
   end
   if abs(t ⋅ cut) > 1e-12
      error("Cut direction is not orthogonal to dislocation tangent!")
   end
   return Dislocation3D(Nquad, b, t, cut, C, remove_singularity)
end

function _autocut_(b::Vec3, t::Vec3)
   # Normalise vectors
   t /= norm(t)
   b /= norm(b)

   # Construct a cut plane "automatically", depending on whethere screw or edge
   if abs(b⋅t) > 0.9
      cut, _ = onb3D(t)
   else
      cut = b-(b⋅t)*t
      cut /= norm(cut)
   end

	return Vec3(cut)
end

Dislocation(b::Array{T,1},t::Array{T,1},C::Array{T, 4}; kwargs...) where {T} =
		Dislocation(Vec3{T}(b),Vec3{T}(t),Ten43{T}(C); kwargs...)

function (Disl::Dislocation3D{T})(x) where T
   if Disl.remove_singularity && norm(x) < 1e-10
      return @SVector zeros(T, 3)
   end
   return eval_dislocation(Vec3(x), Disl.b, Disl.t, Disl.cut, Disl.C, Disl.Nquad)
end

function grad(Disl::Dislocation3D{T}, x) where T
   if Disl.remove_singularity && norm(x) < 1e-10
      return @SMatrix zeros(T, 3, 3)
   end
   return grad_dislocation(Vec3(x), Disl.b, Disl.t, Disl.cut, Disl.C, Disl.Nquad)
end

# ========== Edge dislocation isotropic solid ==============
"""
`u_edge_isotropic(x, y, b, ν) -> u_x, u_y`
compute the displacement field `ux, uy` for an edge dislocation in an
isotropic linearly elastic medium, with core at (0,0),
burgers vector `b * [1.0;0.0]` and Poisson ratio `ν`
This is to be used primarily for comparison, since the exact solution will
not be the isotropic elasticity solution.
"""
function u_edge_isotropic(x, y, b, ν)
    # @warn("This function is to be replaced by new Isotropic Edge Dislocation type")
    r² = x.^2 + y.^2
    ux = b/(2*π) * ( angle.(x + im*y) + (x .* y) ./ (2*(1-ν) * r²) )
    uy = -b/(2*π) * ( (1-2*ν)/(4*(1-ν)) * log.(r²) - 2 * y.^2 ./ (4*(1-ν) * r²) )
    return ux, uy
end

# ========== The Stroh / Hirth&Lothe Horror! =======================

include("sextic.jl")

# ========== Implementation of the BBS79 formula for a dislocation =============

function QSB(C, m0::Vec3{TT}, n0::Vec3{TT}, Nquad) where TT
   Q, S, B = zero(Mat3{TT}), zero(Mat3{TT}), zero(Mat3{TT})
   nn, nm, mm = zero(MMat3{TT}), zero(MMat3{TT}), zero(MMat3{TT})

	# Since we integrate a periodic function, use the trapezium rule.
   for ω in range(0, step=pi/Nquad, length=Nquad)
      m = cos(ω) * m0 + sin(ω) * n0
      n = -sin(ω) * m0 + cos(ω) * n0
      @einsum nn[i,j] = n[α] * C[i,α,j,β] * n[β]
      @einsum nm[i,j] = n[α] * C[i,α,j,β] * m[β]
      @einsum mm[i,j] = m[α] * C[i,α,j,β] * m[β]
      nn⁻¹ = inv(nn)
      Q += nn⁻¹                           # (3.6.4)
      S += nn⁻¹ * nm                      # (3.6.6)
      B += mm - nm' * nn⁻¹ * nm           # (3.6.9) and using  mn = nm'
   end

   return Q * (-1/Nquad), S * (-1/Nquad), B * (1/4/Nquad/pi)
end

function eval_dislocation(x::AbstractVector{TT}, b, t, cut, C, Nquad=10) where TT
   # normalise dislocation tangent direction
   t /= norm(t)
   # fixed coordinate system w.r.t cut, from which we compute ω
   m0 = cut
	n0 = t × m0

   # project x into the plane normal to t, this will not change the value of u
   x -= (t ⋅ x) * t

   # construct a right--handed ONB (t, m, n)
   mω = x / norm(x)   # p.145, l.7
   nω = t × mω

   # compute x ⤅ (r, ω)
   r = norm(x)
   ω = angle(mω ⋅ m0 , mω ⋅ n0 )
	# Ensure this is positive
   if ω < 0.0
      ω += 2*π
   end

	# ------------ Implement components of (4.1.25) ------------
   # NB: 4.1.25 is incorrect!!! Instead integrate 4.1.24 correctly (Q,S,B)
   # depend upon angle in general.

   nn, nm, mm = zero(MMat3{TT}), zero(MMat3{TT}), zero(MMat3{TT})
   T2, T3 = zero(Mat3{TT}), zero(Mat3{TT})

   # compute S for Term 1
   _, Sω, _ = QSB(C, mω, nω, Nquad)

   # compute Terms 2 and 3 via Legendre quadrature
   Xquad, Wquad = legendre(Float64, 2*Nquad)
   Xquad = ω * (1.0 .+ Xquad) / 2.0     # now Xquad ranges from 0.0 to ω
   Wquad = Wquad * (ω / sum(Wquad))
   for (ξ, dξ) in zip(Xquad, Wquad)
      m = cos(ξ) * m0 + sin(ξ) * n0
      n = -sin(ξ) * m0 + cos(ξ) * n0
      # get the S, B tensors
      _, Sξ, Bξ = QSB(C, m, n, Nquad)
      # compute nn, nm
      @einsum nn[i,j] = n[α] * C[i,α,j,β] * n[β]
      @einsum nm[i,j] = n[α] * C[i,α,j,β] * m[β]
      nn⁻¹ = inv(nn)
      T2 += dξ * (nn⁻¹ * Bξ)
      T3 += dξ * (nn⁻¹ * nm * Sξ)
   end
   #---------- put everything together -------------  (4.1.25)
   u = (- Sω * log(r) + 4*π * T2 + T3) * b / (2*π)
   return Vec3(u)
end

function grad_dislocation(x::AbstractVector{TT}, b, t, cut, C, Nquad=10) where TT
   #x, b, t, C = Vec3(x), Vec3(b), Vec3(t), Ten43(C)
   t /= norm(t)
	m0 = cut
	n0 = t × m0

   Q, S, B = QSB(C, m0, n0, Nquad)
   # compute displacement gradient
   # project x to the (m0, n0) plane and compute the (m, n) vectors
   x -= (x⋅t) * t
   r = norm(x)
   m = x / norm(x)
   n = t × m
   #  Implement (4.1.16)
   nn, nm, mm = zero(MMat3{TT}), zero(MMat3{TT}), zero(MMat3{TT})
   @einsum nn[i,j] = n[α] * C[i,α,j,β] * n[β]
   @einsum nm[i,j] = n[α] * C[i,α,j,β] * m[β]
   nn⁻¹ = inv(nn)
   Du = 1/(2*π*r) * ( kron( (-S*b) ,m' ) + kron(nn⁻¹*(4*π*B + nm*S)*b,n') )
   return Mat3(Du)
end

# ========== Edge dislocation isotropic solid ==============

struct IsoEdgeDislocation3D{T} <: AbstractDislocation{T}
   λ::T
   μ::T
   b::T
   remove_singularity::Bool
end

IsoEdgeDislocation3D(λ, μ, b; remove_singularity = true) =
   IsoEdgeDislocation3D(λ, μ, b, remove_singularity)

function (Disl::IsoEdgeDislocation3D{T})(x) where T
   if Disl.remove_singularity && norm(x) < 1e-10
      return @SVector zeros(T, 3)
   end
   return eval_isoedge(Vec3(x), Disl.b, Disl.λ, Disl.μ)
end

function grad(Disl::IsoEdgeDislocation3D{T}, x) where T
   if Disl.remove_singularity && norm(x) < 1e-10
      return @SMatrix zeros(T, 3, 3)
   end
   return grad_isoedge(Vec3{T}(x), Disl.b, Disl.λ, Disl.μ)
end

"Isotropic CLE Edge dislocation"
function eval_isoedge(x::Vec3{T}, b::Real, λ::Real, μ::Real) where T
   u = zeros(T,3)
   # Compute Poisson ratio
   ν = λ/(2*(λ+μ))
   r² = dot(x[1:2],x[1:2])
	θ = angle.(x[1] + im*x[2])
	if θ < 0
		θ += 2π
	end
   u[1] = b/(2*π) * (θ + (x[1]*x[2])/(2*(1-ν) * r²))
   u[2] = -b/(2*π) * ( (1-2*ν)/(4*(1-ν)) * log(r²) - 2*x[2]^2 /(4*(1-ν)*r²))
   return Vec3(u)
end

"displacement gradient due to isotropic CLE Edge dislocation"
function grad_isoedge(x::Vec3{T}, b::Real, λ::Real, μ::Real) where T
   Du = zeros(T,3,3)
   # Compute Poisson ratio
   ν = λ/(2*(λ+μ))
   r² = dot(x[1:2],x[1:2])
   r⁴ = (r²).^2
   Du[1,1] = b/(2*π) * ( (2*ν-1)/(2*(1-ν)) * x[2] / r² -  x[1]^2 * x[2] / ((1-ν)*r⁴))
   Du[1,2] = b/(2*π) * ( (3-2*ν)/(2*(1-ν)) * x[1] / r² -  x[1]*x[2]^2 / ((1-ν)*r⁴))
   Du[2,1] = -b/(2*π) * ( x[1]/r² + (x[1]*(x[2].^2-x[1]^2))/(2*(1-ν)*r⁴))
   Du[2,2] = b/(2*π) * ( ν*x[2]/((1-ν)*r²) + (x[2]*(x[1]^2-x[2]^2))/(2*(1-ν)*r⁴))
   return Mat3(Du)
end

# ========== Screw dislocation isotropic solid ==============

struct IsoScrewDislocation3D{T} <: AbstractDislocation{T}
   b::T
   remove_singularity::Bool
end

IsoScrewDislocation3D(b; remove_singularity = true) =
   IsoScrewDislocation3D(b, remove_singularity)

function (Disl::IsoScrewDislocation3D{T})(x) where T
   if Disl.remove_singularity && norm(x) < 1e-10
      return @SVector zeros(T, 3)
   end
   return eval_isoscrew(Vec3(x), Disl.b)
end

function grad(Disl::IsoScrewDislocation3D{T}, x) where T
   if Disl.remove_singularity && norm(x) < 1e-10
      return @SMatrix zeros(T, 3, 3)
   end
   return grad_isoscrew(Vec3{T}(x), Disl.b)
end

"Isotropic CLE Screw dislocation"
function eval_isoscrew(x::Vec3{T}, b::Real) where T
   u = zeros(T,3)
	θ = angle( x[1] , x[2] )
	if θ < 0.
		θ += 2π
	end
   u[3] = b/(2*π) * θ
   return Vec3(u)
end

"displacement gradient due to isotropic CLE Screw dislocation"
function grad_isoscrew(x::Vec3{T}, b::Real) where T
   Du = zeros(T,3,3)
   r² = dot(x[1:2],x[1:2])
   Du[3,1] = b/(2*π) * ( -x[2] / r² )
   Du[3,2] =  b/(2*π) * ( x[1] / r² )
   return Mat3(Du)
end


text:
# NOTE: consider this experimental at best... you probably shouldn't use it

# ---------------------------------------------------------------------

type LiquidParams
  l::Int  # column dimensions l x w x h
  w::Int
  h::Int
  neuronType::DataType
  pctInhibitory::Float64
  decayRateDist::Distribution{Univariate,Continuous}
  λ::Float64  # used in probability of synapse connection
  pctInput::Float64
  pctOutput::Float64
  readout::Readout
  dt::Float64         # ms - simulation time step
  # τu::Float64         # ms - membrane potential decay period
  # τq::Float64          # ms - pulse decay period
  uDecayRate::Float64    # 1/ms - equal to 1/τu = membrane potential decay rate
  qDecayRate::Float64    # 1/ms - equal to 1/τq = pulse decay rate
  urest::Float64       # mV - resting membrane potential
  ufire::Float64      # mV - membrane potential immediately after firing a spike
  baseThreshold::Float64
end

function LiquidParams(; l::Int = 3,
                        w::Int = l,
                        h::Int = 5,
                        neuronType::DataType = SRMNeuron,
                        pctInhibitory::Float64 = 0.2,
                        decayRateDist::Distribution{Univariate,Continuous} = Uniform(0.8, 0.99),
                        λ::Float64 = 1.0,
                        pctInput::Float64 = 0.1,
                        pctOutput::Float64 = 0.4,
                        readout::Readout = FireReadout(),
                        dt::Float64 = 0.25,
                        # τu::Float64 = 25.0,
                        # τq::Float64 = 4.0,
                        uDecayRate::Float64 = 1.0 / 25.0,
                        qDecayRate::Float64 = 1.0 / 4.0,
                        urest::Float64 = 0.0,
                        ufire::Float64 = -0.4,
                        baseThreshold::Float64 = 2.0)
  LiquidParams(l, w, h, neuronType, pctInhibitory, decayRateDist, λ, pctInput, pctOutput, readout, dt, uDecayRate, qDecayRate, urest, ufire, baseThreshold)
end

Base.print(io::IO, l::LiquidParams) = print(io, "Params{$(l.l)x$(l.w)x$(l.h), pcti=$(l.pctInhibitory), λ=$(l.λ), pctout=$(l.pctOutput)}")
Base.show(io::IO, l::LiquidParams) = print(io, l)

# ---------------------------------------------------------------------

# maintains a list of GRFs which correspond to several spike trains for each input value
type LiquidInputs{T <: LiquidInput}
  inputs::Vector{T}   # one per input
end

function OnlineStats.fit!(inputs::LiquidInputs, x::VecF, dt::Float64)
  for i in 1:length(inputs.inputs)
    fit!(inputs.inputs[i], x[i], dt)
  end
end

# ---------------------------------------------------------------------

# methods to compute probability of synaptic connection
#   P(connection) = C * exp(-(D(a,b)/λ)^2)

const C_EE = 0.3
const C_EI = 0.2
const C_IE = 0.4
const C_II = 0.1
function C(n1::SpikingNeuron, n2::SpikingNeuron)
  n1.excitatory && (return n2.excitatory ? C_EE : C_EI)
  n2.excitatory ? C_IE : C_II
end

function distance(n1::SpikingNeuron, n2::SpikingNeuron)
  norm(n1.position - n2.position)
end

function probabilityOfConnection(n1::SpikingNeuron, n2::SpikingNeuron, λ::Float64)
  C(n1, n2) * exp(-((distance(n1, n2) / λ) ^ 2))
end

# TODO: make this a parameter
const UNIF_WEIGHT = Uniform(0.1, 1.0)

function weight(n::SpikingNeuron)
  (n.excitatory ? 1.0 : -1.0) * rand(UNIF_WEIGHT)
end

bound{T<:Real}(x::T, l::T, u::T) = max(l, min(x, u))

function delay(n::SpikingNeuron)
  sample(1:MAX_FUTURE)
end

# ---------------------------------------------------------------------

type Liquid{T}
  params::LiquidParams
  neurons::Vector{T}
  outputNeurons::Vector{T}
end

Base.print(io::IO, l::Liquid) = print(io, "Liquid{n=$(length(l.neurons)), nout=$(length(l.outputNeurons))}")
Base.show(io::IO, l::Liquid) = print(io, l)

function Liquid{T}(::Type{T}, params::LiquidParams)

  # create neurons in an (w x w x h) column
  # neurons = DiscreteLeakyIntegrateAndFireNeuron[]
  neurons = T[]
  for i in 1:params.l
    for j in 1:params.w
      for k in 1:params.h
        excitatory = rand() > params.pctInhibitory
        neuron = T([i,j,k], excitatory, params)
        push!(neurons, neuron)
      end
    end
  end

  # output from a random subset of the neuronal column
  outputNeurons = sample(neurons, params.pctOutput)

  # randomly connect the neurons
  for n1 in neurons
    for n2 in neurons
      if rand() <= probabilityOfConnection(n1, n2, params.λ)
        # synapse = DelaySynapse(n2, weight(n1), delay(n1))
        synapse = SRMSynapse(n2, weight(n1))
        push!(n1.synapses, synapse)
      end
    end
    println(n1, ", Synapses:")
    for s in n1.synapses println("     ", s) end
  end

  Liquid{T}(params, neurons, outputNeurons)
end

StatsBase.sample{T<:SpikingNeuron}(neurons::Vector{T}, pct::Float64) = sample(neurons, round(Int, pct * length(neurons)))

function OnlineStats.fit!(liquid::Liquid, dt::Float64)
  for neuron in liquid.neurons
    fit!(neuron, dt)
  end
  foreach(liquid.neurons, fire!)
end

# ---------------------------------------------------------------------

# manages the various layers and flow: input --> liquid --> readout --> readout model

type LiquidStateMachine <: OnlineStat
  nin::Int
  nout::Int
  liquid::Liquid
  inputs::LiquidInputs
  readout::Readout
  readoutModels::Vector{OnlineStat}
  n::Int
end

function LiquidStateMachine(params::LiquidParams, numInputs::Int, numOutputs::Int)
  # initialize liquid
  liquid = Liquid(params.neuronType, params)

  # create input structure
  wgt = ExponentialWeight(10000)
  variances = [Variance(wgt) for i in 1:numInputs]
  inputs = LiquidInputs(liquid, variances)

  # create readout models
  # TODO: make readout model more flexible... param
  # readoutModels = OnlineStat[OnlineFLS(length(liquid.outputNeurons), 0.00001, wgt) for i in 1:numOutputs]
  readoutModels = OnlineStat[Adagrad(length(liquid.outputNeurons); penalty = L2Penalty(0.00001)) for i in 1:numOutputs]

  LiquidStateMachine(numInputs, numOutputs, liquid, inputs, params.readout, readoutModels, 0)
end

# liquidState(lsm::LiquidStateMachine) = Float64[float(neuron.fired) for neuron in lsm.liquid.outputNeurons]
# OnlineStats.statenames(lsm::LiquidStateMachine) = [:liquidState, :nobs]
# OnlineStats.state(lsm::LiquidStateMachine) = Any[liquidState(lsm), nobs(lsm)]

function OnlineStats.fit!(lsm::LiquidStateMachine, x::VecF, y::VecF, dt::Float64 = lsm.liquid.params.dt)
  fit!(lsm.inputs, x, dt)   # update input neurons
  fit!(lsm.liquid, dt)     # update liquid state

  # update readout models
  # TODO: liquidState should be more flexible... multiple models, recent window averages, etc
  fit!(lsm.readout, lsm.liquid.outputNeurons)
  state = liquidState(lsm, lsm.readout)
  for (i,model) in enumerate(lsm.readoutModels)
    fit!(model, state, y[i])
  end

  lsm.n += 1
end

# given the current liquid state and readout model, predict the future
function StatsBase.predict(lsm::LiquidStateMachine)
  state = liquidState(lsm, lsm.readout)
  Float64[predict(model, state) for model in lsm.readoutModels]
end

# ---------------------------------------------------------------------

# ---------------------------------------------------------------------


text:
# Huber loss function

using LinearAlgebra

export HuberLoss

"""
**Huber loss**

    HuberLoss(ρ=1.0, μ=1.0)

Returns the function
```math
f(x) = \\begin{cases}
    \\tfrac{μ}{2}\\|x\\|^2 & \\text{if}\\ \\|x\\| ⩽ ρ \\\\
    ρμ(\\|x\\| - \\tfrac{ρ}{2}) & \\text{otherwise},
\\end{cases}
```
where `ρ` and `μ` are positive parameters.
"""
struct HuberLoss{R <: Real} <: ProximableFunction
    rho::R
    mu::R
    function HuberLoss{R}(rho::R, mu::R) where {R <: Real}
        if rho <= 0.0 || mu <= 0.0
            error("parameters rho and mu must be positive")
        else
            new(rho, mu)
        end
    end
end

is_convex(f::HuberLoss) = true
is_smooth(f::HuberLoss) = true

HuberLoss(rho::R=1.0, mu::R=1.0) where {R <: Real} = HuberLoss{R}(rho, mu)

function (f::HuberLoss)(x::AbstractArray{T}) where T <: Union{Real, Complex}
    normx = norm(x)
    if normx <= f.rho
        return (f.mu/2)*normx^2
    else
        return f.rho*f.mu*(normx-f.rho/2)
    end
end

function gradient!(y::AbstractArray{T}, f::HuberLoss, x::AbstractArray{T}) where T <: Union{Real, Complex}
    normx = norm(x)
    if normx <= f.rho
        y .= f.mu .* x
        v = (f.mu/2)*normx^2
    else
        y .= (f.mu*f.rho)/normx .* x
        v = f.rho*f.mu*(normx-f.rho/2)
    end
    return v
end

function prox!(y::AbstractArray{T}, f::HuberLoss, x::AbstractArray{T}, gamma::Real=1.0) where T <: Union{Real, Complex}
    normx = norm(x)
    mugam = f.mu*gamma
    scal = (1-min(mugam/(1+mugam), mugam*f.rho/(normx)))
    for k in eachindex(y)
        y[k] = scal*x[k]
    end
    normy = scal*normx
    if normy <= f.rho
        return (f.mu/2)*normy^2
    else
        return f.rho*f.mu*(normy-f.rho/2)
    end
end

fun_name(f::HuberLoss) = "Huber loss"
fun_dom(f::HuberLoss) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::HuberLoss) = "x ↦ (μ/2)||x||² if ||x||⩽ρ, μρ(||x||-ρ/2) otherwise"
fun_params(f::HuberLoss) = string("ρ = $(f.rho), μ = $(f.mu)")

function prox_naive(f::HuberLoss, x::AbstractArray{T}, gamma::Real=1.0) where T <: Union{Real, Complex}
    y = (1-min(f.mu*gamma/(1+f.mu*gamma), f.mu*gamma*f.rho/(norm(x))))*x
    if norm(y) <= f.rho
        return y, (f.mu/2)*norm(y)^2
    else
        return y, f.rho*f.mu*(norm(y)-f.rho/2)
    end
end


text:
using NLsolve
using DiffEqBase
using Test
using NLSolversBase
using LineSearches
using LinearAlgebra
import Base.convert
using SparseArrays
using Printf
using IterativeSolvers
using Random
add_jl(x) = endswith(x, ".jl") ? x : x*".jl"

mutable struct WrappedArray{T,N} <: DEDataArray{T,N}
    x::Array{T,N}
end

Base.reshape(A::WrappedArray, dims::Int...) = WrappedArray(reshape(A.x, dims...))
Base.convert(A::Type{WrappedArray{T,N}}, B::Array{T,N}) where {T,N} = WrappedArray(copy(B))

if length(ARGS) > 0
    tests = map(add_jl, ARGS)
else
    tests = ["2by2.jl",
             "linsolve.jl",
             "minpack.jl",
             "singular.jl",
             "finite_difference.jl",
             "iface.jl",
             "incomplete.jl",
             "already_converged.jl",
             "autodiff.jl",
             "josephy.jl",
             "difficult_mcp.jl",
             "sparse.jl",
             "throws.jl",
             "f_g_counts.jl",
             "no_linesearch.jl",
             "abstractarray.jl",
             "complex.jl",
             "interface/caches.jl",
             "fixedpoint/fixedpoint.jl"]
end

println("Running tests:")

for test in tests
    include(test)
end


text:
using Formatting
using JSON, ArgParse

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--config-file"
            help = "Configuration file."
            arg_type = String
            required = true

        "--forcing-file"
            help = "Annual forcing file. It should contain: TAUX, TAUY, SWFLX, NSWFLX, VSFLX"
            arg_type = String
            required = true

        "--nudging-timescale-day"
            help = "Nudging timescale used to replace weak restoring timescale when finding Q-flux in days. Default is 5 days."
            arg_type = Float64
            default = 5.0

        "--stop-n"
            help = "Core of the model."
            arg_type = Int64

        "--time-unit"
            help = "."
            arg_type = String

    end

    return parse_args(s)
end

parsed = parse_commandline()

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
function _checksquare(A::AbstractMatrix)
    size(A,1) != size(A,2) && throw(ArgumentError("Matrix should be square"))
end

function _check_selfloops(A::AbstractMatrix)
    s = min(size(A)...)
    any(!iszero, A[i,i] for i=1:s) && throw(ArgumentError("Self loops are not allowed"))
end

# vector that maps edge i->j to edge j->i in A.nzval
function inverse_edge_idx(A::SparseMatrixCSC)
    issymmetric(A) || throw(ArgumentError("Matrix should be symmetric"))
    n = size(A, 2)
    X = zeros(Int, nnz(A))
    rowcnt = zeros(Int, n)
    for j in 1:n
        for k in nzrange(A, j)
            r = A.rowval[k]
            rowcnt[r] += 1
            if r > j
                z = nzrange(A, r)[rowcnt[r]]
                X[k], X[z] = z, k
            end
        end
    end
    X
end

# return the index in `A.nzval` of element (`i`,`j`) if it's nonzero,
# otherwise raise an error
#
# copied from julia/stdlib/v1.7/SparseArrays/src/sparsematrix.jl:2158
function nzindex(A::SparseMatrixCSC, i::Integer, j::Integer)
    if !(1 <= i <= size(A, 1) && 1 <= j <= size(A, 2)); throw(BoundsError()); end
    k1 = Int(A.colptr[j])
    k2 = Int(A.colptr[j+1]-1)
    (k1 > k2) && throw(ArgumentError("Matrix element ($i,$j) is zero"))
    k = searchsortedfirst(rowvals(A), i, k1, k2, Base.Order.Forward)
    if k > k2 || A.rowval[k] != i
        throw(ArgumentError("Matrix element ($i,$j) is zero"))
    end
    return k
    # ???? Propagate errors to higher-level functions? E.g. do not throw errors
    #  here but return some Int, then from `edge_idx` print something like
    #  "Graph does not contain edge (i,j)" ?????
end

# return the indices `(i,j)` of the `k`th element in `A.nzval`
function nzindex(A::SparseMatrixCSC, k::Integer)
    j = searchsortedfirst(A.colptr, k+1) - 1
    i = rowvals(A)[k]
    i, j
end

# Void type used to fill sparse arrays w/o using memory for the nonzero values,
#  only for their positions
struct NullNumber <: Number; end
Base.zero(::NullNumber) = false
Base.zero(::Type{NullNumber}) = false
Base.iszero(::NullNumber) = false


text:
"""
    abstract type to include all layers
"""
abstract type Layer end

export Layer

### FCLayer

@doc raw"""
    FCLayer(channels=0, actFun=:noAct, [layerInput = nothing; keepProb = 1.0])

Fully-connected layer (equivalent to Dense in TensorFlow etc.)

# Arguments

- `channels` := (`Integer`) is the number of nodes in the layer

- `actFun` := (`Symbol`) is the activation function of this layer

- `layerInput` := (`Layer` or `Array`) the input of this array (optional don't need to assign it)

- `keepProb` := (`AbstractFloat`) the keep probability (1 - prob of the dropout rate)

---------

# Summary

    mutable struct FCLayer <: Layer

-------

# Fields

- `channels::Integer` := is the number of nodes in the layer

- `actFun::Symbol` := the activation function of this layer

- `inputS::Tuple{Integer, Integer}` := input size of the layer, of the shape (channels of the previous layer, size of mini-batch)

- `outputS::Tuple{Integer, Integer}` := output size of the layer, of the shape (channels of this layer, size of mini-batch)

- `keepProb::AbstractFloat` := the keep probability (rate) of the drop-out operation `<1.0`

- `W::Array{T,2} where {T}` := the scaling parameters of this layer `W * X`, of the shape (channels of this layer, channels of the previous layer)

- `B::Array{T,2} where {T}` := the bias of this layer `W * X .+ B`, of the shape (channels of this layer, 1)

- `dW::Array{T,2} where {T}` := the derivative of the loss function to the W parameters $\frac{dJ}{dW}$

- `dB::Array{T,2} where {T}` := the derivative of the loss function to the B parameters $\frac{dJ}{dB}$

- `forwCount::Integer` := forward propagation counter

- `backCount::Integer` := backward propagation counter

- `updateCount::Integer` := update parameters counter

- `prevLayer::L where {L<:Union{Layer,Nothing}}` := the previous layer which is
the input of this layer

- `nextLayers::Array{Layer,1}` := An array of the next `layer`(s)

---------

# Supertype Hierarchy

    FCLayer <: Layer <: Any

---------

# Examples

```julia
X_Input = Input(X_train)
X = FCLayer(20, :relu)(X_Input)

```

In the previous example the variable `X_Input` is a pointer to the `Input` layer, and `X` is an pointer to the `FCLayer(20, :relu)` layer.
**Note** that the layer instance can be used as a connecting function.

"""
mutable struct FCLayer <: Layer
    channels::Integer
    actFun::Symbol

    inputS::Tuple{Integer, Integer}
    outputS::Tuple{Integer, Integer}

    """
        drop-out keeping node probability
    """
    keepProb::AbstractFloat
    W::Array{T,2} where {T}
    B::Array{T,2} where {T}
    dW::Array{T,2} where {T}
    dB::Array{T,2} where {T}
    # ### adding Z & A place holder for recursive calling
    # ### and a counter for how many it was called
    # Z::Array{T,2} where {T}
    # dA::Array{T,2} where {T}
    # A::Array{T,2} where {T}

    V::Dict{Symbol,Array{T,2} where {T,N}}
    S::Dict{Symbol,Array{T,2} where {T,N}}

    forwCount::Integer
    backCount::Integer
    updateCount::Integer
    """
        pointer to previous layer
    """
    prevLayer::L where {L<:Union{Layer,Nothing}}
    nextLayers::Array{Layer,1}
    function FCLayer(channels=0, actFun=:noAct, layerInput = nothing; keepProb = 1.0)
        # W, B
        if isa(layerInput, Layer)
            T = eltype(layerInput.W)
            nl = channels
            nl_1 = size(layerInput.W)[1]
            prevLayer = layerInput
        elseif isa(layerInput, Array)
            T = eltype(layerInput)
            nl = channels
            nl_1 = size(layerInput)[1]
            prevLayer = nothing
        else
            T = Any
            nl = 0
            nl_1 = 0
            prevLayer = nothing
        end
        new(
            channels,
            actFun,
            (0,0), #inputS
            (0,0), #outputS
            keepProb,
            Matrix{T}(undef, 0, 0),
            Matrix{T}(undef, 0, 0),
            Matrix{T}(undef, 0, 0),
            Matrix{T}(undef, 0, 0),
            # Matrix{T}(undef, 0, 0),
            # Matrix{T}(undef, 0, 0),
            # Matrix{T}(undef, 0, 0),
            Dict(:dw => Matrix{T}(undef, 0, 0), :db => Matrix{T}(undef, 0, 0)),
            Dict(:dw => Matrix{T}(undef, 0, 0), :db => Matrix{T}(undef, 0, 0)),
            0,
            0,
            0,
            prevLayer,
            Array{Layer,1}(undef,0)
        )#
    end #FCLayer
end #struct FCLayer

Base.show(io::IO, l::FCLayer) = print(io, "FCLayer($(l.channels), $(l.actFun), $(l.inputS), $(l.outputS), ...)")

export FCLayer

### Mulit-Input Layer MILayer

export MILayer

abstract type MILayer <: Layer end

### AddLayer

@doc raw"""
    AddLayer(; [channels = 0])

Layer performs and addition of multiple previous layers

# Arguments

- `channels` := (`Integer`) number of channels/nodes of this array which equals to the same of the previous layer(s)

---------

# Summary

    mutable struct AddLayer <: MILayer

# Fields

- `channels::Integer` := is the number of nodes or `channels` in the layer

- `inputS::Tuple` := input size of the layer

- `outputS::Tuple` := output size of the layer

- `forwCount::Integer` := forward propagation counter

- `backCount::Integer` := backward propagation counter

- `updateCount::Integer` := update parameters counter

- `nextLayers::Array{Layer,1}` := An array of the next `layer`(s)

- `prevLayer::Array{Layer,1}` := An array of the previous `layer`(s) to be added

---------

# Supertype Hierarchy

    AddLayer <: MILayer <: Layer <: Any

-------

# Examples

```julia
XIn1 = Input(X_train)
X1 = FCLayer(10, :relu)(XIn1)
XIn2 = Input(X_train)
X2 = FCLayer(10, :tanh)(XIn2)

Xa = AddLayer()([X1,X2])
```

"""
mutable struct AddLayer <: MILayer
    channels::Integer

    inputS::Tuple
    outputS::Tuple

    # A::Array{T,N} where {T,N}
    # dA::Array{T,N} where {T,N}

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    nextLayers::Array{Layer,1}
    prevLayer::Array{Layer,1}
    function AddLayer(; channels = 0)
        # channels = l1.channels
        # T = eltype(l1.W)
        new(
            channels,
            (0,), #inputS
            (0,), #outputS
            # Matrix{Nothing}(undef, 0, 0),
            # Matrix{Nothing}(undef, 0, 0),
            0,
            0,
            0,
            Array{Layer,1}(undef,0), #nextLayers
            Array{Layer,1}(undef,0), #prevLayer
            )
    end #function AddLayer
end

Base.show(io::IO, l::AddLayer) = print(io, "AddLayer($(l.channels), $(l.inputS), $(l.outputS), ...)")

export AddLayer

### ConcatLayer

export ConcatLayer

@doc raw"""
    ConcatLayer(; channels = 0)

Perform concatenation of group of previous `Layer`s

# Summary

    mutable struct ConcatLayer <: MILayer

# Fields

    channels    :: Integer
    inputS      :: Tuple
    outputS     :: Tuple
    forwCount   :: Integer
    backCount   :: Integer
    updateCount :: Integer
    nextLayers  :: Array{Layer,1}
    prevLayer   :: Array{Layer,1}

# Supertype Hierarchy

    ConcatLayer <: MILayer <: Layer <: Any
"""
mutable struct ConcatLayer <: MILayer
    channels::Integer

    inputS::Tuple
    outputS::Tuple

    # A::Array{T,N} where {T,N}
    # dA::Array{T,N} where {T,N}

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    nextLayers::Array{Layer,1}
    prevLayer::Array{Layer,1}

    LSlice::Dict{Layer,UnitRange{Integer}}
    function ConcatLayer(; channels = 0)
        # channels = l1.channels
        # T = eltype(l1.W)
        new(
            channels,
            (0,), #inputS
            (0,), #outputS
            # Matrix{Nothing}(undef, 0, 0),
            # Matrix{Nothing}(undef, 0, 0),
            0,
            0,
            0,
            Array{Layer,1}(undef,0), #nextLayers
            Array{Layer,1}(undef,0), #prevLayer
            Dict{Layer,UnitRange{Integer}}(),
            )
    end #function ConcatLayer
end #mutable struct ConcatLayer

Base.show(io::IO, l::ConcatLayer) = print(io, "ConcatLayer($(l.channels), $(l.inputS), $(l.outputS), ...)")

### Activation

@doc raw"""
    Activation(actFun)

# Arguments

- `actFun::Symbol` := the activation function of this layer

---------

# Summary

mutable struct Activation <: Layer

--------

# Fields

- `actFun::Symbol` := the activation function of this layer

- `channels::Integer` := is the number of nodes or `channels` in the layer

- `inputS::Tuple` := input size of the layer

- `outputS::Tuple` := output size of the layer

- `forwCount::Integer` := forward propagation counter

- `backCount::Integer` := backward propagation counter

- `updateCount::Integer` := update parameters counter

- `nextLayers::Array{Layer,1}` := An array of the next `layer`(s)

- `prevLayer::Array{Layer,1}` := An array of the previous `layer`(s) to be added

---------

# Supertype Hierarchy

    Activation <: Layer <: An

---------

# Examples

```julia
X_Input = Input(X_train)
X = FCLayer(10, :noAct)(X_Input)
X = Activation(:relu)(X)
```

"""
mutable struct Activation <: Layer
    actFun::Symbol
    channels::Integer

    inputS::Tuple
    outputS::Tuple

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    nextLayers::Array{Layer,1}
    prevLayer::L where {L<:Union{Layer,Nothing}}

    function Activation(actFun=:relu)
        new(
            actFun,
            0, #channels
            (0,), #inputS
            (0,), #outputS
            0,
            0,
            0,
            Array{Layer,1}(undef,0),
            nothing,
            )
    end #function Activation
end #mutable struct Activation

Base.show(io::IO, l::Activation) = print(io, "Activation($(l.channels), $(l.actFun), $(l.inputS), $(l.outputS), ...)")

export Activation

### Input

@doc raw"""
    Input(X_shape::Tuple)

`Input` `Layer` that is used as a pointer to the input array(s).

# Arguments

- `X_shape::Tuple` := shape of the input Array

-------

# Summary

    mutable struct Input <: Layer

# Fields

- `channels::Integer` := is the number of nodes or `channels` in the layer

- `inputS::Tuple` := input size of the layer

- `outputS::Tuple` := output size of the layer

- `forwCount::Integer` := forward propagation counter

- `backCount::Integer` := backward propagation counter

- `updateCount::Integer` := update parameters counter

- `nextLayers::Array{Layer,1}` := An array of the next `layer`(s)

- `prevLayer::Array{Layer,1}` := An array of the previous `layer`(s) to be added

------

# Supertype Hierarchy

    Input <: Layer <: Any

--------

# Examples

```julia
X_Input = Input(size(X_train))
X = FCLayer(10, :relu)(X_Input)
```

It is possible to use the Array instead of its size `NumNN` will take care of the rest

```julia
X_Input = Input(X_train)
X = FCLayer(10, :relu)(X_Input)
```

"""
mutable struct Input <: Layer
    channels::Integer

    inputS::Tuple
    outputS::Tuple

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    nextLayers::Array{Layer,1}
    prevLayer::L where {L<:Union{Layer,Nothing}}

    function Input(X_shape::Tuple)
        N = length(X_shape)
        if N==2
            channels = X_shape[1]
        elseif N==3
            channels = X_shape[2]
        elseif N==4
            channels = X_shape[3]
        elseif N==5
            channels = X_shape[4]
        end

        new(
            channels,
            X_shape, #inputS
            X_shape, #outputS
            0,
            0,
            0,
            Array{Layer,1}(undef,0),
            nothing,
            )
    end #function Layer
end #mutable struct Input

function Input(X::Array{T,N}) where {T,N}
    X_shape = size(X)
    Input(X_shape)
end #function Input(X::Array{T,N}) where {T,N}

Base.show(io::IO, l::Input) = print(io, "Input($(l.channels), $(l.inputS), $(l.outputS), ...)")

export Input

### BatchNorm

@doc raw"""
    BatchNorm(;dim=1, ϵ=1e-10)

Batch Normalization `Layer` that is used ot normalize across the dimensions specified by the argument `dim`.

# Arguments

- `dim::Integer` := is the dimension to normalize across

- `ϵ::AbstractFloat` := is a backup constant that is used to prevent from division on zero when ``σ^2`` is zero

---------

# Summary

    mutable struct BatchNorm <: Layer

-------

# Fields

- `channels::Integer` := is the number of nodes in the layer

- `inputS::Tuple{Integer, Integer}` := input size of the layer, of the shape (channels of the previous layer, size of mini-batch)

- `outputS::Tuple{Integer, Integer}` := output size of the layer, of the shape (channels of this layer, size of mini-batch)

- `dim::Integer` := the dimension to normalize across

- `ϵ::AbstractFloat` := backup constant to protect from dividing on zero when ``σ^2 = 0``

- `W::Array{T,2} where {T}` := the scaling parameters of this layer `W * X`, same shape of the mean `μ`

- `B::Array{T,2} where {T}` := the bias of this layer `W * X .+ B`, same shape of the variance ``σ^2``

- `dW::Array{T,2} where {T}` := the derivative of the loss function to the W parameters $\frac{dJ}{dW}$

- `dB::Array{T,2} where {T}` := the derivative of the loss function to the B parameters $\frac{dJ}{dB}$

- `forwCount::Integer` := forward propagation counter

- `backCount::Integer` := backward propagation counter

- `updateCount::Integer` := update parameters counter

- `prevLayer::L where {L<:Union{Layer,Nothing}}` := the previous layer which is
the input of this layer

- `nextLayers::Array{Layer,1}` := An array of the next `layer`(s)

---------

# Supertype Hierarchy

    BatchNorm <: Layer <: Any

---------

# Examples

```julia
X_train = rand(14,14,3,32) #input of shape `14×14` with channels of `3` and mini-batch size `32`

X_Input = Input(X_train)
X = Conv2D(10, (3,3))(X_Input)
X = BatchNorm(dim=3) #to normalize across the channels dimension
X = Activation(:relu)
```

```julia
X_train = rand(128,5,32) #input of shape `128` with channels of `5` and mini-batch size `32`

X_Input = Input(X_train)
X = Conv1D(10, 5)(X_Input)
X = BatchNorm(dim=2) #to normalize across the channels dimension
X = Activation(:relu)
``

```julia
X_train = rand(64*64,32) #input of shape `64*64` and mini-batch size `32`

X_Input = Input(X_train)
X = FCLayer(10, :noAct)(X_Input)
X = BatchNorm(dim=1) #to normalize across the features dimension
X = Activation(:relu)
````

"""
mutable struct BatchNorm <: Layer
    channels::Integer

    inputS::Tuple
    outputS::Tuple

    dim::Integer

    ϵ::AbstractFloat

    W::Array{T, N} where {T,N}
    dW::Array{T, N} where {T,N}

    B::Array{T, N} where {T,N}
    dB::Array{T, N} where {T,N}

    V::Dict{Symbol,Array}
    S::Dict{Symbol,Array}

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    nextLayers::Array{Layer,1}
    prevLayer::L where {L<:Union{Layer,Nothing}}

    function BatchNorm(;dim=1, ϵ=1e-7)
        new(
            0, #channels
            (0,), #inputS
            (0,), #outputS
            dim,
            ϵ,
            Array{Any,1}(undef,0), #W
            Array{Any,1}(undef,0), #dW
            Array{Any,1}(undef,0), #B
            Array{Any,1}(undef,0), #dB
            Dict(:dw=>Array{Any,1}(undef,0),
                 :db=>Array{Any,1}(undef,0)), #V
            Dict(:dw=>Array{Any,1}(undef,0),
                 :db=>Array{Any,1}(undef,0)), #S
            0, #forwCount
            0, #backCount
            0, #updateCount
            Array{Layer,1}(undef,0), #nextLayers
            nothing, #prevLayer
            )
    end #function BatchNorm

end #mutable struct BatchNorm

Base.show(io::IO, l::BatchNorm) = print(io, "BatchNorm($(l.channels), dim = $(l.dim), $(l.inputS), $(l.outputS), ...)")

export BatchNorm

### Flatten

export Flatten

@doc raw"""
    Flatten()

Flatten the input into 2D `Array`

# Summary

    mutable struct Flatten <: Layer

# Fields

    channels    :: Integer
    inputS      :: Tuple
    outputS     :: Tuple
    forwCount   :: Integer
    backCount   :: Integer
    updateCount :: Integer
    nextLayers  :: Array{Layer,1}
    prevLayer   :: Union{Nothing, Layer}

# Supertype Hierarchy

    Flatten <: Layer <: Any
"""
mutable struct Flatten <: Layer
    channels::Integer

    inputS::Tuple
    outputS::Tuple

    forwCount::Integer
    backCount::Integer
    updateCount::Integer

    nextLayers::Array{Layer,1}
    prevLayer::L where {L<:Union{Layer,Nothing}}

    function Flatten()

        return new(
            0, #channels
            (0,), #inputS
            (0,), #outputS
            0, #forwCount
            0, #backCount
            0, #updateCount
            Array{Layer,1}(undef,0), #nextLayers
            nothing,
        )
    end

end

Base.show(io::IO, l::Flatten) = print(io, "Flatten($(l.channels), $(l.inputS), $(l.outputS), ...)")

### Model

@doc raw"""

    function Model(
        X,
        Y,
        inLayer::Layer,
        outLayer::Layer,
        α;
        optimizer = :gds,
        β1 = 0.9,
        β2 = 0.999,
        ϵAdam = 1e-8,
        regulization = 0,
        λ = 1.0,
        lossFun = :categoricalCrossentropy,
        paramsDtype::DataType = Float64,
    )

# Summary

    mutable struct Model <: Any

# Fields

    inLayer      :: Layer
    outLayer     :: Layer
    lossFun      :: Symbol
    paramsDtype  :: DataType
    regulization :: Integer
    λ            :: AbstractFloat
    α            :: AbstractFloat
    optimizer    :: Symbol
    ϵAdam        :: AbstractFloat
    β1           :: AbstractFloat
    β2           :: AbstractFloat
"""
mutable struct Model
    # layers::Array{Layer,1}
    inLayer::Layer
    outLayer::Layer
    lossFun::Symbol

    paramsDtype::DataType

    """
        regulization type
            0 : means no regulization
            1 : L1 regulization
            2 : L2 regulization
    """
    regulization::Integer

    """
        regulization constant
    """
    λ::AbstractFloat

    """
        learning rate
    """
    α::AbstractFloat

    optimizer::Symbol
    ϵAdam::AbstractFloat
    β1::AbstractFloat
    β2::AbstractFloat

    function Model(
        X,
        Y,
        inLayer::Layer,
        outLayer::Layer,
        α;
        optimizer = :gds,
        β1 = 0.9,
        β2 = 0.999,
        ϵAdam = 1e-8,
        regulization = 0,
        λ = 1.0,
        lossFun = :categoricalCrossentropy,
        paramsDtype::DataType = eltype(X),
    )

        deepInitWB!(outLayer; dtype = paramsDtype)
        resetCount!(outLayer, :forwCount)
        deepInitVS!(outLayer, optimizer)
        resetCount!(outLayer, :forwCount)
        @assert regulization in [0, 1, 2]
        return new(
            inLayer,
            outLayer,
            lossFun,
            paramsDtype,
            regulization,
            λ,
            α,
            optimizer,
            ϵAdam,
            β1,
            β2,
        )
    end #inner-constructor

end #Model

export Model


text:
using LinearAlgebra

"""
    angularWaveVectorComponent(n, L)

Returns the angular wave vector component for a rectangular room given the
mode index and room length along the required component cartesian direction, in
meters.
# Example
```jldoctest
julia> nx = 7
7

julia> Lx = 5
5

julia> kx = angularWaveVectorComponent(nx, Lx)
4.39822971502571
```
"""
function angularWaveVectorComponent(n::Integer, L::Real)
    return n * π / L
end

"""
    angularWaveVector(nx, ny, nz, Lx, Ly, Lz)

Returns the angular wave vector of a rectangular room mode. Mode numbers along
the cartesian directions are nx ny and nz, the room sizes along the cartesian
directions are Lx, Ly and Lz meters.
# Example
```jldoctest
julia> angularWaveVector(1, 4, 7, 5.0, 4.0, 3.0)
3-element Array{Float64,1}:
 0.6283185307179586
 3.141592653589793
 7.330382858376184
```
"""
function angularWaveVector(
    nx::Integer, ny::Integer, nz::Integer,
    Lx::Real, Ly::Real, Lz::Real
    )

    return [
        angularWaveVectorComponent(nx, Lx);
        angularWaveVectorComponent(ny, Ly);
        angularWaveVectorComponent(nz, Lz)
        ]

end

"""
    modeFrequency(nx, ny, nz, Lx, Ly, Lz, c)

Returns the modal frequency of a mode characterized by indeces nx, ny and nz
along the cartesian directions for a rectangular room of sizes Lx, Ly and Lz
meters along the cartesian directions. c is the speed of sound in air, by
default 343.0 meters per second (speed at standard conditions).
# Example
```jldoctest
julia> modeFrequency(1, 4, 7, 5.0, 4.0, 3.0)
436.7174156260672

julia> modeFrequency(1, 4, 7, 5.0, 4.0, 3.0, 340.0)
432.89772977511035
```
"""
function modeFrequency(
    nx::Integer, ny::Integer, nz::Integer,
    Lx::Real, Ly::Real, Lz::Real,
    c::Real = 343.0
    )

    return norm(angularWaveVector(nx, ny, nz, Lx, Ly, Lz)) * c / 2π

end

"""
    mode(x, y, z, nx, ny, nz, Lx, Ly, Lz)

Returns the modal shape at coordinates in x, y, z for mode characterized by
indeces nx, ny and nz along the cartesian directions for a rectangular room of
sizes Lx, Ly and Lz meters along the cartesian directions. Result returned
normalized between 1 and -1, arbitrary units.
# Example
```jldoctest
julia> mode(2.5, 2.0, 1.5, 1, 4, 7, 5.0, 4.0, 3.0)
0.0

julia> mode.([3.0; 2.5], [1.0; 2.0], [0.75; 1.5], 1, 4, 7, 5.0, 4.0, 3.0)
2-element Array{Float64,1}:
 0.2185080122244105
 0.0
```
"""
function mode(
    x::Real, y::Real, z::Real,
    nx::Integer, ny::Integer, nz::Integer,
    Lx::Real, Ly::Real, Lz::Real
    )

return cospi(nx * x / Lx) * cospi(ny * y / Ly) * cospi(nz * z / Lz)

end

"""
    meshGrid(x, y, z)

Returns the 3-D grid coordinates defined by the vectors x, y, and z.
The grid is represented by three matrices of size length(y) by length(x) by
length(z). This is similar to MATLAB's meshgrid.
# Example
```jldoctest
julia> x = [1; 2]
2-element Array{Int64,1}:
 1
 2

julia> y = [3; 4]
2-element Array{Int64,1}:
 3
 4

julia> z = [5; 6]
2-element Array{Int64,1}:
 5
 6

julia> X, Y, Z = meshGrid(x, y, z)
([1 2; 1 2]

[1 2; 1 2], [3 3; 4 4]

[3 3; 4 4], [5 5; 5 5]

[6 6; 6 6])
```
"""
function meshGrid(x::AbstractArray, y::AbstractArray, z::AbstractArray)

    X = zeros(eltype(x), length(y), length(x), length(z))
    fX(i, j, k) = X[i, j, k] = x[j]
    [fX(i, j, k) for i in 1:length(y), j in 1:length(x), k in 1:length(z)]

    Y = zeros(eltype(y), length(y), length(x), length(z))
    fY(i, j, k) = Y[i, j, k] = y[i]
    [fY(i, j, k) for i in 1:length(y), j in 1:length(x), k in 1:length(z)]

    Z = zeros(eltype(z), length(y), length(x), length(z))
    fZ(i, j, k) = Z[i, j, k] = z[k]
    [fZ(i, j, k) for i in 1:length(y), j in 1:length(x), k in 1:length(z)]

    return X, Y, Z

end

"""
    indexGrid(Nx, Ny, Nz)

Returns a 3-D grid up to indeces Nx, Ny and Nz along the 3 cartesian
directions. The grid is represented by three matrices of size Ny by Nx by Nz.
# Example
```jldoctest
julia> Gx, Gy, Gz = indexGrid(3, 4, 5)
([1 2 3; 1 2 3; 1 2 3; 1 2 3]

[1 2 3; 1 2 3; 1 2 3; 1 2 3]

[1 2 3; 1 2 3; 1 2 3; 1 2 3]

[1 2 3; 1 2 3; 1 2 3; 1 2 3]

[1 2 3; 1 2 3; 1 2 3; 1 2 3], [1 1 1; 2 2 2; 3 3 3; 4 4 4]

[1 1 1; 2 2 2; 3 3 3; 4 4 4]

[1 1 1; 2 2 2; 3 3 3; 4 4 4]

[1 1 1; 2 2 2; 3 3 3; 4 4 4]

[1 1 1; 2 2 2; 3 3 3; 4 4 4], [1 1 1; 1 1 1; 1 1 1; 1 1 1]

[2 2 2; 2 2 2; 2 2 2; 2 2 2]

[3 3 3; 3 3 3; 3 3 3; 3 3 3]

[4 4 4; 4 4 4; 4 4 4; 4 4 4]

[5 5 5; 5 5 5; 5 5 5; 5 5 5])

```
"""
function indexGrid(Nx::Integer, Ny::Integer, Nz::Integer)
    return meshGrid(0:Nx, 0:Ny, 0:Nz)
end

"""
    roomAxes(rx, ry, rz, Lx, Ly, Lz)

Returns the 1-D axes coordinates along the sides of a rectangular room with
reference to a corner. The room has sizes Lx, Ly and Lz meters along the
cartesian directions, while rx, ry and rz express the step of the axes along
each cartesian direction.
"""
function roomAxes(
    rx::Real, ry::Real, rz::Real,
    Lx::Real, Ly::Real, Lz::Real
    )

    x = range(0.0, stop = Lx, length = round(Integer, 1 + Lx / rx))
    y = range(0.0, stop = Ly, length = round(Integer, 1 + Ly / ry))
    z = range(0.0, stop = Lz, length = round(Integer, 1 + Lz / rz))

    return x, y, z

end

"""
    roomGrid(rx, ry, rz, Lx, Ly, Lz)

Returns the 3-D grid coordinates inside a rectangular room with reference to a
corner. The room has sizes Lx, Ly and Lz meters along the cartesian directions,
while rx, ry and rz express the step of the grid along each cartesian direction.
"""
function roomGrid(
    rx::Real, ry::Real, rz::Real,
    Lx::Real, Ly::Real, Lz::Real
    )

    x, y, z = roomAxes(rx, ry, rz, Lx, Ly, Lz)

    return meshGrid(x, y, z)

end

max_stars_repo_name: CrocoDuckoDucks/RectangularRoom

max_issues_repo_name: CrocoDuckoDucks/RectangularRoom

max_forks_repo_name: CrocoDuckoDucks/RectangularRoom


text:
half = Float16(4.71);
float = Float32(4.71);
double = Float64(4.71);

w1(x) = x^3 - convert(typeof(x),6.0)*x^2 + convert(typeof(x),3.0)*x - 0.149;
w2(x) = ((x-convert(typeof(x),6.0))*x + convert(typeof(x),3.0)) * x -0.149;
w1h = ( w1(half) - -14.636489) /  -14.636489;
w1f = ( w1(float) - -14.636489) /  -14.636489;
w1d = ( w1(double) - -14.636489) /  -14.636489;
w2h = ( w1(half) - -14.636489) /  -14.636489;
w2f = ( w1(float) - -14.636489) /  -14.636489;
w2d = ( w1(double) - -14.636489) /  -14.636489;
print(w1(half)); print("\twynik half\n");

print(w1h); print("\t<-w1\n");
print(w2h); print("\n");
print(w1f); print("\t<-w1\n");
print(w2f); print("\n");
print(w1d); print("\t<-w1\n");
print(w2d); print("\n");
print("aletypv23\n");

# cwiczenia

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
using Random
using GymSpaces: Box, Discrete

struct HotterColderEnv <: AbstractEnv
    range
    bounds
    action_space
    observation_space
    number
    guess_count
    guess_max
    observation
    seed::MersenneTwister
end

function HotterColderEnv()
    seed = MersenneTwister()
    range = 1000
    bounds = 2000

    action_space = Box([-bounds], [bounds], Float32)
    observation_space = Discrete(4)

    number = 2range * rand(Float32) - range
    guess_count = 0
    guess_max = 200
    observation = 0
    HotterColderEnv(range, bounds, action_space, observation_space, number, guess_count,
                    guess_max, observation, seed)
end

seed!(env::HotterColderEnv) = (env.seed = MersenneTwister())
seed!(env::HotterColderEnv, int) = (env.seed = MersenneTwister(int))

function reset!(env::HotterColderEnv)
    env.number = 2range * rand(Float32) - range
    env.guess_count = 0
    env.observation = 0
    return env.observation
end

function step!(env::HotterColderEnv, action)
    @assert action ∈ env.action_space

    if action < env.number
        env.observation = 1
    elseif action == env.number
        env.observation = 2
    elseif action == env.number
        env.observation = 3
    end

    reward = ((min(action, env.number) + env.bounds) / (max(action, env.number) + env.bounds)) ^ 2
    env.guess_count += 1
    done = env.guess_count ≥ env.guess_max

    return env.observation, reward[1], done Dict(:number => env.number, :guess => env.guess_count)
end


text:
using Test
using YaoBlocks.AD
using YaoBlocks
using YaoArrayRegister

@testset "NoParams" begin
    c = chain(put(3, 2 => Ry(0.2)), put(3, 1 => Rx(0.5)))
    np = chain(NoParams(put(3, 2 => Ry(0.2))), put(3, 1 => Rx(0.5)))
    reg = rand_state(3)
    @test content(np[1]) == put(3, 2 => Ry(0.2))
    @test parameters(np) == [0.5]
    @test nparameters(np) == 1
    @test getiparams(np[1]) == ()
    @test niparams(np[1]) == 0
    @test mat(c) ≈ mat(np)
    @test apply!(copy(reg), c) ≈ apply!(copy(reg), np)
    @test np' == chain(put(3, 1 => Rx(-0.5)), NoParams(put(3, 2 => Ry(-0.2))))
    @test dispatch!(np, [0.3]) ==
          chain(NoParams(put(3, 2 => Ry(0.2))), put(3, 1 => Rx(0.3)))

    u = chain(1, NoParams(Ry(1.0)), Rx(1.0))
    u2 = dispatch(u, [2.0]) # returns the following error
    @test u2 == chain(1, NoParams(Ry(1.0)), Rx(2.0))
    @test chsubblocks(u2[1], X) == NoParams(X)
end


text:
filename() = "day06.txt"

# read data
read_data() = split(String(read(filename())), "\n\n")

# Part 1
function part1()
    # how many unique poeple are in this set
    count_unique(s) = length(unique(collect(replace(s, "\n"=>""))))
    return count_unique.(read_data()) |> sum
end

# Part 2
# Count everyone who answered 'yes', so we just need to take the intersection
# of all people.
function count_everyone(line)
    people = split(line, "\n")
    set = Set(collect('a':'z'))
    for p in people
        new_set = collect(p)
        intersect!(set, new_set)
    end
    return length(set)
end

part2() = count_everyone.(read_data()) |> sum

# part 2 (refactored using splatting)
count_everyone_splat(s) = length(intersect(collect.(split(s, "\n"))...))
part2b() = count_everyone_splat.(read_data()) |> sum


text:
"""
    binary_tournament(population)
Apply binary tournament to obtain a solution from from population.
"""
function binary_tournament(population)
    a, b = rand(1:length(P), 2)

    return is_better(population[a], population[b]) ? population[a] : population[b]
end


text:
module UnitJuMP

using Reexport
@reexport using JuMP
@reexport using Unitful

import MutableArithmetics
_MA = MutableArithmetics

include("units.jl")
include("mutable_arithmetics.jl")
include("operators.jl")

export UnitAffExpr

end