text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## Diagonal matrices

struct Diagonal{T,V<:AbstractVector{T}} <: AbstractMatrix{T}
    diag::V

    function Diagonal{T,V}(diag) where {T,V<:AbstractVector{T}}
        require_one_based_indexing(diag)
        new{T,V}(diag)
    end
end
Diagonal(v::AbstractVector{T}) where {T} = Diagonal{T,typeof(v)}(v)
Diagonal{T}(v::AbstractVector) where {T} = Diagonal(convert(AbstractVector{T}, v)::AbstractVector{T})

"""
    Diagonal(A::AbstractMatrix)

Construct a matrix from the diagonal of `A`.

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

julia> Diagonal(A)
3×3 Diagonal{Int64,Array{Int64,1}}:
 1  ⋅  ⋅
 ⋅  5  ⋅
 ⋅  ⋅  9
```
"""
Diagonal(A::AbstractMatrix) = Diagonal(diag(A))

"""
    Diagonal(V::AbstractVector)

Construct a matrix with `V` as its diagonal.

# Examples
```jldoctest
julia> V = [1, 2]
2-element Array{Int64,1}:
 1
 2

julia> Diagonal(V)
2×2 Diagonal{Int64,Array{Int64,1}}:
 1  ⋅
 ⋅  2
```
"""
Diagonal(V::AbstractVector)

Diagonal(D::Diagonal) = D
Diagonal{T}(D::Diagonal{T}) where {T} = D
Diagonal{T}(D::Diagonal) where {T} = Diagonal{T}(D.diag)

AbstractMatrix{T}(D::Diagonal) where {T} = Diagonal{T}(D)
Matrix(D::Diagonal) = diagm(0 => D.diag)
Array(D::Diagonal) = Matrix(D)

# For D<:Diagonal, similar(D[, neweltype]) should yield a Diagonal matrix.
# On the other hand, similar(D, [neweltype,] shape...) should yield a sparse matrix.
# The first method below effects the former, and the second the latter.
similar(D::Diagonal, ::Type{T}) where {T} = Diagonal(similar(D.diag, T))
# The method below is moved to SparseArrays for now
# similar(D::Diagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)

copyto!(D1::Diagonal, D2::Diagonal) = (copyto!(D1.diag, D2.diag); D1)

size(D::Diagonal) = (length(D.diag),length(D.diag))

function size(D::Diagonal,d::Integer)
    if d<1
        throw(ArgumentError("dimension must be ≥ 1, got $d"))
    end
    return d<=2 ? length(D.diag) : 1
end

@inline function getindex(D::Diagonal, i::Int, j::Int)
    @boundscheck checkbounds(D, i, j)
    if i == j
        @inbounds r = D.diag[i]
    else
        r = diagzero(D, i, j)
    end
    r
end
diagzero(::Diagonal{T},i,j) where {T} = zero(T)
diagzero(D::Diagonal{Matrix{T}},i,j) where {T} = zeros(T, size(D.diag[i], 1), size(D.diag[j], 2))

function setindex!(D::Diagonal, v, i::Int, j::Int)
    @boundscheck checkbounds(D, i, j)
    if i == j
        @inbounds D.diag[i] = v
    elseif !iszero(v)
        throw(ArgumentError("cannot set off-diagonal entry ($i, $j) to a nonzero value ($v)"))
    end
    return v
end

## structured matrix methods ##
function Base.replace_in_print_matrix(A::Diagonal,i::Integer,j::Integer,s::AbstractString)
    i==j ? s : Base.replace_with_centered_mark(s)
end

parent(D::Diagonal) = D.diag

ishermitian(D::Diagonal{<:Real}) = true
ishermitian(D::Diagonal{<:Number}) = isreal(D.diag)
ishermitian(D::Diagonal) = all(ishermitian, D.diag)
issymmetric(D::Diagonal{<:Number}) = true
issymmetric(D::Diagonal) = all(issymmetric, D.diag)
isposdef(D::Diagonal) = all(isposdef, D.diag)

factorize(D::Diagonal) = D

real(D::Diagonal) = Diagonal(real(D.diag))
imag(D::Diagonal) = Diagonal(imag(D.diag))

iszero(D::Diagonal) = all(iszero, D.diag)
isone(D::Diagonal) = all(isone, D.diag)
isdiag(D::Diagonal) = all(isdiag, D.diag)
isdiag(D::Diagonal{<:Number}) = true
istriu(D::Diagonal) = true
istril(D::Diagonal) = true
function triu!(D::Diagonal,k::Integer=0)
    n = size(D,1)
    if !(-n + 1 <= k <= n + 1)
        throw(ArgumentError(string("the requested diagonal, $k, must be at least ",
            "$(-n + 1) and at most $(n + 1) in an $n-by-$n matrix")))
    elseif k > 0
        fill!(D.diag,0)
    end
    return D
end

function tril!(D::Diagonal,k::Integer=0)
    n = size(D,1)
    if !(-n - 1 <= k <= n - 1)
        throw(ArgumentError(string("the requested diagonal, $k, must be at least ",
            "$(-n - 1) and at most $(n - 1) in an $n-by-$n matrix")))
    elseif k < 0
        fill!(D.diag,0)
    end
    return D
end

(==)(Da::Diagonal, Db::Diagonal) = Da.diag == Db.diag
(-)(A::Diagonal) = Diagonal(-A.diag)
(+)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag + Db.diag)
(-)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag - Db.diag)

(*)(x::Number, D::Diagonal) = Diagonal(x * D.diag)
(*)(D::Diagonal, x::Number) = Diagonal(D.diag * x)
(/)(D::Diagonal, x::Number) = Diagonal(D.diag / x)
(*)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag .* Db.diag)
(*)(D::Diagonal, V::AbstractVector) = D.diag .* V

(*)(A::AbstractTriangular, D::Diagonal) = rmul!(copy(A), D)
(*)(D::Diagonal, B::AbstractTriangular) = lmul!(D, copy(B))

(*)(A::AbstractMatrix, D::Diagonal) =
    rmul!(copyto!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), A), D)
(*)(D::Diagonal, A::AbstractMatrix) =
    lmul!(D, copyto!(similar(A, promote_op(*, eltype(A), eltype(D.diag)), size(A)), A))

function rmul!(A::AbstractMatrix, D::Diagonal)
    require_one_based_indexing(A)
    A .= A .* permutedims(D.diag)
    return A
end

function lmul!(D::Diagonal, B::AbstractMatrix)
    require_one_based_indexing(B)
    B .= D.diag .* B
    return B
end

rmul!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal) = typeof(A)(rmul!(A.data, D))
function rmul!(A::UnitLowerTriangular, D::Diagonal)
    rmul!(A.data, D)
    for i = 1:size(A, 1)
        A.data[i,i] = D.diag[i]
    end
    LowerTriangular(A.data)
end
function rmul!(A::UnitUpperTriangular, D::Diagonal)
    rmul!(A.data, D)
    for i = 1:size(A, 1)
        A.data[i,i] = D.diag[i]
    end
    UpperTriangular(A.data)
end

function lmul!(D::Diagonal, B::UnitLowerTriangular)
    lmul!(D, B.data)
    for i = 1:size(B, 1)
        B.data[i,i] = D.diag[i]
    end
    LowerTriangular(B.data)
end
function lmul!(D::Diagonal, B::UnitUpperTriangular)
    lmul!(D, B.data)
    for i = 1:size(B, 1)
        B.data[i,i] = D.diag[i]
    end
    UpperTriangular(B.data)
end

*(D::Adjoint{<:Any,<:Diagonal}, B::Diagonal) = Diagonal(adjoint.(D.parent.diag) .* B.diag)
*(A::Adjoint{<:Any,<:AbstractTriangular}, D::Diagonal) = rmul!(copy(A), D)
function *(adjA::Adjoint{<:Any,<:AbstractMatrix}, D::Diagonal)
    A = adjA.parent
    Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
    adjoint!(Ac, A)
    rmul!(Ac, D)
end

*(D::Transpose{<:Any,<:Diagonal}, B::Diagonal) = Diagonal(transpose.(D.parent.diag) .* B.diag)
*(A::Transpose{<:Any,<:AbstractTriangular}, D::Diagonal) = rmul!(copy(A), D)
function *(transA::Transpose{<:Any,<:AbstractMatrix}, D::Diagonal)
    A = transA.parent
    At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
    transpose!(At, A)
    rmul!(At, D)
end

*(D::Diagonal, B::Adjoint{<:Any,<:Diagonal}) = Diagonal(D.diag .* adjoint.(B.parent.diag))
*(D::Diagonal, B::Adjoint{<:Any,<:AbstractTriangular}) = lmul!(D, collect(B))
*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = (Q = adjQ.parent; rmul!(Array(D), adjoint(Q)))
function *(D::Diagonal, adjA::Adjoint{<:Any,<:AbstractMatrix})
    A = adjA.parent
    Ac = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
    adjoint!(Ac, A)
    lmul!(D, Ac)
end

*(D::Diagonal, B::Transpose{<:Any,<:Diagonal}) = Diagonal(D.diag .* transpose.(B.parent.diag))
*(D::Diagonal, B::Transpose{<:Any,<:AbstractTriangular}) = lmul!(D, copy(B))
function *(D::Diagonal, transA::Transpose{<:Any,<:AbstractMatrix})
    A = transA.parent
    At = similar(A, promote_op(*, eltype(A), eltype(D.diag)), (size(A, 2), size(A, 1)))
    transpose!(At, A)
    lmul!(D, At)
end

*(D::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:Diagonal}) =
    Diagonal(adjoint.(D.parent.diag) .* adjoint.(B.parent.diag))
*(D::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:Diagonal}) =
    Diagonal(transpose.(D.parent.diag) .* transpose.(B.parent.diag))

rmul!(A::Diagonal, B::Diagonal) = Diagonal(A.diag .*= B.diag)
lmul!(A::Diagonal, B::Diagonal) = Diagonal(B.diag .= A.diag .* B.diag)

function lmul!(adjA::Adjoint{<:Any,<:Diagonal}, B::AbstractMatrix)
    A = adjA.parent
    return lmul!(adjoint(A), B)
end
function lmul!(transA::Transpose{<:Any,<:Diagonal}, B::AbstractMatrix)
    A = transA.parent
    return lmul!(transpose(A), B)
end

function rmul!(A::AbstractMatrix, adjB::Adjoint{<:Any,<:Diagonal})
    B = adjB.parent
    return rmul!(A, adjoint(B))
end
function rmul!(A::AbstractMatrix, transB::Transpose{<:Any,<:Diagonal})
    B = transB.parent
    return rmul!(A, transpose(B))
end

# Get ambiguous method if try to unify AbstractVector/AbstractMatrix here using AbstractVecOrMat
mul!(out::AbstractVector, A::Diagonal, in::AbstractVector) = out .= A.diag .* in
mul!(out::AbstractVector, A::Adjoint{<:Any,<:Diagonal}, in::AbstractVector) = out .= adjoint.(A.parent.diag) .* in
mul!(out::AbstractVector, A::Transpose{<:Any,<:Diagonal}, in::AbstractVector) = out .= transpose.(A.parent.diag) .* in

mul!(out::AbstractMatrix, A::Diagonal, in::StridedMatrix) = out .= A.diag .* in
mul!(out::AbstractMatrix, A::Adjoint{<:Any,<:Diagonal}, in::StridedMatrix) = out .= adjoint.(A.parent.diag) .* in
mul!(out::AbstractMatrix, A::Transpose{<:Any,<:Diagonal}, in::StridedMatrix) = out .= transpose.(A.parent.diag) .* in

# ambiguities with Symmetric/Hermitian
# RealHermSymComplex[Sym]/[Herm] only include Number; invariant to [c]transpose
*(A::Diagonal, transB::Transpose{<:Any,<:RealHermSymComplexSym}) = A * transB.parent
*(transA::Transpose{<:Any,<:RealHermSymComplexSym}, B::Diagonal) = transA.parent * B
*(A::Diagonal, adjB::Adjoint{<:Any,<:RealHermSymComplexHerm}) = A * adjB.parent
*(adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}, B::Diagonal) = adjA.parent * B
*(transA::Transpose{<:Any,<:RealHermSymComplexSym}, transD::Transpose{<:Any,<:Diagonal}) = transA.parent * transD
*(transD::Transpose{<:Any,<:Diagonal}, transA::Transpose{<:Any,<:RealHermSymComplexSym}) = transD * transA.parent
*(adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}, adjD::Adjoint{<:Any,<:Diagonal}) = adjA.parent * adjD
*(adjD::Adjoint{<:Any,<:Diagonal}, adjA::Adjoint{<:Any,<:RealHermSymComplexHerm}) = adjD * adjA.parent
mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:RealHermSym}) = mul!(C, A, B.parent)
mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:RealHermSymComplexHerm}) = mul!(C, A, B.parent)
mul!(C::AbstractMatrix, A::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:RealHermSym}) = mul!(C, A, B.parent)
mul!(C::AbstractMatrix, A::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:RealHermSymComplexSym}) = mul!(C, A, B.parent)
mul!(C::AbstractMatrix, A::Adjoint{<:Any,<:Diagonal}, B::Adjoint{<:Any,<:RealHermSymComplexSym}) = C .= adjoint.(A.parent.diag) .* B
mul!(C::AbstractMatrix, A::Transpose{<:Any,<:Diagonal}, B::Transpose{<:Any,<:RealHermSymComplexHerm}) = C .= transpose.(A.parent.diag) .* B

(/)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag ./ Db.diag)

function ldiv!(D::Diagonal{T}, v::AbstractVector{T}) where {T}
    if length(v) != length(D.diag)
        throw(DimensionMismatch("diagonal matrix is $(length(D.diag)) by $(length(D.diag)) but right hand side has $(length(v)) rows"))
    end
    for i = 1:length(D.diag)
        d = D.diag[i]
        if iszero(d)
            throw(SingularException(i))
        end
        v[i] = d\v[i]
    end
    v
end
function ldiv!(D::Diagonal{T}, V::AbstractMatrix{T}) where {T}
    require_one_based_indexing(V)
    if size(V,1) != length(D.diag)
        throw(DimensionMismatch("diagonal matrix is $(length(D.diag)) by $(length(D.diag)) but right hand side has $(size(V,1)) rows"))
    end
    for i = 1:length(D.diag)
        d = D.diag[i]
        if iszero(d)
            throw(SingularException(i))
        end
        for j = 1:size(V,2)
            @inbounds V[i,j] = d\V[i,j]
        end
    end
    V
end

ldiv!(adjD::Adjoint{<:Any,<:Diagonal{T}}, B::AbstractVecOrMat{T}) where {T} =
    (D = adjD.parent; ldiv!(conj(D), B))
ldiv!(transD::Transpose{<:Any,<:Diagonal{T}}, B::AbstractVecOrMat{T}) where {T} =
    (D = transD.parent; ldiv!(D, B))

function ldiv!(D::Diagonal, A::Union{LowerTriangular,UpperTriangular})
    broadcast!(\, parent(A), D.diag, parent(A))
    A
end

function rdiv!(A::AbstractMatrix{T}, D::Diagonal{T}) where {T}
    require_one_based_indexing(A)
    dd = D.diag
    m, n = size(A)
    if (k = length(dd)) ≠ n
        throw(DimensionMismatch("left hand side has $n columns but D is $k by $k"))
    end
    @inbounds for j in 1:n
        ddj = dd[j]
        if iszero(ddj)
            throw(SingularException(j))
        end
        for i in 1:m
            A[i, j] /= ddj
        end
    end
    A
end

function rdiv!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal)
    broadcast!(/, parent(A), parent(A), permutedims(D.diag))
    A
end

rdiv!(A::AbstractMatrix{T}, adjD::Adjoint{<:Any,<:Diagonal{T}}) where {T} =
    (D = adjD.parent; rdiv!(A, conj(D)))
rdiv!(A::AbstractMatrix{T}, transD::Transpose{<:Any,<:Diagonal{T}}) where {T} =
    (D = transD.parent; rdiv!(A, D))

(/)(A::Union{StridedMatrix, AbstractTriangular}, D::Diagonal) =
    rdiv!((typeof(oneunit(eltype(D))/oneunit(eltype(A)))).(A), D)

(\)(F::Factorization, D::Diagonal) =
    ldiv!(F, Matrix{typeof(oneunit(eltype(D))/oneunit(eltype(F)))}(D))
\(adjF::Adjoint{<:Any,<:Factorization}, D::Diagonal) =
    (F = adjF.parent; ldiv!(adjoint(F), Matrix{typeof(oneunit(eltype(D))/oneunit(eltype(F)))}(D)))
(\)(A::Union{QR,QRCompactWY,QRPivoted}, B::Diagonal) =
    invoke(\, Tuple{Union{QR,QRCompactWY,QRPivoted}, AbstractVecOrMat}, A, B)

function kron(A::Diagonal{T1}, B::Diagonal{T2}) where {T1<:Number, T2<:Number}
    valA = A.diag; nA = length(valA)
    valB = B.diag; nB = length(valB)
    valC = Vector{typeof(zero(T1)*zero(T2))}(undef,nA*nB)
    @inbounds for i = 1:nA, j = 1:nB
        valC[(i-1)*nB+j] = valA[i] * valB[j]
    end
    return Diagonal(valC)
end

function kron(A::Diagonal{T}, B::AbstractMatrix{S}) where {T<:Number, S<:Number}
    Base.require_one_based_indexing(B)
    (mA, nA) = size(A); (mB, nB) = size(B)
    R = zeros(Base.promote_op(*, T, S), mA * mB, nA * nB)
    m = 1
    for j = 1:nA
        A_jj = A[j,j]
        for k = 1:nB
            for l = 1:mB
                R[m] = A_jj * B[l,k]
                m += 1
            end
            m += (nA - 1) * mB
        end
        m += mB
    end
    return R
end

function kron(A::AbstractMatrix{T}, B::Diagonal{S}) where {T<:Number, S<:Number}
    require_one_based_indexing(A)
    (mA, nA) = size(A); (mB, nB) = size(B)
    R = zeros(promote_op(*, T, S), mA * mB, nA * nB)
    m = 1
    for j = 1:nA
        for l = 1:mB
            Bll = B[l,l]
            for k = 1:mA
                R[m] = A[k,j] * Bll
                m += nB
            end
            m += 1
        end
        m -= nB
    end
    return R
end

conj(D::Diagonal) = Diagonal(conj(D.diag))
transpose(D::Diagonal{<:Number}) = D
transpose(D::Diagonal) = Diagonal(transpose.(D.diag))
adjoint(D::Diagonal{<:Number}) = conj(D)
adjoint(D::Diagonal) = Diagonal(adjoint.(D.diag))

function diag(D::Diagonal, k::Integer=0)
    # every branch call similar(..., ::Int) to make sure the
    # same vector type is returned independent of k
    if k == 0
        return copyto!(similar(D.diag, length(D.diag)), D.diag)
    elseif -size(D,1) <= k <= size(D,1)
        return fill!(similar(D.diag, size(D,1)-abs(k)), 0)
    else
        throw(ArgumentError(string("requested diagonal, $k, must be at least $(-size(D, 1)) ",
            "and at most $(size(D, 2)) for an $(size(D, 1))-by-$(size(D, 2)) matrix")))
    end
end
tr(D::Diagonal) = sum(tr, D.diag)
det(D::Diagonal) = prod(det, D.diag)
logdet(D::Diagonal{<:Real}) = sum(log, D.diag)
function logdet(D::Diagonal{<:Complex}) # make sure branch cut is correct
    z = sum(log, D.diag)
    complex(real(z), rem2pi(imag(z), RoundNearest))
end

# Matrix functions
for f in (:exp, :log, :sqrt,
          :cos, :sin, :tan, :csc, :sec, :cot,
          :cosh, :sinh, :tanh, :csch, :sech, :coth,
          :acos, :asin, :atan, :acsc, :asec, :acot,
          :acosh, :asinh, :atanh, :acsch, :asech, :acoth)
    @eval $f(D::Diagonal) = Diagonal($f.(D.diag))
end

#Linear solver
function ldiv!(D::Diagonal, B::StridedVecOrMat)
    m, n = size(B, 1), size(B, 2)
    if m != length(D.diag)
        throw(DimensionMismatch("diagonal matrix is $(length(D.diag)) by $(length(D.diag)) but right hand side has $m rows"))
    end
    (m == 0 || n == 0) && return B
    for j = 1:n
        for i = 1:m
            di = D.diag[i]
            if di == 0
                throw(SingularException(i))
            end
            B[i,j] = di \ B[i,j]
        end
    end
    return B
end
(\)(D::Diagonal, A::AbstractMatrix) =
    ldiv!(D, (typeof(oneunit(eltype(D))/oneunit(eltype(A)))).(A))

(\)(D::Diagonal, b::AbstractVector) = D.diag .\ b
(\)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag .\ Db.diag)

function inv(D::Diagonal{T}) where T
    Di = similar(D.diag, typeof(inv(zero(T))))
    for i = 1:length(D.diag)
        if D.diag[i] == zero(T)
            throw(SingularException(i))
        end
        Di[i] = inv(D.diag[i])
    end
    Diagonal(Di)
end

function pinv(D::Diagonal{T}) where T
    Di = similar(D.diag, typeof(inv(zero(T))))
    for i = 1:length(D.diag)
        isfinite(inv(D.diag[i])) ? Di[i]=inv(D.diag[i]) : Di[i]=zero(T)
    end
    Diagonal(Di)
end
function pinv(D::Diagonal{T}, tol::Real) where T
    Di = similar(D.diag, typeof(inv(zero(T))))
    if( !isempty(D.diag) ) maxabsD = maximum(abs.(D.diag)) end
    for i = 1:length(D.diag)
        if( abs(D.diag[i]) > tol*maxabsD && isfinite(inv(D.diag[i])) )
            Di[i]=inv(D.diag[i])
        else
            Di[i]=zero(T)
        end
    end
    Diagonal(Di)
end

#Eigensystem
eigvals(D::Diagonal{<:Number}; permute::Bool=true, scale::Bool=true) = D.diag
eigvals(D::Diagonal; permute::Bool=true, scale::Bool=true) =
    [eigvals(x) for x in D.diag] #For block matrices, etc.
eigvecs(D::Diagonal) = Matrix{eltype(D)}(I, size(D))
function eigen(D::Diagonal; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=nothing)
    if any(!isfinite, D.diag)
        throw(ArgumentError("matrix contains Infs or NaNs"))
    end
    Eigen(sorteig!(eigvals(D), eigvecs(D), sortby)...)
end

#Singular system
svdvals(D::Diagonal{<:Number}) = sort!(abs.(D.diag), rev = true)
svdvals(D::Diagonal) = [svdvals(v) for v in D.diag]
function svd(D::Diagonal{<:Number})
    S   = abs.(D.diag)
    piv = sortperm(S, rev = true)
    U   = Diagonal(D.diag ./ S)
    Up  = hcat([U[:,i] for i = 1:length(D.diag)][piv]...)
    V   = Diagonal(fill!(similar(D.diag), one(eltype(D.diag))))
    Vp  = hcat([V[:,i] for i = 1:length(D.diag)][piv]...)
    return SVD(Up, S[piv], copy(Vp'))
end

# disambiguation methods: * of Diagonal and Adj/Trans AbsVec
*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal) = Adjoint(map((t,s) -> t'*s, D.diag, parent(x)))
*(x::Adjoint{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) =
    mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y))
*(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal) = Transpose(map((t,s) -> transpose(t)*s, D.diag, parent(x)))
*(x::Transpose{<:Any,<:AbstractVector}, D::Diagonal, y::AbstractVector) =
    mapreduce(t -> t[1]*t[2]*t[3], +, zip(x, D.diag, y))

function cholesky!(A::Diagonal, ::Val{false} = Val(false); check::Bool = true)
    info = 0
    diagonal = A.diag
    for i in axes(diagonal, 1)
        d = diagonal[i]
        if !(d == 0 || (isreal(d) && d < 0))
            diagonal[i] = √d
        elseif check
            throw(PosDefException(i))
        else
            info = i
            break
        end
    end
    Cholesky(A, 'U', convert(BlasInt, info))
end

cholesky(A::Diagonal, ::Val{false} = Val(false); check::Bool = true) =
    cholesky!(cholcopy(A), Val(false); check = check)

max_stars_repo_name: gasagna/julia

max_issues_repo_name: gasagna/julia

max_forks_repo_name: gasagna/julia


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -23.82106505746807 <= q <= 20.87497636142748)

                     add_NL_constraint(m, :(log(1 + exp(-0.7432570628691102 + 0.9028765283094948*log(1 + exp(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))) + -0.9177729976456397*log(1 + exp(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))) + -0.6091404444372248*log(1 + exp(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))) + 0.7075296204914552*log(1 + exp(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))) + -0.9115635329002854*log(1 + exp(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))))) + log(1 + exp(-0.2671327260424561 + -0.018499351907995898*log(1 + exp(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))) + 0.12529725503678657*log(1 + exp(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))) + 0.9839396030143157*log(1 + exp(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))) + 0.28500560463058955*log(1 + exp(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))) + 0.970334521625404*log(1 + exp(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))))) + log(1 + exp(0.40040910594674495 + -0.7454849743358074*log(1 + exp(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))) + 0.7108986625228169*log(1 + exp(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))) + -0.30349624151173904*log(1 + exp(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))) + -0.26057109207128715*log(1 + exp(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))) + 0.7858005505148937*log(1 + exp(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))))) + log(1 + exp(-0.7120753920031948 + -0.005433419997351141*log(1 + exp(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))) + -0.8839514190404207*log(1 + exp(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))) + -0.2933551592941672*log(1 + exp(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))) + -0.5418193939730438*log(1 + exp(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))) + 0.8978172795430157*log(1 + exp(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))))) + log(1 + exp(0.7950909223530527 + -0.6380481687930124*log(1 + exp(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3]))) + -0.8467501205760102*log(1 + exp(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3]))) + 0.5001198683847332*log(1 + exp(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3]))) + -0.03855272462979853*log(1 + exp(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3]))) + 0.3301524112226937*log(1 + exp(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

ext: jl

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "eccodes"
version = v"2.19.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://confluence.ecmwf.int/download/attachments/45757960/eccodes-$version-Source.tar.gz", "9964bed5058e873d514bd4920951122a95963128b12f55aa199d9afbafdd5d4b"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/eccodes-*-Source
if [[ ${target} = *-mingw* ]] ; then
    chmod +x cmake/ecbuild_windows_replace_symlinks.sh
    atomic_patch -p1 /workspace/srcdir/patches/windows.patch
else
    atomic_patch -p1 /workspace/srcdir/patches/unix.patch
fi
mkdir build
cd build
export CFLAGS="-I${includedir}"
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DENABLE_NETCDF=OFF \
    -DENABLE_PNG=ON \
    -DENABLE_PYTHON=OFF \
    -DENABLE_FORTRAN=OFF \
    -DENABLE_ECCODES_THREADS=ON \
    ..
make -j${nproc}
make install
install_license ../LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libeccodes", :eccodes)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="libpng_jll", uuid="b53b4c65-9356-5827-b1ea-8c7a1a84506f"))
    Dependency(PackageSpec(name="OpenJpeg_jll", uuid="643b3616-a352-519d-856d-80112ee9badc"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
cd(@__DIR__) # changes the directory to the current directory, the default I guess is the HOME
using Pkg; Pkg.activate("."); Pkg.instantiate()
#=
Pkg is Julia's built-in package manager, and handles operations such as installing, updating and removing packages.
Just like cargo it creates a toml-file that describes the versions of the software used & dependencies.
=#
Pkg.add("OrdinaryDiffEq")
Pkg.add("ModelingToolkit")
Pkg.add("DataDrivenDiffEq")
Pkg.add("LinearAlgebra")
Pkg.add("DiffEqSensitivity")
Pkg.add("Optim")
Pkg.add("DiffEqFlux")
Pkg.add("Flux")
Pkg.add("Plots")
using OrdinaryDiffEq
using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra, DiffEqSensitivity, Optim
using DiffEqFlux, Flux
using Plots
gr() # function from Plots

function corona!(du,u,p,t)
    S,E,I,R,N,D,C = u
    F, β0,α,κ,μ,σ,γ,d,λ = p
    dS = -β0*S*F/N - β(t,β0,D,N,κ,α)*S*I/N -μ*S # susceptible
    dE = β0*S*F/N + β(t,β0,D,N,κ,α)*S*I/N -(σ+μ)*E # exposed
    dI = σ*E - (γ+μ)*I # infected
    dR = γ*I - μ*R # removed (recovered + dead)
    dN = -μ*N # total population
    dD = d*γ*I - λ*D # severe, critical cases, and deaths
    dC = σ*E # +cumulative cases

    du[1] = dS; du[2] = dE; du[3] = dI; du[4] = dR
    du[5] = dN; du[6] = dD; du[7] = dC
end
β(t,β0,D,N,κ,α) = β0*(1-α)*(1-D/N)^κ
S0 = 14e6
u0 = [0.9*S0, 0.0, 0.0, 0.0, S0, 0.0, 0.0]
p_ = [10.0, 0.5944, 0.4239, 1117.3, 0.02, 1/3, 1/5,0.2, 1/11.2]
R0 = p_[2]/p_[7]*p_[6]/(p_[6]+p_[5])
tspan = (0.0, 21.0)
prob = ODEProblem(corona!, u0, tspan, p_)
solution = solve(prob, Vern7(), abstol=1e-12, reltol=1e-12, saveat = 1)

tspan2 = (0.0,60.0)
prob = ODEProblem(corona!, u0, tspan2, p_)
solution_extrapolate = solve(prob, Vern7(), abstol=1e-12, reltol=1e-12, saveat = 1)

# Ideal data
tsdata = Array(solution)
# Add noise to the data
noisy_data = tsdata + Float32(1e-5)*randn(eltype(tsdata), size(tsdata))

plot(abs.(tsdata-noisy_data)')

### Neural ODE

ann_node = FastChain(FastDense(7, 64, tanh),FastDense(64, 64, tanh), FastDense(64, 64, tanh), FastDense(64, 7))
p = Float64.(initial_params(ann_node))

function dudt_node(u,p,t)
    S,E,I,R,N,D,C = u
    F,β0,α,κ,μ,σ,γ,d,λ = p_
    dS,dE,dI,dR,dD = ann_node([S/N,E,I,R,N,D/N,C],p)

    dN = -μ*N # total population
    dC = σ*E # +cumulative cases

    [dS,dE,dI,dR,dN,dD,dC]
end
prob_node = ODEProblem(dudt_node, u0, tspan, p)
s = concrete_solve(prob_node, Tsit5(), u0, p, saveat = solution.t)

function predict(θ)
    Array(concrete_solve(prob_node, Vern7(), u0, θ, saveat = 1,
                         abstol=1e-6, reltol=1e-6,
                         sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
end

# No regularisation right now
function loss(θ)
    pred = predict(θ)
    sum(abs2, (noisy_data[2:4,:] .- pred[2:4,:])), pred # + 1e-5*sum(sum.(abs, params(ann)))
end

loss(p)

const losses = []
callback(θ,l,pred) = begin
    push!(losses, l)
    if length(losses)%50==0
        println(losses[end])
    end
    false
end

res1_node = DiffEqFlux.sciml_train(loss, p, ADAM(0.01), cb=callback, maxiters = 500)
res2_node = DiffEqFlux.sciml_train(loss, res1_node.minimizer, BFGS(initial_stepnorm=0.01), cb=callback, maxiters = 10000)

prob_node2 = ODEProblem(dudt_node, u0, tspan, res2_node.minimizer)
s = solve(prob_node2, Tsit5(), saveat = 1)
scatter(solution, vars=[2,3,4], label=["True Exposed" "True Infected" "True Recovered"])
plot!(s, vars=[2,3,4], label=["Estimated Exposed" "Estimated Infected" "Estimated Recovered"])

# Plot the losses
plot(losses, yaxis = :log, xaxis = :log, xlabel = "Iterations", ylabel = "Loss")

# Extrapolate out
prob_node_extrapolate = ODEProblem(dudt_node,u0, tspan2, res2_node.minimizer)
_sol_node = solve(prob_node_extrapolate, Vern7(), abstol=1e-12, reltol=1e-12, saveat = 1)
p_node = scatter(solution_extrapolate, vars=[2,3,4], legend = :topleft, label=["True Exposed" "True Infected" "True Recovered"], title="Neural ODE Extrapolation")
plot!(p_node,_sol_node, lw=5, vars=[2,3,4], label=["Estimated Exposed" "Estimated Infected" "Estimated Recovered"])
plot!(p_node,[20.99,21.01],[0.0,maximum(hcat(Array(solution_extrapolate[2:4,:]),Array(_sol_node[2:4,:])))],lw=5,color=:black,label="Training Data End")

savefig("neuralode_extrapolation.png")
savefig("neuralode_extrapolation.pdf")

### Universal ODE Part 1

ann = FastChain(FastDense(3, 64, tanh),FastDense(64, 64, tanh), FastDense(64, 1))
p = Float64.(initial_params(ann))

function dudt_(u,p,t)
    S,E,I,R,N,D,C = u
    F, β0,α,κ,μ,σ,γ,d,λ = p_
    z = ann([S/N,I,D/N],p) # Exposure does not depend on exposed, removed, or cumulative!
    dS = -β0*S*F/N - z[1] -μ*S # susceptible
    dE = β0*S*F/N + z[1] -(σ+μ)*E # exposed
    dI = σ*E - (γ+μ)*I # infected
    dR = γ*I - μ*R # removed (recovered + dead)
    dN = -μ*N # total population
    dD = d*γ*I - λ*D # severe, critical cases, and deaths
    dC = σ*E # +cumulative cases

    [dS,dE,dI,dR,dN,dD,dC]
end
prob_nn = ODEProblem(dudt_,u0, tspan, p)
s = concrete_solve(prob_nn, Tsit5(), u0, p, saveat = 1)

plot(solution, vars=[2,3,4])
plot!(s[2:4,:]')

function predict(θ)
    Array(concrete_solve(prob_nn, Vern7(), u0, θ, saveat = solution.t,
                         abstol=1e-6, reltol=1e-6,
                         sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
end

# No regularisation right now
function loss(θ)
    pred = predict(θ)
    sum(abs2, noisy_data[2:4,:] .- pred[2:4,:]), pred # + 1e-5*sum(sum.(abs, params(ann)))
end

loss(p)

const losses = []
callback(θ,l,pred) = begin
    push!(losses, l)
    if length(losses)%50==0
        println(losses[end])
    end
    false
end

res1_uode = DiffEqFlux.sciml_train(loss, p, ADAM(0.01), cb=callback, maxiters = 500)
res2_uode = DiffEqFlux.sciml_train(loss, res1_uode.minimizer, BFGS(initial_stepnorm=0.01), cb=callback, maxiters = 10000)

loss(res2_uode.minimizer)

prob_nn2 = ODEProblem(dudt_,u0, tspan, res2_uode.minimizer)
uode_sol = solve(prob_nn2, Tsit5(), saveat = 1)
plot(solution, vars=[2,3,4])
plot!(uode_sol, vars=[2,3,4])

# Plot the losses
plot(losses, yaxis = :log, xaxis = :log, xlabel = "Iterations", ylabel = "Loss")

# Collect the state trajectory and the derivatives
X = noisy_data
# Ideal derivatives
DX = Array(solution(solution.t, Val{1}))

# Extrapolate out
prob_nn2 = ODEProblem(dudt_,u0, tspan2, res2_uode.minimizer)
_sol_uode = solve(prob_nn2, Vern7(), abstol=1e-12, reltol=1e-12, saveat = 1)
p_uode = scatter(solution_extrapolate, vars=[2,3,4], legend = :topleft, label=["True Exposed" "True Infected" "True Recovered"], title="Universal ODE Extrapolation")
plot!(p_uode,_sol_uode, lw = 5, vars=[2,3,4], label=["Estimated Exposed" "Estimated Infected" "Estimated Recovered"])
plot!(p_uode,[20.99,21.01],[0.0,maximum(hcat(Array(solution_extrapolate[2:4,:]),Array(_sol_uode[2:4,:])))],lw=5,color=:black,label="Training Data End")

savefig("universalode_extrapolation.png")
savefig("universalode_extrapolation.pdf")

### Universal ODE Part 2: SInDy to Equations

# Create a Basis
@variables u[1:3]
# Lots of polynomials
polys = Operation[]
for i ∈ 0:2, j ∈ 0:2, k ∈ 0:2
    push!(polys, u[1]^i * u[2]^j * u[3]^k)
end

# And some other stuff
h = [cos.(u)...; sin.(u)...; unique(polys)...]
basis = Basis(h, u)

X = noisy_data
# Ideal derivatives
DX = Array(solution(solution.t, Val{1}))
S,E,I,R,N,D,C = eachrow(X)
F,β0,α,κ,μ,_,γ,d,λ = p_
L = β.(0:tspan[end],β0,D,N,κ,α).*S.*I./N
L̂ = vec(ann([S./N I D./N]',res2_uode.minimizer))
X̂ = [S./N I D./N]'

scatter(L,title="Estimated vs Expected Exposure Term",label="True Exposure")
plot!(L̂,label="Estimated Exposure")
savefig("estimated_exposure.png")
savefig("estimated_exposure.pdf")

# Create an optimizer for the SINDY problem
opt = SR3()
# Create the thresholds which should be used in the search process
thresholds = exp10.(-6:0.1:1)

# Test on original data and without further knowledge
Ψ_direct = SInDy(X[2:4, :], DX[2:4, :], basis, thresholds, opt = opt, maxiter = 50000) # Fail
println(Ψ_direct.basis)
# Test on ideal derivative data ( not available )
Ψ_ideal = SInDy(X[2:4, 5:end], L[5:end], basis, thresholds, opt = opt, maxiter = 50000) # Succeed
println(Ψ_ideal.basis)
# Test on uode derivative data
Ψ = SInDy(X̂[:, 2:end], L̂[2:end], basis, thresholds,  opt = opt, maxiter = 10000, normalize = true, denoise = true) # Succeed
println(Ψ.basis)

# Build a ODE for the estimated system
function approx(u,p,t)
    S,E,I,R,N,D,C = u
    F, β0,α,κ,μ,σ,γ,d,λ = p_
    z = Ψ([S/N,I,D/N]) # Exposure does not depend on exposed, removed, or cumulative!
    dS = -β0*S*F/N - z[1] -μ*S # susceptible
    dE = β0*S*F/N + z[1] -(σ+μ)*E # exposed
    dI = σ*E - (γ+μ)*I # infected
    dR = γ*I - μ*R # removed (recovered + dead)
    dN = -μ*N # total population
    dD = d*γ*I - λ*D # severe, critical cases, and deaths
    dC = σ*E # +cumulative cases

    [dS,dE,dI,dR,dN,dD,dC]
end

# Create the approximated problem and solution
a_prob = ODEProblem{false}(approx, u0, tspan2, p_)
a_solution = solve(a_prob, Tsit5())

p_uodesindy = scatter(solution_extrapolate, vars=[2,3,4], legend = :topleft, label=["True Exposed" "True Infected" "True Recovered"])
plot!(p_uodesindy,a_solution, lw = 5, vars=[2,3,4], label=["Estimated Exposed" "Estimated Infected" "Estimated Recovered"])
plot!(p_uodesindy,[20.99,21.01],[0.0,maximum(hcat(Array(solution_extrapolate[2:4,:]),Array(_sol_uode[2:4,:])))],lw=5,color=:black,label="Training Data End")

savefig("universalodesindy_extrapolation.png")
savefig("universalodesindy_extrapolation.pdf")


text:
using Distributions, Statistics, LinearAlgebra

import Base: size, length

import Statistics: mean, var, std

export initialize, initialize!

function initialize(N::Int, Distx0::MultivariateDistribution)
    NS = length(Distx0)

    # pre-allocate space
    ENS = EnsembleState(N, zeros(NS))

    ENS.S .= [rand(Distx0) for i = 1:N]
    return ENS
end
#
function initialize(N::Int, NS::Int)
    # pre-allocate space
    ENS = EnsembleState(N, zeros(NS))

    Dist = MvNormal(zeros(NS), I)
    ENS.S .= [rand(Dist) for i = 1:N]
    return ENS
end
#
function initialize!(ENS::EnsembleState{N, TS}, u::TS)  where {N, TS}
    NS = length(u)
    Dist = MvNormal(zeros(NS), I)
    ENS.S .= [rand(Dist) for i = 1:N]
    return ENS
end

max_stars_repo_name: mleprovost/EnKF

max_issues_repo_name: mleprovost/EnKF

max_forks_repo_name: mleprovost/EnKF


text:
####
#### CLI settings
####

const ALLOWED_FILE_SUFFIXES = (".mat", ".nii", ".nii.gz", ".par", ".xml", ".rec")
const ALLOWED_FILE_SUFFIXES_STRING = join(ALLOWED_FILE_SUFFIXES, ", ", ", and ")

const CLI_SETTINGS = ArgParseSettings(
    prog = "",
    fromfile_prefix_chars = "@",
    error_on_conflict = false,
    exit_after_help = false,
)

@add_arg_table! CLI_SETTINGS begin
    "input"
        nargs = '+' # At least one input is required
        arg_type = String
        required = true
        help = "one or more input filenames. Valid file types are limited to: $ALLOWED_FILE_SUFFIXES_STRING"
    "--mask", "-m"
        nargs = '+' # If --mask is passed, at least one input is required
        arg_type = String
        help = "one or more mask filenames. Masks are loaded and subsequently applied to the corresponding input files via elementwise multiplication. The number of mask files must equal the number of input files. Valid file types are the same as for input files, and are limited to: $ALLOWED_FILE_SUFFIXES_STRING"
    "--output", "-o"
        arg_type = String
        help = "output directory. If not specified, output file(s) will be stored in the same location as the corresponding input file(s). Outputs are stored with the same basename as the input files with additional suffixes; see --T2map and --T2part"
    "--T2map"
        action = :store_true
        help = "call T2mapSEcorr to compute T2 distributions from 4D multi spin-echo input images. T2 distributions and T2 maps produced by T2mapSEcorr are saved as MAT files with extensions .t2dist.mat and .t2maps.mat"
    "--T2part"
        action = :store_true
        help = "call T2partSEcorr to analyze 4D T2 distributions to produce parameter maps. If --T2map is also passed, input 4D arrays are interpreted as multi spin-echo images and T2 distributions are first computed by T2mapSEcorr. If only --T2part is passed, input 4D arrays are interpreted as T2 distributions and only T2partSEcorr is called. Output T2 parts are saved as a MAT file with extension .t2parts.mat"
    "--quiet", "-q"
        action = :store_true
        help = "suppress printing to the terminal. Note: all terminal outputs, including errors and warnings, are still printed to the log file"
    "--dry"
        action = :store_true
        help = "execute dry run of processing without saving any results"
    "--legacy"
        action = :store_true
        help = "use legacy settings and algorithms from the original MATLAB pipeline. This ensures that the same T2-distributions and T2-parts will be produced as those from MATLAB. Note that execution time will be much slower, and less robust algorithms will be used."
end

add_arg_group!(CLI_SETTINGS,
    "T2mapSEcorr/T2partSEcorr arguments",
    :t2_map_part,
)

@add_arg_table! CLI_SETTINGS begin
    "--MatrixSize"
        nargs = 3
        arg_type = Int
        help = "first three dimensions of input 4D image. Inferred automatically"
        group = :t2_map_part
    "--nTE"
        arg_type = Int
        help = "number of echoes in input signal. Inferred automatically when --T2map is passed"
        group = :t2_map_part
    "--TE"
        arg_type = Float64
        help = "inter-echo spacing (Units: seconds). Required when --T2map is passed"
        group = :t2_map_part
    "--nT2"
        arg_type = Int
        help = "number of T2 components used in the multi-exponential analysis. Required when --T2map is passed. Inferred from fourth dimension of input image if only --T2part is passed"
        group = :t2_map_part
    "--T2Range"
        nargs = 2
        arg_type = Float64
        help = "minimum and maximum T2 values (Units: seconds). T2 components are logarithmically spaced between these bounds. Required parameter."
        group = :t2_map_part
    "--SPWin"
        nargs = 2
        arg_type = Float64
        help = "minimum and maximum T2 values of the short peak window (Units: seconds). Required parameter when --T2part is passed"
        group = :t2_map_part
    "--MPWin"
        nargs = 2
        arg_type = Float64
        help = "minimum and maximum T2 values of the middle peak window (Units: seconds). Required parameter when --T2part is passed"
        group = :t2_map_part
    "--T1"
        arg_type = Float64
        help = "assumed value of longitudinal T1 relaxation (Units: seconds)."
        group = :t2_map_part
    "--Reg"
        arg_type = String
        help = "routine used for choosing regularization parameter. One of \"none\", \"chi2\", \"gcv\", or \"lcurve\", representing no regularization, --Chi2Factor based Tikhonov regularization, generalized cross-validation based regularization, and L-curve based regularization, respectively."
        group = :t2_map_part
    "--Chi2Factor"
        arg_type = Float64
        help = "if --Reg=\"chi2\", the T2 distribution is regularized such that the chi^2 goodness of fit is increased by a multiplicative factor --Chi2Factor relative to the unregularized solution"
        group = :t2_map_part
    "--Sigmoid"
        arg_type = Float64
        help = "replace the hard upper limit cutoff time of the short peak window, --SPWin[2], with a smoothed sigmoidal cutoff function, f, scaled and shifted such that f(--SPWin[2] +/- --Sigmoid) = 0.5 -/+ 0.4. --Sigmoid is the time scale of the smoothing (Units: seconds)"
        group = :t2_map_part
    "--Threshold"
        arg_type = Float64
        help = "first echo intensity cutoff for empty voxels. Processing is skipped for voxels with intensity <= --Threshold"
        group = :t2_map_part
end

add_arg_group!(CLI_SETTINGS,
    "B1 correction and stimulated echo correction",
    :B1_SE_corr,
)

@add_arg_table! CLI_SETTINGS begin
    "--nRefAngles"
        arg_type = Int
        help = "in estimating the local refocusing flip angle to correct for B1 inhomogeneities, up to --nRefAngles angles in the range [--MinRefAngle, 180] are explicitly checked. The optimal angle is then determined through interpolation from these --nRefAngles observations."
        group = :B1_SE_corr
    "--nRefAnglesMin"
        arg_type = Int
        help = "initial number of angles to check during flip angle estimation before refinement near likely optima. Setting --nRefAnglesMin equal to --nRefAngles forces all angles to be checked."
        group = :B1_SE_corr
    "--MinRefAngle"
        arg_type = Float64
        help = "minimum refocusing angle for flip angle estimation (Units: degrees)."
        group = :B1_SE_corr
    "--SetFlipAngle"
        arg_type = Float64
        help = "to skip B1 inhomogeneity correction, use --SetFlipAngle to assume a fixed refocusing flip angle for all voxels (Units: degrees)."
        group = :B1_SE_corr
    "--SetRefConAngle"
        arg_type = Float64
        help = "refocusing pulse control angle for stimulated echo correction. Unlike B1 inhomogeneity correction, stimulated echo correction must be performed manually. By default, --SetRefConAngle is set to 180 degrees, equivalent to no stimulated echo correction (Units: degrees)."
        group = :B1_SE_corr
    "--RefConAngle"
        arg_type = Float64
        help = "Deprecated flag; see --SetRefConAngle."
        group = :B1_SE_corr
end

add_arg_group!(CLI_SETTINGS,
    "Additional save options",
    :save_opts,
)

@add_arg_table! CLI_SETTINGS begin
    "--SaveDecayCurve"
        action = :store_true
        help = "include a 4D array of the time domain decay curves resulting from the NNLS fits in the output maps dictionary"
        group = :save_opts
    "--SaveNNLSBasis"
        action = :store_true
        help = "include a 5D (or 2D if --SetFlipAngle is used) array of NNLS basis matrices in the output maps dictionary. Note: this 5D array is extremely large for typical image sizes; in most cases, this flag should only be set when debugging small datasets"
        group = :save_opts
    "--SaveRegParam"
        action = :store_true
        help = "include 3D arrays of the regularization parameters and resulting chi^2-factors in the output maps dictionary"
        group = :save_opts
    "--SaveResidualNorm"
        action = :store_true
        help = "include a 3D array of the l2-norms of the residuals from the NNLS fits in the output maps dictionary"
        group = :save_opts
end

add_arg_group!(CLI_SETTINGS,
    "BET arguments",
    :bet_args,
)

@add_arg_table! CLI_SETTINGS begin
    "--bet"
        action = :store_true
        help = "use the BET brain extraction tool from the FSL library of analyis tools to automatically create a binary brain mask. Only voxels within the binary mask will be analyzed. Note that if a mask is passed explicitly with the --mask flag, this mask will be used and the --bet flag will be ignored"
        group = :bet_args
    "--betargs"
        arg_type = String
        default = "-m -n -f 0.25 -R"
        help = "BET command line interface arguments. Must be passed as a single string with arguments separated by spaces, e.g. '-m -n'. The flag '-m' indicates that a binary mask should be computed, and therefore will be added to the list of arguments if not provided"
        group = :bet_args
    "--betpath"
        arg_type = String
        default = "bet"
        help = "path to BET executable."
        group = :bet_args
end

"""
    main(command_line_args = ARGS)

Entry point function for command line interface, parsing the command line arguments `ARGS` and subsequently calling one or both of `T2mapSEcorr` and `T2partSEcorr` with the parsed settings.
See the [Arguments](@ref) section for available options.

See also:
* [`T2mapSEcorr`](@ref)
* [`T2partSEcorr`](@ref)
"""
function main(command_line_args::Vector{String} = ARGS)

    # Parse command line arguments
    opts = parse_args(command_line_args, CLI_SETTINGS; as_symbols = true)
    opts === nothing && return # Help message was triggered. Return nothing instead of exit(0)
    opts = handle_cli_deprecations!(opts)

    # Get input file list and output folder list
    for file_info in get_file_infos(opts)
        # Make output path
        if !opts[:dry]
            mkpath(file_info[:outputfolder])
        end

        # Save settings files
        if !opts[:dry]
            map(filter(s -> startswith(s, "@"), command_line_args)) do settingsfile
                src = settingsfile[2:end] # drop "@" character
                dst = joinpath(file_info[:outputfolder], file_info[:choppedinputfile] * "." * basename(src))
                cp(src, dst; force = true)
            end
        end

        # Main processing
        tee_capture(
            logfile = joinpath(file_info[:outputfolder], file_info[:choppedinputfile] * ".log"),
            suppress_terminal = opts[:quiet],
            suppress_logfile = opts[:dry],
        ) do io
            try
                main_(io, file_info, opts)
            catch e
                @warn "Error during processing of file: $(file_info[:inputfile])"
                @warn sprint(showerror, e, catch_backtrace())
            end
        end
    end

    return nothing
end

function main_(io::IO, file_info::Dict{Symbol,Any}, opts::Dict{Symbol,Any})

    # Starting message/starting time
    t_start = tic()
    printheader(io, "Starting DECAES with $(Threads.nthreads()) threads")

    # Load image(s)
    image = @showtime(io,
        "Loading input file: $(file_info[:inputfile])",
        load_image(file_info[:inputfile]),
    )

    # Apply mask
    if file_info[:maskfile] !== nothing
        @showtime(io,
            "Applying mask from file: '$(file_info[:maskfile])'",
            try_apply_maskfile!(image, file_info[:maskfile]),
        )
    elseif opts[:bet]
        @showtime(io,
            "Making and applying BET mask with args: '$(opts[:betargs])'",
            try_apply_bet!(image, opts[:betpath], opts[:betargs]),
        )
    end

    # Compute T2 distribution from input 4D multi-echo image
    if opts[:T2map]
        maps, dist = @showtime(io,
            "Running T2mapSEcorr on file: $(file_info[:inputfile])",
            T2mapSEcorr(io, image, t2map_options(image, opts)),
        )

        # Save T2-distribution to .mat file
        savefile = joinpath(file_info[:outputfolder], file_info[:choppedinputfile] * ".t2dist.mat")
        if !opts[:dry]
            @showtime(io,
                "Saving T2 distribution to file: $savefile",
                MAT.matwrite(savefile, Dict{String,Any}("dist" => dist)),
            )
        end

        # Save T2-maps to .mat file
        savefile = joinpath(file_info[:outputfolder], file_info[:choppedinputfile] * ".t2maps.mat")
        if !opts[:dry]
            @showtime(io,
                "Saving T2 parameter maps to file: $savefile",
                MAT.matwrite(savefile, maps),
            )
        end
    else
        # Input image is the T2 distribution
        dist = image
    end

    # Analyze T2 distribution to produce parameter maps
    if opts[:T2part]
        parts = @showtime(io,
            "Running T2partSEcorr",
            T2partSEcorr(io, dist, t2part_options(dist, opts)),
        )

        # Save T2-parts to .mat file
        savefile = joinpath(file_info[:outputfolder], file_info[:choppedinputfile] * ".t2parts.mat")
        if !opts[:dry]
            @showtime(io,
                "Saving T2 parts maps to file: $savefile",
                MAT.matwrite(savefile, parts),
            )
        end
    end

    # Done message
    printheader(io, "Finished ($(round(toc(t_start); digits = 2)) seconds)")

    return nothing
end

"""
Entrypoint function for compiling DECAES into an executable [app](https://julialang.github.io/PackageCompiler.jl/dev/apps/).
"""
function julia_main()::Cint
    try
        main(ARGS)
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

####
#### Helper functions
####

function handle_cli_deprecations!(opts)
    handle_renamed_cli_flag!(opts, :RefConAngle => :SetRefConAngle)
end

function handle_renamed_cli_flag!(opts, oldnew::Pair{Symbol, Symbol})
    oldflag, newflag = oldnew
    if opts[oldflag] !== nothing
        if opts[newflag] !== nothing
            error("The flag --$newflag and the deprecated flag --$oldflag were both passed; use --$newflag only.")
        else
            @warn "The flag --$oldflag is deprecated and will be removed in future releases; use --$newflag instead."
        end
        opts[newflag] = opts[oldflag]
        delete!(opts, oldflag)
    end
    return opts
end

function t2map_options(image::Array, opts::Dict{Symbol,Any})
    fields = fieldsof(T2mapOptions, Set)
    kwargs = Dict{Symbol,Any}()
    for (k,v) in opts
        (v === nothing) && continue # filter unset cli args
        (v isa AbstractVector && isempty(v)) && continue # filter unset cli args (empty vectors are unset cli varargs)
        (k ∉ fields) && continue # filter T2mapOptions fields
        kwargs[k] = v isa AbstractVector ? tuple(v...) : v
    end
    T2mapOptions(image; kwargs...)
end

function t2part_options(dist::Array, opts::Dict{Symbol,Any})
    fields = fieldsof(T2partOptions, Set)
    kwargs = Dict{Symbol,Any}()
    for (k,v) in opts
        (v === nothing) && continue # filter unset cli args
        (v isa AbstractVector && isempty(v)) && continue # filter unset cli args (empty vectors are unset cli varargs)
        (k === :nT2 && !opts[:T2map]) && continue # nT2 must be explicitly passed, unless not performing T2-mapping, in which case it is inferred from `dist`
        (k ∉ fields) && continue # filter T2mapOptions fields
        kwargs[k] = v isa AbstractVector ? tuple(v...) : v
    end
    T2partOptions(dist; kwargs...)
end

function get_file_infos(opts::Dict{Symbol,Any})
    @unpack input, output, mask = opts

    # Read in input files
    inputfiles = String[path for path in input if is_allowed_suffix(path)]

    if isempty(inputfiles)
        msg = if !isempty(input) && isfile(input[1])
            "No valid file types were found for processing, but a file name was passed.\n" *
            "Perhaps you meant to prepend an '@' character to a settings file, e.g. @$(input[1])?\n" *
            "If not, note that only $ALLOWED_FILE_SUFFIXES_STRING file types are supported"
        else
            "No valid files were found for processing. Note that currently only $ALLOWED_FILE_SUFFIXES_STRING file types are supported"
        end
        error(msg)
    end

    # Get output folders
    outputfolders = if output === nothing
        dirname.(inputfiles)
    else
        [output for _ in 1:length(inputfiles)]
    end

    # Get mask files
    maskfiles = if isempty(mask)
        fill(nothing, length(inputfiles))
    elseif length(mask) == length(inputfiles)
        String.(mask)
    else
        error("Number of mask files passed does not equal the number of input image files passed")
    end

    # Create file_info dictionaries
    file_info = Dict{Symbol,Any}[]
    for (inputfile, outputfolder, maskfile) in zip(inputfiles, outputfolders, maskfiles)
        d = eltype(file_info)(
            :inputfile => inputfile,
            :outputfolder => outputfolder,
            :maskfile => maskfile,
            :choppedinputfile => chop_allowed_suffix(basename(inputfile)),
        )
        push!(file_info, d)
    end

    return file_info
end

function load_image(filename, ::Val{N}) where {N}
    image = if maybe_get_suffix(filename) == ".mat"
        # Load first `N`-dimensional array which is found, or throw an error if none are found
        data = MAT.matread(filename)
        array_keys = findall(x -> x isa AbstractArray{T,N} where {T}, data)
        if isempty(array_keys)
            error("No 4D array was found in the input file: $filename")
        end
        if length(array_keys) > 1
            array_keys = sort(array_keys)
            @warn "Multiple possible images found in file: $(filename)\nChoosing field $(repr(array_keys[1])) out of the following options: $(join(repr.(array_keys), ", "))"
        end
        data[array_keys[1]]

    elseif maybe_get_suffix(filename) ∈ (".nii", ".nii.gz")
        # Check slope field; if scl_slope == 0, data is not scaled and raw data should be returned
        #   See e.g. https://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/scl_slopeinter.html
        # Loaded data is coerced to a `N`-dimensional array
        data = NIfTI.niread(filename)
        if data.header.scl_slope == 0
            data.raw[ntuple(_ -> Colon(), N)...] # Return raw data
        else
            data[ntuple(_ -> Colon(), N)...] # Return scaled data (getindex from the NIfTI package handles scaling)
        end

    elseif maybe_get_suffix(filename) ∈ (".par", ".xml", ".rec")
        # Load PAR/REC or XML/REC file, coercing resulting data into a `N`-dimensional array
        _, data = ParXRec.parxrec(filename)
        data[ntuple(_ -> Colon(), N)...]

    else
        error("Currently, only $ALLOWED_FILE_SUFFIXES_STRING files are supported")
    end

    # Currently, the pipeline is ~twice as fast on Float64 arrays than Float32 arrays (unclear why).
    # However, the MATLAB toolbox converts images to double as well, so here we simply do the same
    image = convert(Array{Float64,N}, image)

    return image
end
load_image(filename; ndims::Int = 4) = load_image(filename, Val(ndims))

function try_apply_maskfile!(image, maskfile)
    try
        image .*= load_image(maskfile, Val(3))
    catch e
        @warn "Error while loading mask file: $maskfile"
        @warn sprint(showerror, e, catch_backtrace())
    end
    return image
end

function try_apply_bet!(image, betpath, betargs)
    try
        image .*= make_bet_mask(image, betpath, betargs)
    catch e
        @warn "Error while making mask using BET"
        @warn sprint(showerror, e, catch_backtrace())
    end
    return image
end

function make_bet_mask(image::Array{T,3}, betpath, betargs) where {T}
    # Split betargs, and ensure that "-m" (make binary mask) is among args
    args = convert(Vector{String}, filter!(!isempty, split(betargs, " ")))
    if "-m" ∉ args
        push!(args, "-m")
    end

    # Create mask using BET and return mask
    mask = mktempdir() do temppath
        tempbase = basename(tempname())
        nifti_imagefile = joinpath(temppath, tempbase * ".nii")
        nifti_maskfile = joinpath(temppath, tempbase * ".bet")
        NIfTI.niwrite(nifti_imagefile, NIfTI.NIVolume(image)) # create nifti file for bet
        run(Cmd([
            betpath;
            nifti_imagefile;
            nifti_maskfile;
            args
        ]))
        # BET appends "_mask" and ".nii.gz" to output file name.
        # Find this file, ensure it is unique, then load and return it
        bet_maskfiles = filter!(file -> startswith(file, tempbase * ".bet_mask"), readdir(temppath))
        @assert length(bet_maskfiles) == 1 # ensure unique; this should never be false using a temp filename
        load_image(joinpath(temppath, bet_maskfiles[1]), Val(3))
    end

    return mask
end
make_bet_mask(image::Array{T,4}, args...; kwargs...) where {T} = make_bet_mask(image[:,:,:,1], args...; kwargs...) # use first echo

maybe_get_first(f, xs) = findfirst(f, xs) |> I -> I === nothing ? nothing : xs[I]
maybe_get_suffix(filename) = maybe_get_first(ext -> endswith(lowercase(filename), ext), ALLOWED_FILE_SUFFIXES) # case-insensitive
is_allowed_suffix(filename) = maybe_get_suffix(filename) !== nothing

function chop_allowed_suffix(filename::AbstractString)
    ext = maybe_get_suffix(filename)
    if ext !== nothing
        return filename[1:end-length(ext)]
    else
        error("Currently only $ALLOWED_FILE_SUFFIXES_STRING file types are supported")
    end
end


text:
using OrnsteinZernike
using Documenter

DocMeta.setdocmeta!(
    OrnsteinZernike, :DocTestSetup, :(using OrnsteinZernike); recursive=true
)

makedocs(;
    modules=[OrnsteinZernike],
    authors="Edwin Bedolla",
    repo="https://github.com/edwinb-ai/OrnsteinZernike.jl/blob/{commit}{path}#{line}",
    sitename="OrnsteinZernike.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://edwinb-ai.github.io/OrnsteinZernike.jl",
        assets=String[],
    ),
    pages=["Home" => "index.md"],
)

deploydocs(; repo="github.com/edwinb-ai/OrnsteinZernike.jl", devbranch="main")


text:
"""
	trend1d(cmd0::String="", arg1=[], kwargs...)

Fit a [weighted] [robust] polynomial/Fourier model for y = f(x) to xy[w] data.

Full option list at [`trend1d`](http://gmt.soest.hawaii.edu/doc/latest/trend1d.html)

Parameters
----------

- **F** : **output** : -- Str --   Flags = xymrw|p|P|c

    Specify up to five letters from the set {x y m r w} in any order to create columns of output.
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/trend1d.html#f)
- **N** : **model_components** : -- Str --      Flags = [p|P|f|F|c|C|s|S|x]n[,…][+llength][+oorigin][+r]

    Specify Specify the number of terms in the model, n_model, and append +r to do a robust fit. E.g., a robust bilinear model is -N4+r.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/trend1d.html#n)

- **C** : **condition_number** : -- Number --   Flags = condition_number

    Set the maximum allowed condition number for the matrix solution.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/trend1d.html#c)
- **I** : **confidence_level** : -- Number or [] --   Flags = [confidence_level]

    Iteratively increase the number of model parameters, starting at one, until n_model is reachedx
    or the reduction in variance of the model is not significant at the confidence_level level.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/trend1d.html#i)
- **W** : **weights** : -- Str or [] --     Flags = [+s]

    Weights are supplied in input column 3. Do a weighted least squares fit [or start with
    these weights when doing the iterative robust fit].
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/trend1d.html#w)
- $(GMT.opt_V)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_swap_xy)
"""
function trend1d(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("trend1d", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd = parse_V_params("", d)
	cmd, = parse_b(cmd, d)
	cmd, = parse_d(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'C', d, [:C :condition_number])
	cmd = add_opt(cmd, 'I', d, [:I :confidence_level])
	cmd = add_opt(cmd, 'F', d, [:F :output])
	cmd = add_opt(cmd, 'N', d, [:N :model_components])
	cmd = add_opt(cmd, 'W', d, [:W :weights])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "trend1d", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
trend1d(arg1=[], cmd0::String=""; kw...) = trend1d(cmd0, arg1; kw...)


text:
module XIntsTests
using ReTest
include("tests.jl")
end


text:
###############################################
#FEniCS tutorial demo program: Incompressible Navier-Stokes equations
#for flow around a cylinder using the Incremental Pressure Correction
#Scheme (IPCS).
#  u' + u . nabla(u)) - div(σ(u, p)) = f
#                                 div(u) = 0
###############################################

module ft08

using FenicsPy

# 参数：ρ,μ,Δt
T = 5
num_steps = 5000
Δt = T / num_steps
μ = 0.001
ρ = 1

# Create mesh
channel = Rectangle(Point([0.0, 0.0]), Point([2.2, 0.41]))
cylinder = Circle(Point([0.2, 0.2]), 0.05)
Ω = channel - cylinder
mesh = generate_mesh(Ω, 64)

# Define function spaces
V = VectorFunctionSpace(mesh, "P", 2)
Q = FunctionSpace(mesh, "P", 1)

# Define boundaries
inflow   = "near(x[0], 0)"
outflow  = "near(x[0], 2.2)"
walls    = "near(x[1], 0) || near(x[1], 0.41)"
cylinder = "on_boundary && x[0]>0.1 && x[0]<0.3 && x[1]>0.1 && x[1]<0.3"

# Define inflow profile
inflow_profile = ("4.0*1.5*x[1]*(0.41 - x[1]) / pow(0.41, 2)", "0")

# Define boundary conditions
bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree=2), inflow)
bcu_walls = DirichletBC(V, Constant((0, 0)), walls)
bcu_cylinder = DirichletBC(V, Constant((0, 0)), cylinder)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_inflow,bcu_walls,bcu_cylinder]
bcp = [bcp_outflow]

u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = FeFunction(V)
u_ = FeFunction(V)
p_n = FeFunction(Q)
p_ = FeFunction(Q)

U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 0))

# Define symmetric gradient
function ϵ(u)
    return sym(nabla_grad(u))
end

# Define stress tensor
function σ(u, p)
    return 2*μ*ϵ(u) - p*Identity(len(u))
end

# Define variational problem for step 1
F1 = ρ*dot((u - u_n) / Δt, v)*dx  + ρ*dot(dot(u_n, nabla_grad(u_n)), v)*dx +
            inner(σ(U, p_n), ϵ(v))*dx  + dot(p_n*n, v)*ds -
            dot(2*μ*nabla_grad(U)*n, v)*ds   - ρ*dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (ρ/Δt)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - (Δt/ρ)*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files for visualization output
xdmffile_u = XDMFFile("navier_stokes_cylinder/velocity.xdmf")
xdmffile_p = XDMFFile("navier_stokes_cylinder/pressure.xdmf")

# Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries("navier_stokes_cylinder/velocity_series")
timeseries_p = TimeSeries("navier_stokes_cylinder/pressure_series")

# Save mesh to file (for use in reaction_system.py)
File("navier_stokes_cylinder/cylinder.xml.gz") << mesh

t = 0
for i in 1:num_steps

    # Update current time
    global t += Δt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1, "bicgstab", "hypre_amg")

    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2, "bicgstab", "hypre_amg")

    b3 = assemble(L3)
    solve(A3, u_.vector(), b3, "cg", "sor")

    #update values
    u_n.assign(u_)
    p_n.assign(p_)

    if (i-1)%25 == 0 || i == num_steps

        xdmffile_u.write(u_, t)
        xdmffile_p.write(p_, t)

        # Save nodal values to file
        timeseries_u.store(u_.vector(), t)
        timeseries_p.store(p_.vector(), t)

        println(i, "/", num_steps, " \t max:", max(array(u_)...))

    end

    # Plot solution
    #plot(u_)
    #plot(p_)
end

xdmffile_u.close()
xdmffile_p.close()

end # module ft08


text:
using Test
using CSV
using DataFrames
using StatsBase
using BenchmarkTools
using PrettyTables
using Observables
using ScottishTaxBenefitModel
using ScottishTaxBenefitModel.GeneralTaxComponents
using ScottishTaxBenefitModel.STBParameters
using ScottishTaxBenefitModel.Runner: do_one_run
using ScottishTaxBenefitModel.RunSettings
using .Utils
using .Monitor: Progress
using .ExampleHelpers
using .STBOutput: make_poverty_line, summarise_inc_frame,
    dump_frames, summarise_frames, make_gain_lose

settings = Settings()

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 120
BenchmarkTools.DEFAULT_PARAMETERS.samples = 2

tot = 0

# observer = Observer(Progress("",0,0,0))
obs = Observable( Progress(settings.uuid,"",0,0,0,0))
of = on(obs) do p
    global tot
    println(p)

    tot += p.step
    println(tot)
end

function basic_run( ; print_test :: Bool, mtrouting :: MT_Routing )
    settings.means_tested_routing = mtrouting
    settings.run_name="run-$(mtrouting)-$(date_string())"
    sys = [get_system(scotland=false), get_system( scotland=true )]
    tot = 0
    results = do_one_run( settings, sys, obs )
    h1 = results.hh[1]
    pretty_table( h1[:,[:weighted_people,:bhc_net_income,:eq_bhc_net_income,:ahc_net_income,:eq_ahc_net_income]] )
    settings.poverty_line = make_poverty_line( results.hh[1], settings )
    dump_frames( settings, results )
    println( "poverty line = $(settings.poverty_line)")
    outf = summarise_frames( results, settings )
    println( outf )
    gl = make_gain_lose( results.hh[1], results.hh[2], settings )
    println(gl)
end

@testset "MR test" begin
    @time begin
        settings.means_tested_routing = modelled_phase_in
        settings.do_marginal_rates = true
        settings.dump_frames = true
        sys = [get_system(scotland=false), get_system( scotland=true )]
        results = do_one_run( settings, sys, obs )
        settings.poverty_line = make_poverty_line( results.hh[1], settings )
        outf = summarise_frames( results, settings )
        println( outf.metrs[1] )
    end
end

@testset "basic run timing" begin
    for mt in instances( MT_Routing )
        println( "starting run using $mt routing")
        @time basic_run( print_test=true, mtrouting = mt )
    end
    # @benchmark frames =
    # print(t)
end

@testset "thread test" begin
    # 42.531500 53 secs
    settings.requested_threads = 4
    @time basic_run( print_test=true, mtrouting = modelled_phase_in )

    settings.do_marginal_rates = false
    @time basic_run( print_test=true, mtrouting = modelled_phase_in )
end

#=
if print_test
    summary_output = summarise_results!( results=results, base_results=base_results )
    print( "   deciles = $( summary_output.deciles)\n\n" )
    print( "   poverty_line = $(summary_output.poverty_line)\n\n" )
    print( "   inequality = $(summary_output.inequality)\n\n" )
    print( "   poverty = $(summary_output.poverty)\n\n" )
    print( "   gainlose_by_sex = $(summary_output.gainlose_by_sex)\n\n" )
    print( "   gainlose_by_thing = $(summary_output.gainlose_by_thing)\n\n" )
    print( "   metr_histogram= $(summary_output.metr_histogram)\n\n")
    println( "SUMMARY OUTPUT")
    println( summary_output )
    println( "as JSON")
    println( JSON.json( summary_output ))
end

=#

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
using Documenter, OrdinaryDiffEq

makedocs(
    sitename="OrdinaryDiffEq.jl",
    authors="Chris Rackauckas et al.",
    clean=true,
    doctest=false,
    modules=[OrdinaryDiffEq],

    format=Documenter.HTML(assets=["assets/favicon.ico"],
                           canonical="https://ordinarydiffeq.sciml.ai/stable/"),

    pages=[
        "OrdinaryDiffEq.jl: ODE solvers and utilities" => "index.md",
        "Usage" => "usage.md"
    ]
)

deploydocs(
    repo="github.com/SciML/OrdinaryDiffEq.jl";
    push_preview=true
)


text:
using MyExample
using Test

using LinearAlgebra: Eigen
import LinearAlgebra: eigen, \, det, logdet, inv

#test

@testset "MyExample.jl" begin
    # Write your own tests here.

    @test NSKRRegressor(randn(4,4), [2,2], "zeros").λ_ == [2,2]

    @test NSKRRegressor([2,2], "zeros").λ_ == [2,2]

    @test length(Mps(randn(2,2), randn(4,4)).matrices) == 2

    model = NSKRRegressor([2,2], "zeros")
    fit(model, eigen.([randn(2,2), randn(5,5)]), randn(2,5), "traininds")

end

max_stars_repo_name: padwulf/MyExample

max_issues_repo_name: padwulf/MyExample

max_forks_repo_name: padwulf/MyExample


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestConstraintSOCtoPSD

using Test

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_SOCtoPSD()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.SOCtoPSD{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_VectorOfVariabless_SecondOrderCone",
            "test_basic_VectorAffineFunction_SecondOrderCone",
            "test_basic_VectorQuadraticFunction_SecondOrderCone",
        ],
    )
    MOI.empty!(bridged_mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1.0, 1 / √2, 1 / √2],
            (
                MOI.VectorAffineFunction{Float64},
                MOI.PositiveSemidefiniteConeTriangle,
            ) => [[√2 / 2, -1 / 2, √2 / 4, -1 / 2, √2 / 4, √2 / 4]],
            (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-√2]],
        )
    MOI.Test.test_conic_SecondOrderCone_VectorOfVariables(bridged_mock, config)
    MOI.empty!(bridged_mock)
    MOI.Test.test_conic_SecondOrderCone_VectorAffineFunction(
        bridged_mock,
        config,
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.SecondOrderCone,
            }(),
        ),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = [√2, -1.0, -1.0]
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        3,
        ((
            MOI.VectorAffineFunction{Float64},
            MOI.PositiveSemidefiniteConeTriangle,
            0,
        ),),
    )
    return
end

function test_RSOCtoPSD()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.RSOCtoPSD{Float64}(mock)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_VectorOfVariabless_RotatedSecondOrderCone",
            "test_basic_VectorAffineFunction_RotatedSecondOrderCone",
            "test_basic_VectorQuadraticFunction_RotatedSecondOrderCone",
        ],
    )
    MOI.empty!(bridged_mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [0.5, 1.0, 1 / √2, 1 / √2],
            (MOI.VariableIndex, MOI.EqualTo{Float64}) => [-√2, -1 / √2],
            (
                MOI.VectorAffineFunction{Float64},
                MOI.PositiveSemidefiniteConeTriangle,
            ) => [[√2, -1 / 2, √2 / 8, -1 / 2, √2 / 8, √2 / 8]],
        )
    MOI.Test.test_conic_RotatedSecondOrderCone_VectorOfVariables(
        bridged_mock,
        config,
    )
    MOI.empty!(bridged_mock)
    mock.optimize! =
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [1 / √2, 1 / √2],
            (
                MOI.VectorAffineFunction{Float64},
                MOI.PositiveSemidefiniteConeTriangle,
            ) => [[√2, -1 / 2, √2 / 8, -1 / 2, √2 / 8, √2 / 8]],
        )
    MOI.Test.test_conic_RotatedSecondOrderCone_VectorAffineFunction(
        bridged_mock,
        config,
    )
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.RotatedSecondOrderCone,
            }(),
        ),
    )
    for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = [√2, √2 / 4, -1.0, -1.0]
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        ((
            MOI.VectorAffineFunction{Float64},
            MOI.PositiveSemidefiniteConeTriangle,
            0,
        ),),
    )
    return
end

# TODO(odow): we should fix this so it does not error.
function test_rsoc_to_psd_dimension_2()
    inner = MOI.Utilities.Model{Float64}()
    model = MOI.Bridges.Constraint.RSOCtoPSD{Float64}(inner)
    x = MOI.add_variables(model, 2)
    @test_throws(
        ErrorException,
        MOI.add_constraint(
            model,
            MOI.VectorOfVariables(x),
            MOI.RotatedSecondOrderCone(2),
        ),
    )
    return
end

end  # module

TestConstraintSOCtoPSD.runtests()

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
using StatGeochemBase
using Documenter

DocMeta.setdocmeta!(StatGeochemBase, :DocTestSetup, :(using StatGeochemBase); recursive=true)

makedocs(;
    modules=[StatGeochemBase],
    authors="C. Brenhin Keller",
    repo="https://github.com/brenhinkeller/StatGeochemBase.jl/blob/{commit}{path}#{line}",
    sitename="StatGeochemBase.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://brenhinkeller.github.io/StatGeochemBase.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/brenhinkeller/StatGeochemBase.jl",
    devbranch = "main",
)

max_stars_repo_name: brenhinkeller/StatGeochemBase

max_issues_repo_name: brenhinkeller/StatGeochemBase

max_forks_repo_name: brenhinkeller/StatGeochemBase


text:
# test if the PooldDataArray is constructed the way I think it is
# f.pool
ft = dtmtable.features[2]
ft.pool
unique(ft)

for j = 1:size(dtmtable.features, 2)
    ft = dtmtable.features[j]
    for i = 1:length(ft.pool)
        @assert all(ft[ft.refs .== UInt8(i)] .== ft.pool[i])
    end
end

max_stars_repo_path: testing/pda order critical property.jl

max_issues_repo_path: testing/pda order critical property.jl

max_forks_repo_path: testing/pda order critical property.jl


text:
using Base.Test
using LibPNG
using ColorTypes
using FixedPointNumbers
using ImageCore
using MicroLogging

configure_logging(min_level=MicroLogging.Debug)

include("images.jl")

tmpdir = joinpath(tempdir(), "LibPNG")
try
    rm(tmpdir, recursive = true)
end


text:
# This file extends the FileIO.jl interface for filenames ending with ".mid".
# See https://juliaio.github.io/FileIO.jl/stable/implementing/#Implementing-loaders/savers

using FileIO
export load, save

function fileio_load(f::File{format"MIDI"})
    open(f) do s
        skipmagic(s)
        midifile = load(s)
    end
end

function fileio_load(s::Stream{format"MIDI"})
    midifile = MIDIFile()

    # Skip the next four bytes - this is the header size, and it's always equal to 6.
    skip(s, 4)

    # Read the format code. 0 = single track, 1 = multiple tracks, 2 = multiple songs
    # Remember - MIDI files store data in big-endian format, which is why ntoh is used
    midifile.format = ntoh(read(s, UInt16))

    # Get the number of tracks and time division
    numberoftracks = ntoh(read(s, UInt16))
    midifile.tpq = ntoh(read(s, Int16))
    midifile.tracks = [readtrack(s.io) for x in 1:numberoftracks]

    return midifile
end

function fileio_save(f::File{format"MIDI"}, data::MIDIFile)
    open(f, "w") do s
        fileio_save(s, data)
    end
end

function fileio_save(s::Stream{format"MIDI"}, data::MIDIFile)
    write(s, magic(format"MIDI"))

    write(s, hton(convert(UInt32, 6))) # Header length
    write(s, hton(data.format))
    write(s, hton(convert(UInt16, length(data.tracks))))
    write(s, hton(data.tpq))

    map(track -> writetrack(s.io, track), data.tracks)
    return data
end

function fileio_save(s::Union{File{format"MIDI"}, Stream{format"MIDI"}}, notes::Notes)
    track = MIDITrack()
    addnotes!(track, notes)
    midi = MIDIFile(1, notes.tpq, [track])
    fileio_save(s, midi)
    return midi
end