text:
abstract type AbstractDiagTriangularMatrix{P,T,L} <: AbstractMatrix{T} end

abstract type AbstractTriangularMatrix{P,T,L} <: AbstractDiagTriangularMatrix{P,T,L} end
abstract type AbstractLowerTriangularMatrix{P,T,L} <: AbstractTriangularMatrix{P,T,L} end
abstract type AbstractUpperTriangularMatrix{P,T,L} <: AbstractTriangularMatrix{P,T,L} end
abstract type AbstractMutableLowerTriangularMatrix{P,T,L} <: AbstractLowerTriangularMatrix{P,T,L} end
abstract type AbstractMutableUpperTriangularMatrix{P,T,L} <: AbstractUpperTriangularMatrix{P,T,L} end

struct LowerTriangularMatrix{P,T,L} <: AbstractLowerTriangularMatrix{P,T,L}
    data::NTuple{L,T}
end
mutable struct MutableLowerTriangularMatrix{P,T,L} <: AbstractMutableLowerTriangularMatrix{P,T,L}
    data::NTuple{L,T}
    MutableLowerTriangularMatrix{P,T,L}(::UndefInitializer) where {P,T,L} = new{P,T,L}()
end
@inline LowerTriangularMatrix(M::MutableLowerTriangularMatrix{P,T,L}) where {P,T,L} = LowerTriangularMatrix{P,T,L}(M.data)
struct UpperTriangularMatrix{P,T,L} <: AbstractUpperTriangularMatrix{P,T,L}
    data::NTuple{L,T}
end
mutable struct MutableUpperTriangularMatrix{P,T,L} <: AbstractMutableUpperTriangularMatrix{P,T,L}
    data::NTuple{L,T}
    MutableUpperTriangularMatrix{P,T,L}(::UndefInitializer) where {P,T,L} = new{P,T,L}()
end
struct PtrLowerTriangularMatrix{P,T,L} <: AbstractMutableLowerTriangularMatrix{P,T,L}
    ptr::Ptr{T}
end
struct PtrUpperTriangularMatrix{P,T,L} <: AbstractMutableUpperTriangularMatrix{P,T,L}
    ptr::Ptr{T}
end

struct UninitializedLowerTriangularMatrix{P,T,L} <: AbstractMutableLowerTriangularMatrix{P,T,L}
    ptr::Ptr{T}
end
struct UninitializedUpperTriangularMatrix{P,T,L} <: AbstractMutableUpperTriangularMatrix{P,T,L}
    ptr::Ptr{T}
end
@inline function ReverseDiffExpressionsBase.uninitialized(L::AbstractMutableLowerTriangularMatrix{P,T,N}) where {P,T,N}
    UninitializedLowerTriangularMatrix{P,T,N}(pointer(L))
end
@inline function ReverseDiffExpressionsBase.uninitialized(U::AbstractMutableUpperTriangularMatrix{P,T,N}) where {P,T,N}
    UninitializedUpperTriangularMatrix{P,T,N}(pointer(U))
end
@inline Base.pointer(L::UninitializedLowerTriangularMatrix) = L.ptr
@inline Base.pointer(U::UninitializedUpperTriangularMatrix) = U.ptr
ReverseDiffExpressionsBase.isinitialized(::Type{<:UninitializedLowerTriangularMatrix}) = false
ReverseDiffExpressionsBase.isinitialized(::Type{<:UninitializedUpperTriangularMatrix}) = false
function PtrLowerTriangularMatrix{P,T,L}(sp::StackPointer) where {P,T,L}
    sp + VectorizationBase.align(sizeof(T)*L), PtrLowerTriangularMatrix{P,T,L}(pointer(sp, T))
end

@generated function PtrLowerTriangularMatrix{P,T}(sp::StackPointer) where {P,T}
    L = PaddedMatrices.calc_padding(binomial2(P+1),T)
    quote
        $(Expr(:meta,:inline))
        A = PtrLowerTriangularMatrix{$P,$T,$L}(pointer(sp, $T))
        sp + $(VectorizationBase.align(L*sizeof(T))), A
    end
end
@generated function PtrUpperTriangularMatrix{P,T}(sp::StackPointer) where {P,T}
    L = PaddedMatrices.calc_padding(binomial2(P+1),T)
    quote
        $(Expr(:meta,:inline))
        A = PtrUpperTriangularMatrix{$P,$T,$L}(pointer(sp, $T))
        sp + $(VectorizationBase.align(L*sizeof(T))), A
    end
end
@generated function ReverseDiffExpressionsBase.alloc_adjoint(
    L::AbstractMutableLowerTriangularMatrix{P,T}
) where {P,T}
    N = PaddedMatrices.calc_padding(binomial2(N+1), T)
    quote
        $(Expr(:meta,:inline))
        MutableLowerTriangularMatrix{$P,$T,$N}(undef)
    end
end
@generated function ReverseDiffExpressionsBase.alloc_adjoint(
    L::AbstractMutableUpperTriangularMatrix{P,T}
) where {P,T}
    N = PaddedMatrices.calc_padding(binomial2(N+1), T)
    quote
        $(Expr(:meta,:inline))
        MutableUpperTriangularMatrix{$P,$T,$N}(undef)
    end
end
@generated function ReverseDiffExpressionsBase.alloc_adjoint(
    ptr::Ptr{T}, L::AbstractMutableLowerTriangularMatrix{P,T}
) where {P,T}
    N = binomial2(P+1)
    quote
        $(Expr(:meta,:inline))
        PtrLowerTriangularMatrix{$P,$T,$N}(ptr)
    end
end
@generated function ReverseDiffExpressionsBase.alloc_adjoint(
    ptr::Ptr{T}, L::AbstractMutableUpperTriangularMatrix{P,T}
) where {P,T}
    N = binomial2(P+1)
    quote
        $(Expr(:meta,:inline))
        PtrUpperTriangularMatrix{$P,$T,$N}(ptr)
    end
end
@generated function ReverseDiffExpressionsBase.alloc_adjoint(
    ptr::StackPointer, L::AbstractMutableLowerTriangularMatrix{P,T}
) where {P,T}
    N = PaddedMatrices.calc_padding(binomial2(P+1), T)
    quote
        $(Expr(:meta,:inline))
        PtrLowerTriangularMatrix{$P,$T,$N}(ptr)
    end
end
@generated function ReverseDiffExpressionsBase.alloc_adjoint(
    ptr::StackPointer, L::AbstractMutableUpperTriangularMatrix{P,T}
) where {P,T}
    N = PaddedMatrices.calc_padding(binomial2(P+1), T)
    quote
        $(Expr(:meta,:inline))
        PtrUpperTriangularMatrix{$P,$T,$N}(ptr)
    end
end

@inline lt_sub2ind_fast(P, i, j) = i == j ? i : j * P - binomial2(j) + i - j
@inline ut_sub2ind_fast(P, i, j) = i == j ? i : 1 + P + binomial2(j) + i - j
function lt_sub2ind(P, i, j)
    i == j && return i
    j, i = minmax(i, j)
    j * P - binomial2(j) + i - j
end
function ut_sub2ind(P, i, j)
    i == j && return i
    i, j = minmax(i, j)
    1 + P + binomial2(j) + i - j
end

@inline UpperTriangularMatrix(M::MutableUpperTriangularMatrix{P,T,L}) where {P,T,L} = UpperTriangularMatrix{P,T,L}(M.data)
@generated function MutableLowerTriangularMatrix{P,T}(undef) where {P,T}
    Lbase = binomial2(P+1)
    Wm1 = VectorizationBase.pick_vector_width(Lbase,T) - 1
    L = (Lbase + Wm1) & ~Wm1
    quote
        $(Expr(:meta,:inline))
        MutableLowerTriangularMatrix{$P,$T,$L}(undef)
    end
end
@generated function MutableUpperTriangularMatrix{P,T}(undef) where {P,T}
    Lbase = binomial2(P+1)
    Wm1 = VectorizationBase.pick_vector_width(Lbase,T) - 1
    L = (Lbase + Wm1) & ~Wm1
    quote
        $(Expr(:meta,:inline))
        MutableUpperTriangularMatrix{$P,$T,$L}(undef)
    end
end
function MutableLowerTriangularMatrix(
    B::Union{A,LinearAlgebra.Adjoint{T,A}}
) where {T,M,A <: AbstractFixedSizeMatrix{M,M,T}}
    L = MutableLowerTriangularMatrix{M,T}(undef)
    @inbounds for mc ∈ 1:M, mr ∈ mc:M
        L[mr,mc] = B[mr,mc]
    end
    L
end
function MutableUpperTriangularMatrix(
    B::Union{A,LinearAlgebra.Adjoint{T,A}}
) where {T,M,A <: AbstractFixedSizeMatrix{M,M,T}}
    U = MutableLowerTriangularMatrix{M,T}(undef)
    @inbounds for mc ∈ 1:M, mr ∈ 1:mc
        U[mr,mc] = B[mr,mc]
    end
    U
end

function copyto!(L::MutableLowerTriangularMatrix{M,T}, A::Matrix) where {M,T}
    for m ∈ 1:M
        L[m] = A[m,m]
    end
    for mc ∈ 1:M-1
        for mr ∈ mc+1:M
            L[mr,mc] = A[mr,mc]
        end
    end
end
function MutableLowerTriangularMatrix(A::Matrix{T}) where {T}
    M, N = size(A)
    @assert M == N
    L = MutableLowerTriangularMatrix{M,T}(undef)
    copyto!(L, A)
    L
end

function copyto!(U::MutableUpperTriangularMatrix{M,T}, A::Matrix) where {M,T}
    for m ∈ 1:M
        U[m] = A[m,m]
    end
    for mc ∈ 2:M
        for mr ∈ 1:mc-1
            U[mr,mc] = A[mr,mc]
        end
    end
end
function MutableUpperTriangularMatrix(A::Matrix{T}) where {T}
    M, N = size(A)
    @assert M == N
    U = MutableUpperTriangularMatrix{M,T}(undef)
    copyto!(U, A)
    U
end

@inline Base.similar(A::AT) where {AT <: Union{<:AbstractMutableLowerTriangularMatrix,<:AbstractMutableUpperTriangularMatrix}} = AT(undef)
@inline function Base.similar(A::PtrLowerTriangularMatrix{P,T,L}) where {P,T,L}
    MutableLowerTriangularMatrix{P,T,L}(undef)
end
@inline function Base.similar(sp::StackPointer, A::AbstractLowerTriangularMatrix{P,T,L}) where {P,T,L}
    PtrLowerTriangularMatrix{P,T,L}(sp)
end
@inline function Base.similar(A::PtrUpperTriangularMatrix{P,T,L}) where {P,T,L}
    MutableUpperTriangularMatrix{P,T,L}(undef)
end
@inline function Base.similar(sp::StackPointer, A::AbstractUpperTriangularMatrix{P,T,L}) where {P,T,L}
    PtrUpperTriangularMatrix{P,T,L}(sp)
end

abstract type AbstractSymmetricMatrix{P,T,L} <: AbstractDiagTriangularMatrix{P,T,L} end
abstract type AbstractSymmetricMatrixL{P,T,L} <: AbstractSymmetricMatrix{P,T,L} end
abstract type AbstractSymmetricMatrixU{P,T,L} <: AbstractSymmetricMatrix{P,T,L} end
abstract type AbstractMutableSymmetricMatrixL{P,T,L} <: AbstractSymmetricMatrixL{P,T,L} end
abstract type AbstractMutableSymmetricMatrixU{P,T,L} <: AbstractSymmetricMatrixU{P,T,L} end

struct SymmetricMatrixL{P,T,L} <: AbstractSymmetricMatrixL{P,T,L}
    data::NTuple{L,T}
end
mutable struct MutableSymmetricMatrixL{P,T,L} <: AbstractMutableSymmetricMatrixL{P,T,L}
    data::NTuple{L,T}
    MutableSymmetricMatrixL{P,T,L}(::UndefInitializer) where {P,T,L} = new{P,T,L}()
end
function MutableSymmetricMatrixL(S::SymmetricMatrixL{M,T,L}) where {M,T,L}
    Sm = MutableSymmetricMatrixL{M,T,L}(undef)
    Sm.data = S.data
    Sm
end

struct PtrSymmetricMatrixL{P,T,L} <: AbstractMutableSymmetricMatrixL{P,T,L}
    ptr::Ptr{T}
end
# @generated function SymmetricMatrixL(S::PaddedMatrices.AbstractFixedSizeMatrix{P,P,T,R}) where {P,T,R}
#     q = quote end
#     qa = q.args
#     PaddedMatrices.load_L_quote!(qa, P, R, :Σ, :Σ)
#     L = binomial2(P+1)
#     Wm1 = VectorizationBase.pick_vector_width(L, T)-1
#     L = (L + Wm1) & ~Wm1
#     lq = store_packed_L_quote!(qa, P, :Σ, T, L)
#     quote
#         # $(Expr(:meta,:inline))
#         @inbounds begin
#             # begin
#             $q
#             $lq
#         end
#     end
# end

struct SymmetricMatrixU{P,T,L} <: AbstractSymmetricMatrix{P,T,L}
    data::NTuple{L,T}
end
struct PtrSymmetricMatrixU{P,T,L} <: AbstractMutableSymmetricMatrixU{P,T,L}
    ptr::Ptr{T}
end
mutable struct MutableSymmetricMatrixU{P,T,L} <: AbstractMutableSymmetricMatrixU{P,T,L}
    data::NTuple{L,T}
    MutableSymmetricMatrixU{P,T,L}(::UndefInitializer) where {P,T,L} = new{P,T,L}()
end
@generated function MutableSymmetricMatrixL{P,T}(undef) where {P,T}
    Lbase = binomial2(P+1)
    W = VectorizationBase.pick_vector_width(Lbase,T)
    Wm1 = W - 1
    L = (Lbase + Wm1) & ~Wm1
    quote
        $(Expr(:meta,:inline))
        MutableSymmetricMatrixL{$P,$T,$L}(undef)
    end
end
@generated function MutableSymmetricMatrixU{P,T}(undef) where {P,T}
    Lbase = binomial2(P+1)
    W = VectorizationBase.pick_vector_width(Lbase,T)
    Wm1 = W - 1
    L = (Lbase + Wm1) & ~Wm1
    quote
        $(Expr(:meta,:inline))
        MutableSymmetricMatrixU{$P,$T,$L}(undef)
    end
end
function Base.copyto!(L::AbstractLowerTriangularMatrix{M,T}, A::AbstractMatrix{T}) where {M,T}
    @boundscheck begin
        N, P = size(A)
        @assert M == N == P
    end
    @inbounds for m in 1:M
        L[m] = A[m,m]
    end
    ind = M
    for c in 1:M
        for r in c+1:M
            ind += 1
            L[ind] = A[r, c]
        end
    end
    L
end
function MutableLowerTriangularMatrix(A::AbstractMatrix{T}) where {T}
    @assert size(A,1) == size(A,2) "Matrix of size $(size(A)) is not square."
    M = size(A,1)
    L = MutableLowerTriangularMatrix{M,T}(undef)
    copyto!(L, A)
    L
end
function MutableLowerTriangularMatrix{M}(A::AbstractMatrix{T}) where {M,T}
    @assert M == size(A,1) == size(A,2) "Matrix of size $(size(A)) is not a square $M x $M matrix."
    L = MutableLowerTriangularMatrix{M,T}(undef)
    copyto!(L, A)
    L
end
function MutableLowerTriangularMatrix{M,T}(A::AbstractMatrix) where {M,T}
    @assert M == size(A,1) == size(A,2) "Matrix of size $(size(A)) is not a square $M x $M matrix."
    L = MutableLowerTriangularMatrix{M,T}(undef)
    copyto!(L, A)
    L
end

const AbstractMutableDiagMatrix{P,T,L} = Union{
    AbstractMutableLowerTriangularMatrix{P,T,L},
    AbstractMutableUpperTriangularMatrix{P,T,L},
    AbstractMutableSymmetricMatrixL{P,T,L},
    AbstractMutableSymmetricMatrixU{P,T,L}
}

@inline function Base.getindex(A::AbstractDiagTriangularMatrix{P,T,L}, i::Integer) where {P,T,L}
    @boundscheck i > L && ThrowBoundsError("i = $i > L = $L")
    @inbounds A.data[i]
end
@inline function Base.getindex(A::LinearAlgebra.Adjoint{T,<:AbstractDiagTriangularMatrix{P,T,L}}, i::Integer) where {P,T,L}
    @boundscheck i > L && ThrowBoundsError("i = $i > L = $L")
    @inbounds A.parent.data[i]
end
@inline function Base.getindex(A::LinearAlgebra.Adjoint{Union{},<:AbstractDiagTriangularMatrix{P,Vec{W,T},L}}, i::Integer) where {P,T,L,W}
    @boundscheck i > L && ThrowBoundsError("i = $i > L = $L")
    @inbounds A.parent.data[i]
end

const MutableDiagTriangle{P,T,L} = Union{
    MutableLowerTriangularMatrix{P,T,L},
    MutableUpperTriangularMatrix{P,T,L},
    MutableSymmetricMatrixL{P,T,L},
    MutableSymmetricMatrixU{P,T,L}
}

@inline Base.pointer(A::MutableDiagTriangle{P,T,L}) where {P,T,L} = Base.unsafe_convert(Ptr{T}, pointer_from_objref(A))
@inline function Base.pointer(A::MutableDiagTriangle{P,NTuple{W,Core.VecElement{T}},L}) where {P,T,L,W}
    Base.unsafe_convert(Ptr{T}, pointer_from_objref(A))
end
@inline Base.pointer(A::PtrLowerTriangularMatrix) = A.ptr
@inline Base.pointer(A::PtrUpperTriangularMatrix) = A.ptr
@inline Base.pointer(A::PtrSymmetricMatrixL) = A.ptr
@inline Base.pointer(A::PtrSymmetricMatrixU) = A.ptr

@inline Base.pointer(A::PtrLowerTriangularMatrix{P,Vec{W,T}}) where {P,W,T} = Base.unsafe_convert(Ptr{T},A.ptr)
@inline Base.pointer(A::PtrUpperTriangularMatrix{P,Vec{W,T}}) where {P,W,T} = Base.unsafe_convert(Ptr{T},A.ptr)
@inline Base.pointer(A::PtrSymmetricMatrixL{P,Vec{W,T}}) where {P,W,T} = Base.unsafe_convert(Ptr{T},A.ptr)
@inline Base.pointer(A::PtrSymmetricMatrixU{P,Vec{W,T}}) where {P,W,T} = Base.unsafe_convert(Ptr{T},A.ptr)

@inline function Base.getindex(A::AbstractMutableDiagMatrix{P,T,L}, i::Integer) where {P,T,L}
    @boundscheck i > L && ThrowBoundsError("i = $i > L = $L")
    VectorizationBase.load(pointer(A) + sizeof(T) * (i-1))
end
@inline function Base.getindex(A::LinearAlgebra.Adjoint{T,<:AbstractMutableDiagMatrix{P,T,L}}, i::Integer) where {P,T,L}
    @boundscheck i > L && ThrowBoundsError("i = $i > L = $L")
    VectorizationBase.load(pointer(A.parent) + sizeof(T) * (i-1))
end
@inline function Base.getindex(A::AbstractMutableDiagMatrix{P,Vec{W,T},L}, i::Integer) where {P,W,T,L}
    @boundscheck i > L && ThrowBoundsError("i = $i > L = $L")
    SIMDPirates.vload(Vec{W,T}, pointer(A) + W*sizeof(T)*(i-1))
end
@inline function Base.getindex(A::LinearAlgebra.Adjoint{Union{},<:AbstractMutableDiagMatrix{P,Vec{W,T},L}}, i::Integer) where {P,W,T,L}
    @boundscheck i > L && ThrowBoundsError("i = $i > L = $L")
    SIMDPirates.vload(Vec{W,T}, pointer(A.parent) + W*sizeof(T)*(i-1))
end

@inline function Base.setindex!(A::AbstractMutableDiagMatrix{P,T,L}, v, i::Integer) where {P,T,L}
    @boundscheck i > L && ThrowBoundsError("i > $L.")
    VectorizationBase.store!(pointer(A) + (i-1) * sizeof(T), convert(T,v))
    v
end
@inline function Base.setindex!(A::AbstractMutableDiagMatrix{P,NTuple{W,Core.VecElement{T}},L}, v::NTuple{W,Core.VecElement{T}}, i::Integer) where {P,T,L,W}
    @boundscheck i > L && ThrowBoundsError("i > $L.")
    SIMDPirates.vstore!(pointer(A) + (i-1) * sizeof(NTuple{W,Core.VecElement{T}}), v)
    v
end

@inline function Base.setindex!(A::AbstractMutableLowerTriangularMatrix{P,T,L}, v, i::Integer, j::Integer) where {P,T,L}
    ind = lt_sub2ind_fast(P, i, j)
    @boundscheck ind > L && ThrowBoundsError("i > $L.")
    VectorizationBase.store!(pointer(A) + (ind-1) * sizeof(T), convert(T,v))
    v
end
@inline function Base.setindex!(A::AbstractMutableLowerTriangularMatrix{P,NTuple{W,Core.VecElement{T}},L}, v::NTuple{W,Core.VecElement{T}}, i::Integer, j::Integer) where {P,T,L,W}
    ind = lt_sub2ind_fast(P, i, j)
    @boundscheck ind > L && ThrowBoundsError("i > $L.")
    SIMDPirates.vstore!(pointer(A) + (ind-1) * sizeof(NTuple{W,Core.VecElement{T}}), v)
    v
end
@inline function Base.setindex!(A::AbstractMutableUpperTriangularMatrix{P,T,L}, v, i::Integer, j::Integer) where {P,T,L}
    ind = ut_sub2ind_fast(P, i, j)
    @boundscheck ind > L && ThrowBoundsError("i > $L.")
    VectorizationBase.store!(pointer(A) + (ind-1) * sizeof(T), convert(T,v))
    v
end
@inline function Base.setindex!(A::AbstractMutableUpperTriangularMatrix{P,NTuple{W,Core.VecElement{T}},L}, v::NTuple{W,Core.VecElement{T}}, i::Integer, j::Integer) where {P,T,L,W}
    ind = ut_sub2ind_fast(P, i, j)
    @boundscheck ind > L && ThrowBoundsError("i > $L.")
    SIMDPirates.vstore!(pointer(A) + (ind-1) * sizeof(NTuple{W,Core.VecElement{T}}), v)
    v
end

#
# struct LowerTriangleColumnIterator{MT}
#     num_col::Int
#     W::Int
#     first_mask::MT
#     last_mask::MT
#     max_mask::MT
# end
# function LowerTriangleColumnIterator(::Val{P}, ::Type{T}) where {P,T}
#     W, Wshift = pick_vector_width_shift(P, T)
#     Wm1 = W - 1
#     MT = VectorizationBase.mask_type(Val(P), T)
#
#     max_mask = MT((1 << Wshift) - 1)
#     LowerTriangleColumnIterator(P, W, first_mask, last_mask, max_mask)
# end
# # Base.eltype(::Type{LowerTriangleColumnIterator}) =
# Base.length(i::LowerTriangleColumnIterator) = i.num_col - 1
#
# function Base.iterate(t::LowerTriangleColumnIterator)
#
# end
# function Base.iterate(t::LowerTriangleColumnIterator, state)
#
# end

@generated type_length(::Type{<:AbstractDiagTriangularMatrix{M}}) where {M} = binomial2(M+1)
@generated type_length(::AbstractDiagTriangularMatrix{P}) where {P} = binomial2(P+1)
@generated param_type_length(::Type{<:AbstractDiagTriangularMatrix{P}}) where {P} = binomial2(P+1)
@generated param_type_length(::AbstractDiagTriangularMatrix{P}) where {P} = binomial2(P+1)

Base.size(::AbstractDiagTriangularMatrix{P}) where {P} = (P,P)

@inline VectorizationBase.vectorizable(A::AbstractMutableDiagMatrix) = VectorizationBase.Pointer(pointer(A))
@inline VectorizationBase.vectorizable(A::AbstractDiagTriangularMatrix) = PaddedMatrices.vStaticPaddedArray(A)

# @inline binom2(n::Int) = (nu = reinterpret(UInt, n); reinterpret(Int, (nu*(nu-one(UInt))) >> one(UInt)))

# function padded_diagonal_length(P, T)
#     Wm1 = VectorizationBase.pick_vector_width(P, T) - 1
#     (P + Wm1) & ~Wm1
# end
# function num_complete_blocks(P, T)
#     W, Wshift = VectorizationBase.pick_vector_width_shift(P, T)
#     (P - 1) >> Wshift
# end
# function calculate_L(P, T)
#     W, Wshift = VectorizationBase.pick_vector_width_shift(P, T)
#     W² = W*W
#     Wm1 = W - 1
#     rem = P & Wm1
#     padded_diag_length = (P + Wm1) & ~Wm1
#     L = padded_diag_length
#
#     Pm1 = P - 1
#     num_complete_blocks = Pm1 >> Wshift
#     L += W² * binomial(1+num_complete_blocks, 2)
#     rem_block = Pm1 & Wm1
#     L += (1+num_complete_blocks)*W*rem_block
#     L
# end
#
# function lower_triangle_sub2ind_quote(P, T)
#     # assume inputs are i, j
#     # for row, column
#     W, Wshift = VectorizationBase.pick_vector_width_shift(P, T)
#     W² = W * W
#     Wm1 = W - 1
#     rem = P & Wm1
#     if rem == 0
#         pad = 0
#         diag_length = P
#         num_rows = P
#         rem_row = Wm1
#     else
#         pad = W - rem
#         diag_length = P + pad
#         if rem == 1
#             num_rows = P - rem
#             rem_row = W
#         else
#             num_rows = diag_length
#             rem_row = rem - 1
#         end
#         # num_rows = rem == 1 ? P - rem : diag_length
#         # rem_row = rem - 1
#     end
#     # @show diag_length, num_rows, rem_row
#     quote
#         # We want i >= j, for a lower triangular matrix.
#         j, i = minmax(j, i)
#         # @boundscheck i > $P && ThrowBoundsError(str)
#         i == j && return i
#         ind = $(P - num_rows + Wm1 - rem_row) + j*$num_rows + i
#         if j > $rem_row
#             j2 = j - $rem_row
#             jrem = j2 & $Wm1
#             number_blocks = 1+(j2 >> $Wshift)
#             ind -= jrem * $W * number_blocks + binomial2(number_blocks) << $(2Wshift)
#         end
#         ind
#     end
# end
# @generated function lower_triangle_sub2ind(::Val{P}, ::Type{T}, i, j) where {T,P}
#     quote
#         # $(Expr(:meta, :inline))
#         $(lower_triangle_sub2ind_quote(P, T))
#     end
# end
#
# function number_lower_triangular_vector_loads_for_column(P, T, j)
#     W, Wshift = VectorizationBase.pick_vector_width_shift(P, T)
# end
#
# function upper_triangle_sub2ind_quote(P, T)
#     # assume inputs are i, j
#     # for row, column
#     W, Wshift = VectorizationBase.pick_vector_width_shift(P, T)
#     Wm1 = W - 1
#     rem = P & Wm1
#     diag_length = (P + Wm1) & ~Wm1
#     quote
#         # We want i <= j, for a upper triangular matrix.
#         i, j = minmax(i, j)
#         i == j && return i
#         num_blocks = (j-2) >> $Wshift
#         ind = $diag_length + ((num_blocks*(num_blocks+1)) << $(2Wshift - 1))
#         jrem = (j-2) & $(W-1)
#         ind += ((num_blocks+1) * jrem) << $Wshift
#         ind + i
#     end
# end
# @generated function upper_triangle_sub2ind(::Val{P}, ::Type{T}, i, j) where {T,P}
#     quote
#         # $(Expr(:meta, :inline))
#         $(upper_triangle_sub2ind_quote(P, T))
#     end
# end

@generated function Base.:+(A::AbstractLowerTriangularMatrix{M,T,L}, B::AbstractLowerTriangularMatrix{M,T,L}) where {M,T,L}
    quote
        $(Expr(:meta,:inline))
        C = MutableLowerTriangularMatrix{$M,$T,$L}(undef)
        @vectorize $T for l ∈ 1:$L
            C[l] = A[l] + B[l]
        end
        C
    end
end
@generated function Base.:+(A::AbstractUpperTriangularMatrix{M,T,L}, B::AbstractUpperTriangularMatrix{M,T,L}) where {M,T,L}
    quote
        $(Expr(:meta,:inline))
        C = MutableUpperTriangularMatrix{$M,$T,$L}(undef)
        @vectorize $T for l ∈ 1:$L
            C[l] = A[l] + B[l]
        end
        C
    end
end
@generated function Base.:+(A::UpperTriangularMatrix{M,T,L}, B::UpperTriangularMatrix{M,T,L}) where {M,T,L}
    quote
        $(Expr(:meta,:inline))
        C = MutableUpperTriangularMatrix{$M,$T,$L}(undef)
        @vectorize $T for l ∈ 1:$L
            C[l] = A[l] + B[l]
        end
        UpperTriangularMatrix(C)
    end
end

@inline SIMDPirates.lifetime_start!(A::MutableLowerTriangularMatrix) = nothing
@inline SIMDPirates.lifetime_end!(A::MutableLowerTriangularMatrix) = nothing
@inline function SIMDPirates.lifetime_start!(A::AbstractMutableLowerTriangularMatrix{M,T,L}) where {M,T,L}
    SIMDPirates.lifetime_start!(pointer(A), Val{L}())
end
@inline function SIMDPirates.lifetime_end!(A::AbstractMutableLowerTriangularMatrix{M,T,L}) where {M,T,L}
    SIMDPirates.lifetime_end!(pointer(A), Val{L}())
end


text:
abstract type AbstractCosGaussian end

"""
    CosGaussianParams( dims, k, α, σ, μ, δ )

Parameters of a distribution with is a product of a Cosine
distribution along x and a Normal distribution along v.

- `n_gaussians` : Number of Gaussians
- `n_cos`       : Number of cosines
- `normal`      : Normalization constant of each Gaussian
"""
struct CosGaussianParams
    dims::Tuple{Int64,Int64}
    n_cos::Int64
    n_gaussians::Int64
    k::Array{Vector{Float64},1}
    α::Vector{Float64}
    σ::Array{Vector{Float64},1}
    μ::Array{Vector{Float64},1}
    normal::Vector{Float64}
    δ::Vector{Float64}

    function CosGaussianParams(
        dims::Tuple{Int64,Int64},
        k::Array{Vector{Float64},1},
        α::Vector{Float64},
        σ::Array{Vector{Float64},1},
        μ::Array{Vector{Float64},1},
        δ::Vector{Float64}=[1.0],
    )
        n_cos = length(k)
        @assert n_cos == length(α)
        for i in 1:n_cos
            @assert length(k[i]) == dims[1]
        end
        n_gaussians = length(σ)
        @assert n_gaussians == length(μ)
        @assert all(σ .!= 0.0)
        normal = zeros(n_gaussians)
        for j in 1:n_gaussians
            @assert length(σ[j]) == dims[2]
            normal[j] = 1.0 / ((2π)^(0.5 * dims[2]) * prod(σ[j]))
        end
        @assert sum(δ) == 1.0

        return new(dims, n_cos, n_gaussians, k, α, σ, μ, normal, δ)
    end
end

export CosSumGaussian

"""
    CosSumGaussian{D,V}( n_cos, n_gaussians, k, α, σ, μ, δ )

Data type for parameters of initial distribution

```math
(1+ \\cos( \\sum^{n_{cos}}_{i=1} k_i x))
\\cdot
\\sum_{j=1}^{n_{gaussians}}
\\delta_j
\\exp \\big( -\\frac{1}{2}
\\frac{(v-\\mu_j)^2}{\\sigma_j^2} \\big)
```

## Parameters

- `k` : values of the wave numbers (one array for each cosines)
- `α` : strength of perturbations
- `σ` : variance of the Gaussian (one velocity vector for each gaussian).
- `μ` : mean value of the Gaussian (one velocity vector for each gaussian).
- `δ` : portion of each Gaussian

# Example

```math
f(x,v_1,v_2)=\\frac{1}{2\\pi\\sigma_1\\sigma_2} \\exp \\Big( - \\frac{1}{2}
\\big( \\frac{v_1^2}{\\sigma_1^2} + \\frac{v_2^2}{\\sigma_2^2} \\big)
\\Big) ( 1 + \\alpha \\cos(kx)),
```

```julia
df = CosSumGaussian{1,2}([[k]],[α], [[σ₁,σ₂]], [[μ₁,μ₂]])
```

"""
struct CosSumGaussian{D,V} <: AbstractCosGaussian
    dims::Tuple{Int64,Int64}
    params::CosGaussianParams

    function CosSumGaussian{D,V}(
        k::Array{Vector{Float64},1},
        α::Vector{Float64},
        σ::Array{Vector{Float64},1},
        μ::Array{Vector{Float64},1},
        δ::Vector{Float64}=[1.0],
    ) where {D,V}
        dims = (D, V)
        params = CosGaussianParams(dims, k, α, σ, μ, δ)

        return new(dims, params)
    end
end

export SumCosGaussian

"""
    SumCosGaussian( dims, n_cos, n_gaussians, k, α, σ, μ, δ )
Data type for parameters of initial distribution
```math
(1+ \\sum_{i=1}^{n_{cos}} \\alpha_i \\cos(  k_i \\mathbf{x}))
\\cdot
\\sum_{j=1}^{n_{gaussians}}
\\delta_j \\exp
\\big( -\\frac{1}{2}
\\frac{(\\mathbf{v}-\\mu_j)^2}{\\sigma_j^2} \\big)
```
## Parameters
- `k` : values of the wave numbers (Array of vectors for multiple cosines)
- `α` : strength of perturbations
- `σ` : variance of the Gaussian ( Array of vectors for multiple Gaussians)
- `μ` : mean value of the Gaussian ( Array multiple Gaussians)
- `normal` : Normalization constant of each Gaussian
- `n_gaussians` : Number of Gaussians
- `n_cos` : Number of cosines
- `δ` : portion of each Gaussian

# Example

```math
f(x,v_1,v_2) = \\frac{1}{2\\pi\\sigma_1\\sigma_2}
\\exp \\Big( - \\frac{1}{2} \\big( \\frac{v_1^2}{\\sigma_1^2}
 + \\frac{v_2^2}{\\sigma_2^2} \\big) \\Big)
( 1 + \\alpha_1 \\cos(k_1 x) + \\alpha_2 \\cos(k_2 x) ),
```
```julia
df = SumCosGaussian{1,2}([[k₁],[k₂]], [α₁, α₂], [[σ₁,σ₂]], [[0.0,0.0]])

```
"""
struct SumCosGaussian{D,V} <: AbstractCosGaussian
    dims::Tuple{Int64,Int64}
    params::CosGaussianParams

    function SumCosGaussian{D,V}(
        k::Array{Array{Float64,1},1},
        α::Vector{Float64},
        σ::Array{Array{Float64,1},1},
        μ::Array{Array{Float64,1},1},
        δ::Array{Float64,1}=[1.0],
    ) where {D,V}
        dims = (D, V)
        params = CosGaussianParams(dims, k, α, σ, μ, δ)
        return new(dims, params)
    end
end

function eval_x_density(f::CosSumGaussian, x::Union{Float64,Vector{Float64}})
    fval = 1.0
    for j in 1:(f.params.n_cos)
        fval += f.params.α[j] * cos(sum(f.params.k[j] .* x))
    end
    return fval
end

"""
    eval_x_density( f, x )

evaluate the cosine part of the distribution function
"""
function eval_x_density(f::SumCosGaussian, x::Union{Float64,Vector{Float64}})
    fval = 1.0
    for j in 1:(f.params.n_cos)
        fval += f.params.α[j] * cos(sum(f.params.k[j] .* x))
    end
    return fval
end

"""
    eval_v_density( f, v )

evaluate the normal part of the distribution function
"""
function eval_v_density(f::AbstractCosGaussian, v::Union{Float64,Vector{Float64}})
    fval = 0.0
    for j in 1:(f.params.n_gaussians)
        fval +=
            f.params.normal[j] * f.params.δ[j] .*
            exp(-0.5 * sum(((v .- f.params.μ[j]) ./ f.params.σ[j]) .^ 2))
    end
    return fval
end

function (f::CosSumGaussian)(x, v)
    return eval_x_density(f, x) * eval_v_density(f, v)
end

function (f::SumCosGaussian)(x, v)
    return eval_x_density(f, x) * eval_v_density(f, v)
end

ext: jl


text:
module SpecialExponentialFamilies

const SEF = SpecialExponentialFamilies

using Reexport

@reexport using Distributions
using StatsBase
import StatsBase: confint, fit

using DiffResults
using Empirikos
using Expectations

using ForwardDiff
using GLM
using LinearAlgebra

using OffsetArrays
using Optim
using LinearAlgebra
using Statistics

using Random
using Requires

import Base: eltype
import Base.Broadcast: broadcastable
import Base: minimum, maximum, extrema
import Distributions: probs, logpdf, pdf, cf, mgf, insupport, support, sampler
import StatsBase: fit, moment
import Statistics: mean, var
import Random: rand, rand!, AbstractRNG

using UnPack

include("splines.jl")
include("continuous_exp_families.jl")
include("lindsey.jl")
include("logspline_deconvolution.jl")
include("datasets.jl")

function __init__()
    @require ApproxFun = "28f2ccd6-bb30-5033-b560-165f7b14dc2f" include(
        "sample_continuous_exp_families.jl"
    )
end

export ExponentialFamily, ExponentialFamilyDistribution, LindseyMethod, SEF

export fit, confint

end # module


text:
import Base: convert,getindex,setindex!,+,-,*,/,==,length,norm

# Vector3
struct Vec3{T}
  x::T
  y::T
  z::T
end
Vec3(x,y)=Vec3(x,y,y)
Vec3(x)=Vec3(x,x,x)
Vec3()=Vec3(0.0)
function getindex(V::Vec3,i::Int)
  if i==1
    X=V.x
  elseif i==2
    X=V.y
  elseif i==3
    X=V.z
  else
    error("Out of Bound of Vector3.")
  end
end
function setindex!(V::Vec3,X,i::Int)
  if i==1
    V.x=X
  elseif i==2
    V.y=X
  elseif i==3
    V.z=X
  else
    error("Out of Bound of Vector3.")
  end
end
function convert(::Type{Vec3},a::AbstractVector)
  l=length(a)
  if l==0
    error("Empty Vector.")
  elseif l==1
    v = Vec3(a[1])
  elseif l==2
    v = Vec3(a[1],a[2])
  else
    v = Vec3(a[1],a[2],a[3])
  end
end
function convert(::Type{Vector},v::Vec3)
  a = [v.x,v.y,v.z]
end
function convert{T<:Vec3}(::Type{Matrix},vs::AbstractVector{T})
  l = length(vs)
  if l==0
    error("Empty Vector.")
  else
    a = Array(typeof(vs[1].x),3,l)
    for i in 1:l
      a[1,i]=vs[i].x
      a[2,i]=vs[i].y
      a[3,i]=vs[i].z
    end
  end
  return a
end
function convert{T<:Vec3}(::Type{Vector{T}},a::Matrix)
  s1,s2 = size(a)
  if s1 != 4
    error("Size of Dim 1 of Matrix Doesn't Match Vec4")
  elseif s2 == 0
    return []
  else
    vs = Array(T,s2)
    for i in 1:s2
      vs[i] = Vec3(a[1,i],a[2,i],a[3,i])
    end
    return vs
  end
end
+(a::Vec3,b::Vec3)=Vec3(a.x+b.x,a.y+b.y,a.z+b.z)
-(a::Vec3,b::Vec3)=Vec3(a.x-b.x,a.y-b.y,a.z-b.z)
*(a::Vec3,b::Real)=Vec3(a.x*b,a.y*b,a.z*b)
*(b::Real,a::Vec3)=*(a,b)
/(a::Vec3,b::Real)=Vec3(a.x/b,a.y/b,a.z/b)
/(b::Real,a::Vec3)=/(a,b)
==(a::Vec3,b::Vec3)=(a.x==b.x) && (a.y==b.y) && (a.z==b.z)

# Vector4
struct Vec4{T}
  x::T
  y::T
  z::T
  w::T
end
Vec4(x,y,z)=Vec4(x,y,z,z)
Vec4(x,y)=Vec4(x,y,y,y)
Vec4(x)=Vec4(x,x,x,x)
Vec4()=Vec4(0.0,0.0,0.0,1.0)
function getindex(V::Vec4,i::Int)
  if i==1
    X=V.x
  elseif i==2
    X=V.y
  elseif i==3
    X=V.z
  elseif i==4
    X=V.w
  else
    error("Out of Bound of Vector4.")
  end
end
function setindex!(V::Vec4,X,i::Int)
  if i==1
    V.x=X
  elseif i==2
    V.y=X
  elseif i==3
    V.z=X
  elseif i==4
    V.w=X
  else
    error("Out of Bound of Vector4.")
  end
end
function convert(::Type{Vec4},a::AbstractVector)
  l=length(a)
  if l==0
    error("Empty Vector.")
  elseif l==1
    v = Vec4(a[1])
  elseif l==2
    v = Vec4(a[1],a[2])
  elseif l==3
    v = Vec4(a[1],a[2],a[3])
  else
    v = Vec4(a[1],a[2],a[3],a[4])
  end
end
function convert(::Type{Vector},v::Vec4)
  a = [v.x,v.y,v.z,v.w]
end
function convert{T<:Vec4}(::Type{Matrix},vs::AbstractVector{T})
  l = length(vs)
  if l==0
    error("Empty Vector.")
  else
    a = Array(typeof(vs[1].x),4,l)
    for i in 1:l
      a[1,i]=vs[i].x
      a[2,i]=vs[i].y
      a[3,i]=vs[i].z
      a[4,i]=vs[i].w
    end
  end
  return a
end
function convert{T<:Vec4}(::Type{Vector{T}},a::Matrix)
  s1,s2 = size(a)
  if s1 != 4
    error("Size of Dim 1 of Matrix Doesn't Match Vec4")
  elseif s2 == 0
    return []
  else
    vs = Array(T,s2)
    for i in 1:s2
      vs[i] = Vec4(a[1,i],a[2,i],a[3,i],a[4,i])
    end
    return vs
  end
end
convert(::Type{Vec3},v::Vec4)=Vec3(v.x,v.y,v.z)
convert(::Type{Vec4},v::Vec3)=Vec4(v.x,v.y,v.z,1.0)
+(a::Vec4,b::Vec4)=Vec4(a.x+b.x,a.y+b.y,a.z+b.z,a.w+b.w)
-(a::Vec4,b::Vec4)=Vec4(a.x-b.x,a.y-b.y,a.z-b.z,a.w-b.w)
*(a::Vec4,b::Real)=Vec4(a.x*b,a.y*b,a.z*b,a.w)
*(b::Real,a::Vec4)=*(a,b)
/(a::Vec4,b::Real)=Vec4(a.x/b,a.y/b,a.z/b,a.w)
/(b::Real,a::Vec4)=/(a,b)
==(a::Vec4,b::Vec4)=(a.x==b.x) && (a.y==b.y) && (a.z==b.z) && (a.w==b.w)

function length(V::Union(Vec3,Vec4))
  sqrt(V.x^2+V.y^2+V.z^2)
end
function norm(V::Union(Vec3,Vec4))
  V/length(V)
end

# Transformation
function translation(x,y,z)
  tm =  [1.0 0.0 0.0   x;
         0.0 1.0 0.0   y;
         0.0 0.0 1.0   z;
         0.0 0.0 0.0 1.0]
end
translation(t::Union(Vec3,Vec4))=translation(t.x,t.y,t.z)
translation(;x=0.0,y=0.0,z=0.0)=translation(x,y,z)

function rotationxyz(theta::Real,axis::Int)
  if (3<axis) || (axis <1)
    error("Only X:1, Y:2, Z:3 Axis Allowed.")
  else
    r=Vec3()
    r[axis]=theta
    rotationxyz(r)
  end
end
function rotationxyz(r::Vec3)
  ar = convert(Vector,r)
  ri = find(ar)
  rn = length(ri)
  if rn == 0
    return eye(4)
  elseif rn == 1
    axis = ri[1]
  else
    error("Only One Axis Allowed.")
  end

  c = cos(ar[axis])
  s = sin(ar[axis])

  if axis==1
    rm = [1.0 0.0 0.0 0.0;
          0.0   c  -s 0.0;
          0.0   s   c 0.0;
          0.0 0.0 0.0 1.0]
  elseif axis==2
    rm = [  c 0.0   s 0.0;
          0.0 1.0 0.0 0.0;
          -s 0.0   c 0.0;
          0.0 0.0 0.0 1.0]
  else
    rm = [  c  -s 0.0 0.0;
          s   c 0.0 0.0;
          0.0 0.0 1.0 0.0;
          0.0 0.0 0.0 1.0]
  end
end
rotationxyz(ar::AbstractVector)=rotationxyz(convert(Vec3,ar))

function reflection(r::Vec3)
  ar = convert(Vector,r)
  ri = find(ar)
  rn = length(ri)
  if rn != 2
    error("Only XY, XZ, YZ Plane Allowed.")
  end
  pi = sum(ri)
  if pi==3 # xy
    rm = [1.0  0.0  0.0  0.0;
          0.0  1.0  0.0  0.0;
          0.0  0.0 -1.0  0.0;
          0.0  0.0  0.0  1.0]
  elseif pi==4 # xz
    rm = [1.0  0.0  0.0  0.0;
          0.0 -1.0  0.0  0.0;
          0.0  0.0  1.0  0.0;
          0.0  0.0  0.0  1.0]
  else # yz
    rm = [-1.0  0.0  0.0  0.0;
          0.0  1.0  0.0  0.0;
          0.0  0.0  1.0  0.0;
          0.0  0.0  0.0  1.0]
  end
end
reflection(ar::AbstractVector)=reflection(convert(Vec3,ar))

function transform(v::Vec4,m::Matrix)
  s1,s2 = size(m)
  if (s1 != s2) || (s1 != 4)
    error("Invalid Transform Matrix.")
  end
  convert(Vec4,m*convert(Vector,v))
end
function transform{T<:Vec4}(vs::AbstractVector{T},m::Matrix)
  s1,s2 = size(m)
  if (s1 != s2) || (s1 != 4)
    error("Invalid Transform Matrix.")
  end
  convert(Vector{T},m*convert(Matrix,vs))
end
translate(v::Vec4,t::Vec3)=transform(v,translation(t))
translate{T<:Vec4}(vs::AbstractVector{T},t::Vec3)=transform(vs,translation(t))
rotatexyz(v::Vec4,r::Vec3)=transform(v,rotationxyz(r))
rotatexyz{T<:Vec4}(vs::AbstractVector{T},r::Vec3)=transform(vs,rotationxyz(r))
function transrotatez(v::Vec4,t::Vec3,theta::Real)
  tm = translation(t)
  rm = rotationxyz(Vec3(0.0,0.0,theta))
  convert(Vec4,rm*tm*convert(Vector,v))
end
function transrotatez{T<:Vec4}(vs::AbstractVector{T},t::Vec3,theta::Real)
  tm = translation(t)
  rm = rotationxyz(Vec3(0.0,0.0,theta))
  convert(Vector{T},rm*tm*convert(Matrix,vs))
end


text:
using Combinatorics: length
using DataFrames: Matrix
using MAT, DelimitedFiles, JLD, JSON, CSV, LinearAlgebra, Combinatorics, DataFrames, Statistics, GlobalEnergyGIS

include("data_generation.jl")

# name the instance that you're creating
sspscenario_input = "ssp2-26"
sspyear_input = 2050
era_year_input = 2018

# select the folder in which you want to write the files, "folder" and "subfolder" should not change
folder = "examples"
subfolder = "structures"
instance = "Europe"

# Path to the output folder of GlobalEnergyGIS
inputdata =  "/scratch/work/beranej1/Copernicus/output"       #"D:\\Eigene Dateien\\Studium\\Master\\RA\\Copernicus\\output"

# Write the countries/regions for which you want to create the dataset. Use either country Names or predifined regions: Nordics, Eastern, Western, Central, Mediterranian
Regions = ["Finland", "Norway", "Sweden", "Denmark", "France", "Germany", "United Kingdom", "Spain", "Poland", "Italy"]

# Desired time period in hours
T = 8760
# number of technologies
t = 8       # 8 since hydro is not included here
## Value of Fmin
Fmin = 0.05

Dataset = get_countries()

create_data_sets(inputdata, sspscenario_input, sspyear_input, era_year_input, Dataset, folder, subfolder, instance, T, t, Fmin)


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using CSV
using DataFrames
using Dates
using Plots
using BSON
using StatsBase
using Underscores
import Downloads
using ColorSchemes
using PrettyTables
colors = ColorSchemes.tableau_10.colors
# using Makie
# using GLMakie
# GLMakie.activate!()
plotlyjs()
##

ps = let
    regions2key =  Dict(
        "Hungary" => "HU",
        "UK" => "GB",
        "France" => "FR",
        "Sweden" => "SE",
        "Denmark" => "DK",
        "Germany" => "DE",
        "Italy" => "IT",
        "USA" => "US")
    key2regions = Dict( c=>k for (k,c) in pairs(regions2key) )
    ks = collect(keys(key2regions))
    regions = collect(values(key2regions))
    (; regions2key, key2regions, keys=ks, regions)
end

function epidata(ps)
    @_ "https://storage.googleapis.com/covid19-open-data/v3/epidemiology.csv" |>
        Downloads.download |>
        CSV.read(__, DataFrame) |>
        filter(:location_key => ∈(ps.keys), __) |>
        select(__, :date, :new_deceased => :deaths, :new_confirmed => :cases,
            :location_key =>ByRow( x->ps.key2regions[x] )=> :region)
end

function popdata(ps)
    @_ "https://storage.googleapis.com/covid19-open-data/v3/demographics.csv" |>
        Downloads.download |>
        CSV.read(__, DataFrame) |>
        filter(:location_key => ∈(ps.keys), __) |>
        select(__, :population,
            :location_key =>ByRow( x->ps.key2regions[x] )=> :region)
end

function vacdata(ps)
    @_ "https://storage.googleapis.com/covid19-open-data/v3/vaccinations.csv" |>
        Downloads.download |>
        CSV.read(__, DataFrame) |>
        filter(:location_key => ∈(ps.keys), __) |>
        select(__, :date, :cumulative_persons_vaccinated => :vaccinated,
            :location_key =>ByRow( x->ps.key2regions[x] )=> :region)
end

function disallownegative!(df)
    for c in [:deaths, :cases]
        idx = findall(<(0), df[!,c])
        if !isempty(idx)
            @info "replace negative counts with zero"
            pretty_table(df[idx,:])
            df[idx,c] .= 0
        end
    end
    return df
end

function rollingmean(xs, window)
    n  = length(xs)-window+1
    ys = Vector{Float64}(undef, n)
    for i in 1:n
        ys[i] = mean(xs[i:i+window-1])
    end
    return ys
end

function smoothing(df, window)
    @assert isodd(window)
    Δ = window ÷ 2
    smoothed = DataFrame(
        :date => df.date[1+Δ:end-Δ],
        :cases => df.cases[1+Δ:end-Δ],
        :deaths => df.deaths[1+Δ:end-Δ], )
    for col in ["cases", "deaths"]
        smoothed[!,"$(col)_smoothed"] = rollingmean(df[:,col], window)
    end
    smoothed
end

function replacemissing!(df)
    cols = filter(x->nonmissingtype( eltype( df[!,x] ) ) <: Number, names(df))
    for col in cols
        for i in findall(ismissing, df[!,col])
            df[i, col] = i == 1 ? 0 : df[i-1, col]
        end
    end
    df
end

function addpadding(df, dates)
    cols = filter(x->nonmissingtype( eltype( df[!,x] ) ) <: Number, names(df))
    padded = DataFrame( :date => dates )
    for col in cols
        df[2:end,col] = diff(df[:,col])
    end
    padded =@_ leftjoin(padded, df; on=:date) |>
        sort(__, :date) |>
        transform(__, cols .=> x->cumsum(coalesce.(x,0)), renamecols=false)
    padded
end

function main(epi, vac, pop)
    epidata = Dict{String,DataFrame}()
    for region in ps.regions
        @info region
        regional = @_ epi |>
            filter(:region => ==(region), __) |>
            sort(__, :date) |>
            filter(:date => <=(Date("2021-09-01")), __) |>
            disallowmissing |>
            disallownegative! |>
            smoothing(__, 7)

        vaccination = @_ vac |>
            filter(:region => ==(region), __) |>
            select(__, Not(:region)) |>
            sort(__, :date) |>
            filter(:date => <(Date("2021-09-01")), __) |>
            replacemissing!(__) |>
            disallowmissing |>
            addpadding(__,regional.date)

        population = @_ pop |>
            filter(:region => ==(region), __) |>
            only(__.population)

        regional[!,:population] .= population
        regional[!,:cases_per_100000] = regional[!,"cases_smoothed"] / population * 100_000
        regional[!,:deaths_per_100000] = regional[!,"deaths_smoothed"] / population * 100_000
        vaccination[!,:fraction_vaccinated] = vaccination[!,"vaccinated"] / population
        regional = leftjoin(regional, vaccination; on=:date)
        disallowmissing!(regional)

        @assert all(regional.date .== regional.date[1]:Day(1):regional.date[end])
        CSV.write(projectdir("data/epidata/", "epidata_region=$region.csv"), regional)
        epidata[region] = regional
    end
    bson(projectdir("data/epidata/", "epidata.bson"), epidata)
    epidata
end

## ========================================================================
data = main(
    epidata(ps),
    vacdata(ps),
    popdata(ps)
    )

## =========================================================================
#          merge data for epidemia
function dict2df(data)
    regions = collect(keys(data))
    r = regions[1]
    epi = @_ data[r] |>
        select(__, All(),
            :date => ( x -> div.(0:length(x)-1, 7) .+ 1 ) => :week)
    epi[!, :region] .= r
    for r in regions[2:end]
        df = @_ data[r] |>
            select(__, All(),
                :date => ( x -> div.(0:length(x)-1, 7) .+ 1 ) => :week)
        df[!, :region] .= r
        epi = vcat(epi, df)
    end
    epi
end
epi = dict2df(data)
CSV.write(projectdir("data/epidata/", "epidata_longformat.csv"), epi)

## ==============================================================
#     DEPRICATED WHO DATA

who_params = (
    regions = ["Hungary", "The United Kingdom", "France", "Sweden", "Denmark", "Germany", "Italy", "United States of America"]
    , old2new = Dict(
        "Hungary" => "Hungary",
        "The United Kingdom" => "UK",
        "France" => "France",
        "Sweden" => "Sweden",
        "Denmark" => "Denmark",
        "Germany" => "Germany",
        "Italy" => "Italy",
        "United States of America" => "USA"
    )
    , fname = projectdir("data/epinow/WHO-COVID-19-global-data.csv")
)

let
    p1 = bar(rt.date, rt.new_cases, fc=colors[1], lc=colors[1])
    p2 = plot(rt.date, rt.mean, ribbon=(rt.lower_80, rt.upper_80))
    p = plot(p1, p2, layout=(2,1), size=(1000, 2*250), sharex=true, link=:x)
end
who = CSV.read(who_params.fname, DataFrame)

function load_data(regions, params)
    who = CSV.read(params.fname, DataFrame)
    observations = Dict{String, DataFrame}()
    for region in params.regions
        regional = @_ who |>
                filter(:Country => ==(region), __) |>
                select(__, :Date_reported => :date, :New_cases => :cases, :New_deaths => :deaths) #|>
        #__[Date("2020-08-01") <. __.date <. Date("2021-02-01") ]
        observations[params.old2new[region]] = regional
    end
    observations
end

obs = load_data(regions, who_params)

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
__precompile__(false)

"""

`RAFF.jl` is a Jula package.

"""
module RAFF

# Dependencies
using Distributed
using ForwardDiff
using LinearAlgebra
using Statistics
using Printf
using Random
using SharedArrays
using Logging

export lmlovo, raff, praff

# Set RAFF logger
raff_logger = ConsoleLogger(stdout, Logging.Error)

lm_logger = ConsoleLogger(stdout, Logging.Error)

# Load code

include("raffoutput.jl")

include("utils.jl")

include("dutils.jl")

include("generator.jl")

"""
    lmlovo(model::Function [, θ::Vector{Float64} = zeros(n)], data::Array{Float64, 2},
           n::Int, p::Int [; kwargs...])

    lmlovo(model::Function, gmodel!::Function [, θ::Vector{Float64} = zeros(n)],
           data::Array{Float64,2}, n::Int, p::Int [; MAXITER::Int=200,
           ε::Float64=10.0^-4])

Fit the `n`-parameter model `model` to the data given by matrix
`data`. The strategy is based on the LOVO function, which means that
only `p` (0 < `p` <= rows of `data`) points are trusted. The
Levenberg-Marquardt algorithm is implemented in this version.

Matriz `data` is the data to be fit. This matrix should be in the form

    x11 x12 ... x1N y1
    x21 x22 ... x2N y2
    :

where `N` is the dimension of the argument of the model
(i.e. dimension of `x`).

If `θ` is provided, then it is used as the starting point.

The signature of function `model` should be given by

    model(x::Union{Vector{Float64}, SubArray}, θ::Vector{Float64})

where `x` are the variables and `θ` is a `n`-dimensional vector of
parameters. If the gradient of the model `gmodel!`

    gmodel! = (g::SubArray, x::Union{Vector{Float64}, SubArray},
               θ::Vector{Float64})

is not provided, then the function ForwardDiff.gradient! is called to
compute it.  **Note** that this choice has an impact in the
computational performance of the algorithm. In addition, if
`ForwardDiff.jl` is being used, then one **MUST** remove the signature
of vector `θ` from function `model`.

The optional arguments are

  - `MAXITER`: maximum number of iterations
  - `ε`: tolerance for the gradient of the function

Returns a [`RAFFOutput`](@ref) object.

"""
function lmlovo(model::Function, gmodel!::Function, θ::Vector{Float64},
                data::Array{Float64,2}, n::Int, p::Int;
                MAXITER::Int=200, ε::Float64=10.0^-4)

    @assert(n > 0, "Dimension should be positive.")
    @assert(p >= 0, "Trusted points should be nonnegative.")

    npun, = size(data)

    with_logger(lm_logger) do

        @debug("Size of data matrix ", size(data))

    end

    # Counters for calls to F and its Jacobian
    nf = 0

    nj = 0

    (p == 0) && return RAFFOutput(1, θ, 0, p, 0.0, nf, nj, [1:npun;])

    # Main function - the LOVO function
    LovoFun = let

        npun_::Int = npun

        ind::Vector{Int} = Vector{Int}(undef, npun_)

        F::Vector{Float64} = Vector{Float64}(undef, npun_)

        p_::Int = p

        # Return a ordered set index and lovo value
        (θ) -> begin

            nf += 1

            @views for i = 1:npun_
                F[i] = (model(data[i,1:(end - 1)], θ) - data[i, end])^2
            end

            indF, orderedF = sort_fun!(F, ind, p_)

            return indF, sum(orderedF)
        end

    end

    # Residue and Jacobian of residue
    val_res::Vector{Float64} = Vector{Float64}(undef, p)

    jac_res::Array{Float64, 2} = Array{Float64}(undef, p, n)

    # This function returns the residue and Jacobian of residue
    ResFun!(θ::Vector{Float64}, ind, r::Vector{Float64},
            rJ::Array{Float64, 2}) = begin

       nj += 1

       for (k, i) in enumerate(ind)

            x = data[i, 1:(end - 1)]

            r[k] = model(x, θ) - data[i, end]

            v = @view(rJ[k, :])

            gmodel!(v, x, θ)

        end

    end

    # Levenberg-Marquardt algorithm
    # -----------------------------

    Id = Matrix(1.0I, n, n)

    # Status = 1 means success
    status = 1

    # Parameters

    λ_up      = 2.0
    λ_down    = 2.0
    λ         = 1.0
    maxoutind = min(p, 5)

    # Allocation
    θnew = Vector{Float64}(undef, n)
    d = Vector{Float64}(undef, n)
    y = Vector{Float64}(undef, n)
    G = Array{Float64, 2}(undef, n, n)
    grad_lovo = Vector{Float64}(undef, n)

    # Initial steps

    ind_lovo, best_val_lovo = LovoFun(θ)

    ResFun!(θ, ind_lovo, val_res, jac_res)

    BLAS.gemv!('T', 1.0, jac_res, val_res, 0.0, grad_lovo)

    ngrad_lovo = norm(grad_lovo, 2)

    safecount = 1

    # Main loop

    while (ngrad_lovo >= ε) && (safecount < MAXITER)

        with_logger(lm_logger) do

            @info("Iteration $(safecount)")
            @info("  Current value:   $(best_val_lovo)")
            @info("  ||grad_lovo||_2: $(ngrad_lovo)")
            @info("  Current iterate: $(θ)")
            @info("  Best indices (first $(maxoutind)): $(ind_lovo[1:maxoutind])")
            @info("  lambda: $(λ)")

        end

        G .= Id

        BLAS.gemm!('T', 'N', 1.0, jac_res, jac_res, λ, G)

        F = qr!(G)
        #F = cholesky!(G, Val(true))

        ad = try

            ldiv!(d, F, grad_lovo)

            d .*= -1.0

        catch
            "error"
        end
        if ad == "error" #restarting if lapack fails
            with_logger(lm_logger) do

                @warn "Failed to solve the linear system. Will try new point."
                d .= - 1.0 .* grad_lovo
                θ .= rand(n)

            end
        else
            d .= ad
        end

        θnew .= θ .+ d

        ind_lovo, val_lovo = LovoFun(θnew)

        if  val_lovo <= best_val_lovo

            θ .= θnew

            best_val_lovo = val_lovo

            λ = λ / λ_down

            ResFun!(θ, ind_lovo, val_res, jac_res)

            BLAS.gemv!('T', 1.0, jac_res, val_res, 0.0, grad_lovo)

            ngrad_lovo = norm(grad_lovo, 2)

            with_logger(lm_logger) do

                @info("  Better function value found, lambda changed to $(λ).")

            end

        else

            λ = λ * λ_up

            with_logger(lm_logger) do

                @info("  No improvement, lambda changed to $(λ).")

            end

        end

        safecount += 1

    end

    if safecount == MAXITER
        with_logger(lm_logger) do

            @info("No solution was found in $(safecount) iterations.")

        end
        status = 0
    end

    # TODO: Create a test for this case
    if isnan(ngrad_lovo)
        with_logger(lm_logger) do

            @info("Incorrect value for gradient norm $(ngrad_lovo).")

        end
        status = 0
    end

    outliers = [1:npun;]
    setdiff!(outliers, ind_lovo)

    with_logger(lm_logger) do

        @info("""

        Final iteration (STATUS=$(status))
          Solution found:       $(θ)
          ||grad_lovo||_2:      $(ngrad_lovo)
          Function value:       $(best_val_lovo)
          Number of iterations: $(safecount)
          Outliers:             $(outliers)

        """)

    end

    return RAFFOutput(status, θ, safecount, p, best_val_lovo, nf, nj, outliers)

end

function lmlovo(model::Function, θ::Vector{Float64}, data::Array{Float64,2},
                n::Int, p::Int; kwargs...)

    # Define closures for derivative and initializations

    # 'x' is considered as global parameter for this function
    model_cl(θ) = model(x, θ)

    grad_model!(g, x_, θ) = begin

        global x = x_

        ForwardDiff.gradient!(g, model_cl, θ)

    end

    return lmlovo(model, grad_model!, θ, data, n, p; kwargs...)

end

lmlovo(model::Function, gmodel!::Function, data::Array{Float64,2}, n::Int, p::Int; kwargs...) =
    lmlovo(model, gmodel!, zeros(Float64, n), data, n, p; kwargs...)

lmlovo(model::Function, data::Array{Float64,2}, n::Int, p::Int; kwargs...) =
    lmlovo(model, zeros(Float64, n), data, n, p; kwargs...)

"""

    raff(model::Function, data::Array{Float64, 2}, n::Int; kwargs...)

    raff(model::Function, gmodel!::Function, data::Array{Float64, 2},
        n::Int; MAXMS::Int=1, SEEDMS::Int=123456789,
        initguess::Vector{Float64}=zeros(Float64, n),
        ε::Float64=1.0e-4, noutliers::Int=-1, ftrusted::Union{Float64,
        Tuple{Float64, Float64}}=0.5)

Robust Algebric Fitting Function (RAFF) algorithm. This function uses
a voting system to automatically find the number of trusted data
points to fit the `model`.

  - `model`: function to fit data. Its signature should be given by

        model(x, θ)

    where `x` is the multidimensional argument and `θ` is the
    `n`-dimensional vector of parameters

  - `gmodel!`: gradient of the model function. Its signature should be
    given by

        gmodel!(g, x, θ)

    where `x` is the multidimensional argument, `θ` is the
    `n`-dimensional vector of parameters and the gradient is written
    in `g`.

  - `data`: data to be fit. This matrix should be in the form

        x11 x12 ... x1N y1
        x21 x22 ... x2N y2
        :

    where `N` is the dimension of the argument of the model
    (i.e. dimension of `x`).

  - `n`: dimension of the parameter vector in the model function

The optional arguments are

  - `MAXMS`: number of multistart points to be used
  - `SEEDMS`: integer seed for random multistart points
  - `initialguess`: a good guess for the starting point and for
    generating random points in the multistart strategy
  - `ε`: gradient stopping criteria to `lmlovo`
  - `noutliers`: integer describing the maximum expected number of
    outliers. The default is *half*. *Deprecated*.
  - `ftrusted`: float describing the minimum expected percentage of
    trusted points. The default is *half* (0.5). Can also be a
    Tuple of the form `(fmin, fmax)` percentages of trusted points.

Returns a [`RAFFOutput`](@ref) object with the best parameter found.

"""
function raff(model::Function, gmodel!::Function,
    data::Array{Float64, 2}, n::Int; MAXMS::Int=1,
    SEEDMS::Int=123456789, initguess::Vector{Float64}=zeros(Float64,
    n), ε::Float64=1.0e-4, noutliers::Int=-1,
    ftrusted::Union{Float64, Tuple{Float64, Float64}}=0.5)

    np, = size(data)

    if noutliers != -1

        Base.depwarn("Optional argument `noutliers::Int` is deprecated and will be removed from future releases. Use `ftrusted::Float` or `itrusted::Tuple{Float, Float}` instead.", :raff)

        ftrusted = (noutliers >= 0) ? max(0, np - noutliers) / np : 0.5

    end

    # Initialize random generator
    seedMS = MersenneTwister(SEEDMS)

    # Define interval for trusted points

    pliminf, plimsup = try

        check_ftrusted(ftrusted, np)

    catch e

        with_logger(raff_logger) do

            @error("Error in optional parameter `ftrusted`.", e)

        end

        return RAFFOutput()

    end

    lv = plimsup - pliminf + 1

    sols = Vector{RAFFOutput}(undef, lv)

    for i = pliminf:plimsup

        vbest = RAFFOutput(initguess, i)

        ind = i - pliminf + 1

        for j = 1:MAXMS

            with_logger(raff_logger) do

                @debug("Running lmlovo for p = $(i). Repetition $(j).")

            end

            # Starting point
            θ = randn(seedMS, Float64, n)
            θ .= θ .+ initguess

            # Call function and store results
            sols[ind] = lmlovo(model, gmodel!, θ, data, n, i; ε=ε, MAXITER=400)

            # Update the best point and functional value
            (sols[ind].status == 1) && (sols[ind].f < vbest.f) && (vbest = sols[ind])

        end

        sols[ind] = vbest

        with_logger(raff_logger) do

            @debug("Best solution for p = $(i).", vbest.solution)

        end

    end

    # Count the total number of iterations, and function and Jacobian
    # evaluations.

    nf = 0
    nj = 0
    ni = 0

    for s in sols

        nf += s.nf
        nj += s.nj
        ni += s.iter

    end

    # Apply the filter and the voting strategy to all the solutions
    # found

    votsis = with_logger(raff_logger) do

        voting_strategy(model, data, sols, pliminf, plimsup)

    end

    mainind = findlast(x->x == maximum(votsis), votsis)

    s = sols[mainind]

    return RAFFOutput(s.status, s.solution, ni, s.p, s.f, nf, nj, s.outliers)

end

function raff(model::Function, data::Array{Float64, 2}, n::Int; kwargs...)

    # Define closures for derivative and initializations

    # 'x' is considered as global for this function
    model_cl(θ) = model(x, θ)

    grad_model!(g, x_, θ) = begin

        global x = x_

        return ForwardDiff.gradient!(g, model_cl, θ)

    end

    return raff(model, grad_model!, data, n; kwargs...)

end

"""

    praff(model::Function, data::Array{Float64, 2}, n::Int; kwargs...)

    praff(model::Function, gmodel!::Function, data::Array{Float64, 2},
        n::Int; MAXMS::Int=1, SEEDMS::Int=123456789, batches::Int=1,
        initguess::Vector{Float64}=zeros(Float64, n),
        ε::Float64=1.0e-4, noutliers::Int=-1, ftrusted::Union{Float64,
        Tuple{Float64, Float64}}=0.5)

Multicore distributed version of RAFF. See the description of the
[`raff`](@ref) function for the main (non-optional) arguments. All the
communication is performed by channels.

This function uses all available **local** workers to run RAFF
algorithm. Note that this function does not use *Tasks*, so all the
parallelism is based on the [Distributed](https://docs.julialang.org/en/latest/manual/parallel-computing/#Multi-Core-or-Distributed-Processing-1) package.

The optional arguments are

  - `MAXMS`: number of multistart points to be used
  - `SEEDMS`: integer seed for random multistart points
  - `batches`: size of batches to be send to each worker
  - `initguess`: starting point to be used in the multistart procedure
  - `ε`: stopping tolerance
  - `noutliers`: integer describing the maximum expected number of
    outliers. The default is *half*. *Deprecated*.
  - `ftrusted`: float describing the minimum expected percentage of
    trusted points. The default is *half* (0.5). Can also be a
    Tuple of the form `(fmin, fmax)` percentages of trusted points.

Returns a [`RAFFOutput`](@ref) object containing the solution.

"""
function praff(model::Function, gmodel!::Function,
               data::Array{Float64, 2}, n::Int; MAXMS::Int=1,
               SEEDMS::Int=123456789, batches::Int=1,
               initguess::Vector{Float64}=zeros(Float64, n),
               ε::Float64=1.0e-4, noutliers::Int=-1,
               ftrusted::Union{Float64, Tuple{Float64, Float64}}=0.5)

    np, = size(data)

    if noutliers != -1

        Base.depwarn("Optional argument `noutliers::Int` is deprecated and will be removed from future releases. Use `ftrusted::Float` or `itrusted::Tuple{Float, Float}` instead.", :raff)

        ftrusted = (noutliers >= 0) ? max(0, np - noutliers) / np : 0.5

    end

    # Initialize random generator
    seedMS = MersenneTwister(SEEDMS)

    # Define interval for trusted points

    pliminf, plimsup = try

        check_ftrusted(ftrusted, np)

    catch e

        with_logger(raff_logger) do

            @error("Error in optional parameter `ftrusted`.", e)

        end

        return RAFFOutput()

    end

    lv = plimsup - pliminf + 1

    # Create a RemoteChannel to receive solutions
    bqueue = RemoteChannel(() -> Channel{Vector{Float64}}(div(lv, 2)))
    # TODO: Check a smart way for not creating a large channel
    squeue = RemoteChannel(() -> Channel{RAFFOutput}(lv))
    # Create another channel to assign tasks
    tqueue = RemoteChannel(() -> Channel{UnitRange{Int}}(0))

    # This command selects only nodes which are local to myid()
    curr_workers = workers()

    futures = Vector{Future}(undef, length(curr_workers))

    # Start updater Task
    # This task is not needed up to now.
    # @async with_logger(()-> update_best(bqueue, bestx), raff_logger)

    # Start workers Tasks (CPU intensive)
    with_logger(raff_logger) do

        @debug("Workers", curr_workers)

    end

    for (i, t) in enumerate(curr_workers)

        futures[i] = @spawnat(t, with_logger( ()-> try

            @debug("Creating worker $(t).")

            consume_tqueue(bqueue, tqueue, squeue,
                           model, gmodel!, data, n, pliminf,
                           plimsup, MAXMS, seedMS, initguess)
            catch e

               @error("Unable to start worker $(t).", e)

            end, raff_logger
        ))

        with_logger(raff_logger) do

            @debug("Created worker $(t).")

        end
    end

    # Check asynchronously if there is at least one live worker
    @async with_logger(
        () -> check_and_close(bqueue, tqueue, squeue, futures),
        raff_logger)

    # Populate the task queue with jobs

    for p = pliminf:batches:plimsup

        try

            put!(tqueue, p:min(plimsup, p + batches - 1))

        catch e

            with_logger(raff_logger) do

                @warn("Tasks queue prematurely closed while inserting tasks. Will exit.")

            end

            break

        end

        with_logger(raff_logger) do

            @debug("Added problem $(p) to tasks queue.")

        end

    end

    # The task queue can be closed, since all the problems have been
    # read, due to the size 0 of this channel
    close(tqueue)

    with_logger(raff_logger) do

        @debug("Waiting for workers to finish.")

    end

    # Create a vector of solutions to store the results from workers
    sols = Vector{RAFFOutput}(undef, lv)

    for i in 1:lv

        try

            rout = take!(squeue)

            sols[rout.p - pliminf + 1] = rout

            with_logger(raff_logger) do

                @debug("Stored solution for p=$(rout.p).")

            end

        catch e

            with_logger(raff_logger) do

                @error("Error when retrieving solutions.", e)

            end

        end

    end

    close(bqueue)

    close(squeue)

    # Count the total number of iterations, and function and Jacobian
    # evaluations.

    nf = 0
    nj = 0
    ni = 0

    for s in sols

        nf += s.nf
        nj += s.nj
        ni += s.iter

    end

    # Apply the filter and the voting strategy to all the solutions
    # found

    votsis = with_logger(raff_logger) do

        voting_strategy(model, data, sols, pliminf, plimsup)

    end

    mainind = findlast(x->x == maximum(votsis), votsis)

    s = sols[mainind]

    return RAFFOutput(s.status, s.solution, ni, s.p, s.f, nf, nj, s.outliers)

end

function praff(model::Function, data::Array{Float64, 2}, n::Int; kwargs...)

    # Define closures for derivative and initializations

    # 'x' is considered as global parameter for this function
    model_cl(θ) = model(x, θ)

    grad_model!(g, x_, θ) = begin

        global x = x_

        return ForwardDiff.gradient!(g, model_cl, θ)

    end

    return praff(model, grad_model!, data, n; kwargs...)

end

end


text:
function _update_parameter_values!(
    ::AbstractArray{T},
    ::NoAttributes,
    args...,
) where {T <: Union{Float64, PJ.ParameterRef}} end

######################## Methods to update Parameters from Time Series #####################
function _set_param_value!(param::JuMPParamArray, value::Float64, name::String, t::Int)
    JuMP.set_value(param[name, t], value)
    return
end

function _set_param_value!(
    param::DenseAxisArray{Vector{NTuple{2, Float64}}},
    value::Vector{NTuple{2, Float64}},
    name::String,
    t::Int,
)
    param[name, t] = value
    return
end

function _set_param_value!(param::JuMPFloatArray, value::Float64, name::String, t::Int)
    param[name, t] = value
    return
end

function _set_param_value!(
    param::SparseAxisArray,
    value::Float64,
    name::String,
    subcomp::String,
    t::Int,
)
    _set_parameter_value_sparse_array!(param[name, subcomp, t], value)
    return
end

function _set_parameter_value_sparse_array!(parameter::Float64, value::Float64)
    parameter = value
    return
end

function _set_parameter_value_sparse_array!(parameter::PJ.ParameterRef, value::Float64)
    JuMP.set_value(parameter, value)
    return
end

function _update_parameter_values!(
    param_array::AbstractArray{T},
    attributes::TimeSeriesAttributes{U},
    ::Type{V},
    model::DecisionModel,
    ::DatasetContainer{DataFrameDataset},
) where {
    T <: Union{PJ.ParameterRef, Float64},
    U <: PSY.AbstractDeterministic,
    V <: PSY.Component,
}
    initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels
    horizon = get_time_steps(get_optimization_container(model))[end]
    components = get_available_components(V, get_system(model))
    for component in components
        name = PSY.get_name(component)
        ts_vector = get_time_series_values!(
            U,
            model,
            component,
            get_time_series_name(attributes),
            get_time_series_multiplier_id(attributes),
            initial_forecast_time,
            horizon,
        )
        for (t, value) in enumerate(ts_vector)
            _set_param_value!(param_array, value, name, t)
        end
    end
end

function _update_parameter_values!(
    param_array::SparseAxisArray,
    attributes::TimeSeriesAttributes{U},
    ::Type{V},
    model::DecisionModel,
    ::DatasetContainer{DataFrameDataset},
) where {U <: PSY.AbstractDeterministic, V <: PSY.HybridSystem}
    initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels
    horizon = get_time_steps(get_optimization_container(model))[end]
    components = get_available_components(V, get_system(model))
    for component in components, subcomp_type in [PSY.RenewableGen, PSY.ElectricLoad]
        name = PSY.get_name(component)
        !does_subcomponent_exist(component, subcomp_type) && continue
        subcomponent = get_subcomponent(component, subcomp_type)
        ts_vector = get_time_series_values!(
            U,
            model,
            component,
            make_subsystem_time_series_name(subcomponent, get_time_series_name(attributes)),
            get_time_series_multiplier_id(attributes),
            initial_forecast_time,
            horizon,
        )
        for (t, value) in enumerate(ts_vector)
            _set_param_value!(param_array, value, name, string(subcomp_type), t)
        end
    end
end

function _update_parameter_values!(
    param_array::AbstractArray{T},
    attributes::TimeSeriesAttributes{U},
    service::V,
    model::DecisionModel,
    ::DatasetContainer{DataFrameDataset},
) where {
    T <: Union{PJ.ParameterRef, Float64},
    U <: PSY.AbstractDeterministic,
    V <: PSY.Service,
}
    initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels
    horizon = get_time_steps(get_optimization_container(model))[end]
    ts_vector = get_time_series_values!(
        U,
        model,
        service,
        get_time_series_name(attributes),
        get_time_series_multiplier_id(attributes),
        initial_forecast_time,
        horizon,
    )
    service_name = PSY.get_name(service)
    for (t, value) in enumerate(ts_vector)
        _set_param_value!(param_array, value, service_name, t)
    end
end

function _update_parameter_values!(
    param_array::AbstractArray{T},
    attributes::TimeSeriesAttributes{U},
    ::Type{V},
    model::EmulationModel,
    ::DatasetContainer{DataFrameDataset},
) where {T <: Union{PJ.ParameterRef, Float64}, U <: PSY.SingleTimeSeries, V <: PSY.Device}
    initial_forecast_time = get_current_time(model)
    components = get_available_components(V, get_system(model))
    for component in components
        # Note: This interface reads one single value per component at a time.
        ts_vector = get_time_series_values!(
            U,
            model,
            component,
            get_time_series_name(attributes),
            get_time_series_multiplier_id(attributes),
            initial_forecast_time,
        )
        _set_param_value!(param_array, ts_vector[1], PSY.get_name(component), 1)
    end
    return
end

function _update_parameter_values!(
    param_array::AbstractArray{T},
    attributes::VariableValueAttributes,
    ::Type{<:PSY.Component},
    model::DecisionModel,
    state::DatasetContainer{DataFrameDataset},
) where {T <: Union{PJ.ParameterRef, Float64}}
    current_time = get_current_time(model)
    state_values = get_dataset_values(state, get_attribute_key(attributes))
    component_names, time = axes(param_array)
    resolution = get_resolution(model)

    state_data = get_dataset(state, get_attribute_key(attributes))
    state_timestamps = state_data.timestamps
    max_state_index = length(state_data)

    state_data_index = find_timestamp_index(state_timestamps, current_time)

    sim_timestamps = range(current_time, step=resolution, length=time[end])
    for t in time
        timestamp_ix = min(max_state_index, state_data_index + 1)
        @debug "parameter horizon is over the step" max_state_index > state_data_index + 1
        if state_timestamps[timestamp_ix] <= sim_timestamps[t]
            state_data_index = timestamp_ix
        end
        for name in component_names
            # Pass indices in this way since JuMP DenseAxisArray don't support view()
            _set_param_value!(param_array, state_values[state_data_index, name], name, t)
        end
    end
    return
end

function _update_parameter_values!(
    param_array::AbstractArray{T},
    attributes::VariableValueAttributes{VariableKey{OnVariable, U}},
    ::Type{U},
    model::DecisionModel,
    state::DatasetContainer{DataFrameDataset},
) where {T <: Union{PJ.ParameterRef, Float64}, U <: PSY.Component}
    current_time = get_current_time(model)
    state_values = get_dataset_values(state, get_attribute_key(attributes))
    component_names, time = axes(param_array)
    resolution = get_resolution(model)

    state_data = get_dataset(state, get_attribute_key(attributes))
    state_timestamps = state_data.timestamps
    max_state_index = length(state_data)

    state_data_index = find_timestamp_index(state_timestamps, current_time)

    sim_timestamps = range(current_time, step=resolution, length=time[end])
    for t in time
        timestamp_ix = min(max_state_index, state_data_index + 1)
        @debug "parameter horizon is over the step" max_state_index > state_data_index + 1
        if state_timestamps[timestamp_ix] <= sim_timestamps[t]
            state_data_index = timestamp_ix
        end
        for name in component_names
            # Pass indices in this way since JuMP DenseAxisArray don't support view()
            val = round(state_values[state_data_index, name])
            @assert 0.0 <= val <= 1.0
            _set_param_value!(param_array, val, name, t)
        end
    end
    return
end

function _update_parameter_values!(
    param_array::AbstractArray{T},
    attributes::VariableValueAttributes,
    ::Type{<:PSY.Component},
    model::EmulationModel,
    state::DatasetContainer{DataFrameDataset},
) where {T <: Union{PJ.ParameterRef, Float64}}
    current_time = get_current_time(model)
    state_values = get_dataset_values(state, get_attribute_key(attributes))
    component_names, _ = axes(param_array)
    state_data = get_dataset(state, get_attribute_key(attributes))
    state_timestamps = state_data.timestamps
    state_data_index = find_timestamp_index(state_timestamps, current_time)
    for name in component_names
        # Pass indices in this way since JuMP DenseAxisArray don't support view()
        _set_param_value!(param_array, state_values[state_data_index, name], name, 1)
    end
    return
end

function _update_parameter_values!(
    param_array::AbstractArray{T},
    attributes::VariableValueAttributes{VariableKey{OnVariable, U}},
    ::Type{<:PSY.Component},
    model::EmulationModel,
    state::DatasetContainer{DataFrameDataset},
) where {T <: Union{PJ.ParameterRef, Float64}, U <: PSY.Component}
    current_time = get_current_time(model)
    state_values = get_dataset_values(state, get_attribute_key(attributes))
    component_names, _ = axes(param_array)
    state_data = get_dataset(state, get_attribute_key(attributes))
    state_timestamps = state_data.timestamps
    state_data_index = find_timestamp_index(state_timestamps, current_time)
    for name in component_names
        # Pass indices in this way since JuMP DenseAxisArray don't support view()
        val = round(state_values[state_data_index, name])
        @assert 0.0 <= val <= 1.0
        _set_param_value!(param_array, val, name, 1)
    end
    return
end

function _update_parameter_values!(
    ::AbstractArray{T},
    ::VariableValueAttributes,
    ::Type{<:PSY.Component},
    ::EmulationModel,
    ::EmulationModelStore,
) where {T <: Union{PJ.ParameterRef, Float64}}
    error("The emulation model has parameters that can't be updated from its results")
    return
end

"""
Update parameter function an OperationModel
"""
function update_parameter_values!(
    model::OperationModel,
    key::ParameterKey{T, U},
    input::DatasetContainer{DataFrameDataset},
) where {T <: ParameterType, U <: PSY.Component}
    # Enable again for detailed debugging
    # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin
    optimization_container = get_optimization_container(model)
    # Note: Do not instantite a new key here because it might not match the param keys in the container
    # if the keys have strings in the meta fields
    parameter_array = get_parameter_array(optimization_container, key)
    parameter_attributes = get_parameter_attributes(optimization_container, key)
    _update_parameter_values!(parameter_array, parameter_attributes, U, model, input)
    IS.@record :execution ParameterUpdateEvent(
        T,
        U,
        parameter_attributes,
        get_current_timestamp(model),
        get_name(model),
    )
    # end
    return
end

function update_parameter_values!(
    model::OperationModel,
    key::ParameterKey{T, U},
    input::DatasetContainer{DataFrameDataset},
) where {T <: EnergyLimitParameter, U <: PSY.Generator}
    # Enable again for detailed debugging
    # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin
    optimization_container = get_optimization_container(model)
    # Note: Do not instantite a new key here because it might not match the param keys in the container
    # if the keys have strings in the meta fields
    parameter_array = get_parameter_array(optimization_container, key)
    parameter_attributes = get_parameter_attributes(optimization_container, key)
    internal = get_internal(model)
    execution_count = internal.execution_count
    current_time = get_current_time(model)
    state_values = get_dataset_values(input, get_attribute_key(parameter_attributes))
    component_names, time = axes(parameter_array)
    resolution = get_resolution(model)
    interval_time_steps = Int(get_interval(model.internal.store_parameters) / resolution)
    state_data = get_dataset(input, get_attribute_key(parameter_attributes))
    state_timestamps = state_data.timestamps
    max_state_index = length(state_data)

    state_data_index = find_timestamp_index(state_timestamps, current_time)
    sim_timestamps = range(current_time, step=resolution, length=time[end])
    old_parameter_values = JuMP.value.(parameter_array)
    # The current method uses older parameter values because when passing the energy output from one stage
    # to the next, the aux variable values gets over-written by the lower level model after its solve.
    # This approach is a temporary hack and will be replaced in future versions.
    for t in time
        timestamp_ix = min(max_state_index, state_data_index + 1)
        @debug "parameter horizon is over the step" max_state_index > state_data_index + 1
        if state_timestamps[timestamp_ix] <= sim_timestamps[t]
            state_data_index = timestamp_ix
        end
        for name in component_names
            # the if statement checks if its the first solve of the model and uses the values stored in the state
            # and for subsequent solves uses the state data to update the parameter values for the last set of time periods
            # that are equal to the length of the interval i.e. the time periods that dont overlap between each solves.
            if execution_count == 0 || t > time[end] - interval_time_steps
                # Pass indices in this way since JuMP DenseAxisArray don't support view()
                _set_param_value!(
                    parameter_array,
                    state_values[state_data_index, name],
                    name,
                    t,
                )
            else
                # Currently the update method relies on using older parameter values of the EnergyLimitParameter
                # to update the parameter for overlapping periods between solves i.e. we ingoring the parameter values
                # in the model interval time periods.
                _set_param_value!(
                    parameter_array,
                    old_parameter_values[name, t + interval_time_steps],
                    name,
                    t,
                )
            end
        end
    end

    IS.@record :execution ParameterUpdateEvent(
        T,
        U,
        parameter_attributes,
        get_current_timestamp(model),
        get_name(model),
    )
    return
end

"""
Update parameter function an OperationModel
"""
function update_parameter_values!(
    model::OperationModel,
    key::ParameterKey{T, U},
    input::DatasetContainer{DataFrameDataset},
) where {T <: ParameterType, U <: PSY.Service}
    # Enable again for detailed debugging
    # TimerOutputs.@timeit RUN_SIMULATION_TIMER "$T $U Parameter Update" begin
    optimization_container = get_optimization_container(model)
    # Note: Do not instantite a new key here because it might not match the param keys in the container
    # if the keys have strings in the meta fields
    param_array = get_parameter_array(optimization_container, key)
    parameter_attributes = get_parameter_attributes(optimization_container, key)
    service = PSY.get_component(U, get_system(model), key.meta)
    _update_parameter_values!(param_array, parameter_attributes, service, model, input)
    IS.@record :execution ParameterUpdateEvent(
        T,
        U,
        parameter_attributes,
        get_current_timestamp(model),
        get_name(model),
    )
    #end
    return
end

function _update_parameter_values!(
    param_array,
    attributes::CostFunctionAttributes,
    ::Type{V},
    model::DecisionModel,
    ::DatasetContainer{DataFrameDataset},
) where {V <: PSY.Component}
    initial_forecast_time = get_current_time(model) # Function not well defined for DecisionModels
    time_steps = get_time_steps(get_optimization_container(model))
    horizon = time_steps[end]
    container = get_optimization_container(model)
    if is_synchronized(container)
        obj_func = get_objective_function(container)
        set_synchronized_status(obj_func, false)
        reset_variant_terms(obj_func)
    end
    components = get_available_components(V, get_system(model))

    for component in components
        if _has_variable_cost_parameter(component)
            name = PSY.get_name(component)
            ts_vector = PSY.get_variable_cost(
                component,
                PSY.get_operation_cost(component);
                start_time=initial_forecast_time,
                len=horizon,
            )
            variable_cost_forecast_values = TimeSeries.values(ts_vector)
            for (t, value) in enumerate(variable_cost_forecast_values)
                if attributes.uses_compact_power
                    value, _ = _convert_variable_cost(value)
                end
                _set_param_value!(param_array, PSY.get_cost(value), name, t)
                update_variable_cost!(container, param_array, attributes, component, t)
            end
        end
    end
    return
end

_has_variable_cost_parameter(component::PSY.Component) =
    _has_variable_cost_parameter(PSY.get_operation_cost(component))
_has_variable_cost_parameter(::PSY.MarketBidCost) = true
_has_variable_cost_parameter(::T) where {T <: PSY.OperationalCost} = false

function _update_pwl_cost_expression(
    container::OptimizationContainer,
    ::Type{T},
    component_name::String,
    time_period::Int,
    cost_data::Vector{NTuple{2, Float64}},
) where {T <: PSY.Component}
    pwl_var_container = get_variable(container, PieceWiseLinearCostVariable(), T)
    resolution = get_resolution(container)
    dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR
    gen_cost = JuMP.AffExpr(0.0)
    slopes = PSY.get_slopes(cost_data)
    upb = PSY.get_breakpoint_upperbounds(cost_data)
    for i in 1:length(cost_data)
        JuMP.add_to_expression!(
            gen_cost,
            slopes[i] * upb[i] * dt * pwl_var_container[(component_name, i, time_period)],
        )
    end
    return gen_cost
end

function update_variable_cost!(
    container::OptimizationContainer,
    param_array::JuMPFloatArray,
    attributes::CostFunctionAttributes{Float64},
    component::T,
    time_period::Int,
) where {T <: PSY.Component}
    resolution = get_resolution(container)
    dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR
    base_power = get_base_power(container)
    component_name = PSY.get_name(component)
    cost_data = param_array[component_name, time_period]
    if iszero(cost_data)
        return
    end
    variable = get_variable(container, get_variable_type(attributes)(), T)
    gen_cost = variable[component_name, time_period] * cost_data * base_power * dt
    # Attribute doesn't have multiplier
    # gen_cost = attributes.multiplier * gen_cost_
    add_to_objective_variant_expression!(container, gen_cost)
    set_expression!(container, ProductionCostExpression, gen_cost, component, time_period)
    return
end

function update_variable_cost!(
    container::OptimizationContainer,
    param_array::DenseAxisArray{Vector{NTuple{2, Float64}}},
    ::CostFunctionAttributes{Vector{NTuple{2, Float64}}},
    component::T,
    time_period::Int,
) where {T <: PSY.Component}
    component_name = PSY.get_name(component)
    cost_data = param_array[component_name, time_period]
    if all(iszero.(last.(cost_data)))
        return
    end

    gen_cost =
        _update_pwl_cost_expression(container, T, component_name, time_period, cost_data)
    # Attribute doesn't have multiplier
    # gen_cost = attributes.multiplier * gen_cost_
    add_to_objective_variant_expression!(container, gen_cost)
    set_expression!(container, ProductionCostExpression, gen_cost, component, time_period)
    return
end


text:
"""
    `get_αβ_from_ϕ(ϕ,REG_STRUCT)`

    Get full α and β vectors from vector ϕ=[a,b,c] using Nl, ρ̄l, and dl, assuming energy function
    from Jenkinson et al 2017 and 2018.

        α_l = (a + b⋅ρ̄_l)⋅N_l , β_l = c/d_l

    # Examples
    ```julia-repl
    julia> rs=CpelNano.RegStruct(); rs.L=10; rs.Nl=fill(5.0,rs.L); rs.ρl=fill(0.1,rs.L); rs.dl=fill(10.0,rs.L-1);
    julia> CpelNano.get_αβ_from_ϕ([0.0,0.0,0.0],rs)
    ([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    ```
"""
function get_αβ_from_ϕ(ϕ::Vector{Float64}, rs::RegStruct)::NTuple{2,Vector{Float64}}

    # Obtain full α & β vectors
    α = [rs.Nl[l] * (ϕ[1] + ϕ[2] * rs.ρl[l]) for l = 1:rs.L]
    β = [ϕ[3] / rs.dl[l] for l = 1:(rs.L - 1)]

    # Return α and β vector
    return α, β

end
"""
    `get_α_from_ex(exs)`

    Get α vector from E[X] vector. Only used for marginal model

        α_n = 0.5 * [log(E[X]+1) - log(1-E[X])]

    # Examples
    ```julia-repl
    julia> CpelNano.get_α_from_ex(zeros(5))

    ```
"""
function get_α_from_ex(exs::Vector{Float64})::Vector{Float64}

    # Obtain full α vector
    α = 0.5 * [log((exs[l] + 1) / (1 - exs[l])) for l = 1:length(exs)]
    α = max.(-AMAX, α)
    α = min.(AMAX, α)

    # Return α vector
    return α

end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

module AngularPower
export Lambertian, Cosine, Gaussian

using ....OpticSim
using ...Emitters
using ...Geometry
using LinearAlgebra

abstract type AbstractAngularPowerDistribution{T<:Real} end

"""
    Lambertian{T} <: AbstractAngularPowerDistribution{T}

Ray power is unaffected by angle.
"""
struct Lambertian{T} <: AbstractAngularPowerDistribution{T}
    Lambertian(::Type{T} = Float64) where {T<:Real} = new{T}()
end

Emitters.apply(::Lambertian, ::Transform, power::Real, ::OpticSim.Ray{<:Real,3}) = power

"""
    Cosine{T} <: AbstractAngularPowerDistribution{T}

Cosine power distribution. Ray power is calculated by:

`power = power * (cosine_angle ^ cosine_exp)`
"""
struct Cosine{T} <: AbstractAngularPowerDistribution{T}
    cosine_exp::T

    function Cosine(cosine_exp::T = one(T)) where {T<:Real}
        new{T}(cosine_exp)
    end
end

function Emitters.apply(d::Cosine, tr::Transform, power::Real, ray::OpticSim.Ray{<:Real,3})
    cosanglebetween = dot(OpticSim.direction(ray), forward(tr))
    return power * cosanglebetween ^ d.cosine_exp
end

"""
    Gaussian{T} <: AbstractAngularPowerDistribution{T}

GGaussian power distribution. Ray power is calculated by:

`power = power * exp(-(gaussianu * l^2 + gaussianv * m^2))`
where l and m are the cos_angles between the two axes respectivly.
"""
struct Gaussian{T} <: AbstractAngularPowerDistribution{T}
    gaussianu::T
    gaussianv::T

    function Gaussian(gaussianu::T, gaussianv::T) where {T<:Real}
        new{T}(gaussianu, gaussianv)
    end
end

function Emitters.apply(d::Gaussian, tr::Transform, power::Real, ray::OpticSim.Ray{<:Real,3})
    l = dot(OpticSim.direction(ray), right(tr))
    m = dot(OpticSim.direction(ray), up(tr))
    return power * exp(-(d.gaussianu * l^2 + d.gaussianv * m^2))
end

end # module Angular Power


text:
module TranscodingStreams

export
    TranscodingStream,
    Noop,
    NoopStream

const ByteData = Union{Vector{UInt8},Base.CodeUnits{UInt8}}

include("memory_unsafe.jl")
include("buffer.jl")
include("_data.jl")
include("codec.jl")
include("stream.jl")
include("io_unsafe.jl")
include("noop.jl")
include("transcode.jl")
include("testtools.jl")

end # module


text:
module APMAE4001

using Reexport

include("Matrix.jl")
@reexport using .Matrix

include("FiniteDifference/FiniteDifference.jl")

include("InverseMatrices/InverseMatrices.jl")
@reexport using .InverseMatrices

include("LUFactorization/LUFactorization.jl")
@reexport using .LUFactorization

end # module


text:
iround(x) = round(Int64, x)

triangle(n) = iround.(exp(diagm(-1=> 1:n)))

function pascal(n)
   t=triangle(n)
   println.(join.([filter(!iszero, t[i,:]) for i in 1:(n+1)], " "))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export mkvar, mkpoly

const PolyT = Polynomial{true,Float64}

mkpoly(x::Expr) = eval(symbol_walk(y->mkvar(y), x)) * one(PolyT)
mkpoly(x::Symbol) = mkvar(x) * one(PolyT)
mkpoly(x::Number) = iszero(x) ? zero(PolyT) : one(PolyT) * x
mkpoly(x::T) where {T<:AbstractVariable} = x * one(PolyT)

_varmap = Dict{Symbol, AbstractVariable}()

mkvar(x::String) = mkvar(Symbol(x))
mkvar(x::T) where {T <: AbstractVariable} = x
function mkvar(x::Symbol)
    if haskey(_varmap, x)
        return _varmap[x]
    end
    v = PolyVar{true}(string(x))
    push!(_varmap, x=>v)
    v
end

MultivariatePolynomials.isconstant(p::AbstractPolynomialLike) = all(isconstant, terms(p))


text:
struct Callbacks
    cbs::Vector
    executor::CallbackExecutor
    graph::SimpleDiGraph
end
Callbacks(cbs, executor = LinearExecutor()) = Callbacks(cbs, executor, callbackgraph(cbs))

getmetrics(cbs::Callbacks) = [cb for cb in cbs.cbs if cb isa AbstractMetric]
getloss(cbs::Callbacks) = only([cb for cb in cbs.cbs if cb isa AverageLoss])

"""
    $TYPEDEF

Stores data of the last processed batch.

$FIELDS

(!) If used in callbacks, some fields may be `nothing` as
they are reset after every step.
"""
# TODO: maybe make this dependent on the current `Phase`?
@with_kw mutable struct BatchState
    xs = nothing
    ys = nothing
    ŷs = nothing
    loss = nothing
    grads = nothing
end

"""
    Learner(model, data, opt, lossfn)

Central object for training that holds all necessary state.

$(TYPEDFIELDS)

# Fields
- `device`: device to train on, usually `Flux.cpu` or `Flux.gpu`
- `params`: model parameters, i.e. `Flux.params(model)`
- `phase::`[`AbstractFittingPhase`](@ref): current fitting phase
- `metrics::AbstractVector{`[`AbstractMetric`](@ref)`}`: metric callbacks to
    evaluate performance during training
- `config`
- `callbacks::AbstractVector{`[`AbstractCallback`](@ref)`}`: callbacks for training loop
- `recorder::`[`Recorder`](@ref): special callback that records metrics and hyperparameters
- `scheduler::`[`Scheduler`](@ref): special callback that records metrics and hyperparameters
- `batch::`[`BatchState`](@ref): Holds training state during batch
"""
@with_kw mutable struct Learner
    model
    data
    opt
    lossfn
    params
    batch::BatchState
    callbacks::Callbacks
    cbstate::Dict
end

function Learner(
        model, data, opt, lossfn;
        callbacks = [], metrics = [], schedule = Schedules(),
        usedefaultcallbacks = true, config = Dict(), cbexecutor = LinearExecutor()
    )
    if usedefaultcallbacks
        callbacks = vcat(defaultcallbacks(), callbacks)
    end
    callbacks = vcat(callbacks, metrics)
    Callbacks(callbacks)

    return Learner(
        model, data, opt, lossfn,
        Flux.params(model),
        BatchState(),
        Callbacks(callbacks, cbexecutor),
        Dict{Symbol, Any}())
end

defaultcallbacks()::Vector{AbstractCallback} = [
    ProgressBarLogger(),
    MetricsLogger(),
    StopOnNaNLoss(),
    Recorder(),
    AverageLoss(),
]

#  Callback handling
handle(event, learner, phase) = handle(learner.callbacks.executor, event, phase, learner)

# Other

getdataloader(phase::AbstractTrainingPhase, learner) = learner.data[1]
getdataloader(phase::ValidationPhase, learner) = learner.data[2]

function model!(learner, model)
    learner.model = model
    learner.params = Flux.params(model)
end

numsteps(learner::Protected, phase) = numsteps(getfield(learner, :data), phase)
numsteps(learner, phase) = length(getdataloader(phase, learner))

# TODO: fix

#=
function setschedule!(learner, schedule)
    learner.scheduler = ParamScheduler(
        delayschedule(schedule, learner.recorder.epoch)
    )
end
=#

# TODO: fix
#=
function numsteps(
        learner::Learner,
        phase::AbstractFittingPhase = learner.state.phase)
    return length(getdataloader(phase, learner))
end
numsteps(p::Protected{Learner}, phase) = numsteps(getfield(p, :data), phase)
=#

# TODO: move to Callback

#=
function artifactpath(learner)
    p = joinpath(pwd(), ".artifacts", string(learner.runid))
    if !isdir(p)
        mkpath(p)
    end
    return p
end
=#