text:
"""
`PairwiseListMatrix{T, diagonal, VT}` is a (squared) symmetric matrix that stores a `list`
of type `VT` with values of type `T` for the pairwise comparison/evaluation of `nelements`.
If `diagonal` is `true` the first element of the list is `1, 1` otherwise is `1, 2`.
If `diagonal` is `false` the diagonal values are stored in a vector on the `diag` field.
"""
mutable struct PairwiseListMatrix{T,diagonal,VT} <: AbstractArray{T, 2}
    list::VT
    diag::VT
    nelements::Int

    function PairwiseListMatrix{T,diagonal,VT}(list::VT,
                                               diag::VT,
                                               N::Int) where {T,diagonal,VT}
        @assert eltype(VT) === T "Field vectors must have the same element type than the PairwiseListMatrix"
        if !diagonal
            n_diag = length(diag)
            n_list = _nelements(length(list))
            @assert n_diag == N "There are $n_diag elements in the diag field instead of $N"
        else
            n_list = _nelements_with_diagonal(length(list))
            @assert n_list == N "There are $n_list elements in the list field instead of $N"
        end
        new(list, diag, N)
    end
end

"Returns `true` if the list has diagonal values."
@inline hasdiagonal(::PairwiseListMatrix{T,diagonal,VT}) where {T,diagonal,VT} = diagonal
@inline hasdiagonal(nplm::NamedArray{T,2,PairwiseListMatrix{T,D,TV},DN}) where {T,D,TV,DN} = D

"Retuns the `list` vector."
@inline getlist(plm::PairwiseListMatrix) = plm.list

"Retuns the `diag` vector (which contains the diagonal values if `diagonal` is `false`)."
@inline getdiag(plm::PairwiseListMatrix) = plm.diag

# listlength
# ==========

"Returns the length of the `list` field"
lengthlist(plm::PairwiseListMatrix) = length(getlist(plm))

function lengthlist(nelements::Int, ::Type{Val{diagonal}}) where diagonal
    diagonal ? div(nelements*(nelements+1),2) : div(nelements*(nelements-1),2)
end

"""
Returns the list length needed for a pairwise measures or comparisons of `nelements`.
If `diagonal` is `true`, diagonal values are included in the list.

```jldoctest
julia> using PairwiseListMatrices

julia> plm = PairwiseListMatrix([1, 2, 3, 4, 5, 6], false)
4×4 PairwiseListMatrix{Int64,false,Array{Int64,1}}:
 0  1  2  3
 1  0  4  5
 2  4  0  6
 3  5  6  0

julia> lengthlist(4, false)
6

julia> lengthlist(plm)
6

```
"""
lengthlist(nelements::Int, diagonal::Bool) = lengthlist(nelements, Val{diagonal})

# Indexes: ij2k
# =============

@inline @fastmath function ij2k(i::Int, j::Int, nelements::Int,
                                ::Type{Val{diagonal}}) where diagonal
    div(diagonal ?  (nelements*(nelements+1))-((nelements-i)*(nelements-i+1)) :
    (nelements*(nelements-1))-((nelements-i)*(nelements-i-1)), 2) - nelements + j
end

"""
Returns the `k` index of the `list` from the indixes `i` and `j` with `i<j` from a matrix
of `nelements` by `nelements`. `diagonal` should be `true` or `Val{true}` if the diagonal
values are on the `list`. You must not use it with `i>j`.

```jldoctest
julia> using PairwiseListMatrices

julia> plm = PairwiseListMatrix([10,20,30,40,50,60], true)
3×3 PairwiseListMatrix{Int64,true,Array{Int64,1}}:
 10  20  30
 20  40  50
 30  50  60

julia> ij2k(1, 2, 3, true)
2

julia> getlist(plm)[2]
20

```
"""
@inline function ij2k(i::Int, j::Int, nelements::Int, diagonal::Bool)
    ij2k(i, j, nelements, Val{diagonal})
end

# Empties
# -------

function _diagonal_value(::Type{T}) where T
    if hasmethod(zero, (T,))
        return zero(T)
    else
        if T isa Union
            T_types = Base.uniontypes(T)
            if Missing ∈ T_types
                return missing
            end
            if Nothing ∈ T_types
                return nothing
            end
        elseif T === Any
            return nothing
        end
    end
    throw(ErrorException(
       "Please use the last argument to fill the diagonal. It should be of type $T."))
end

"""
An empty `PairwiseListMatrix` can be created for a given `Type` and a number of elements
`nelements`. The `diagonal` (default to `false`) can be declared as `true` to indicate that
the list needs space for the diagonal elements. If `diagonal` is `false` the diagonal
values are stored in a vector on the `diag` field instead of being on the list.
The `diag` vector can be filled with the optional `diagonalvalue` argument (default to `0`).

```julia
PairwiseListMatrix(Int, 3)
PairwiseListMatrix(Int, 3, true)
```
"""
function PairwiseListMatrix(::Type{T},
                           nelements::Int,
                           diagonal::Bool = false,
                           diagonalvalue::T = _diagonal_value(T)) where T
    if diagonal
        return PairwiseListMatrix{T, true, Vector{T}}(
            Array{T}(undef, lengthlist(nelements, Val{true})),
            T[],
            nelements)
    else
        return PairwiseListMatrix{T, false, Vector{T}}(
            Array{T}(undef, lengthlist(nelements, Val{false})),
            fill!(Array{T}(undef, nelements), diagonalvalue),
            nelements)
    end
end

"""
```julia
PairwiseListMatrix{Int,false,Vector{Int}}(3)
```
"""
function PairwiseListMatrix{T,diagonal,VT}(nelements::Int) where {T,diagonal,VT}
    if diagonal
        return PairwiseListMatrix{T,true,VT}(
            convert(VT, Array{T}(undef, lengthlist(nelements, Val{true}))),
            convert(VT, T[]),
            nelements)
    else
        return PairwiseListMatrix{T, false, VT}(
            Array{T}(undef, lengthlist(nelements, Val{false})),
            convert(VT, Array{T}(undef, nelements)),
            nelements)
    end
end

# Convert
# -------

function Base.convert(::Type{PairwiseListMatrix{T,D,VT}},
                     mat::PairwiseListMatrix{F,D,VF}) where {T,F,D,VT,VF}
    PairwiseListMatrix{T,D,VT}(convert(VT, mat.list),
        convert(VT, mat.diag),
        copy(mat.nelements))
end

function Base.convert(::Type{PairwiseListMatrix{T,true,VT}},
                     mat::PairwiseListMatrix{F,false,VF}) where {T,F,VT,VF}
    N = copy(mat.nelements)
    list = convert(VT, Array{T}(undef, lengthlist(N, Val{true})))
    k = 0
    @inbounds for i in 1:N
        for j in i:N
            list[k += 1] = mat[i,j]
        end
    end
    PairwiseListMatrix{T, true, VT}(list, convert(VT,T[]), N)
end

function Base.convert(::Type{PairwiseListMatrix{T, false, VT}},
                     mat::PairwiseListMatrix{F, true, VF}) where {T,F,VT,VF}
    N = copy(mat.nelements)
    diag = convert(VT, Array{T}(undef, N))
    matlist = mat.list
    list = convert(VT, Array{T}(undef, length(matlist) - N))
    l = 0
    d = 0
    m = 0
    @inbounds for i in 1:N
        for j in i:N
            if i != j
                @inbounds list[l += 1] = matlist[m += 1]
            else
                @inbounds diag[d += 1] = matlist[m += 1]
            end
        end
    end
    PairwiseListMatrix{T, false, VT}(list, diag, N)
end

function Base.convert(::Type{PairwiseListMatrix{T,D,VT}}, mat::Matrix{F}) where {T,F,D,VT}
    N = size(mat, 1)
    if N != size(mat, 2)
        throw(ErrorException("It should be a squared matrix."))
    end
    plm = PairwiseListMatrix{T,D,VT}(N)
    @inbounds for i in 1:N
        for j in i:N
            value = mat[i,j]
            if i != j && value != mat[j,i]
                throw(ErrorException("It should be a symmetric matrix."))
            end
            plm[i,j] = value
        end
    end
    plm
end

# This is faster than list comprehension (2.4 x)
function Base.convert(::Type{Array{F,2}}, lm::PairwiseListMatrix{T, true, VT}) where {T,F,VT}
    N = lm.nelements
    complete = Array{F}(undef, N, N)
    list = lm.list
    k = 0
    @inbounds for col in 1:N
        for row in col:N
            complete[row, col] = complete[col, row] = F(list[k += 1])
        end
    end
    complete
end

function Base.convert(::Type{Array{F,2}}, lm::PairwiseListMatrix{T, false, VT}) where {T,F,VT}
    N = lm.nelements
    complete = Array{F}(undef, N, N)
    list = lm.list
    diag = lm.diag
    k = 0
    l = 0
    @inbounds for col in 1:(N-1)
        for row in (col+1):N
            complete[row, col] = complete[col, row] = F(list[k += 1])
        end
    end
    @inbounds @simd for i in 1:N
        complete[i, i] = F(diag[i])
    end
    complete
end

# Matrix(m::Symmetric{T, PairwiseListMatrix{T,D,VT}}) where {T,D,VT} = Matrix(m.data)

# From a list
# -----------

@inline _nelements(len::Int) = @fastmath div(1+Int(sqrt(1+8*len)),2)
@inline _nelements_with_diagonal(len::Int) = @fastmath div(Int(sqrt(1+8*len)-1),2)

"""
A `PairwiseListMatrix` can be created from a `list`. The `diagonal` (default to `false`)
could be declared as `true` to indicate that the list has the diagonal elements. If
`diagonal` is `false`, the diagonal values are stored in a vector on the `diag` field
instead of being on the list. The `diag` vector can be filled with the optional
`diagonalvalue` argument (default to `0`).
"""
function PairwiseListMatrix(list::AbstractVector{T},
                           diagonal::Bool = false,
                           diagonalvalue::T = _diagonal_value(T)) where T
    VT = typeof(list)
    if diagonal
        nelements = _nelements_with_diagonal(length(list))
        return PairwiseListMatrix{T, true, VT}(list, convert(VT,T[]), nelements)
    else
        nelements = _nelements(length(list))
        return PairwiseListMatrix{T, false, VT}(list,
            fill!(convert(VT, Array{T}(undef, nelements)), diagonalvalue),
            nelements)
    end
end

# AbstractArray methods
# =====================

Base.size(m::PairwiseListMatrix) = (m.nelements, m.nelements)
Base.length(m::PairwiseListMatrix) = m.nelements * m.nelements

Base.eltype(m::PairwiseListMatrix{T, diagonal, VT}) where {T, diagonal, VT} = T

for F in [:similar, :copy, :zeros, :ones]
    @eval begin
        function Base.$F(m::PairwiseListMatrix{T, diagonal, VT}) where {T, diagonal, VT}
            list = $F(m.list)
            VOUT = typeof(list)
            diag = convert(VOUT, $F(m.diag))
            PairwiseListMatrix{T, diagonal, VOUT}(list, diag, copy(m.nelements))
        end
    end
end

function Base.similar(m::PairwiseListMatrix{T, diagonal, VT},
                      ::Type{S}) where {T, diagonal, VT, S}
    list = similar(m.list, S)
    VOUT = typeof(list)
    diag = convert(VOUT, similar(m.diag, S))
    PairwiseListMatrix{S, diagonal, VOUT}(list, diag, copy(m.nelements))
end

# Indexing (getindex)
# ===================

function Base.getindex(lm::PairwiseListMatrix{T, true, VT}, i::Int, j::Int) where {T, VT}
    if i <= j
        return(lm.list[ij2k(i, j, lm.nelements, Val{true})])
    else
        return(lm.list[ij2k(j, i, lm.nelements, Val{true})])
    end
end

function Base.getindex(lm::PairwiseListMatrix{T, false, VT}, i::Int, j::Int) where {T, VT}
    if i < j
        return(lm.list[ij2k(i, j, lm.nelements, Val{false})])
    elseif i > j
        return(lm.list[ij2k(j, i, lm.nelements, Val{false})])
    else
        return(lm.diag[i])
    end
end

function Base.getindex(lm::PairwiseListMatrix{T, true, VT}, i::Int) where {T, VT}
    n = lm.nelements
    row = Int(ceil(i/n))
    col = i - n*(row-1)
    if row <= col
        @inbounds return lm.list[ij2k(row, col, n, Val{true})]
    else
        @inbounds return lm.list[ij2k(col, row, n, Val{true})]
    end
end

function Base.getindex(lm::PairwiseListMatrix{T, false, VT}, i::Int) where {T, VT}
    n = lm.nelements
    row = Int(ceil(i/n))
    col = i - n*(row-1)
    if row < col
        @inbounds return lm.list[ij2k(row, col, n, Val{false})]
    elseif row > col
        @inbounds return lm.list[ij2k(col, row, n, Val{false})]
    else
        @inbounds return lm.diag[row]
    end
end

# Set values (setindex!)
# ======================

function Base.setindex!(lm::PairwiseListMatrix{T, true, VT}, v, i::Int, j::Int) where {T, VT}
    if i <= j
        return setindex!(lm.list, v, ij2k(i, j, lm.nelements, Val{true}))
    else
        return setindex!(lm.list, v, ij2k(j, i, lm.nelements, Val{true}))
    end
end

function Base.setindex!(lm::PairwiseListMatrix{T, false, VT}, v, i::Int, j::Int) where {T, VT}
    if i < j
        return setindex!(lm.list, v, ij2k(i, j, lm.nelements, Val{false}))
    elseif i > j
        return setindex!(lm.list, v, ij2k(j, i, lm.nelements, Val{false}))
    else
        return setindex!(lm.diag, v, i)
    end
end

function Base.setindex!(lm::PairwiseListMatrix{T, true, VT}, v, i::Int) where {T, VT}
    n = lm.nelements
    row = Int(ceil(i/n))
    col = i - n*(row-1)
    if row <= col
        return setindex!(lm.list, v, ij2k(row, col, n, Val{true}))
    else
        return setindex!(lm.list, v, ij2k(col, row, n, Val{true}))
    end
end

function Base.setindex!(lm::PairwiseListMatrix{T, false, VT}, v, i::Int) where {T, VT}
    n = lm.nelements
    row = Int(ceil(i/n))
    col = i - n*(row-1)
    if row < col
        return setindex!(lm.list, v, ij2k(row, col, n, Val{false}))
    elseif row > col
        return setindex!(lm.list, v, ij2k(col, row, n, Val{false}))
    else
        return setindex!(lm.diag, v, row)
    end
end

# Transpose
# =========

LinearAlgebra.transpose(lm::PairwiseListMatrix)  = lm
LinearAlgebra.transpose!(lm::PairwiseListMatrix) = lm

LinearAlgebra.adjoint(lm::PairwiseListMatrix)  = lm
LinearAlgebra.adjoint!(lm::PairwiseListMatrix) = lm

# lu
# ==

function LinearAlgebra.lu!(lm::PairwiseListMatrix, args...; kargs...)
    LinearAlgebra.lu!(Matrix(lm), args...; kargs...)
end

# Diagonal
# ========

"""
Returns a vector of type `VT` from a `PairwiseListMatrix{T, false, VT}` that has
the diagonal values.

```jldoctest
julia> using PairwiseListMatrices

julia> plm = PairwiseListMatrix([10,20,30,40,50,60], true)
3×3 PairwiseListMatrix{Int64,true,Array{Int64,1}}:
 10  20  30
 20  40  50
 30  50  60

julia> diagonal(plm)
3-element Array{Int64,1}:
 10
 40
 60

```
"""
diagonal(lm::PairwiseListMatrix{T, false, VT}) where {T,VT} = lm.diag

function diagonal(lm::PairwiseListMatrix{T, true, VT}) where {T,VT}
    convert(VT, T[ lm[i,i] for i in 1:lm.nelements ])
end

# Customizing broadcasting
# ========================

function Base.BroadcastStyle(::Type{PairwiseListMatrix{T, D, VT}}) where {T, D, VT}
    Broadcast.DefaultArrayStyle{2}()
end

function Base.Broadcast.broadcasted(::Broadcast.DefaultArrayStyle{2},
                                    fun::F,
                                    plm::PairwiseListMatrix{T,D,V}) where {F<:Function,T,D,V}
    list = fun.(plm.list)
    E = eltype(list)
    if D
        diag = Array{E, 1}()
    else
        _diag = fun.(plm.diag)
        if eltype(_diag) !== E
            diag = copyto!(similar(_diag, E), _diag)
        else
            diag = _diag
        end
    end
    PairwiseListMatrix{E, D, typeof(list)}(list, diag, plm.nelements)
end

# Slow for multiple operations (it doesn't fuse them), e.g.: 0.5 .* $plm ./ 10.0
function Base.Broadcast.broadcasted(::Broadcast.DefaultArrayStyle{2},
                                    fun::F,
                                    plm::PairwiseListMatrix{T,D,V},
                                    n::Number) where {F<:Function,T,D,V}
    list = broadcast(fun, plm.list, n)
    E = eltype(list)
    if D
        diag = Array{E, 1}()
    else
        _diag = broadcast(fun, plm.diag, n)
        if eltype(_diag) !== E
            diag = copyto!(similar(_diag, E), _diag)
        else
            diag = _diag
        end
    end
    PairwiseListMatrix{E, D, typeof(list)}(list, diag, plm.nelements)
end

function Base.Broadcast.broadcasted(::Broadcast.DefaultArrayStyle{2},
                                    fun::F,
                                    n::Number,
                                    plm::PairwiseListMatrix{T,D,V}) where {F<:Function,T,D,V}
    list = broadcast(fun, n, plm.list)
    E = eltype(list)
    if D
        diag = Array{E, 1}()
    else
        _diag = broadcast(fun, n, plm.diag)
        if eltype(_diag) !== E
            diag = copyto!(similar(_diag, E), _diag)
        else
            diag = _diag
        end
    end
    PairwiseListMatrix{E, D, typeof(list)}(list, diag, plm.nelements)
end

@inline function Base.copyto!(dest::PairwiseListMatrix{T, D, VT},
                              bc::Base.Broadcast.Broadcasted{Nothing}) where {T, D, VT}
    axes(dest) == axes(bc) || Base.Broadcast.throwdm(axes(dest), axes(bc))
    bc_ = Base.Broadcast.preprocess(dest, bc)
    @iterateupper dest true list[k] = bc_[CartesianIndex(i,j)] # slow: bc_ has a plm
    return dest
end

# Unary operations (faster than default methods)
# ==============================================

for F in (:abs, :-)
    @eval begin
        function Base.$F(lm::PairwiseListMatrix{T,true,VT}) where {T,VT}
            list = $F(lm.list)
            VOUT = typeof(list)
            diag = convert(VOUT, lm.diag)
            PairwiseListMatrix{eltype(list), true, VOUT}(list, diag, lm.nelements)
        end
        function Base.$F(lm::PairwiseListMatrix{T,false,VT}) where {T,VT}
            list = $F(lm.list)
            VOUT = typeof(list)
            diag = convert(VOUT, $F(lm.diag))
            PairwiseListMatrix{eltype(list), false, VOUT}(list, diag, lm.nelements)
        end
    end
end

for F in (:map, ) # , :broadcast)
    @eval begin
        function Base.$F(f, plm::PairwiseListMatrix{T,true,VT}) where {T,VT}
            list = $F(f, plm.list)
            VOUT = typeof(list)
            diag = convert(VOUT, plm.diag)
            PairwiseListMatrix{eltype(list), true, VOUT}(list, diag, plm.nelements)
        end

        function Base.$F(f, plm::PairwiseListMatrix{T,false,VT}) where {T,VT}
            list = $F(f, plm.list)
            VOUT = typeof(list)
            diag = convert(VOUT, $F(f, plm.diag))
            PairwiseListMatrix{eltype(list), false, VOUT}(list, diag, plm.nelements)
        end
    end
end

for F in (:map!, ) # , :broadcast!)
    @eval begin
        function Base.$F(f, plm::PairwiseListMatrix{T,true,VT}) where {T,VT}
            $F(f, plm.list)
            plm
        end

        function Base.$F(f, plm::PairwiseListMatrix{T,false,VT}) where {T,VT}
            $F(f, plm.list)
            $F(f, plm.diag)
            plm
        end
    end
end

LinearAlgebra.svd(m::PairwiseListMatrix; kargs...) = LinearAlgebra.svd(Matrix(m), kargs...)

function LinearAlgebra.svd(m::UpperTriangular{T,PairwiseListMatrix{T,D,VT}}) where {T,D,VT}
    svd(Matrix(m))
end

function LinearAlgebra.svd(m::LowerTriangular{T,PairwiseListMatrix{T,D,VT}}) where {T,D,VT}
    svd(Matrix(m))
end

# Faster mean and sum
# ===================

function _sum_kernel!(sum_i, list, N)
    k = 0
    l = 0
    for i in 1:N
        l += 1
        @inbounds for j in i:N
            sum_i[i] += list[k += 1]
        end
        @inbounds for j in (i+1):N
            sum_i[j] += list[l += 1]
        end
    end
    sum_i
end

function _sum_kernel!(sum_i, diag, list, N)
    k = 0
    l = 0
    for i in 1:(N-1)
        @inbounds for j in (i+1):N
            sum_i[i] += list[k += 1]
        end
        @inbounds for j in (i+1):N
            sum_i[j] += list[l += 1]
        end
    end
    @inbounds @simd for i in 1:N
        sum_i[i] += diag[i]
    end
    sum_i
end

function _test_and_setsum(::Type{T}, dims, N) where T
    if dims == 1
        sum_i = zeros(T, 1, N)
    elseif dims == 2
        sum_i = zeros(T, N, 1)
    else
        throw(ErrorException("dims should be 1 or 2"))
    end
    sum_i
end

_sum(m::PairwiseListMatrix{T, false}) where {T} = T(2)*sum(m.list) + sum(m.diag)
_sum(m::PairwiseListMatrix{T, true}) where {T} = T(2)*sum(m.list) - sum(diagonal(m))

function Base.sum(lm::PairwiseListMatrix{T,true,VT}; dims::Union{Int,Colon}=:) where {T,VT}
    if isa(dims, Colon)
        return _sum(lm)
    end
    N = lm.nelements
    sum_i = _test_and_setsum(T, dims, N)
    _sum_kernel!(sum_i, lm.list, N)
end

function Base.sum(lm::PairwiseListMatrix{T,false,VT};
                  dims::Union{Int,Colon}=:) where {T, VT}
    if isa(dims, Colon)
        return _sum(lm)
    end
    N = lm.nelements
    sum_i = _test_and_setsum(T, dims, N)
    _sum_kernel!(sum_i, lm.diag, lm.list, N)
end

function Statistics.mean(m::PairwiseListMatrix; dims::Union{Int,Colon}=:)
    if isa(dims, Colon)
        sum(m)/length(m)
    else
        sum(m, dims=dims) ./ m.nelements
    end
end

# Sum/Mean without diagonal: sum/mean_nodiag
# ------------------------------------------

function _sum_nodiag_kernel!(sum_i, lm::PairwiseListMatrix{T,true,VT}, N) where {T,VT}
    list = lm.list
    k = 0
    @inbounds for i in 1:N
        for j in i:N
            k += 1
            if i != j
                value = list[k]
                sum_i[i] += value
                sum_i[j] += value
            end
        end
    end
    sum_i
end

function _sum_nodiag_kernel!(sum_i, lm::PairwiseListMatrix{T,false,VT}, N) where {T,VT}
    list = lm.list
    k = 0
    l = 0
    for i in 1:(N-1)
        @inbounds for j in (i+1):N
            sum_i[i] += list[k += 1]
        end
        @inbounds for j in (i+1):N
            sum_i[j] += list[l += 1]
        end
    end
    sum_i
end

_sum_nodiag(m::PairwiseListMatrix{T,false,VT}) where {T,VT} = T(2) * sum(m.list)
function _sum_nodiag(m::PairwiseListMatrix{T,true,VT}) where {T,VT}
    T(2) * sum(m.list) - sum(diagonal(m))
end

"Sum the values outside the diagonal"
function sum_nodiag(lm::PairwiseListMatrix{T,diagonal,VT};
                    dims::Union{Int,Colon}=:) where {T,diagonal,VT}
    if isa(dims, Colon)
        return _sum_nodiag(lm)
    end
    N = lm.nelements
    sum_i = _test_and_setsum(T, dims, N)
    _sum_nodiag_kernel!(sum_i, lm, N)
end

"Mean of the values outside the diagonal"
function mean_nodiag(m::PairwiseListMatrix; dims::Union{Int,Colon}=:)
    if isa(dims, Colon)
        return sum_nodiag(m) / (length(m) - m.nelements)
    end
    sum_nodiag(m, dims=dims) ./ (m.nelements-1)
end

# Operations on Vector{PairwiseListMatrix}
# ========================================

# sum
# ---

@inline _has_diagonal(x::PairwiseListMatrix{T,diagonal,VT}) where {T,diagonal,VT} = diagonal

function Base.sum(list::Vector{PairwiseListMatrix{T,diagonal,VT}}) where {T,diagonal,VT}
    samples = length(list)
    if samples == 0
        throw(ErrorException("Empty list"))
    end
    if samples == 1
        return list[1]
    else
        start = copy(list[1])
        i = 2
        start_list = start.list
        N = length(start_list)
        i = 2
        while i <= samples
            @inbounds ylist = list[i].list
            if length(ylist) != N
                throw(ErrorException("Different number of elements"))
            end
            @inbounds @simd for k in 1:N
                start_list[k] += ylist[k]
            end
            i += 1
        end
        if !diagonal
            i = 2
            start_diag = start.diag
            while i <= samples
                @inbounds ydiag = list[i].diag
                @inbounds @simd for k in 1:length(start_diag)
                    start_diag[k] += ydiag[k]
                end
                i += 1
            end
        end
        return start
    end
end

# mean
# ----

function Statistics.mean(list::Vector{PairwiseListMatrix{T,diagonal,VT}}) where {T,diagonal,VT}
    sum(list) ./ length(list)
end

# std
# ---

function Statistics.varm(list::Vector{PairwiseListMatrix{T,D,VT}},
                   mean::PairwiseListMatrix{M,D,VM}) where {T,D,VT,M,VM}
    samples = length(list)
    if samples < 2
        throw(ErrorException("You need at least 2 samples."))
    end
    mean_list = mean.list
    N = length(mean_list)
    OutType = promote_type(T, M)
    out = fill!(PairwiseListMatrix{OutType,D,Vector{OutType}}(mean.nelements),zero(OutType))
    out_list = out.list
    @inbounds for sample in list
        sample_list = sample.list
        if length(sample_list) != N
            throw(ErrorException("Different number of elements"))
        end
        @inbounds @simd for k in 1:N
            out_list[k] += abs2( sample_list[k] - mean_list[k])
        end
    end
    if !D
        out_diag = out.diag
        mean_diag = mean.diag
        @inbounds for sample in list
            sample_diag = sample.diag
            for k in 1:length(out_diag)
                @inbounds out_diag[k] += abs2( sample_diag[k] - mean_diag[k])
            end
        end
    end
    out ./ samples
end

function Statistics.var(list::Vector{PairwiseListMatrix{T,D,VT}}; mean=nothing) where {T,D,VT}
    if mean === nothing
        varm(list, Statistics.mean(list))
    else
        varm(list, mean)
    end
end

function Statistics.std(list::Vector{PairwiseListMatrix{T,D,VT}}; mean=nothing) where {T,D,VT}
    sqrt.(var(list, mean=mean))
end

# zscore
# ======

"""
This function takes a vector of `PairwiseListMatrix` objects and a `PairwiseListMatrix` and
fill the matrix with the zscore value using the median and std of the vector.
"""
function zscore!(list::Vector{PairwiseListMatrix{L,D,VL}},
                 mat::PairwiseListMatrix{E,D,VE}) where {L,D,VL,E,VE}
    list_mean = mean(list)
    if size(mat) != size(list_mean)
        throw(ErrorException("PairwiseListMatrices must have the same size"))
    end
    list_std  = std(list, mean=list_mean)

    mat_list = mat.list
    std_list = list_std.list
    mean_list = list_mean.list

    @inbounds for i in 1:length(mat_list)
        s = std_list[i]
        m = mean_list[i]
        value = mat_list[i]
        if (m - value + one(E)) ≉ one(E)
            if (s + one(E)) ≉ one(E)
                mat_list[i] = (value - m)/s
            else
                mat_list[i] = NaN
            end
        else
            mat_list[i] = zero(E)
        end
    end

    if D
        mat_diag = mat.diag
        std_diag = list_std.diag
        mean_diag = list_mean.diag
        @inbounds for i in 1:length(mat_diag)
            s = std_diag[i]
            m = mean_diag[i]
            value = mat_diag[i]
            if (m - value + one(E)) ≉ one(E)
                if (s + one(E)) ≉ one(E)
                    mat_diag[i] = (value - m)/s
                else
                    mat_diag[i] = NaN
                end
            else
                mat_diag[i] = zero(E)
            end
        end
    end

    mat
end

"It's like `zscore!` but without modifying the `PairwiseListMatrix`"
zscore(list, mat) = zscore!(list, copy(mat))

# NamedArrays
# ===========

# labels
# ------

"""
It gets the labels of a PairwiseListMatrix.

```jldoctest
julia> using PairwiseListMatrices

julia> plm  = PairwiseListMatrix(ones(3), false)
3×3 PairwiseListMatrix{Float64,false,Array{Float64,1}}:
 0.0  1.0  1.0
 1.0  0.0  1.0
 1.0  1.0  0.0

julia> getlabels(plm)
3-element Array{String,1}:
 "1"
 "2"
 "3"

julia> nplm  = setlabels(plm, ["a","b","c"])
3×3 Named PairwiseListMatrix{Float64,false,Array{Float64,1}}
A ╲ B │   a    b    c
──────┼──────────────
a     │ 0.0  1.0  1.0
b     │ 1.0  0.0  1.0
c     │ 1.0  1.0  0.0

julia> getlabels(nplm)
3-element Array{String,1}:
 "a"
 "b"
 "c"

```
"""
function getlabels(plm::PairwiseListMatrix{T,D,TV}) where {T,D,TV}
    String[ string(i) for i in 1:(plm.nelements) ]
end

getlabels(nplm::NamedArray{T,2,PairwiseListMatrix{T,D,TV},DN}) where {T,D,TV,DN} = names(nplm)[1]

"It changes the labels of a Named PairwiseListMatrix"
function setlabels!(nplm::NamedArray{T,2,PairwiseListMatrix{T,D,TV},DN},
                    labels::Vector{String}) where {T,D,TV,DN}
    NE = nplm.array.nelements
    NL = length(labels)
    @assert NL == NE "The number of elements, $NE, must be equal to the number of labels, $NL"
    setnames!(nplm, labels, 1)
    setnames!(nplm, labels, 2)
    nplm
end

"""
Creates a Named PairwiseListMatrix.

```jldoctest
julia> using PairwiseListMatrices

julia> plm  = PairwiseListMatrix(ones(3), false)
3×3 PairwiseListMatrix{Float64,false,Array{Float64,1}}:
 0.0  1.0  1.0
 1.0  0.0  1.0
 1.0  1.0  0.0

julia> nplm  = setlabels(plm, ["a","b","c"])
3×3 Named PairwiseListMatrix{Float64,false,Array{Float64,1}}
A ╲ B │   a    b    c
──────┼──────────────
a     │ 0.0  1.0  1.0
b     │ 1.0  0.0  1.0
c     │ 1.0  1.0  0.0

```
"""
function setlabels(nplm::NamedArray{T,2,PairwiseListMatrix{T,D,TV},DN},
                   labels::Vector{String}) where {T,D,TV,DN}
    setlabels!(copy(nplm), labels)
end

function setlabels(plm::PairwiseListMatrix{T,D,TV}, labels::Vector{String}) where {T,D,TV}
    nplm = NamedArray(plm)
    setlabels!(nplm, labels)
end

# Delegate to PairwiseListMatrix
# ------------------------------

for F in (:getlist, :getdiag, :lengthlist, :sum_nodiag, :mean_nodiag,
          :diagonal)
    @eval begin
        function $F(m::NamedArray{T,2,PairwiseListMatrix{T,D,TV},DN}) where {T,D,TV,DN}
            $F(m.array)
        end
    end
end

# Keep names
for F in (:map, :map!, :broadcast, :broadcast!)
    @eval begin
        function Base.$F(f, m::NamedArray{T,2,PairwiseListMatrix{T,D,TV},DN}) where {T,D,TV,DN}
            NamedArray(($F)(f, m.array), m.dicts, m.dimnames)
        end
    end
end

# Tables
# ======

"""
Creates a `Matrix{Any}`, labels are stored in the columns 1 and 2, and the
values in the column 3. Diagonal values are included by default.

```jldoctest
julia> using PairwiseListMatrices

julia> plm = PairwiseListMatrix([10,20,30], false)
3×3 PairwiseListMatrix{Int64,false,Array{Int64,1}}:
  0  10  20
 10   0  30
 20  30   0

julia> to_table(plm)
6×3 Array{Any,2}:
 "1"  "1"   0
 "1"  "2"  10
 "1"  "3"  20
 "2"  "2"   0
 "2"  "3"  30
 "3"  "3"   0

julia> to_table(plm, diagonal=false)
3×3 Array{Any,2}:
 "1"  "2"  10
 "1"  "3"  20
 "2"  "3"  30

```
"""
function to_table(plm::PairwiseListMatrix; diagonal::Bool = true, labels = getlabels(plm))
    N = plm.nelements
    table = Array{Any}(undef, diagonal ? div(N*(N+1),2) : div(N*(N-1),2), 3)
    t = 0
    @iterateupper plm diagonal begin
        t += 1
        table[t, 1] = labels[i]
        table[t, 2] = labels[j]
        table[t, 3] = list[k]
    end
    table
end

function to_table(nplm::NamedArray{T,2,PairwiseListMatrix{T,D,TV},DN};
                 diagonal::Bool = true,
                 labels::Vector{String} = getlabels(nplm)) where {T,D,TV,DN}
    to_table(nplm.array, diagonal=diagonal, labels=labels)
end

"""
It takes a `PairwiseListMatrix` and converts it to a `Dict` of `Symbol`s to arrays.
The returned dictionary can be easily converted into a `DataFrame`.

```jldoctest
julia> using PairwiseListMatrices, DataFrames

julia> nplm = setlabels(PairwiseListMatrix([10,20,30], false), ["a","b","c"])
3×3 Named PairwiseListMatrix{Int64,false,Array{Int64,1}}
A ╲ B │  a   b   c
──────┼───────────
a     │  0  10  20
b     │ 10   0  30
c     │ 20  30   0

julia> dict = to_dict(nplm, diagonal=false)
Dict{Symbol,Array{T,1} where T} with 3 entries:
  :values => [10, 20, 30]
  :j      => ["b", "c", "c"]
  :i      => ["a", "a", "b"]

julia> DataFrame(dict)
3×3 DataFrames.DataFrame
│ Row │ i      │ j      │ values │
│     │ String │ String │ Int64  │
├─────┼────────┼────────┼────────┤
│ 1   │ a      │ b      │ 10     │
│ 2   │ a      │ c      │ 20     │
│ 3   │ b      │ c      │ 30     │
```
"""
function to_dict(plm::PairwiseListMatrix{T,D,TV};
                 diagonal::Bool=true,
                 labels::Vector{String} = getlabels(plm)) where {T,D,TV}
    N = plm.nelements
    L = diagonal ? div(N*(N+1),2) : div(N*(N-1),2)
    I = Array{String}(undef, L)
    J = Array{String}(undef, L)
    K = Array{T}(undef, L)
    t = 0
    @iterateupper plm diagonal begin
        t += 1
        I[t] = labels[i]
        J[t] = labels[j]
        K[t] = list[k]
    end
    Dict(:i => I, :j => J, :values => K)
end

function to_dict(nplm::NamedArray{T,2,PairwiseListMatrix{T,D,TV},DN};
                 diagonal::Bool = true,
                 labels::Vector{String} = getlabels(nplm)) where {T,D,TV,DN}
    to_dict(nplm.array, diagonal=diagonal, labels=labels)
end

"""
Creation of a `PairwiseListMatrix` from a `Matrix`, `DataFrame` or similar structure.
By default the columns with the labels for i (slow) and j (fast) are 1 and 2.
Values are taken from the column 3 by default.

```jldoctest
julia> using PairwiseListMatrices, Pkg, DelimitedFiles

julia> import PairwiseListMatrices

julia> filename = joinpath(dirname(pathof(PairwiseListMatrices)), "..", "test", "example.csv");

julia> dat = readdlm(filename, ',')
3×3 Array{Any,2}:
 "A"  "B"  10
 "A"  "C"  20
 "B"  "C"  30

julia> from_table(dat, false)
3×3 Named PairwiseListMatrix{Any,false,Array{Any,1}}
A ╲ B │       A        B        C
──────┼──────────────────────────
A     │ nothing       10       20
B     │      10  nothing       30
C     │      20       30  nothing
```

This is also useful to create a `PairwiseListMatrix` from a `DataFrame`:

```
julia> using PairwiseListMatrices, DataFrames, CSV, Pkg

julia> import PairwiseListMatrices

julia> filename = joinpath(dirname(pathof(PairwiseListMatrices)), "..", "test", "example.csv");

julia> df = CSV.read(filename, header=false)
3×3 DataFrames.DataFrame
│ Row │ Column1 │ Column2 │ Column3 │
│     │ String⍰ │ String⍰ │ Int64⍰  │
├─────┼─────────┼─────────┼─────────┤
│ 1   │ A       │ B       │ 10      │
│ 2   │ A       │ C       │ 20      │
│ 3   │ B       │ C       │ 30      │

julia> from_table(df, false)
3×3 Named PairwiseListMatrix{Union{Missing, Int64},false,Array{Union{Missing, Int64},1}}
A ╲ B │       A        B        C
──────┼──────────────────────────
A     │ missing       10       20
B     │      10  missing       30
C     │      20       30  missing
```
"""
function from_table(table,
                    diagonal::Bool;
                    labelcols::Vector{Int} = [1,2],
                    valuecol::Int = 3,
                    diagonalvalue = :default)
    @assert size(table,2) >= 3
    values = table[:,valuecol]
    if diagonalvalue == :default
        diagonalvalue = _diagonal_value(eltype(values))
    end
    plm  = PairwiseListMatrix(values, diagonal, diagonalvalue)
    nplm = NamedArray(plm)
    if length(labelcols) == 2
        unique_labels = unique(vcat(table[:, labelcols[1]],
                                    table[:, labelcols[2]]))
        setlabels!(nplm, String[ string(lab) for lab in unique_labels ])
    end
    nplm
end

# write...
# ========

"""
This function takes the filename as first argument and a `PairwiseListMatrix` as second
argument. If the `diagonal` keyword argument is `true` (default), the diagonal is included
in the output. The keyword argument `delim` (by default is `'\t'`) allows to modified the
character used as delimiter.

```jldoctest
julia> using PairwiseListMatrices, DelimitedFiles

julia> plm  = PairwiseListMatrix(trues(3), false)
3×3 PairwiseListMatrix{Bool,false,BitArray{1}}:
 false   true   true
  true  false   true
  true   true  false

julia> writedlm("example.csv", plm, diagonal=false, delim=',')

julia> println(read("example.csv", String))
1,2,true
1,3,true
2,3,true

```
"""
function DelimitedFiles.writedlm(filename::String,
                       plm::PairwiseListMatrix{T,D,TV};
                       diagonal::Bool = true,
                       delim::Char = '\t',
                       labels::Vector{String} = getlabels(plm)) where {T,D,TV}
    open(filename, "w") do fh
        @iterateupper plm diagonal println(fh, labels[i], delim, labels[j], delim, list[k])
    end
end

function DelimitedFiles.writedlm(filename::String,
                       nplm::NamedArray{T,2,PairwiseListMatrix{T,D,TV},DN};
                       diagonal::Bool = true,
                       delim::Char = '\t',
                       labels::Vector{String} = getlabels(nplm)) where {T,D,TV,DN}
    writedlm(filename, nplm.array, diagonal=diagonal, delim=delim, labels=labels)
end

# triu! & triu
# ============

function LinearAlgebra.triu!(mat::PairwiseListMatrix)
    throw(ErrorException("PairwiseListMatrix must be Symmetric, use triu instead of triu!"))
end
LinearAlgebra.triu!(mat::PairwiseListMatrix, k::Int) = triu!(mat)
LinearAlgebra.triu(mat::PairwiseListMatrix) = triu(Matrix(mat))
LinearAlgebra.triu(mat::PairwiseListMatrix, k::Int) = triu(Matrix(mat), k::Int)

# join
# ====

"""
This function join two PairwiseListMatrices by their labels, returning two
PairwiseListMatrices with same size and labels. There are 4 `kind`s of joins:

- `:inner` : Intersect. The output matrices only include the labels that are in both PairwiseListMatrices
- `:outer` : Union. Include the labels of the two PairwiseListMatrices.
- `:left` : Only use labels from the first argument.
- `:right` : Only use labels from the second argument.

`NaN`s are filled in where needed to complete joins. The default value for missing values
can be changed passing a tuple to `missing`.

```jldoctest
julia> using PairwiseListMatrices

julia> l = setlabels(PairwiseListMatrix([1.,2.,3.], false), ["a","b","c"]) # a b c
3×3 Named PairwiseListMatrix{Float64,false,Array{Float64,1}}
A ╲ B │   a    b    c
──────┼──────────────
a     │ 0.0  1.0  2.0
b     │ 1.0  0.0  3.0
c     │ 2.0  3.0  0.0

julia> r = setlabels(PairwiseListMatrix([1.,2.,3.], false), ["b","c","d"]) # b c d
3×3 Named PairwiseListMatrix{Float64,false,Array{Float64,1}}
A ╲ B │   b    c    d
──────┼──────────────
b     │ 0.0  1.0  2.0
c     │ 1.0  0.0  3.0
d     │ 2.0  3.0  0.0

julia> join(l, r, kind=:inner) # b c
(2×2 Named PairwiseListMatrix{Float64,false,Array{Float64,1}}
A ╲ B │   b    c
──────┼─────────
b     │ 0.0  3.0
c     │ 3.0  0.0, 2×2 Named PairwiseListMatrix{Float64,false,Array{Float64,1}}
A ╲ B │   b    c
──────┼─────────
b     │ 0.0  1.0
c     │ 1.0  0.0)

julia> join(l, r, kind=:outer) # a b c d
(4×4 Named PairwiseListMatrix{Float64,false,Array{Float64,1}}
A ╲ B │   a    b    c    d
──────┼───────────────────
a     │ 0.0  1.0  2.0  NaN
b     │ 1.0  0.0  3.0  NaN
c     │ 2.0  3.0  0.0  NaN
d     │ NaN  NaN  NaN  NaN, 4×4 Named PairwiseListMatrix{Float64,false,Array{Float64,1}}
A ╲ B │   a    b    c    d
──────┼───────────────────
a     │ NaN  NaN  NaN  NaN
b     │ NaN  0.0  1.0  2.0
c     │ NaN  1.0  0.0  3.0
d     │ NaN  2.0  3.0  0.0)

```
"""
function Base.join(
        left::NamedArray{L,2,PairwiseListMatrix{L,DL,VL},NL},
        right::NamedArray{R,2,PairwiseListMatrix{L,DR,VR},NR};
        kind::Symbol = :inner,
        missing::Tuple{L,R} = (L(NaN), R(NaN))
        ) where {L <: AbstractFloat, R <: AbstractFloat,DL,DR,VL,VR,NL,NR}

    labels_left = getlabels(left)
    labels_right = getlabels(right)

    if labels_left == labels_right
        return(left, right)
    end

    if kind == :inner
        out_labels = intersect(labels_left, labels_right)
        N = length(out_labels)
        out_L = PairwiseListMatrix(L, N, DL)
        out_R = PairwiseListMatrix(R, N, DR)
        for i in 1:N
            li = out_labels[i]
            for j in i:N
                lj = out_labels[j]
                out_L[i,j] = left[li, lj]
                out_R[i,j] = right[li, lj]
            end
        end
        return(setlabels(out_L, out_labels), setlabels(out_R, out_labels))
    elseif kind == :left
        out_labels = labels_left
        N = length(out_labels)
        out_R = PairwiseListMatrix(R, N, DR)
        for i in 1:N
            li = out_labels[i]
            flag_i_r = li in labels_right
            for j in i:N
                lj = out_labels[j]
                out_R[i,j] = (flag_i_r && (lj in labels_right)) ? right[li, lj] : missing[2]
            end
        end
        return(left, setlabels(out_R, out_labels))
    elseif kind == :right
        out_labels = labels_right
        N = length(out_labels)
        out_L = PairwiseListMatrix(L, N, DL)
        for i in 1:N
            li = out_labels[i]
            flag_i_l = li in labels_left
            for j in i:N
                lj = out_labels[j]
                out_L[i,j] = (flag_i_l  && (lj in labels_left)) ? left[li, lj] : missing[1]
            end
        end
        return(setlabels(out_L, out_labels), right)
    elseif kind == :outer
        out_labels = union(labels_left, labels_right)
        N = length(out_labels)
        out_L = PairwiseListMatrix(L, N, DL)
        out_R = PairwiseListMatrix(R, N, DR)
        for i in 1:N
            li = out_labels[i]
            flag_i_l = li in labels_left
            flag_i_r = li in labels_right
            for j in i:N
                lj = out_labels[j]
                out_L[i,j] = (flag_i_l && (lj in labels_left))  ? left[li, lj] : missing[1]
                out_R[i,j] = (flag_i_r && (lj in labels_right)) ? right[li, lj] : missing[2]
            end
        end
        return(setlabels(out_L , out_labels), setlabels(out_R, out_labels))
    else
        throw(ArgumentError("Unknown kind of join requested: use :inner, :left, :right or :outer"))
    end

end


text:
using HttpServer
using WebSockets
using JSON

type Token
    i::Int
    j::Int
    label::String
end

function Token(i::Int, j::Int)
    h = i + j
    id = Int(rem(h,length(posconf))) + 1
    label = posconf[id][1]
    Token(i, j, label)
end

function readconf(path::String)
    lines = open(readlines, path)
    map(lines) do line
        items = split(chomp(line), '\t')
        (items[1], items[2])
    end
end

const filepath = dirname(@__FILE__)
const clients = Dict()
const posconf = readconf("$(filepath)/pos.conf")
const entityconf = readconf("$(filepath)/entity.conf")

function create_error(msg::String)
    dict = Dict(
        "error": msg,
    )
    JSON.json(dict)
end

wsh = WebSocketHandler() do req, client
    println("Client: $(client.id) is connected.")
    while true
        #println("Request from $(client.id) recieved.")
        msg = String(read(client))
        if length(msg) > 2000
            write(client, create_error("Error: the mssage length exceeds 2000."))
            continue
        end

        json = try
            JSON.parse(msg)
        catch
            write(client, create_error("Error: invalid json format: $(msg)"))
            continue
        end
        if !haskey(json, "text")
            write(client, create_error("Error: key: text does not exist. Check the JSON format."))
            continue
        end
        text = json["text"]
        doc = Vector{Token}[]
        sents = split(text, '\n', keep=false)
        for sent in sents
            tokens = Token[]
            index = 1
            for i = 1:length(sent)
                c = sent[i]
                if c == ' '
                    index < i && push!(tokens, Token(index,i-1))
                    index = i + 1
                end
            end
            index <= length(sent) && push!(tokens, Token(index,length(sent)))
            length(tokens) > 0 && push!(doc, tokens)
        end

        sentences = []
        #trans_ja, trans_en, trans_cn = [], [], []
        for i = 1:length(doc)
            tokens = doc[i]
            #str = "I have a pen."
            annos = []
            for t in tokens
                push!(annos, ["pos", t.i-1,t.j-1,t.label])
                if rem(t.i,4) == 0
                    id = Int(rem(t.i,length(entityconf))) + 1
                    label = entityconf[id][1]
                    push!(annos, ["entity",t.i-1,t.j-1,label])
                end
            end
            push!(sentences, Dict("text"=>sents[i], "anno"=>annos))
            #push!(trans_ja, "日本語の翻訳結果です。")
            #push!(trans_en, "This is an example of English translation.")
            #push!(trans_cn, "这是一个中国的翻译结果。")
        end
        res = JSON.json(sentences)
        #println(res)
        write(client, res)
    end
end

onepage = readstring(joinpath(dirname(@__FILE__), "index.html"))
httph = HttpHandler() do req::Request, res::Response
    Response(onepage)
end
server = Server(httph, wsh)
run(server, 3000)

ext: jl

max_stars_repo_name: hshindo/JukaiNLP-web

max_issues_repo_name: hshindo/JukaiNLP-web

max_forks_repo_name: hshindo/JukaiNLP-web


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(domain::Domain)
  collect(domain)
end

@recipe function f(domain::Domain, data::AbstractVector)
  points = [centroid(domain, i) for i in 1:nelements(domain)]
  PointSet(points), data
end


text:
using Distributed
addprocs(16)
@everywhere include("test_rijke.jl")
collect_sensitivities()

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
# From https://github.com/MichaelHatherly/Lexicon.jl/issues/134#issue-121962949
module Lexicon134

export f, g

"
Prints stuff
"
f(a) = println(a)

"
a
b
c
d
e
"
g(x) = x

end


text:
using Distributions: Gaussian
using Test

@testset "Events" begin include("events.jl") end
@testset "TrajectoryTests" begin include("trajectories.jl") end
@testset "TrajectoryProbability" begin include("probability.jl") end
@testset "SimpleSystem" begin include("SimpleSystem.jl") end
@testset "ComplexSystem" begin include("ComplexSystem.jl") end
@testset "SMCEstimate" begin include("SMCEstimate.jl") end
@testset "MCMCMoves" begin include("mcmc_moves.jl") end
@testset "Gaussian" begin include("gaussian.jl") end


text:
#!/bin/julia

# Path to potential REQUIRE file is passed on the command line.
require_file = ARGS[1]

if VERSION < v"0.7-"
  pkg_dir = "$(ENV["JULIA_PKGDIR"])/v$(VERSION.major).$(VERSION.minor)"
  open("$pkg_dir/REQUIRE", "a") do io
    write(io, read(require_file))
  end
  Pkg.resolve();
  Reqs = Pkg.Reqs
else
  using Pkg
  Reqs = Pkg.Pkg2.Reqs
  emptyversionlower = v"0.0.0-"
  for reqline in Reqs.read(require_file)
    if reqline isa Reqs.Requirement
        pkg = String(reqline.package)
        if pkg == "julia" continue end
        version = try; reqline.versions.intervals[1].lower; catch; emptyversionlower; end
        if version != emptyversionlower
          Pkg.add(PackageSpec(name=pkg, version=version))
        else
          Pkg.add(pkg)
        end
    end
  end
end
# Precompile the packages
for reqline in Reqs.read(require_file)
  if reqline isa Reqs.Requirement
      pkg = reqline.package
      pkg != "julia" && eval(:(using $(Symbol(pkg))))
  end
end


text:
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: May 2020
#

ftol = 1f-6

function test_model(ndim; tti=false)
    n = Tuple(121 for i=1:ndim)
    o = Tuple(0. for i=1:ndim)
    d = Tuple(10. for i=1:ndim)
    m = .5f0 .+ rand(Float32, n...)
    if !tti
        model = Model(n, d, o, m; nb=10)
        model2 = Model(n, d, o, m, nothing; nb=10)
        @test model.n == model2.n
        @test model.d == model2.d
        @test model.o == model2.o
        @test model.m == model2.m
        @test [keys(model.params)...] == [:m]
        @test [keys(model2.params)...] == [:m]
    else
        epsilon = .1f0 * m
        delta = .1f0 * m
        theta = .1f0 * m
        phi = .1f0 * m
        model = Model(n, d, o, m; nb=10, epsilon=epsilon, delta=delta, theta=theta, phi=phi)
        @test all(k in keys(model.params) for k in [:m, :epsilon, :delta, :theta, :phi])
    end
    return model
end

# Test padding and its adjoint
function test_padding(ndim)

    model = test_model(ndim; tti=false)
    modelPy = devito_model(model, Options())

    m0 = model.m
    mcut = remove_padding(deepcopy(modelPy.m.data), modelPy.padsizes; true_adjoint=true)
    mdata = deepcopy(modelPy.m.data)

    @show dot(m0, mcut)/dot(mdata, mdata)
    @test isapprox(dot(m0, mcut), dot(mdata, mdata))
end

function test_limit_m(ndim, tti)
    model = test_model(ndim; tti=tti)

    @test get_dt(model) == calculate_dt(model)
    @test model[:m] == model.m
    mloc = model.m
    model.m .*= 0f0
    @test norm(model.m) == 0
    model.m = mloc
    @test model.m == mloc

    srcGeometry = example_src_geometry()
    recGeometry = example_rec_geometry(cut=true)
    dm = rand(Float32, model.n...)
    new_mod, dm_n = limit_model_to_receiver_area(srcGeometry, recGeometry, deepcopy(model), 100f0; pert=dm)

    # check inds
    inds = ndim == 3 ? [6:116, 1:11, 1:121] : [6:116, 1:121]
    @test new_mod.n[1] == 111
    @test new_mod.n[end] == 121
    if ndim == 3
        @test new_mod.n[2] == 11
    end

    @test isapprox(new_mod.m, model.m[inds...]; rtol=ftol)
    @test isapprox(dm_n, vec(dm[inds...]); rtol=ftol)
    if tti
        @test isapprox(new_mod.epsilon, model.epsilon[inds...]; rtol=ftol)
        @test isapprox(new_mod.delta, model.delta[inds...]; rtol=ftol)
        @test isapprox(new_mod.theta, model.theta[inds...]; rtol=ftol)
        @test isapprox(new_mod.phi, model.phi[inds...]; rtol=ftol)
    end

    ex_dm = extend_gradient(model, new_mod, dm_n)
    @test size(ex_dm) == size(dm)
    dmt = PhysicalParameter(dm_n, new_mod.n, model.d, new_mod.o)
    ex_dm = extend_gradient(model, new_mod, dmt)
    @test size(ex_dm.data) == size(dm)
    @test ex_dm.o == model.o
    @test ex_dm.n == model.n

end

function setup_3d()
    xrec = Array{Any}(undef, 2)
    yrec = Array{Any}(undef, 2)
    zrec = Array{Any}(undef, 2)
    for i=1:2
        xrec[i] = i .+ collect(0:10)
        yrec[i] = i .+ collect(0:10)
        zrec[i] = i
    end
    x3d, y3d, z3d = setup_3D_grid(xrec[1],yrec[1],zrec[1])
    for k=1:11
        @test all(y3d[(11*(k-1))+1:(11*(k-1))+11] .== k)
        @test all(x3d[k:11:end] .== k)
    end
    @test all(z3d[:] .== 1)
    x3d, y3d, z3d = setup_3D_grid(xrec,yrec,zrec)
    for i=1:2
        for k=1:11
            @test all(y3d[i][(11*(k-1))+1:(11*(k-1))+11] .== k + i -1)
            @test all(x3d[i][k:11:end] .== k + i - 1)
        end
        @test all(z3d[i][:] .== i)
    end
end

@testset "Test basic utilities" begin
    @timeit TIMEROUTPUT "Basic setup utilities" begin
        setup_3d()
        for ndim=[2, 3]
            test_padding(ndim)
        end
        opt = Options(frequencies=[[2.5, 4.5], [3.5, 5.5]])
        @test subsample(opt, 1).frequencies == [2.5, 4.5]
        @test subsample(opt, 2).frequencies == [3.5, 5.5]

        for ndim=[2, 3]
            for tti=[true, false]
                test_limit_m(ndim, tti)
            end
        end

        # Test model
        for ndim=[2, 3]
            for tti=[true, false]
                model =  test_model(ndim; tti=tti)

                # Default dt
                modelPy = devito_model(model, Options())
                @test get_dt(model) == calculate_dt(model)
                @test isapprox(modelPy.critical_dt, calculate_dt(model))
                @test isapprox(calculate_dt(model; dt=.5f0), .5f0)

                # Input dt
                modelPy = devito_model(model, Options(dt_comp=.5f0))
                @test modelPy.critical_dt == .5f0

                # Verify nt
                srcGeometry = example_src_geometry()
                recGeometry = example_rec_geometry(cut=true)
                nt1 = get_computational_nt(srcGeometry, recGeometry, model)
                nt2 = get_computational_nt(srcGeometry, model)
                nt3 = get_computational_nt(srcGeometry, recGeometry, model; dt=1f0)
                nt4 = get_computational_nt(srcGeometry, model; dt=1f0)
                @test all(nt1 .== (trunc(Int64, 1000f0 / calculate_dt(model)) + 1))
                @test all(nt1 .== nt2)
                @test all(nt3 .== 1001)
                @test all(nt4 .== 1001)
            end
        end
    end
end


text:
struct PoissonData
    X::Matrix{Int}
    x::Vector{Int}
    y::Vector{Int}
    Xunique::Vector{Vector{Int}}
    function PoissonData(X::Matrix{Int}, y::Vector{Int})
        x = denserank([X[i, :] for i in 1:size(X, 1)])
        Xunique = sort(unique([X[i, :] for i in 1:size(X, 1)]))
        new(X, x, y, Xunique)
    end
end

PoissonData(x::Vector{Int}, y::Vector{Int}) = PoissonData(x[:, :], y)

struct PoissonDDP <: AbstractDPM
    parent::DPM
    a0::Float64
    b0::Float64
    a1::Vector{Vector{Int}}
    b1::Vector{Vector{Int}}
    sumlogu::Vector{Vector{Float64}}
    gammaprior::Womack
    gamma::Vector{Bool}
    G::Int
    function PoissonDDP(
        rng::AbstractRNG,
        N::Int,
        G::Int;
        K0::Int = 1,
        αa0::Float64 = 2.0,
        αb0::Float64 = 4.0,
        a0::Float64 = 2.0,
        b0::Float64 = 4.0,
        ζ0::Float64 = 1.0,
    )
        parent = DPM(rng, N; K0, a0 = αa0, b0 = αb0)
        a1 = [a0 * ones(Int, G)]
        b1 = [b0 * ones(Int, G)]
        sumlogu = [zeros(G)]
        gammaprior = Womack(G - 1, ζ0)
        gamma = ones(Bool, G)
        new(parent, a0, b0, a1, b1, sumlogu, gammaprior, gamma, G)
    end
end

function parent_dpm(m::PoissonDDP)
    m.parent
end

function add_cluster!(m::PoissonDDP)
    @extract m : G a0 b0 a1 b1 sumlogu
    push!(a1, a0 * ones(G))
    push!(b1, b0 * ones(G))
    push!(sumlogu, zeros(G))
end

function update_suffstats!(m::PoissonDDP, data)
    @extract data : y x
    @extract m : a0 b0 a1 b1 sumlogu gamma
    d = cluster_labels(m)
    while length(a1) < cluster_capacity(m)
        add_cluster!(m)
    end
    for k in active_clusters(m)
        a1[k] .= a0
        b1[k] .= b0
        sumlogu[k] .= 0.0
    end
    for i = 1:length(y)
        di = d[i]
        zi = iszero(gamma[x[i]]) ? 1 : x[i]
        sumlogu[di][zi] += logfactorial(y[i])
        a1[di][zi] += y[i]
        b1[di][zi] += 1
    end
end

function update_suffstats!(m::PoissonDDP, data, i::Int, k1::Int, k2::Int)
    @extract data : y x
    @extract m : a0 b0 a1 b1 sumlogu gamma
    while length(a1) < cluster_capacity(m)
        add_cluster!(m)
    end
    zi = iszero(gamma[x[i]]) ? 1 : x[i]

    # Modify cluster/group k2/zi
    a1[k2][zi] += y[i]
    b1[k2][zi] += 1
    sumlogu[k2][zi] += logfactorial(y[i])

    # Modify cluster/group k1/zi
    a1[k1][zi] -= y[i]
    b1[k1][zi] -= 1
    sumlogu[k1][zi] -= logfactorial(y[i])
end

function logpredlik(m::PoissonDDP, data, i::Int, k::Int)
    d = cluster_labels(m)
    @extract m : a1 b1 gamma
    @extract data : y x
    j = iszero(gamma[x[i]]) ? 1 : x[i]
    a1kj = a1[k][j] - (d[i] == k) * y[i]
    b1kj = b1[k][j] - (d[i] == k)
    return logpdf(NegativeBinomial(a1kj, b1kj / (b1kj + 1)), y[i])
end

function logpredlik(m::PoissonDDP, train, predict, i::Int, k::Int)
    @extract m : a1 b1 gamma
    @extract predict : y x
    j = iszero(gamma[x[i]]) ? 1 : x[i]
    a1kj = a1[k][j]
    b1kj = b1[k][j]
    return logpdf(NegativeBinomial(a1kj, b1kj / (b1kj + 1)), y[i])
end

function logmglik(m::PoissonDDP, j::Int, k::Int)
    @extract m : a0 b0 a1 b1 sumlogu
    return (
        a0 * log(b0) - a1[k][j] * log(b1[k][j]) +
        loggamma(a1[k][j]) - loggamma(a0) -
        sumlogu[k][j]
    )
end

function update_gamma!(rng::AbstractRNG, m::PoissonDDP, data)
    @extract m : gammaprior gamma
    A = active_clusters(m)

    # Resample gamma[g], given the other gamma's
    for g = 2:length(gamma)
        # log-odds (numerator)
        gamma[g] = 1
        update_suffstats!(m, data)
        log_num = logpdf(gammaprior, gamma[2:end])
        for k ∈ A, j ∈ (1, g)
            log_num += logmglik(m, j, k)
        end

        # log-odds (denominator)
        gamma[g] = 0
        update_suffstats!(m, data)
        log_den = logpdf(gammaprior, gamma[2:end])
        for k ∈ A, j ∈ (1)
            log_den += logmglik(m, j, k)
        end

        # log-odds and new gamma[g]
        log_odds = log_num - log_den
        gamma[g] = rand(rng) <= 1 / (1 + exp(-log_odds))
    end
end

function update_hyperpars!(rng::AbstractRNG, m::PoissonDDP, data)
    update_gamma!(rng, m, data)
end

# Tõnu Kollo (tonu.kollo@ut.ee) University of Tartu, Tartu, Estonia


text:
function test_lrn_layer(backend::Backend, mode::LRNModeType, tensor_dim, T, eps)
  println("-- Testing LRN($(typeof(mode))) on $(typeof(backend)){$T}...")

  dims = tuple(rand(6:11, tensor_dim)...)
  op_dim = max(abs(rand(Int)) % tensor_dim, 1)

  println("    > Setup with dims $dims")

  input = rand(T, dims)
  input_blobs = Blob[make_blob(backend, input)]
  diff_blobs = Blob[make_blob(backend, input)]

  layer = LRNLayer(tops=[:output],bottoms=[:input], mode=mode, channel_dim=op_dim)
  state = setup(backend, layer, input_blobs, diff_blobs)

  @test size(state.blobs[1]) == size(input)

  println("    > Forward")
  forward(backend, state, input_blobs)
  expected_output = lrn_forward(input, state, op_dim)
  got_output = similar(input)
  copy!(got_output, state.blobs[1])
  @test all(abs.(got_output - expected_output) .< eps)

  println("    > Backward")
  top_diff = rand(T, size(input))
  copy!(state.blobs_diff[1], top_diff)
  backward(backend, state, input_blobs, diff_blobs)
  got_grad = zeros(T, size(input))
  copy!(got_grad, diff_blobs[1])
  expected_grad = lrn_backward(input, top_diff, state, op_dim)
  @test all(abs.(got_grad - expected_grad) .< eps)

  shutdown(backend, state)
end

function lrn_forward_across_channel{T}(input::Array{T}, state, op_dim)
  output = similar(input)
  pre_dim, chann_dim, post_dim = split_dims(input, op_dim)
  pre_pad = div(state.layer.kernel-1,2)
  post_pad = state.layer.kernel - pre_pad - 1

  canonical_input = reshape(input, (pre_dim, chann_dim, post_dim))
  canonical_output = reshape(output, (pre_dim, chann_dim, post_dim))

  for n = 1:post_dim
    for c = 1:chann_dim
      cstart = c-pre_pad
      cend   = min(c + post_pad, chann_dim)
      cstart = max(1, cstart)

      tmp = canonical_input[:,cstart:cend,n].^2 * (state.layer.scale / state.layer.kernel)
      tmp = (sum(tmp, 2) + state.layer.shift) .^ state.layer.power
      canonical_output[:,c,n] = canonical_input[:,c,n] ./ tmp
    end
  end

  return output
end
function lrn_forward_within_channel{T}(input::Array{T}, state)
  output = similar(input)
  width, height, channels, num = size(input)
  pooled_width = width; pooled_height = height
  kernel_size = state.layer.kernel^2
  pre_pad = div(state.layer.kernel-1,2)

  for n = 1:num
    for c = 1:channels
      for ph = 1:pooled_height
        for pw = 1:pooled_width
          hstart = ph - pre_pad
          wstart = pw - pre_pad
          hend = min(hstart + state.layer.kernel - 1, height)
          wend = min(wstart + state.layer.kernel - 1, width)
          hstart = max(1, hstart)
          wstart = max(1, wstart)

          tmp = (input[wstart:wend,hstart:hend,c,n]).^2 * state.layer.scale / kernel_size
          tmp = (sum(tmp) + state.layer.shift) .^ state.layer.power
          output[pw,ph,c,n] = input[pw,ph,c,n] / tmp
        end
      end
    end
  end

  return output
end
function lrn_forward{T}(input::Array{T}, state, op_dim)
  if isa(state.layer.mode, LRNMode.AcrossChannel)
    lrn_forward_across_channel(input, state, op_dim)
  elseif isa(state.layer.mode, LRNMode.WithinChannel)
    lrn_forward_within_channel(input, state)
  else
    error("Unknown LRN-mode $(state.layer.mode)")
  end
end

function lrn_backward_across_channel{T}(input::Array{T}, top_diff::Array{T}, state, op_dim)
  output = zeros(T, size(input))
  pre_dim, chann_dim, post_dim = split_dims(input, op_dim)
  pre_pad = div(state.layer.kernel-1,2)
  post_pad = state.layer.kernel - pre_pad - 1

  canonical_input = reshape(input, (pre_dim, chann_dim, post_dim))
  canonical_output = reshape(output, (pre_dim, chann_dim, post_dim))
  canonical_diff = reshape(top_diff, (pre_dim, chann_dim, post_dim))

  for n = 1:post_dim
    for c = 1:chann_dim
      cstart = c-pre_pad
      cend   = min(c + post_pad, chann_dim)
      cstart = max(1, cstart)

      tmp = canonical_input[:,cstart:cend,n].^2 * (state.layer.scale / state.layer.kernel)
      tmp = (sum(tmp, 2) + state.layer.shift)

      canonical_output[:,c,n] += tmp .^ (-state.layer.power) .* canonical_diff[:,c,n]

      tmp = -state.layer.power * tmp .^ (-state.layer.power - 1)
      tmp = 2 * state.layer.scale / state.layer.kernel * tmp
      canonical_output[:,cstart:cend,n] += tmp .* canonical_input[:,cstart:cend,n] .*
          canonical_input[:,c,n] .* canonical_diff[:,c,n]
    end
  end

  return output
end
function lrn_backward_within_channel{T}(input::Array{T}, top_diff::Array{T}, state)
  output = zeros(T, size(input))
  width, height, channels, num = size(input)
  pooled_width = width; pooled_height = height
  kernel_size = state.layer.kernel^2
  pre_pad = div(state.layer.kernel-1,2)

  for n = 1:num
    for c = 1:channels
      for ph = 1:pooled_height
        for pw = 1:pooled_width
          hstart = ph - pre_pad
          wstart = pw - pre_pad
          hend = min(hstart + state.layer.kernel - 1, height)
          wend = min(wstart + state.layer.kernel - 1, width)
          hstart = max(1, hstart)
          wstart = max(1, wstart)

          tmp = (input[wstart:wend,hstart:hend,c,n]).^2 * state.layer.scale / kernel_size
          tmp = (sum(tmp) + state.layer.shift)
          output[pw,ph,c,n] += tmp .^ (-state.layer.power) * top_diff[pw,ph,c,n]

          tmp = -state.layer.power * tmp .^ (-state.layer.power-1)
          tmp = 2 * state.layer.scale / kernel_size * tmp
          output[wstart:wend, hstart:hend, c, n] +=
              tmp * input[wstart:wend, hstart:hend, c, n] * input[pw,ph,c,n] * top_diff[pw,ph,c,n]
        end
      end
    end
  end

  return output
end

function lrn_backward{T}(input::Array{T}, top_diff::Array{T}, state, op_dim)
  if isa(state.layer.mode, LRNMode.AcrossChannel)
    lrn_backward_across_channel(input, top_diff, state, op_dim)
  elseif isa(state.layer.mode, LRNMode.WithinChannel)
    lrn_backward_within_channel(input, top_diff, state)
  else
    error("Unknown LRN-mode $(state.layer.mode)")
  end
end

function test_lrn_layer(backend::Backend, mode::LRNModeType, T, eps)
  test_lrn_layer(backend, mode, 4, T, eps)
end

function test_lrn_layer(backend::Backend, T, eps)
  test_lrn_layer(backend, LRNMode.AcrossChannel(), T, eps)
  test_lrn_layer(backend, LRNMode.WithinChannel(), T, eps)
end

function test_lrn_layer(backend::Backend)
  test_lrn_layer(backend, Float32, 1e-3)
  test_lrn_layer(backend, Float64, 1e-9)
end

if test_cpu
  test_lrn_layer(backend_cpu)
end
if test_gpu
  test_lrn_layer(backend_gpu)
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#
"""
$(TYPEDEF)

Epigraph of Euclidean norm (AKA second-order) cone of dimension `dim`.

    $(FUNCTIONNAME){T}(dim::Int)
"""
mutable struct EpiNormEucl{T <: Real} <: Cone{T}
    dim::Int

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}

    dist::T

    function EpiNormEucl{T}(dim::Int) where {T <: Real}
        @assert dim >= 2
        cone = new{T}()
        cone.dim = dim
        return cone
    end
end

use_dual_barrier(::EpiNormEucl) = false

reset_data(cone::EpiNormEucl) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = false)

use_sqrt_hess_oracles(::Int, cone::EpiNormEucl) = true

get_nu(cone::EpiNormEucl) = 2

function set_initial_point!(
    arr::AbstractVector,
    cone::EpiNormEucl{T},
    ) where {T <: Real}
    arr .= 0
    arr[1] = sqrt(T(get_nu(cone)))
    return arr
end

# TODO refac with dual feas check
function update_feas(cone::EpiNormEucl{T}) where T
    @assert !cone.feas_updated
    u = cone.point[1]

    if u > eps(T)
        @views w = cone.point[2:end]
        cone.dist = (abs2(u) - sum(abs2, w)) / 2
        cone.is_feas = (cone.dist > eps(T))
    else
        cone.is_feas = false
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::EpiNormEucl{T}) where T
    u = cone.dual_point[1]

    if u > eps(T)
        w = view(cone.dual_point, 2:cone.dim)
        @views dual_dist = abs2(u) - sum(abs2, cone.dual_point[2:end])
        return (dual_dist > 2 * eps(T))
    end

    return false
end

function update_grad(cone::EpiNormEucl)
    @assert cone.is_feas

    @. cone.grad = cone.point / cone.dist
    cone.grad[1] *= -1

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::EpiNormEucl)
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)

    mul!(cone.hess.data, cone.grad, cone.grad')
    inv_dist = inv(cone.dist)
    @inbounds for j in eachindex(cone.grad)
        cone.hess[j, j] += inv_dist
    end
    cone.hess[1, 1] -= inv_dist + inv_dist

    cone.hess_updated = true
    return cone.hess
end

function update_inv_hess(cone::EpiNormEucl)
    @assert cone.is_feas
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)

    mul!(cone.inv_hess.data, cone.point, cone.point')
    @inbounds for j in eachindex(cone.grad)
        cone.inv_hess[j, j] += cone.dist
    end
    cone.inv_hess[1, 1] -= cone.dist + cone.dist

    cone.inv_hess_updated = true
    return cone.inv_hess
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiNormEucl,
    )
    @assert cone.is_feas
    u = cone.point[1]
    w = @view cone.point[2:end]

    @inbounds for j in 1:size(prod, 2)
        uj = arr[1, j]
        wj = @view arr[2:end, j]
        ga = (dot(w, wj) - u * uj) / cone.dist
        prod[1, j] = -ga * u - uj
        @. @views prod[2:end, j] = ga * w + wj
    end
    @. prod ./= cone.dist

    return prod
end

function inv_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiNormEucl,
    )
    @assert cone.is_feas

    @inbounds for j in 1:size(prod, 2)
        @views pa = dot(cone.point, arr[:, j])
        @. @views prod[:, j] = pa * cone.point
    end
    @. @views prod[1, :] -= cone.dist * arr[1, :]
    @. @views prod[2:end, :] += cone.dist * arr[2:end, :]

    return prod
end

function sqrt_hess_prod!(
    prod::AbstractVecOrMat{T},
    arr::AbstractVecOrMat{T},
    cone::EpiNormEucl{T},
    ) where {T <: Real}
    @assert cone.is_feas
    u = cone.point[1]
    w = @view cone.point[2:end]

    rt2 = sqrt(T(2))
    distrt2 = cone.dist * rt2
    rtdist = sqrt(cone.dist)
    urtdist = u + rtdist * rt2
    @inbounds for j in 1:size(arr, 2)
        uj = arr[1, j]
        @views wj = arr[2:end, j]
        dotwwj = dot(w, wj)
        prod[1, j] = (u * uj - dotwwj) / distrt2
        wmulj = (dotwwj / urtdist - uj) / distrt2
        @. @views prod[2:end, j] = w * wmulj + wj / rtdist
    end

    return prod
end

function inv_sqrt_hess_prod!(
    prod::AbstractVecOrMat{T},
    arr::AbstractVecOrMat{T},
    cone::EpiNormEucl{T},
    ) where {T <: Real}
    @assert cone.is_feas
    u = cone.point[1]
    w = @view cone.point[2:end]

    rt2 = sqrt(T(2))
    rtdist = sqrt(cone.dist)
    urtdist = u + rtdist * rt2
    @inbounds for j in 1:size(arr, 2)
        uj = arr[1, j]
        @views wj = arr[2:end, j]
        dotwwj = dot(w, wj)
        prod[1, j] = (u * uj + dotwwj) / rt2
        wmulj = (dotwwj / urtdist + uj) / rt2
        @. @views prod[2:end, j] = w * wmulj + wj * rtdist
    end

    return prod
end

function dder3(cone::EpiNormEucl, dir::AbstractVector)
    @assert cone.grad_updated
    dim = cone.dim
    dder3 = cone.dder3
    point = cone.point
    u = point[1]
    u_dir = dir[1]
    @views w = point[2:end]
    @views w_dir = dir[2:end]

    jdotpd = u * u_dir - dot(w, w_dir)
    hess_prod!(dder3, dir, cone)
    dotdHd = -dot(dir, dder3)
    dotpHd = dot(point, dder3)
    dder3 .*= jdotpd
    @. @views dder3[2:end] += dotdHd * w + dotpHd * w_dir
    dder3[1] += -dotdHd * u - dotpHd * u_dir
    dder3 ./= 2 * cone.dist

    return dder3
end


text:
using IntervalMatrices, Test, LinearAlgebra

@testset "Interval arithmetic" begin
    a = -1.5 ± 0.5
    b = -1..1
    @test a + b == -3.0..0.0
    @test a * b == -2.0..2.0
    @test a * b + a == -4.0..1.0
    @test a * (b + 1) == -4.0..0.0
end

@testset "Interval matrix arithmetic" begin
    a = 1.0..1.3; b = 2.0..3.4; c = -0.5 ± 0.5; d = 0.0 ± 0.1
    A = IntervalMatrix([a b; c d])

    # can add, TODO: add test
    B = A + A
    @test B isa IntervalMatrix

    # can multiply, TODO: add test
    B = A * A
    @test B isa IntervalMatrix

    # multiply interval and interval matrix
    x = 0.0..2.0
    @test x * A == A * x ==
          IntervalMatrix([0.0..2.6 0.0..6.8; -2.0..0.0 -0.2..0.2])
    # multiply scalar and interval matrix
    x = 1.0
    for A2 in [x * A, A * x]
        @test A2 == A && typeof(A2) == typeof(A)
    end
end

@testset "Interval matrix methods" begin
    m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.8..4.2 0.0..0.9])
    @test opnorm(m) == opnorm(m, Inf) ≈ 5.2
    @test opnorm(m, 1) ≈ 5.3
    l = inf(m)
    r = sup(m)
    m2 = copy(m)
    @test m2 isa IntervalMatrix && m.mat == m2.mat
end

@testset "Interval matrix exponential" begin
    m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9])
    a_over = expm_overapproximation(m, 1.0, 4)
    a_under = expm_underapproximation(m, 1.0, 4)
    @test a_over isa IntervalMatrix && a_under isa IntervalMatrix
end

@testset "Interval matrix split" begin
    m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9])
    c, s = split(m)
    @test c ≈ [-0.1 -4.; 4. -0.1]
    @test s ≈ [1. 0.1; 0.1 1.]
end

@testset "Interval matrix membership" begin
    m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9])
    a1 = [0. -4.;  4. 0.]
    a2 = [0. -3.; 4. 0.]
    @test a1 ∈ m
    @test a2 ∉ m
end

@testset "Interval matrix rand" begin
    m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9])
    for i in 1:3
        @test rand(m) ∈ m
    end
end

@testset "Interval matrix rand" begin
    m = IntervalMatrix([-1.1..0.9 -4.1.. -3.9; 3.9..4.1 -1.1..0.9])
    f = IntervalMatrices.correction_hull(m, 1e-3, 5)
end


text:
function read_from_documentation(raw_code)
    raw_code_without_spaces = replace(replace(raw_code, " " => ""), "\n" => "")
    code_hex_array = map(x -> parse(UInt16, x, base = 16), collect(raw_code_without_spaces))
    code_bit_string = string(string.(code_hex_array, base = 2, pad = 4)...)
    map(x -> parse(Int8, x, base = 2), collect(code_bit_string)) .* Int8(2) .- Int8(1)
end

const GALILEO_E1B_CODES = extend_front_and_back(read_in_codes(
    joinpath(dirname(pathof(GNSSSignals)), "..", "data", "codes_galileo_e1b.bin"),
    50,
    4092
))

"""
$(SIGNATURES)

Get codes of type GalileoE1B as a Matrix where each column
represents a PRN.
```julia-repl
julia> get_code(GalileoE1B)
```
"""
function get_codes(::Type{GalileoE1B})
    @view GALILEO_E1B_CODES[2:end - 2, :]
end

"""
$(SIGNATURES)

Get code length of GNSS system GalileoE1B.
```julia-repl
julia> get_code_length(GalileoE1B)
```
"""
@inline function get_code_length(::Type{GalileoE1B})
    4092
end

"""
$(SIGNATURES)

Get shortest code length of GNSS system GalileoE1B.
```julia-repl
julia> get_shortest_code_length(GalileoE1B)
```
"""
@inline function get_secondary_code_length(::Type{GalileoE1B})
    1
end

"""
$(SIGNATURES)

Get center frequency of GNSS system GalileoE1B.
```julia-repl
julia> get_center_frequency(GalileoE1B)
```
"""
@inline function get_center_frequency(::Type{GalileoE1B})
    1_575_420_000Hz
end

"""
$(SIGNATURES)

Get code frequency of GNSS system GalileoE1B.
```julia-repl
julia> get_code_frequency(GalileoE1B)
```
"""
@inline function get_code_frequency(::Type{GalileoE1B})
    1023_000Hz
end

"""
$(SIGNATURES)

Get data frequency of GNSS system GalileoE1B.
```julia-repl
julia> get_data_frequency(GalileoE1B)
```
"""
@inline function get_data_frequency(::Type{GalileoE1B})
    250Hz
end

"""
$(SIGNATURES)

Get code of type GalileoE1B at phase `phase` of PRN `prn`. Includes only BOC(1,1) at the
moment.
```julia-repl
julia> get_code(GalileoE1B, 10.3, 1)
```
"""
Base.@propagate_inbounds function get_code(::Type{GalileoE1B}, phase, prn::Integer)
    floored_2phase = floor(Int, 2 * phase)
    get_code_unsafe(
        GalileoE1B,
        mod(
            floored_2phase >> 1,
            get_code_length(GalileoE1B) * get_secondary_code_length(GalileoE1B)
        ),
        prn
    ) * (iseven(floored_2phase) * 2 - 1)
end

"""
$(SIGNATURES)

Get code of type GalileoE1B at phase `phase` of PRN `prn`.
The phase will not be wrapped by the code length. The phase has to smaller
than the code length. Includes only BOC(1,1) at the moment.
```julia-repl
julia> get_code_unsafe(GalileoE1B, 10.3, 1)
```
"""
Base.@propagate_inbounds function get_code_unsafe(::Type{GalileoE1B}, phase, prn::Integer)
    floored_2phase = floor(Int, 2 * phase)
    get_code_unsafe(GalileoE1B, floored_2phase >> 1, prn) *
        (iseven(floored_2phase) * 2 - 1)
end

"""
$(SIGNATURES)

Get code of GNSS system GalileoE1B at phase `phase` of PRN `prn`.
The phase will not be wrapped by the code length. The phase has to smaller
than the code length and must be an integer.
```julia-repl
julia> get_code_unsafe(GalileoE1B, 10, 1)
```
"""
Base.@propagate_inbounds function get_code_unsafe(
    ::Type{GalileoE1B},
    phase::Integer,
    prn::Integer
)
    GALILEO_E1B_CODES[2 + phase, prn]
end


text:
module SemiparametricAUC

# package code goes here
using Distributions
using DataStructures
using Iterators
using DataArrays
using DataFrames
import Base: log, mean, var, zeros, length, isa, values, sort, inv, diagm, sqrt, diag,hcat

include("calculate-auc.jl")
include("coef-table.jl")
include("sauc.jl")
include("calculate-auc-simulation.jl")
include("simulation-one-predictor.jl")

export calculate_auc
export coefs_table
export semiparametricAUC
export calculate_auc_simulation
export simulate_one_predictor

end # module


text:
import CSV

# This file implements using Cullen's Evaluated Atomic Data Library 1997 for emission probabilities.
# The file relax.csv contains five columns, the atomic number, 3 sub-shell indices and a probability.
# The first two columns indicate the atom and the sub-shell ionized.  The third and fourth columns
# contain the inner and outer sub-shells involved in the X-ray transition and the fifth contains
# the probability of an ionization in column 2 producing an X-ray in columns 3 and 4.  This is what
# Perkins would call the Total yield - An ionization in column 2 will initiate a cascade of transitions.
# When the atom finally returns to the ground state, it could have emitted zero or more x-rays and
# zero or more Augers.  A K shell vacancy could relax via L3 which could relax via M5 and so on.

struct CullenEADL <: NeXLAlgorithm end

let transitions_data, xrayweights_data

    function loadAltWeights()
        @info "Loading EADL transition rate data."
        nn = (
            1, # Shell index
            (2 for _ in 1:3)...,
            (3 for _ in 1:5)...,
            (4 for _ in 1:7)...,
            (5 for _ in 1:9)...,
            (6 for _ in 1:11)...,
            (7 for _ in 1:13)...
        ) # = collect(Iterators.flatten(collect(i for _ in 1:(2i-1)) for i in 1:7 ))
        trans = Dict{Tuple{Int,Int},Int}()
        xrw = Dict{Tuple{Int,Int},Dict{Tuple{Int,Int},Float64}}()
        for row in CSV.File(joinpath(dirname(pathof(@__MODULE__)), "..", "data", "relax.csv"))
            z, ionized, inner, outer, weight = row.ZZ, row.II, row.NN, row.OO, row.PP
            if (z <= 99) && hasedge(z, inner, FFASTDB) && hasedge(z, outer, FFASTDB) &&
            (nn[inner] != nn[outer])
                # There seems to be a problem with the L2-M1 and L3-M1 weights which I resolve with this ad-hoc fix.
                if (outer == 5) && ((inner == 4) || (inner == 3))
                    if z >= 29
                        weight *= max(0.1, 0.1 + ((0.9 * (z - 29.0)) / (79.0 - 29.0)))
                    else
                        weight *= max(0.1, 0.2 - ((0.1 * (z - 22.0)) / (29.0 - 22.0)))
                    end
                end
                get!(xrw, (z, ionized), Dict{Tuple{Int,Int}, Float64}())[(inner, outer)] = weight
                trans[(inner, outer)] = get(trans, (inner, outer), 0) + 1
            end
        end
        # Add these which aren't in Cullen (The weights are WAGs)
        for x in ((3, 1, 1, 2, 0.00001), (4, 1, 1, 2, 0.00005), (5, 1, 1, 3, 0.0002))
            xrw[(x[1], x[2])] = Dict( (x[3], x[4]) =>x[5] )
        end
        return ( trans, xrw )
    end
    transitions_data, xrayweights_data = loadAltWeights()
    global transitions()::Dict{Tuple{Int,Int},Int} = transitions_data
    global xrayweights()::Dict{Tuple{Int,Int},Dict{Tuple{Int,Int},Float64}} = xrayweights_data
end

"""
    totalWeight(z::Int, ionized::Int, inner::Int, outer::Int, ::Type{CullenEADL})

The line weight for the transition `(inner,outer)` which results from an ionization of `ionized`.
"""
function totalWeight(z::Int, ionized::Int, inner::Int, outer::Int, ::Type{CullenEADL})
    trs = get(xrayweights(), (z, ionized), nothing)
    return isnothing(trs) ? 0.0 : get(trs, (inner,outer), 0.0)
end

"""
    allTotalWeights(z::Int, ionized::Int)

Returns a Vector containing tuples `(inner, outer, weight)` for each transition which could
result when the specified shell is ionized.
"""
allTotalWeights(z::Int, ionized::Int, ::Type{CullenEADL})::Dict{Tuple{Int,Int},Float64} =
    get(xrayweights(), (z, ionized), Dict{Tuple{Int,Int},Float64}())

isAvailable(z::Int, inner::Int, outer::Int, ::Type{CullenEADL}) =
    totalWeight(z, inner, inner, outer, CullenEADL) > 0.0


text:
"""
    scattering_field(args)

Returns a function which gives the average scattering coefficients for any vector `x` inside the material. This field is defined by Equation (3.13) in [AL Gower and G Kristensson, "Effective waves for random three-dimensional particulate materials", (2021)](https://arxiv.org/pdf/2010.00934.pdf)
"""
scattering_field

"Calculates the effective wavenumbers and return Vector{EffectivePlaneWaveMode}."
function WaveModes(ω::T, source::AbstractSource, material::Material{Dim,S,Sps}; kws...) where {T,Dim,S<:Shape{Dim},Sps<:Species{T,Dim}} # without the parametric types we get a "Unreachable reached" error

    # The wavenumbers are calculated without knowledge of the materail symmetry. This is because the plane-wave symmetry leads to all possible wavenumbers and is simple to calculate.
    k_effs = wavenumbers(ω, source.medium, material.species; numberofparticles = material.numberofparticles, kws... )

    # The wavemodes need to know the material symmetry as the eigenvectors do depend on material shape and symetry.
    wave_effs = [
        WaveMode(ω, k_eff, source, material; kws...)
    for k_eff in k_effs]

    return wave_effs
end

"""
    WaveMode(ω::T, wavenumber::Complex{T}, eigenvectors::Array{Complex{T}}, ::SetupSymmetry; kws...)

Returns a concrete subtype of AbstractWaveMode depending on the SetupSymmetry. The returned type should have all the necessary fields to calculate scattered waves (currently not true for EffectivePlanarWaves).
"""
function WaveMode(ω::T, wavenumber::Complex{T}, source::AbstractSource, material::Material{Dim}; kws...) where {T,Dim}

    eigvectors = eigenvectors(ω, wavenumber, source, material; kws...)

    α = solve_boundary_condition(ω, wavenumber, eigvectors, source, material; kws...)

    # After this normalisation, sum(eigvectors, dims = 3) will satisfy the boundary conditions
    eigvectors = [eigvectors[i] * α[i[3]] for i in CartesianIndices(eigvectors)]

    return EffectiveRegularWaveMode(ω, wavenumber, source, material, eigvectors; kws...)
end

function WaveMode(ω::T, wavenumber::Complex{T}, psource::PlaneSource{T,Dim,1}, material::Material{Dim,Halfspace{T,Dim}};
    tol::T = 1e-6, kws...) where {T,Dim}

    direction = transmission_direction(wavenumber, (ω / psource.medium.c) * psource.direction, material.shape.normal)
    eigvectors = eigenvectors(ω, wavenumber, psource, material; direction_eff = direction, kws...)

    α = solve_boundary_condition(ω, wavenumber, eigvectors, psource, material; kws...)

    # After this normalisation, sum(eigvectors, dims = 3) will satisfy the boundary conditions
    eigvectors = [eigvectors[i] * α[i[3]] for i in CartesianIndices(eigvectors)]

    return EffectivePlaneWaveMode(ω, wavenumber, direction, eigvectors)
end

function WaveMode(ω::T, wavenumber::Complex{T}, psource::PlaneSource{T,Dim,1}, material::Material{Dim,Plate{T,Dim}}; kws...) where {T,Dim}

    # First we calculate the outward pointing normal
    n = material.shape.normal;
    n = - n .* sign(real(dot(n,psource.direction)));

    k = ω / psource.medium.c

    direction1 = transmission_direction(wavenumber, k * psource.direction, n)
    eigvectors1 = eigenvectors(ω, wavenumber, psource, material; direction_eff = direction1, kws...)

    # we choose direction2 so that k2 .* direction2 = - k1 .* direction1, where k1 = wavenumber, and k2 = - wavenumber
    direction2 = direction1
    eigvectors2 = eigenvectors(ω, - wavenumber, psource, material; direction_eff = direction2, kws...)

    α = solve_boundary_condition(ω, wavenumber, eigvectors1, eigvectors2, psource, material; kws...)

    # apply normalisation
    eigvectors1 = eigvectors1 .* α[1]
    eigvectors2 = eigvectors2 .* α[2]

    mode1 = EffectivePlaneWaveMode(ω, wavenumber, direction1, eigvectors1)
    mode2 = EffectivePlaneWaveMode(ω, - wavenumber, direction2, eigvectors2)

    return [mode1,mode2]
end

# eigensystem(ω::T, source::AbstractSource, material::Material; kws...) where T<:AbstractFloat = eigensystem(ω, source.medium, material.species, Symmetry(source,material); numberofparticles = material.numberofparticles, kws...)

function solve_boundary_condition(ω::T, wavenumber::Complex{T}, eigvectors::Array, source::AbstractSource, material::Material; kws...) where T
    return solve_boundary_condition(ω, wavenumber, eigvectors, source, material, Symmetry(source,material); kws...)
end

function solve_boundary_condition(ω::T, wavenumber::Complex{T}, eigvectors1::Array, eigvectors2::Array, source::AbstractSource, material::Material; kws...) where T
    return solve_boundary_condition(ω, wavenumber, eigvectors1, eigvectors2, source, material, Symmetry(source,material); kws...)
end

"""
    convert_eigenvector_basis(medium::PhysicalMedium,sym::AbstractSymmetry,eigvecs::Array)

The eigenvectors from high symmetric scenarios are smaller then the more general scenarios. This function pads with zeros the more symmetric cases to match more general cases, so that we can use the functions for both.
"""
convert_eigenvector_basis(medium::PhysicalMedium,sym::AbstractSymmetry,eigvecs::Array) = eigvecs

eigenvectors(ω::T, k_eff::Complex{T}, source::AbstractSource, material::Material; kws...) where T<:AbstractFloat = eigenvectors(ω, k_eff::Complex{T}, source.medium, material.species, Symmetry(source,material); numberofparticles = material.numberofparticles, kws...)

# For plane waves, it is simpler to write all cases in the format for the most general case. For example, for PlanarAzimuthalSymmetry the eignvectors are much smaller. So we will turn these into the more general eigvector case by padding it with zeros.
# function eigenvectors(ω::T, k_eff::Complex{T}, source::PlaneSource{T}, material::Material{Dim,S}; kws...) where {T<:AbstractFloat,Dim,S<:Union{Plate,Halfspace}}
#
#     eigvecs = eigenvectors(ω, k_eff, source.medium, material.species, Symmetry(source,material); kws...)
#
#     if Symmetry(source,material) == PlanarAzimuthalSymmetry{Dim}()
#         eigvecs = azimuthal_to_planar_eigenvector(typeof(source.medium),eigvecs)
#     end
#
#     return eigvecs
#
# end

function eigenvectors(ω::T, k_eff::Complex{T}, medium::PhysicalMedium, species::Vector{Sp}, symmetry::AbstractSymmetry;
        tol::T = 1e-4, kws...
    ) where {T<:AbstractFloat, Sp<:Specie{T}}

    MM = eigensystem(ω, medium, species, symmetry; kws...)

    # calculate eigenvectors
    MM_svd = svd(MM(k_eff))
    inds = findall(MM_svd.S .< tol)

    if isempty(inds)
        @warn("No eigenvectors found with the tolerance tol = $tol. Will use only one eigenvector with the eigenvalue $(MM_svd.S[end]), which should be less than tol.")
        inds = [length(MM_svd.S)]
    end

    #NOTE: MM(k_eff) ≈ MM_svd.U * diagm(0 => MM_svd.S) * MM_svd.Vt
    eigvectors = MM_svd.V[:,inds]

    # Reshape to separate different species and eigenvectors
    S = length(species)

    # pads with zeros if necessary to match the more general case with less symmetry
    eigvectors = convert_eigenvector_basis(medium,symmetry,reshape(eigvectors,(:,S,size(eigvectors,2))))

    return eigvectors
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
using BuiltUtils
using Base.Test

@testset "rpath" begin
    # macOS
    @test rpath_origin(:Darwin) == "-Wl,-rpath,'@loader_path/'"
    @test rpath_origin(:Darwin, escape=true) == rpath_origin(:Darwin, escape=false)

    # Windows
    @test rpath_origin(:NT, escape=true) == rpath_origin(:NT, escape=false) == ""

    # Everything else
    @test rpath_origin(:Linux) == "-Wl,-rpath,'\$ORIGIN' -Wl,-z,origin"
    @test rpath_origin(:Linux) == rpath_origin(:FreeBSD)
    @test rpath_origin(:Solaris, escape=true) == "-Wl,-rpath,'\$\$ORIGIN' -Wl,-z,origin"
end

@testset "whole archive" begin
    lib = "/path/to/nowhere.so"
    @test whole_archive(:Darwin, lib) == "-Xlinker -all_load " * lib
    @test whole_archive(:NT, lib) == "-Wl,--whole-archive " * lib * " -Wl,--no-whole-archive"
end

@testset "is executable" begin
    @test isexecutable("sh")  # Everyone should have sh, so though should be safe
    @test !isexecutable("./wacky-whatever.txt")
end

@testset "ldd" begin
    d = object_dependencies(Libdl.dlpath("libjulia"))
    @test d isa Dict{String,String}
    @test any(k->contains(k, "libLLVM"), keys(d))
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

doc"""
    RMSProp(; kwargs...)

Scale learning rates by dividing with the moving average of the root mean
squared (RMS) gradients. See [1] for further description.

### Arguments

* `η`: default `0.1`, learning rate.
* `ρ`: default `0.9`, gradient moving average decay factor.
* `ϵ`: default `1e-8`, small value added for numerical stability.
* `clip`: default `0`, gradient clipping.
  If positive, will clip the gradient into the range `[-clip, clip]`.
* `scale`: default `0`, gradient rescaling.
  If != 0, multiply the gradient with `scale` before updating.
  Often choose to be `1.0 / batch_size`.
  If leave it default, high-level API like `fit!` will set it to
  `1.0 / batch_size`, since `fit!` knows the `batch_size`.
* `λ`: default `0.00001`, weight decay is equivalent
  to adding a global l2 regularizer for all the parameters.

### Notes
`ρ` should be between 0 and 1. A value of `ρ` close to 1 will decay the
moving average slowly and a value close to 0 will decay the moving average
fast.

Using the step size `η` and a decay factor `ρ the
learning rate `ηₜ` is calculated as:

```math
\begin{align*}
  r_t &= ρ r_{t-1} + (1 - ρ)g^2 \\
  η_t &= \frac{η}{\sqrt{r_t + ϵ}}
\end{align*}
```

### References
1. Tieleman, T. and Hinton, G. (2012):
   Neural Networks for Machine Learning, Lecture 6.5 - rmsprop.
   Coursera. [http://www.youtube.com/watch?v=O3sxAc4hxZU]
   (http://www.youtube.com/watch?v=O3sxAc4hxZU) (formula @5:20)
"""
RMSProp

@defstruct RMSProp <: AbstractOptimizer (
  (η      :: Real = 0.001, η > 0),
  (ρ      :: Real = 0.9,   0 < ρ < 1),
  (ϵ      :: Real = 1e-8,  ϵ > 0),
  (clip   :: Real = 0,     clip >= 0),
   scale  :: Real = 0,
  (λ      :: Real = 1e-5,  λ >= 0),
  η_sched :: Any  = initlrsched(η)
)

create_state(::RMSProp, ::Int, W::NDArray) = zeros(size(W), context(W))

function update!(rms::RMSProp, ::Int, W::NDArray, ∇::NDArray, s::NDArray)
  η = get(rms.η_sched)
  ρ = rms.ρ
  ϵ = rms.ϵ

  normgrad!(rms, W, ∇)

  @inplace s .*= ρ
  @inplace s .+= (1 - ρ) .* (∇.^2)
  @inplace W .+= -η .* ∇ ./ sqrt(s .+ ϵ)  # FIXME: sqrt should be dot-call
end

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet