text:
select(pars::AbstractArray{T,N}, elements_per_dim...; newfun = randoutzeroin) where {T, N} = NaiveNASlib.parselect(pars, elements_per_dim...; newfun)
select(::Missing, args...;kwargs...) = missing
select(::Flux.Zeros, args...;kwargs...) = Flux.Zeros()

struct WeightParam end
struct BiasParam end
struct RecurrentWeightParam end
struct RecurrentState end

"""
    neuroninsert(lt, partype)

Return a function which creates new parameters for layers of type `lt` to use for [`select`](@Ref).
"""
neuroninsert(t, partype) = randoutzeroin
neuroninsert(t, parname::Symbol) = neuroninsert(t, Val(parname))
neuroninsert(t, name::Val) = randoutzeroin

neuroninsert(lt::FluxParNorm, t::Val) = norminsert(lt, t)
norminsert(::FluxParNorm, ::Union{Val{:β},Val{:μ}}) = (args...) -> 0
norminsert(::FluxParNorm, ::Union{Val{:γ},Val{:σ²}}) = (args...) -> 1

# Coupling between input and output weights make it difficult to do anything else?
neuroninsert(::FluxDepthwiseConv, partype) = (args...) -> 0

randoutzeroin(T, d, s...) = _randoutzeroin(T,d,s)
_randoutzeroin(T, d, s) = 0
_randoutzeroin(T, d, s::NTuple{2, Int}) = d == indim(FluxDense()) ? 0 : randn(T, s) ./ prod(s)
_randoutzeroin(T, d, s::NTuple{N, Int}) where N = d == indim(FluxConv{N-2}()) ? 0 : randn(T, s) ./ prod(s)

"""
    KernelSizeAligned(Δsize; pad)
    KernelSizeAligned(Δs::Integer...;pad)

Strategy for changing kernel size of convolutional layers where filters remain phase aligned. In other words, the same
element indices are removed/added for all filters and only 'outer' elements are dropped or added.

Call with vertex as input to change weights.

### Examples

```jldoctest
julia> using NaiveNASflux, Flux

julia> cv = fluxvertex(Conv((3,3), 1=>1;pad=SamePad()), conv2dinputvertex("in", 1));

julia> cv(ones(Float32, 4,4,1,1)) |> size
(4, 4, 1, 1)

julia> layer(cv).weight |> size
(3, 3, 1, 1)

julia> cv |> KernelSizeAligned(-1, 1; pad=SamePad());

julia> cv(ones(Float32, 4,4,1,1)) |> size
(4, 4, 1, 1)

julia> layer(cv).weight |> size
(2, 4, 1, 1)
```
"""
struct KernelSizeAligned{T, P}
    Δsize::T
    pad::P
end
KernelSizeAligned(Δs::Integer...;pad = ntuple(i -> 0, length(Δs))) = KernelSizeAligned(Δs, pad)

(s::KernelSizeAligned)(l) = selectfilters(layertype(l), l, s)
(s::KernelSizeAligned)(v::AbstractVertex) = mutate_weights(v, s)

otherpars(s::KernelSizeAligned, l) = paddingfor(layertype(l), l, s)
paddingfor(lt, l, s) = ()
paddingfor(::FluxConvolutional{N}, l, s) where N = (;pad = Flux.calc_padding(typeof(l), s.pad, size(weights(l))[1:N] .+ s.Δsize, l.dilation, l.stride))

selectfilters(t, l, s) = ()
selectfilters(::FluxConvolutional, l, s) = selectfilters(s, weights(l))

function selectfilters(s::KernelSizeAligned, pars)
    csize = size(pars)
    mids = csize ./ 2 .+ 0.5
    byfun(dim) = i -> abs(mids[dim] - i) * 1 / (1 + mean(abs.(selectdim(pars, dim, i))))
    ps = Pair{Int, Vector{Int}}[]
    for (dim, Δsize) in enumerate(s.Δsize)
        if Δsize < 0
            els = sort(partialsort(1:csize[dim], 1:csize[dim]+Δsize, by=byfun(dim)))
            push!(ps, dim => els)
        elseif Δsize > 0
            els = -ones(Int, csize[dim] + Δsize)
            offs = Δsize ÷ 2
            els[offs+1:csize[dim]+offs] = 1:csize[dim]
            push!(ps, dim => els)
        end
    end
    return ps
end


text:
# We reverse the order of comparisons here so that the result
# of x < y is equal to the result of Monomial(x) < Monomial(y)
Base.@pure Base.isless(v1::AbstractVariable, v2::AbstractVariable) = name(v1) > name(v2)
Base.isless(m1::AbstractTermLike, m2::AbstractTermLike) = isless(promote(m1, m2)...)

# Implement this to make coefficients be compared with terms.
function isless_coefficient(a::Real, b::Real)
    return a < b
end
function isless_coefficient(a::Number, b::Number)
    return abs(a) < abs(b)
end
# By default, coefficients are not comparable so `a` is not strictly
# less than `b`, they are considered sort of equal.
isless_coefficient(a, b) = false

function Base.isless(t1::AbstractTerm, t2::AbstractTerm)
    if monomial(t1) < monomial(t2)
        return true
    elseif monomial(t1) == monomial(t2)
        return isless_coefficient(coefficient(t1), coefficient(t2))
    else
        return false
    end
end

# promoting multiplication is not a good idea
# For example a polynomial of Float64 * a polynomial of JuMP affine expression
# is a polynomial of JuMP affine expression but if we promote it would be a
# polynomial of quadratic expression
for op in [:+, :-, :(==)]
    @eval Base.$op(p1::APL, p2::APL) = $op(promote(p1, p2)...)
end
# Promotion between `I` and `1` is `Any`.
# Promotion between `I` and `2I` is `UniformScaling`.
for op in [:+, :-]
    @eval Base.$op(p1::APL, p2::APL{<:LinearAlgebra.UniformScaling}) = $op(p1, mapcoefficientsnz(J -> J.λ, p2))
    @eval Base.$op(p1::APL{<:LinearAlgebra.UniformScaling}, p2::APL) = $op(mapcoefficientsnz(J -> J.λ, p1), p2)
    @eval Base.$op(p1::APL{<:LinearAlgebra.UniformScaling}, p2::APL{<:LinearAlgebra.UniformScaling}) = $op(mapcoefficientsnz(J -> J.λ, p1), p2)
end
Base.isapprox(p1::APL, p2::APL; kwargs...) = isapprox(promote(p1, p2)...; kwargs...)

# @eval $op(p::APL, α) = $op(promote(p, α)...) would be less efficient
for (op, fun) in [(:+, :plusconstant), (:-, :minusconstant), (:*, :multconstant), (:(==), :eqconstant)]
    @eval Base.$op(p::APL, α) = $fun(p, α)
    @eval Base.$op(α, p::APL) = $fun(α, p)
end
# Fix ambiguity between above methods and methods in MA
Base.:+(::MA.Zero, p::APL) = MA.copy_if_mutable(p)
Base.:+(p::APL, ::MA.Zero) = MA.copy_if_mutable(p)
Base.:-(::MA.Zero, p::APL) = MA.operate(-, p)
Base.:-(p::APL, ::MA.Zero) = MA.copy_if_mutable(p)

# Special case AbstractArrays of APLs
# We add these instead of relying on the broadcasting API since the above method definitions are very wide.
# In particular, there is support for Matrices as coefficents. In order to avoid issues like #104 we therefore
# explicitly define this (instead of implictly getting unexpected results).
for op in [:+, :-]
    @eval Base.$op(p::APL, A::AbstractArray{<:APL}) = map(f -> $op(p, f), A)
    @eval Base.$op(A::AbstractArray{<:APL}, p::APL) = map(f -> $op(f, p), A)
end
Base.:*(p::APL, A::AbstractArray) = map(f -> p * f, A)
Base.:*(A::AbstractArray, p::APL) = map(f -> f * p, A)
Base.:/(A::AbstractArray, p::APL) = map(f -> f / p, A)

constant_function(::typeof(+)) = plusconstant
constant_function(::typeof(-)) = minusconstant
MA.operate!(op::Union{typeof(+), typeof(-)}, p::APL, α) = MA.operate!(constant_function(op), p, α)
MA.operate_to!(output::AbstractPolynomial, op::typeof(*), α, p::APL) = MA.operate_to!(output, multconstant, α, p)
MA.operate_to!(output::AbstractPolynomial, op::typeof(*), p::APL, α) = MA.operate_to!(output, multconstant, p, α)

function polynomial_merge!(
    n1::Int, n2::Int, get1::Function, get2::Function,
    set::Function, push::Function, compare_monomials::Function,
    combine::Function, keep::Function, resize::Function)
    buffer = nothing
    i = j = k = 1
    # Invariant:
    # The terms p[0] -> p[k-1] are sorted and are smaller than the remaining terms.
    # The terms p[k] -> p[i-1] are garbage.
    # The terms p[i] -> p[end] are sorted and still need to be added.
    # The terms q[j] -> p[end] are sorted and still need to be added.
    # If `buffer` is not empty:
    #   The terms in `buffer` are sorted and still need to be added.
    #   Moreover, they are smaller than the terms p[i] -> p[end].
    while i <= n1 && j <= n2
        @assert buffer === nothing || isempty(buffer)
        comp = compare_monomials(i, j)
        if comp > 0
            if k == i
                t0 = get1(i)
                if buffer === nothing
                    buffer = DataStructures.Queue{typeof(t0)}()
                end
                DataStructures.enqueue!(buffer, t0)
                i += 1
            end
            set(k, get2(j))
            j += 1
            k += 1
        elseif iszero(comp)
            combine(i, j)
            if keep(i)
                if k != i
                    @assert k < i
                    set(k, get1(i))
                end
                k += 1
            end
            i += 1
            j += 1
        else
            if k != i
                set(k, get1(i))
            end
            i += 1
            k += 1
        end
        while buffer !== nothing && !isempty(buffer) && j <= n2
            @assert i == k
            t = first(buffer)
            comp = compare_monomials(t, j)
            if comp >= 0
                if comp > 0
                    t = get2(j)
                else
                    t = combine(t, j)
                end
                j += 1
            end
            if comp <= 0
                DataStructures.dequeue!(buffer)
            end
            # if `comp` is zero, we called `combine` so `t`
            # might not be kept. If `comp` is not zero, we
            # skip the `keep` call that might be costly.
            if iszero(comp) && !keep(t)
                continue
            end
            if k <= n1
                DataStructures.enqueue!(buffer, get1(i))
                set(k, t)
            else
                push(t)
                n1 += 1
            end
            i += 1
            k += 1
        end
    end
    if buffer !== nothing && !isempty(buffer)
        @assert j == n2 + 1
        @assert i == k
        n = length(buffer)
        resize(n1 + n)
        for k in n1:-1:i
            set(k + n, get1(k))
        end
        for k in i:(i + n - 1)
            set(k, DataStructures.dequeue!(buffer))
        end
        n1 += n
    else
        len = (k - 1) + (n2 - (j - 1)) + (n1 - (i - 1))
        if n1 < len
            resize(len)
        end
        while j <= n2
            set(k, get2(j))
            j += 1
            k += 1
        end
        @assert j == n2 + 1
        while i <= n1
            set(k, get1(i))
            i += 1
            k += 1
        end
        if len < n1
            resize(len)
        end
        @assert i == n1 + 1
        @assert k == len + 1
    end
    return
end

#MA.operate!(op::Union{typeof(+), typeof(-)}, p::AbstractPolynomial, q::AbstractPolynomial) = MA.operate_to!(p, op, p, q)
MA.operate!(op::Union{typeof(+), typeof(-)}, p::AbstractPolynomial, q::AbstractPolynomialLike) = MA.operate!(op, p, polynomial(q))

function mul_to_terms!(ts::Vector{<:AbstractTerm}, p1::APL, p2::APL)
    for t1 in terms(p1)
        for t2 in terms(p2)
            push!(ts, t1 * t2)
        end
    end
    return ts
end
function Base.:*(p::AbstractPolynomial, q::AbstractPolynomial)
    polynomial!(mul_to_terms!(MA.promote_operation(*, termtype(p), termtype(q))[], p, q))
end

Base.isapprox(p::APL, α; kwargs...) = isapprox(promote(p, α)...; kwargs...)
Base.isapprox(α, p::APL; kwargs...) = isapprox(promote(p, α)...; kwargs...)

# `MA.operate(-, p)` redirects to `-p` as it assumes that `-p` can be modified
# through the MA API without modifying `p`. We should either copy the monomial
# here or implement a `MA.operate(-, p)` that copies it. We choose the first
# option.
Base.:-(m::AbstractMonomialLike) = _term(-1, MA.copy_if_mutable(m))
Base.:-(t::AbstractTermLike) = _term(MA.operate(-, coefficient(t)), monomial(t))
Base.:-(p::APL) = polynomial!((-).(terms(p)))
Base.:+(p::Union{APL, RationalPoly}) = p
Base.:*(p::Union{APL, RationalPoly}) = p

# Avoid adding a zero constant that might artificially increase the Newton polytope
# Need to add polynomial conversion for type stability
plusconstant(p::APL{S}, α::T)  where {S, T} = iszero(α) ? polynomial( p, MA.promote_operation(+, S, T)) : p + constantterm(α, p)
plusconstant(α::S, p::APL{T})  where {S, T} = iszero(α) ? polynomial( p, MA.promote_operation(+, S, T)) : constantterm(α, p) + p
function MA.operate!(::typeof(plusconstant), p::APL, α)
    if !iszero(α)
        MA.operate!(+, p, constantterm(α, p))
    end
    return p
end
minusconstant(p::APL{S}, α::T) where {S, T} = iszero(α) ? polynomial( p, MA.promote_operation(-, S, T)) : p - constantterm(α, p)
minusconstant(α::S, p::APL{T}) where {S, T} = iszero(α) ? polynomial(-p, MA.promote_operation(-, S, T)) : constantterm(α, p) - p
function MA.operate!(::typeof(minusconstant), p::APL, α)
    if !iszero(α)
        MA.operate!(-, p, constantterm(α, p))
    end
    return p
end

# Coefficients and variables commute
multconstant(α, v::AbstractVariable) = multconstant(α, monomial(v)) # TODO linear term
multconstant(m::AbstractMonomialLike, α) = multconstant(α, m)

_multconstant(α, f, t::AbstractTermLike) = mapcoefficientsnz(f, t)
function _multconstant(α::T, f, p::AbstractPolynomial{S}) where {S, T}
    if iszero(α)
        zero(polynomialtype(p, MA.promote_operation(*, T, S)))
    else
        mapcoefficientsnz(f, p)
    end
end
_multconstant(α, f, p::AbstractPolynomialLike) = _multconstant(α, f, polynomial(p))

multconstant(α, p::AbstractPolynomialLike) = _multconstant(α, β -> α*β, p)
multconstant(p::AbstractPolynomialLike, α) = _multconstant(α, β -> β*α, p)

function mapcoefficientsnz_to! end

function _multconstant_to!(output, α, f, p)
    if iszero(α)
        MA.operate!(zero, output)
    else
        mapcoefficientsnz_to!(output, f, p)
    end
end
function MA.operate_to!(output, ::typeof(multconstant), p::APL, α)
    _multconstant_to!(output, α, β -> β*α, p)
end
function MA.operate_to!(output, ::typeof(multconstant), α, p::APL)
    _multconstant_to!(output, α, β -> α*β, p)
end

MA.operate_to!(output::AbstractMonomial, ::typeof(*), m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents_to!(output, +, m1, m2)
MA.operate!(::typeof(*), m1::AbstractMonomial, m2::AbstractMonomialLike) = mapexponents!(+, m1, m2)
Base.:*(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(+, m1, m2)
#Base.:*(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = *(monomial(m1), monomial(m2))

Base.:*(m::AbstractMonomialLike, t::AbstractTermLike) = term(coefficient(t), m * monomial(t))
Base.:*(t::AbstractTermLike, m::AbstractMonomialLike) = term(coefficient(t), monomial(t) * m)
Base.:*(t1::AbstractTermLike, t2::AbstractTermLike) = term(coefficient(t1) * coefficient(t2), monomial(t1) * monomial(t2))

Base.:*(t::AbstractTermLike, p::APL) = polynomial!(map(te -> t * te, terms(p)))
Base.:*(p::APL, t::AbstractTermLike) = polynomial!(map(te -> te * t, terms(p)))
Base.:*(p::APL, q::APL) = polynomial(p) * polynomial(q)

# guaranteed that monomial(t1) > monomial(t2)
function _polynomial_2terms(t1::TT, t2::TT, ::Type{T}) where {TT<:AbstractTerm, T}
    if iszero(t1)
        polynomial(t2, T)
    elseif iszero(t2)
        polynomial(t1, T)
    else
        # not `polynomial!` because we `t1` and `t2` cannot be modified
        polynomial(termtype(TT, T)[t1, t2], SortedUniqState())
    end
end

_term(α, mono) = term(α, MA.copy_if_mutable(mono))

for op in [:+, :-]
    @eval begin
        Base.$op(t1::AbstractTermLike, t2::AbstractTermLike) = $op(term(t1), term(t2))
        Base.$op(t1::AbstractTerm, t2::AbstractTerm) = $op(_promote_terms(t1, t2)...)
        function Base.$op(t1::TT, t2::TT) where {T, TT <: AbstractTerm{T}}
            S = MA.promote_operation($op, T, T)
            # t1 > t2 would compare the coefficient in case the monomials are equal
            # and it will throw a MethodError in case the coefficients are not comparable
            if monomial(t1) == monomial(t2)
                polynomial(_term($op(coefficient(t1), coefficient(t2)), monomial(t1)), S)
            elseif monomial(t1) > monomial(t2)
                ts = _polynomial_2terms(t1, $op(t2), S)
            else
                ts = _polynomial_2terms($op(t2), t1, S)
            end
        end
    end
end
_promote_terms(t1, t2) = promote(t1, t2)
# Promotion between `I` and `1` is `Any`.
_promote_terms(t1::AbstractTerm, t2::AbstractTerm{<:LinearAlgebra.UniformScaling}) = _promote_terms(t1, coefficient(t2).λ * monomial(t2))
_promote_terms(t1::AbstractTerm{<:LinearAlgebra.UniformScaling}, t2::AbstractTerm) = _promote_terms(coefficient(t1).λ * monomial(t1), t2)
# Promotion between `I` and `2I` is `UniformScaling`, not `UniformScaling{Int}`.
function _promote_terms(t1::AbstractTerm{LinearAlgebra.UniformScaling{S}}, t2::AbstractTerm{LinearAlgebra.UniformScaling{T}}) where {S<:Number, T<:Number}
    U = LinearAlgebra.UniformScaling{promote_type(S, T)}
    _promote_terms(MA.scaling_convert(U, coefficient(t1)) * monomial(t1), MA.scaling_convert(U, coefficient(t2)) * monomial(t2))
end
function _promote_terms(t1::AbstractTerm{LinearAlgebra.UniformScaling{T}}, t2::AbstractTerm{LinearAlgebra.UniformScaling{T}}) where T<:Number
    return promote(t1, t2)
end

LinearAlgebra.adjoint(v::AbstractVariable) = v
LinearAlgebra.adjoint(m::AbstractMonomial) = m
LinearAlgebra.adjoint(t::AbstractTerm) = _term(LinearAlgebra.adjoint(coefficient(t)), monomial(t))
LinearAlgebra.adjoint(p::AbstractPolynomialLike) = polynomial(map(LinearAlgebra.adjoint, terms(p)))
LinearAlgebra.adjoint(r::RationalPoly) = adjoint(numerator(r)) / adjoint(denominator(r))

LinearAlgebra.transpose(v::AbstractVariable) = v
LinearAlgebra.transpose(m::AbstractMonomial) = m
LinearAlgebra.transpose(t::AbstractTerm) = _term(LinearAlgebra.transpose(coefficient(t)), monomial(t))
LinearAlgebra.transpose(p::AbstractPolynomialLike) = polynomial(map(LinearAlgebra.transpose, terms(p)))
LinearAlgebra.transpose(r::RationalPoly) = transpose(numerator(r)) / transpose(denominator(r))

LinearAlgebra.dot(p1::AbstractPolynomialLike, p2::AbstractPolynomialLike) = p1' * p2
LinearAlgebra.dot(x, p::AbstractPolynomialLike) = x' * p
LinearAlgebra.dot(p::AbstractPolynomialLike, x) = p' * x

LinearAlgebra.symmetric_type(PT::Type{<:APL}) = PT
LinearAlgebra.symmetric(p::APL, ::Symbol) = p

# Amazingly, this works! Thanks, StaticArrays.jl!
"""
Convert a tuple of variables into a static vector to allow array-like usage.
The element type of the vector will be Monomial{vars, length(vars)}.
"""
Base.vec(vars::Tuple{Vararg{AbstractVariable}}) = [vars...]
# vec(vars::Tuple{Vararg{<:TypedVariable}}) = SVector(vars)

# https://github.com/JuliaLang/julia/pull/23332
Base.:^(x::AbstractPolynomialLike, p::Integer) = Base.power_by_squaring(x, p)

function MA.operate_to!(output::AbstractPolynomial, op::MA.AddSubMul, x, args::Vararg{Any, N}) where N
    return MA.operate_to!(output, MA.add_sub_op(op), x, *(args...))
end
function MA.operate!(op::MA.AddSubMul, x::AbstractPolynomial, y, z, args::Vararg{Any, N}) where N
    return MA.operate!(MA.add_sub_op(op), x, *(y, z, args...))
end
MA.buffer_for(::MA.AddSubMul, ::Type{<:AbstractPolynomial}, args::Vararg{Type, N}) where {N} = zero(MA.promote_operation(*, args...))
function MA.buffered_operate_to!(buffer::AbstractPolynomial, output::AbstractPolynomial, op::MA.AddSubMul, x::AbstractPolynomial, y, z, args::Vararg{Any, N}) where N
    product = MA.operate_to!!(buffer, *, y, z, args...)
    return MA.operate_to!(output, MA.add_sub_op(op), x, product)
end
function MA.buffered_operate!(buffer::AbstractPolynomial, op::MA.AddSubMul, x::AbstractPolynomial, y, z, args::Vararg{Any, N}) where N
    product = MA.operate_to!!(buffer, *, y, z, args...)
    return MA.operate!(MA.add_sub_op(op), x, product)
end

ext: jl


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Implementation of the Camera abstract type and the derivating concrete types

###############
# Camera types

"""
    Camera

An abstract type representing an observer.

See also: [`OrthogonalCamera`](@ref), [`PerspectiveCamera`](@ref)
"""
abstract type Camera end

"""
    OrthogonalCamera <: Camera

A [`Camera`](@ref) implementing an orthogonal 3D → 2D projection.

# Fields

- `aspect_ratio::Float32`: defines how larger than the height is the image (16/9, 4/3, ...).
- `transformation::Transformation`: define the [`Transformation`](@ref) applied to the rays generated by the camera.

See also: [`fire_ray(::OrthogonalCamera, ::Float32, ::Float32)`](@ref)
"""
Base.@kwdef struct OrthogonalCamera <: Camera
    aspect_ratio::Float32 = 1f0
    transformation::Transformation = Transformation()
end

@doc """
    OrthogonalCamera(aspect_ratio::Float32, transformation::Transformation)

Constructor for an [`OrthogonalCamera`](@ref) instance.
""" OrthogonalCamera(::Float32, ::Transformation)

@doc """
    OrthogonalCamera(; aspect_ratio::Float32 = 1f0,
                       transformation::Transformation = Transformation())

Keyword-based constructor for an [`OrthogonalCamera`](@ref) instance.

If no parameter is specified, it return a camera with square aspect ratio and an identity transformation.
""" OrthogonalCamera(; ::Float32, ::Transformation)

"""
    PerspectiveCamera <: Camera

A [`Camera`](@ref) implementing a perspective 3D → 2D projection.

# Fields

- `aspect_ratio::Float32`: defines how larger than the height is the image (16/9, 4/3, ...).
- `transformation::Transformation`: define the [`Transformation`](@ref) applied to the rays generated by the camera.
- `screen_distance::Float32`: tells how much far from the eye of the observer is the screen and it influences the FOV (field-of-view).

See also: [`fire_ray(::PerspectiveCamera, ::Float32, ::Float32)`](@ref), [`aperture_deg`](@ref)
"""
Base.@kwdef struct PerspectiveCamera <: Camera
    aspect_ratio::Float32 = 1f0
    transformation::Transformation = Transformation()
    screen_distance::Float32 = 1f0
end

@doc """
    PerspectiveCamera(aspect_ratio::Float32, transformation::Transformation, screen_distance::Float32)

Constructor for an [`PerspectiveCamera`](@ref) instance.
""" PerspectiveCamera(::Float32, ::Transformation, ::Float32)

@doc """
    PerspectiveCamera(; aspect_ratio::Float32 = 1f0,
                        transformation::Transformation = Transformation(),
                        screen_distance::Float32 = 1f0)

Keyword-based constructor for a [`PerspectiveCamera`](@ref) instance.

If no parameter is specified, it return a camera with square aspect ratio, an identity transformation, and a
screen distance of 1, giving a FOV of 90°.
""" PerspectiveCamera(; ::Float32, ::Transformation, ::Float32)

##########
# Methods

"""
    fire_ray(camera::OrthogonalCamera, u::Float32, v::Float32)

Fire a [`Ray`](@ref) through an [`OrthogonalCamera`](@ref) at a position ``(u, v)`` on the screen,
using an orthogonal projection.

Parameters `u` and `v` are bound between 0 and 1:

    (0, 1)                            (1, 1)
        +------------------------------+
        |                              |
        |                              |
        |                              |
        +------------------------------+
    (0, 0)                            (1, 0)
"""
function fire_ray(camera::OrthogonalCamera, u::Float32, v::Float32)
    camera.transformation * Ray(Point(-1f0, (1f0 - 2u) * camera.aspect_ratio, 2v - 1f0),
                                VEC_X,
                                tmin = 1f0)
end

"""
    fire_ray(camera::PerspectiveCamera, u::Float32, v::Float32)

Fire a [`Ray`](@ref) through a [`PerspectiveCamera`](@ref) at a position ``(u, v)`` on the screen,
using a perspective projection.

Parameters `u` and `v` are bound between 0 and 1:

    (0, 1)                            (1, 1)
        +------------------------------+
        |                              |
        |                              |
        |                              |
        +------------------------------+
    (0, 0)                            (1, 0)
"""
function fire_ray(camera::PerspectiveCamera, u::Float32, v::Float32)
    camera.transformation * Ray(Point(-camera.screen_distance, 0f0, 0f0),
                                Vec(camera.screen_distance, (1f0 - 2u) * camera.aspect_ratio, 2v - 1f0),
                                tmin = 1f0)
end

"""
    aperture_deg(camera::PerspectiveCamera)

Compute the FOV of the camera in degrees for a [`PerspectiveCamera`](@ref).

# Examples

FOV for a camera with screen distance of 1 and aspect ratio of 1:

```jldoctest
julia> aperture_deg(PerspectiveCamera())
90.0f0
```

FOV for a camera with screen distance of 1 and aspect ratio of 16/9:

```jldoctest
julia> aperture_deg(PerspectiveCamera(aspect_ratio = 16//9))
58.715508f0
```
"""
aperture_deg(camera::PerspectiveCamera) = 2 * rad2deg(atan(camera.screen_distance, camera.aspect_ratio))


text:
using Test
using FakePkg
using FromFile
using Comonicon.Builder

@testset "FakePkg" begin
    @test FakePkg.CASTED_COMMANDS["main"].root.name == "pkg"

    @test FakePkg.command_main(["-h"]) == 0
    @test FakePkg.command_main(["--help"]) == 0
    @test FakePkg.command_main(["-V"]) == 0
    @test FakePkg.command_main(["--version"]) == 0

    @test FakePkg.command_main(["add", "-h"]) == 0
    @test FakePkg.command_main(["rm", "-h"]) == 0
    @test FakePkg.command_main(["noarguments", "-h"]) == 0

    @test FakePkg.command_main(["rm", "CCC", "-h"]) == 0
    @test FakePkg.command_main(["add", "CCC", "-h"]) == 0

    @test FakePkg.command_main(["add", "ABC"]) == 0
    @test FakePkg.command_main(["rm", "ABC"]) == 0
    @test FakePkg.command_main(["activate", "fake", "-s"]) == 0
    @test FakePkg.command_main(["activate", "fake", "--shared"]) == 0

    @test FakePkg.command_main(["registry", "add", "abc"]) == 0
    @test FakePkg.command_main(["registry", "rm", "abc"]) == 0
end

using Comonicon.Options: read_options, get_path, @asset_str

@testset "test assets path" begin
    options = read_options(FakePkg)
    for asset in options.application.assets
        @test ispath(get_path(FakePkg, asset))
    end
end

@from "../../../test/utils.jl" import with_args

@testset "build package" begin
    @test with_args(["-h"]) do
        FakePkg.comonicon_install()
    end == 0

    @test with_args() do
        FakePkg.comonicon_install()
    end == 0

    @test isfile(joinpath(".julia", "bin", "pkg"))
    if haskey(ENV, "SHELL") && basename(ENV["SHELL"]) == "zsh"
        @test isfile(joinpath(".julia", "completions", "_pkg"))
    end

    @test with_args(["sysimg"]) do
        FakePkg.comonicon_install()
    end == 0

    @test with_args(["sysimg", "tarball"]) do
        FakePkg.comonicon_install()
    end == 0

    @test with_args(["app"]) do
        FakePkg.comonicon_install()
    end == 0

    @test with_args(["app", "tarball"]) do
        FakePkg.comonicon_install()
    end == 0

    @test isdir("build")
    @test isfile(Builder.tarball_name(FakePkg, "pkg", "application"))
end


text:
using ErlangTermFormat
using Test

isdir("data") || error("Run generator.exs before testing")

const etfs = read.(joinpath.("data", readdir("data")), String)

@testset "ErlangTermFormat.jl" begin
    for e in etfs
        @test encode(decode(e)) == e
        # We should also be able to handle extra trailing data.
        @test encode(decode(e * "foo")) == e
    end
end


text:
module VectorSpaces

export VectorSpace, ZeroSpace, dim, span, span_transposed, kernel, image, preimage

using LinearAlgebra
using LinAlgBigFloat

export rank

import Base: show, size, copy, adjoint, union, intersect, in, issubset, (==), (*), (\)
import LinearAlgebra: rank

mutable struct VectorSpace{T<:Number,S<:AbstractMatrix}
  n::Int
  rank::Int
  first::Int
  qrf::Union{QRPivoted{T, S},Nothing}
  VectorSpace{T,S}(n::Int, rank::Int, first::Int, qrf::Union{QRPivoted{T,S},Nothing}) where {T,S<:AbstractMatrix} = new(n, rank, first, qrf)
end

function show(io::IO, vs::VectorSpace)
  println(io, typeof(vs), " rank ", vs.rank, " dim ", vs.n)
end

"""

  `VectorSpace(A::Union{AbstractMatrix,AbstractVector})`

  Create VectorSpace as subspace spanned by the column vectors of A.
  The base space has dimension of the number of rows of A.
  The dimension of the subspace is the (numerical) rank of A.
"""
function VectorSpace(A::AbstractMatrix{T}) where T <: Number
  n, m = size(A)
  r = 0
  if m > 0
    tol = tolerance(A)
    qrf = qr(A, ColumnNorm())
    r = rank(qrf, tol)
  end
  if r <= 0 || r >= n
    qrf = nothing
  end
  VectorSpace{T,typeof(A)}(n, r, 1, qrf)
end

VectorSpace(B::AbstractVector{T}) where T <: Number = VectorSpace(reshape(B, size(B,1), 1))
"""

  `ZeroSpace([T<:Number,] n) -> VectorSpace`

  Create zero space in base space of dimension `n`. Default type is Float64.
"""
ZeroSpace(::Type{T}, n::Int) where T <: Number = VectorSpace{T,Matrix{T}}(n, 0, 1, nothing)
ZeroSpace(n::Int) = ZeroSpace(Float64, n)

"""

  `VectorSpace([T<:Number], n)`

  Create base space of dimension`n`. Default type is Float64.
"""
VectorSpace(::Type{T}, n::Int) where T <: Number = VectorSpace{T,Matrix{T}}(n, n, 1, nothing)
VectorSpace(n::Int) = VectorSpace(Float64, n)

size(vs::VectorSpace) = (vs.n, vs.rank)
rank(vs::VectorSpace) = vs.rank
dim(vs::VectorSpace) = vs.n

copy(vs::VectorSpace{T,S}) where {T, S} =
  VectorSpace{T,S}(vs.n, vs.rank, vs.first, vs.qrf)

"""

  `adjoint(vs::VectorSpace) -> VectorSpace`

  Create VectorSpace representing the orthogonal complement of `vs`.
"""
function adjoint(vs::VectorSpace{T,S}) where {T,S}
  n, r = size(vs)
  if vs.first == 1
    f = r + 1
  else
    f = 1
  end
  r = n - r
  VectorSpace{T,S}(n, r, f, vs.qrf)
end

"""

  `span(vs::VectorSpace) -> Matrix`

  Return unitary matrix containing basis vectors of `vs`.
"""
function span(vs::VectorSpace{T}) where T
  if 0 < vs.rank < vs.n
    r = _range1(vs)
    vs.qrf.Q[:,r]
  else
    Matrix{T}(I, vs.n, vs.rank)
  end
end

"""

  `span_transposed(vs::VectorSpace) -> Matrix`

  Return unitary matrix containing basis vectors of complement of `vs`.
  Generally `span(vs') == span_transposed(vs) is always true`.
"""
function span_transposed(vs::VectorSpace{T}) where T
  if 0 < vs.rank < vs.n
    r = _range2(vs)
    vs.qrf.Q[:,r]
  else
    Matrix{T}(I, vs.n, vs.n - vs.rank)
  end
end

@inline function _range1(vs::VectorSpace)
  f = vs.first
  l = vs.rank + f - 1
  f:l
end

@inline function _range2(vs::VectorSpace)
  n, r = size(vs)
  f = vs.first == 1 ? r + 1 : 1
  l = n - r + f - 1
  f:l
end

"""
  `*(A::AbstractMatrix, vs::VectorSpace) -> VectorSpace`

  Calculate the image of the vector space `vs` under matrix `A`.
  Discard small column vectors in `A * span(vs)`.
  Use `A * vs` rather than `VectorSpace(A * span(A))` for accuracy.
"""
function *(A::AbstractMatrix{T}, vs::VectorSpace{T}) where T <: Number
  n, m = size(A)
  nv, r = size(vs)
  m == nv || error("dimension mismatch")

  if r == 0
    ZeroSpace(T, n)
  elseif r == nv
    VectorSpace(A)
  else
    spana = A * span(vs)
    tol = tolerance(A)
    for k = 1:r
      if norm(spana[:,k], Inf) < tol
        spana[:,k] = 0
      end
    end
    VectorSpace(spana)
  end
end

"""
  `\\(A::AbstractMatrix, vs::VectorSpace) -> VectorSpace`

  Calculate the preimage of the vector space `vs` under matrix `A`.
  Discard small row vectors in `span(vs')' * A`.
  Use `A \\ vs` rather than `VectorSpace(span(vs')'A)` for accuracy.
"""
function \(A::AbstractMatrix{T}, vs::VectorSpace{T}) where T <: Number
  n, m = size(A)
  nv, r = size(vs)
  n == nv || error("dimension mismatch")

  if r == nv
    VectorSpace(T, m)
  elseif r == 0
    kernel(A)
  else
    spana = span_transposed(vs)'A
    tol = tolerance(A)
    for k = 1:nv-r
      if norm(spana[k,:], Inf) < tol
        spana[k,:] = 0
      end
    end
    kernel(spana)
  end
end

"""

`image(A:AbstractMatrix[, vs::VectorSpace] [, k = 1]) -> VectorSpace`

  Create image of `vs` under `A^k`. If k ≠ 1 A must be square matrix.
  If vs is not given, the usual image (range) of `A` is returned.
"""
image(A::AbstractArray{T}) where T <: Number = VectorSpace(A)
function image(A::AbstractArray{T}, k::Int) where T <: Number
  k == 1 ?  VectorSpace(A) : image(A, VectorSpace(T, size(A,2)), k)
end

function image(A::AbstractMatrix, vs::VectorSpace, k::Int = 1)
  k ≥  0 || error("exponent ≥ 0 required")
  n, m = size(A)
  na, ra = size(vs)
  n == m || k == 1 || error("exponent k ≠ 1 only for square matrix")
  na == m || error("matrix dimension 2 does not match source space dimension")

  va = vs
  if k > 0
    va = ra == m ? image(A) : A * va
    r = rank(va)
    while k > 1 && 0 < r < n
      vaa = A * va
      n = r
      r = rank(vaa)
      k -= 1
      if r < n
        va = vaa
      end
    end
  end
  va
end

"""

  `kernel(A::AbstractArray) -> VectorSpace`

  Create VectorSpace representing the kernel of the mapping given by matrix `A`.
"""
kernel(A::AbstractArray) =  VectorSpace(kernel_matrix(A))

"""

  `preimage(A::AbstractArray[, vs::VectorSpace] [, k::Int=1])`

  Create vector space, which is the inverse image of `vs` under the mapping `A^k`.
  `vs` defaults to the zero space; in that case the kernel of `A^k` is returned.
"""
preimage(A::AbstractArray{T}) where T<:Number = kernel(A)
function preimage(A::AbstractArray{T}, k::Int) where T <:Number
  k == 1 ? kernel(A) : preimage(A, ZeroSpace(T,size(A,1)), k)
end

function preimage1(A::AbstractArray{T}, vs::VectorSpace, k::Int = 1) where T <: Number
  k ≥  0 || error("exponent ≥ 0 required")
  n, m = size(A)
  na, ra = size(vs)
  n == m || k == 1 || error("exponent k ≠ 1 only for square matrix")
  na == n || error("matrix dimension 1 does not match target space dimension")

  va = vs
  qra = qrfact(A, Val{true})
  tol = tolerance(A)
  r = rank(qra)
  va = zeros(T, m, m-r+ra)
  R = view(qra[:R], 1:r, 1:r)
  va[1:r,1:m-r] = -R \ qra[:R][1:r,r+1:m]
  va[1:r,m-r+1:m-r+ra] = R \ (qra[:Q] * span(va))[1:r,:]
  va[r+1:m,1:m-r] = 0

end

function preimage(A::AbstractArray, vs::VectorSpace, k::Int = 1)
  k ≥  0 || error("exponent ≥ 0 required")
  n, m = size(A)
  na, ra = size(vs)
  n == m || k == 1 || error("exponent k ≠ 1 only for square matrix")
  na == n || error("matrix dimension 1 does not match target space dimension")

  va = vs
  if k > 0
    va = ra == 0 ? kernel(A) : A \ va
    rp = 0
    r = va.rank
    while k > 1 && r > rp
      vaa = A \ va
      rp = r
      r = rank(vaa)
      k -= 1
      if r > rp
        va = vaa
      end
    end
  end
  va
end

"""

  `union(va::VectorSpace, vb::VectorSpace) -> VectorSpace`

  Create VectorSpace representing `va ∪ vb`.
"""
function union(va::VectorSpace, vb::VectorSpace)
  na, ra = size(va)
  nb, rb = size(vb)
  na == nb || error("dimension mismatch")

  if ra >= na || rb >= nb
    VectorSpace{eltype(va)}(va.na)
  elseif ra == 0
    vb
  elseif rb == 0
    va
  else
    A = [span(va) span(vb)]
    VectorSpace(A)
  end
end

"""

  `intersect(va::VectorSpace, vb::VectorSpace) -> VectorSpace`

  Create VectorSpace representing `va ∩ vb`.
"""
function intersect(va::VectorSpace, vb::VectorSpace)
  na, ra = size(va)
  nb, rb = size(vb)
  na == nb || error("dimension mismatch")

  if ra == 0 || rb == 0
    ZeroSpace(typeof(va), va.na)
  elseif ra == na
    vb
  elseif rb == nb
    va
  else
    U = span(va)
    V = span(vb)
    K = kernel_matrix([U -V])
    VectorSpace(U * K[1:rank(va),:])
  end
end

"""

  `in(A::AbstractArray, vb::VectorSpace) -> Bool`

  Determine if `A[:,k] ∈ vb ∀ k`.
"""
function in(A::AbstractArray, vs::VectorSpace)
  issubset(VectorSpace(A), vs)
end

"""

  `issubset(va::VectorSpace, vb::VectorSpace) -> Bool`

  Determine if `va ⊆ vb`.
"""
function issubset(va::VectorSpace{T}, vb::VectorSpace{T}) where T <: Number
  na, ra = size(va)
  nb, rb = size(vb)
  na == nb || error("vector spaces not in same base space")

  rb >= nb || ra <= 0 && return true
  ra > rb && return false
  spanba = span_transposed(vb)'span(va)
  norm(vec(spanba), Inf) <= 4eps(T) * na
end

"""

  ` va == vb `

  Determine if VectorSpaces are (numerically) equal.
  Equivalent to `va ⊆ vb and vb ⊆ va`.
"""
function ==(va::VectorSpace, vb::VectorSpace)
  va.rank == vb.rank &&
  issubset(va, vb)
end

#########################################################################################
function rank(A::AbstractMatrix{T}) where T <: Number
  tol = tolerance(A)
  rank(qrfact(A, Val{true}), tol)
end

function rank(QR::LinearAlgebra.QRPivoted{T, S}, tol::AbstractFloat) where S <: AbstractMatrix{T} where T <: Number
  minimum(size(QR)) == 0 && ( return 0 )
  R = QR.R
  sv = sort(abs.(diag(R)), rev = true)
  n = length(sv)
  while n > 0 && sv[n] <= tol
    n -= 1
  end
  n
end

"""
  `tolerance(A::AbstractMatrix)`

  standard tolerance for a matrix.
"""
@inline tolerance(A::AbstractMatrix) = eps(norm(A,1)) * maximum(size(A))

"""

  `kernel_matrix(A::AbstractMatrix) -> Matrix`

Produce unitary matrix, whose vectors are basis of the kernel of Matrix A.
If the matrix is injective, the a matrix with zero columns is returned.
If the matrix has rank zero, the unit matrix is returned.
"""
function kernel_matrix(A::AbstractMatrix{T}) where T <: Number

  tol = tolerance(A)
  QR = qr(A, ColumnNorm())
  r = rank(QR, tol)
  n, m = size(A)
  K = zeros(T, m, m - r)
  for k = r+1:m
    K[k,k-r] = one(T)
  end
  if 0 < r < m
    r1r = 1:r
    rrm = r+1:m
    R = QR[:R]
    RV = view(R, r1r, r1r)
    K[r1r,:] = -RV \ R[r1r, rrm]
    p = sortperm(QR[:p]) # inverse of permutation
    K, R = qr(K[p,:])
  end
  K
end

function kernel_matrix(A::AbstractVector{T}) where T <: Number
  kernel_matrix(reshape(A, size(A,1), 1))
end

end # module

max_stars_repo_name: KlausC/LinAlgBigFloat

max_issues_repo_name: KlausC/LinAlgBigFloat

max_forks_repo_name: KlausC/LinAlgBigFloat


text:
using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--part2"
    help = "Do part 2"
    action = :store_true
    "--test"
    help = "Run test data"
    action = :store_true
end
parsed_args = parse_args(ARGS, s)

test_suffix = parsed_args["test"] ? "_test" : ""
chunks = readlines("input/day_10$test_suffix.txt")
const opens = Dict('(' => ')', '{' => '}', '<' => '>', '[' => ']')
const scores = Dict(')' => 3, ']' => 57, '}' => 1197, '>' => 25137)

function assess_chunk(c)
    stack = []
    for chr in c
        if haskey(opens, chr)
            push!(stack, opens[chr])
        elseif isempty(stack) || pop!(stack) != chr
            return (score = scores[chr], stack = [])
        end
    end
    (score = 0, stack = stack)
end

function part1()
    sum(map(c -> assess_chunk(c).score, chunks))
end

const completion_scores = Dict(')' => 1, ']' => 2, '}' => 3, '>' => 4)
function part2()
    stacks = map(ass -> ass.stack, filter(ass -> ass.score == 0, map(assess_chunk, chunks)))
    function score_stack(stack)
        score = 0
        while !isempty(stack)
            score = score * 5 + completion_scores[pop!(stack)]
        end
        score
    end
    scores = sort(map(score_stack, stacks))
    mid = div(length(scores) + 1, 2)
    scores[mid]
end

p = parsed_args["part2"] ? part2 : part1
answer = p()
println("The answer is: $answer")


text:
using Oceananigans.Diagnostics: AdvectiveCFL, DiffusiveCFL
using Oceananigans.Simulations: TimeStepWizard
using Oceananigans.Utils: prettytime
using Oceananigans: iteration, time
using Printf

export SimpleProgressMessenger, SingleLineProgressMessenger, TimedProgressMessenger

function SimpleProgressMessenger_function(simulation; LES=false, SI_units=true,
                                           initial_wall_time_seconds=1e-9*time_ns())
    model = simulation.model
    Δt = simulation.Δt

    iter, t = iteration(simulation), time(simulation)

    progress = 100 * (t / simulation.stop_time)

    current_wall_time = 1e-9 * time_ns() - initial_wall_time_seconds

    u_max = maximum(abs, model.velocities.u)
    v_max = maximum(abs, model.velocities.v)
    w_max = maximum(abs, model.velocities.w)

    if SI_units
        @info @sprintf("[%06.2f%%] iter: % 6d,     time: % 10s,     Δt: % 10s,     wall time: % 8s",
                        progress, iter, prettytime(t), prettytime(Δt), prettytime(current_wall_time),)
        if LES
            ν_max = maximum(abs, model.diffusivity_fields.νₑ)
            @info @sprintf("          └── max(|u⃗|): [%.2e, %.2e, %.2e] m/s,     adv CFL: %.2e,     diff CFL: %.2e,     νₘₐₓ: %.2e m²/s",
                           u_max, v_max, w_max, AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model), ν_max)
        else
            @info @sprintf("          └── max(|u⃗|): [%.2e, %.2e, %.2e] m/s,     adv CFL: %.2e,     diff CFL: %.2e",
                           u_max, v_max, w_max, AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model))
        end

    else
        @info @sprintf("[%06.2f%%] iter: % 6d,     time: %13.2f,     Δt: %13.2f,     wall time: % 8s",
                       progress, iter, t, Δt, prettytime(current_wall_time))
        if LES
            ν_max = maximum(abs, model.diffusivity_fields.νₑ)
            @info @sprintf("          └── max(|u⃗|): [%.2e, %.2e, %.2e],     adv CFL: %.2e,     diff CFL: %.2e,     νₘₐₓ: %.2e",
                           u_max, v_max, w_max, AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model), ν_max)
        else
            @info @sprintf("          └── max(|u⃗|): [%.2e, %.2e, %.2e],     adv CFL: %.2e,     diff CFL: %.2e",
                           u_max, v_max, w_max, AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model))
        end
    end

    @info ""

    return nothing
end

function SimpleProgressMessenger(; kwargs...)
    return simulation -> SimpleProgressMessenger_function(simulation; kwargs...)
end

function SingleLineProgressMessenger_func(simulation; LES=false, SI_units=true,
                                          initial_wall_time_seconds=1e-9*time_ns())
    model = simulation.model
    Δt = simulation.Δt

    iter, t = iteration(simulation), time(simulation)

    progress = 100 * (t / simulation.stop_time)

    current_wall_time = 1e-9 * time_ns() - initial_wall_time_seconds

    if SI_units
        if LES
            ν_max = maximum(abs, model.diffusivity_fields.νₑ)
            @info @sprintf("[%06.2f%%] iter: % 6d,     time: %10s,     Δt: %10s,     wall time: %8s,     adv CFL: %.2e,     diff CFL: %.2e,     νₘₐₓ: %.2e m²/s",
                           progress, iter, prettytime(t), prettytime(Δt), prettytime(current_wall_time),
                           AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model), ν_max)
        else
            @info @sprintf("[%06.2f%%] iter: % 6d,     time: %10s,     Δt: %10s,     wall time: %8s,     adv CFL: %.2e,     diff CFL: %.2e",
                           progress, iter, prettytime(t), prettytime(Δt), prettytime(current_wall_time),
                           AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model))
        end

    else
        if LES
            ν_max = maximum(abs, model.diffusivity_fields.νₑ)
            @info @sprintf("[%06.2f%%] iter: % 6d,     time: %13.2f,     Δt: %13.2f,     wall time: % 8s,     adv CFL: %.2e,     diff CFL: %.2e,     νₘₐₓ: %.2e",
                           progress, iter, t, Δt, prettytime(current_wall_time), AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model), ν_max)
        else
            @info @sprintf("[%06.2f%%] iter: % 6d,     time: %13.2f,     Δt: %13.2f,     wall time: % 8s,     adv CFL: %.2e,     diff CFL: %.2e",
                           progress, iter, t, Δt, prettytime(current_wall_time), AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model))
        end
    end

    return nothing
end

function SingleLineProgressMessenger(; kwargs...)
    return simulation -> SingleLineProgressMessenger_func(simulation; kwargs...)
end

mutable struct TimedProgressMessenger{T, I, L} <: Function
    wall_time₀ :: T  # Wall time at simulation start
    wall_time⁻ :: T  # Wall time at previous calback
    iteration⁻ :: I  # Iteration at previous calback
           LES :: L
end

function TimedProgressMessenger(; LES=false,
                                wall_time₀=1e-9*time_ns(),
                                wall_time⁻=1e-9*time_ns(),
                                iteration⁻=0,
    )

    return TimedProgressMessenger(wall_time₀, wall_time⁻, iteration⁻, LES)
end

function (pm::TimedProgressMessenger)(simulation)
    model = simulation.model
    Δt = simulation.Δt
    adv_cfl = AdvectiveCFL(simulation.Δt)(model)
    dif_cfl = DiffusiveCFL(simulation.Δt)(model)

    iter, t = iteration(simulation), time(simulation)

    progress = 100 * (t / simulation.stop_time)

    current_wall_time = 1e-9 * time_ns() - pm.wall_time₀
    time_since_last_callback = 1e-9 * time_ns() - pm.wall_time⁻
    iterations_since_last_callback = iter==0 ? Inf : iter - pm.iteration⁻
    wall_time_per_step = time_since_last_callback / iterations_since_last_callback
    pm.wall_time⁻ = 1e-9 * time_ns()
    pm.iteration⁻ = iter

    u_max = maximum(abs, model.velocities.u)
    v_max = maximum(abs, model.velocities.v)
    w_max = maximum(abs, model.velocities.w)

    @info @sprintf("[%06.2f%%] iteration: % 6d, time: % 10s, Δt: % 10s, wall time: % 8s (% 8s / time step)",
                    progress, iter, prettytime(t), prettytime(Δt), prettytime(current_wall_time), prettytime(wall_time_per_step))

    if pm.LES
        ν_max = maximum(abs, model.diffusivity_fields.νₑ)
        @info @sprintf("          └── max(|u⃗|): [%.2e, %.2e, %.2e] m/s, CFL: %.2e, νCFL: %.2e, ν_max: %.2e m²/s",
                        u_max, v_max, w_max, adv_cfl, dif_cfl, ν_max)
    else
        @info @sprintf("          └── max(|u⃗|): [%.2e, %.2e, %.2e] m/s, CFL: %.2e, νCFL: %.2e",
                        u_max, v_max, w_max, adv_cfl, dif_cfl)
    end

    @info ""

    return nothing
end


text:
"""
    fixedStateCV(
        states :: Vector{<:AbstractMatrix},
        kraus_ops :: Vector{<:AbstractMatrix}
    ) :: Tuple{Float64, Vector{Matrix{ComplexF64}}}

For a fixed ensemble of `states` and quantum channel described by `kraus_ops`,
the communication value (CV) is evaluated by maximizing over all POVM measurements.
This optimization is expressed in primal form as the semidefinite program:

```math
\\begin{matrix}
    \\text{fixedStateCV}(\\mathcal{N})&= \\max_{\\{\\Pi_x\\}_{x}} \\sum_x \\text{Tr}\\left[\\Pi_x\\mathcal{N}(\\rho_x)\\right] \\\\
    & \\\\
    & \\text{s.t.} \\quad \\sum_x\\Pi_x = \\mathbb{I} \\;\\; \\text{and} \\;\\; \\Pi_x \\geq 0
\\end{matrix}
```

where each state ``\\rho_x`` satisfies ``\\text{Tr}[\\rho_x] = 1`` and ``\\rho_x \\geq 0``.
The channel ``\\mathcal{N}`` is applied to each state as
``\\mathcal{N}(\\rho_x) = \\sum_j k \\rho_x k^\\dagger``.

# Returns

A `Tuple`, `(cv, opt_povm)` where `cv` is the evaluated communication value and
`opt_povm` is the optimal POVM measurement.
"""
function fixedStateCV(
    states :: Vector{<:AbstractMatrix},
    kraus_ops :: Vector{<:AbstractMatrix}
) :: Tuple{Float64, Vector{Matrix{ComplexF64}}}
    d = size(kraus_ops[1],1)
    n = length(states)

    # add povm variables and constraints
    povm_vars = map(i -> HermitianSemidefinite(d), 1:n)
    constraints = sum(map(Πy -> real(Πy), povm_vars)) == Matrix{Float64}(I, d, d)
    constraints += sum(map(Πy -> imag(Πy), povm_vars)) == zeros(Float64, d, d)

    # apply channel to states
    evolved_states = map(ρ ->  sum(k -> k * ρ * k' , kraus_ops), states)

    # maximize CV over POVM measurements
    objective = maximize(real(sum(tr.(povm_vars .* evolved_states))), constraints)
    qsolve!(objective)

    cv = objective.optval
    opt_povm = map(Πy -> Πy.value, povm_vars)

    return cv, opt_povm
end

"""
    fixedMeasurementCV(
        povm :: Vector{<:AbstractMatrix},
        kraus_ops :: Vector{<:AbstractMatrix}
    ) :: Tuple{Float64, Vector{Matrix{ComplexF64}}}

For a fixed `povm` measurement and quantum channel described by `kraus_ops`, the
communication value (CV) and optimal state encodings are computed.
The fixed measurement CV is evaluated as

```math
\\text{fixedMeasurementCV}(\\mathcal{N}) = \\sum_y ||\\mathcal{N}^{\\dagger}(\\Pi_y)||_{\\infty}
```

where ``||\\mathcal{N}^{\\dagger}(\\Pi_y)||_{\\infty}`` is the largest eigenvalue
of the POVM element ``\\Pi_y`` evolved by the adjoint channel,
``\\mathcal{N}^{\\dagger}(\\Pi_y) = \\sum_j k^{\\dagger}_j \\Pi_y k_j``.
The states which maximize the CV are simply the eigenvectors corresponding to the
largest eigenvalue of each respective POVM element.

# Returns

A `Tuple`, `(cv, opt_states)` where `cv` is the communication value and
`opt_states` is the set of optimal states.
"""
function fixedMeasurementCV(
    povm :: Vector{<:AbstractMatrix},
    kraus_ops :: Vector{<:AbstractMatrix}
) :: Tuple{Float64, Vector{Matrix{ComplexF64}}}
    d = size(kraus_ops[1],2)
    n = length(povm)

    # apply adjoint channel to POVM
    evolved_povm = map(Π ->  sum(k -> k' * Π * k , kraus_ops), povm)

    opt_states = Vector{Matrix{ComplexF64}}(undef, n)
    cv = 0
    for i in 1:n
        povm_el = evolved_povm[i]
        vals, vecs = eigen(povm_el)

        opt_states[i] = vecs[:,end] * vecs[:,end]'
        cv += real(vals[end])
    end

    return cv, opt_states
end

"""
    seesawCV(
        init_states :: Vector{<:AbstractMatrix},
        kraus_ops :: Vector{<:AbstractMatrix},
        num_steps :: Int64;
        verbose :: Bool = false
    )

Performs the see-saw optimization technique to maximize the communication
value (CV) of the channel described by `kraus_ops` over all states and measurements.
This iterative and biconvex optimization technique combines coordinate ascent
maximization with semidefinite programming.
The number of iterations is determined by `num_steps` where each iteration
consists of a two-step procedure:

1. The POVM measurement is optimized with respect to a fixed state ensemble
   using the [`fixedStateCV`](@ref) function.
2. The state ensemble is optimized with respect to a fixed povm state using
   the [`fixedMeasurementCV`](@ref) function.

This procedure is initialized with `init_states` and after many iterations, a
local maximum of the CV is found.
The `verbose` keyword argument can be used to print out the CV evaluated in each
step.

The see-saw method has shown success in similar encoding/decoding
optimization problems in quantum information, *e.g.*,
[https://arxiv.org/abs/quant-ph/0307138v2](https://arxiv.org/abs/quant-ph/0307138v2)
and [https://arxiv.org/abs/quant-ph/0606078v1](https://arxiv.org/abs/quant-ph/0606078v1).
We note that our implementation is quite distinct from previous works, however,
the core iterative approach remains the same.

# Returns

A `Tuple` containing the following data in order:

1. `max_cv_tuple :: Tuple`, `(max_cv, opt_states, opt_povm)` A 3-tuple containing the maximal
   communication value and the optimal states/POVM that achieve this value.
2. `cvs :: Vector{Float64}`, A list of each evaluated CV. Since states and measurements
   are optimized in each iteration, we have `length(cvs) == 2 * num_steps`.
3. `opt_ensembles :: Vector{Vector{Matrix{ComplexF64}}}`, A list of state ensembles
   optimized in each step, where `length(opt_ensembles) == num_steps`.
4. `opt_povms :: Vector{Vector{Matrix{ComplexF64}}}`, A list of POVM measurements
   optimized in each step, where `length(opt_povms) == num_steps`.

!!! warning "Optimum Not Guaranteed"
    This function is not guaranteed to find a global or local optima. However,
    `seesawCV` will always provide a lower bound on the communication value.
"""
function seesawCV(
    init_states :: Vector{<:AbstractMatrix},
    kraus_ops :: Vector{<:AbstractMatrix},
    num_steps :: Int64;
    verbose :: Bool = false
) :: Tuple
    opt_ensembles = Vector{Vector{Matrix{ComplexF64}}}(undef, num_steps)
    opt_povms = Vector{Vector{Matrix{ComplexF64}}}(undef, num_steps)
    cvs = zeros(Float64, 2*num_steps)
    cv_id = 1

    max_cv_tuple = (1, [], [])

    for i in 1:num_steps
        # maximizing CV over POVMs
        opt_states = (i == 1) ? init_states : opt_ensembles[i-1]

        fixed_state_cv, opt_povm = fixedStateCV(opt_states, kraus_ops)

        cvs[cv_id] = fixed_state_cv
        cv_id += 1
        opt_povms[i] = opt_povm

        if fixed_state_cv > max_cv_tuple[1]
            max_cv_tuple = (fixed_state_cv, opt_states, opt_povm)
        end

        # maximizing CV over states
        fixed_povm_cv, opt_states = fixedMeasurementCV(opt_povm, kraus_ops)

        cvs[cv_id] = fixed_povm_cv
        cv_id += 1
        opt_ensembles[i] = opt_states

        if fixed_povm_cv > max_cv_tuple[1]
            max_cv_tuple = (fixed_povm_cv, opt_states, opt_povm)
        end

        if verbose
            println("i = ", i)
            println("fixed_state_cv = ", fixed_state_cv)
            println("fixed_povm_cv = ", fixed_povm_cv)
            println("max_cv = ", max_cv_tuple[1])
        end
    end

    return max_cv_tuple, cvs, opt_ensembles, opt_povms
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
"""
    DistanceMethod

Abstract supertype for methods of calculating distances and discretised dispersal
probabilities between cells in a grid.
between cells in a grid.

Distance calculation methods include:
- [`CentroidToCentroid`](@ref)
- [`AreaToArea`](@ref)
- [`AreaToCentroid`](@ref)

Which are adapted from: "On the approximation of continuous dispersal kernels
in discrete-space models, Joseph D. Chipperfield et al 2011"

The `CentroidToArea` method has not been implemented.
"""
abstract type DistanceMethod end

subsampling(method::DistanceMethod) = method.subsampling

"""
    CentroidToCentroid <: DistanceMethod

    CentroidToCentroid()

Calculates the discrete probability of dispersal between source and destination cell
centroids. This is the naive method, but it will not handle low grid resolution well
due to severe truncation.

See: "On the approximation of continuous dispersal kernels in discrete-space models,
Joseph D. Chipperfield et al 2011"
"""
struct CentroidToCentroid <: DistanceMethod end

dispersalprob(f, ::CentroidToCentroid, x, y, cellsize) = f(sqrt(x^2 + y^2) * cellsize)

"""
    AreaToCentroid <: DistanceMethod

    AreaToCentroid(subsampling)
    AreaToCentroid(; subsampling=10.0)

Calculates the discrete probability of dispersal between source cell area and destination
centroid.

See: "On the approximation of continuous dispersal kernels in discrete-space models,
Joseph D. Chipperfield et al 2011"
"""
Base.@kwdef struct AreaToCentroid{SS<:Number} <: DistanceMethod
    "Subsampling for numerical integration"
    subsampling::SS = Param(10.0; bounds=(2.0, 40.0))
end

@inline function dispersalprob(f, dm::AreaToCentroid, x, y, cellsize)
    prob = zero(cellsize)
    centerfirst = 1 / subsampling(dm) / 2 - 0.5
    centerlast = centerfirst * -1
    range = LinRange(centerfirst, centerlast, round(Int, subsampling(dm)))
    for j in range, i in range
        prob += sqrt((x + i)^2 + (y + j)^2) * cellsize |> f
    end
    prob / subsampling(dm)^2
end

"""
    AreaToArea <: DistanceMethod

    AreaToArea(subsampling)
    AreaToArea(; subsampling=10.0)

Calculates the discrete probability of dispersal between source and destination based on
cell areas.

See: "On the approximation of continuous dispersal kernels in discrete-space models,
Joseph D. Chipperfield et al 2011"
"""
Base.@kwdef struct AreaToArea{SS<:Number} <: DistanceMethod
    subsampling::SS = Param(10.0; bounds=(2.0, 40.0))
end

@inline function dispersalprob(f, dm::AreaToArea, x, y, cellsize)
    prob = zero(cellsize)
    # Get the center point of the first cell (for both dimensions)
    centerfirst = 1 / subsampling(dm) / 2 - 0.5
    centerlast = centerfirst * -1
    range = LinRange(centerfirst, centerlast, round(Int, subsampling(dm)))
    for i in range, j in range
        for a in range, b in range
            prob += sqrt((x + i + a)^2 + (y + j + b)^2) * cellsize |> f
        end
    end
    prob / subsampling(dm)^4
end


text:
function mean(iterable)
    state = start(iterable)
    if done(iterable, state)
        error("mean of empty collection undefined: $(repr(iterable))")
    end
    count = 1
    total, state = next(iterable, state)
    while !done(iterable, state)
        value, state = next(iterable, state)
        total += value
        count += 1
    end
    return total/count
end
mean(v::AbstractArray) = sum(v) / length(v)

function mean!{T}(r::AbstractArray{T}, v::AbstractArray)
    sum!(r, v; init=true)
    rs = convert(T, length(v) / length(r))
    if rs != 1
        for i = 1:length(r)
            @inbounds r[i] /= rs
        end
    end
    return r
end

meantype{T}(::Type{T}) = typeof((zero(T) + zero(T)) / 2)
mean{T}(v::AbstractArray{T}, region) =
    mean!(Array(meantype(T), reduced_dims(size(v), region)), v)

function median!{T<:Real}(v::AbstractVector{T}; checknan::Bool=true)
    isempty(v) && error("median of an empty array is undefined")
    if checknan
        for x in v
            isnan(x) && return x
        end
    end
    n = length(v)
    if isodd(n)
        return select!(v,div(n+1,2))
    else
        m = select!(v, div(n,2):div(n,2)+1)
        return (m[1] + m[2])/2
    end
end
median{T<:Real}(v::AbstractArray{T}; checknan::Bool=true) =
    median!(vec(copy(v)), checknan=checknan)
median{T}(v::AbstractArray{T}, region; checknan::Bool=true) =
    mapslices( x->median(x; checknan=checknan), v, region )

## variances

function sumabs2_pairwise{T<:Base.LinAlg.BlasFloat}(A::StridedArray{T}, i1::Int, n::Int)
    if n <= 2048
        BLAS.dot(n, pointer(A, i1), stride(A, 1), pointer(A, i1), stride(A, 1))
    else
        n2 = div(n,2)
        sumabs2_pairwise(A, i1, n2) + sumabs2_pairwise(A, i1+n2, n-n2)
    end
end

function sumabs2_pairwise(A::AbstractArray, i1::Int, n::Int)
    if n < 256
        @inbounds s = abs2(A[i1])
        for i=i1+1:i1+n-1
            @inbounds s += abs2(A[i])
        end
        return s
    else
        n2 = div(n,2)
        return sumabs2_pairwise(A, i1, n2) + sumabs2_pairwise(A, i1+n2, n-n2)
    end
end

sumabs2(v::AbstractArray) = sumabs2_pairwise(v, 1, length(v))

plusabs2(x, y) = x + abs2(y)
eval(ngenerate(:N, :(typeof(R)), :(_sumabs2!{T,N}(R::AbstractArray, A::AbstractArray{T,N})), N->gen_reduction_body(N, plusabs2)))
sumabs2!{R}(r::AbstractArray{R}, A::AbstractArray; init::Bool=true) = _sumabs2!(initarray!(r, zero(R), init), A)
sumabs2{T}(A::AbstractArray{T}, region) = _sumabs2!(reduction_init(A, region, abs2(zero(T))+abs2(zero(T))), A)

function varzm(v::AbstractArray; corrected::Bool=true)
    n = length(v)
    n == 0 && return NaN
    return sumabs2(v) / (n - int(corrected))
end

function varzm(v::AbstractArray, region; corrected::Bool=true)
    cn = regionsize(v, region) - int(corrected)
    sumabs2(v, region) / cn
end

function varm_pairwise(A::AbstractArray, m::Number, i1::Int, n::Int) # see sum_pairwise
    if n < 256
        @inbounds s = abs2(A[i1] - m)
        for i = i1+1:i1+n-1
            @inbounds s += abs2(A[i] - m)
        end
        return s
    else
        n2 = div(n,2)
        return varm_pairwise(A, m, i1, n2) + varm_pairwise(A, m, i1+n2, n-n2)
    end
end

function varm(v::AbstractArray, m::Number; corrected::Bool=true)
    n = length(v)
    n == 0 && return NaN
    return varm_pairwise(v, m, 1, n) / (n - int(corrected))
end

@ngenerate N Array{typeof((abs2(zero(T))+abs2(zero(T)))/1), N} function _varm{S,T,N}(v::AbstractArray{S,N}, m::AbstractArray{T,N}, region, corrected::Bool)
    rdims = reduced_dims(v, region)
    rdims == size(m) || error(DimensionMismatch("size of mean does not match reduced dimensions"))

    R = fill!(similar(v, typeof((abs2(zero(T))+abs2(zero(T)))/1), rdims), 0)
    @nextract N sizeR d->size(R,d)
    @nloops N i v d->(j_d = sizeR_d==1 ? 1 : i_d) begin
        @inbounds (@nref N R j) += abs2((@nref N v i) - (@nref N m j))
    end
    scale!(R, 1/(regionsize(v, region) - int(corrected)))
end
varm{S,T,N}(v::AbstractArray{S,N}, m::AbstractArray{T,N}, region; corrected::Bool=true) = _varm(v, m, region, corrected)

function var(v::AbstractArray; corrected::Bool=true, mean=nothing)
    mean == 0 ? varzm(v; corrected=corrected) :
    mean == nothing ? varm(v, Base.mean(v); corrected=corrected) :
    isa(mean, Number) ? varm(v, mean; corrected=corrected) :
    error("Invalid value of mean.")
end

function var(v::AbstractArray, region; corrected::Bool=true, mean=nothing)
    mean == 0 ? varzm(v, region; corrected=corrected) :
    mean == nothing ? varm(v, Base.mean(v, region), region; corrected=corrected) :
    isa(mean, AbstractArray) ? varm(v, mean, region; corrected=corrected) :
    error("Invalid value of mean.")
end

function var(iterable; corrected::Bool=true, mean=nothing)
    state = start(iterable)
    if done(iterable, state)
        error("variance of empty collection undefined: $(repr(iterable))")
    end
    count = 1
    value, state = next(iterable, state)
    if mean == nothing
        # Use Welford algorithm as seen in (among other places)
        # Knuth's TAOCP, Vol 2, page 232, 3rd edition.
        M = value / 1
        S = zero(M)
        while !done(iterable, state)
            value, state = next(iterable, state)
            count += 1
            new_M = M + (value - M) / count
            S = S + (value - M) * (value - new_M)
            M = new_M
        end
        return S / (count - int(corrected))
    else # mean provided
        # Cannot use a compensated version, e.g. the one from
        # "Updating Formulae and a Pairwise Algorithm for Computing Sample Variances."
        # by Chan, Golub, and LeVeque, Technical Report STAN-CS-79-773,
        # Department of Computer Science, Stanford University,
        # because user can provide mean value that is different to mean(iterable)
        sum2 = (value - mean)^2
        while !done(iterable, state)
            value, state = next(iterable, state)
            count += 1
            sum2 += (value - mean)^2
        end
        return sum2 / (count - int(corrected))
    end
end

varm(iterable, m::Number; corrected::Bool=true) =
    var(iterable, corrected=corrected, mean=m)

## variances over ranges

varm(v::Range, m::Number) = var(v)

function var(v::Range)
    s = step(v)
    l = length(v)
    if l == 0 || l == 1
        return NaN
    end
    return abs2(s) * (l + 1) * l / 12
end

## standard deviation

function sqrt!(v::AbstractArray)
    for i = 1:length(v)
        v[i] = sqrt(v[i])
    end
    v
end

stdm(v::AbstractArray, m::Number; corrected::Bool=true) =
    sqrt(varm(v, m; corrected=corrected))

std(v::AbstractArray; corrected::Bool=true, mean=nothing) =
    sqrt(var(v; corrected=corrected, mean=mean))

std(v::AbstractArray, region; corrected::Bool=true, mean=nothing) =
    sqrt!(var(v, region; corrected=corrected, mean=mean))

std(iterable; corrected::Bool=true, mean=nothing) =
    sqrt(var(iterable, corrected=corrected, mean=mean))

stdm(iterable, m::Number; corrected::Bool=true) =
    std(iterable, corrected=corrected, mean=m)

## pearson covariance functions ##

# auxiliary functions

_conj{T<:Real}(x::AbstractArray{T}) = x
_conj(x::AbstractArray) = conj(x)

_getnobs(x::AbstractVector, vardim::Int) = length(x)
_getnobs(x::AbstractMatrix, vardim::Int) = size(x, vardim)

function _getnobs(x::AbstractVecOrMat, y::AbstractVecOrMat, vardim::Int)
    n = _getnobs(x, vardim)
    _getnobs(y, vardim) == n || throw(DimensionMismatch("Dimensions of x and y mismatch."))
    return n
end

_vmean(x::AbstractVector, vardim::Int) = mean(x)
_vmean(x::AbstractMatrix, vardim::Int) = mean(x, vardim)

# core functions

unscaled_covzm(x::AbstractVector) = dot(x, x)
unscaled_covzm(x::AbstractMatrix, vardim::Int) = (vardim == 1 ? _conj(x'x) : x * x')

unscaled_covzm(x::AbstractVector, y::AbstractVector) = dot(x, y)
unscaled_covzm(x::AbstractVector, y::AbstractMatrix, vardim::Int) =
    (vardim == 1 ? At_mul_B(x, _conj(y)) : At_mul_Bt(x, _conj(y)))
unscaled_covzm(x::AbstractMatrix, y::AbstractVector, vardim::Int) =
    (c = vardim == 1 ? At_mul_B(x, _conj(y)) :  x * _conj(y); reshape(c, length(c), 1))
unscaled_covzm(x::AbstractMatrix, y::AbstractMatrix, vardim::Int) =
    (vardim == 1 ? At_mul_B(x, _conj(y)) : A_mul_Bc(x, y))

# covzm (with centered data)

covzm(x::AbstractVector; corrected::Bool=true) = unscaled_covzm(x, x) / (length(x) - int(corrected))

covzm(x::AbstractMatrix; vardim::Int=1, corrected::Bool=true) =
    scale!(unscaled_covzm(x, vardim), inv(size(x,vardim) - int(corrected)))

covzm(x::AbstractVector, y::AbstractVector; corrected::Bool=true) =
    unscaled_covzm(x, y) / (length(x) - int(corrected))

covzm(x::AbstractVecOrMat, y::AbstractVecOrMat; vardim::Int=1, corrected::Bool=true) =
    scale!(unscaled_covzm(x, y, vardim), inv(_getnobs(x, y, vardim) - int(corrected)))

# covm (with provided mean)

covm(x::AbstractVector, xmean; corrected::Bool=true) =
    covzm(x .- xmean; corrected=corrected)

covm(x::AbstractMatrix, xmean; vardim::Int=1, corrected::Bool=true) =
    covzm(x .- xmean; vardim=vardim, corrected=corrected)

covm(x::AbstractVector, xmean, y::AbstractVector, ymean; corrected::Bool=true) =
    covzm(x .- xmean, y .- ymean; corrected=corrected)

covm(x::AbstractVecOrMat, xmean, y::AbstractVecOrMat, ymean; vardim::Int=1, corrected::Bool=true) =
    covzm(x .- xmean, y .- ymean; vardim=vardim, corrected=corrected)

# cov (API)

function cov(x::AbstractVector; corrected::Bool=true, mean=nothing)
    mean == 0 ? covzm(x; corrected=corrected) :
    mean == nothing ? covm(x, Base.mean(x); corrected=corrected) :
    isa(mean, Number) ? covm(x, mean; corrected=corrected) :
    error("Invalid value of mean.")
end

function cov(x::AbstractMatrix; vardim::Int=1, corrected::Bool=true, mean=nothing)
    mean == 0 ? covzm(x; vardim=vardim, corrected=corrected) :
    mean == nothing ? covm(x, _vmean(x, vardim); vardim=vardim, corrected=corrected) :
    isa(mean, AbstractArray) ? covm(x, mean; vardim=vardim, corrected=corrected) :
    error("Invalid value of mean.")
end

function cov(x::AbstractVector, y::AbstractVector; corrected::Bool=true, mean=nothing)
    mean == 0 ? covzm(x, y; corrected=corrected) :
    mean == nothing ? covm(x, Base.mean(x), y, Base.mean(y); corrected=corrected) :
    isa(mean, (Number,Number)) ? covm(x, mean[1], y, mean[2]; corrected=corrected) :
    error("Invalid value of mean.")
end

function cov(x::AbstractVecOrMat, y::AbstractVecOrMat; vardim::Int=1, corrected::Bool=true, mean=nothing)
    if mean == 0
        covzm(x, y; vardim=vardim, corrected=corrected)
    elseif mean == nothing
        covm(x, _vmean(x, vardim), y, _vmean(y, vardim); vardim=vardim, corrected=corrected)
    elseif isa(mean, (Any,Any))
        covm(x, mean[1], y, mean[2]; vardim=vardim, corrected=corrected)
    else
        error("Invalid value of mean.")
    end
end

# cov2cor!

function cov2cor!{T}(C::AbstractMatrix{T}, xsd::AbstractArray)
    nx = length(xsd)
    size(C) == (nx, nx) || throw(DimensionMismatch("Inconsistent dimensions."))
    for j = 1:nx
        for i = 1:j-1
            C[i,j] = C[j,i]
        end
        C[j,j] = one(T)
        for i = j+1:nx
            C[i,j] /= (xsd[i] * xsd[j])
        end
    end
    return C
end

function cov2cor!(C::AbstractMatrix, xsd::Number, ysd::AbstractArray)
    nx, ny = size(C)
    length(ysd) == ny || throw(DimensionMismatch("Inconsistent dimensions."))
    for j = 1:ny
        for i = 1:nx
            C[i,j] /= (xsd * ysd[j])
        end
    end
    return C
end

function cov2cor!(C::AbstractMatrix, xsd::AbstractArray, ysd::Number)
    nx, ny = size(C)
    length(xsd) == nx || throw(DimensionMismatch("Inconsistent dimensions."))
    for j = 1:ny
        for i = 1:nx
            C[i,j] /= (xsd[i] * ysd)
        end
    end
    return C
end

function cov2cor!(C::AbstractMatrix, xsd::AbstractArray, ysd::AbstractArray)
    nx, ny = size(C)
    (length(xsd) == nx && length(ysd) == ny) ||
        throw(DimensionMismatch("Inconsistent dimensions."))
    for j = 1:ny
        for i = 1:nx
            C[i,j] /= (xsd[i] * ysd[j])
        end
    end
    return C
end

# # corzm (non-exported, with centered data)

corzm{T}(x::AbstractVector{T}) = float(one(T) * one(T))

corzm(x::AbstractMatrix; vardim::Int=1) =
    (c = unscaled_covzm(x, vardim); cov2cor!(c, sqrt!(diag(c))))

function corzm(x::AbstractVector, y::AbstractVector)
    n = length(x)
    length(y) == n || throw(DimensionMismatch("Inconsistent lengths."))
    x1 = x[1]
    y1 = y[1]
    xx = abs2(x1)
    yy = abs2(y1)
    xy = x1 * conj(y1)
    i = 1
    while i < n
        i += 1
        @inbounds xi = x[i]
        @inbounds yi = y[i]
        xx += abs2(xi)
        yy += abs2(yi)
        xy += xi * conj(yi)
    end
    return xy / (sqrt(xx) * sqrt(yy))
end

corzm(x::AbstractVector, y::AbstractMatrix; vardim::Int=1) =
    cov2cor!(unscaled_covzm(x, y, vardim), sqrt(sumabs2(x)), sqrt!(sumabs2(y, vardim)))

corzm(x::AbstractMatrix, y::AbstractVector; vardim::Int=1) =
    cov2cor!(unscaled_covzm(x, y, vardim), sqrt!(sumabs2(x, vardim)), sqrt(sumabs2(y)))

corzm(x::AbstractMatrix, y::AbstractMatrix; vardim::Int=1) =
    cov2cor!(unscaled_covzm(x, y, vardim), sqrt!(sumabs2(x, vardim)), sqrt!(sumabs2(y, vardim)))

# corm

corm(x::AbstractVector, xmean) = corzm(x .- xmean)

corm(x::AbstractMatrix, xmean; vardim::Int=1) = corzm(x .- xmean; vardim=vardim)

corm(x::AbstractVector, xmean, y::AbstractVector, ymean) = corzm(x .- xmean, y .- ymean)

corm(x::AbstractVecOrMat, xmean, y::AbstractVecOrMat, ymean; vardim::Int=1) =
    corzm(x .- xmean, y .- ymean; vardim=vardim)

# cor

function cor(x::AbstractVector; mean=nothing)
    mean == 0 ? corzm(x) :
    mean == nothing ? corm(x, Base.mean(x)) :
    isa(mean, Number) ? corm(x, mean) :
    error("Invalid value of mean.")
end

function cor(x::AbstractMatrix; vardim::Int=1, mean=nothing)
    mean == 0 ? corzm(x; vardim=vardim) :
    mean == nothing ? corm(x, _vmean(x, vardim); vardim=vardim) :
    isa(mean, AbstractArray) ? corm(x, mean; vardim=vardim) :
    error("Invalid value of mean.")
end

function cor(x::AbstractVector, y::AbstractVector; mean=nothing)
    mean == 0 ? corzm(x, y) :
    mean == nothing ? corm(x, Base.mean(x), y, Base.mean(y)) :
    isa(mean, (Number,Number)) ? corm(x, mean[1], y, mean[2]) :
    error("Invalid value of mean.")
end

function cor(x::AbstractVecOrMat, y::AbstractVecOrMat; vardim::Int=1, mean=nothing)
    if mean == 0
        corzm(x, y; vardim=vardim)
    elseif mean == nothing
        corm(x, _vmean(x, vardim), y, _vmean(y, vardim); vardim=vardim)
    elseif isa(mean, (Any,Any))
        corm(x, mean[1], y, mean[2]; vardim=vardim)
    else
        error("Invalid value of mean.")
    end
end

## nice-valued ranges for histograms

function histrange{T<:FloatingPoint,N}(v::AbstractArray{T,N}, n::Integer)
    if length(v) == 0
        return 0.0:1.0:0.0
    end
    lo, hi = extrema(v)
    if hi == lo
        step = 1.0
    else
        bw = (hi - lo) / n
        e = 10.0^floor(log10(bw))
        r = bw / e
        if r <= 2
            step = 2*e
        elseif r <= 5
            step = 5*e
        else
            step = 10*e
        end
    end
    start = step*(ceil(lo/step)-1)
    nm1 = iceil((hi - start)/step)
    start:step:(start + nm1*step)
end

function histrange{T<:Integer,N}(v::AbstractArray{T,N}, n::Integer)
    if length(v) == 0
        return 0:1:0
    end
    lo, hi = extrema(v)
    if hi == lo
        step = 1
    else
        bw = (hi - lo) / n
        e = 10^max(0,ifloor(log10(bw)))
        r = bw / e
        if r <= 1
            step = e
        elseif r <= 2
            step = 2*e
        elseif r <= 5
            step = 5*e
        else
            step = 10*e
        end
    end
    start = step*(ceil(lo/step)-1)
    nm1 = iceil((hi - start)/step)
    start:step:(start + nm1*step)
end

## midpoints of intervals
midpoints(r::Range) = r[1:length(r)-1] + 0.5*step(r)
midpoints(v::AbstractVector) = [0.5*(v[i] + v[i+1]) for i in 1:length(v)-1]

## hist ##
function sturges(n)  # Sturges' formula
    n==0 && return one(n)
    iceil(log2(n))+1
end

function hist!{HT}(h::AbstractArray{HT}, v::AbstractVector, edg::AbstractVector; init::Bool=true)
    n = length(edg) - 1
    length(h) == n || error("length(h) must equal length(edg) - 1.")
    if init
        fill!(h, zero(HT))
    end
    for x in v
        i = searchsortedfirst(edg, x)-1
        if 1 <= i <= n
            h[i] += 1
        end
    end
    edg, h
end

hist(v::AbstractVector, edg::AbstractVector) = hist!(Array(Int, length(edg)-1), v, edg)
hist(v::AbstractVector, n::Integer) = hist(v,histrange(v,n))
hist(v::AbstractVector) = hist(v,sturges(length(v)))

function hist!{HT}(H::AbstractArray{HT,2}, A::AbstractMatrix, edg::AbstractVector; init::Bool=true)
    m, n = size(A)
    size(H) == (length(edg)-1, n) || error("Incorrect size of H.")
    if init
        fill!(H, zero(HT))
    end
    for j = 1:n
        hist!(sub(H, :, j), sub(A, :, j), edg)
    end
    edg, H
end

hist(A::AbstractMatrix, edg::AbstractVector) = hist!(Array(Int, length(edg)-1, size(A,2)), A, edg)
hist(A::AbstractMatrix, n::Integer) = hist(A,histrange(A,n))
hist(A::AbstractMatrix) = hist(A,sturges(size(A,1)))

## hist2d
function hist2d!{HT}(H::AbstractArray{HT,2}, v::AbstractMatrix,
                     edg1::AbstractVector, edg2::AbstractVector; init::Bool=true)
    size(v,2) == 2 || error("hist2d requires an Nx2 matrix.")
    n = length(edg1) - 1
    m = length(edg2) - 1
    size(H) == (n, m) || error("Incorrect size of H.")
    if init
        fill!(H, zero(HT))
    end
    for i = 1:size(v,1)
        x = searchsortedfirst(edg1, v[i,1]) - 1
        y = searchsortedfirst(edg2, v[i,2]) - 1
        if 1 <= x <= n && 1 <= y <= m
            @inbounds H[x,y] += 1
        end
    end
    edg1, edg2, H
end

hist2d(v::AbstractMatrix, edg1::AbstractVector, edg2::AbstractVector) =
    hist2d!(Array(Int, length(edg1)-1, length(edg2)-1), v, edg1, edg2)

hist2d(v::AbstractMatrix, edg::AbstractVector) = hist2d(v, edg, edg)

hist2d(v::AbstractMatrix, n1::Integer, n2::Integer) =
    hist2d(v, histrange(sub(v,:,1),n1), histrange(sub(v,:,2),n2))
hist2d(v::AbstractMatrix, n::Integer) = hist2d(v, n, n)
hist2d(v::AbstractMatrix) = hist2d(v, sturges(size(v,1)))

## quantiles ##

# for now, use the R/S definition of quantile; may want variants later
# see ?quantile in R -- this is type 7
function quantile!(v::AbstractVector, q::AbstractVector)
    isempty(v) && error("empty data array")
    isempty(q) && error("empty quantile array")

    # make sure the quantiles are in [0,1]
    q = bound_quantiles(q)

    lv = length(v)
    lq = length(q)

    index = 1 .+ (lv-1)*q
    lo = ifloor(index)
    hi = iceil(index)
    sort!(v)
    isnan(v[end]) && error("quantiles are undefined in presence of NaNs")
    i = find(index .> lo)
    r = float(v[lo])
    h = (index.-lo)[i]
    r[i] = (1.-h).*r[i] + h.*v[hi[i]]
    return r
end
quantile(v::AbstractVector, q::AbstractVector) = quantile!(copy(v),q)
quantile(v::AbstractVector, q::Number) = quantile(v,[q])[1]

function bound_quantiles(qs::AbstractVector)
    epsilon = 100*eps()
    if (any(qs .< -epsilon) || any(qs .> 1+epsilon))
        error("quantiles out of [0,1] range")
    end
    [min(1,max(0,q)) for q = qs]
end

max_stars_repo_name: aplund/julia

max_issues_repo_name: aplund/julia

max_forks_repo_name: aplund/julia


text:
# Definitions borrow from https://github.com/FluxML/NNlib.jl (names used to
# standardize package). TODO: Decide how/if to incorporate NNlib depedency into
# McCormick.jl/EAGO.jl.
oftf(x, y) = oftype(float(x), y)
leakyrelu(x, a=oftf(x, 0.01)) = max(a * x, x)
swish1(x) = x * sigmoid(x)

"""
    AtomType
"""
@enum(AtomType, VAR_ATOM, PARAM_ATOM, CONST_ATOM, SELECT_ATOM,
                ABS, ABS2, INV, RAD2DEG, DEG2RAD, ONE, ZERO,
                MULT, DIV, PLUS, MINUS, POW, MIN, MAX, STEP, SIGN,
                LOG, LOG2, LOG10, LOG1P, EXP, EXP2, EXP10, EXPM1,
                SIN, COS, TAN, CSC, SEC, COT,
                ASIN, ACOS, ATAN, ACSC, ASEC, ACOT,
                SINH, COSH, TANH, CSCH, SECH, COTH,
                ASINH, ACOSH, ATANH, ACSCH, ASECH, ACOTH,
                ERF, ERFC, ERFINV, ERFCINV, SQRT, CBRT,
                # ML specific functions
                RELU, LEAKYRELU, SOFTPLUS, SOFTSIGN, GELU, SILU, SIGMOID,
                # EAGO modeling functions
                XLOGX, ARH,
                # EAGO bound enforcing functions
                POS, NEG, LOWER_BND, UPPER_BND, BND,
                # Atoms for nonexpression types
                USER, USERN, SUBEXPR)

# Define univariate functions that correspond to the AtomType
const UNIVARIATE_ATOM_DICT = Dict{AtomType, Symbol}()

UNIVARIATE_ATOM_DICT[PLUS]      = :+
UNIVARIATE_ATOM_DICT[MINUS]     = :-
UNIVARIATE_ATOM_DICT[INV]       = :inv
UNIVARIATE_ATOM_DICT[ONE]       = :one
UNIVARIATE_ATOM_DICT[ZERO]      = :zero
UNIVARIATE_ATOM_DICT[ABS]       = :abs
UNIVARIATE_ATOM_DICT[ABS2]      = :abs2
UNIVARIATE_ATOM_DICT[RAD2DEG]   = :rad2deg
UNIVARIATE_ATOM_DICT[DEG2RAD]   = :deg2rad
UNIVARIATE_ATOM_DICT[STEP]      = :step
UNIVARIATE_ATOM_DICT[SIGN]      = :sign

UNIVARIATE_ATOM_DICT[LOG]       = :log
UNIVARIATE_ATOM_DICT[LOG2]      = :log2
UNIVARIATE_ATOM_DICT[LOG10]     = :log10
UNIVARIATE_ATOM_DICT[LOG1P]     = :log1p
UNIVARIATE_ATOM_DICT[EXP]       = :exp
UNIVARIATE_ATOM_DICT[EXP2]      = :exp2
UNIVARIATE_ATOM_DICT[EXP10]     = :exp10
UNIVARIATE_ATOM_DICT[EXPM1]     = :expm1

UNIVARIATE_ATOM_DICT[SIN]       = :sin
UNIVARIATE_ATOM_DICT[COS]       = :cos
UNIVARIATE_ATOM_DICT[TAN]       = :tan
UNIVARIATE_ATOM_DICT[CSC]       = :csc
UNIVARIATE_ATOM_DICT[SEC]       = :sec
UNIVARIATE_ATOM_DICT[COT]       = :cot
UNIVARIATE_ATOM_DICT[ASIN]      = :asin
UNIVARIATE_ATOM_DICT[ACOS]      = :acos
UNIVARIATE_ATOM_DICT[ATAN]      = :atan
UNIVARIATE_ATOM_DICT[ACSC]      = :acsc
UNIVARIATE_ATOM_DICT[ASEC]      = :asec
UNIVARIATE_ATOM_DICT[ACOT]      = :acot
UNIVARIATE_ATOM_DICT[SINH]      = :sinh
UNIVARIATE_ATOM_DICT[COSH]      = :cosh
UNIVARIATE_ATOM_DICT[TANH]      = :tanh
UNIVARIATE_ATOM_DICT[CSCH]      = :csch
UNIVARIATE_ATOM_DICT[SECH]      = :sech
UNIVARIATE_ATOM_DICT[COTH]      = :coth
UNIVARIATE_ATOM_DICT[ASINH]     = :asinh
UNIVARIATE_ATOM_DICT[ACOSH]     = :acosh
UNIVARIATE_ATOM_DICT[ATANH]     = :atanh
UNIVARIATE_ATOM_DICT[ACSCH]     = :acsch
UNIVARIATE_ATOM_DICT[ASECH]     = :asech
UNIVARIATE_ATOM_DICT[ACOTH]     = :acoth

UNIVARIATE_ATOM_DICT[ERF]       = :erf
UNIVARIATE_ATOM_DICT[ERFC]      = :erfc
UNIVARIATE_ATOM_DICT[ERFINV]    = :erfinv
UNIVARIATE_ATOM_DICT[ERFCINV]   = :erfcinv

UNIVARIATE_ATOM_DICT[SQRT]      = :sqrt
UNIVARIATE_ATOM_DICT[CBRT]      = :cbrt

UNIVARIATE_ATOM_DICT[POS]       = :positive
UNIVARIATE_ATOM_DICT[NEG]       = :negative
UNIVARIATE_ATOM_DICT[USER]      = :user

UNIVARIATE_ATOM_DICT[RELU]       = :relu
UNIVARIATE_ATOM_DICT[LEAKYRELU]  = :leakyrelu
UNIVARIATE_ATOM_DICT[SIGMOID]    = :sigmoid
UNIVARIATE_ATOM_DICT[GELU]       = :gelu
UNIVARIATE_ATOM_DICT[SOFTPLUS]   = :softplus
UNIVARIATE_ATOM_DICT[SOFTSIGN]   = :softsign
UNIVARIATE_ATOM_DICT[SILU]       = :swish1

UNIVARIATE_ATOM_DICT[XLOGX]      = :xlogx

# Define bivariate functions that correspond to the AtomType
const BIVARIATE_ATOM_DICT = Dict{AtomType, Symbol}()
BIVARIATE_ATOM_DICT[DIV]       = :/
BIVARIATE_ATOM_DICT[POW]       = :^
BIVARIATE_ATOM_DICT[LOWER_BND] = :lower_bnd
BIVARIATE_ATOM_DICT[UPPER_BND] = :upper_bnd
BIVARIATE_ATOM_DICT[ARH]       = :arh
BIVARIATE_ATOM_DICT[MINUS]     = :-
BIVARIATE_ATOM_DICT[MIN]       = :min
BIVARIATE_ATOM_DICT[MAX]       = :max
BIVARIATE_ATOM_DICT[PLUS]      = :+
BIVARIATE_ATOM_DICT[MULT]      = :*
#BIVARIATE_ATOM_DICT[EXPAX]     = :lower_bnd
#BIVARIATE_ATOM_DICT[EXPXY]     = :lower_bnd

# Define n-arity functions that correspond to the AtomType
const NARITY_ATOM_DICT = Dict{AtomType, Symbol}()
NARITY_ATOM_DICT[MIN]  = :min
NARITY_ATOM_DICT[MAX]  = :max
NARITY_ATOM_DICT[PLUS] = :+
NARITY_ATOM_DICT[MULT] = :*
NARITY_ATOM_DICT[USERN]   = :usern

const ALL_ATOM_DICT = Dict{AtomType, Symbol}()
foreach(x -> setindex!(ALL_ATOM_DICT, x[2], x[1]), UNIVARIATE_ATOM_DICT)
foreach(x -> setindex!(ALL_ATOM_DICT, x[2], x[1]), BIVARIATE_ATOM_DICT)
foreach(x -> setindex!(ALL_ATOM_DICT, x[2], x[1]), NARITY_ATOM_DICT)

# A functions that may be 1 to n-arity functions that correspond to the AtomType
ALL_ATOM_DICT[BND]       = :bnd
#ATM_EVAL[QUAD1]     = :quad1
#ATM_EVAL[QUAD2]     = :quad2

# List of keys only
const UNIVARIATE_ATOM_TYPES = AtomType[k for k in keys(UNIVARIATE_ATOM_DICT)]
const BIVARIATE_ATOM_TYPES = AtomType[k for k in keys(BIVARIATE_ATOM_DICT)]
const NARITY_ATOM_TYPES = AtomType[k for k in keys(NARITY_ATOM_DICT)]
const ALL_ATOM_TYPES = AtomType[k for k in keys(ALL_ATOM_DICT)]

# Reverse lookup dicts
const REV_UNIVARIATE_ATOM_DICT = Dict(UNIVARIATE_ATOM_DICT[k] => k for k in keys(UNIVARIATE_ATOM_DICT))
const REV_BIVARIATE_ATOM_DICT = Dict(BIVARIATE_ATOM_DICT[k] => k for k in keys(BIVARIATE_ATOM_DICT))
const REV_NARITY_ATOM_DICT = Dict(NARITY_ATOM_DICT[k] => k for k in keys(NARITY_ATOM_DICT))

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     cell_metadata_json: true
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# # Meridional Overturning Circulation
#
# Global Ocean transport depictions often involve concepts like the [AMOC](https://en.wikipedia.org/wiki/Atlantic_meridional_overturning_circulation), [Thermohaline Circulation](https://en.wikipedia.org/wiki/Thermohaline_circulation), or [conveyor Belt](http://oceanmotion.org/html/background/ocean-conveyor-belt.html). To apply these concepts to gridded ocean and climate models, one can compute an _overturning streamfunction_ as shown below. For more detail, please refer to [Forget et al, 2015](https://doi.org/10.5194/gmd-8-3071-2015) and the other [GlobalOceanNotebooks](https://github.com/JuliaClimate/GlobalOceanNotebooks).

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ### Read grid & velocities from file
#
# 1. pre-requisites
# 2. read variables
# 3. conversion to transports

# +
include("helper_functions.jl")

pth=MeshArrays.GRID_LLC90
γ=GridSpec("LatLonCap",pth)
Γ=GridLoad(γ;option="full")
LC=LatitudeCircles(-89.0:89.0,Γ);

# +
using OceanStateEstimation
OceanStateEstimation.get_ecco_velocity_if_needed();
#OceanStateEstimation.get_occa_velocity_if_needed();

"""
    read_velocities(γ::gcmgrid,t::Int,pth::String)

Read velocity components `u,v` from files in `pth`for time `t`
"""
function read_velocities(γ::gcmgrid,t::Int,pth::String)
    u=read_nctiles("$pth"*"UVELMASS/UVELMASS","UVELMASS",γ,I=(:,:,:,t))
    v=read_nctiles("$pth"*"VVELMASS/VVELMASS","VVELMASS",γ,I=(:,:,:,t))
    return u,v
end

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ### Compute Overturning Streamfunction
#
# 1. integrate along latitude circles
# 2. integrate from the bottom

# + {"slideshow": {"slide_type": "-"}}
nz=size(Γ.hFacC,2); nt=12; nl=length(LC)
ov=Array{Float64,3}(undef,nl,nz,nt)

#integrate across latitude circles
for t=1:nt
    (U,V)=read_velocities(Γ.XC.grid,t,ECCOclim_path)
    (U,V)=convert_velocities(U,V,Γ)
    for z=1:nz
        UV=Dict("U"=>U[:,z],"V"=>V[:,z],"dimensions"=>["x","y"])
        [ov[l,z,t]=ThroughFlow(UV,LC[l],Γ) for l=1:nl]
    end
end

#integrate from bottom
ov=reverse(cumsum(reverse(ov,dims=2),dims=2),dims=2);

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ### Plot Annual Mean And Variability
# -

x=vec(-89.0:89.0); y=reverse(vec(Γ.RF[1:end-1])); #coordinate variables

# + {"slideshow": {"slide_type": "fragment"}, "cell_style": "split"}
tmp=dropdims(mean(ov,dims=3),dims=3)
z=reverse(tmp,dims=2); z[z.==0.0].=NaN

contourf(x,y,1e-6*transpose(z),clims=(-40,40),
    title="Overturning mean (Eulerian; in Sv)",titlefontsize=10)
#savefig("MOC_mean.png")
# + {"slideshow": {"slide_type": "fragment"}, "cell_style": "split"}
tmp=dropdims(std(ov,dims=3),dims=3)
z=reverse(tmp,dims=2); z[z.==0.0].=NaN

contourf(x,y,1e-6*transpose(z),clims=(-40,40),
    title="Overturning standard deviation (Eulerian; in Sv)",titlefontsize=10)
#savefig("MOC_std.png")
# -

max_stars_repo_name: gaelforget/Notebooks

max_issues_repo_name: gaelforget/Notebooks

max_forks_repo_name: gaelforget/Notebooks


text:
using Decomp
using Base.Test

a = zeros(3,3)
for i=1:100
  a[1,1] = rand(-2.:10e-8:2.)
  a[2,2] = rand(-2.:10e-8:2.)
  a[3,3] = rand(-2.:10e-8:2.)
  a[1,2] = rand(-2.:10e-8:2.)
  a[1,3] = rand(-2.:10e-8:2.)
  a[2,3] = rand(-2.:10e-8:2.)
  a[2,1] = a[1,2]
  a[3,1] = a[1,3]
  a[3,2] = a[2,3]
  eigv,eigvec1,eigvec2,eigvec3 = eigen(a)
  eigvc,eigvec123 = eig(a)
  @test eigv ≈ sort!(eigvc,by=abs,rev=true)
  @test eigv[1] * eigvec1*eigvec1' .+ eigv[2] * eigvec2*eigvec2' .+ eigv[3] * eigvec3*eigvec3' ≈ a
end


text:
@testset "eternity2.jl" begin
    @testset "fill_with_generator!" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        n=6

        generator = SeaPearl.Eternity2Generator(6,6,6)

        for _ in 1:3
            SeaPearl.fill_with_generator!(model, generator)
            @test length(model.variables) == 120
            @test length(model.constraints) == 61
            empty!(model)
        end
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Represents a Quobyte mount that lasts the lifetime of a pod. Quobyte volumes do not support ownership management or SELinux relabeling.

    IoK8sApiCoreV1QuobyteVolumeSource(;
        group=nothing,
        readOnly=nothing,
        registry=nothing,
        tenant=nothing,
        user=nothing,
        volume=nothing,
    )

    - group::String : Group to map volume access to Default is no group
    - readOnly::Bool : ReadOnly here will force the Quobyte volume to be mounted with read-only permissions. Defaults to false.
    - registry::String : Registry represents a single or multiple Quobyte Registry services specified as a string as host:port pair (multiple entries are separated with commas) which acts as the central registry for volumes
    - tenant::String : Tenant owning the given Quobyte volume in the Backend Used with dynamically provisioned Quobyte volumes, value is set by the plugin
    - user::String : User to map volume access to Defaults to serivceaccount user
    - volume::String : Volume is a string that references an already created Quobyte volume by name.
"""
mutable struct IoK8sApiCoreV1QuobyteVolumeSource <: SwaggerModel
    group::Any # spec type: Union{ Nothing, String } # spec name: group
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    registry::Any # spec type: Union{ Nothing, String } # spec name: registry
    tenant::Any # spec type: Union{ Nothing, String } # spec name: tenant
    user::Any # spec type: Union{ Nothing, String } # spec name: user
    volume::Any # spec type: Union{ Nothing, String } # spec name: volume

    function IoK8sApiCoreV1QuobyteVolumeSource(;group=nothing, readOnly=nothing, registry=nothing, tenant=nothing, user=nothing, volume=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1QuobyteVolumeSource, Symbol("group"), group)
        setfield!(o, Symbol("group"), group)
        validate_property(IoK8sApiCoreV1QuobyteVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1QuobyteVolumeSource, Symbol("registry"), registry)
        setfield!(o, Symbol("registry"), registry)
        validate_property(IoK8sApiCoreV1QuobyteVolumeSource, Symbol("tenant"), tenant)
        setfield!(o, Symbol("tenant"), tenant)
        validate_property(IoK8sApiCoreV1QuobyteVolumeSource, Symbol("user"), user)
        setfield!(o, Symbol("user"), user)
        validate_property(IoK8sApiCoreV1QuobyteVolumeSource, Symbol("volume"), volume)
        setfield!(o, Symbol("volume"), volume)
        o
    end
end # type IoK8sApiCoreV1QuobyteVolumeSource

const _property_map_IoK8sApiCoreV1QuobyteVolumeSource = Dict{Symbol,Symbol}(Symbol("group")=>Symbol("group"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("registry")=>Symbol("registry"), Symbol("tenant")=>Symbol("tenant"), Symbol("user")=>Symbol("user"), Symbol("volume")=>Symbol("volume"))
const _property_types_IoK8sApiCoreV1QuobyteVolumeSource = Dict{Symbol,String}(Symbol("group")=>"String", Symbol("readOnly")=>"Bool", Symbol("registry")=>"String", Symbol("tenant")=>"String", Symbol("user")=>"String", Symbol("volume")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1QuobyteVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1QuobyteVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1QuobyteVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1QuobyteVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1QuobyteVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1QuobyteVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1QuobyteVolumeSource)
    (getproperty(o, Symbol("registry")) === nothing) && (return false)
    (getproperty(o, Symbol("volume")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1QuobyteVolumeSource }, name::Symbol, val)
end


text:
"""
Created on Wed Sep 23 2020

@author: Yoann Pradat

Multiplicative update algorithms for α divergence.
"""

"""
    _update_nmf_mu_α!(V, WH, W, H, α, ϵ)

"""
function _update_nmf_mu_α!(V::Matrix{T}, WH::Matrix{T}, W::Matrix{T}, H::Matrix{T}, α::Real, ϵ::Real) where T <: Real
    if α == 0
        #### dual Kullback-Leibler

        #### update H
        #### NOTE: if W has columns normalized to sum to 1, the normalization matrix
        #### (transpose(W) * ones(size(V))) .^ (-1) is just ones(K, N)
        H .= H .* exp.((transpose(W) * ones(size(V))) .^ (-1)  .* (transpose(W) * log.(V ./ (WH .+ ϵ))))
        WH .= W * H

        #### update W
        #### NOTE: if H has columns normalized to sum to 1, the normalization matrix
        #### (transpose(H) * ones(size(V))) .^ (-1) is just ones(K, N)
        W .= W .* exp.((ones(size(V)) * transpose(H)) .^ (-1) .* (log.(V ./ (WH .+ ϵ)) * transpose(H)))
        WH .= W * H

    else
        #### general α divergence updates for α != 0

        #### A.Cichocki, S.Amari, R.Zdunek, R.Kompass, G.Hori, and Z.He,
        #### "Extended SMART algorithms for non-negative matrix factorization,"
        #### Artificial Intelligence and Soft Computing, 2006.
        #### Equation (44)

        #### note: it is not clear in eqs (45) and (46) how the normalization
        #### made the denominators disappear in (44).

        # #### case of eqs (45), (46)
        #
        # # update W
        # W .= W .* (((V .+ ϵ) ./ (WH .+ ϵ)) .^ α * transpose(H)) .^ (1/α)

        # #### normalize to unit sums the columns of W if specified
        # #### note: order not as in (45-46)
        # cs_W = sum(W, dims=1)
        # W .= W .* repeat(cs_W .^ -1, size(W, 1), 1)

        # WH .= W * H

        # # update H
        # H .= H .* ((transpose(W) * ((V .+ ϵ)./(WH .+ ϵ)) .^ α)) .^ (1/α)
        # WH .= W * H

        #### eqs (44)

        # update W
        W .= W .* ((((V .+ ϵ) ./ (WH .+ ϵ)) .^ α * transpose(H)) ./ (ones(size(V)) * transpose(H)) .+ ϵ) .^ (1/α)
        WH .= W * H

        # update H
        H .= H .* (((transpose(W) * ((V .+ ϵ) ./ (WH .+ ϵ)) .^ α)) ./ (transpose(W) * ones(size(V))) .+ ϵ) .^ (1/α)
        WH .= W * H

    end
end


text:
"""
Riesz sense fractional derivative
"""
abstract type Riesz <: FracDiffAlg end

"""
# Riesz sense symmetric fractional derivative algorithm.

    fracdiff(f, α, end_point, h, RieszSymmetric())

Compute fractional derivative of Riesz sense using Triangular Strip Matrix algorithm.

### Example

```julia-repl
julia> fracdiff(x->x, 0.5, 1, 0.01, RieszSymmetric())
```
"""
struct RieszSymmetric <: Riesz end

"""
# Riesz sense Ortigueira definition fractional derivative

    fracdiff(f, α, end_point, h, RieszOrtigueira())

Ortigueira's definition of the symmetric Riesz derivative via centred differences

### Usage

```julia-repl
julia> fracdiff(x->x, 0.5, 1, 0.01, RieszOrtigueira())
```
"""
struct RieszOrtigueira <: Riesz end

################################################################
###                    Type definition done                  ###
################################################################

function fracdiff(f::FunctionAndNumber, α, end_point, h, ::RieszSymmetric)
    N=floor(Int, end_point/h)

    mat = RieszMatrix(α, N+1, h)
    return mat*f.(collect(0:h:end_point))
end
function RieszMatrix(α, N, h)
    caputo = B(N+1, α)
    caputo = caputo[2:(N+1), 1:N]
    result = 1/2*(caputo+caputo')
    result = h^(-α)*result
    return result
end

#FIXME: Values are changing when h become smaller?
function fracdiff(f, α, point, h, ::RieszOrtigueira)
    N = round(Int, point/h)
    t = collect(0:h:point)
    return ranort(α, N+1, h)*f.(t)
end

function ranort(alpha, N)
    k=collect(0:N-1)
    rc = ((-1)*ones(size(k))).^k.*gamma.(alpha+1).*(gamma.(alpha*0.5 .-k.+1).*gamma.(alpha*0.5 .+ k.+1)).^(-1)
    rc = rc*(cos(alpha*π*0.5))
    R = zeros(N, N)

    for m=1:N
        R[m, m:N] = rc[1:N-m+1]
    end

    for i=1:N-1
        for j=i:N
            R[j, i] = R[i, j]
        end
    end
    return R
end
# Dispatch for ranort
function ranort(alpha, N, h)
    R = ranort(alpha, N)
    R = R*h^(-alpha)
    return R
end


text:
module HighFrequencyCovariance

using DataFrames
using Dates
using Distributions
using LinearAlgebra
using Random
using Statistics
using StochasticIntegrals # Used for Monte Carlo Functions
using UnivariateFunctions # Used for Monte Carlo Functions

# Preliminary functions
include("00_Structs.jl")
export SortedDataFrame, CovarianceMatrix, +, -
export get_correlation, get_volatility
export make_nan_covariance_matrix, duration
export subset_to_tick, subset_to_time, calculate_mean_abs_distance
export to_dataframe, valid_correlation_matrix, is_psd_matrix
export ticks_per_asset, get_assets
export show
include("01_helpers.jl")
export simple_differencing
include("02_subsample_ticks.jl")
export cov2cor, cor2cov, cov2cor_and_vol, covariance, construct_matrix_from_eigen, get_returns
export combine_covariance_matrices, rearrange
export squared_frobenius, squared_frobenius_distance
export get_all_refresh_times, latest_value, time_between_refreshes, random_value_in_interval
include("03_MonteCarlo.jl")
export generate_random_path, ItoSet
include("04_Serialisation.jl")
export to_dataframe, dataframe_to_covariancematrix

# Volatility estimation techniques
include("10_volatility_simple.jl")
export simple_volatility
include("11_volatility_two_scales.jl")
export two_scales_volatility, default_num_grids

# Covariance Estimation techniques
include("20_covariance_simple.jl")
export simple_covariance, simple_covariance_given_returns, default_spacing
include("21_covariance_bnhls.jl")
export HFC_Kernel, parzen, quadratic_spectral, fejer, tukey_hanning, bnhls_2008
export bnhls_covariance, preaveraging_end_returns
include("22_covariance_spectral.jl")
export spectral_covariance
include("23_covariance_preaveragingHY.jl")
export preaveraged_covariance
include("24_covariance_two_scales.jl")
export two_scales_covariance

# Regularisation Techniques
include("30_Regularisation_identity.jl")
export identity_regularisation
include("31_Regularisation_NearestCorrelationMatrix.jl")
export project_to_U, project_to_S, iterate_higham, nearest_correlation_matrix, nearest_psd_matrix
include("32_Regularisation_RMT.jl")
export eigenvalue_clean

# Blocking and regularisation.
include("40_Blocking.jl")
export blockwise_estimation, put_assets_into_blocks_by_trading_frequency

# Convenience Wrappers
include("50_ConvenienceFunctions.jl")
export estimate_volatility, estimate_covariance, regularise
export estimate_microstructure_noise

end