text:
import Temporal.acf  # used for running autocorrelation function

using Statistics

"""
```
runmean(x::Array{T}; n::Int=10, cumulative::Bool=true)::Array{T}
runmean(X::Matrix; n::Int=10, cumulative::Bool=true)::Matrix
```

Compute a running or rolling arithmetic mean of an array.
"""
function runmean(x::Vector{T}; n::Int=10, cumulative::Bool=true)::Vector{T} where {T<:AbstractFloat}
    @assert n<size(x,1) && n>1 "Argument n is out of bounds."
    out = zeros(size(x,1))
    out[1:n-1] .= NaN
    if cumulative
        fi = 1.0:size(x,1)
        @inbounds for i in n:size(x,1)
            out[i] = sum(x[1:i])/fi[i]
        end
    else
        # original shorter but slower version of the code below
        # @inbounds for i = n:size(x,1)
        #     out[i] = mean(x[i-n+1:i])
        # end
        n_current = 0
        s::T = 0
        @inbounds for i = 1:n
            if isfinite(x[i])
                s += x[i]
                n_current += 1
            end
        end
        out[n] = n_current == n ? s / n : NaN
        @inbounds for i = n+1:size(x,1)
            if isfinite(x[i])
                s += x[i]
                n_current += 1
            end
            if isfinite(x[i-n])
                s -= x[i-n]
                n_current -= 1
            end
            out[i] = (n_current == n) ? s / n_current : NaN
        end
    end
    return out
end
runmean(X::Matrix; n::Int=10, cumulative::Bool=true)::Matrix = hcat((runmean(X[:,j], n=n, cumulative=cumulative) for j in 1:size(X,2))...)

"""
```
runsum(x::Vector{T}; n::Int=10, cumulative::Bool=true)::Vector{T}
runsum(X::Matrix; n::Int=10, cumulative::Bool=true)::Matrix
```

Compute a running or rolling summation of an array.
"""
function runsum(x::Vector{T}; n::Int=10, cumulative::Bool=true)::Vector{T} where {T<:AbstractFloat}
    @assert n<size(x,1) && n>1 "Argument n is out of bounds."
    if cumulative
        out = cumsum(x, dims=1)
        out[1:n-1] .= NaN
    else
        out = zeros(size(x))
        out[1:n-1] .= NaN
        # original shorter but slower version of the code below
        # @inbounds for i = n:size(x,1)
        #     out[i] = sum(x[i-n+1:i])
        # end
        n_current = 0
        s::T = 0
        @inbounds for i = 1:n
            if isfinite(x[i])
                s += x[i]
                n_current += 1
            end
        end
        out[n] = n_current == n ? s : NaN
        @inbounds for i = n+1:size(x,1)
            if isfinite(x[i])
                s += x[i]
                n_current += 1
            end
            if isfinite(x[i-n])
                s -= x[i-n]
                n_current -= 1
            end
            out[i] = (n_current == n) ? s : NaN
        end
    end
    return out
end
runsum(X::Matrix; n::Int=10, cumulative::Bool=true)::Matrix = hcat((runsum(X[:,j], n=n, cumulative=cumulative) for j in 1:size(X,2))...)

"""
```
runmad(x::Vector{T}; n::Int=10, cumulative::Bool=true, fun::Function=median)::Vector{T}
runmad(X::Matrix; n::Int=10, cumulative::Bool=true)::Matrix
```

Compute the running or rolling mean absolute deviation of an array
"""
function runmad(x::Vector{T}; n::Int=10, cumulative::Bool=true, fun::Function=median)::Vector{T} where {T<:AbstractFloat}
    @assert n<size(x,1) && n>1 "Argument n is out of bounds."
    out = zeros(size(x))
    out[1:n-1] .= NaN
    center = 0.0
    if cumulative
        fi = collect(1.0:size(x,1))
        @inbounds for i = n:size(x,1)
            center = fun(x[1:i])
            out[i] = sum(abs.(x[1:i].-center)) / fi[i]
        end
    else
        fn = float(n)
        @inbounds for i = n:size(x,1)
            center = fun(x[i-n+1:i])
            out[i] = sum(abs.(x[i-n+1:i].-center)) / fn
        end
    end
    return out
end
runmad(X::Matrix; n::Int=10, cumulative::Bool=true)::Matrix = hcat((runmad(X[:,j], n=n, cumulative=cumulative) for j in 1:size(X,2))...)

"""
```
runvar(x::Vector{T}; n::Int=10, cumulative=true)::Vector{T}
runvar(X::Matrix; n::Int=10, cumulative::Bool=true)::Matrix
```

Compute the running or rolling variance of an array
"""
function runvar(x::Vector{T}; n::Int=10, cumulative=true)::Vector{T} where {T<:AbstractFloat}
    @assert n<size(x,1) && n>1 "Argument n is out of bounds."
    out = zeros(size(x))
    out[1:n-1] .= NaN
    if cumulative
        @inbounds for i = n:size(x,1)
            out[i] = var(x[1:i])
        end
    else
        @inbounds for i = n:size(x,1)
            out[i] = var(x[i-n+1:i])
        end
    end
    return out
end
runvar(X::Matrix; n::Int=10, cumulative::Bool=true)::Matrix = hcat((runvar(X[:,j], n=n, cumulative=cumulative) for j in 1:size(X,2))...)

"""
```
runsd(x::Vector{T}; n::Int=10, cumulative::Bool=true)::Vector{T}
runsd(X::Matrix; n::Int=10, cumulative::Bool=true)::Matrix
```

Compute the running or rolling standard deviation of an array
"""
function runsd(x::Vector{T}; n::Int=10, cumulative::Bool=true)::Vector{T} where {T<:AbstractFloat}
    return sqrt.(runvar(x, n=n, cumulative=cumulative))
end
runsd(X::Matrix; n::Int=10, cumulative::Bool=true)::Matrix = sqrt.(hcat((runvar(X[:,j], n=n, cumulative=cumulative) for j in 1:size(X,2))...))

"""
```
runcov(x::Vector{T}, y::Vector{T}; n::Int=10, cumulative::Bool=true)::Vector{T}
runcov(X::Matrix, Y::Matrix; n::Int=10, cumulative::Bool=true)::Matrix
```

Compute the running or rolling covariance of two arrays
"""
function runcov(x::Vector{T}, y::Vector{T}; n::Int=10, cumulative::Bool=true)::Vector{T} where {T<:AbstractFloat}
    @assert length(x) == length(y) "Dimension mismatch: length of `x` not equal to length of `y`."
    @assert n<size(x,1) && n>1 "Argument n is out of bounds."
    out = zeros(size(x))
    out[1:n-1] .= NaN
    if cumulative
        @inbounds for i = n:length(x)
            out[i] = cov(x[1:i], y[1:i])
        end
    else
        @inbounds for i = n:length(x)
            out[i] = cov(x[i-n+1:i], y[i-n+1:i])
        end
    end
    return out
end
runcov(X::Matrix, Y::Matrix; n::Int=10, cumulative::Bool=true)::Matrix = hcat((runcov(X[:,j], Y[:,j], n=n, cumulative=cumulative) for j in 1:size(X,2))...)

"""
```
runcor(x::Vector{T}, y::Vector{T}; n::Int=10, cumulative::Bool=true)::Vector{T}
runcor(X::Matrix, y::Vector; n::Int=10, cumulative::Bool=true)::Matrix
runcor(X::Matrix, Y::Matrix; n::Int=10, cumulative::Bool=true)::Matrix
```

Compute the running or rolling correlation of two arrays
"""
function runcor(x::Vector{T}, y::Vector{T}; n::Int=10, cumulative::Bool=true)::Vector{T} where {T<:Float64}
    @assert length(x) == length(y) "Dimension mismatch: length of `x` not equal to length of `y`."
    @assert n<size(x,1) && n>1 "Argument n is out of bounds."
    out = zeros(size(x))
    out[1:n-1] .= NaN
    if cumulative
        @inbounds for i = n:length(x)
            out[i] = cor(x[1:i], y[1:i])
        end
    else
        @inbounds for i = n:length(x)
            out[i] = cor(x[i-n+1:i], y[i-n+1:i])
        end
    end
    return out
end
runcor(X::Matrix, Y::Matrix; n::Int=10, cumulative::Bool=true)::Matrix = hcat((runcor(X[:,j], Y[:,j], n=n, cumulative=cumulative) for j in 1:size(X,2))...)
runcor(X::Matrix, y::Vector; n::Int=10, cumulative::Bool=true)::Matrix = hcat((runcor(X[:,j], y, n=n, cumulative=cumulative) for j in 1:size(X,2))...)

"""
```
runmax(x::Vector{T}; n::Int=10, cumulative::Bool=true, inclusive::Bool=true)::Vector{T}
runmax(X::Matrix; n::Int=10, cumulative::Bool=true, inclusive::Bool=true)::Matrix
```

Compute the running or rolling maximum of an array
"""
function runmax(x::Vector{T}; n::Int=10, cumulative::Bool=true, inclusive::Bool=true)::Vector{T} where {T<:AbstractFloat}
    @assert n<size(x,1) && n>1 "Argument n is out of bounds."
    out = zeros(size(x))
    if inclusive
        if cumulative
            out[n] = maximum(x[1:n])
            @inbounds for i = n+1:size(x,1)
                out[i] = max(out[i-1], x[i])
            end
        else
            @inbounds for i = n:size(x,1)
                out[i] = maximum(x[i-n+1:i])
            end
        end
        out[1:n-1] .= NaN
        return out
    else
        if cumulative
            out[n+1] = maximum(x[1:n])
            @inbounds for i = n+1:size(x,1)-1
                out[i+1] = max(out[i-1], x[i-1])
            end
        else
            @inbounds for i = n:size(x,1)-1
                out[i+1] = maximum(x[i-n+1:i])
            end
        end
        out[1:n] .= NaN
        return out
    end
end
runmax(X::Matrix; n::Int=10, cumulative::Bool=true, inclusive::Bool=true)::Matrix = hcat((runmax(X[:,j], n=n, cumulative=cumulative, inclusive=inclusive) for j in 1:size(X,2))...)

"""
```
runmin(x::Vector{T}; n::Int=10, cumulative::Bool=true, inclusive::Bool=true)::Vector{T}
runmin(X::Matrix; n::Int=10, cumulative::Bool=true, inclusive::Bool=true)::Matrix
```

Compute the running or rolling minimum of an array
"""
function runmin(x::Vector{T}; n::Int=10, cumulative::Bool=true, inclusive::Bool=true)::Vector{T} where {T<:AbstractFloat}
    @assert n<size(x,1) && n>1 "Argument n is out of bounds."
    out = zeros(size(x))
    if inclusive
        if cumulative
            out[n] = minimum(x[1:n])
            @inbounds for i = n+1:size(x,1)
                out[i] = min(out[i-1], x[i])
            end
        else
            @inbounds for i = n:size(x,1)
                out[i] = minimum(x[i-n+1:i])
            end
        end
        out[1:n-1] .= NaN
        return out
    else
        if cumulative
            out[n+1] = minimum(x[1:n])
            @inbounds for i = n+1:size(x,1)-1
                out[i+1] = min(out[i-1], x[i-1])
            end
        else
            @inbounds for i = n:size(x,1)-1
                out[i+1] = minimum(x[i-n+1:i])
            end
        end
        out[1:n] .= NaN
        return out
    end
end
runmin(X::Matrix; n::Int=10, cumulative::Bool=true, inclusive::Bool=true)::Matrix = hcat((runmin(X[:,j], n=n, cumulative=cumulative, inclusive=inclusive) for j in 1:size(X,2))...)

"""
```
runquantile(x::Vector{T}; p::T=0.05, n::Int=10, cumulative::Bool=true)::Vector{T}
runquantile(X::Matrix; n::Int=10, cumulative::Bool=true, p::Real=0.05)::Matrix
```

Compute the running/rolling quantile of an array
"""
function runquantile(x::Vector{T}; p::T=0.05, n::Int=10, cumulative::Bool=true)::Vector{T} where {T<:AbstractFloat}
    @assert n<size(x,1) && n>1 "Argument n is out of bounds."
    out = zeros(T, size(x,1))
    if cumulative
        @inbounds for i in 2:size(x,1)
            out[i] = quantile(x[1:i], p)
        end
        out[1] = NaN
    else
        @inbounds for i in n:size(x,1)
            out[i] = quantile(x[i-n+1:i], p)
        end
        out[1:n-1] .= NaN
    end
    return out
end
runquantile(X::Matrix; n::Int=10, cumulative::Bool=true, p::Real=0.05)::Matrix = hcat((runquantile(X[:,j], n=n, cumulative=cumulative, p=p) for j in 1:size(X,2))...)

"""
```
function runacf(x::Vector{T};
                n::Int = 10,
                maxlag::Int = n-3,
                lags::AbstractVector{Int,1} = 0:maxlag,
                cumulative::Bool = true)::Matrix{T} where {T<:AbstractFloat}
runacf(X::Matrix; n::Int=10, cumulative::Bool=true, maxlag::Int=n-3, lags::AbstractVector{Int}=0:maxlag)::Matrix
```

Compute the running/rolling autocorrelation of a vector.
"""
function runacf(x::Vector{T};
                n::Int = 10,
                maxlag::Int = n-3,
                lags::AbstractVector{Int} = 0:maxlag,
                cumulative::Bool = true)::Matrix{T} where {T<:AbstractFloat}
    @assert size(x, 2) == 1 "Autocorrelation input array must be one-dimensional"
    N = size(x, 1)
    @assert n < N && n > 0
    if length(lags) == 1 && lags[1] == 0
        return ones((N, 1))
    end
    out = zeros((N, length(lags))) * NaN
    if cumulative
        @inbounds for i in n:N
            out[i,:] = acf(x[1:i], lags=lags)
        end
    else
        @inbounds for i in n:N
            out[i,:] = acf(x[i-n+1:i], lags=lags)
        end
    end
    return out
end
runacf(X::Matrix; n::Int=10, cumulative::Bool=true, maxlag::Int=n-3, lags::AbstractVector{Int}=0:maxlag)::Matrix = hcat((runacf(X[:,j], n=n, cumulative=cumulative, maxlag=maxlag, lags=lags) for j in 1:size(X,2))...)

"""
```
runfun(x::Vector{T}, f::Function; n::Int = 10, cumulative::Bool=false, args...)::Vector{T} where {T<:AbstractFloat}
runfun(X::Matrix{T}, f::Function; n::Int=10, cumulative::Bool=false, args...)::Matrix{T} where {T<:AbstractFloat}
```

Apply a general function `f` that returns a scalar over an array
"""
function runfun(x::Vector{T}, f::Function; n::Int = 10, cumulative::Bool=false, args...)::Vector{T} where {T<:AbstractFloat}
    N = size(x,1)
    out = zeros(T, N) .* NaN
    if cumulative
        for i in n:N
            result::T = f(x[1:i]; args...)
            out[i] = result
        end
    else
        for i in n:N
            result::T = f(x[i-n+1:i]; args...)
            out[i] = result
        end
    end
    return out
end

function runfun(X::Matrix{T}, f::Function; n::Int=10, cumulative::Bool=false, args...)::Matrix{T} where {T<:AbstractFloat}
    return hcat((runfun(X[:,j], f, n=n, cumulative=cumulative; args...) for j in 1:size(X,2))...)
end


text:
using BPPLib
using Test

@testset "BPPLoadSymbol" begin
    data = loadBPP(:PM_u010_01)
    @test length(data.weights) == 10
    @test data.capacity == 150
    @test data.lb == 0
    @test data.ub == typemax(Int64)
    @test data.name == "PM_u010_01"
end

@testset "BPPErrorSymbol" begin
    data = loadBPP(:PM_u010_0)
    @test data === nothing
end

@testset "BPPLoadString" begin
    data = loadBPP(joinpath(pkgdir(BPPLib), "test/data/bpp/Falkenauer_u120_00.txt"))
    @test length(data.weights) == 120
    @test data.capacity == 150
    @test data.lb == 48
    @test data.ub == 48
    @test data.name == "Falkenauer_u120_00"
end

@testset "BPPErrorString" begin
    data = loadBPP(joinpath(pkgdir(BPPLib), "not a file"))
    @test data === nothing
end

@testset "BPPShow" begin
    data = loadBPP(:PM_u010_01)
    @test_nowarn println(data)
    data = loadBPP(joinpath(pkgdir(BPPLib), "test/data/bpp/Falkenauer_u120_00.txt"))
    @test_nowarn println(data)
end

@testset "CSPLoadSymbol" begin
    data = loadCSP(:PM_u010_01)
    @test length(data.weights) == 10
    @test length(data.demands) == 10
    @test data.capacity == 150
    @test data.lb == 0
    @test data.ub == typemax(Int64)
    @test data.name == "PM_u010_01"
end

@testset "CSPErrorSymbol" begin
    data = loadCSP(:PM_u010_0)
    @test data === nothing
end

@testset "CSPLoadString" begin
    data = loadCSP(joinpath(pkgdir(BPPLib), "test/data/csp/Falkenauer_u120_00.txt"))
    @test length(data.weights) == 58
    @test length(data.demands) == 58
    @test data.capacity == 150
    @test data.lb == 48
    @test data.ub == 48
    @test data.name == "Falkenauer_u120_00"
end

@testset "CSPErrorString" begin
    data = loadCSP(joinpath(pkgdir(BPPLib), "not a file"))
    @test data === nothing
end

@testset "CSPShow" begin
    data = loadCSP(:PM_u010_01)
    @test_nowarn println(data)
    data = loadCSP(joinpath(pkgdir(BPPLib), "test/data/csp/Falkenauer_u120_00.txt"))
    @test_nowarn println(data)
end

@testset "toBPP" begin
    data = loadCSP(:PM_u010_01)
    data2 = toBPP(data)
    @test length(data2.weights) == 10
    @test data2.capacity == 150
    @test data2.name == "PM_u010_01"
end

@testset "toCSP" begin
    data = loadBPP(:PM_u010_01)
    data2 = toCSP(data)
    @test length(data2.weights) == 10
    @test length(data2.demands) == 10
    @test data2.capacity == 150
    @test data2.name == "PM_u010_01"
end

ext: jl


text:
struct KroneckerMap{T, As<:LinearMapTuple} <: LinearMap{T}
    maps::As
    function KroneckerMap{T}(maps::LinearMapTuple) where {T}
        for TA in Base.Generator(eltype, maps)
            promote_type(T, TA) == T ||
                error("eltype $TA cannot be promoted to $T in KroneckerMap constructor")
        end
        return new{T,typeof(maps)}(maps)
    end
end

"""
    kron(A::LinearMap, B::LinearMap)::KroneckerMap
    kron(A, B, Cs...)::KroneckerMap

Construct a (lazy) representation of the Kronecker product `A⊗B`. One of the two factors
can be an `AbstractMatrix`, which is then promoted to a `LinearMap` automatically.

To avoid fallback to the generic [`Base.kron`](@ref) in the multi-map case,
there must be a `LinearMap` object among the first 8 arguments in usage like
`kron(A, B, Cs...)`.

For convenience, one can also use `A ⊗ B` or `⊗(A, B, Cs...)` (typed as `\\otimes+TAB`) to construct the
`KroneckerMap`, even when all arguments are of `AbstractMatrix` type.

If `A`, `B`, `C` and `D` are linear maps of such size that one can form the matrix
products `A*C` and `B*D`, then the mixed-product property `(A⊗B)*(C⊗D) = (A*C)⊗(B*D)`
holds. Upon vector multiplication, this rule is checked for applicability.

# Examples
```jldoctest; setup=(using LinearAlgebra, SparseArrays, LinearMaps)
julia> J = LinearMap(I, 2) # 2×2 identity map
2×2 LinearMaps.UniformScalingMap{Bool} with scaling factor: true

julia> E = spdiagm(-1 => trues(1)); D = E + E' - 2I;

julia> Δ = kron(D, J) + kron(J, D); # discrete 2D-Laplace operator

julia> Matrix(Δ)
4×4 Array{Int64,2}:
 -4   1   1   0
  1  -4   0   1
  1   0  -4   1
  0   1   1  -4
```
"""
Base.kron(A::LinearMap, B::LinearMap) =
    KroneckerMap{promote_type(eltype(A), eltype(B))}((A, B))
Base.kron(A::LinearMap, B::KroneckerMap) =
    KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A, B.maps...))
Base.kron(A::KroneckerMap, B::LinearMap) =
    KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A.maps..., B))
Base.kron(A::KroneckerMap, B::KroneckerMap) =
    KroneckerMap{promote_type(eltype(A), eltype(B))}(tuple(A.maps..., B.maps...))
# hoist out scalings
Base.kron(A::ScaledMap, B::LinearMap) = A.λ * kron(A.lmap, B)
Base.kron(A::LinearMap, B::ScaledMap) = kron(A, B.lmap) * B.λ
Base.kron(A::ScaledMap, B::ScaledMap) = (A.λ * B.λ) * kron(A.lmap, B.lmap)
# reduce UniformScalingMaps
Base.kron(A::UniformScalingMap, B::UniformScalingMap) = UniformScalingMap(A.λ * B.λ, A.M * B.M)
# disambiguation
Base.kron(A::ScaledMap, B::KroneckerMap) = A.λ * kron(A.lmap, B)
Base.kron(A::KroneckerMap, B::ScaledMap) = kron(A, B.lmap) * B.λ
# generic definitions
Base.kron(A::LinearMap, B::LinearMap, C::LinearMap, Ds::LinearMap...) =
    kron(kron(A, B), C, Ds...)
Base.kron(A::AbstractVecOrMat, B::LinearMap) = kron(LinearMap(A), B)
Base.kron(A::LinearMap, B::AbstractVecOrMat) = kron(A, LinearMap(B))
# promote AbstractMatrix arguments to LinearMaps, then take LinearMap-Kronecker product
for k in 3:8 # is 8 sufficient?
    Is = ntuple(n->:($(Symbol(:A, n))::AbstractVecOrMat), Val(k-1))
    # yields (:A1, :A2, :A3, ..., :A(k-1))
    L = :($(Symbol(:A, k))::LinearMap)
    # yields :Ak::LinearMap
    mapargs = ntuple(n -> :(LinearMap($(Symbol(:A, n)))), Val(k-1))
    # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1)))

    @eval Base.kron($(Is...), $L, As::MapOrVecOrMat...) =
        kron($(mapargs...), $(Symbol(:A, k)), convert_to_lmaps(As...)...)
end

struct KronPower{p}
    function KronPower(p::Integer)
        p > 1 || throw(ArgumentError("the Kronecker power is only defined for exponents larger than 1, got $k"))
        return new{p}()
    end
end

"""
    ⊗(k::Integer)

Construct a lazy representation of the `k`-th Kronecker power
`A^⊗(k) = A ⊗ A ⊗ ... ⊗ A`, where `A` can be an `AbstractMatrix` or a `LinearMap`.
"""
⊗(k::Integer) = KronPower(k)

⊗(A, B, Cs...) = kron(convert_to_lmaps(A, B, Cs...)...)

Base.:(^)(A::MapOrMatrix, ::KronPower{p}) where {p} =
    kron(ntuple(n -> convert_to_lmaps_(A), Val(p))...)

Base.size(A::KroneckerMap) = map(*, size.(A.maps)...)

LinearAlgebra.issymmetric(A::KroneckerMap) = all(issymmetric, A.maps)
LinearAlgebra.ishermitian(A::KroneckerMap) = all(ishermitian, A.maps)

LinearAlgebra.adjoint(A::KroneckerMap) = KroneckerMap{eltype(A)}(map(adjoint, A.maps))
LinearAlgebra.transpose(A::KroneckerMap) = KroneckerMap{eltype(A)}(map(transpose, A.maps))

Base.:(==)(A::KroneckerMap, B::KroneckerMap) = (eltype(A) == eltype(B) && A.maps == B.maps)

#################
# multiplication helper functions
#################

@inline function _kronmul!(y, B, x, A, T)
    ma, na = size(A)
    mb, nb = size(B)
    X = reshape(x, (nb, na))
    Y = reshape(y, (mb, ma))
    if B isa UniformScalingMap
        _unsafe_mul!(Y, X, transpose(A))
        lmul!(B.λ, y)
    else
        temp = similar(Y, (ma, nb))
        _unsafe_mul!(temp, A, copy(transpose(X)))
        _unsafe_mul!(Y, B, transpose(temp))
    end
    return y
end
@inline function _kronmul!(y, B, x, A::UniformScalingMap, _)
    ma, na = size(A)
    mb, nb = size(B)
    iszero(A.λ) && return fill!(y, zero(eltype(y)))
    X = reshape(x, (nb, na))
    Y = reshape(y, (mb, ma))
    _unsafe_mul!(Y, B, X)
    !isone(A.λ) && rmul!(y, A.λ)
    return y
end
@inline function _kronmul!(y, B, x, A::VecOrMatMap, _)
    ma, na = size(A)
    mb, nb = size(B)
    X = reshape(x, (nb, na))
    Y = reshape(y, (mb, ma))
    At = transpose(A.lmap)
    if B isa UniformScalingMap
        # the following is (perhaps due to the reshape?) faster than
        # _unsafe_mul!(Y, B * X, At)
        _unsafe_mul!(Y, X, At)
        lmul!(B.λ, y)
    elseif nb*ma <= mb*na
        _unsafe_mul!(Y, B, X * At)
    else
        _unsafe_mul!(Y, Matrix(B * X), At)
    end
    return y
end
const VectorMap{T} = WrappedMap{T,<:AbstractVector}
const AdjOrTransVectorMap{T} = WrappedMap{T,<:LinearAlgebra.AdjOrTransAbsVec}
@inline _kronmul!(y, B::AdjOrTransVectorMap, x, a::VectorMap, _) = mul!(y, a.lmap, B.lmap * x)

#################
# multiplication with vectors
#################

const KroneckerMap2{T} = KroneckerMap{T, <:Tuple{LinearMap, LinearMap}}

function _unsafe_mul!(y::AbstractVecOrMat, L::KroneckerMap2, x::AbstractVector)
    require_one_based_indexing(y)
    A, B = L.maps
    _kronmul!(y, B, x, A, eltype(L))
    return y
end
function _unsafe_mul!(y::AbstractVecOrMat, L::KroneckerMap, x::AbstractVector)
    require_one_based_indexing(y)
    A = first(L.maps)
    B = kron(Base.tail(L.maps)...)
    _kronmul!(y, B, x, A, eltype(L))
    return y
end
# mixed-product rule, prefer the right if possible
# (A₁ ⊗ A₂ ⊗ ... ⊗ Aᵣ) * (B₁ ⊗ B₂ ⊗ ... ⊗ Bᵣ) = (A₁B₁) ⊗ (A₂B₂) ⊗ ... ⊗ (AᵣBᵣ)
function _unsafe_mul!(y::AbstractVecOrMat,
                        L::CompositeMap{<:Any,<:Tuple{KroneckerMap,KroneckerMap}},
                        x::AbstractVector)
    require_one_based_indexing(y)
    B, A = L.maps
    if length(A.maps) == length(B.maps) && all(_iscompatible, zip(A.maps, B.maps))
        _unsafe_mul!(y, kron(map(*, A.maps, B.maps)...), x)
    else
        _unsafe_mul!(y, LinearMap(A)*B, x)
    end
    return y
end
# mixed-product rule, prefer the right if possible
# (A₁⊗B₁) * (A₂⊗B₂) * ... * (Aᵣ⊗Bᵣ) = (A₁*A₂*...*Aᵣ) ⊗ (B₁*B₂*...*Bᵣ)
function _unsafe_mul!(y::AbstractVecOrMat,
                        L::CompositeMap{T, <:Tuple{Vararg{KroneckerMap2}}},
                        x::AbstractVector) where {T}
    require_one_based_indexing(y)
    As = map(AB -> AB.maps[1], L.maps)
    Bs = map(AB -> AB.maps[2], L.maps)
    As1, As2 = Base.front(As), Base.tail(As)
    Bs1, Bs2 = Base.front(Bs), Base.tail(Bs)
    apply = all(_iscompatible, zip(As1, As2)) && all(_iscompatible, zip(Bs1, Bs2))
    if apply
        _unsafe_mul!(y, kron(prod(As), prod(Bs)), x)
    else
        _unsafe_mul!(y, CompositeMap{T}(map(LinearMap, L.maps)), x)
    end
    return y
end

###############
# KroneckerSumMap
###############
struct KroneckerSumMap{T, As<:Tuple{LinearMap, LinearMap}} <: LinearMap{T}
    maps::As
    function KroneckerSumMap{T}(maps::Tuple{LinearMap,LinearMap}) where {T}
        A1, A2 = maps
        (size(A1, 1) == size(A1, 2) && size(A2, 1) == size(A2, 2)) ||
            throw(ArgumentError("operators need to be square in Kronecker sums"))
        for TA in Base.Generator(eltype, maps)
            promote_type(T, TA) == T ||
                error("eltype $TA cannot be promoted to $T in KroneckerSumMap constructor")
        end
        return new{T, typeof(maps)}(maps)
    end
end

"""
    kronsum(A, B)::KroneckerSumMap
    kronsum(A, B, Cs...)::KroneckerSumMap

Construct a (lazy) representation of the Kronecker sum `A⊕B = A ⊗ Ib + Ia ⊗ B`
of two square linear maps of type `LinearMap` or `AbstractMatrix`. Here, `Ia`
and `Ib` are identity operators of the size of `A` and `B`, respectively.
Arguments of type `AbstractMatrix` are automatically promoted to `LinearMap`.

For convenience, one can also use `A ⊕ B` or `⊕(A, B, Cs...)` (typed as
`\\oplus+TAB`) to construct the `KroneckerSumMap`.

# Examples
```jldoctest; setup=(using LinearAlgebra, SparseArrays, LinearMaps)
julia> J = LinearMap(I, 2) # 2×2 identity map
2×2 LinearMaps.UniformScalingMap{Bool} with scaling factor: true

julia> E = spdiagm(-1 => trues(1)); D = LinearMap(E + E' - 2I);

julia> Δ₁ = kron(D, J) + kron(J, D); # discrete 2D-Laplace operator, Kronecker sum

julia> Δ₂ = kronsum(D, D);

julia> Δ₃ = D^⊕(2);

julia> Matrix(Δ₁) == Matrix(Δ₂) == Matrix(Δ₃)
true
```
"""
kronsum(A::MapOrMatrix, B::MapOrMatrix) =
    KroneckerSumMap{promote_type(eltype(A), eltype(B))}(convert_to_lmaps(A, B))
kronsum(A::MapOrMatrix, B::MapOrMatrix, C::MapOrMatrix, Ds::MapOrMatrix...) =
    kronsum(A, kronsum(B, C, Ds...))

struct KronSumPower{p}
    function KronSumPower(p::Integer)
        p > 1 || throw(ArgumentError("the Kronecker sum power is only defined for exponents larger than 1, got $k"))
        return new{p}()
    end
end

"""
    ⊕(k::Integer)

Construct a lazy representation of the `k`-th Kronecker sum power `A^⊕(k) = A ⊕ A ⊕ ... ⊕ A`,
where `A` can be a square `AbstractMatrix` or a `LinearMap`.
"""
⊕(k::Integer) = KronSumPower(k)

⊕(a, b, c...) = kronsum(a, b, c...)

Base.:(^)(A::MapOrMatrix, ::KronSumPower{p}) where {p} =
    kronsum(ntuple(n->convert_to_lmaps_(A), Val(p))...)

Base.size(A::KroneckerSumMap, i) = prod(size.(A.maps, i))
Base.size(A::KroneckerSumMap) = (size(A, 1), size(A, 2))

LinearAlgebra.issymmetric(A::KroneckerSumMap) = all(issymmetric, A.maps)
LinearAlgebra.ishermitian(A::KroneckerSumMap) = all(ishermitian, A.maps)

LinearAlgebra.adjoint(A::KroneckerSumMap) =
    KroneckerSumMap{eltype(A)}(map(adjoint, A.maps))
LinearAlgebra.transpose(A::KroneckerSumMap) =
    KroneckerSumMap{eltype(A)}(map(transpose, A.maps))

Base.:(==)(A::KroneckerSumMap, B::KroneckerSumMap) =
    (eltype(A) == eltype(B) && A.maps == B.maps)

function _unsafe_mul!(y::AbstractVecOrMat, L::KroneckerSumMap, x::AbstractVector)
    A, B = L.maps
    ma, na = size(A)
    mb, nb = size(B)
    X = reshape(x, (nb, na))
    Y = reshape(y, (nb, na))
    _unsafe_mul!(Y, X, transpose(A))
    _unsafe_mul!(Y, B, X, true, true)
    return y
end


text:
module TestMulticlassLGBM

using MLJBase
using Test
using Random: seed!

import CategoricalArrays
import LightGBM

## CLASSIFIER -- shamelessly adapted from binary_classification.jl

# test multiclass case:
N = 5
Nsamples = 3000
seed!(0)

model = LightGBM.MLJInterface.LGBMClassifier(num_iterations=100)

X       = (x1=rand(Nsamples), x2=rand(Nsamples), x3=rand(Nsamples))
ycat    = string.(mod.(round.(Int, X.x1 * 10), N)) |> MLJBase.categorical
weights = Float64.(MLJBase.int(ycat)) # just use the 1's/2's directly as multipliers

y = MLJBase.identity.(ycat) # make plain Vector with categ. elements (actually not sure what this is doing)

# fit once, without weights
train, test              = MLJBase.partition(MLJBase.eachindex(y), 0.6)
fitresult, cache, report = MLJBase.fit(model, 0, MLJBase.selectrows(X, train), y[train];)
yhat                     = MLJBase.mode.(MLJBase.predict(model, fitresult, MLJBase.selectrows(X, test)))

# fit again with weights
fitresult, cache, report = MLJBase.fit(model, 0, MLJBase.selectrows(X, train), y[train], weights[train])
yhat_with_weights        = MLJBase.mode.(MLJBase.predict(model, fitresult, MLJBase.selectrows(X, test)))

misclassification_rate   = sum(yhat .!= y[test])/length(test)

@test misclassification_rate < 0.015

# All we can really say about fitting with/without weights for this example is that the solutions shouldn't be identical
@test yhat_with_weights != yhat

# Cache contains iterations counts history
@test cache isa NamedTuple
@test cache.num_boostings_done == [100]

@test isa(report, NamedTuple)

expected_return_type = Tuple{
    LightGBM.LGBMClassification,
    CategoricalArrays.CategoricalArray,
    LightGBM.MLJInterface.LGBMClassifier,
}

@test isa(fitresult, expected_return_type)

# Provided by Anthony Blaom as a simple integration test
X, y = @load_iris;
model = LightGBM.MLJInterface.LGBMClassifier()
yhat = fit!(machine(model, X, y); verbosity=0) |> predict;
@test scitype(mode.(yhat)) <: AbstractVector{Multiclass{3}}
@test mean(cross_entropy(yhat, y)) < 0.6  # or do proper out-of-sample test

end # module
true


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module PendulumDAE

DAE-model of a mass point attached via a rod to a revolute joint (GGL formulation).

Starting equations (der(xx) = dxx/dt):

     der(x) = vx
     der(y) = vy
  m*der(vx) = -x*lambda
  m*der(vy) = -y*lambda-m*g
  x*x + y*y = L*L

Gear/Gupta/Leimkuhler equations:
   # lambda = der(lambda_int)
   # mue    = der(mue_int)
   # g      = x*x + y*y - L*L
   # G      = dg/d[x;y] = [2*x 2*y]

   der([x,y]) = [vx,vy] - G'*mue
   0 = m*der(vx) + lambda*x
   0 = m*der(vy) + lambda*y + m*g
   0 = x*x + y*y - L*L
   0 = 2*x*vx + 2*y*vy

Arguments of getResidues! function:
   _x    = [x; y; vx; vy; lambda_int; mue_int]
   _derx = [der(x); der(y); der(vx); der(vy); lambda; mue]

Modelica model: ModelicaReferenceModels.ODAEs.Pendulum
"""
module PendulumDAE

import ..ModiaMath

mutable struct Model <: ModiaMath.AbstractSimulationModel
    simulationState::ModiaMath.SimulationState

    # Parameters
    L::Float64
    m::Float64
    g::Float64

    function Model(;L=1.0, m=1.0, g=9.81, x0=L / 2.0, y0=-0.5, x_fixed=false)   #y0=-sqrt(L*L - x0*x0))
        @assert(L > 0.0)
        @assert(m > 0.0)
        @assert(-L <= x0 <= L)
        @assert(-L <= y0 <= L)
        simulationState = ModiaMath.SimulationState("PendulumDAE", getModelResidues!, [x0,y0,1.0,1.0,0.0,0.0], getVariableName;
                                x_fixed=[x_fixed, false, x_fixed, false, false, false],
                                is_constraint = [false,false,false,false,true,true],
								has_constraintDerivatives = true)
        new(simulationState, L, m, g)
    end
end

getVariableName(model, vcat, vindex) = ModiaMath.getVariableName(model, vcat, vindex;
                                                         xNames=["x", "y", "vx", "vy", "lambda_int", "mue_int"],
                                                         derxNames=["der(x)", "der(y)", "der(vx)", "der(vy)", "lambda", "mue"])

function getModelResidues!(m::Model, t::Float64, _x::Vector{Float64}, _derx::Vector{Float64}, _r::Vector{Float64}, _w::Vector{Float64})
    x      = _x[1]
    y      = _x[2]
    vx     = _x[3]
    vy     = _x[4]
    dervx  = _derx[3]
    dervy  = _derx[4]
    lambda = _derx[5]
    mue    = _derx[6]

    _r[1] = _derx[1] - vx + x * mue
    _r[2] = _derx[2] - vy + y * mue
    _r[3] = m.m * dervx + lambda * x
    _r[4] = m.m * dervy + lambda * y + m.m * m.g

	if ModiaMath.compute_der_fc(m)
		_r[5] = x * _derx[1] + y * _derx[2]
		_r[6] = x * dervx + y * dervy + vx*vx + vy*vy
	else
		_r[5] = (x * x + y * y - m.L * m.L) / 2.0
		_r[6] = x * vx + y * vy
    end
    return nothing
end

end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module AndreNumbers
using Nemo

export ModuleAndreNumbers
export André, C000111, V000111, V178963, V178964, V181936, V250283

"""

Generalized André numbers count the ``m``-alternating permutations of length ``n``, cf. A181937.

```
[  SEQ  ] n|k [0][1][2][3][4] [5] [6]  [7]   [8]   [9]  [10]
[V000012] [1]  1, 1, 1, 1, 1,  1,  1,   1,    1,    1,     1
[V000111] [2]  1, 1, 1, 2, 5, 16, 61, 272, 1385, 7936, 50521
[V178963] [3]  1, 1, 1, 1, 3,  9, 19,  99,  477, 1513, 11259
[V178964] [4]  1, 1, 1, 1, 1,  4, 14,  34,   69,  496,  2896
[V181936] [5]  1, 1, 1, 1, 1,  1,  5,  20,   55,  125,   251
[V250283] [6]  1, 1, 1, 1, 1,  1,  1,   6,   27,   83,   209
```

* André, C000111, V000111, V178963, V178964, V181936, V250283.
"""
const ModuleAndreNumbers = ""

"""

Return the generalized André numbers which are the ``m``-alternating permutations of length ``n``, cf. A181937.
"""
function André(m::Int, n::Int)
    haskey(CacheAndré, (m, n)) && return CacheAndré[(m, n)]
    n ≤ 0 && return ZZ(1)
    r = range(0, step = m, stop = n - 1)
    S = sum(binomial(ZZ(n), ZZ(k)) * André(m, k) for k ∈ r)
    return CacheAndré[(m, n)] = n % m == 0 ? -S : S
end
const CacheAndré = Dict{Tuple{Int,Int},fmpz}()

"""

Return the generalized André numbers which are the ``m``-alternating permutations of length ``n``.
"""
V181937(m::Int, n::Int) = abs(André(m, n))

"""

Return the up-down numbers (2-alternating permutations).
"""
V000111(n::Int) = abs(André(2, n))

"""

Return the number of 3-alternating permutations.
"""
V178963(n::Int) = abs(André(3, n))

"""

Return the number of 4-alternating permutations.
"""
V178964(n::Int) = abs(André(4, n))

"""

Return the number of 5-alternating permutations.
"""
V181936(n::Int) = abs(André(5, n))

"""

Return the number of 6-alternating permutations.
"""
V250283(n::Int) = abs(André(6, n))

"""

Generate the André numbers (a.k.a. Euler-up-down numbers A000111). Don't confuse with the Euler numbers A122045.
"""
C000111() =
    Channel(csize = 2) do c
        D = Dict{Int,fmpz}(0 => 1, -1 => 0)
        i = k = 0
        s = 1

        while true
            A = 0
            D[k+s] = 0
            s = -s
            for j ∈ 0:i
                A += D[k]
                D[k] = A
                k += s
            end
            put!(c, A)
            i += 1
        end
    end

#START-TEST-########################################################

using Test, SeqTests, SeqUtils

function test()
    @testset "André" begin

        @test isa(André(2, 10), fmpz)

        @test André(2, 10) == fmpz(-50521)
        @test André(2, 50) == fmpz(-6053285248188621896314383785111649088103498225146815121)
        @test V178963(30) == fmpz(2716778010767155313771539)
        @test V178964(40) == fmpz(11289082167259099068433198467575829)

        if data_installed()
            V = [V000111, V178963, V178964, V181936, V250283]
            for v ∈ V
                SeqTest(v, 'V')
            end
        end

        V = [1, 1, 1, 2, 5, 16, 61, 272, 1385, 7936, 50521]
        generator = C000111()
        for n ∈ 1:10
            v = take!(generator)
            @test V[n] == v
        end

        v = take!(generator)
        close(generator)

        @test isa(v, fmpz)
    end
end

function demo()
    for m ∈ 1:8
        Println([André(m, n) for n ∈ 0:11])
    end

    println()
    generator = C000111()
    for n ∈ 0:10
        v = take!(generator)
        println(n, " ↦ ", v)
    end
    close(generator)
end

"""

for m ∈ 1:20, n ∈ 0:100 André(m,n) end
    0.012252 seconds (66.73 k allocations: 1.341 MiB)
"""
function perf()
    GC.gc()
    @time (for m ∈ 1:20, n ∈ 0:100
        André(m, n)
    end)
end

function main()
    test()
    demo()
    perf()
end

main()

for n ∈ 0:32
    println(André(3, n))
end

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
"""
This function adds the power limits of  hydro generators when there are no CommitmentVariables
"""
function powerconstraints(m::JuMP.AbstractModel, devices::Array{T,1}, time_periods::Int64) where T <: PowerSystems.HydroCurtailment

    phy = m[:phy]
    time_index = m[:phy].axes[2]
    name_index = m[:phy].axes[1]

    (length(phy.axes[2]) != time_periods) ? @error("Length of time dimension inconsistent") : true

    pmax_th = JuMP.JuMPArray(Array{ConstraintRef}(length.(JuMP.axes(phy))), name_index, time_index)
    pmin_th = JuMP.JuMPArray(Array{ConstraintRef}(length.(JuMP.axes(phy))), name_index, time_index)

    for t in phy.axes[2], (ix, name) in enumerate(phy.axes[1])
        if name == devices[ix].name
            pmin_hg[name, t] = @constraint(m, phy[name, t] >= 0.0)
            pmax_hg[name, t] = @constraint(m, phy[name, t] <= devices[ix].tech.activepowerlimits.max * values(devices[ix].scalingfactor)[t])
        else
            @error "Bus name in Array and variable do not match"
        end
    end

    JuMP.register_object(m, :pmax_hydro, pmax_hg)
    JuMP.register_object(m, :pmin_hydro, pmin_hg)

    return m

end


text:
using Test, Pkg

# no need to change pwd for running tests

using AdventOfCode2020:day25
@testset "Day 25" begin
    @testset "Day 25, Unit tests" begin
        @test day25.part1(day25.ex) == 14897079
        @test_broken day25.part2(day25.ex) == nothing
    end

    @testset "Day 24, Puzzle tests (post hoc)" begin
        @test day25.part1(day25.input) == 4968512
        @test_broken day25.part2(day25.input) == nothing
    end
end

using AdventOfCode2020:day24
@testset "Day 24" begin
    @testset "Day 24, Unit tests" begin
        @test day24.part1(day24.ex) == 10
        @test day24.part2(day24.ex) == 2208
    end

    @testset "Day 24, Puzzle tests (post hoc)" begin
        @test day24.part1(day24.input) == 354
        @test day24.part2(day24.input) == 3608
    end
end

using AdventOfCode2020:day23
@testset "Day 23" begin
    @testset "Day 23, Unit tests" begin
        @test day23.part1(day23.ex) == "67384529"
    end

    @testset "Day 23, Puzzle tests (post hoc)" begin
        @test day23.part1(day23.input) == "69852437"
    end
end

using AdventOfCode2020:day22
@testset "Day 22" begin
    @testset "Day 22, Unit tests" begin
        @test day22.part1(day22.ex) == 306
        @test day22.part2(day22.ex) == 291
    end

    @testset "Day 22, Puzzle tests (post hoc)" begin
        @test day22.part1(day22.input) == 33010
        @test day22.part2(day22.input) == 32769
    end
end

using AdventOfCode2020:day21
@testset "Day 21" begin
    @testset "Day 21, Unit tests" begin
        @test day21.part1(day21.ex) == 5
        @test day21.part2(day21.ex) == "mxmxvkd,sqjhc,fvjkl"
    end

    @testset "Day 21, Puzzle tests (post hoc)" begin
        @test day21.part1(day21.input) == 2098
        @test day21.part2(day21.input) == "ppdplc,gkcplx,ktlh,msfmt,dqsbql,mvqkdj,ggsz,hbhsx"
    end
end

using AdventOfCode2020:day20
@testset "Day 20" begin
    @testset "Day 20, Puzzle tests (post hoc)" begin
        @test_skip day20.part1(day20.input) == nothing
        @test_skip day20.part2(day20.input) == nothing
    end
end

using AdventOfCode2020:day19
@testset "Day 19" begin
    @testset "Day 19, Puzzle tests (post hoc)" begin
        @test day19.part1(day19.input) == 182
        @test day19.part2(day19.input) == 334
    end
end

using AdventOfCode2020:day18
@testset "Day 18" begin
    @testset "Day 18, Puzzle tests (post hoc)" begin
        @test day18.part1(day18.input) == 800602729153
        @test day18.part2(day18.input) == 92173009047076
    end
end

using AdventOfCode2020:day17
@testset "Day 17" begin
    @testset "Day 17, Unit tests" begin
        @test day17.part1(day17.ex) == 112
        @test day17.part2(day17.ex) == 848
    end

    @testset "Day 17, Puzzle tests (post hoc)" begin
        @test day17.part1(day17.input) == 280
        @test day17.part2(day17.input) == 1696
    end
end

using AdventOfCode2020:day16
@testset "Day 16" begin
    @testset "Day 16, Puzzle tests (post hoc)" begin
        @test day16.part1(day16.input) == 21081
        @test day16.part2(day16.input) == 314360510573
    end
end

using AdventOfCode2020:day15
@testset "Day 15" begin
    @testset "Day 15, Unit tests" begin
        @test day15.part1([3,1,2]) == 1836
        @test day15.part2([3,1,2]) == 362
    end

    @testset "Day 15, Puzzle tests (post hoc)" begin
        @test day15.part1([6,4,12,1,20,0,16]) == 475
        @test day15.part2([6,4,12,1,20,0,16]) == 11261
    end
end

using AdventOfCode2020:day14
@testset "Day 14" begin
    @testset "Day 14, Unit tests" begin
        @test day14.part1(day14.ex1) == 165
        @test day14.part2(day14.ex2) == 208
    end

    @testset "Day 14, Puzzle tests (post hoc)" begin
        @test day14.part1(day14.input) == 9967721333886
        @test day14.part2(day14.input) == 4355897790573
    end
end

using AdventOfCode2020:day13
@testset "Day 13" begin
    @testset "Day 13, Unit tests" begin
        @test day13.part1(day13.ex) == 295
        @test day13.part2(day13.ex) == 1068781
    end

    @testset "Day 13, Puzzle tests (post hoc)" begin
        @test day13.part1(day13.input) == 246
        @test day13.part2(day13.input) == 939490236001473
    end
end

using AdventOfCode2020:day12
@testset "Day 12" begin
    @testset "Day 12, Unit tests" begin
        @test day12.part1(day12.ex) == 25
        @test day12.part2(day12.ex) == 286
    end

    @testset "Day 12, Puzzle tests (post hoc)" begin
        @test day12.part1(day12.input) == 415
        @test day12.part2(day12.input) == 29401
    end
end

using AdventOfCode2020:day11
@testset "Day 11" begin
    @testset "Day 11, Unit tests" begin
        @test day11.part1(day11.ex) == 37
        @test day11.part2(day11.ex) == 26
    end

    @testset "Day 11, Puzzle tests (post hoc)" begin
        @test day11.part1(day11.input) == 2273
        @test day11.part2(day11.input) == 2064
    end
end

using AdventOfCode2020:day10
@testset "Day 10" begin
    @testset "Day 10, Unit tests" begin
        @test day10.part1(day10.ex1) == 220
        @test day10.part2(day10.ex0) == 8
        @test day10.part2(day10.ex1) == 19208
    end

    @testset "Day 10, Puzzle tests (post hoc)" begin
        @test day10.part1(day10.input) == 2176
        @test day10.part2(day10.input) == 18512297918464
    end
end

using AdventOfCode2020:day09
@testset "Day 9" begin
    @testset "Day 9, Unit tests" begin
        @test day09.first_invalid(day09.ex, 5) == 127
        @test day09.part2(day09.ex, 5) == 62
    end

    @testset "Day 9, Puzzle tests (post hoc)" begin
    @test day09.first_invalid(day09.input) == 556543474
    @test day09.part2(day09.input, 25) == 76096372
    end
end

using AdventOfCode2020:day08
@testset "Day 8" begin
    @testset "Day 8, Unit tests" begin
        @test day08.part1(day08.ex) == 5
        @test day08.part2(day08.ex) == 8
    end

    @testset "Day 8, Puzzle tests (post hoc)" begin
        @test day08.part1(day08.input) == 1930
        @test day08.part2(day08.input) == 1688
    end
end

using AdventOfCode2020:day07
@testset "Day 7" begin
    @testset "Day 7, Unit tests" begin
        @test day07.part1(day07.ex) == 4
        @test day07.part2(day07.ex) == 32
    end

    @testset "Day 7, Puzzle tests (post hoc)" begin
        @test day07.part1(day07.input) == 211
        @test day07.part2(day07.input) == 12414
    end
end

using AdventOfCode2020:day06
@testset "Day 6" begin
    @testset "Day 6, Unit tests" begin
        @test day06.count_yesses(day06.ex1)[1] == 6
        @test day06.count_yesses(day06.ex2)[1] == 11
        @test day06.count_yesses(day06.ex2)[2] == 6
    end

    @testset "Day 6, Puzzle tests (post hoc)" begin
        @test day06.count_yesses(day06.input)[1] == 6506
        @test day06.count_yesses(day06.input)[2] == 3243
    end
end

using AdventOfCode2020:day05
@testset "Day 5" begin
    @testset "Day 5, Unit tests" begin
        @test day05.parse_bp("BFFFBBFRRR") == 567
        @test day05.parse_bp("FFFBBBFRRR") == 119
        @test day05.parse_bp("BBFFBBFRLL") == 820
    end

    @testset "Day 5, Puzzle tests (post hoc)" begin
        @test day05.highest_id(day05.input) == 922
        @test day05.find_seat(day05.input)[1] == 747
    end
end

using AdventOfCode2020:day04
@testset "Day 4" begin
    @testset "Day 4, Unit tests" begin
        @test day04.count_valid_pp(day04.ex)[1] == 2
        @test day04.count_valid_pp(day04.invalid)[2] == 0
        @test day04.count_valid_pp(day04.valid)[2] == 4
    end

    @testset "Day 4, Puzzle tests (post hoc)" begin
        @test day04.count_valid_pp(day04.input)[1] == 206
        @test day04.count_valid_pp(day04.input)[2] == 123
    end
end

using AdventOfCode2020:day03
@testset "Day 3" begin
    @testset "Day 3, Unit tests" begin
        @test day03.count_trees(day03.ex) == 7
        @test day03.assess_slopes(day03.ex) == 336
    end

    @testset "Day 3, Puzzle tests (post hoc)" begin
        @test day03.count_trees(day03.input) == 189
        @test day03.assess_slopes(day03.input) == 1718180100
    end
end

using AdventOfCode2020:day02
@testset "Day 2" begin
    @testset "Day 2, Unit tests" begin
        @test day02.count_valids(day02.ex) == 2
        @test day02.count_valids_p2(day02.ex) == 1
    end

    @testset "Day 2, Puzzle tests (post hoc)" begin
        @test day02.count_valids(day02.input) == 416
        @test day02.count_valids_p2(day02.input) == 688
    end
end

using AdventOfCode2020:day01
@testset "Day 1" begin
    @testset "Day 1, Unit tests" begin
        @test day01.sum_two_to_2020(day01.ex) == 514579
        @test day01.sum_three_to_2020(day01.ex) == 241861950
    end

    @testset "Day 1, Puzzle tests (post hoc)" begin
        @test day01.sum_two_to_2020(day01.input) == 996996
        @test day01.sum_three_to_2020(day01.input) == 9210402
    end
end


text:
@testset "Diversity Measures" begin

    # A simple test to ensure the math works out.
    testP = [0.0 0.006 0.002 0.004;
             0.006 0.0 0.008 0.010;
             0.002 0.008 0.0 0.002;
             0.004 0.010 0.002 0.0;]

    testQ = [0.40, 0.20, 0.20, 0.20]

    @test NL79(testP, testQ) ≈ 0.00352

    # A simple test to ensure the math works out AND the allele
    # frequencies and so on are properly computed from sequence data.
    testSeqs = [dna"AAAACTTTTACCCCCGGGGG",
                dna"AAAACTTTTACCCCCGGGGG",
                dna"AAAACTTTTACCCCCGGGGG",
                dna"AAAACTTTTACCCCCGGGGG",
                dna"AAAAATTTTACCCCCGTGGG",
                dna"AAAAATTTTACCCCCGTGGG",
                dna"AAAACTTTTTCCCCCGTAGG",
                dna"AAAACTTTTTCCCCCGTAGG",
                dna"AAAAATTTTTCCCCCGGAGG",
                dna"AAAAATTTTTCCCCCGGAGG"]

    testP = [0.0 0.1 0.15 0.15;
             0.1 0.0 0.15 0.15;
             0.15 0.15 0.0 0.1;
             0.15 0.15 0.1 0.0;]

    @test NL79(testSeqs) ≈ NL79(testP, testQ)

    testSeqs = [dna"ATAATAAAAAAATAATAAAAAAATAAAAAAAATAAAAAAAA",
                dna"AAAAAAAATAAATAATAAAAAAATAAAAAAAAAAAAAAAAA",
                dna"AAAATAAAAATATAATAAAAAAATATAAAAAAAAAAAAAAA",
                dna"AAAAAAAAAAAATAATAAAAAAATAAATAAATAAAAAAAAA",
                dna"AAAATAAAAAAAATATAAAAAAATAAAAAAAAAAAAAAAAA",
                dna"AAAATAAAAAAAAAATAAAAAAAAAAAAAAAAAAATAAAAA",
                dna"AAAAAATAAAAATAATAAAAAAATAAAAAAAAAAAAAAAAA",
                dna"AAAAAAAAAAAAAAATAAAAAAATAAAAAAAAAAAAAAATA",
                dna"AAAAAAAAAAAAAAAAAAAAAAATAAAAAAAAAAAAAAAAA",
                dna"AAAAAAAAAAAAAAATAAAAAAATAATAAAAAAAAAAAAAA"]

    @test avg_mut(testSeqs) ≈ 3.888888 atol=10e-5
end


text:
using WAVEFORM, LinearAlgebra

g1 = range(0, stop=1, length=100)
g2 = range(0, stop=1, length=200)
I = joLagrangeInterp1D(g1, g2, Float64)

x = randn(100)
y = randn(200)

# Dot test
@test dot(y, I*x) - dot(x, I'*y) < 1e-9


text:
using LightGraphs
using LinearAlgebra
using StatsBase
using GraphPlot
using Optim
using ReverseDiff
using ForwardDiff

include("src/estimating_functions.jl")
include("src/generators.jl")

n = 50
beta = zeros(n)
for i in sample(collect(1:n), Int(floor(n / 10)))
    beta[i] = sqrt(log(n))
end
mu = -log(n)

g = sparse_beta_graph(beta, mu)

layout = (args...) -> spring_layout(args...; C=9)
gplot(g, layout=layout)

# Parameter Estimation

function sparse_beta_loglik(g::Graph, mu, beta)
    beta = exp.(log.(beta))
    tot = 0
    for i in 1:n
        tot -= degree(g)[i] * beta[i]
        for j in i:n
            tot += log(1 + exp(mu + beta[i] + beta[j]))
        end
    end
    return -(g.ne * mu) + tot
end

function sparse_beta_loglik_rp(g::Graph, mu, gamma)
    beta = exp.(gamma)
    tot = 0
    for i in 1:n
        tot -= degree(g)[i] * beta[i]
        for j in i:n
            tot += log(1 + exp(mu + beta[i] + beta[j]))
        end
    end
    return -(g.ne * mu) + tot
end

f = mubeta -> sparse_beta_loglik_rp(g, mubeta[1], mubeta[2:end])

tape = ReverseDiff.GradientTape(f, vcat(-1.5, repeat([1], n)))
(storage, params) -> ReverseDiff.gradient!(storage, tape, params)

optimize(
    mubeta -> sparse_beta_loglik_rp(g, mubeta[1], mubeta[2:end]),
    (storage, params) -> ReverseDiff.gradient!(storage, tape, params),
    # repeat([0], n+1),
    # repeat([Inf], n+1),
    vcat(-1.5, repeat([1], n)),
    LBFGS(),
    Optim.Options(show_every=1))

# # optimize(mubeta -> sparse_beta_loglik(g, mubeta[1], mubeta[2:end]), vcat(-1.5, repeat([0], n)))

# tape = ReverseDiff.GradientTape(mubeta -> sparse_beta_loglik(g, mubeta[1], mubeta[2:end]), vcat(-1.5, repeat([0], n)))
# f(mubeta) = sparse_beta_loglik(g, mubeta[1], mubeta[2:end])
# f_tape = ReverseDiff.GradientTape(f, vcat(-1.5, repeat([0], n)))
# g!(storage, params) = ReverseDiff.gradient!(storage, f_tape, params)

# optimize(
#     f,
#     g!(storage, params),
#     vcat(-1.5, repeat([0], n)),
#     LBFGS(),
#     Optim.Options(
#         show_every=1))

max_stars_repo_name: jacobusmmsmit/sparse-beta

max_issues_repo_name: jacobusmmsmit/sparse-beta

max_forks_repo_name: jacobusmmsmit/sparse-beta


text:
module LeftChildRightSiblingTrees

# See `abstracttrees.jl` for the only dependency of this package

export Node,
    addchild,
    addsibling,
    depth,
    graftchildren!,
    isroot,
    isleaf,
    lastsibling,
    prunebranch!

mutable struct Node{T}
    data::T
    parent::Node{T}
    child::Node{T}
    sibling::Node{T}

    """
        root = Node(data)

    Construct a disconnected node, which can serve as the root of a new tree.
    `root.data` holds `data`.
    """
    function Node{T}(data) where T
        n = new{T}(data)
        n.parent = n
        n.child = n
        n.sibling = n
        n
    end

    """
        node = Node(data, parent::Node)

    Construct a `node` with `parent` as its parent node. `node.data` stores `data`.
    Node that this does *not* update any links in, e.g., `parent`'s other children.

    For a higher-level interface, see [`addchild`](@ref).
    """
    function Node{T}(data, parent::Node) where T
        n = new{T}(data, parent)
        n.child = n
        n.sibling = n
        n
    end
end
Node(data::T) where {T} = Node{T}(data)
Node(data, parent::Node{T}) where {T} = Node{T}(data, parent)

"""
    node = lastsibling(child)

Return the last sibling of `child`.
"""
function lastsibling(sib::Node)
    newsib = sib.sibling
    while sib != newsib
        sib = newsib
        newsib = sib.sibling
    end
    sib
end

"""
    node = addsibling(oldersib, data)

Append a new "youngest" sibling, storing `data` in `node`. `oldersib` must be
the previously-youngest sibling (see [`lastsibling`](@ref)).
"""
function addsibling(oldersib::Node{T}, data) where T
    if oldersib.sibling != oldersib
        error("Truncation of sibling list")
    end
    youngersib = Node(data, oldersib.parent)
    oldersib.sibling = youngersib
    youngersib
end

"""
    node = addchild(parent::Node, data)

Create a new child of `parent`, storing `data` in `node.data`.
This adjusts all links to ensure the integrity of the tree.
"""
function addchild(parent::Node{T}, data) where T
    newc = Node(data, parent)
    prevc = parent.child
    if prevc == parent
        parent.child = newc
    else
        prevc = lastsibling(prevc)
        prevc.sibling = newc
    end
    newc
end

"""
    isroot(node)

Returns `true` if `node` is the root of a tree (meaning, it is its own parent).
"""
isroot(n::Node) = n == n.parent

"""
    isleaf(node)

Returns `true` if `node` has no children.
"""
isleaf(n::Node) = n == n.child

makeleaf!(n::Node) = n.child = n

Base.show(io::IO, n::Node) = print(io, "Node(", n.data, ')')

# Iteration over children
# for c in parent
#     # do something
# end
Base.IteratorSize(::Type{<:Node}) = Base.SizeUnknown()
Base.eltype(::Type{Node{T}}) where T = Node{T}

function Base.iterate(n::Node, state::Node = n.child)
    n === state && return nothing
    return state, state === state.sibling ? n : state.sibling
end

# To support Base.pairs
struct PairIterator{T}
    parent::Node{T}
end
Base.pairs(node::Node) = PairIterator(node)
Base.IteratorSize(::Type{<:PairIterator}) = Base.SizeUnknown()

function Base.iterate(iter::PairIterator, state::Node=iter.parent.child)
    iter.parent === state && return nothing
    return state=>state, state === state.sibling ? iter.parent : state.sibling
end

function showedges(io::IO, parent::Node, printfunc = identity)
    str = printfunc(parent.data)
    if str != nothing
        if isleaf(parent)
            println(io, str, " has no children")
        else
            print(io, str, " has the following children: ")
            for c in parent
                print(io, printfunc(c.data), "    ")
            end
            print(io, "\n")
            for c in parent
                showedges(io, c, printfunc)
            end
        end
    end
end
showedges(parent::Node) = showedges(stdout, parent)

depth(node::Node) = depth(node, 1)
function depth(node::Node, d)
    childd = d + 1
    for c in node
        d = max(d, depth(c, childd))
    end
    return d
end

"""
    graftchildren!(dest, src)

Move the children of `src` to become children of `dest`.
`src` becomes a leaf node.
"""
function graftchildren!(dest, src)
    for c in src
        c.parent = dest
    end
    if isleaf(dest)
        dest.child = src.child
    else
        lastsib = lastsibling(dest.child)
        lastsib.sibling = src.child
    end
    makeleaf!(src)
    return dest
end

"""
    prunebranch!(node)

Eliminate `node` and all its children from the tree.
"""
function prunebranch!(node)
    isroot(node) && error("cannot prune the root")
    p = node.parent
    if p.child == node
        # `node` is the first child of p
        if node.sibling === node
            p.child = p   # p is now a leaf
        else
            p.child = node.sibling
        end
    else
        # `node` is a middle or last child of p
        child = p.child
        sib = child.sibling
        while sib != node
            @assert sib != child
            child = sib
            sib = child.sibling
        end
        if sib.sibling === sib
            # node is the last child of p, just truncate
            child.sibling = child
        else
            # skip over node
            child.sibling = sib.sibling
        end
    end
    return p
end

include("abstracttrees.jl")

end # module


text:
using .Plots

function preprocess(name, plot::Plots.Plot, data)
#=     pb = PipeBuffer()
    show(pb, MIME("image/png"), plot)
    arr = FileIO.load(pb)
 =#    push!(data, name=>plot)
    return data
end

function preprocess(name, plots::AbstractArray{<:Plots.Plot}, data)
    for (i, plot)=enumerate(plots)
        preprocess(name*"/$i", plot, data)
    end
    return data
end

function log_plot(lg::MLFLogger, key::AbstractString, obj::Plots.Plot; step=nothing)
    mktempdir(;prefix="jlmlf_") do dirname
        fname = joinpath(dirname, key)
        mkpath(fname)
        fname = joinpath(fname, "plot.png")
        savefig(obj, fname)
        log_image(lg, key, fname; step=step)
    end

end

process(lg::MLFLogger, name::AbstractString, obj::Plots.Plot, step::Int) = log_plot(lg, name, obj; step=step)


text:
function shock_loading(m::RealBondMkup, TTT_jump::Matrix{Float64})
    exo = m.exogenous_shocks
    endo = m.endogenous_states

    _RRR = zeros(n_backward_looking_states(m), n_shocks_exogenous(m))
    RRR  = zeros(n_model_states(m), n_shocks_exogenous(m))

    # Technology shock loading on states
    _RRR[endo[:z′_t], exo[:z_sh]] .= 1.

    # Monetary policy shock loading on states
    _RRR[endo[:mon′_t], exo[:mon_sh]] .= 1.

    # Markup shock loading on states
    _RRR[endo[:mkp′_t], exo[:mkp_sh]] .= 1.

    # Loading on states and jumps
    RRR[1:n_backward_looking_states(m), :] = _RRR
    RRR[n_backward_looking_states(m)+1:end, :] = TTT_jump*_RRR

    return RRR
end


text:
struct Layers
    layers::Vector{Layer}
end

Base.convert(::Type{Layers}, s::Layer) = Layers([s])
Base.convert(::Type{Layers}, s::Layers) = s

Base.getindex(v::Layers, i::Int) = v.layers[i]
Base.length(v::Layers) = length(v.layers)
Base.eltype(::Type{Layers}) = Layer
Base.iterate(v::Layers, args...) = iterate(v.layers, args...)

const OneOrMoreLayers = Union{Layers, Layer}

function Base.:+(s1::OneOrMoreLayers, s2::OneOrMoreLayers)
    l1::Layers, l2::Layers = s1, s2
    return Layers(vcat(l1.layers, l2.layers))
end

function Base.:*(s1::OneOrMoreLayers, s2::OneOrMoreLayers)
    l1::Layers, l2::Layers = s1, s2
    return Layers([el1 * el2 for el1 in l1 for el2 in l2])
end

function summary(e::Entry)
    summaries = Dict{KeyType, Any}()
    for (i, tup) in enumerate((e.primary, e.positional, e.named))
        for (key, val) in pairs(tup)
            summaries[key] = if i > 1 && iscontinuous(val)
                AbstractPlotting.extrema_nan(val)
            elseif i == 2 && isgeometry(val)
                nothing
            else
                collect(uniquesorted(vec(val)))
            end
        end
    end
    return Entry(e; summaries)
end

mergesummaries(s1::AbstractVector, s2::AbstractVector) = mergesorted(s1, s2)
mergesummaries(s1::Tuple, s2::Tuple) = extend_extrema(s1, s2)
mergesummaries(::Nothing, ::Nothing) = nothing

mergelabels(a, b) = a

function compute_axes_grid(fig, s::OneOrMoreLayers;
                           axis=NamedTuple(), palettes=NamedTuple())
    layers::Layers = s
    entries = [summary(e) for layer in layers for e in process(layer)]
    summaries = mapfoldl(
        entry -> entry.summaries,
        mergewith!(mergesummaries),
        entries,
        init=Dict{KeyType, Any}()
    )
    labels = mapfoldl(
        entry -> entry.labels,
        mergewith!(mergelabels),
        entries,
        init=Dict{KeyType, Any}()
    )

    palettes = merge(default_palettes(), palettes)
    scales = default_scales(summaries, palettes)

    e = (; entries, scales, labels)

    rowcol = (:row, :col)

    layout_scale, scales... = map((:layout, rowcol...)) do sym
        return get(e.scales, sym, nothing)
    end

    grid_size = map(scales, (first, last)) do scale, f
        isnothing(scale) || return maximum(scale.plot)
        isnothing(layout_scale) || return maximum(f, layout_scale.plot)
        return 1
    end

    axes_grid = map(CartesianIndices(grid_size)) do c
        type = get(axis, :type, Axis)
        options = Base.structdiff(axis, (; type))
        ax = type(fig[Tuple(c)...]; options...)
        return AxisEntries(ax, Entry[], e.scales, e.labels)
    end

    for entry in e.entries
        rows, cols = map(rowcol, scales, (first, last)) do sym, scale, f
            v = get(entry.primary, sym, nothing)
            layout_v = get(entry.primary, :layout, nothing)
            # without layout info, plot on all axes
            # all values in `v` and `layout_v` are equal
            isnothing(v) || return rescale(v[1:1], scale)
            isnothing(layout_v) || return map(f, rescale(layout_v[1:1], layout_scale))
            return 1:f(grid_size)
        end
        for i in rows, j in cols
            ae = axes_grid[i, j]
            push!(ae.entries, entry)
        end
    end

    # Link colors
    labeledcolorbar = getlabeledcolorbar(axes_grid)
    if !isnothing(labeledcolorbar)
        _, colorbar = labeledcolorbar
        colorrange = colorbar.extrema
        for entry in AlgebraOfGraphics.entries(axes_grid)
            entry.attributes[:colorrange] = colorrange
        end
    end

    return axes_grid

end

function AbstractPlotting.plot!(fig, s::OneOrMoreLayers;
                                axis=NamedTuple(), palettes=NamedTuple())
    grid = compute_axes_grid(fig, s; axis, palettes)
    foreach(plot!, grid)
    return grid
end

function AbstractPlotting.plot(s::OneOrMoreLayers;
                               axis=NamedTuple(), figure=NamedTuple(), palettes=NamedTuple())
    fig = Figure(; figure...)
    grid = plot!(fig, s; axis, palettes)
    return FigureGrid(fig, grid)
end

# Convenience function, which may become superfluous if `plot` also calls `facet!`
function draw(s::OneOrMoreLayers;
              axis = NamedTuple(), figure=NamedTuple(), palettes=NamedTuple())
    fg = plot(s; axis, figure, palettes)
    facet!(fg)
    Legend(fg)
    resizetocontent!(fg)
    return fg
end

function draw!(fig, s::OneOrMoreLayers;
               axis=NamedTuple(), palettes=NamedTuple())
    ag = plot!(fig, s; axis, palettes)
    facet!(fig, ag)
    return ag
end


text:
"""
Contains all found `Branch`es
"""
type Project
	branches::Vector{Branch}	# unordered
	activeSolution::Union{Solution, Vector{Float64}, Type{Void}}
	O::Observer
	# lock::ReentrantLock
	Project() = new(Branch[], Void, Observer())
	Project(B::Vector{Branch}) = new(B, Void, Observer())
end
@relay(Project, :branches)

#prevent enormous amounts of data being printed
Base.show(io::Base.IO, P::Project) = write(io, "***PROJECT***")

#TODO lots of functions missing... add? pop e.g. really unlikely op... -> add on demand

# Base.lock(P::Project) = lock(P.lock)
# Base.unlock(P::Project) = unlock(P.lock)

mbObserve.observe(cb::Function, P::Project, s) = observe(cb, P.O, s)
mbObserve.notify(P::Project, s, v...) = mbObserve.notify(P.O, s, v...)

function Base.push!(P::Project, branches...)
	for b in branches
		push!(P.branches, b)
		mbObserve.notify(P, :pushBranch, b)
	end
	return P
end

function Base.deleteat!(P::Project, i::Int)
	B = P[i]
	P.activeSolution in B && setActiveSolution(P, P.activeSolution.data)
	deleteat!(P.branches, i)
	mbObserve.notify(P, :delBranch, B)
end

function setActiveSolution(P::Project, V::Vector{Float64})
	P.activeSolution = V
	mbObserve.notify(P, :activeSolutionChanged)
	return Void
end

function setActiveSolution(P::Project, S::Solution)
	B = S.parent
	P = B.parent
	P.activeSolution = S
	mbObserve.notify(P, :activeSolutionChanged)
	return Void
end


text:
function layoutable(::Type{Box}, fig_or_scene; bbox = nothing, kwargs...)

    topscene = get_topscene(fig_or_scene)

    default_attrs = default_attributes(Box, topscene).attributes
    theme_attrs = subtheme(topscene, :Box)
    attrs = merge!(merge!(Attributes(kwargs), theme_attrs), default_attrs)

    @extract attrs (color, visible, valign, halign, padding, strokewidth,
        strokevisible, strokecolor)

    layoutobservables = LayoutObservables{Box}(attrs.width, attrs.height, attrs.tellwidth, attrs.tellheight,
        halign, valign, attrs.alignmode; suggestedbbox = bbox)

    strokecolor_with_visibility = lift(strokecolor, strokevisible) do col, vis
        vis ? col : RGBAf0(0, 0, 0, 0)
    end

    ibbox = @lift(round_to_IRect2D($(layoutobservables.computedbbox)))

    r = poly!(topscene, ibbox, color = color, visible = visible, raw = true,
        strokecolor = strokecolor_with_visibility, strokewidth = strokewidth,
        inspectable = false)

    elements = Dict(:rect => r)

    # trigger bbox
    layoutobservables.suggestedbbox[] = layoutobservables.suggestedbbox[]

    Box(fig_or_scene, layoutobservables, attrs, elements)
end


text:
struct RRVector
end

struct RRStringArray
end

struct RRDoubleMatrix
end

struct RRComplex
end

struct RRComplexVector
end

struct RRComplexMatrix
end

struct RRCData
end

struct RRList
end

struct RRListItem
end


text:
"""
   Traits for the NEID spectrograph (data products from DRP v0.6).
   https://neid.psu.edu/
Author: Eric Ford and collaborators
Created: August 2020
"""

""" Delegates loading of code specifying types essential to the package.  """

import RvSpectMLBase: min_order, max_order, min_pixel_in_order, max_pixel_in_order, min_pixel, max_pixel

#= Pre-ship
min_order(::NEID2D) = 1
max_order(::NEID2D) = 90
=#
# At KPNO
min_order(::NEID2D) = 1
#max_order(::NEID2D) = 126 # DRS 0.5
max_order(::NEID2D) = 122 # DRS 0.6
min_pixel_in_order(inst::NEID2D) = 1
max_pixel_in_order(inst::NEID2D) = 9216

min_pixel(::NEID1D) = 1
max_pixel(::NEID1D) = (max_order(NEID1D())-min_order(NEID1D())+1)*9216 # TODO: Update once know size of NEID's 1d extracted spectra

bad_col_ranges(inst::NEID2D, ord::Int) = [439:449, 1934:1943, 6714:6714]

import RvSpectMLBase: orders_to_use_default, min_col_default, max_col_default

#= Pre-ship
#orders_to_use_default(inst::NEID2D) = 1:52   # Avoiding redder orders due to tellurics
orders_to_use_default(inst::NEID2D) = 1:71   # Avoiding 72 because of NaNs in solar data
min_col_default(::NEID2D, ord::Integer) = 451              # Avoiding smaller columns due to NaNs
#min_col_default(::NEID2D) = 2000              # Avoiding smaller columns due to lower flux and distortions
max_col_default(::NEID2D, ord::Integer) = 9216 - (min_col_default(NEID2D(),ord)-1)   # Avoiding larger columns for symmetry
=#
# At KPNO
#orders_to_use_default(::NEID2D) = 55:113 # everything plaussibly usable as of DRS 0.5
#orders_to_use_default(::NEID2D) = 56:111 # everything plaussibly usable as of DRS 0.6
#orders_to_use_default(::NEID2D) = 60:90  # relatively safe
orders_to_use_default(::NEID2D) = 56:98   # avoid worst of tellurics
function min_col_default(::NEID2D, ord::Integer)
#    return 1445 +500 # DRS 0.6, avoiding NaN in cols 1435-1444
    if ord == 55
        return 787
    else
        return 500
    end
end
#max_col_default(::NEID2D, ord::Integer) = 8429  # DRS 0.5
#max_col_default(::NEID2D, ord::Integer) = 6214  # DRS 0.6, avoiding NaN in col 6215
max_col_default(::NEID2D, ord::Integer) = 9215  # DRS 0.6

import RvSpectMLBase: get_pixel_range
function get_pixel_range(inst::NEID2D, ord::Integer)
    minc = max(min_col_default(inst, ord)) #, min_col_excalibur(inst,order), min_col_nonnan(inst,order))
    maxc = min(max_col_default(inst, ord)) #, max_col_excalibur(inst,order), max_col_nonnan(inst,order))
    return minc:maxc
end

import RvSpectMLBase: metadata_symbols_default, metadata_strings_default
#= Pre-ship headers
#metadata_symbols_default(::AnyNEID) = Symbol[:bjd, :target, :ssbz]
#metadata_strings_default(::AnyNEID) = String["OBSJD", "SKY-OBJ", "SSBZ000"]
=#
# Headers at KPNO
metadata_symbols_default(::AnyNEID) = Symbol[:bjd, :target, :exptime, :airmass, :ssbz]
metadata_strings_default(::AnyNEID) = String["OBSJD", "OBJECT", "EXPTIME", "AIRMASS",  "SSBZ100"]

import RvSpectMLBase: default_ccf_mask_v_width
default_ccf_mask_v_width(::AnyNEID) = 620.953

import RvSpectMLBase: get_inst_module
get_inst_module(::AnyNEID) = NEID

import RvSpectMLBase: get_λ_range
function get_λ_range(data::CLT) where { T1<:Real, T2<:Real, T3<:Real, A1<:AbstractArray{T1,2}, A2<:AbstractArray{T2,2}, A3<:AbstractArray{T3,2},
                                       IT<:AnyNEID, CLT<:Spectra2DBasic{T1,T2,T3,A1,A2,A3,IT} }
   (λmin, λmax) = NaNMath.extrema(data.λ)
   return (min=λmin, max=λmax)
end

#=
# Pre-ship
default_λmin = 3950.0  # Based on HD solar data from PSU, should generalize
default_λmax = 9500.0  #
=#
# At KPNO
#default_λmin = 5040.04045429369  # Based on one FITS file, should generalize
#default_λmax = 9810.539033567005 # Based on one FITS file, should generalize
default_λmin = 5145.0
default_λmax = 9822.0

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
include("../jlproto/RaftRPC.jl")

###
# Types for state stored on servers
###

# Lots of different indexes in this format
typealias CandidateId UInt64
typealias LogIndex UInt64
typealias Term UInt64
typealias Server UInt64 # not sure what to do with this?

# Remember to save to disk before responding to RPCs
type PersistentState
  currentTerm::Term
  votedFor::CandidateId
  log::Array{Tuple{Term, RaftRPC.LogEntry}}
end

type VolatileState
  commitIndex::LogIndex
  lastApplied::LogIndex
end

type LeaderVolatileState
  nextIndex::Dict{Server, LogIndex}
  matchIndex::Dict{Server, LogIndex}
end

# Configuration
type Address
  ipaddress::IPv4
  port::UInt
end

type ServerConfig
  id::Server # my id
  port::UInt # port to listen on
  servers::Dict{Server, Address}
end


text:
#using .Bio3DView

function get_residues(ta::TrjArray)
    natom = ta.natom

    properties = Dict()
    if isempty(ta.chainid)
        properties["chainid"] = [nothing]
    else
        properties["chainid"] = unique(ta.chainid)
    end
    if isempty(ta.chainname)
        properties["chainname"] = [nothing]
    else
        properties["chainname"] = unique(ta.chainname)
    end
    if isempty(ta.resid)
        properties["resid"] = [nothing]
    else
        properties["resid"] = unique(ta.resid)
    end
    if isempty(ta.resname)
        properties["resname"] = [nothing]
    else
        properties["resname"] = unique(ta.resname)
    end

    res_array = []
    id = BitArray(undef, natom)
    id_c1 = BitArray(undef, natom)
    id_c2 = BitArray(undef, natom)
    id_r1 = BitArray(undef, natom)
    id_r2 = BitArray(undef, natom)
    for c1 in properties["chainid"]
        id .= id_c1 .= id_c2 .= id_r1 .= id_r2 .= 1
        if !isnothing(c1)
            id_c1 .= (id .& Bool(ta.chainid .== c1))
        end
        for c2 in properties["chainname"]
            if !isnothing(c2)
                id_c2 .= (id_c1 .& (ta.chainname .== c2))
            end
            for r1 in properties["resid"]
                if !isnothing(r1)
                    id_r1 .= (id_c2 .& (ta.resid .== r1))
                end
                for r2 in properties["resname"]
                    if !isnothing(r2)
                        id_r2 .= (id_r1 .& (ta.resname .== r2))
                    end
                    if any(id_r2)
                        push!(res_array, ta[:, id_r2])
                    end
                end
            end
        end
    end

    return res_array
end

#"""
#viewstruct
#
#based on Bio3Dview.jl
#"""
#function viewstruc(ta::TrjArray; kwargs...)
#    io = IOBuffer()
#    writepdb(io, ta)
#    return Bio3DView.view("data-type='pdb'", String(take!(io)); style=Bio3DView.defaultstyle("pdb"), kwargs...)
#end

#cpu(m) = fmap(x -> adapt(Array, x), m)
#gpu(x) = use_cuda[] ? fmap(CUDA.cu, x) : x

function gpu(ta::TrjArray{T, U}) where {T, U}
    TrjArray{T, U}(CuArray(ta.xyz), CuArray(ta.boxsize),
      ta.chainname, CuArray(ta.chainid),
      ta.resname, CuArray(ta.resid),
      ta.atomname, CuArray(ta.atomid),
      CuArray(ta.mass), CuArray(ta.radius), CuArray(ta.charge), CuArray(ta.sasa),
      CuArray(ta.list_bond), CuArray(ta.list_angle), CuArray(ta.list_dihedral),
      CuArray(ta.list_improper), CuArray(ta.list_cmap))
end

function gpu32(ta::TrjArray)
    TrjArray{Float32, Int64}(CUDA.cu(ta.xyz), CUDA.cu(ta.boxsize),
      ta.chainname, CUDA.cu(ta.chainid),
      ta.resname, CUDA.cu(ta.resid),
      ta.atomname, CUDA.cu(ta.atomid),
      CUDA.cu(ta.mass), CUDA.cu(ta.radius), CUDA.cu(ta.charge), CUDA.cu(ta.sasa),
      CUDA.cu(ta.list_bond), CUDA.cu(ta.list_angle), CUDA.cu(ta.list_dihedral),
      CUDA.cu(ta.list_improper), CUDA.cu(ta.list_cmap))
end

function cpu(ta::TrjArray{T, U}) where {T, U}
    TrjArray{T, U}(Array(ta.xyz), Array(ta.boxsize),
      ta.chainname, Array(ta.chainid),
      ta.resname, Array(ta.resid),
      ta.atomname, Array(ta.atomid),
      Array(ta.mass), Array(ta.radius), Array(ta.charge), Array(ta.sasa),
      Array(ta.list_bond), Array(ta.list_angle), Array(ta.list_dihedral),
      Array(ta.list_improper), Array(ta.list_cmap))
end

function cpu64(ta::TrjArray)
    TrjArray{Float64, Int64}(Array{Float64}(ta.xyz), Array{Float64}(ta.boxsize),
      ta.chainname, Array{Int64}(ta.chainid),
      ta.resname, Array{Int64}(ta.resid),
      ta.atomname, Array{Int64}(ta.atomid),
      Array{Float64}(ta.mass), Array{Float64}(ta.radius), Array{Float64}(ta.charge), Array{Float64}(ta.sasa),
      Array{Int64}(ta.list_bond), Array{Int64}(ta.list_angle), Array{Int64}(ta.list_dihedral),
      Array{Int64}(ta.list_improper), Array{Int64}(ta.list_cmap))
end

"""
    logsumexp(X)

Compute `log(sum(exp, X))` avoiding possible overflows/undeflows.
`X` should be an iterator of real numbers.
"""
function logsumexp(X)
    isempty(X) && return log(sum(X))
    reduce(logaddexp, X)
end

function logsumexp(X::AbstractArray{T}; dims=:) where {T<:Real}
    # Do not use log(zero(T)) directly to avoid issues with ForwardDiff (#82)
    u = reduce(max, X, dims=dims, init=oftype(log(zero(T)), -Inf))
    u isa AbstractArray || isfinite(u) || return float(u)
    let u=u # avoid https://github.com/JuliaLang/julia/issues/15276
        # TODO: remove the branch when JuliaLang/julia#31020 is merged.
        if u isa AbstractArray
            u .+ log.(sum(exp.(X .- u); dims=dims))
        else
            u + log(sum(x -> exp(x-u), X))
        end
    end
end


text:
using Revise, YAML, ApproxOperator

config = YAML.load_file("test.yml")

# elements,nodes = importmsh("./msh/bar.msh",config)
#
# sp = RegularGrid(nodes[:x],nodes[:y],nodes[:z],n = 2,γ = 1)
#
# sp(elements["Ω"],elements["Γᵍ"])

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
"""
```
function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matrix{S};
             kwargs...) where {S<:AbstractFloat, U<:Number}
```

### Arguments:

- `loglikelihood::Function`: Log-likelihood function of model being estimated. Takes `parameters`
    and `data` as arguments.
- `parameters::ParameterVector{U}`: Model parameter vector, which stores parameter values,
    prior dists, and bounds.
- `data`: A matrix or dataframe containing the time series of the observables used in
    the calculation of the posterior/loglikelihood
- `old_data`: A matrix containing the time series of observables of previous data
    (with `data` being the new data) for the purposes of a time tempered estimation
    (that is, using the posterior draws from a previous estimation as the initial set
    of draws for an estimation with new data)

### Keyword Arguments:
- `verbose::Symbol`: Desired frequency of function progress messages printed to standard out.
	- `:none`: No status updates will be reported.
	- `:low`: Status updates for SMC initialization and recursion will be included.
	- `:high`: Status updates for every iteration of SMC is output, which includes
    the mean and standard deviation of each parameter draw after each iteration,
    as well as calculated acceptance rate, ESS, and number of times resampled.

- `parallel::Bool`: Flag for running algorithm in parallel.
- `n_parts::Int`: Number of particles.
- `n_blocks::Int`: Number of parameter blocks in mutation step.
- `n_mh_steps::Int`: Number of Metropolis Hastings steps to attempt during the mutation step.
- `λ::S`: The 'bending coefficient' λ in Φ(n) = (n/N(Φ))^λ
- `n_Φ::Int`: Number of stages in the tempering schedule.
- `resampling_method::Symbol`: Which resampling method to use.
    - `:systematic`: Will use sytematic resampling.
    - `:multinomial`: Will use multinomial resampling.
    - `:polyalgo`: Samples using a polyalgorithm.
- `threshold_ratio::S`: Threshold s.t. particles will be resampled when the population
    drops below threshold * N.
- `c::S`: Scaling factor for covariance of the particles. Controls size of steps in mutation step.
- `α::S`: The mixture proportion for the mutation step's proposal distribution.
- `target::S`: The initial target acceptance rate for new particles during mutation.

- `use_chand_recursion::Bool`: Flag for using Chandrasekhar Recursions in Kalman filter.
- `use_fixed_schedule::Bool`: Flag for whether or not to use a fixed tempering (ϕ) schedule.
- `tempering_target::S`: Coefficient of the sample size metric to be targeted when solving
    for an endogenous ϕ.

- `old_data::Matrix{S}`: data from vintage of last SMC estimation.
- `old_cloud::Cloud`: associated cloud borne of old data in previous SMC estimation.
- `old_vintage::String`: String for vintage date of old data
- `smc_iteration::Int`: The iteration index for the number of times SMC has been run on the
     same data vintage. Primarily for numerical accuracy/testing purposes.

- `run_test::Bool`: Flag for when testing accuracy of program
- `filestring_addl::Vector{String}`: Additional file string extension for loading old cloud.
- `continue_intermediate::Bool`: Flag to indicate whether one is continuing SMC from an
    intermediate stage.
- `intermediate_stage_start::Int`: Intermediate stage at which one wishes to begin the estimation.
- `save_intermediate::Bool`: Flag for whether one wants to save intermediate Cloud objects
- `intermediate_stage_increment::Int`: Save Clouds at every increment
    (1 = each stage, 10 = every 10th stage, etc.)

### Outputs

- `cloud`: The Cloud object containing all of the information about the
    parameter values from the sample, their respective log-likelihoods, the ESS
    schedule, tempering schedule etc., which is saved in the saveroot.

### Overview

Sequential Monte Carlo can be used in lieu of Random Walk Metropolis Hastings to
    generate parameter samples from high-dimensional parameter spaces using
    sequentially constructed proposal densities to be used in iterative importance
    sampling.

This implementation is based on Edward Herbst and Frank Schorfheide's 2014 paper
    'Sequential Monte Carlo Sampling for DSGE Models' and the code accompanying their
    book 'Bayesian Estimation of DSGE Models'.

SMC is broken up into three main steps:

- `Correction`: Reweight the particles from stage n-1 by defining incremental weights,
    which gradually "temper in" the loglikelihood function p(Y|θ)^(ϕ_n - ϕ_n-1) into the
    normalized particle weights.
- `Selection`: Resample the particles if the distribution of particles begins to
    degenerate, according to a tolerance level for the ESS.
- `Mutation`: Propagate particles {θ(i), W(n)} via N(MH) steps of a Metropolis
    Hastings algorithm.
"""
function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matrix{S};
             verbose::Symbol = :low,
             testing::Bool   = false,
             data_vintage::String = Dates.format(today(), "yymmdd"),

             parallel::Bool  = false,
             n_parts::Int    = 5_000,
             n_blocks::Int   = 1,
             n_mh_steps::Int = 1,

             λ::S = 2.1,
             n_Φ::Int64 = 300,

             resampling_method::Symbol = :systematic,
             threshold_ratio::S = 0.5,

             # Mutation Settings
             c::S = 0.5,       # step size
             α::S = 1.0,       # mixture proportion
             target::S = 0.25, # target accept rate

             use_fixed_schedule::Bool = true,
             tempering_target::S = 0.97,

             old_data::Matrix{S} = Matrix{S}(undef, size(data, 1), 0),
             old_cloud::Cloud = Cloud(0, 0),
             old_vintage::String = "",
             smc_iteration::Int = 1,

             run_test::Bool = false,
             filestring_addl::Vector{String} = Vector{String}(),
             loadpath::String = "",
             savepath::String = "smc_cloud.jld2",
             particle_store_path::String = "smcsave.h5",
             continue_intermediate::Bool = false,
             intermediate_stage_start::Int = 0,
             save_intermediate::Bool = false,
             intermediate_stage_increment::Int = 10) where {S<:AbstractFloat, U<:Number}

    ########################################################################################
    ### Settings
    ########################################################################################

    # Construct closure of mutation function so as to avoid issues with serialization
    # across workers with different Julia system images
    sendto(workers(), parameters = parameters)
    sendto(workers(), data = data)

    function mutation_closure(p::Vector{S}, d_μ::Vector{S}, d_Σ::Matrix{S},
                 blocks_free::Vector{Vector{Int64}}, blocks_all::Vector{Vector{Int64}},
                 ϕ_n::S, ϕ_n1::S; c::S = 1.0, α::S = 1.0, n_mh_steps::Int = 1,
                 old_data::T = Matrix{S}(undef, size(data, 1), 0)) where {S<:Float64, T<:Matrix}
        return mutation(loglikelihood, parameters, data, p, d_μ, d_Σ, blocks_free, blocks_all,
                        ϕ_n, ϕ_n1; c = c, α = α, n_mh_steps = n_mh_steps, old_data = old_data)
    end
    @everywhere function mutation_closure(p::Vector{S}, d_μ::Vector{S}, d_Σ::Matrix{S},
                 blocks_free::Vector{Vector{Int64}}, blocks_all::Vector{Vector{Int64}},
                 ϕ_n::S, ϕ_n1::S; c::S = 1.0, α::S = 1.0, n_mh_steps::Int = 1,
                 old_data::T = Matrix{S}(undef, size(data, 1), 0)) where {S<:Float64, T<:Matrix}
        return mutation(loglikelihood, parameters, data, p, d_μ, d_Σ, blocks_free, blocks_all,
                        ϕ_n, ϕ_n1; c = c, α = α, n_mh_steps = n_mh_steps, old_data = old_data)
    end

    # Check that if there's a tempered update, old and current vintages are different
    tempered_update = !isempty(old_data) # Time tempering
    @assert !(tempered_update&(old_vintage==data_vintage)) "Old & current vintages the same!"

    # General
    i   = 1             # Index tracking the stage of the algorithm
    j   = 2             # Index tracking the fixed_schedule entry ϕ_prop
    ϕ_n = ϕ_prop = 0.   # Instantiate ϕ_n and ϕ_prop variables

    resampled_last_period = false # Ensures proper resetting of ESS_bar after resample
    #use_fixed_schedule = (tempering_target == 0.0)
    threshold          = threshold_ratio * n_parts

    fixed_para_inds = findall([ θ.fixed for θ in parameters])
    free_para_inds  = findall([!θ.fixed for θ in parameters])
    para_symbols    = [θ.key for θ in parameters]

    n_para          = length(parameters)
    n_free_para     = length(free_para_inds)
    @assert n_free_para > 0 "All model parameters are fixed!"

    # Initialization of Particle Array Cloud
    cloud = Cloud(n_para, n_parts)

    #################################################################################
    ### Initialize Algorithm: Draws from prior
    #################################################################################
    println(verbose, :low, "\n\n SMC " * (testing ? "testing " : "") * "starts ....\n\n")

    if tempered_update
        # If user does not input Cloud object themselves, looks for cloud in loadpath.
        cloud = cloud_isempty(old_cloud) ? load(loadpath, "cloud") : old_cloud

        initialize_cloud_settings!(cloud; tempered_update = tempered_update,
                                   n_parts = n_parts, n_Φ = n_Φ, c = c, accept = target)
        initialize_likelihoods!(loglikelihood, parameters, data, cloud; parallel = parallel)

    elseif continue_intermediate
        cloud = load(loadpath, "cloud")
    else
        # Instantiating Cloud object, update draws, loglh, & logpost
        initial_draw!(loglikelihood, parameters, data, cloud; parallel = parallel)
        initialize_cloud_settings!(cloud; tempered_update = tempered_update,
                                   n_parts = n_parts, n_Φ = n_Φ, c = c, accept = target)
    end

    # Fixed schedule for construction of ϕ_prop
    if use_fixed_schedule
        cloud.tempering_schedule = ((collect(1:n_Φ) .- 1) / (n_Φ-1)) .^ λ
    else
        proposed_fixed_schedule  = ((collect(1:n_Φ) .- 1) / (n_Φ-1)) .^ λ
    end

    # Instantiate incremental and normalized weight matrices for logMDD calculation
    if continue_intermediate
        w_matrix = load(loadpath, "w")
        W_matrix = load(loadpath, "W")
        j        = load(loadpath, "j")
        i        = cloud.stage_index
        c        = cloud.c
        ϕ_prop   = proposed_fixed_schedule[j]
    else
        w_matrix = zeros(n_parts, 1)
        W_matrix = tempered_update ? get_weights(cloud) : fill(1/n_parts,(n_parts, 1))
    end

    # Printing
    init_stage_print(cloud, para_symbols; verbose = verbose,
                     use_fixed_schedule = use_fixed_schedule)
    println(verbose, :low, "\n\n SMC recursion starts... \n\n")

    #################################################################################
    ### Recursion
    #################################################################################
    while ϕ_n < 1.

        start_time = time_ns()
        cloud.stage_index = i += 1

        #############################################################################
        ### Setting ϕ_n (either adaptively or by the fixed schedule)
        #############################################################################
        ϕ_n1 = cloud.tempering_schedule[i-1]

        if use_fixed_schedule
            ϕ_n = cloud.tempering_schedule[i]
        else
            ϕ_n, resampled_last_period, j, ϕ_prop = solve_adaptive_ϕ(cloud,
                                                       proposed_fixed_schedule,
                                                       i, j, ϕ_prop, ϕ_n1,
                                                       tempering_target,
                                                       resampled_last_period)
        end

        #############################################################################
        ### Step 1: Correction
        #############################################################################

        # Calculate incremental weights (if no old data, get_old_loglh(cloud) = 0)
        incremental_weights = exp.((ϕ_n1 - ϕ_n) * get_old_loglh(cloud) +
                                   (ϕ_n - ϕ_n1) * get_loglh(cloud))

        # Update weights
        update_weights!(cloud, incremental_weights)
        mult_weights = get_weights(cloud)

        # Normalize weights
        normalize_weights!(cloud)
        normalized_weights = get_weights(cloud)

        w_matrix = hcat(w_matrix, incremental_weights)
        W_matrix = hcat(W_matrix, normalized_weights)

        ##############################################################################
        ### Step 2: Selection
        ##############################################################################

        # Calculate the degeneracy/effective sample size metric
        push!(cloud.ESS, 1 / sum(normalized_weights .^ 2))

        # If this assertion does not hold then there are likely too few particles
        @assert !isnan(cloud.ESS[i]) "No particles have non-zero weight."

        # Resample if degeneracy/ESS metric falls below the accepted threshold
        if (cloud.ESS[i] < threshold)

            # Resample according to particle weights, uniformly reset weights to 1/n_parts
            new_inds = resample(normalized_weights; method = resampling_method)
            cloud.particles = [deepcopy(cloud.particles[k,j]) for k in new_inds,
                               j=1:size(cloud.particles, 2)]
            reset_weights!(cloud)
            cloud.resamples += 1
            resampled_last_period = true
            W_matrix[:, i] .= 1/n_parts
        end

        ##############################################################################
        ### Step 3: Mutation
        ##############################################################################

        # Calculate adaptive c-step for use as scaling coefficient in mutation MH step
        c = c * (0.95 + 0.10 * exp(16.0 * (cloud.accept - target)) /
                 (1.0 + exp(16.0 * (cloud.accept - target))))
        cloud.c = c

        θ_bar = weighted_mean(cloud)
        R     = weighted_cov(cloud)

        # Ensures marix is positive semi-definite symmetric
        # (not off due to numerical error) and values haven't changed
        R_fr = (R[free_para_inds, free_para_inds] + R[free_para_inds, free_para_inds]') / 2.

        # MvNormal centered at ̄θ with var-cov ̄Σ, subsetting out the fixed parameters
        θ_bar_fr = θ_bar[free_para_inds]

        # Generate random parameter blocks
        blocks_free = generate_free_blocks(n_free_para, n_blocks)
        blocks_all  = generate_all_blocks(blocks_free, free_para_inds)

        new_particles = if parallel
            @distributed (hcat) for k in 1:n_parts
                mutation_closure(cloud.particles[k, :], θ_bar_fr, R_fr, blocks_free,
                                 blocks_all, ϕ_n, ϕ_n1; c = c, α = α,
                                 n_mh_steps = n_mh_steps, old_data = old_data)
            end
        else
            hcat([mutation_closure(cloud.particles[k, :], θ_bar_fr, R_fr,
                                   blocks_free, blocks_all, ϕ_n, ϕ_n1; c = c,
                                   α = α, n_mh_steps = n_mh_steps,
                                   old_data = old_data) for k=1:n_parts]...)
        end
        update_cloud!(cloud, new_particles)
        update_acceptance_rate!(cloud)

        ##############################################################################
        ### Timekeeping and Output Generation
        ##############################################################################
        cloud.total_sampling_time += Float64((time_ns() - start_time) * 1e-9)

        end_stage_print(cloud, para_symbols; verbose = verbose,
                        use_fixed_schedule = use_fixed_schedule)

        if run_test && (i == 3)
            break
        end

        if mod(cloud.stage_index, intermediate_stage_increment) == 0 && save_intermediate
            jldopen(replace(savepath, ".jld2" => "_stage=$(cloud.stage_index).jld2"),
                    true, true, true, IOStream) do file
                write(file, "cloud", cloud)
                write(file, "w", w_matrix)
                write(file, "W", W_matrix)
                write(file, "j", j)
            end
        end
    end

    ##################################################################################
    ### Saving data
    ##################################################################################
    if !testing
        simfile = h5open(particle_store_path, "w")
        particle_store = d_create(simfile, "smcparams", datatype(Float64),
                                  dataspace(n_parts, n_para))
        for k in 1:n_parts; particle_store[k,:] = cloud.particles[k, 1:n_para] end
        close(simfile)
        jldopen(savepath, true, true, true, IOStream) do file
            write(file, "cloud", cloud)
            write(file, "w", w_matrix)
            write(file, "W", W_matrix)
        end
    end
end


text:
using Oceananigans: znodes

"""
    animate_training_data(ds, ds_coarse; filepath, frameskip=1, fps=15)

Create an animation of the variable `v` and coarse-grained variable `v̄`.
`xlabel` should be specified for the plot. A `filepath should also be specified`.
A `.gif` and `.mp4` will be saved. `frameskip` > 1 can be used to skip some frames
to produce the animation more quickly. The output animation will use `fps` frames
per second.
"""
function animate_training_data(ds, ds_coarse; filepath, frameskip=1, fps=15)
    T = ds["T"]
    wT = ds["wT"]
    zc = znodes(T)
    zf = znodes(wT)
    times = T.times

    T_coarse = ds_coarse["T"]
    wT_coarse = ds_coarse["wT"]
    z̄c = znodes(T_coarse)
    z̄f = znodes(wT_coarse)

    kwargs = (linewidth=3, linealpha=0.8, ylims=extrema(zf),
              grid=false, legend=:bottomright, framestyle=:box,
              foreground_color_legend=nothing, background_color_legend=nothing)

    T_lims = interior(T) |> extrema
    wT_lims = interior(wT) |> extrema

    anim = @animate for n=1:frameskip:length(times)
        @info "Plotting free convection data [frame $n/$(length(times))]: $filepath..."

        time_str = @sprintf("%.2f days", times[n] / days)

        wT_profile = interior(wT)[1, 1, :, n]
        wT_coarse_profile = interior(wT_coarse)[1, 1, :, n]

        wT_plot = plot()
        plot!(wT_plot, wT_profile, zf, label="LES (Nz=$(length(zf)))", xlabel="Heat flux wT (m/s °C)",
              xlims=wT_lims, ylabel="Depth z (meters)", title="Free convection: $time_str"; kwargs...)
        plot!(wT_plot, wT_coarse_profile, z̄f, label="coarse (Nz=$(length(z̄f)))"; kwargs...)

        T_profile = interior(T)[1, 1, :, n]
        T_coarse_profile = interior(T_coarse)[1, 1, :, n]

        T_plot = plot()
        plot!(T_plot, T_profile, zc, label="LES (Nz=$(length(zc)))", xlabel="Temperature T (°C)",
              xlims=T_lims; kwargs...)
        plot!(T_plot, T_coarse_profile, z̄c, label="coarse (Nz=$(length(z̄c)))"; kwargs...)

        plot(wT_plot, T_plot, dpi=200)
    end

    @info "Saving $filepath..."
    mp4(anim, filepath * ".mp4", fps=fps)
    gif(anim, filepath * ".gif", fps=fps)

    return nothing
end

function animate_learned_free_convection(ds, NN, NN_function, NDEType, algorithm, T_scaling, wT_scaling; filepath, frameskip=1, fps=15)
    T = ds["T"]
    wT = ds["wT"]
    Nz = size(T, 3)
    Nt = size(T, 4)
    zc = znodes(T)
    zf = znodes(wT)
    times = T.times

    nde_params = FreeConvectionNDEParameters(ds, T_scaling, wT_scaling)
    T₀ = T_scaling.(interior(T)[1, 1, :, 1])
    nde = NDEType(NN, ds)
    nde_sol = solve_nde(nde, NN, T₀, algorithm, nde_params)

    kwargs = (linewidth=3, linealpha=0.8, ylims=extrema(zf),
              grid=false, legend=:bottomright, framestyle=:box,
              foreground_color_legend=nothing, background_color_legend=nothing)

    T_lims = interior(T) |> extrema
    wT_lims = interior(wT) |> extrema

    anim = @animate for n=1:frameskip:Nt
        @info "Plotting learned free convection [frame $n/$Nt]: $filepath ..."

        time_str = @sprintf("%.2f days", times[n] / days)

        wT_plot = plot()
        T_plot = plot()

        T_profile = interior(T)[1, 1, :, n]
        wT_profile = interior(wT)[1, 1, :, n]

        plot!(wT_plot, wT_profile, zf, label="Oceananigans", xlabel="Heat flux w'T' (m/s K)",
              xlims=wT_lims, ylabel="Depth z (meters)", title="Free convection: $time_str"; kwargs...)

        plot!(T_plot, T_profile, zc[:], label="Oceananigans", xlabel="Temperature T (°C)",
              xlims=T_lims; kwargs...)

        T_NN = nde_sol[n]
        bottom_flux = wT_scaling(interior(wT)[1, 1, 1, n])
        top_flux = wT_scaling(interior(wT)[1, 1, Nz+1, n])
        input = FreeConvectionTrainingDataInput(T_NN, bottom_flux, top_flux)
        wT_NN = unscale.(NN_function(input), Ref(wT_scaling))
        T_NN = unscale.(T_NN, Ref(T_scaling))

        plot!(wT_plot, wT_NN, zf, label="neural network"; kwargs...)
        plot!(T_plot, T_NN, zc, label="NDE"; kwargs...)

        plot(wT_plot, T_plot, dpi=200)
    end

    @info "Saving $filepath..."
    mp4(anim, filepath * ".mp4", fps=fps)
    gif(anim, filepath * ".gif", fps=fps)

    return nothing
end