text:
# This file is based on JuliaLang/LinearAlgebra/src/lu.jl, a part of Julia.
# License is MIT: https://julialang.org/license

####################
# UL Factorization #
####################
"""
    UL <: Factorization

Matrix factorization type of the `UL` factorization of a square matrix `A`. This
is the return type of [`ul`](@ref), the corresponding matrix factorization function.

The individual components of the factorization `F::UL` can be accessed via [`getproperty`](@ref):

| Component | Description                              |
|:----------|:-----------------------------------------|
| `F.U`     | `U` (upper triangular) part of `UL`      |
| `F.L`     | `L` (unit lower triangular) part of `UL` |
| `F.p`     | (right) permutation `Vector`             |
| `F.P`     | (right) permutation `Matrix`             |

Iterating the factorization produces the components `F.U`, `F.L`, and `F.p`.

# Examples
```jldoctest
julia> A = [4 3; 6 3]
2×2 Array{Int64,2}:
 4  3
 6  3

julia> F = ul(A)
UL{Float64,Array{Float64,2}}
L factor:
2×2 Array{Float64,2}:
 1.0       0.0
 0.666667  1.0
U factor:
2×2 Array{Float64,2}:
 6.0  3.0
 0.0  1.0

julia> F.L * F.U == A[F.p, :]
true

julia> l, u, p = ul(A); # destructuring via iteration

julia> l == F.L && u == F.U && p == F.p
true
```
"""
struct UL{T,S<:AbstractMatrix{T},IPIV<:AbstractVector{<:Integer}} <: Factorization{T}
    factors::S
    ipiv::IPIV
    info::BlasInt

    function UL{T,S,IPIV}(factors, ipiv, info) where {T,S<:AbstractMatrix{T},IPIV<:AbstractVector{<:Integer}}
        require_one_based_indexing(factors)
        new{T,S,IPIV}(factors, ipiv, info)
    end
end
function UL(factors::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer}, info::Integer) where {T}
    UL{T,typeof(factors),typeof(ipiv)}(factors, ipiv, BlasInt(info))
end
function UL{T}(factors::AbstractMatrix, ipiv::AbstractVector{<:Integer}, info::Integer) where {T}
    UL(convert(AbstractMatrix{T}, factors),
       ipiv,
       info)
end

function UL{T,S}(factors::AbstractMatrix, ipiv::AbstractVector{<:Integer}, info::Integer) where {T,S<:AbstractMatrix{T}}
    UL(convert(S, factors),
       ipiv,
       info)
end

# iteration for destructuring into components
Base.iterate(S::UL) = (S.U, Val(:L))
Base.iterate(S::UL, ::Val{:L}) = (S.L, Val(:p))
Base.iterate(S::UL, ::Val{:p}) = (S.p, Val(:done))
Base.iterate(S::UL, ::Val{:done}) = nothing

adjoint(F::UL) = Adjoint(F)
transpose(F::UL) = Transpose(F)

# AbstractMatrix
function ul!(A::AbstractMatrix{T}, pivot::Union{Val{false}, Val{true}} = Val(true);
             check::Bool = true) where T<:BlasFloat
    return generic_ulfact!(A, pivot; check = check)
end
# function ul!(A::HermOrSym, pivot::Union{Val{false}, Val{true}} = Val(true); check::Bool = true)
#     copytri!(A.data, A.uplo, isa(A, Hermitian))
#     ul!(A.data, pivot; check = check)
# end

"""
    ul!(A, pivot=Val(true); check = true) -> UL

`ul!` is the same as [`ul`](@ref), but saves space by overwriting the
input `A`, instead of creating a copy. An [`InexactError`](@ref)
exception is thrown if the factorization produces a number not representable by the
element type of `A`, e.g. for integer types.

# Examples
```jldoctest
julia> A = [4. 3.; 6. 3.]
2×2 Array{Float64,2}:
 4.0  3.0
 6.0  3.0

julia> F = ul!(A)
UL{Float64,Array{Float64,2},Vector{Int}}
L factor:
2×2 Array{Float64,2}:
 1.0       0.0
 0.666667  1.0
U factor:
2×2 Array{Float64,2}:
 6.0  3.0
 0.0  1.0

julia> iA = [4 3; 6 3]
2×2 Array{Int64,2}:
 4  3
 6  3

julia> ul!(iA)
ERROR: InexactError: Int64(0.6666666666666666)
Stacktrace:
[...]
```
"""
ul!(A::AbstractMatrix, pivot::Union{Val{false}, Val{true}} = Val(true); check::Bool = true) =
    generic_ulfact!(A, pivot; check = check)
function generic_ulfact!(A::AbstractMatrix{T}, ::Val{Pivot} = Val(true);
                         check::Bool = true) where {T,Pivot}
    m, n = size(A)
    minmn = min(m,n)
    info = 0
    ipiv = Vector{BlasInt}(undef, minmn)
    @inbounds begin
        for k = minmn:-1:1
            # find index max
            kp = k
            if Pivot
                amax = abs(zero(T))
                for i = 1:k
                    absi = abs(A[i,k])
                    if absi > amax
                        kp = i
                        amax = absi
                    end
                end
            end
            ipiv[k] = kp
            if !iszero(A[kp,k])
                if k != kp
                    # Interchange
                    for i = 1:n
                        tmp = A[k,i]
                        A[k,i] = A[kp,i]
                        A[kp,i] = tmp
                    end
                end
                # Scale first coulmn
                Akkinv = inv(A[k,k])
                for i = 1:k-1
                    A[i,k] *= Akkinv
                end
            elseif info == 0
                info = k
            end
            # Update the rest
            for j = 1:k-1
                for i = 1:k-1
                    A[i,j] -= A[i,k]*A[k,j]
                end
            end
        end
    end
    check && checknonsingular(info, Val{Pivot}())
    return UL{T}(A, ipiv, convert(BlasInt, info))
end

function ultype(T::Type)
    # In generic_ulfact!, the elements of the lower part of the matrix are
    # obtained using the division of two matrix elements. Hence their type can
    # be different (e.g. the division of two types with the same unit is a type
    # without unit).
    # The elements of the upper part are obtained by U - U * L
    # where U is an upper part element and L is a lower part element.
    # Therefore, the types LT, UT should be invariant under the map:
    # (LT, UT) -> begin
    #     L = oneunit(UT) / oneunit(UT)
    #     U = oneunit(UT) - oneunit(UT) * L
    #     typeof(L), typeof(U)
    # end
    # The following should handle most cases
    UT = typeof(oneunit(T) - oneunit(T) * (oneunit(T) / (oneunit(T) + zero(T))))
    LT = typeof(oneunit(UT) / oneunit(UT))
    S = promote_type(T, LT, UT)
end

# for all other types we must promote to a type which is stable under division
"""
    ul(A, pivot=Val(true); check = true) -> F::UL

Compute the UL factorization of `A`.

When `check = true`, an error is thrown if the decomposition fails.
When `check = false`, responsibility for checking the decomposition's
validity (via [`issuccess`](@ref)) lies with the user.

In most cases, if `A` is a subtype `S` of `AbstractMatrix{T}` with an element
type `T` supporting `+`, `-`, `*` and `/`, the return type is `UL{T,S{T}}`. If
pivoting is chosen (default) the element type should also support [`abs`](@ref) and
[`<`](@ref).

The individual components of the factorization `F` can be accessed via [`getproperty`](@ref):

| Component | Description                         |
|:----------|:------------------------------------|
| `F.U`     | `U` (upper triangular) part of `UL` |
| `F.L`     | `L` (lower triangular) part of `UL` |
| `F.p`     | (right) permutation `Vector`        |
| `F.P`     | (right) permutation `Matrix`        |

Iterating the factorization produces the components `F.L`, `F.U`, and `F.p`.

The relationship between `F` and `A` is

`F.L*F.U == A[F.p, :]`

`F` further supports the following functions:

| Supported function               | `UL` | `UL{T,Tridiagonal{T}}` |
|:---------------------------------|:-----|:-----------------------|
| [`/`](@ref)                      | ✓    |                        |
| [`\\`](@ref)                     | ✓    | ✓                      |
| [`inv`](@ref)                    | ✓    | ✓                      |
| [`det`](@ref)                    | ✓    | ✓                      |
| [`logdet`](@ref)                 | ✓    | ✓                      |
| [`logabsdet`](@ref)              | ✓    | ✓                      |
| [`size`](@ref)                   | ✓    | ✓                      |

# Examples
```jldoctest
julia> A = [4 3; 6 3]
2×2 Array{Int64,2}:
 4  3
 6  3

julia> F = ul(A)
UL{Float64,Array{Float64,2},Vector{Int}}
L factor:
2×2 Array{Float64,2}:
 1.0       0.0
 0.666667  1.0
U factor:
2×2 Array{Float64,2}:
 6.0  3.0
 0.0  1.0

julia> F.L * F.U == A[F.p, :]
true

julia> l, u, p = ul(A); # destructuring via iteration

julia> l == F.L && u == F.U && p == F.p
true
```
"""
function _ul(layout, A::AbstractMatrix{T}, pivot::Union{Val{false}, Val{true}}=Val(true);
            check::Bool = true) where T
    S = ultype(T)
    ul!(copy_oftype(A, S), pivot; check = check)
end

ul(A::AbstractMatrix{T}, pivot::Union{Val{false}, Val{true}}=Val(true); check::Bool = true) where T =
    _ul(MemoryLayout(A), A, pivot; check=check)

ul(S::UL) = S
function ul(x::Number; check::Bool=true)
    info = x == 0 ? one(BlasInt) : zero(BlasInt)
    check && checknonsingular(info)
    return UL(fill(x, 1, 1), BlasInt[1], info)
end

function UL{T}(F::UL) where T
    M = convert(AbstractMatrix{T}, F.factors)
    UL{T,typeof(M)}(M, F.ipiv, F.info)
end
UL{T,S}(F::UL) where {T,S} = UL{T,S}(convert(S, F.factors), F.ipiv, F.info)
Factorization{T}(F::UL{T}) where {T} = F
Factorization{T}(F::UL) where {T} = UL{T}(F)

copy(A::UL{T,S}) where {T,S} = UL{T,S}(copy(A.factors), copy(A.ipiv), A.info)

size(A::UL)    = size(getfield(A, :factors))
size(A::UL, i) = size(getfield(A, :factors), i)

getU(F::UL) = getU(F, size(F.factors))
function getU(F::UL{T}, _) where T
    m, n = size(F)
    U = triu!(layout_getindex(getfield(F, :factors),1:m,1:min(m,n)))
    for i = 1:min(m,n); U[i,i] = one(T); end
    return U
end

getL(F::UL) = getL(F, size(F.factors))
function getL(F::UL, _)
    m, n = size(F)
    tril!(layout_getindex(getfield(F, :factors),1:min(m,n),1:n))
end

function getproperty(F::UL{T,<:AbstractMatrix}, d::Symbol) where T
    m, n = size(F)
    if d === :L
        return getL(F)
    elseif d === :U
        return getU(F)
    elseif d === :p
        return invperm(ipiv2perm(getfield(F, :ipiv), m))
    elseif d === :P
        return Matrix{T}(I, m, m)[:,invperm(F.p)]
    else
        getfield(F, d)
    end
end

Base.propertynames(F::UL, private::Bool=false) =
    (:L, :U, :p, :P, (private ? fieldnames(typeof(F)) : ())...)

issuccess(F::UL) = F.info == 0

function show(io::IO, mime::MIME{Symbol("text/plain")}, F::UL)
    if issuccess(F)
        summary(io, F); println(io)
        println(io, "U factor:")
        show(io, mime, F.U)
        println(io, "\nL factor:")
        show(io, mime, F.L)
    else
        print(io, "Failed factorization of type $(typeof(F))")
    end
end

_apply_inverse_ipiv_rows!(A::UL, B::AbstractVecOrMat) = _ipiv_rows!(A, axes(A.ipiv,1), B)
_apply_ipiv_rows!(A::UL, ipiv, B::AbstractVecOrMat) = _ipiv_rows!(A, reverse(axes(ipiv,1)), B)
_apply_ipiv_rows!(A::UL, ::AbstractUnitRange, B::AbstractVecOrMat) = B
_apply_ipiv_rows!(A::UL, B::AbstractVecOrMat) = _apply_ipiv_rows!(A, A.ipiv, B)

function _ipiv_rows!(A::UL, order::OrdinalRange, B::AbstractVecOrMat)
    for i = order
        if i != A.ipiv[i]
            _swap_rows!(B, i, A.ipiv[i])
        end
    end
    B
end

function _swap_rows!(B::AbstractVector, i::Integer, j::Integer)
    B[i], B[j] = B[j], B[i]
    B
end

function _swap_rows!(B::AbstractMatrix, i::Integer, j::Integer)
    for col = 1 : size(B, 2)
        B[i,col], B[j,col] = B[j,col], B[i,col]
    end
    B
end

function ldiv!(A::UL{<:Any,<:AbstractMatrix}, B::AbstractVecOrMat)
    _apply_ipiv_rows!(A, B)
    ArrayLayouts.ldiv!(LowerTriangular(A.factors), ArrayLayouts.ldiv!(UnitUpperTriangular(A.factors), B))
end

ldiv!(A::UL{<:Any,<:AbstractMatrix}, B::LayoutVector) = Base.invoke(ldiv!, Tuple{UL{<:Any,<:AbstractMatrix},AbstractVecOrMat}, A, B)
ldiv!(A::UL{<:Any,<:AbstractMatrix}, B::LayoutMatrix) = Base.invoke(ldiv!, Tuple{UL{<:Any,<:AbstractMatrix},AbstractVecOrMat}, A, B)

function ldiv!(transA::Transpose{<:Any,<:UL{<:Any,<:AbstractMatrix}}, B::AbstractVecOrMat)
    A = transA.parent
    ArrayLayouts.ldiv!(transpose(UnitUpperTriangular(A.factors)), ArrayLayouts.ldiv!(transpose(LowerTriangular(A.factors)), B))
    _apply_inverse_ipiv_rows!(A, B)
end

ldiv!(adjF::Adjoint{T,<:UL{T,<:AbstractMatrix}}, B::AbstractVecOrMat{T}) where {T<:Real} =
    (F = adjF.parent; ldiv!(transpose(F), B))

function ldiv!(adjA::Adjoint{<:Any,<:UL{<:Any,<:AbstractMatrix}}, B::AbstractVecOrMat)
    A = adjA.parent
    ArrayLayouts.ldiv!(adjoint(UnitUpperTriangular(A.factors)), ArrayLayouts.ldiv!(adjoint(LowerTriangular(A.factors)), B))
    _apply_inverse_ipiv_rows!(A, B)
end

(\)(A::Adjoint{<:Any,<:UL}, B::Adjoint{<:Any,<:AbstractVecOrMat}) = A \ copy(B)
(\)(A::Transpose{<:Any,<:UL}, B::Transpose{<:Any,<:AbstractVecOrMat}) = A \ copy(B)

function (/)(A::AbstractMatrix, F::Adjoint{<:Any,<:UL})
    T = promote_type(eltype(A), eltype(F))
    return adjoint(ldiv!(F.parent, copy_oftype(adjoint(A), T)))
end
# To avoid ambiguities with definitions in adjtrans.jl and factorizations.jl
(/)(adjA::Adjoint{<:Any,<:AbstractVector}, F::Adjoint{<:Any,<:UL}) = adjoint(F.parent \ adjA.parent)
(/)(adjA::Adjoint{<:Any,<:AbstractMatrix}, F::Adjoint{<:Any,<:UL}) = adjoint(F.parent \ adjA.parent)
function (/)(trA::Transpose{<:Any,<:AbstractVector}, F::Adjoint{<:Any,<:UL})
    T = promote_type(eltype(trA), eltype(F))
    return adjoint(ldiv!(F.parent, conj!(AbstractVector{T}(trA.parent))))
end
function (/)(trA::Transpose{<:Any,<:AbstractMatrix}, F::Adjoint{<:Any,<:UL})
    T = promote_type(eltype(trA), eltype(F))
    return adjoint(ldiv!(F.parent, conj!(AbstractMatrix{T}(trA.parent))))
end

function det(F::UL{T}) where T
    n = checksquare(F)
    issuccess(F) || return zero(T)
    P = one(T)
    c = 0
    @inbounds for i = 1:n
        P *= F.factors[i,i]
        if F.ipiv[i] != i
            c += 1
        end
    end
    s = (isodd(c) ? -one(T) : one(T))
    return P * s
end

function logabsdet(F::UL{T}) where T  # return log(abs(det)) and sign(det)
    n = checksquare(F)
    issuccess(F) || return log(zero(real(T))), log(one(T))
    c = 0
    P = one(T)
    abs_det = zero(real(T))
    @inbounds for i = 1:n
        dg_ii = F.factors[i,i]
        P *= sign(dg_ii)
        if F.ipiv[i] != i
            c += 1
        end
        abs_det += log(abs(dg_ii))
    end
    s = ifelse(isodd(c), -one(real(T)), one(real(T))) * P
    abs_det, s
end

inv!(A::UL{T,<:AbstractMatrix}) where {T} =
    ldiv!(A.factors, copy(A), Matrix{T}(I, size(A, 1), size(A, 1)))
inv(A::UL{<:BlasFloat,<:AbstractMatrix}) = inv!(copy(A))

# Tridiagonal

# # See dgttrf.f
# function ul!(A::Tridiagonal{T,V}, pivot::Union{Val{false}, Val{true}} = Val(true);
#              check::Bool = true) where {T,V}
#     n = size(A, 1)
#     info = 0
#     ipiv = Vector{BlasInt}(undef, n)
#     dl = A.dl
#     d = A.d
#     du = A.du
#     if dl === du
#         throw(ArgumentError("off-diagonals of `A` must not alias"))
#     end
#     dl2 = fill!(similar(d, n-2), 0)::V

#     @inbounds begin
#         for i = 1:n
#             ipiv[i] = i
#         end
#         for i = n:-1:3
#             # pivot or not?
#             if pivot === Val(false) || abs(d[i]) >= abs(dl[i])
#                 # No interchange
#                 if d[i] != 0
#                     fact = du[i-1]/d[i]
#                     du[i-1] = fact
#                     d[i-1] -= fact*dl[i-1]
#                     dl2[i] = 0
#                 end
#             else
#                 # Interchange
#                 fact = d[i]/du[i-1]
#                 d[i] = du[i-1]
#                 du[i-1] = fact
#                 tmp = dl[i-1]
#                 dl[i-1] = d[i-1]
#                 d[i-1] = tmp - fact*d[i-1]
#                 dl2[i] = dl[i-1]
#                 dl[i-1] = -fact*dl[i-1]
#                 ipiv[i] = i-1
#             end
#         end
#         if n > 1
#             i = 2
#             if pivot === Val(false) || abs(d[i]) >= abs(dl[i])
#                 if d[i] != 0
#                     fact = du[i-1]/d[i]
#                     du[i-1] = fact
#                     d[i-1] -= fact*dl[i-1]
#                 end
#             else
#                 fact = d[i]/du[i-1]
#                 d[i] = du[i-1]
#                 du[i-1] = fact
#                 tmp = dl[i-1]
#                 dl[i-1] = d[i-1]
#                 d[i-1] = tmp - fact*d[i-1]
#                 ipiv[i] = i-1
#             end
#         end
#         # check for a zero on the diagonal of U
#         for i = 1:n
#             if d[i] == 0
#                 info = i
#                 break
#             end
#         end
#     end
#     B = Tridiagonal{T,V}(dl, d, du, dl2)
#     check && checknonsingular(info, pivot)
#     return UL{T,Tridiagonal{T,V}}(B, ipiv, convert(BlasInt, info))
# end

function getU(F::UL{T,<:Tridiagonal}, _) where T
    m, n = size(F)
    @assert m == n
    data = getfield(F, :factors)
    Bidiagonal(ones(T,m), copy(data.du), :U)
end
function getL(F::UL{T,<:Tridiagonal}, _) where T
    m, n = size(F)
    @assert m == n
    data = getfield(F, :factors)
    Bidiagonal(copy(data.d), copy(data.dl), :L)
end

# function getproperty(F::UL{T,Tridiagonal{T,V}}, d::Symbol) where {T,V}
#     m, n = size(F)
#     if d === :U
#         du = getfield(getfield(F, :factors), :du)
#         U = Array(Bidiagonal(fill!(similar(du, n), one(T)), du, d))
#         for i = 2:n
#             tmp = U[getfield(F, :ipiv)[i], i+1:end]
#             U[getfield(F, :ipiv)[i], i+1:end] = U[i, i+1:end]
#             U[i, i+1:end] = tmp
#         end
#         return U
#     elseif d === :L
#         L = Array(Bidiagonal(getfield(getfield(F, :factors), :d), getfield(getfield(F, :factors), :dl), d))
#         for i = 1:n - 2
#             L[i,i + 2] = getfield(getfield(F, :factors), :du2)[i]
#         end
#         return L
#     elseif d === :p
#         return ipiv2perm(getfield(F, :ipiv), m)
#     elseif d === :P
#         return Matrix{T}(I, m, m)[:,invperm(F.p)]
#     end
#     return getfield(F, d)
# end

# # See dgtts2.f
# function ldiv!(A::UL{T,Tridiagonal{T,V}}, B::AbstractVecOrMat) where {T,V}
#     require_one_based_indexing(B)
#     n = size(A,1)
#     if n != size(B,1)
#         throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
#     end
#     nrhs = size(B,2)
#     dl = A.factors.dl
#     d = A.factors.d
#     du = A.factors.du
#     du2 = A.factors.du2
#     ipiv = A.ipiv
#     @inbounds begin
#         for j = 1:nrhs
#             for i = 1:n-1
#                 ip = ipiv[i]
#                 tmp = B[i+1-ip+i,j] - dl[i]*B[ip,j]
#                 B[i,j] = B[ip,j]
#                 B[i+1,j] = tmp
#             end
#             B[n,j] /= d[n]
#             if n > 1
#                 B[n-1,j] = (B[n-1,j] - du[n-1]*B[n,j])/d[n-1]
#             end
#             for i = n-2:-1:1
#                 B[i,j] = (B[i,j] - du[i]*B[i+1,j] - du2[i]*B[i+2,j])/d[i]
#             end
#         end
#     end
#     return B
# end

# function ldiv!(transA::Transpose{<:Any,<:UL{T,Tridiagonal{T,V}}}, B::AbstractVecOrMat) where {T,V}
#     require_one_based_indexing(B)
#     A = transA.parent
#     n = size(A,1)
#     if n != size(B,1)
#         throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
#     end
#     nrhs = size(B,2)
#     dl = A.factors.dl
#     d = A.factors.d
#     du = A.factors.du
#     du2 = A.factors.du2
#     ipiv = A.ipiv
#     @inbounds begin
#         for j = 1:nrhs
#             B[1,j] /= d[1]
#             if n > 1
#                 B[2,j] = (B[2,j] - du[1]*B[1,j])/d[2]
#             end
#             for i = 3:n
#                 B[i,j] = (B[i,j] - du[i-1]*B[i-1,j] - du2[i-2]*B[i-2,j])/d[i]
#             end
#             for i = n-1:-1:1
#                 if ipiv[i] == i
#                     B[i,j] = B[i,j] - dl[i]*B[i+1,j]
#                 else
#                     tmp = B[i+1,j]
#                     B[i+1,j] = B[i,j] - dl[i]*tmp
#                     B[i,j] = tmp
#                 end
#             end
#         end
#     end
#     return B
# end

# # Ac_ldiv_B!(A::UL{T,Tridiagonal{T}}, B::AbstractVecOrMat) where {T<:Real} = At_ldiv_B!(A,B)
# function ldiv!(adjA::Adjoint{<:Any,UL{T,Tridiagonal{T,V}}}, B::AbstractVecOrMat) where {T,V}
#     require_one_based_indexing(B)
#     A = adjA.parent
#     n = size(A,1)
#     if n != size(B,1)
#         throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
#     end
#     nrhs = size(B,2)
#     dl = A.factors.dl
#     d = A.factors.d
#     du = A.factors.du
#     du2 = A.factors.du2
#     ipiv = A.ipiv
#     @inbounds begin
#         for j = 1:nrhs
#             B[1,j] /= conj(d[1])
#             if n > 1
#                 B[2,j] = (B[2,j] - conj(du[1])*B[1,j])/conj(d[2])
#             end
#             for i = 3:n
#                 B[i,j] = (B[i,j] - conj(du[i-1])*B[i-1,j] - conj(du2[i-2])*B[i-2,j])/conj(d[i])
#             end
#             for i = n-1:-1:1
#                 if ipiv[i] == i
#                     B[i,j] = B[i,j] - conj(dl[i])*B[i+1,j]
#                 else
#                     tmp = B[i+1,j]
#                     B[i+1,j] = B[i,j] - conj(dl[i])*tmp
#                     B[i,j] = tmp
#                 end
#             end
#         end
#     end
#     return B
# end

rdiv!(B::AbstractMatrix, A::UL) = transpose(ldiv!(transpose(A), transpose(B)))
rdiv!(B::AbstractMatrix, A::Transpose{<:Any,<:UL}) = transpose(ldiv!(A.parent, transpose(B)))
rdiv!(B::AbstractMatrix, A::Adjoint{<:Any,<:UL}) = adjoint(ldiv!(A.parent, adjoint(B)))

# Conversions
AbstractMatrix(F::UL) = (F.L * F.U)[invperm(F.p),:]
AbstractArray(F::UL) = AbstractMatrix(F)
Matrix(F::UL) = Array(AbstractArray(F))
Array(F::UL) = Matrix(F)

# function Tridiagonal(F::UL{T,Tridiagonal{T,V}}) where {T,V}
#     n = size(F, 1)

#     dl  = copy(F.factors.dl)
#     d   = copy(F.factors.d)
#     du  = copy(F.factors.du)
#     du2 = copy(F.factors.du2)

#     for i = n - 1:-1:1
#         li         = dl[i]
#         dl[i]      = li*d[i]
#         d[i + 1]  += li*du[i]
#         if i < n - 1
#             du[i + 1] += li*du2[i]
#         end

#         if F.ipiv[i] != i
#             tmp   = dl[i]
#             dl[i] = d[i]
#             d[i]  = tmp

#             tmp      = d[i + 1]
#             d[i + 1] = du[i]
#             du[i]    = tmp

#             if i < n - 1
#                 tmp       = du[i + 1]
#                 du[i + 1] = du2[i]
#                 du2[i]    = tmp
#             end
#         end
#     end
#     return Tridiagonal(dl, d, du)
# end
# AbstractMatrix(F::UL{T,Tridiagonal{T,V}}) where {T,V} = Tridiagonal(F)
# AbstractArray(F::UL{T,Tridiagonal{T,V}}) where {T,V} = AbstractMatrix(F)
# Matrix(F::UL{T,Tridiagonal{T,V}}) where {T,V} = Array(AbstractArray(F))
# Array(F::UL{T,Tridiagonal{T,V}}) where {T,V} = Matrix(F)


text:
export get_guild_audit_log

"""
    get_guild_audit_log(c::Client, guild::Integer; kwargs...) -> AuditLog

Get a [`Guild`](@ref)'s [`AuditLog`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/audit-log#get-guild-audit-log).
"""
function get_guild_audit_log(c::Client, guild::Integer; kwargs...)
    return Response{AuditLog}(c, :GET, "/guilds/$guild/audit-logs"; kwargs...)
end

ext: jl


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

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

    IoK8sApiCoreV1CephFSVolumeSource(;
        monitors=nothing,
        path=nothing,
        readOnly=nothing,
        secretFile=nothing,
        secretRef=nothing,
        user=nothing,
    )

    - monitors::Vector{String} : Required: Monitors is a collection of Ceph monitors More info: https://examples.k8s.io/volumes/cephfs/README.md#how-to-use-it
    - path::String : Optional: Used as the mounted root, rather than the full Ceph tree, default is /
    - readOnly::Bool : Optional: Defaults to false (read/write). ReadOnly here will force the ReadOnly setting in VolumeMounts. More info: https://examples.k8s.io/volumes/cephfs/README.md#how-to-use-it
    - secretFile::String : Optional: SecretFile is the path to key ring for User, default is /etc/ceph/user.secret More info: https://examples.k8s.io/volumes/cephfs/README.md#how-to-use-it
    - secretRef::IoK8sApiCoreV1LocalObjectReference : Optional: SecretRef is reference to the authentication secret for User, default is empty. More info: https://examples.k8s.io/volumes/cephfs/README.md#how-to-use-it
    - user::String : Optional: User is the rados user name, default is admin More info: https://examples.k8s.io/volumes/cephfs/README.md#how-to-use-it
"""
mutable struct IoK8sApiCoreV1CephFSVolumeSource <: SwaggerModel
    monitors::Any # spec type: Union{ Nothing, Vector{String} } # spec name: monitors
    path::Any # spec type: Union{ Nothing, String } # spec name: path
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    secretFile::Any # spec type: Union{ Nothing, String } # spec name: secretFile
    secretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1LocalObjectReference } # spec name: secretRef
    user::Any # spec type: Union{ Nothing, String } # spec name: user

    function IoK8sApiCoreV1CephFSVolumeSource(;monitors=nothing, path=nothing, readOnly=nothing, secretFile=nothing, secretRef=nothing, user=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1CephFSVolumeSource, Symbol("monitors"), monitors)
        setfield!(o, Symbol("monitors"), monitors)
        validate_property(IoK8sApiCoreV1CephFSVolumeSource, Symbol("path"), path)
        setfield!(o, Symbol("path"), path)
        validate_property(IoK8sApiCoreV1CephFSVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1CephFSVolumeSource, Symbol("secretFile"), secretFile)
        setfield!(o, Symbol("secretFile"), secretFile)
        validate_property(IoK8sApiCoreV1CephFSVolumeSource, Symbol("secretRef"), secretRef)
        setfield!(o, Symbol("secretRef"), secretRef)
        validate_property(IoK8sApiCoreV1CephFSVolumeSource, Symbol("user"), user)
        setfield!(o, Symbol("user"), user)
        o
    end
end # type IoK8sApiCoreV1CephFSVolumeSource

const _property_map_IoK8sApiCoreV1CephFSVolumeSource = Dict{Symbol,Symbol}(Symbol("monitors")=>Symbol("monitors"), Symbol("path")=>Symbol("path"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("secretFile")=>Symbol("secretFile"), Symbol("secretRef")=>Symbol("secretRef"), Symbol("user")=>Symbol("user"))
const _property_types_IoK8sApiCoreV1CephFSVolumeSource = Dict{Symbol,String}(Symbol("monitors")=>"Vector{String}", Symbol("path")=>"String", Symbol("readOnly")=>"Bool", Symbol("secretFile")=>"String", Symbol("secretRef")=>"IoK8sApiCoreV1LocalObjectReference", Symbol("user")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1CephFSVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1CephFSVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1CephFSVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1CephFSVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1CephFSVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1CephFSVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1CephFSVolumeSource)
    (getproperty(o, Symbol("monitors")) === nothing) && (return false)
    true
end

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


text:
function ec_pie_init()
    opt = ec_init()
    opt["tooltip"]["trigger"] = "item"
    opt["calculable"] = true

    return opt
end

function ec_pie{S<:String}(legendlist::Vector{S}, data::Vector; radius="55%", center = ["50%", "55%"])
    opt = ec_pie_init()
    opt["legend"]["data"] = legendlist

    m = length(data)
    data_dict = Dict(
    	"name" => "",
    	"type" => "pie",
    	"radius" => radius,
    	"center" => center,
		"data" => Array(Any, m)
		)

    for i = 1:m
        data_dict["data"][i] = Dict(
            "name" => legendlist[i],
            "value" => data[i]
            )
    end

    push!(opt["series"], data_dict)

    return opt
end

function ec_ring{S<:String}(legendlist::Vector{S}, data::Vector; radius=["50%", "70%"], center = ["50%", "55%"])
    ec_pie(legendlist, data; radius = radius, center = center)
end

function ec_rose_base{S<:String}(legendlist::Vector{S}, data::Vector; radius=["20%", "100%"], center = ["50%", "55%"])
	opt = ec_pie_init()
    opt["legend"]["data"] = legendlist

    m = length(data)
    data_dict = Dict(
    	"name" => "",
    	"type" => "pie",
    	"radius" => radius,
    	"center" => center,
		"data" => Array(Any, m)
		)

    for i = 1:m
        data_dict["data"][i] = Dict(
            "name" => legendlist[i],
            "value" => data[i]
            )
    end

    push!(opt["series"], data_dict)

    return opt
end

function ec_rose_area{S<:String}(legendlist::Vector{S}, data::Vector; radius=["20%", "70%"], center = ["50%", "55%"])
	opt = ec_rose_base(legendlist, data; radius = radius, center = center)
	opt["series"][1]["roseType"] = "area"

	return opt
end

function ec_rose_radius{S<:String}(legendlist::Vector{S}, data::Vector; radius=["20%", "70%"], center = ["50%", "55%"])
	opt = ec_rose_base(legendlist, data; radius = radius, center = center)
	opt["series"][1]["roseType"] = "radius"

	return opt
end


text:
function highpass_filter(X::Array{Float64,1},sampling_rate::Float64, cutoff::Float64=300.0;method=Butterworth(4))
  filtfilt(digitalfilter(Highpass(cutoff;fs=sampling_rate), method), X)
end

function get_threshold(X::Array{Float64,1},θ=6.0)
  n = length(X)
  μ = mean(X)
  σ = std(X)::Float64 #because of type instability in std
  l = μ - θ*σ
  u = μ + θ*σ
  σ2 = 0.0
  for x in X
    if l < x < u
      σ2 +=  x*x
    end
  end
  σ2 = sqrt(σ2/n - μ*μ)
  μ, σ2
end

function extract_spikes(X::Array{Float64,1};μ::Float64=NaN, σ::Float64=NaN,nq::Tuple{Int64, Int64}=(20,40),θ=6.0,only_negative=true)
  if isnan(σ) || isnan(μ)
    μ, σ = get_threshold(X,θ)
  end
  n = length(X)
  l = μ - θ*σ
  u = μ + θ*σ
  pidx = Array{Int64}(0)
  xmin = Inf
  xmax = -Inf
  j = 1
  in_peak = false
  _idx = 0
  while j <= n
    #check whether points are outside μ ± σ
    #fast forward points that are within noise levels
    x = X[j]
    if x < l
      if x < xmin
        xmin = x
        _idx = j
      end
      in_peak = true
  elseif (x > u ) && !only_negative
      if x > xmax
        xmax = x
        _idx = j
      end
      in_peak = true
    else
      if in_peak == true #just exited a peak; reset
        push!(pidx, _idx)
      end
      xmin = Inf
      xmax = -Inf
      in_peak = false
    end
    j += 1
  end
  pidx, waveforms = extract_spikes(X, pidx,nq)
  pidx, waveforms
end

function extract_spikes(X::Array{Float64,1}, idx::Array{Int64,1}, n::Tuple{Int64, Int64}=(20,40))
  nt = sum(n)
  nb,na = n
  idxl = find(x->x>nb, idx)
  nx = length(idxl)
  waveforms = zeros(nt, nx)
  for i in 1:nx
     _idx = idx[idxl[i]]
    for j in -nb:na-1
      waveforms[j+nb+1,i] = X[_idx + j]
    end
  end
  idx[idxl], waveforms
end


text:
immutable Logistic <: ContinuousUnivariateDistribution
    location::Real
    scale::Real
    function Logistic(l::Real, s::Real)
    	if s > 0.0
    		new(float64(l), float64(s))
    	else
    		error("scale must be positive")
    	end
	end
end

Logistic(l::Real) = Logistic(l, 1.0)
Logistic()  = Logistic(0.0, 1.0)

@_jl_dist_2p Logistic logis

entropy(d::Logistic) = log(d.scale) + 2.0

insupport(d::Logistic, x::Number) = isreal(x) && isfinite(x)

kurtosis(d::Logistic) = 1.2

mean(d::Logistic) = d.location

median(d::Logistic) = d.location

modes(d::Logistic) = [d.location]

skewness(d::Logistic) = 0.0

std(d::Logistic) = pi * d.scale / sqrt(3.0)

var(d::Logistic) = (pi * d.scale)^2 / 3.0


text:
module Vis

using GLMakie: GLMakie
using ColorTypes: RGB
using FromFile: @from
@from "../Geometry/Geometry.jl" import Geometry:
    Surface,
    CSGTree,
    CSGGenerator,
    BoundingBox,
    TriangleMesh,
    makemesh,
    makiemesh

global current_main_scene = nothing
global current_layout_scene = nothing
global current_3d_scene = nothing
global current_mode = nothing
# modes:    nothing, :default  -> Original Vis beheviour
#           :pluto             -> support pluto notebooks
#           :docs              -> support documenter figures

# added the following 2 functions to allow us to hack the drawing mechanisim while in a pluto notebook
set_current_main_scene(scene) = (global current_main_scene = scene)
set_current_3d_scene(lscene) = (global current_3d_scene = lscene)

get_current_mode() = begin
    global current_mode
    return current_mode
end
set_current_mode(mode) = (global current_mode = mode)

"""
    scene(resolution = (1000, 1000))

Create a new GLMakie scene with the given resolution including control buttons.
"""
function scene(resolution = (1000, 1000))
    @assert resolution[1] > 0 && resolution[2] > 0

    scene, layout = GLMakie.layoutscene(resolution = resolution)
    global current_main_scene = scene
    global current_layout_scene = layout
    lscene =
        layout[1, 1] = GLMakie.LScene(
            scene,
            scenekw = (
                camera = GLMakie.cam3d_cad!,
                axis_type = GLMakie.axis3d!,
                raw = false,
            ),
        )
    global current_3d_scene = lscene

    # in these modes we want to skip the creation of the utility buttons as these modes are not interactive
    if (get_current_mode() == :pluto || get_current_mode() == :docs)
        return scene, lscene
    end

    threedbutton = GLMakie.Button(
        scene,
        label = "3D",
        buttoncolor = RGB(0.8, 0.8, 0.8),
        height = 40,
        width = 80,
    )
    twodxbutton = GLMakie.Button(
        scene,
        label = "2D-x",
        buttoncolor = RGB(0.8, 0.8, 0.8),
        height = 40,
        width = 80,
    )
    twodybutton = GLMakie.Button(
        scene,
        label = "2D-y",
        buttoncolor = RGB(0.8, 0.8, 0.8),
        height = 40,
        width = 80,
    )
    savebutton = GLMakie.Button(
        scene,
        label = "Screenshot",
        buttoncolor = RGB(0.8, 0.8, 0.8),
        height = 40,
        width = 160,
    )

    GLMakie.on(threedbutton.clicks) do nclicks
        make3d(lscene)
        return yield()
    end

    GLMakie.on(twodybutton.clicks) do nclicks
        make2dy(lscene)
        return yield()
    end

    GLMakie.on(twodxbutton.clicks) do nclicks
        make2dx(lscene)
        return yield()
    end

    GLMakie.on(savebutton.clicks) do nclicks
        Vis.save("screenshot.png")
        return yield()
    end

    layout[2, 1] = GLMakie.grid!(
        hcat(threedbutton, twodxbutton, twodybutton, savebutton),
        tellwidth = false,
        tellheight = true,
    )
    return scene, lscene
end

function make3d(scene::GLMakie.LScene = current_3d_scene)
    s = scene.scene
    # use 3d camera
    GLMakie.cam3d_cad!(s)
    # reset scene rotation
    s.transformation.rotation[] = GLMakie.Quaternion(0.0, 0.0, 0.0, 1.0)
    # show all the axis ticks
    s[GLMakie.OldAxis].attributes.showticks[] = (true, true, true)
    # reset tick and axis label rotation and position
    rot = (
        GLMakie.Quaternion(0.0, 0.0, -0.7071067811865476, -0.7071067811865475),
        GLMakie.Quaternion(0.0, 0.0, 1.0, 0.0),
        GLMakie.Quaternion(0.0, 0.7071067811865475, 0.7071067811865476, 0.0),
    )
    s[GLMakie.OldAxis].attributes.ticks.rotation[] = rot
    s[GLMakie.OldAxis].attributes.names.rotation[] = rot
    s[GLMakie.OldAxis].attributes.ticks.align =
        ((:left, :center), (:right, :center), (:right, :center))
    s[GLMakie.OldAxis].attributes.names.align =
        ((:left, :center), (:right, :center), (:right, :center))
    # reset scene limits to automatic
    s.limits[] = GLMakie.Makie.Automatic()
    return GLMakie.update!(s)
end

function make2dy(scene::GLMakie.LScene = current_3d_scene)
    s = scene.scene
    # use 2d camera
    GLMakie.cam2d!(s)

    scene_transform = GLMakie.qrotation(GLMakie.Vec3f(0, 1, 0), 0.5pi)
    scene_transform_inv = GLMakie.qrotation(GLMakie.Vec3f(0, 1, 0), -0.5pi)    # to use with the ticks and names

    # set rotation to look onto yz plane
    s.transformation.rotation[] = scene_transform
    # hide x ticks

    # there is a bug in GLMakie 0.14.2 which cause an exception setting the X showticks to false.
    # we work wround it by making sure the labels we want to turn off are ortogonal to the view direction
    # s[GLMakie.OldAxis].attributes.showticks[] = (false, true, true)
    s[GLMakie.OldAxis].attributes.showticks[] = (true, true, true)

    # set tick and axis label rotation and position
    s[GLMakie.OldAxis].attributes.ticks.rotation[] =
        (0.0, scene_transform_inv, scene_transform_inv)
    s[GLMakie.OldAxis].attributes.names.rotation[] =
        s[GLMakie.OldAxis].attributes.ticks.rotation[]
    s[GLMakie.OldAxis].attributes.ticks.align =
        ((:right, :center), (:right, :center), (:center, :right))
    s[GLMakie.OldAxis].attributes.names.align =
        ((:left, :center), (:left, :center), (:center, :left))
    # update the scene limits automatically to get true reference values
    s.limits[] = GLMakie.Makie.Automatic()
    GLMakie.update_limits!(s)
    # manually set the scene limits to draw the axes correctly
    o, w = GLMakie.origin(s.data_limits[]), GLMakie.widths(s.data_limits[])
    s.limits[] = GLMakie.FRect3D((1000.0f0, o[2], o[3]), (w[2], w[2], w[3]))
    # set the eye (i.e. light) position to behind the camera
    s.camera.eyeposition[] = (0, 0, -100)
    return GLMakie.update!(s)
end

function make2dx(scene::GLMakie.LScene = current_3d_scene)
    s = scene.scene
    # use 2d camera
    GLMakie.cam2d!(s)

    scene_transform =
        GLMakie.qrotation(GLMakie.Vec3f(0, 0, 1), 0.5pi) *
        GLMakie.qrotation(GLMakie.Vec3f(1, 0, 0), 0.5pi)
    scene_transform_inv =
        GLMakie.qrotation(GLMakie.Vec3f(1, 0, 0), -0.5pi) *
        GLMakie.qrotation(GLMakie.Vec3f(0, 0, 1), -0.5pi)

    # set rotation to look onto yz plane
    s.transformation.rotation[] = scene_transform
    # hide y ticks

    # there is a bug in GLMakie 0.14.2 which cause an exception setting the X showticks to false.
    # we work wround it by making sure the labels we want to turn off are ortogonal to the view direction
    s[GLMakie.OldAxis].attributes.showticks[] = (true, false, true)

    # set tick and axis label rotation and position
    s[GLMakie.OldAxis].attributes.ticks.rotation[] =
        (scene_transform_inv, 0.0, scene_transform_inv)
    s[GLMakie.OldAxis].attributes.names.rotation[] =
        s[GLMakie.OldAxis].attributes.ticks.rotation[]
    s[GLMakie.OldAxis].attributes.ticks.align =
        ((:right, :center), (:right, :center), (:center, :center))
    s[GLMakie.OldAxis].attributes.names.align =
        ((:left, :center), (:right, :center), (:center, :center))
    # update the scene limits automatically to get true reference values
    s.limits[] = GLMakie.Makie.Automatic()
    GLMakie.update_limits!(s)
    # manually set the scene limits to draw the axes correctly
    o, w = GLMakie.origin(s.data_limits[]), GLMakie.widths(s.data_limits[])
    s.limits[] = GLMakie.FRect3D((o[1], -1000.0f0, o[3]), (w[1], w[1], w[3]))
    # set the eye (i.e. light) position to behind the camera
    s.camera.eyeposition[] = (0, 0, -100)
    return GLMakie.update!(s)
end

"""
    draw(ob; resolution = (1000, 1000), kwargs...)

Draw an object in a new scene.
`kwargs` depends on the object type.
"""
function draw(ob; resolution = (1000, 1000), kwargs...)
    scene, lscene = Vis.scene(resolution)
    draw!(lscene, ob; kwargs...)
    display(scene)

    if (get_current_mode() == :pluto || get_current_mode() == :docs)
        return scene
    end
end

"""
    draw!([scene = currentscene], ob; kwargs...)

Draw an object in an existing scene.
`kwargs` depends on the object type.
"""
function draw!(ob; kwargs...)
    if current_3d_scene === nothing
        scene, lscene = Vis.scene()
    else
        scene = current_main_scene
        lscene = current_3d_scene
    end

    GLMakie.cameracontrols(lscene.scene).attributes.mouse_rotationspeed[] =
        0.0001f0 #doesn't seem to do anything
    draw!(lscene, ob; kwargs...)
    display(scene)

    if (get_current_mode() == :pluto || get_current_mode() == :docs)
        return scene
    end
end

## GEOMETRY

"""
    draw!(scene::GLMakie.LScene, surf::Surface{T}; numdivisions = 20, normals = false, normalcolor = :blue, kwargs...)

Transforms `surf` into a mesh using [`makemesh`](@ref) and draws the result.
`normals` of the surface can be drawn at evenly sampled points with provided `normalcolor`.
`numdivisions` determines the resolution with which the mesh is triangulated.
`kwargs` is passed on to the [`TriangleMesh`](@ref) drawing function.
"""
function draw!(
    scene::GLMakie.LScene,
    surf::Surface{T};
    numdivisions::Int = 30,
    normals::Bool = false,
    normalcolor = :blue,
    kwargs...,
) where {T <: Real}
    mesh = makemesh(surf, numdivisions)
    if nothing === mesh
        return
    end
    draw!(scene, mesh; kwargs..., normals = false)
    if normals
        ndirs = GLMakie.Point3f.(samplesurface(surf, normal, numdivisions ÷ 10))
        norigins =
            GLMakie.Point3f.(samplesurface(surf, point, numdivisions ÷ 10))
        GLMakie.arrows!(
            scene,
            norigins,
            ndirs,
            arrowsize = 0.2,
            arrowcolor = normalcolor,
            linecolor = normalcolor,
            linewidth = 2,
        )
    end
end

"""
    draw!(scene::GLMakie.LScene, tmesh::TriangleMesh{T}; linewidth = 3, shaded = true, wireframe = false, color = :orange, normals = false, normalcolor = :blue, transparency = false, kwargs...)

Draw a [`TriangleMesh`](@ref), optionially with a visible `wireframe`. `kwargs` are passed on to [`GLMakie.mesh`](http://makie.juliaplots.org/stable/plotting_functions.html#mesh).
"""
function draw!(
    scene::GLMakie.LScene,
    tmesh::TriangleMesh{T};
    linewidth = 3,
    shaded::Bool = true,
    wireframe::Bool = false,
    color = :orange,
    normals::Bool = false,
    normalcolor = :blue,
    transparency::Bool = false,
    kwargs...,
) where {T <: Real}
    points, indices = makiemesh(tmesh)
    if length(points) > 0 && length(indices) > 0
        GLMakie.mesh!(
            scene,
            points,
            indices;
            kwargs...,
            color = color,
            shading = shaded,
            transparency = transparency,
            visible = shaded,
        )
        if wireframe
            mesh = scene.scene[end][1]
            if shaded
                GLMakie.wireframe!(
                    scene,
                    mesh,
                    color = (:black, 0.1),
                    linewidth = linewidth,
                )
            else
                GLMakie.wireframe!(
                    scene,
                    mesh,
                    color = color,
                    linewidth = linewidth,
                )
            end
        end
    end
    if normals
        @warn "Normals being drawn from triangulated mesh, precision may be low"
        norigins =
            [GLMakie.Point3f(centroid(t)) for t in tmesh.triangles[1:10:end]]
        ndirs = [GLMakie.Point3f(normal(t)) for t in tmesh.triangles[1:10:end]]
        if length(norigins) > 0
            GLMakie.arrows!(
                scene,
                norigins,
                ndirs,
                arrowsize = 0.2,
                arrowcolor = normalcolor,
                linecolor = normalcolor,
                linewidth = 2,
            )
        end
    end
end

"""
    draw!(scene::GLMakie.LScene, meshes::Vararg{S}; colors::Bool = false, kwargs...) where {T<:Real,S<:Union{TriangleMesh{T},Surface{T}}}

Draw a series of [`TriangleMesh`](@ref) or [`Surface`](@ref) objects, if `colors` is true then each mesh will be colored automatically with a diverse series of colors.
`kwargs` are is passed on to the drawing function for each element.
"""
function draw!(
    scene::GLMakie.LScene,
    meshes::Vararg{S};
    colors::Bool = false,
    kwargs...,
) where {T <: Real, S <: Union{TriangleMesh{T}, Surface{T}}}
    for i in 1:length(meshes)
        if colors
            col = indexedcolor2(i)
        else
            col = :orange
        end
        draw!(scene, meshes[i]; kwargs..., color = col)
    end
end

function draw(
    meshes::Vararg{S};
    kwargs...,
) where {T <: Real, S <: Union{TriangleMesh{T}, Surface{T}}}
    scene, lscene = Vis.scene()
    draw!(lscene, meshes...; kwargs...)
    GLMakie.display(scene)
    if (get_current_mode() == :pluto || get_current_mode() == :docs)
        return scene
    end
end

"""
    draw!(scene::GLMakie.LScene, csg::Union{CSGTree,CSGGenerator}; numdivisions::Int = 20, kwargs...)

Convert a CSG object ([`CSGTree`](@ref) or [`CSGGenerator`](@ref)) to a mesh using [`makemesh`](@ref) with resolution set by `numdivisions` and draw the resulting [`TriangleMesh`](@ref).
"""
function draw!(
    scene::GLMakie.LScene,
    csg::CSGTree{T};
    numdivisions::Int = 30,
    kwargs...,
) where {T <: Real}
    return draw!(scene, makemesh(csg, numdivisions); kwargs...)
end
function draw!(
    scene::GLMakie.LScene,
    csg::CSGGenerator{T};
    kwargs...,
) where {T <: Real}
    return draw!(scene, csg(); kwargs...)
end

"""
    draw!(scene::GLMakie.LScene, bbox::BoundingBox{T}; kwargs...)

Draw a [`BoundingBox`](@ref) as a wireframe, ie series of lines.
"""
function draw!(
    scene::GLMakie.LScene,
    bbox::BoundingBox{T};
    kwargs...,
) where {T <: Real}
    p1 = SVector{3, T}(bbox.xmin, bbox.ymin, bbox.zmin)
    p2 = SVector{3, T}(bbox.xmin, bbox.ymax, bbox.zmin)
    p3 = SVector{3, T}(bbox.xmin, bbox.ymax, bbox.zmax)
    p4 = SVector{3, T}(bbox.xmin, bbox.ymin, bbox.zmax)
    p5 = SVector{3, T}(bbox.xmax, bbox.ymin, bbox.zmin)
    p6 = SVector{3, T}(bbox.xmax, bbox.ymax, bbox.zmin)
    p7 = SVector{3, T}(bbox.xmax, bbox.ymax, bbox.zmax)
    p8 = SVector{3, T}(bbox.xmax, bbox.ymin, bbox.zmax)
    return GLMakie.linesegments!(
        scene,
        [
            p1,
            p2,
            p2,
            p3,
            p3,
            p4,
            p4,
            p1,
            p1,
            p5,
            p2,
            p6,
            p3,
            p7,
            p4,
            p8,
            p5,
            p6,
            p6,
            p7,
            p7,
            p8,
            p8,
            p5,
        ];
        kwargs...,
    )
end

"""
    save(path::String)

Save the current GLMakie scene to an image file.
"""
function save(path::String)
    # save closes the window so just display it again as a work-around
    # for some reason the size isn't maintained automatically so we just reset it manually
    size = GLMakie.size(current_main_scene)
    GLMakie.save(path, current_3d_scene.scene)
    GLMakie.resize!(current_main_scene, size)
    return display(current_main_scene)
end
function save(::Nothing) end

end # module


text:
"""
Abstract supertype for parameters. Theses are wrappers for model parameter values and
metadata that are returned from [`params`](@ref), and used in
`getfield/setfield/getpropery/setproperty` methods and to generate the Tables.jl interface.
They are stripped from the model with [`stripparams`](@ref).

An `AbstractParam` must define a `Base.parent` method that returns a `NamedTuple`, and a
constructor that accepts a `NamedTuple`. It must have a `val` property, and should use
`checkhasval` in its constructor.
"""
abstract type AbstractParam{T} <: AbstractNumbers.AbstractNumber{T} end

function ConstructionBase.setproperties(p::P, patch::NamedTuple) where P <: AbstractParam
    fields = ConstructionBase.setproperties(parent(p), patch)
    P.name.wrapper(fields)
end

@inline withunits(m, args...) = map(p -> withunits(p, args...), params(m))
@inline function withunits(p::AbstractParam, fn::Symbol=:val)
    _applyunits(*, getproperty(p, fn), get(p, :units, nothing))
end

@inline stripunits(m, xs) = map(stripunits, params(m), xs)
@inline function stripunits(p::AbstractParam, x)
    _applyunits(/, x, get(p, :units, nothing))
end

# Param might have `nothing` for units
@inline _applyunits(f, x, units) = f(x, units)
@inline _applyunits(f, x, ::Nothing) = x
@inline _applyunits(f, xs::Tuple, units) = map(x -> f(x, units), xs)
@inline _applyunits(f, xs::Tuple, units::Nothing) = xs
@inline _applyunits(f, ::Nothing, units) = nothing
@inline _applyunits(f, ::Nothing, ::Nothing) = nothing

# Base NamedTuple-like interface
Base.keys(p::AbstractParam) = keys(parent(p))
# Base.values has the potential to be confusing, as we
# have a val field in Param.  Not sure what to do about this.
Base.values(p::AbstractParam) = values(parent(p))
@inline Base.propertynames(p::AbstractParam) = propertynames(parent(p))
@inline Base.getproperty(p::AbstractParam, x::Symbol) = getproperty(parent(p), x)
@inline Base.get(p::AbstractParam, key::Symbol, default) = get(parent(p), key, default)
@inline Base.getindex(p::AbstractParam, i) = getindex(parent(p), i)

# AbstractNumber interface
Base.convert(::Type{Number}, x::AbstractParam) = AbstractNumbers.number(x)
Base.convert(::Type{P}, x::P) where {P<:AbstractParam} = x
AbstractNumbers.number(p::AbstractParam) = withunits(p)
AbstractNumbers.basetype(::Type{<:AbstractParam{T}}) where T = T
AbstractNumbers.like(::Type{<:AbstractParam}, x) = x

# Flatten.jl defaults defined here: AbstractParam needs to be defined first
const SELECT = AbstractParam
const IGNORE = AbstractDict # What else to blacklist?

# Concrete implementation

"""
    Param(p::NamedTuple)
    Param(; kw...)
    Param(val)

A wrapper type that lets you extract model parameters and metadata about the model like
bounding val, units priors, or anything else you want to attach.

The first argument is assigned to the `val` field, and if only keyword arguments are used,
`val`, must be one of them. `val` is used as the number val if the model us run
without stripping out the `Param` fields. `stripparams` also takes only the `:val` field.
"""
struct Param{T,P<:NamedTuple} <: AbstractParam{T}
    parent::P
end
Param(nt::NT) where {NT<:NamedTuple} = begin
    _checkhasval(nt)
    Param{typeof(nt.val),NT}(nt)
end
Param(val; kwargs...) = Param((; val=val, kwargs...))
Param(; kwargs...) = Param((; kwargs...))

Base.parent(p::Param) = getfield(p, :parent)

# Methods for objects that hold params
params(x) = Flatten.flatten(x, SELECT, IGNORE)
stripparams(x) = hasparam(x) ? Flatten.reconstruct(x, withunits(x), SELECT, IGNORE) : x

# Utils
hasparam(obj) = length(params(obj)) > 0

_checkhasval(nt::NamedTuple{Keys}) where {Keys} = first(Keys) == :val || _novalerror(nt)
# @noinline avoids allocations unless there is actually an error
@noinline _novalerror(nt) = throw(ArgumentError("First field of Param must be :val"))


text:
# a simple interface to convert to and from DataFrames.
# DataFrames directly corresponds to DictArray of dimension 1 with symbol field names.

import DataFrames: DataFrame, DataArray, PooledDataArray, Index

dataframe_to_dictarray(df::DataFrame) = begin
  colindex = df.colindex
  columns = df.columns
  newcolumns = map(columns) do column
    if isa(column, DataArray)
      simplify_array(create_naarray(column))
    elseif isa(column, PooledDataArray)
      convert(EnumerationArray, column)
    else
      simplify_array(nalift(column))
    end
  end
  lookup = colindex.lookup
  newpairs = map(colindex.names) do name
    name => newcolumns[lookup[name]]
  end
  DictArray(newpairs...)
end

create_naarray{T}(arr::DataArray{T}) = begin
  data = arr.data
  na = arr.na
  map(data, na) do elem, isna
    isna ? Nullable{T}() : Nullable(elem)
  end
end
create_dataarray{T}(arr::AbstractArray{Nullable{T}}) = begin
  isna_array0 = isna(arr)
  isna_array = if isa(isna_array0, BitArray)
    isna_array0
  else
    result = BitArray(size(isna_array0))
    copy!(result, isna_array0)
    result
  end
  value_array = Array{T}(size(arr))
  arbitrary_array = similar([],T,1)
  anyelem = if isassigned(arbitrary_array,1)
    arbitrary_array[1]
  else
    find_nonnull(arr)
  end
  if !isempty(value_array)
    map!(value_array, arr) do elem
      if isnull(elem)
        anyelem
      else
        elem.value
      end
    end
  end
  DataArray(value_array, isna_array)
end

# promote to Nullable{Any} if the element type is just Nullable.
create_dataarray(arr::AbstractArray{Nullable}) = create_dataarray(convert(AbstractArray{Nullable{Any}}, arr))

find_nonnull(x::AbstractArray) = begin
  for elem in x
    if !isnull(elem)
      return elem.value
    end
  end
end

dictarray_to_dataframe{T}(arr::DictArray{T,1}) = begin
  keys = [Symbol(k) for k in arr.data.keys]
  colindex = Index(Dict(map(ik->ik[2]=>ik[1], enumerate(keys))), keys)
  columns = Any[create_dataarray(v) for v in arr.data.values]
  DataFrame(columns, colindex)
end

"""

`convert(::Type{DictArray}, df::DataFrame)` converts a `DataFrame` into `DictArray`.

"""
Base.convert(::Type{DictArray}, df::DataFrame) = dataframe_to_dictarray(df)

"""

`convert(::Type{LabeledArray}, df::DataFrame)` converts a `DataFrame` into `LabeledArray` simply by wrapping `convert(DictArray, df)` by `LabeledArray`.

"""
Base.convert(::Type{LabeledArray}, df::DataFrame) = LabeledArray(dataframe_to_dictarray(df))

"""

`convert(::Type{DataFrame}, arr::DictArray)` converts a `DictArray` into a `DataFrame`. If the dimensions of `arr` are greater than 1, `arr` is first flattend into 1 dimension using `collapse_axes`, and then converted into a `DataFrame`.

"""
Base.convert(::Type{DataFrame}, arr::DictArray) = dictarray_to_dataframe(collapse_axes(arr, 1, ndims(arr)))

"""

`convert(::Type{DataFrame}, arr::LabeledArray)` converts a `LabeledArray` into a `DataFrame` by first creating a `DictArray` by broadcasting all axes, and then convert that `DictArray` into a `DataFrame`.

"""
Base.convert(::Type{DataFrame}, arr::LabeledArray) = dictarray_to_dataframe((darr=selectfields(arr, allfieldnames(arr)...);collapse_axes(darr,1,ndims(darr))))

"""

`convert(::Type{EnumerationArray}, arr::PooledDataArray)` converts a `PooledDataArray` into an `EnumerationArray`.

"""
Base.convert(::Type{EnumerationArray}, arr::PooledDataArray) = EnumerationArray((arr.refs, arr.pool))


text:
#
# timedAcquisition.jl -
#
# This example demonstrates the usage of the `readTimed` method which is
# returning the time at which the exposition ended.
#
# Please refer to the "simpleAcquisition.jl" example if you're looking for a
# simpler example to acquire images.
#

module TimedAcquisition

# In this example, we want to use the low-level interface as it is the closest
# to the C API.
import NuvuCameras: NC

function main()
    # Opens the acquisition channel using the automatic detection and 4 loop
    # buffers (recommended).
    myCam = NC.open(NC.Cam, NC.AUTO_UNIT, NC.AUTO_CHANNEL, 4)

    # For this example we will not be checking for the available readout modes
    # on the camera and simply assume that the first one is valid.
    NC.setReadoutMode(myCam, 1)

    print("The camera just finished it's initialization, we'll now reset the ",
          "timer and grab a 'timed' image.\n\n")
    readline(STDIN)

    # Open the shutter for the acquisition.
    NC.setShutterMode(myCam, NC.OPEN)

    # Start the timer.
    NC.resetTimer(myCam, 0)

    # Launches an acquisition on the frame grabber and requests an image from
    # the camera (this function does not wait for the acquisition to be
    # complete before returning).
    NC.start(myCam, 1)

    # Reads the image when it's received and records the time; if a timeout
    # occurs an error code will be returned.
    myImage, imageTime = NC.readTimed(myCam) # FIXME: deprecated

    # Close the shutter, now that the acquistion is complete.
    NC.setShutterMode(myCam, NC.CLOSE)

    @printf("The time at which the exposure stopped is %f ms ", imageTime)
    print("after the reset of the timer (done just before requesting an ",
          "image from the camera).\n\n")

    NC.saveImage(myCam, myNcImage, "FirstImage", NC.TIF,
                 "This is the first image grabbed in the timed example", true)

    printf("You just acquired one image, by pressing Enter you will ",
           "acquired 15 new images\n")
    readline(STDIN)

    # Open the shutter for the acquisition.
    NC.setShutterMode(myCam, NC.OPEN)

    # Launches 15 acquisitions on the frame grabber and requests 15 images from
    # the camera (this function does not wait for the acquisition to be
    # complete before returning).
    NC.start(myCam, 15)

    # Loop in which the acquired images are read.
    for i in 1:15

 	# Reads the image received, if a timeout occurs an error code will be
 	# returned.
        myImage, imageTime = NC.readTimed(myCam) # FIXME: deprecated

 	# Checks if an overrun occured on the last image (implying that the
 	# buffer we are reading has been overwritten prior to this "read"
 	# call).
        overrun = NC.getOverrun(myCam)

 	@printf("The time at which the exposition stopped is of %fms ",
                imageTime)
        printf("after the reset of the timer (done just before requesting ",
               "the first image from the camera).\n\n")

 	# If you wish to proceed with some processing on each image read, it
 	# should be done here (in this example we will only save the image).

 	# Saves each image acquired, at the end of each name the loop index
 	# will be added.
 	NC.saveImage(myCam, myNcImage, "Image_$i", NC.TIF,
                     "This is one of the images grabbed in the timed example",
                     true)
    end

    # Close the shutter, now that the acquistion is complete.
    NC.setShutterMode(myCam, NC.CLOSE)

    # Closes the acquisition channel no longer in use.
    NC.CamClose(myCam)

    if is_windows()
	print("\nHit the Return key to close the command window.")
        flush(STDOUT)
        readline(STDIN)
    else
	println("The program ended.")
    end
end

end # module


text:
using SHIPs, Printf, Test, LinearAlgebra
using SHIPs: PolyTransform, rbasis, eval_basis, eval_basis_d
using SHIPs.JacobiPolys: Jacobi
using PyPlot

trans = PolyTransform(1, 1)
B2 = rbasis(4, trans, 2, 0.5, 2.0)

rr = range(0, 2.5, length=100)
BB = zeros(100, length(B2))
for (i, r) in enumerate(rr)
   BB[i, :] = eval_basis(B2, r)
end

##
clf()
for i = 1:length(B2)
   plot(rr, BB[:,i], label = "P_$i")
end
legend()
xlabel("r")
ylabel("P_i(r) = J_i(x) * w(x)")

PyPlot.display_figs()

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
using NeighbourLists
using JuLIP: Atoms, cutoff, neighbourlist, JVecF, positions, set_positions!
using JuLIP.MLIPs: IPCollection, IPSuperBasis
using JuLIP.Potentials: site_energy, site_energy_d
import JuLIP
import SHIPs
using SHIPs: SHIPBasis, eval_basis, eval_basis_d
using LinearAlgebra: qr, norm, cond
using LowRankApprox: pqrfact
using PrettyTables: pretty_table, ft_printf

_hcat(X::AbstractVector) = hcat(X...)

_site_energy_d(args...) = hcat( site_energy_d(args...)... )

"""
l0 : site that we are perturbing
 l : site at which the site energy is evaluated
 i : direction of the perturbation (E1, E2, E3)
 h : finite-difference step
"""
function _site_energy_d2h(basis, at::Atoms,
                           l0::Integer, l::Integer, i::Integer, h::Real)
   X = positions(at)
   X[l0] += h * evec(i)
   set_positions!(at, X)
   dBp = _site_energy_d(basis, at, l)
   X[l0] -= 2*h * evec(i)
   set_positions!(at, X)
   dBm = _site_energy_d(basis, at, l)
   X[l0] += h * evec(i)
   set_positions!(at, X)
   return (dBp - dBm) / (2*h)
end

"""
l0    : first site that we are perturbing
 l    : site at which the site energy is evaluated
 k    : second site that we are perturbin
 i_l0 : direction of the perturbation of l0
 i_k  : direction of the perturbation of k
 h    : finite-difference step
"""
function _site_energy_d3h(basis, at::Atoms,
                           l0::Integer, l::Integer, k::Integer,
                           i_l0::Integer, i_k::Integer, h::Real)
   X = positions(at)
   # -------------------------
   X[l0] += h * evec(i_l0)
   X[k] += h * evec(i_k)
   set_positions!(at, X)
   dBpp = _site_energy_d(basis, at, l)
   X[l0] -= h * evec(i_l0)
   X[k] -= h * evec(i_k)
   # -------------------------
   X[l0] += h * evec(i_l0)
   X[k] -= h * evec(i_k)
   set_positions!(at, X)
   dBpm = _site_energy_d(basis, at, l)
   X[l0] -= h * evec(i_l0)
   X[k] += h * evec(i_k)
   # -------------------------
   X[l0] -= h * evec(i_l0)
   X[k] += h * evec(i_k)
   set_positions!(at, X)
   dBmp = _site_energy_d(basis, at, l)
   X[l0] += h * evec(i_l0)
   X[k] -= h * evec(i_k)
   # -------------------------
   X[l0] -= h * evec(i_l0)
   X[k] -= h * evec(i_k)
   set_positions!(at, X)
   dBmm = _site_energy_d(basis, at, l)
   X[l0] += h * evec(i_l0)
   X[k] += h * evec(i_k)
   # -------------------------
   set_positions!(at, X)
   return (dBpp - dBpm - dBmp + dBmm) / (2*h)^2
end

function _asm_(djEs, djB, w, nn)
   Y = zeros(length(nn) * 3)
   A = zeros(3*length(nn), size(djB, 2))
   for in = 1:length(nn)
      rows = 3*(in-1) .+ (1:3)
      Y[3*(in-1) .+ (1:3)] .= w * djEs[in]
      for iB = 1:size(djB, 2)
         A[rows, iB] .= w * djB[in, iB]
      end
   end
   return A, Y
end

function assemble_lsq(::Val{:dEs}, basis, config, at, w, key)
   # import data
   dEs = config[key]
   l0 = 1
   @assert length(dEs) == length(at)
   # assemble basis
   dB = _site_energy_d(basis, at, l0)
   @assert size(dB) == (length(at), length(basis))
   # -------
   return _asm_(dEs, dB, w, 1:length(at))
end

function assemble_lsq(::Val{:d2Esh}, basis, config, at, w, key)
   # import data
   d2Es = config[key]
   l = config["l"]    # the atom where we evaluate the site energy deriv
   i = config["i"]    # the direction of the perturbation of l0 = 1
   h = config["h"]
   l0 = 1
   @assert length(d2Es) == length(at)
   # assemble basis
   d2B = _site_energy_d2h(basis, at, l0, l, i, h)
   @assert size(d2B) == (length(at), length(basis))
   return _asm_(d2Es, d2B, w, 1:length(at))
end

function assemble_lsq(::Val{:d3Esh}, basis, config, at, w, key)
   # import data
   d3Es = config[key]
   h = config["h"]
   l = config["l"]    # the atom where we evaluate the site energy deriv
   l0 = 1              # first perturbed atom
   i_l0 = config["i"]    # the direction of the perturbation of l0 = 1
   k = config["k"]       # second perturbed atom
   i_k = config["j"]     # direction of perturbation of k
   @assert length(d3Es) == length(at)
   # assemble basis
   d3B = _site_energy_d3h(basis, at, l0, l, k, i_l0, i_k, h)
   @assert size(d3B) == (length(at), length(basis))
   # ----------------
   return _asm_(d3Es, d3B, w, 1:length(at))
end

function assemble_lsq(::Val{:Es}, basis, config, at, w, key)
   Y = [ w * config[key] ]
   A = Matrix( w * site_energy(basis, at, 1)' )
   @assert size(Y) == (1,)
   @assert size(A) == (1, length(basis))
   return A, Y
end

function assemble_lsq(basis, D::Dict, weights::Dict, key="train")

   at = Atoms(D["at"])::Atoms
   data = D["data"]

   # --------------- assemble local matrices  ---------------
   AA = Matrix{Float64}[]
   YY = Vector{Float64}[]
   DTs = String[]
   nrows = 0
   for dat in data
      # get some info about the local lsq system
      dt = dat["datatype"]
      if !haskey(weights, dt); continue; end
      w = weights[dt]
      if w == 0.0; continue; end
      # assemble the local lsq system
      A, Y = assemble_lsq(Val(Symbol(dt)), basis, dat, at, w, key)
      # store it in the AA, YY arrays
      if length(Y) != 0
         @assert size(A, 2) == length(basis)
         @assert size(A, 1) == length(Y)
         push!(AA, A)
         push!(YY, Y)
         push!(DTs, dt)   # remember the data-types to assemble the errors!
         nrows += length(Y)
      end
   end

   # --------------- assemble global matrix  ---------------
   A = zeros(nrows, length(basis))
   Y = zeros(nrows)
   DT = Vector{String}(undef, nrows)
   irow = 0
   for (a, y, dt) in zip(AA, YY, DTs)
      rows = irow .+ (1:size(a,1))
      A[rows, :] .= a
      Y[rows] .= y
      DT[rows] .= dt
      irow += length(rows)
   end

   return A, Y, DT
end

function lsqfit(basis, D::Dict, weights::Dict;
                verbose=true, key="train", kwargs...)
   A, Y, DT = assemble_lsq(basis, D, weights, key)
   qrA = pqrfact(A; rtol=1e-5, kwargs...)
   condA = cond(Matrix(qrA[:R]))
   verbose && @show condA
   c = qrA \ Y
   verbose && @show norm(c), norm(c, Inf)
   verbose && @show norm(A * c - Y) / norm(Y)
   return JuLIP.MLIPs.combine(basis, c),
          _fitinfo(basis, D, weights, kwargs, A, Y, DT, c)
end

function _fitinfo(basis, D, weights, kwargs, A, Y, DT, c)
   D = Dict("basis" => Dict(basis),
            "weights" => weights,
            "kwargs" => Dict(kwargs...),
            "rmse" => Dict{String, Any}(),
            "maxe" => Dict{String, Any}()  )
   errs = A * c - Y
   for dt in unique(DT)
      Idt = findall(DT .== dt)
      D["rmse"][dt] = norm(errs[Idt]) / norm(Y[Idt])
      D["maxe"][dt] = norm(errs[Idt], Inf) / norm(Y[Idt], Inf)
   end
   return D
end

function print_errors(fitinfo::Dict; fmt="%.3e")
   dts = ["Es", "dEs", "d2Esh", "d3Esh"]
   maxe = zeros(length(dts))
   rmse = zeros(length(dts))
   for (i, dt) in enumerate(dts)
      if haskey(fitinfo["maxe"], dt)
         maxe[i] = fitinfo["maxe"][dt]
      else
         maxe[i] = NaN
      end
      if haskey(fitinfo["rmse"], dt)
         rmse[i] = fitinfo["rmse"][dt]
      else
         rmse[i] = NaN
      end
   end
   fltfmt = ft_printf(fmt)[0]
   pretty_table( [dts maxe rmse], ["datatype", "maxe", "rmse"],
                 formatter = Dict(1 => (v,i) -> v,
                                  2 => fltfmt,
                                  3 => fltfmt) )
end


text:
# Test routine to test the infeasibility example from the OSQP Paper

# workspace()
# include("../../../src/Solver.jl")

using COSMO, Base.Test
maxIter = 3000
rng = Random.MersenneTwister(1313)

setUnNonAdaptive =   Settings(max_iter=maxIter,verbose=true,check_termination=1,checkInfeasibility=50,scaling = 0 ,scaleFunc=1,adaptive_rho=false,eps_abs=1e-7,eps_rel=1e-7)
setMeanNonAdaptive = Settings(max_iter=maxIter,verbose=false,check_termination=1,checkInfeasibility=50,scaling = 10,scaleFunc=2,adaptive_rho=false,eps_abs=1e-7,eps_rel=1e-7)

setUnAdaptive =   Settings(max_iter=maxIter,verbose=false,check_termination=1,checkInfeasibility=50,scaling = 0 ,scaleFunc=1,adaptive_rho=true,eps_abs=1e-7,eps_rel=1e-7)
setMeanAdaptive = Settings(max_iter=maxIter,verbose=false,check_termination=1,checkInfeasibility=50,scaling = 10,scaleFunc=2,adaptive_rho=true,eps_abs=1e-7,eps_rel=1e-7)

P = sparse([1. 0;0 0])
q = [1.;-1]

# # optimality
# a = 1
# u1 = 5
# u3 = 3
# A = sparse([1 a;-1 -a;1 0;-1 0; 0 1; 0 -1])
# b = [u1;0;3;-1;u3;-1]
# P = sparse([1 0;0 0])
# q = [1;-1]
# Kf = 0
# Kl = 6

# primal infeasibility
# case = "PRIMAL"
# A = sparse([1. 1;-1 -1;1 0;-1 0;0 1;0 -1])
# b = [0.;0;3;-1;3;-1]
# Kf = 0
# Kl = 6

# # dual infeasibility
# case = "DUAL"
# A = sparse([1 0;-1 0;1 0;-1 0; 0 -1])
# b = [2;0;3;-1;-1]
# Kf = 0
# Kl = 5

# primal and dual infeasiblity
case = "BOTH"
A = sparse([1 1;-1 -1;1 0;-1 0;;0 -1])
b = [0;0;3;-1;-1]
Kf = 0
Kl = 5

## primal infeasible test problem from MOSEK doc, i.e. min x1 s.t. x1 <= 5
# P = zeros(2,2)
# q = [1;0]
# A = [0 1;1 0]
# b = [0;5]
# Kf = 1
# Kl = 1
# Ks = []
# Kq = []

Kq = []
Ks = []
# define cone membership
K = COSMO.Cone(Kf,Kl,Kq,Ks)

# Solve with COSMO
res1,ws,δx1,δμ1 = COSMO.solve(P,q,A,b,K,setMeanAdaptive);
print("\n. δx1: $(δx1), δμ1: $(δμ1)")
# res2,nothing,δx2,δμ2 = COSMO.solve(P,q,A,b,K,setMeanNonAdaptive);
# print("\n. δx2: $(δx2), δμ2: $(δμ2)")
# res3,nothing,δx3,δμ3 = COSMO.solve(P,q,A,b,K,setUnAdaptive);
# print("\n. δx3: $(δx3), δμ3: $(δμ3)")
# res4,nothing,δx4,δμ4 = COSMO.solve(P,q,A,b,K,setMeanAdaptive);
# print("\n. δx4: $(δx4), δμ4: $(δμ4)\n")

# eps_pinf = 1e-4
# eps_dinf = 1e-4

# @testset "Test infeasibility conditions" begin
#     if case == "PRIMAL"
#       δy = -δμ1
#       @test norm(A'*δy,Inf)/norm(δy,Inf) <= eps_pinf
#       support_δy = b'*δy/norm(δy,Inf)
#       @test support_δy <= -eps_pinf
#     elseif case == "DUAL"
#       @test norm(P*δx1,Inf)/norm(δx1,Inf) <= eps_dinf
#       @test maximum(A*δx1)/norm(δx1,Inf) <= eps_dinf
#       @test q'*δx1/norm(δx1,Inf) <= -eps_dinf
#     elseif case == "BOTH"
#       @test norm(P*δx1,Inf)/norm(δx1,Inf) <= eps_dinf
#       @test maximum(A*δx1)/norm(δx1,Inf) <= eps_dinf
#       @test q'*δx1/norm(δx1,Inf) <= -eps_dinf
#       @test norm(A'*δy,Inf)/norm(δy,Inf) <= eps_pinf
#       support_δy = b'*δy/norm(δy,Inf)
#       @test support_δy <= -eps_pinf
#     end
# end


text:
"""
   mₑ : Electron rest mass (in eV)
"""
const mₑ = ustrip((ElectronMass * SpeedOfLightInVacuum^2) |> u"eV")

"""
    λₑ(E::Float64)

Wavelength of an electron in cm. (non-relativistic)

  * E in eV
"""
function λₑ(E::Float64)
    return 1.0e2 * ustrip(PlanckConstant / sqrt(2.0 * ElectronMass * ElementaryCharge * E * u"eV")) # cm
end

"""
    kₑ(E::Float64)

Electron wavenumber (inverse wavelength) in rad⋅cm⁻¹.
"""
kₑ(E::Float64) = 2π / λₑ(E)
# = 0.01*sqrt(2.0*9.109e-31*1.602e-19*E)/1.05457e-34


text:
using Oceananigans
using Oceananigans.Units
using ParameterEstimocean
using LinearAlgebra
using GLMakie
using DataDeps
using Distributions

using Oceananigans.TurbulenceClosures: RiBasedVerticalDiffusivity

using Oceananigans.TurbulenceClosures:
    PiecewiseLinearRiDependentTapering,
    ExponentialRiDependentTapering,
    HyperbolicTangentRiDependentTapering

include("utils.jl")

#####
##### Setup
#####

cases = ["free_convection",
         "strong_wind_weak_cooling",
         "weak_wind_strong_cooling",
         "strong_wind",
         "strong_wind_no_rotation"]

paths = [@datadep_str("two_day_suite_1m/$(case)_instantaneous_statistics.jld2") for case in cases]
times = [2hours, 12hours, 24hours] # 44hours
field_names = (:b, :u, :v)
transformation = ZScore()

z = [-256.0,
     -238.3,
     -207.2,
     -182.3,
     -162.5,
     -146.5,
     -133.0,
     -121.3,
     -110.5,
     -100.2,
     - 90.0,
     - 80.0,
     - 70.0,
     - 60.0,
     - 50.0,
     - 40.0,
     - 30.0,
     - 20.0,
     - 10.0,
        0.0]

regrid = RectilinearGrid(size=length(z)-1; z, topology=(Flat, Flat, Bounded))

observations = [SyntheticObservations(paths[i]; field_names, times, transformation, regrid) for i=1:5]

colorcycle = [:black, :red, :darkblue, :orange, :pink1, :seagreen, :magenta2]
markercycle = [:rect, :utriangle, :star5, :circle, :cross, :+, :pentagon]

function make_figure_axes(n=1)
    fig = Figure(resolution=(1200, n*400))
    axs = []
    for i = 1:n
        ax_b = Axis(fig[i, 1], xlabel = "Buoyancy \n[cm s⁻²]", ylabel = "z [m]")
        ax_u = Axis(fig[i, 2], xlabel = "x-velocity, u \n[cm s⁻¹]")
        ax_v = Axis(fig[i, 3], xlabel = "y-velocity, v \n[cm s⁻¹]")
        push!(axs, (ax_b, ax_u, ax_v))
    end
    n == 1 && (axs = first(axs))
    return fig, axs
end

function plot_fields!(axs, b, u, v, label, color, grid=first(observations).grid; kwargs...)
    z = znodes(Center, grid)
    ## Note unit conversions below, eg m s⁻² -> cm s⁻²:
    lines!(axs[1], 1e2 * b, z; color, label, kwargs...)
    lines!(axs[2], 1e2 * u, z; color, label, kwargs...)
    lines!(axs[3], 1e2 * v, z; color, label, kwargs...)
    return nothing
end

#####
##### Calibration
#####

tapering = ExponentialRiDependentTapering()
ri_based_closure = RiBasedVerticalDiffusivity(Ri_dependent_tapering=tapering)

simulation = ensemble_column_model_simulation(observations;
                                              Nensemble = 2000,
                                              architecture = CPU(),
                                              tracers = (:b, :e),
                                              closure = ri_based_closure)

Qᵘ = simulation.model.velocities.u.boundary_conditions.top.condition
Qᵇ = simulation.model.tracers.b.boundary_conditions.top.condition
N² = simulation.model.tracers.b.boundary_conditions.bottom.condition

simulation.Δt = 10minutes

for (i, obs) in enumerate(observations)
    view(Qᵘ, :, i) .= obs.metadata.parameters.momentum_flux
    view(Qᵇ, :, i) .= obs.metadata.parameters.buoyancy_flux
    view(N², :, i) .= obs.metadata.parameters.N²_deep
end

priors = (ν₀   = ScaledLogitNormal(bounds=(0, 1)),
          κ₀   = ScaledLogitNormal(bounds=(0, 1)),
          Ri₀ν = Normal(-0.5, 1.0),
          Ri₀κ = Normal(-0.5, 1.0),
          Riᵟν = lognormal(mean=1.0,  std=1.0),
          Riᵟκ = lognormal(mean=1.0,  std=1.0))

free_parameters = FreeParameters(priors)
calibration = InverseProblem(observations, simulation, free_parameters)
eki = EnsembleKalmanInversion(calibration; convergence_rate=0.95)

#####
##### Visualization
#####

modeled_time_serieses = calibration.time_series_collector.field_time_serieses
observed, mean_modeled, best_modeled, worst_modeled = [], [], [], []
Nt = length(first(observations).times)
for (c, obs) in enumerate(observations)
    push!(observed,      map(name -> interior(obs.field_time_serieses[name][Nt], 1, 1, :), field_names))
    push!(mean_modeled,  map(name -> interior(  modeled_time_serieses[name][Nt], 1, c, :), field_names))
    push!(best_modeled,  map(name -> interior(  modeled_time_serieses[name][Nt], 2, c, :), field_names))
    push!(worst_modeled, map(name -> interior(  modeled_time_serieses[name][Nt], 3, c, :), field_names))
end

function compare_model_observations(model_label="modeled")
    fig, axs = make_figure_axes(length(observations))
    mean_model_label = "mean " * model_label
    best_model_label = "best " * model_label
    worst_model_label = "worst " * model_label
    for (c, obs) in enumerate(observations)
        plot_fields!(axs[c], observed[c]..., "observed at t = " * prettytime(times[end]), (:black, 0.6), linewidth=5)
        plot_fields!(axs[c], mean_modeled[c]..., mean_model_label, :darkblue, linewidth=2)
        plot_fields!(axs[c], best_modeled[c]..., best_model_label, :orange, linewidth=3)
        plot_fields!(axs[c], worst_modeled[c]..., worst_model_label, :seagreen, linewidth=1)
        [axislegend(ax, position=:rb, merge=true, fontsize=10) for ax in axs[c]]
    end
    return fig
end

#####
##### Calibration
#####

while true
    #try
        iterate!(eki, iterations=4)
        latest_summary = eki.iteration_summaries[end]
        @show latest_summary

        # Plot best current parameters
        Niter = eki.iteration

        min_error, k_min = finitefindmin(latest_summary.mean_square_errors)
        max_error, k_max = finitefindmax(latest_summary.mean_square_errors)

        mean_parameters = latest_summary.ensemble_mean
        best_parameters = latest_summary.parameters[k_min]
        worst_parameters = latest_summary.parameters[k_max]

        forward_run!(calibration, [mean_parameters, best_parameters, worst_parameters])
        fig = compare_model_observations("($Niter iters)")
        save("multi_case_model_observation_comparison_iteration_$Niter.png", fig)
        display(fig)

    # catch err
    #     err isa InterruptException && break
    #     throw(err)
    # end
end

#=
ensemble_means = NamedTuple(n => map(summary -> summary.ensemble_mean[n], eki.iteration_summaries)
                            for n in calibration.free_parameters.names)

fig = Figure()
ax = Axis(fig[1, 1], xlabel = "Ensemble Kalman iteration", ylabel = "Parameter value")

for (i, name) in enumerate(calibration.free_parameters.names)
    label = string(name)
    marker = markercycle[i]
    color = colorcycle[i]
    scatterlines!(ax, 0:Niter, parent(ensemble_means[name]); marker, color, label)
end
=#

max_stars_repo_name: glwagner/LocalOceanClosureCalibration

max_issues_repo_name: glwagner/LocalOceanClosureCalibration

max_forks_repo_name: glwagner/LocalOceanClosureCalibration


text:
module TurtleEx1
using Serd, Serd.CSerd, Serd.RDF

# https://www.w3.org/TeamSubmission/turtle/#sec-examples
const turtle_path = joinpath(@__DIR__, "turtle_ex1.ttl")
const turtle_alt_path = joinpath(@__DIR__, "turtle_ex1_alt.ttl")
const turtle = read(turtle_path, String)
const turtle_alt = read(turtle_alt_path, String) # Round-tripped through Serd

const statements = Statement[
  Prefix("rdf"),
  Prefix("dc", "http://purl.org/dc/elements/1.1/"),
  Prefix("ex", "http://example.org/stuff/1.0/"),
  Triple(
    Resource("http://www.w3.org/TR/rdf-syntax-grammar"),
    Resource("dc", "title"),
    Literal("RDF/XML Syntax Specification (Revised)")),
  Triple(
    Resource("http://www.w3.org/TR/rdf-syntax-grammar"),
    Resource("ex", "editor"),
    Blank("b1")),
  Triple(
    Blank("b1"),
    Resource("ex", "fullname"),
    Literal("Dave Beckett")),
  Triple(
    Blank("b1"),
    Resource("ex", "homePage"),
    Resource("http://purl.org/net/dajobe/")),
]

const serd_prefixes = [
  (SerdNode("rdf", SERD_LITERAL), SerdNode("http://www.w3.org/1999/02/22-rdf-syntax-ns#", SERD_URI)),
  (SerdNode("dc", SERD_LITERAL), SerdNode("http://purl.org/dc/elements/1.1/", SERD_URI)),
  (SerdNode("ex", SERD_LITERAL), SerdNode("http://example.org/stuff/1.0/", SERD_URI))
]

const serd_triples = [
  SerdStatement(
    0,
    nothing,
    SerdNode("http://www.w3.org/TR/rdf-syntax-grammar", SERD_URI),
    SerdNode("dc:title", SERD_CURIE),
    SerdNode("RDF/XML Syntax Specification (Revised)", SERD_LITERAL),
    nothing,
    nothing,
  ),
  SerdStatement(
    SERD_ANON_O_BEGIN,
    nothing,
    SerdNode("http://www.w3.org/TR/rdf-syntax-grammar", SERD_URI),
    SerdNode("ex:editor", SERD_CURIE),
    SerdNode("b1", SERD_BLANK),
    nothing,
    nothing,
  ),
  SerdStatement(
    SERD_ANON_CONT,
    nothing,
    SerdNode("b1", SERD_BLANK),
    SerdNode("ex:fullname", SERD_CURIE),
    SerdNode("Dave Beckett", SERD_LITERAL),
    nothing,
    nothing,
  ),
  SerdStatement(
    SERD_ANON_CONT,
    nothing,
    SerdNode("b1", SERD_BLANK),
    SerdNode("ex:homePage", SERD_CURIE),
    SerdNode("http://purl.org/net/dajobe/", SERD_URI),
    nothing,
    nothing,
  ),
]

const serd_quads = [
  SerdStatement(
    0,
    SerdNode("ex:graph", SERD_CURIE),
    SerdNode("http://www.w3.org/TR/rdf-syntax-grammar", SERD_URI),
    SerdNode("dc:title", SERD_CURIE),
    SerdNode("RDF/XML Syntax Specification (Revised)", SERD_LITERAL),
    nothing,
    nothing,
  ),
  SerdStatement(
    SERD_ANON_O_BEGIN,
    SerdNode("ex:graph", SERD_CURIE),
    SerdNode("http://www.w3.org/TR/rdf-syntax-grammar", SERD_URI),
    SerdNode("ex:editor", SERD_CURIE),
    SerdNode("b1", SERD_BLANK),
    nothing,
    nothing,
  ),
  SerdStatement(
    SERD_ANON_CONT,
    SerdNode("ex:graph", SERD_CURIE),
    SerdNode("b1", SERD_BLANK),
    SerdNode("ex:fullname", SERD_CURIE),
    SerdNode("Dave Beckett", SERD_LITERAL),
    nothing,
    nothing,
  ),
  SerdStatement(
    SERD_ANON_CONT,
    SerdNode("ex:graph", SERD_CURIE),
    SerdNode("b1", SERD_BLANK),
    SerdNode("ex:homePage", SERD_CURIE),
    SerdNode("http://purl.org/net/dajobe/", SERD_URI),
    nothing,
    nothing,
  ),
]

end


text:
using ITensors,
      Test,
      Suppressor

@testset "Example Codes" begin

  examples_dir = joinpath("..", "examples")

  @testset "Basic Ops" begin
    @test_nowarn begin
      @capture_out begin
        include(joinpath(examples_dir, "basic_ops", "basic_ops.jl"))
      end
    end
  end

  @testset "DMRG with Observer" begin
    @test_nowarn begin
      @capture_out begin
        include(joinpath(examples_dir, "dmrg", "1d_ising_with_observer.jl"))
      end
    end
  end

  @testset "Package Compile Code" begin
    @test_nowarn begin
      @capture_out begin
        include(joinpath("..", "src", "packagecompile", "precompile_itensors.jl"))
      end
    end
  end

end

nothing


text:
function manual()
    firstChain::Union{Chain,Nothing} = comp.first_child_
    resIter::Union{Residue,Nothing} = firstChain.first_child_
    counter::Int = 0
    while !isnothing(firstChain)
        while !isnothing(resIter)
            if resIter.name_ == "ALA"
                counter += 1
            end
        resIter = resIter.next_
        end
    firstChain = firstChain.next_
   !isnothing(firstChain) ? resIter = firstChain.first_child_ : nothing
    end
    return counter
end

function man2()
counter::Int = 0
for x in ResidueIterator(comp)
    if x.name_ == "ALA"
        counter +=1
    end
end
return counter
end

function man3()
    counter::Int = 0
    for x in eachrow(dfs.residues)
        if x.name_ =="ALA"
            counter += 1
        end
    end
return counter
end

function man4()
    counter::Int = 0
    for x in dfs.residues.name_
        if x == "ALA"
            counter += 1
        end
    end
return counter
end

function man5()
return length(filter(x->x == "ALA", dfs.residues.name_))
end


text:
using QCompEmulator
using Test

@testset "QCompEmulator.jl" begin
    # Write your tests here.
   @test my_f(2,1) == 5
end


text:
### Numbering
@inbounds function renumber(data, original, new, numbering)
    if numbering
        lookup = Dict{Int, Int}()
        for (x, y) in zip(original, new)
            lookup[x] = y
        end

        for i = 1:length(data)
            data[i] = lookup[data[i]]
        end
    end

    return data
end

### Save data in h5-file or xlsx-file
@inbounds function savedata(args...; group, header, path, info = "", skip = [])
    extension = match(r"\.[A-Za-z0-9]+$", path).match

    if extension == ".h5"
        h5open(path, "w") do file
            for i in keys(group)
                if group[i] != 0
                    write(file, String(i), getfield(args[group[i]], i))
                    try
                        head = header[i]
                        N = size(getfield(args[group[i]], i), 2)
                        if N == 1
                            N = size(getfield(args[group[i]], i), 1)
                        end
                        atr = Dict(string("row", k) => string(head[1, k], ": ", head[2, k]) for k = 1:N)
                        h5writeattr(path, String(i), atr)
                    catch
                    end
                end
            end
            if !isempty(info)
                write(file, "info", info)
            end
        end
        println("Data is successfully exported to $path")
    elseif extension == ".xlsx"
        k = 1
        XLSX.openxlsx(path, mode = "w") do xf
            for i in keys(group)
                if group[i] != 0
                    if k == 1
                        sheet = xf[k]
                        XLSX.rename!(sheet, String(i))
                    else
                        XLSX.addsheet!(xf, String(i))
                        sheet = xf[k]
                    end
                    N = size(getfield(args[group[i]], i), 2)
                    if N == 1
                        N = size(getfield(args[group[i]], i), 1)
                    end
                    sheet["A1"] = header[i][:, 1:N]
                    sheet["A3"] = getfield(args[group[i]], i)
                    k += 1
                end
            end
            if !isempty(info)
                XLSX.addsheet!(xf, "info")
                sheet = xf[k]
                sheet["A1"] = info
            end
        end
        println("Data is successfully exported to $path")
    else
        throw(ErrorException("the SAVE data format is not supported"))
    end

end

### Start row in xlsx-file
function startxlsx(xf)
    start = 1
    for r in XLSX.eachrow(xf)
        if !isa(r[1], String)
            start = XLSX.row_number(r)
            break
        end
    end
    return start
end

### Weighted least-squares method
@inbounds function wls(G, b, method)
    if method == "lu"
        F = lu(G)
        L, U, p, q, Rs = F.:(:)
        x = U \  (L \ ((Rs .* b)[p]))
        x = x[sortperm(q)]
    else
        x = G \ b
    end

    return x
end

### Least-squares method
@inbounds function ls(A, b, method)
    if method == "lu"
        F = lu(A)
        L, U, p, q, Rs = F.:(:)
        x = U \  (L \ ((Rs .* b)[p]))
        x = x[sortperm(q)]
    else
        x = A \ b
    end

    return x
end

### Least absolute value method
@inbounds function lav(J, z, Nvar, Nmeasure, settings)
    t = [zeros(2 * Nvar); ones(2 * Nmeasure)]
    E = spdiagm(0 => ones(Nmeasure))
    A = [J -J E -E]
    lb = zeros(2 * (Nvar + Nmeasure))

    if settings.lav[:optimize] == 2
        model = Model(Ipopt.Optimizer)
    else
        model = Model(GLPK.Optimizer)
    end

    @variable(model, x[i = 1:(2 * (Nvar + Nmeasure))])
    @constraint(model, lb .<= x)
    @constraint(model, A * x .== z)
    @objective(model, Min, t' * x)
    optimize!(model)
    println("Least Aabsolute value state estimation status: $(JuMP.termination_status(model))")

    return JuMP.value.(x)
end

### Try-catch next element
function nextelement(set, current)
    value = 0.0
    try
        value = set[current + 1]
        if value == "all"
            value = -200.0
        end
        value = Float64(value)
    catch
        throw(ErrorException("the name-value pair setting is missing"))
    end

    return value
end

### Data structure
function datastruct(data, max; var = "")
    Nrow, Ncol = size(data)
    if Ncol < max
        throw(DomainError(var, "dimension mismatch"))
    end

    return Nrow
end

### The sparse inverse subset of a real sparse square matrix: Construct the factorization
function constructfact(A)
    F = lu(A)
    L, U, p, q, Rss = F.:(:)
    U = copy(transpose(U))
    n = size(U, 2)
    d = fill(0.0, n)
    Rs = fill(0.0, n)
    for i = 1:n
        Rs[i] = Rss[p[i]]
        d[i] = U[i,i]
        U[i, i] = 0.0
        L[i, i] = 0.0
        for j = U.colptr[i]:(U.colptr[i + 1] - 1)
            if i != U.rowval[j]
                U.nzval[j] = U.nzval[j] / d[i]
            end
        end
    end
    dropzeros!(U)
    dropzeros!(L)
    S = A[p, q]
    S = S + transpose(S)

    return S, L, U, d, p, q, Rs
end

### The sparse inverse subset of a real sparse square matrix: Compute the elimination tree of a sparse matrix
# Copyright (c) 2013-2014 Viral Shah, Douglas Bates and other contributors
# https://github.com/JuliaPackageMirrors/SuiteSparse.jl/blob/master/src/csparse.jl
# Based on Direct Methods for Sparse Linear Systems, T. A. Davis, SIAM, Philadelphia, Sept. 2006.
function etree(A)
    n = size(A, 2)
    parent = fill(0, n)
    ancestor = fill(0, n)
    for k in 1:n, p in A.colptr[k]:(A.colptr[k + 1] - 1)
        i = A.rowval[p]
        while i != 0 && i < k
            inext = ancestor[i]
            ancestor[i] = k
            if (inext == 0)
                parent[i] = k
            end
            i = inext
        end
    end

    head = fill(0, n)
    next = fill(0, n)
    for j in n:-1:1
        if (parent[j] == 0); continue; end
        next[j] = head[parent[j]]
        head[parent[j]] = j
    end
    stack = Int64[]
    for j in 1:n
        if (parent[j] != 0)
            continue
        end
        push!(stack, j)
        while (!isempty(stack))
            p = stack[end]
            i = head[p]
            if (i == 0)
                pop!(stack)
            else
            head[p] = next[i]
            push!(stack, i)
            end
        end
    end

    return parent
end

### The sparse inverse subset of a real sparse square matrix: Find nonzero pattern of Cholesky
# Copyright (c) 2013-2014 Viral Shah, Douglas Bates and other contributors
# https://github.com/JuliaPackageMirrors/SuiteSparse.jl/blob/master/src/csparse.jl
# Based on Direct Methods for Sparse Linear Systems, T. A. Davis, SIAM, Philadelphia, Sept. 2006.
function symbfact(A, parent)
    m,n = size(A); Ap = A.colptr; Ai = A.rowval
    col = Int64[]; sizehint!(col, n)
    row = Int64[]; sizehint!(row, n)

    visited = falses(n)
    for k = 1:m
        visited = falses(n)
        visited[k] = true
        for p in Ap[k]:(Ap[k + 1] - 1)
            i = Ai[p]
            if i > k continue end
            while !visited[i]
                push!(col, i)
                push!(row, k)
                visited[i] = true
                i = parent[i]
            end
        end
    end
    R = sparse([col; collect(1:n)], [row; collect(1:n)], ones(length(row) + n), n, n)

    return R
end

### The sparse inverse subset of a real sparse square matrix
# Copyright 2011, Timothy A. Davis
# http://www.suitesparse.com
@inbounds function sparseinv(L, U, d, p, q, Rs, R)
    Zpattern = R + R'
    n = size(Zpattern, 1)

    Zdiagp = fill(1, n)
    Lmunch = fill(0, n)

    znz = nnz(Zpattern)
    Zx = fill(0.0, znz)
    z = zeros(n)
    k = 0

    Zcolptr = Zpattern.colptr
    Zrowval = Zpattern.rowval
    flag = true
    for j = 1:n
        pdiag = -1
        for p = Zcolptr[j]:(Zcolptr[j + 1] - 1)
            if pdiag == -1
                if Zrowval[p] == j
                    pdiag = p
                    Zx[p] = 1 / (d[j] / Rs[j])
                end
            end
        end
        if pdiag == -1
            flag = false
          break
        end
        Zdiagp[j] = pdiag
    end

    if flag
        for k = 1:n
            Lmunch[k] = L.colptr[k + 1] - 1
        end

        for j = n:-1:1
            for p = Zdiagp[j]:Zcolptr[j + 1] - 1
                z[Zrowval[p]] = Zx[p]
            end

            for p = Zdiagp[j]-1:-1:Zcolptr[j]
                k = Zrowval[p]
                zkj = 0.0
                for up = U.colptr[k]:U.colptr[k + 1] - 1
                    i = U.rowval[up]
                    if i > k && z[i] != 0.0
                        zkj -= U.nzval[up] * z[i]
                    end
                end
                z[k] = zkj
            end

            for p = Zdiagp[j]-1:-1:Zcolptr[j]
                k = Zrowval[p]
                if (Lmunch[k] < L.colptr[k] || L.rowval[Lmunch[k]] != j)
                    continue
                end
                ljk = L.nzval[Lmunch[k]] * Rs[k] / Rs[L.rowval[Lmunch[k]]]
                Lmunch[k] -= 1
                for zp = Zdiagp[k]:Zcolptr[k + 1] - 1
                    Zx[zp] -= z[Zrowval[zp]] * ljk
                end
            end

            for p = Zcolptr[j]:Zcolptr[j + 1] - 1
                i = Zrowval[p]
                Zx[p] = z[i]
                z[i] = 0.0
            end
        end
    end

    idx = findall(!iszero, Zpattern)
    for (k, i) in enumerate(idx)
        Zpattern[i] = Zx[k]
    end

    Zpattern = Zpattern[invperm(q), invperm(p)]

    return Zpattern
end

### Independent system of equations
function independentset(Jt)
    F = qr(Jt)
    R = F.R
    P = F.pcol

    flag = (10^-5) * R[1,1]
    r = 0
    @inbounds for i = 1:size(R,1)
        if abs(R[i,i]) >= flag
            r = i
        end
    end
    p = P[1:r]
    Ht = Jt[:, p]
    dropzeros!(Ht)

    return Ht
end

### Make factor graph
function factorgraph(Ht, Nvar, Nfac, Nlink, Nind)
    ########## Define indirect factor nodes arrays ##########
    row = fill(0, Nlink); col = similar(row)
    rowptr = fill(0, Nind); colptr = fill(0, Nvar)
    idxi = 1; idxr = 1
    idxT = findall(!iszero, Ht); links = Dict(i => Int64[] for i = 1:Nfac)
    @inbounds for i in idxT
        if (Ht.colptr[i[2] + 1] - Ht.colptr[i[2]]) != 1
            row[idxi] = idxr
            col[idxi] = i[1]
            push!(links[idxr], i[1])

            idxi += 1

            colptr[i[1]] += 1
            rowptr[idxr] = idxi

            if idxT[Ht.colptr[i[2] + 1] - 1] == i
                idxr += 1
            end
        end
    end
    pushfirst!(colptr, 1)
    pushfirst!(rowptr, 1)
    colptr = cumsum(colptr)

    ########## Initialize arrays for messages ##########
    Vvar_fac = fill(1e-5, Nlink)
    Wfac_var = similar(Vvar_fac)

    ########## Indices for sending messages ##########
    temp = collect(1:Nlink)
    sort_col_idx = sortperm(col)
    to_fac = temp[sort_col_idx]

    new_row = row[sort_col_idx]
    sort_row_idx = sortperm(new_row)
    to_var = temp[sort_row_idx]

    return Vvar_fac, Wfac_var, to_fac, to_var, colptr, rowptr, idxT, row, col, links
end

### Matpower parse lines
@inline function parse_line(i, line, flag, str)
    sublines = split(line, "[")[end]
    sublines = split(sublines, ";")

    for k = 1:length(sublines)
        subline = strip(sublines[k])
        if !isempty(subline)
            push!(str, subline)
        end
    end

    if occursin("]", line)
        lastlines = split(line, "]")[1]
        lastlines = split(lastlines, ";")
        for k = 1:length(lastlines)
            last = lastlines[k]
            if isempty(strip(last))
                pop!(str)
            end
        end
        flag = false
    end

    return flag, str
end

max_stars_repo_name: mcosovic/JuliaGrid

max_issues_repo_name: mcosovic/JuliaGrid

max_forks_repo_name: mcosovic/JuliaGrid


text:
using BenchmarkTools
using GeoPhyInv
using Test
#md using Plots
#md gr()

# # Intro

#md # ```@docs
#md # GeoPhyInv.Medium
#md # Medium(::Symbol, ::Real)
#md # ```

# # Examples

# Load a predefined model.
medium = Medium(:elastic_homo3D);

# Get the medium parameters that are stored.
names(medium)

# Look at the reference values.
medium.ref

# Inspect the bounds.
medium.bounds

# Plotting is easy
#md p1=heatmap(medium, :vp)
#md plot(p1,size=(800,400))

# To construct an instance of 2-D `Medium`, we need ranges for `z` and `x`.
mgrid = [range(0.0, stop = 10.0, step = 0.1), range(0.0, stop = 30.0, step = 0.2)];

# For 3-D, it is `z`, `y` and `x`.
mgrid = fill(range(-10, stop = 10.0, step = 0.1), 3)

# Allocate basic medium parameters on the grid.
medium = Medium(mgrid, [:vp, :rho, :vs]);

# Bounds for these parameters should be input for modeling or inversion. Use `update!`
vpb = [2100.0, 2200.0];
vsb = [1500, 1700];
rhob = [2100.0, 2300.0];
update!(medium, [:vp, :vs, :rho], [vpb, vsb, rhob]);

# Just fill `medium` with average (reference) values.
fill!(medium);

# Once the basic medium parameters are input, we can access some other derived parameters.
medium[:Zp];

# Otherwise, we can manually update parameters of `medium`.
medium[:vp] .= 3000.0;
medium[:vs] .= 2000.0;

println(medium)

# A model can be also be updated by adding random noise.
update!(medium, [:vp, :rho], randn_perc = 1.0);

# # Methods
#md # ```@docs
#md md # Base.getindex(::Medium, ::Symbol)
#md # GeoPhyInv.update!(::GeoPhyInv.Medium, ::Vector{Symbol},)
#md # Base.copyto!(x::AbstractArray, medium::Medium, fields::Vector{Symbol})
#md # Base.vec(medium::Medium, ::Symbol)
#md # ```


text:
using SparseArrays, LinearAlgebra

"Check whether a package is installed."
function pkg_installed(name::String, uuid::String)
	pkg_id = Base.PkgId(Base.UUID(uuid), name::String)
	return Base.locate_package(pkg_id) != nothing
end

# generate a random pos def matrix with eigenvalues between 0.1 and 2
generate_pos_def_matrix(rng::MersenneTwister, n::Integer, aMin::Float64 = 0.1, aMax::Float64 = 2.) = generate_pos_def_matrix(rng, n, aMin, aMax, MT = Float64)
function generate_pos_def_matrix(rng::MersenneTwister, n::Integer, aMin::Real = 0.1, aMax::Real = 2.; MT::Type{<:AbstractFloat} = Float64 )
	X = rand(rng, MT, n, n)
	# any real square matrix can be QP decomposed into a orthogonal matrix and an uppertriangular matrix R
	Q, R = qr(X)
	eigs = rand(rng ,MT, n) .* (MT(aMax) .- MT(aMin)) .+ MT(aMin)
	X = Q * Matrix(Diagonal(eigs)) * Q'
	X = MT(0.5) * (X + X')
	return X
end

function is_numerically_pos_sem_def(X, atol)
	X = X ./ 2
	X = X + X'

	F = eigfact(X)
	if size(find( x-> x < -atol, F[:values]), 1) == 0
		return true
	else
		return false
	end
end

function is_numerically_symmetric(X,atol)
	n = size(X, 2)
	for i = 1:n-1, j = i+1:n
		if abs(X[i, j] - X[j, i]) >= atol
			return false
		end
	end
	return true
end

function find_nonsymmetric_component(X)
	for i = 2:size(X, 1), j = 1:(i - 1)
		if abs(X[i, j] - X[j, i]) > 0.0
			return i, j, abs(X[i, j] - X[j, i])
		end
	end
end

function find_different_elements(A, B)
	if size(A) != size(B)
		error("Matrices are not the same size")
	end
	m, n = size(A)
	diff_el = Array[]
	for iii = 1:m, jjj = 1:n
		if A[iii, jjj] != B[iii, jjj]
			push!(diff_el, [iii, jjj])
		end
	end
	return diff_el
end

function duplicate_sparsity_pattern(A)
	m, n = size(A)
	B = zeros(m, n)
	for iii = 1:m, jjj = 1:n
		if A[iii, jjj] != 0
			B[iii, jjj] = 1
		end
	end
	return B
end

function apply_pattern!(A::AbstractMatrix, pattern::AbstractMatrix)
  @assert size(A) == size(pattern) "Matrix A and pattern must have same dimensions."
  m, n = size(A)
  for i = 1:m, j = 1:n
    if pattern[i, j] == 0
      A[i, j] = 0.
    end
  end
end

# create a feasible SDP with one PSDConeTriangle constraint
# min c' x
# s.t. At x + s == bt
#       s ∈ PSDConeTriangle
# choose At and bt in such a way that S has the provided sparsity pattern
function feasible_sdp_with_pattern(rng::MersenneTwister, pattern::AbstractMatrix)
  n = size(pattern, 1)
  d = div(n * (n + 1), 2)

  S = generate_pos_def_matrix(rng, n, 0.1, 2.)
  apply_pattern!(S, pattern)
  S = Symmetric(S, :U)

  A1 = rand(rng, n, n)
  apply_pattern!(A1, pattern)
  A1 = Symmetric(A1, :U)
  A = hcat(A1[:])
  s = S[:]
  x = rand(rng, 1)
  b = A * x + s
  B = reshape(b, n, n)

  Y = generate_pos_def_matrix(rng, n,  0.1, 1.)
  y = vec(Y)
  P = sparse(zeros(1, 1))
  q = -P * x - A' * y

  Ct = [COSMO.PsdConeTriangle(d)];
  At = zeros(d)
  bt = zeros(d)
  COSMO.extract_upper_triangle!(A1, At, sqrt(2))
  COSMO.extract_upper_triangle!(B, bt, sqrt(2))
  At = hcat(At)
  return P, q, At, bt, Ct
end

"Take a vector `svec` representing the `d` upper-triangular entries of a matrix `X`, and return `X`."
function matrixify(svec::AbstractVector)
  n = - 0.5 + sqrt(0.25 + 2 * length(svec))
  n = Int64(n)
  X = zeros(n, n)
  COSMO.populate_upper_triangle!(X, svec, 1 /sqrt(2))
  return Symmetric(X, :U)
end

function recreate_sparse_matrix(A::SparseMatrixCSC)
	rowInd = A.rowval
	colPtr = A.colptr
	val = A.nzval

	#compute column indices
	colInd = zeros(Int64, length(rowInd))
	cval = 1
	for iii = 2:length(colPtr)
		currentPtr = colPtr[iii]
		prevPtr = colPtr[iii - 1]
		colInd[prevPtr:currentPtr - 1] = cval
		cval += 1
	end

	# sort rowInd and vals
	p = sortperm(rowInd)
	return sparse(rowInd, colInd, val, size(A, 1), size(A, 2))
end

# Geometric mean from https://github.com/JuliaStats/StatsBase.jl
function gmean(a::AbstractArray{T}) where T<:Real
	s = 0.0
	n = length(a)
	for i in 1:n
		tmp = a[i]
		if tmp < 0.0
			throw(DomainError())
		elseif tmp == 0.0
			return 0.0
		else
			s += log(tmp)
		end
	end
	return exp(s / n)
end

function check_clique_tree(snd, sep, snd_par, snd_child)
	# c.t. should have one root
	ind = findall( x -> x == 0, snd_par)
	length(ind) == 1 || return false

	for (ch, p) in enumerate(snd_par)
		# removed
		if p == -1
			isempty(snd[ch]) || return false
			isempty(sep[ch]) || return false
		# root
		elseif p == 0
			isempty(sep[ch]) || return false
			# find all children of ch
			children = findall(x -> x == ch, snd_par)
			mapreduce(x -> x ∈ snd_child[ch], *,  children) || return false
		else
			intersect(union(snd[ch], sep[ch]), union(snd[p], sep[p])) == sep[ch] || return false
			children = findall(x -> x == ch, snd_par)
			if !isempty(children)
				mapreduce(x -> x ∈ snd_child[ch], *,  children) || return false
			end
		end
	end
	return true
end


text:
__precompile__()

module Tries

    import Base: start, next, done, getindex, setindex!, get,
                 haskey, keys

    using Compat: uninitialized, Nothing, Cvoid, AbstractDict

    export Trie, subtrie, keys_with_prefix, path

    include("trie.jl")
end