text:
#### Specialized matrix types ####

## Hermitian tridiagonal matrices
type SymTridiagonal{T} <: AbstractMatrix{T}
    dv::Vector{T}                        # diagonal
    ev::Vector{T}                        # subdiagonal
    function SymTridiagonal(dv::Vector{T}, ev::Vector{T})
        length(ev)==length(dv)-1 || throw(DimensionMismatch(""))
        new(dv,ev)
    end
end

SymTridiagonal{T}(dv::Vector{T}, ev::Vector{T}) = SymTridiagonal{T}(copy(dv), copy(ev))

function SymTridiagonal{Td,Te}(dv::Vector{Td}, ev::Vector{Te})
    T = promote_type(Td,Te)
    SymTridiagonal(convert(Vector{T}, dv), convert(Vector{T}, ev))
end

SymTridiagonal(A::AbstractMatrix) = diag(A,1)==diag(A,-1)?SymTridiagonal(diag(A), diag(A,1)):throw(DimensionMismatch("matrix is not symmetric; cannot convert to SymTridiagonal"))
full{T}(M::SymTridiagonal{T}) = convert(Matrix{T}, M)
convert{T}(::Type{Matrix{T}}, M::SymTridiagonal{T})=diagm(M.dv)+diagm(M.ev,-1)+conj(diagm(M.ev,1))

size(m::SymTridiagonal) = (length(m.dv), length(m.dv))
size(m::SymTridiagonal, d::Integer) = d<1 ? error("dimension out of range") : (d<=2 ? length(m.dv) : 1)

#Elementary operations
for func in (:copy, :round, :iround, :conj)
    @eval begin
        ($func)(M::SymTridiagonal) = SymTridiagonal(($func)(M.dv), ($func)(M.ev))
    end
end

transpose(M::SymTridiagonal) = M #Identity operation
ctranspose(M::SymTridiagonal) = conj(M)

function diag{T}(M::SymTridiagonal{T}, n::Integer=0)
    absn = abs(n)
    absn==0 ? M.dv : absn==1 ? M.ev : absn<size(M,1) ? zeros(T,size(M,1)-absn) : throw(BoundsError())
end

+(A::SymTridiagonal, B::SymTridiagonal) = SymTridiagonal(A.dv+B.dv, A.ev+B.ev)
-(A::SymTridiagonal, B::SymTridiagonal) = SymTridiagonal(A.dv-B.dv, A.ev-B.ev)
*(A::SymTridiagonal, B::Number) = SymTridiagonal(A.dv*B, A.ev*B)
*(B::Number, A::SymTridiagonal) = A*B
/(A::SymTridiagonal, B::Number) = SymTridiagonal(A.dv/B, A.ev/B)
==(A::SymTridiagonal, B::SymTridiagonal) = (A.dv==B.dv) && (A.ev==B.ev)

## Solver
function \{T<:BlasFloat}(M::SymTridiagonal{T}, rhs::StridedVecOrMat{T})
    if stride(rhs, 1) == 1
        return LAPACK.gtsv!(copy(M.ev), copy(M.dv), copy(M.ev), copy(rhs))
    end
    solve(Tridiagonal(M), rhs)  # use the Julia "fallback"
end

#Wrap LAPACK DSTE{GR,BZ} to compute eigenvalues
eig{T<:BlasFloat}(m::SymTridiagonal{T}) = LAPACK.stegr!('V', copy(m.dv), copy(m.ev))
eigvals{T<:BlasFloat}(m::SymTridiagonal{T}, il::Int, iu::Int) = LAPACK.stegr!('N', 'I', copy(m.dv), copy(m.ev), 0.0, 0.0, il, iu)[1]
eigvals{T<:BlasFloat}(m::SymTridiagonal{T}, vl::Real, vu::Real) = LAPACK.stegr!('N', 'V', copy(m.dv), copy(m.ev), vl, vu, 0, 0)[1]
eigvals{T<:BlasFloat}(m::SymTridiagonal{T}) = LAPACK.stev!('N', copy(m.dv), copy(m.ev))[1]

#Computes largest and smallest eigenvalue
eigmax(m::SymTridiagonal) = eigvals(m, size(m, 1), size(m, 1))[1]
eigmin(m::SymTridiagonal) = eigvals(m, 1, 1)[1]

#Compute selected eigenvectors only corresponding to particular eigenvalues
eigvecs(m::SymTridiagonal) = eig(m)[2]
eigvecs{T<:BlasFloat,Eigenvalue<:Real}(m::SymTridiagonal{T}, eigvals::Vector{Eigenvalue}) = LAPACK.stein!(m.dv, m.ev, eigvals)

###################
# Generic methods #
###################

#Needed for inv_usmani()
type ZeroOffsetVector
    data::Vector
end
getindex (a::ZeroOffsetVector, i) = a.data[i+1]
setindex!(a::ZeroOffsetVector, x, i) = a.data[i+1]=x

#Implements the inverse using the recurrence relation between principal minors
# a, b, c are assumed to be the subdiagonal, diagonal, and superdiagonal of
# a tridiagonal matrix.
#Ref:
#    R. Usmani, "Inversion of a tridiagonal Jacobi matrix",
#    Linear Algebra and its Applications 212-213 (1994), pp.413-414
#    doi:10.1016/0024-3795(94)90414-6
function inv_usmani{T}(a::Vector{T}, b::Vector{T}, c::Vector{T})
    n = length(b)
    θ = ZeroOffsetVector(zeros(T, n+1)) #principal minors of A
    θ[0] = 1
    n>=1 && (θ[1] = b[1])
    for i=2:n
        θ[i] = b[i]*θ[i-1]-a[i-1]*c[i-1]*θ[i-2]
    end
    φ = zeros(T, n+1)
    φ[n+1] = 1
    n>=1 && (φ[n] = b[n])
    for i=n-1:-1:1
        φ[i] = b[i]*φ[i+1]-a[i]*c[i]*φ[i+2]
    end
    α = Array(T, n, n)
    for i=1:n, j=1:n
        sign = (i+j)%2==0 ? (+) : (-)
        if i<j
            α[i,j]=(sign)(prod(c[i:j-1]))*θ[i-1]*φ[j+1]/θ[n]
        elseif i==j
            α[i,i]=                       θ[i-1]*φ[i+1]/θ[n]
        else #i>j
            α[i,j]=(sign)(prod(a[j:i-1]))*θ[j-1]*φ[i+1]/θ[n]
        end
    end
    α
end

#Implements the determinant using principal minors
#Inputs and reference are as above for inv_usmani()
function det_usmani{T}(a::Vector{T}, b::Vector{T}, c::Vector{T})
    n = length(b)
    θa = one(T)
    n==0 && return θa
    θb = b[1]
    for i=2:n
        θb, θa = b[i]*θb-a[i-1]*c[i-1]*θa, θb
    end
    return θb
end

inv(A::SymTridiagonal) = inv_usmani(A.ev, A.dv, A.ev)
det(A::SymTridiagonal) = det_usmani(A.ev, A.dv, A.ev)

function getindex{T}(A::SymTridiagonal{T}, i::Integer, j::Integer)
    (1<=i<=size(A,2) && 1<=j<=size(A,2)) || throw(BoundsError())
    i==j ? A.dv[i] : i==j+1 ? A.ev[j] : i+1==j ? A.ev[i] : zero(T)
end

## Tridiagonal matrices ##
type Tridiagonal{T} <: AbstractMatrix{T}
    dl::Vector{T}    # sub-diagonal
    d::Vector{T}     # diagonal
    du::Vector{T}    # sup-diagonal
    dutmp::Vector{T} # scratch space for vector RHS solver, sup-diagonal
    rhstmp::Vector{T}# scratch space, rhs

    function Tridiagonal(N::Integer)
        dutmp = Array(T, N-1)
        rhstmp = Array(T, N)
        new(dutmp, rhstmp, similar(dutmp), similar(dutmp), similar(rhstmp))
    end

    function Tridiagonal(dl::Vector{T}, d::Vector{T}, du::Vector{T})
        N = length(d)
        if (length(dl) != N-1 || length(du) != N-1)
            error(string("Cannot make Tridiagonal from incompatible lengths of subdiagonal, diagonal and superdiagonal: (", length(dl), ", ", length(d), ", ", length(du),")"))
        end
        new(copy(dl), copy(d), copy(du), Array(T,N-1), Array(T,N))
    end
end

Tridiagonal{T}(dl::Vector{T}, d::Vector{T}, du::Vector{T}) = Tridiagonal{T}(dl, d, du)

function Tridiagonal{Tl, Td, Tu}(dl::Vector{Tl}, d::Vector{Td}, du::Vector{Tu})
    Tridiagonal(map(v->copy(convert(Vector{promote_type(Tl,Td,Tu)}, v)), (dl, d, du))...)
end

size(M::Tridiagonal) = (length(M.d), length(M.d))
size(M::Tridiagonal, d::Integer) = d<1 ? error("dimension out of range") : (d<=2 ? length(M.d) : 1)

full{T}(M::Tridiagonal{T}) = convert(Matrix{T}, M)
function convert{T}(::Type{Matrix{T}}, M::Tridiagonal{T})
    A = zeros(T, size(M))
    for i = 1:length(M.d)
        A[i,i] = M.d[i]
    end
    for i = 1:length(M.d)-1
        A[i+1,i] = M.dl[i]
        A[i,i+1] = M.du[i]
    end
    A
end
function similar(M::Tridiagonal, T, dims::Dims)
    if length(dims) != 2 || dims[1] != dims[2]
        throw(DimensionMismatch("Tridiagonal matrices must be square"))
    end
    Tridiagonal{T}(dims[1])
end

# Operations on Tridiagonal matrices
copy!(dest::Tridiagonal, src::Tridiagonal) = Tridiagonal(copy!(dest.dl, src.dl), copy!(dest.d, src.d), copy!(dest.du, src.du))

#Elementary operations
for func in (:copy, :round, :iround, :conj)
    @eval begin
        ($func)(M::Tridiagonal) = Tridiagonal(map(($func), (M.dl, M.d, M.du))...)
    end
end

transpose(M::Tridiagonal) = Tridiagonal(M.du, M.d, M.dl)
ctranspose(M::Tridiagonal) = conj(transpose(M))

diag{T}(M::Tridiagonal{T}, n::Integer=0) = n==0 ? M.d : n==-1 ? M.dl : n==1 ? M.du : abs(n)<size(M,1) ? zeros(T,size(M,1)-abs(n)) : throw(BoundsError())

function getindex{T}(A::Tridiagonal{T}, i::Integer, j::Integer)
    (1<=i<=size(A,2) && 1<=j<=size(A,2)) || throw(BoundsError())
    i==j ? A.d[i] : i==j+1 ? A.dl[j] : i+1==j ? A.du[i] : zero(T)
end

###################
# Generic methods #
###################

+(A::Tridiagonal, B::Tridiagonal) = Tridiagonal(A.dl+B.dl, A.d+B.d, A.du+B.du)
-(A::Tridiagonal, B::Tridiagonal) = Tridiagonal(A.dl-B.dl, A.d-B.d, A.du+B.du)
*(A::Tridiagonal, B::Number) = Tridiagonal(A.dl*B, A.d*B, A.du*B)
*(B::Number, A::Tridiagonal) = A*B
/(A::Tridiagonal, B::Number) = Tridiagonal(A.dl/B, A.d/B, A.du/B)

==(A::Tridiagonal, B::Tridiagonal) = (A.dl==B.dl) && (A.d==B.d) && (A.du==B.du)
==(A::Tridiagonal, B::SymTridiagonal) = (A.dl==A.du==B.ev) && (A.d==B.dv)
==(A::SymTridiagonal, B::Tridiagonal) = (B.dl==B.du==A.ev) && (B.d==A.dv)

inv(A::Tridiagonal) = inv_usmani(A.dl, A.d, A.du)
det(A::Tridiagonal) = det_usmani(A.dl, A.d, A.du)

# Elementary operations that mix Tridiagonal and SymTridiagonal matrices
convert(::Type{Tridiagonal}, A::SymTridiagonal) = Tridiagonal(A.ev, A.dv, A.ev)
+(A::Tridiagonal, B::SymTridiagonal) = Tridiagonal(A.dl+B.ev, A.d+B.dv, A.du+B.ev)
+(A::SymTridiagonal, B::Tridiagonal) = Tridiagonal(A.ev+B.dl, A.dv+B.d, A.ev+B.du)
-(A::Tridiagonal, B::SymTridiagonal) = Tridiagonal(A.dl-B.ev, A.d-B.dv, A.du-B.ev)
-(A::SymTridiagonal, B::Tridiagonal) = Tridiagonal(A.ev-B.dl, A.dv-B.d, A.ev-B.du)

convert{T}(::Type{Tridiagonal{T}},M::Tridiagonal) = Tridiagonal(convert(Vector{T}, M.dl), convert(Vector{T}, M.d), convert(Vector{T}, M.du))
convert{T}(::Type{AbstractMatrix{T}},M::Tridiagonal) = Tridiagonal(convert(Vector{T}, M.dl), convert(Vector{T}, M.d), convert(Vector{T}, M.du))
convert{T}(::Type{Tridiagonal{T}}, M::SymTridiagonal{T}) = Tridiagonal(M)
convert{T}(::Type{SymTridiagonal{T}}, M::Tridiagonal) = M.dl==M.du ? (SymTridiagonal(M.dl, M.d)) :
    error("Tridiagonal is not symmetric, cannot convert to SymTridiagonal")
convert{T}(::Type{SymTridiagonal{T}},M::SymTridiagonal) = SymTridiagonal(convert(Vector{T}, M.dv), convert(Vector{T}, M.ev))

## Solvers

#### Tridiagonal matrix routines ####
function \{T<:BlasFloat}(M::Tridiagonal{T}, rhs::StridedVecOrMat{T})
    if stride(rhs, 1) == 1
        return LAPACK.gtsv!(copy(M.dl), copy(M.d), copy(M.du), copy(rhs))
    end
    solve(M, rhs)  # use the Julia "fallback"
end

# This is definitely not going to work
#eig(M::Tridiagonal) = LAPACK.stev!('V', copy(M))

# Allocation-free variants
# Note that solve is non-aliasing, so you can use the same array for
# input and output
function solve!{T<:BlasFloat}(x::AbstractArray{T}, xrng::Range{Int}, M::Tridiagonal{T}, rhs::AbstractArray{T}, rhsrng::Range{Int})
    d = M.d
    N = length(d)
    if length(xrng) != N || length(rhsrng) != N
        throw(DimensionMismatch(""))
    end
    dl = M.dl
    du = M.du
    dutmp = M.dutmp
    rhstmp = M.rhstmp
    xstart = first(xrng)
    xstride = step(xrng)
    rhsstart = first(rhsrng)
    rhsstride = step(rhsrng)
    # Forward sweep
    denom = d[1]
    dulast = du[1] / denom
    dutmp[1] = dulast
    rhslast = rhs[rhsstart] / denom
    rhstmp[1] = rhslast
    irhs = rhsstart+rhsstride
    for i in 2:N-1
        dltmp = dl[i-1]
        denom = d[i] - dltmp*dulast
        dulast = du[i] / denom
        dutmp[i] = dulast
        rhslast = (rhs[irhs] - dltmp*rhslast)/denom
        rhstmp[i] = rhslast
        irhs += rhsstride
    end
    dltmp = dl[N-1]
    denom = d[N] - dltmp*dulast
    xlast = (rhs[irhs] - dltmp*rhslast)/denom
    # Backward sweep
    ix = xstart + (N-2)*xstride
    x[ix+xstride] = xlast
    for i in N-1:-1:1
        xlast = rhstmp[i] - dutmp[i]*xlast
        x[ix] = xlast
        ix -= xstride
    end
    nothing
end

function solve!(x::StridedVector, M::Tridiagonal, rhs::StridedVector)
    solve!(x, 1:length(x), M, rhs, 1:length(rhs))
    x
end
solve{TM<:BlasFloat,TB<:BlasFloat}(M::Tridiagonal{TM}, B::StridedVecOrMat{TB}) = solve!(zeros(typeof(one(TM)/one(TB)), size(B)), M, B)
solve(M::Tridiagonal, B::StridedVecOrMat) = solve(float(M), float(B))
function solve!(X::StridedMatrix, M::Tridiagonal, B::StridedMatrix)
    size(B, 1) == size(M, 1) || throw(DimensionMismatch(""))
    size(X) == size(B) || throw(DimensionMismatch(""))
    m, n = size(B)
    for j = 1:n
        r = ((j-1)*m+1):((j-1)*m+m)
        solve!(X, r, M, B, r)
    end
    X
end

# User-friendly solver
\(M::Tridiagonal, rhs::StridedVecOrMat) = solve(M, rhs)

# Tridiagonal multiplication
function mult(x::AbstractArray, xrng::Range{Int}, M::Tridiagonal, v::AbstractArray, vrng::Range{Int})
    dl = M.dl
    d = M.d
    du = M.du
    N = length(d)
    xi = first(xrng)
    xstride = step(xrng)
    vi = first(vrng)
    vstride = step(vrng)
    x[xi] = d[1]*v[vi] + du[1]*v[vi+vstride]
    xi += xstride
    for i = 2:N-1
        x[xi] = dl[i-1]*v[vi] + d[i]*v[vi+vstride] + du[i]*v[vi+2*vstride]
        xi += xstride
        vi += vstride
    end
    x[xi] = dl[N-1]*v[vi] + d[N]*v[vi+vstride]
    x
end

mult(x::StridedVector, M::Tridiagonal, v::StridedVector) = mult(x, 1:length(x), M, v, 1:length(v))

function mult(X::StridedMatrix, M::Tridiagonal, B::StridedMatrix)
    size(B, 1) == size(M, 1) || throw(DimensionMismatch(""))
    size(X) == size(B) || throw(DimensionMismatch(""))
    m, n = size(B)
    for j = 1:n
        r = ((j-1)*m+1):((j-1)*m+m)
        mult(X, r, M, B, r)
    end
    X
end

mult(X::StridedMatrix, M1::Tridiagonal, M2::Tridiagonal) = mult(X, M1, full(M2))

*(M::Tridiagonal, B::Union(StridedVector,StridedMatrix)) = mult(similar(B), M, B)
*(A::Tridiagonal, B::Tridiagonal) = A*full(B)

#### Factorizations for Tridiagonal ####
type LDLTTridiagonal{T<:BlasFloat,S<:BlasFloat} <: Factorization{T}
    D::Vector{S}
    E::Vector{T}
    function LDLTTridiagonal(D::Vector{S}, E::Vector{T})
        typeof(real(E[1])) == eltype(D) ? new(D, E) : error("element types do not match")
        new(D, E)
    end
end

LDLTTridiagonal{S<:BlasFloat,T<:BlasFloat}(D::Vector{S}, E::Vector{T}) = LDLTTridiagonal{T,S}(D, E)

ldltd!{T<:BlasFloat}(A::SymTridiagonal{T}) = LDLTTridiagonal(LAPACK.pttrf!(real(A.dv),A.ev)...)
ldltd{T<:BlasFloat}(A::SymTridiagonal{T}) = ldltd!(copy(A))
ldltd{T}(A::SymTridiagonal{T}) = (S = promote_type(typeof(sqrt(one(T))),Float32); S != T ? ldltd!(convert(SymTridiagonal{S},A)) : ldltd!(copy(A)))
factorize(A::SymTridiagonal) = ldltd(A)

A_ldiv_B!{T<:BlasReal}(C::LDLTTridiagonal{T}, B::StridedVecOrMat{T}) = LAPACK.pttrs!(C.D, C.E, B)
A_ldiv_B!{T<:BlasComplex}(C::LDLTTridiagonal{T}, B::StridedVecOrMat{T}) = LAPACK.pttrs!('L', C.D, C.E, B)
A_ldiv_B!(C::LDLTTridiagonal, B::StridedVecOrMat) = A_ldiv_B!(C, float(B))

type LUTridiagonal{T} <: Factorization{T}
    dl::Vector{T}
    d::Vector{T}
    du::Vector{T}
    du2::Vector{T}
    ipiv::Vector{BlasInt}
    # function LUTridiagonal(dl::Vector{T}, d::Vector{T}, du::Vector{T},
    #                        du2::Vector{T}, ipiv::Vector{BlasInt})
    #     n = length(d)
    #     if length(dl) != n - 1 || length(du) != n - 1 || length(ipiv) != n || length(du2) != n-2
    #         throw(DimensionMismatch("LUTridiagonal")
    #     end
    #     new(dl, d, du, du2, ipiv)
    # end
end
lufact!{T<:BlasFloat}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(A.dl,A.d,A.du)...)
lufact{T<:BlasFloat}(A::Tridiagonal{T}) = lufact!(copy(A))
lufact{T}(A::Tridiagonal{T}) = (S = promote_type(typeof(sqrt(one(T))),Float32); S != T ? lufact!(convert(Tridiagonal{S},A)) : lufact!(copy(A)))
factorize(A::Tridiagonal) = lufact(A)
#show(io, lu::LUTridiagonal) = print(io, "LU decomposition of ", summary(lu.lu))

function det{T}(lu::LUTridiagonal{T})
    n = length(lu.d)
    prod(lu.d) * (bool(sum(lu.ipiv .!= 1:n) % 2) ? -one(T) : one(T))
end

det{T<:BlasFloat}(A::Tridiagonal{T}) = det(lufact(A))

A_ldiv_B!{T<:BlasFloat}(lu::LUTridiagonal{T}, B::StridedVecOrMat{T}) =
    LAPACK.gttrs!('N', lu.dl, lu.d, lu.du, lu.du2, lu.ipiv, B)
A_ldiv_B!(lu::LUTridiagonal, B::StridedVecOrMat) = A_ldiv_B!(lu, float(B))

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
using DataStructures
using Pkg
using Plots

using GeometricProblems.LotkaVolterra2d
using GeometricProblems.Diagnostics

# select Plots backend
gr()
# pgfplots()
# plotly()
# plotlyjs()

# create output directory
if !isdir(out_dir * "/" * fig_dir)
    mkpath(out_dir * "/" * fig_dir)
end

# collect simulation parameters
function get_params(Δt, ntime, problem)
    params = OrderedDict()
    params["Δt"] = Δt
    params["nt"] = ntime
    params["problem"] = problem
    params["GeometricIntegrators.jl"] = Pkg.installed()["GeometricIntegrators"]
    params["GeometricProblems.jl"] = Pkg.installed()["GeometricProblems"]
    params
end

function run_lotka_volterra_2d(int::Integrator, ode::Equation, Δt, nt, filename; kwargs...)
    sim = Simulation(ode, int, Δt, "Lotka-Volterra 2d", prefix * filename * ".hdf5", nt; kwargs...)
    sol = run!(sim)
    run_diagnostics(sol, Δt, nt, filename)
end

function run_lotka_volterra_2d(tableau::AbstractTableau, ode::Equation, Δt, nt; integrator=Integrator, filename=nothing)
    int = integrator(ode, tableau, Δt)

    println()
    println("Running ", tableau.name, " (", filename, ") for ", nt, " time steps...")

    run_lotka_volterra_2d(int, ode, Δt, nt, filename)
end

function run_lotka_volterra_2d(basis::Basis, quadrature::Quadrature, ode::Equation, Δt, nt; integrator=Integrator, filename=nothing)
    int = IntegratorCGVI(ode, basis, quadrature, Δt)

    println()
    println("Running ", filename, " for ", nt, " time steps of Δt = ", Δt, "...")

    run_lotka_volterra_2d(int, ode, Δt, nt, filename)
end

function run_diagnostics(sol, Δt, nt, filename)
    if filename ≠ nothing
        ndrift = div(nt, 10)

        if nt > 100000
            npl = nplot
        else
            npl = 1
        end

        if nt ≤ 200
            markersize=5
        else
            markersize=.5
        end

        # plot overview
        plotlotkavolterra2d(sol)
        savefig(out_dir * "/" * fig_dir * "/" * prefix * filename * fig_suff)

        # plot solution
        plotlotkavolterra2dsolution(sol)
        # markersize=markersize
        # xlim(XLIM)
        # ylim(YLIM)
        # xlabel(L"$q_{1} (t)$", labelpad=10, fontsize=20)
        # ylabel(L"$q_{2} (t)$", labelpad=10, fontsize=20)
        savefig(out_dir * "/" * fig_dir * "/" * prefix * filename * "_solution" * fig_suff)

        # plot energy error
        H, ΔH = compute_energy_error(sol.t, sol.q);
        plotenergyerror(sol.t, ΔH)
        savefig(out_dir * "/" * fig_dir * "/" * prefix * filename * "_energy_error" * fig_suff)

        # plot energy drift
        Tdrift, Hdrift = compute_error_drift(sol.t, ΔH, ndrift)
        plotenergydrift(Tdrift, Hdrift)
        savefig(out_dir * "/" * fig_dir * "/" * prefix * filename * "_energy_drift" * fig_suff)

        # plot momentum error
        if isdefined(sol, :p)
            Δp = compute_momentum_error(sol.t, sol.q, sol.p)
            plotmomentumerror(sol.t, Δp)
            savefig(out_dir * "/" * fig_dir * "/" * prefix * filename * "_momentum_error" * fig_suff)
        end

        # plot Lagrange multiplier
        if isdefined(sol, :λ)
            plotlagrangemultiplier(sol.t, sol.λ)
            savefig(out_dir * "/" * fig_dir * "/" * prefix * filename * "_lambda" * fig_suff)
        end
    end
end

function run_tableau_list(tableaus, ode, Δt, ntime)
    for tab in tableaus
        if length(tab) ≥ 4
            nt = tab[4]
        else
            nt = ntime
        end

        if length(tab) ≥ 5
            dt = tab[5]
        else
            dt = Δt
        end

        run_lotka_volterra_2d(tab[2], ode, dt, nt; integrator=tab[1], filename=tab[3])
    end
end

function run_tableau_list_cgvi(tableaus, ode, Δt, ntime)
    for tab in tableaus
        if length(tab) ≥ 4
            nt = tab[4]
        else
            nt = ntime
        end

        if length(tab) ≥ 5
            dt = tab[5]
        else
            dt = Δt
        end

        run_lotka_volterra_2d(tab[1], tab[2], ode, dt, nt; filename=tab[3])
    end
end

ext: jl


text:
using BayesianDataFusion
using Base.Test

X = sprand(50, 100, 0.1)
Y1 = rand(50, 3)
Y2 = rand(100, 2)

Z1 = X' * Y1
Z2 = X  * Y2
addprocs(2)
@everywhere using BayesianDataFusion

fp  = psparse(X, workers())
Z1p = At_mul_B(fp, Y1)
Z2p = fp * Y2

@test_approx_eq Z1 Z1p
@test_approx_eq Z2 Z2p

W  = sprand(50, 10, 0.1)
rd = RelationData(W, class_cut = 0.5, feat1 = fp)
assignToTest!(rd.relations[1], 2)

result = macau(rd, burnin = 2, psamples = 2, num_latent=5)

@test size(rd.entities[1].model.beta) == (size(fp,2), 5)

########### parallel sparse with CSR ###########
fp2 = psparse(sparse_csr(X), workers())

Z1r = At_mul_B(fp2, Y1)
Z2r = fp2 * Y2

@test_approx_eq Z1 Z1r
@test_approx_eq Z2 Z2r

rd2 = RelationData(W, class_cut = 0.5, feat1 = fp2)
assignToTest!(rd2.relations[1], 2)

result = macau(rd2, burnin = 2, psamples = 2, num_latent=5)

rmprocs( workers() )


text:
"""
    GammaCount(m, s)

"""
struct GammaCount{T<:Real} <: DiscreteUnivariateDistribution
    m::T
    s::T
    GammaCount{T}(m::Real, s::Real) where {T <: Real} = new{T}(m, s)
end

function GammaCount(m::T, s::T; check_args=true) where {T <: Real}
    check_args && @check_args(GammaCount, m > zero(m) && s > zero(s))
    return GammaCount{T}(m, s)
end

GammaCount(m::Real, s::Real) = GammaCount(promote(m, s)...)
GammaCount(m::Integer, s::Integer) = GammaCount(float(m), float(s))
# Poisson() = Poisson(1.0, check_args=false)

@distr_support GammaCount 0 Inf

#### Conversions
convert(::Type{GammaCount{T}}, m::S, s::S) where {T <: Real, S <: Real} = GammaCount(T(m), T(s))
convert(::Type{GammaCount{T}}, d::GammaCount{S}) where {T <: Real, S <: Real} = GammaCount(T(d.m), T(d.s), check_args=false)

### Parameters

params(d::GammaCount) = (d.m, d.s)
partype(::GammaCount{T}) where {T} = T

# rate(d::Poisson) = d.λ

### Evaluation

# @_delegate_statsfuns Poisson pois λ

function pdf(d::GammaCount, x::T) where {T <: Real} # where {T <: Integer}
    if x < 0 || x != floor(x) || isinf(x)
        return 0.0
    elseif x == 0
        return ccdf.(Gamma((0+1)*d.s, 1), d.m*d.s) # 1-cdf.(Gamma((0+1)*d.s, 1), d.m*d.s)
    else
        return cdf.(Gamma((floor(x))*d.s, 1), d.m*d.s) - cdf.(Gamma((floor(x)+1)*d.s, 1), d.m*d.s)
    end
end

pdf(d::GammaCount, x::T) where {T <: Integer} = pdf(d, convert(Float64, x))

function logpdf(d::GammaCount, x::T) where {T <: Real} # where {T <: Integer}
    if x < 0 || x != floor(x) || isinf(x)
        return -Inf
    elseif x == 0
        return logccdf.(Gamma((0+1)*d.s, 1), d.m*d.s)
    else
        return logcdf.(Gamma((floor(x))*d.s, 1), d.m*d.s) + log1mexp.(logcdf.(Gamma((floor(x)+1)*d.s, 1), d.m*d.s) - logcdf.(Gamma((floor(x))*d.s, 1), d.m*d.s))
    end
end

logpdf(d::GammaCount, x::T) where {T <: Integer} = logpdf(d, convert(Float64, x))

function cdf(d::GammaCount, x::T) where {T <: Real} # where {T <: Integer}
    if isinf(x)
        return 1.0
    elseif x < 0
        return 0.0
    else
        return ccdf.(Gamma((floor(x)+1)*d.s, 1), d.m*d.s)
    end
end

cdf(d::GammaCount, x::T) where {T <: Integer} = cdf(d, convert(Float64, x))

function logcdf(d::GammaCount, x::T) where {T <: Real} # where {T <: Integer}
    if isinf(x)
        return 0.0
    elseif x < 0
        return -Inf
    else
        return log1mexp.(logcdf.(Gamma((floor(x)+1)*d.s, 1), d.m*d.s))
    end
end

logcdf(d::GammaCount, x::T) where {T <: Integer} = logcdf(d, convert(Float64, x))

function quantile(d::GammaCount, q::Real)
    if q <= 0 || Distributions.cdf.(d, 0.0) >= q
        return 0
    elseif q >= 1
        return Inf
    else
        return _solve_discrete_quantile(d, q)
    end
end

function rand(rng::AbstractRNG, d::GammaCount)
    # u = 1 - rand(rng)
    return quantile(d, 1 - Base.rand(rng))
end

function _gc_moments(d::GammaCount, m)
    upper_finite = 250 # convert(Int, quantile(d, 1-1e-10))
    series = 0:upper_finite
    # return sum( (series.^m) .* LRMoE.pdf.(d, series) )[1]
    return sum( (series.^m) .* Distributions.pdf.(d, series) )[1]
end

### Statistics

mean(d::GammaCount) = _gc_moments(d, 1)

# mode(d::GammaCount) = floor(Int,d.λ)

# function modes(d::Poisson)
#     λ = d.λ
#     isinteger(λ) ? [round(Int, λ) - 1, round(Int, λ)] : [floor(Int, λ)]
# end

var(d::GammaCount) = _gc_moments(d, 2) - (_gc_moments(d, 1))^2

function skewness(d::GammaCount)
    m1, m2, m3 = _gc_moments(d, 1), _gc_moments(d, 2), _gc_moments(d, 3)
    return (m3 - 3*m1*m2 - m1^3)/(m2^1.5)
end

function kurtosis(d::GammaCount)
    m1, m2, m3, m4 = _gc_moments(d, 1), _gc_moments(d, 2), _gc_moments(d, 3), _gc_moments(d, 4)
    return (m4 - 4*m1*m3 + 6*(m1^2)*m2 - 3*(m1)^4)/(m2^2)
end

# struct RecursiveGammaCountProbEvaluator <: RecursiveProbabilityEvaluator
#     m::Real
#     s::Real
# end

# RecursiveGammaCountProbEvaluator(d::GammaCount) = RecursiveGammaCountProbEvaluator(d.m, d.s)
# nextpdf(d::RecursiveGammaCountProbEvaluator, p::Float64, x::Integer) = p - ccdf.(Gamma((floor(x)+1)*d.s, 1), d.m*d.s) + ccdf.(Gamma((floor(x)+2)*d.s, 1), d.m*d.s)

# Base.broadcast!(::typeof(pdf), r::AbstractArray, d::GammaCount, rgn::UnitRange) =
#     _pdf!(r, d, rgn, RecursiveGammaCountProbEvaluator(d))
# function Base.broadcast(::typeof(pdf), d::GammaCount, X::UnitRange)
#     r = similar(Array{promote_type(partype(d), eltype(X))}, axes(X))
#     r .= pdf.(Ref(d),X)
# end

# function entropy(d::Poisson{T}) where T<:Real
#     λ = rate(d)
#     if λ == zero(T)
#         return zero(T)
#     elseif λ < 50
#         s = zero(T)
#         λk = one(T)
#         for k = 1:100
#             λk *= λ
#             s += λk * loggamma(k + 1) / gamma(k + 1)
#         end
#         return λ * (1 - log(λ)) + exp(-λ) * s
#     else
#         return log(2 * pi * ℯ * λ)/2 -
#                (1 / (12 * λ)) -
#                (1 / (24 * λ * λ)) -
#                (19 / (360 * λ * λ * λ))
#     end
# end

# function mgf(d::Poisson, t::Real)
#     λ = rate(d)
#     return exp(λ * (exp(t) - 1))
# end

# function cf(d::Poisson, t::Real)
#     λ = rate(d)
#     return exp(λ * (cis(t) - 1))
# end

# ### Fitting

# struct PoissonStats <: SufficientStats
#     sx::Float64   # (weighted) sum of x
#     tw::Float64   # total sample weight
# end

# suffstats(::Type{<:Poisson}, x::AbstractArray{T}) where {T<:Integer} = PoissonStats(sum(x), length(x))

# function suffstats(::Type{<:Poisson}, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Integer
#     n = length(x)
#     n == length(w) || throw(DimensionMismatch("Inconsistent array lengths."))
#     sx = 0.
#     tw = 0.
#     for i = 1 : n
#         @inbounds wi = w[i]
#         @inbounds sx += x[i] * wi
#         tw += wi
#     end
#     PoissonStats(sx, tw)
# end

# fit_mle(::Type{<:Poisson}, ss::PoissonStats) = Poisson(ss.sx / ss.tw)

## samplers

# const poissonsampler_threshold = 6

# function sampler(d::Poisson)
#     if rate(d) < poissonsampler_threshold
#         return PoissonCountSampler(d)
#     else
#         return PoissonADSampler(d)
#     end
# end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.test.behaviour.connection;
#
# import typedb.client.TypeDB;
# import typedb.client.api.TypeDBClient;
# import typedb.common.test.server.TypeDBClusterRunner;
# import typedb.common.test.server.TypeDBSingleton;
# import io.cucumber.java.After;
# import io.cucumber.java.Before;
# import io.cucumber.java.en.Given;
#
# import java.io.IOException;
# import java.util.concurrent.TimeoutException;
#
# public class ConnectionStepsCluster extends ConnectionStepsBase {
#     private TypeDBClusterRunner server;
#
#     @Override
#     void beforeAll() {
#         try {
#             server = new TypeDBClusterRunner();
#         } catch (InterruptedException | TimeoutException | IOException e) {
#             throw new RuntimeException(e);
#         }
#         server.start();
#         TypeDBSingleton.setTypeDBRunner(server);
#     }
#
#     @Before
#     public synchronized void before() {
#         super.before();
#     }
#
#     @After
#     public synchronized void after() {
#         super.after();
#     }
#
#     @Override
#     TypeDBClient createTypeDBClient(String address) {
#         return TypeDB.coreClient(address);
#     }
#
#     @Given("connection has been opened")
#     public void connection_has_been_opened() {
#         super.connection_has_been_opened();
#     }
#
#     @Given("connection does not have any database")
#     public void connection_does_not_have_any_database() {
#         super.connection_does_not_have_any_database();
#     }
#
# }


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

using FEMBase
using FEMBase: SparseMatrixCOO, SparseVectorCOO
using FEMBase: resize_sparse, resize_sparsevec
using FEMBase: get_nonzero_rows, get_nonzero_columns
using Test
using SparseArrays

@testset "Add to SparseMatrixCOO" begin
    A = SparseMatrixCOO()
    A2 = reshape(collect(1:9), 3, 3)
    add!(A, sparse(A2))
    @test isapprox(A, A2)
end

@testset "Add to SparseVectorCOO" begin
    b = SparseVectorCOO()
    b2 = collect(1:3)
    add!(b, sparse(b2))
    @test isapprox(b, b2)
end

@testset "Failure to add data to sparse vector due dimension mismatch" begin
    b = SparseVectorCOO()
    @test_throws ErrorException add!(b, [1, 2], [1.0, 2.0, 3.0])
end

@testset "resize of sparse matrix and sparse vector" begin
    A = sparse(rand(3, 3))
    B = resize_sparse(A, 4, 4)
    @test size(B) == (4, 4)

    a = sparse(rand(3))
    b = resize_sparsevec(a, 4)
    @test size(b) == (4, )
end

@testset "type conversions" begin
    A = reshape(collect(1:9), 3, 3)
    A2 = convert(SparseMatrixCOO, A)
    A3 = convert(SparseMatrixCOO, sparse(A))
    @test isapprox(A, A2)
    @test isapprox(A, A3)
    b = sparse(collect(1:3))
    b2 = convert(SparseVectorCOO, b)
    @test isapprox(b, b2)
end

@testset "construct sparse matrix and empty it" begin
    A = SparseMatrixCOO()
    add!(A, 1, 1, 1.0)
    A_expected = sparse([1], [1], [1.0], 2, 2)
    A2 = sparse(A, 2, 2)
    @test isapprox(A2, A_expected)
    empty!(A)
    @test length(A.I) == length(A.J) == length(A.V) == 0
end

@testset "construct sparse matrix with special function" begin
    A = SparseMatrixCOO()
    add!(A, 1, 1, 1.0)
    add!(A, 1, 1, 2.0)
    B = sparse(A, 2, 2, +)
    @test isapprox(B[1,1], 3.0)
end

@testset "normal addition to sparse vector" begin
    A = SparseMatrixCOO()
    add!(A, [1, 2], [1.0 2.0])
end

@testset "get nonzero rows, columns, size" begin
    A = SparseMatrixCOO()
    add!(A, [1, 3], [1, 4], [1.0 2.0; 3.0 4.0])
    B = sparse(A)
    @test get_nonzero_rows(B) == [1, 3]
    @test get_nonzero_columns(B) == [1, 4]
    @test size(A) == (3, 4)
    @test size(A, 1) == 3
end

@testset "isapprox" begin
    A = SparseMatrixCOO()
    B = SparseMatrixCOO()
    add!(A, [1, 3], [1, 4], [1.0 2.0; 3.0 4.0])
    add!(B, [1, 3], [1, 4], [1.0 2.0; 3.0 4.0])
    @test isapprox(A, B)
end


text:
using LNR

function get_cursor(s)
  ls = Jewel.lines(s)
  for (i, l) in enumerate(ls)
    m = match(r".*?\|", l)
    m == nothing || return LNR.cursor(i, length(m.match))
  end
  error("no cursor found")
end

rem_cursor(s) = replace(s, r"\|", "")

to_cursor(s) = (rem_cursor(s), get_cursor(s))

function get_tokens(s)
  s, c = to_cursor(s)
  Jewel.tokens(s, c)
end

function get_scopes(s)
  s, c = to_cursor(s)
  Jewel.scopes(s, c)
end

function get_scope(s)
  sc = get_scopes(s)[end]
  if sc.kind in (:call, :block)
    sc.name
  else
    sc
  end
end


text:
# A* shortest-path algorithm

#################################################
#
#  A* shortest-path algorithm
#
#################################################

module AStar

# The enqueue! and dequeue! methods defined in Base.Collections (needed for
# PriorityQueues) conflict with those used for queues. Hence we wrap the A*
# code in its own module.

using Graphs
using Base.Collections
using Compat

export shortest_path

function a_star_impl!{V,D}(
    graph::AbstractGraph{V},# the graph
    frontier,               # an initialized heap containing the active vertices
    colormap::Vector{Int},  # an (initialized) color-map to indicate status of vertices
    edge_dists::AbstractEdgePropertyInspector{D},  # cost of each edge
    heuristic::Function,    # heuristic fn (under)estimating distance to target
    t::V)  # the end vertex

    while !isempty(frontier)
        (cost_so_far, path, u) = dequeue!(frontier)
        if u == t
            return path
        end

        for edge in out_edges(u, graph)
            v = target(edge)
            if colormap[v] < 2
                colormap[v] = 1
                new_path = cat(1, path, edge)
                path_cost = cost_so_far + edge_property(edge_dists, edge, graph)
                enqueue!(frontier,
                        (path_cost, new_path, v),
                        path_cost + heuristic(v))
            end
        end
        colormap[u] = 2
    end
    nothing
end

function shortest_path{V,E,D}(
    graph::AbstractGraph{V,E},  # the graph
    edge_dists::AbstractEdgePropertyInspector{D},      # cost of each edge
    s::V,                       # the start vertex
    t::V,                       # the end vertex
    heuristic::Function = n -> 0)
            # heuristic (under)estimating distance to target
    frontier = VERSION < v"0.4-" ? PriorityQueue{@compat(Tuple{D,Array{E,1},V}),D}() : PriorityQueue(@compat(Tuple{D,Array{E,1},V}),D)
    frontier[(zero(D), E[], s)] = zero(D)
    colormap = zeros(Int, num_vertices(graph))
    colormap[s] = 1
    a_star_impl!(graph, frontier, colormap, edge_dists, heuristic, t)
end

function shortest_path{V,E,D}(
    graph::AbstractGraph{V,E},  # the graph
    edge_dists::Vector{D},      # cost of each edge
    s::V,                       # the start vertex
    t::V,                       # the end vertex
    heuristic::Function = n -> 0)
    edge_len::AbstractEdgePropertyInspector{D} = VectorEdgePropertyInspector(edge_dists)
    shortest_path(graph, edge_len, s, t, heuristic)
end

end

using .AStar


text:
import Base: print
export latex

using Printf

print(io::IO,A::Union{BosonDestroy,FermionDestroy}) = print(io,"$(A.name)($(A.inds...))")
print(io::IO,A::Union{BosonCreate,FermionCreate}) = print(io,"$(A.name)†($(A.inds...))")
print(io::IO,A::σ) = print(io,"σ$(A.a)($(A.inds...))")
print(io::IO,A::σminus) = print(io,"σ-($(A.inds...))")
print(io::IO,A::σplus) = print(io,"σ+($(A.inds...))")
print(io::IO,A::scal) = print(io,imag(A.v)==0 ? @sprintf("%g",real(A.v)) : (real(A.v)==0 ? @sprintf("%gi",imag(A.v)) : @sprintf("(%g%+gi)",real(x),imag(x))))
print(io::IO,A::param) = print(io,string(A.name, A.state=='c' ? "'" : "",length(A.inds)==0 ? "" : "($(A.inds...))"))
print(io::IO,A::OpProd) = print(io,"$(A.A) $(A.B)")
print(io::IO,A::OpSum) = print(io,"$(A.A) + $(A.B)")
print(io::IO,A::ExpVal) = print(io,"⟨$(A.A)⟩")
print(io::IO,A::Corr) = print(io,"⟨$(A.A)⟩c")
print(io::IO,A::OpSumAnalytic) = print(io,"∑_$(A.ind) $(A.A)")

mystring(x::Number) = @sprintf "%g" x
mystring(x::Rational) = denominator(x)==1 ? "$(numerator(x))" : "\\frac{$(numerator(x))}{$(denominator(x))}"
mystring(x::Complex) = @sprintf "(%g%+gi)" real(x) imag(x)
mystring(x::Complex{Rational{T}}) where T = @sprintf "\\left(%s%s%si\\right)" mystring(real(x)) (imag(x)>=0 ? "+" : "-") mystring(abs(imag(x)))

Base.show(io::IO, ::MIME"text/latex", A::Operator) = print(io,"\$",latex(A),"\$")
latex(A::σ) = string("\\sigma_{$(A.a)",length(A.inds)>0 ? ",$(A.inds...)}" : "}")
latexindstr(inds::OpIndices) = length(inds)==0 ? "" : "_{$(inds...)}"
latex(A::δ) = "δ" * latexindstr(A.inds)
latex(A::Union{BosonDestroy,FermionDestroy}) = "{$(A.name)}$(latexindstr(A.inds))"
latex(A::Union{BosonCreate,FermionCreate}) = "{$(A.name)}$(latexindstr(A.inds))^\\dagger"
latex(A::σminus) = "\\sigma^-" * latexindstr(A.inds)
latex(A::σplus) = "\\sigma^+" * latexindstr(A.inds)
latex(A::scal) = imag(A.v)==0 ? mystring(real(A.v)) : (real(A.v)==0 ? mystring(imag(A.v))*"i" : mystring(A.v))
latex(A::param)  = string(A.name, latexindstr(A.inds), A.state=='c' ? "^*" : "")
latex(A::OpProd) = string(latex(A.A)," ",latex(A.B))
latex(A::OpSum) = string(latex(A.A)," + ",latex(A.B))
latex(A::ExpVal) = "\\langle $(latex(A.A)) \\rangle"
latex(A::Corr) = "\\langle $(latex(A.A)) \\rangle_{c}"
latex(A::OpSumAnalytic) = string("\\sum_{$(A.ind)}",latex(A.A))


text:
@testset "Handles" begin
    @testset "Definitions" begin
        test(HandleDefinition, handle_by_name, :VkInstance, :(
            mutable struct Instance <: Handle
                vks::VkInstance
                refcount::RefCounter
                destructor
                Instance(vks::VkInstance, refcount::RefCounter) = new(vks, refcount, undef)
            end
        ))

        test(HandleDefinition, handle_by_name, :VkPhysicalDevice, :(
            mutable struct PhysicalDevice <: Handle
                vks::VkPhysicalDevice
                instance::Instance
                refcount::RefCounter
                destructor
                PhysicalDevice(vks::VkPhysicalDevice, instance::Instance, refcount::RefCounter) = new(vks, instance, refcount, undef)
            end
        ))

        test(HandleDefinition, handle_by_name, :VkDevice, :(
            mutable struct Device <: Handle
                vks::VkDevice
                physical_device::PhysicalDevice
                refcount::RefCounter
                destructor
                Device(vks::VkDevice, physical_device::PhysicalDevice, refcount::RefCounter) = new(vks, physical_device, refcount, undef)
            end
        ))
    end

    @testset "Constructors" begin
        test_ex(
            Constructor(
                HandleDefinition(handle_by_name(:VkInstance)),
                APIFunction(create_func(:vkCreateInstance), false),
            ),
            :(
                Instance(enabled_layer_names::AbstractArray, enabled_extension_names::AbstractArray; allocator = C_NULL, next=C_NULL, flags=0, application_info=C_NULL) = unwrap(create_instance(enabled_layer_names, enabled_extension_names; allocator, next, flags, application_info))
            )
        )

        test_ex(
            Constructor(
                HandleDefinition(handle_by_name(:VkDebugReportCallbackEXT)),
                APIFunction(create_func(:vkCreateDebugReportCallbackEXT), true),
            ),
            :(
                DebugReportCallbackEXT(instance, pfn_callback::FunctionPtr, fptr_create::FunctionPtr, fptr_destroy::FunctionPtr; allocator = C_NULL, next = C_NULL, flags = 0, user_data = C_NULL) = unwrap(create_debug_report_callback_ext(instance, pfn_callback, fptr_create, fptr_destroy; allocator, next, flags, user_data))
            )
        )

        test_ex(
            Constructor(
                HandleDefinition(handle_by_name(:VkDeferredOperationKHR)),
                APIFunction(func_by_name(:vkCreateDeferredOperationKHR), false),
            ),
            :(
                DeferredOperationKHR(device; allocator = C_NULL) = unwrap(create_deferred_operation_khr(device; allocator))
            )
        )

        test_ex(
            Constructor(
                HandleDefinition(handle_by_name(:VkRenderPass)),
                promote_hl(APIFunction(create_func(:vkCreateRenderPass), false)),
            ),
            :(
                RenderPass(device, attachments::AbstractArray, subpasses::AbstractArray, dependencies::AbstractArray; allocator = C_NULL, next = C_NULL, flags = 0) = unwrap(create_render_pass(device, attachments, subpasses, dependencies; allocator, next, flags))
            ),
        )

        test_ex(
            Constructor(
                HandleDefinition(handle_by_name(:VkRenderPass)),
                promote_hl(APIFunction(create_func(:vkCreateRenderPass2), false)),
            ),
            :(
                RenderPass(device, attachments::AbstractArray, subpasses::AbstractArray, dependencies::AbstractArray, correlated_view_masks::AbstractArray; allocator = C_NULL, next = C_NULL, flags = 0) = unwrap(create_render_pass_2(device, attachments, subpasses, dependencies, correlated_view_masks; allocator, next, flags))
            ),
        )
    end
end


text:
using NumFormat
using Base.Test

function test_equality()
    println( "test equality...")
    srand(10)
    fmts = ASCIIString[ "%10.4f", "%f", "%e", "%10f", "%.3f", "%.3e" ]
    for fmt in fmts
        l = :( x-> x )
        l.args[2].args[2] = Expr( :macrocall, symbol( "@sprintf" ), fmt, :x )
        mfmtr = eval( l )
        for i in 1:10000
            n = erfinv( rand() * 1.99 - 1.99/2.0 )
            expect = mfmtr( n )
            actual = sprintf1( fmt, n )
            @test expect == actual
        end
    end

    fmts = ASCIIString[ "%d", "%10d", "%010d", "%-10d" ]
    for fmt in fmts
        l = :( x-> x )
        l.args[2].args[2] = Expr( :macrocall, symbol( "@sprintf" ), fmt, :x )
        mfmtr = eval( l )
        for i in 1:10000
            j = int( erfinv( rand() * 1.99 - 1.99/2.0 ) * 100000 )
            expect = mfmtr( j )
            actual = sprintf1( fmt, j )
            @test expect == actual
        end
    end
    println( "...Done" )
end

@time test_equality()

println( "\nTest speed" )

function native_int()
    for i in 1:200000
        @sprintf( "%10d", i )
    end
end
function runtime_int()
    for i in 1:200000
        sprintf1( "%10d", i )
    end
end
function runtime_int_bypass()
    f = generate_formatter( "%10d" )
    for i in 1:200000
        f( i )
    end
end

println( "integer @sprintf speed")
@time native_int()
println( "integer sprintf speed")
@time runtime_int()
println( "integer sprintf speed, bypass repeated lookup")
@time runtime_int_bypass()

function native_float()
    srand( 10 )
    for i in 1:200000
        @sprintf( "%10.4f", erfinv( rand() ) )
    end
end
function runtime_float()
    srand( 10 )
    for i in 1:200000
        sprintf1( "%10.4f", erfinv( rand() ) )
    end
end
function runtime_float_bypass()
    f = generate_formatter( "%10.4f" )
    srand( 10 )
    for i in 1:200000
        f( erfinv( rand() ) )
    end
end

println()
println( "float64 @sprintf speed")
@time native_float()
println( "float64 sprintf speed")
@time runtime_float()
println( "float64 sprintf speed, bypass repeated lookup")
@time runtime_float_bypass()

function test_commas()
    println( "\ntest commas..." )
    @test sprintf1( "%'d", 1000 ) == "1,000"
    @test sprintf1( "%'d", -1000 ) == "-1,000"
    @test sprintf1( "%'d", 100 ) == "100"
    @test sprintf1( "%'d", -100 ) == "-100"
    @test sprintf1( "%'f", Inf ) == "Inf"
    @test sprintf1( "%'f", -Inf ) == "-Inf"
    @test sprintf1( "%'s", 1000.0 ) == "1,000.0"
    @test sprintf1( "%'s", 1234567.0 ) == "1.234567e6"
end

function test_format()
    println( "test format...")
    @test format( 10 ) == "10"
    @test format( 10.0 ) == "10"
    @test format( 10.0, precision=2 ) == "10.00"
    @test format( 111//100, precision=2 ) == "1.11"
    @test format( 111//100 ) == "111/100"
    @test format( 1234, commas=true ) == "1,234"
    @test format( 1234, conversion="f", precision=2 ) == "1234.00"
    @test format( 1.23, precision=3 ) == "1.230"
    @test format( 1.23, precision=3, stripzeros=true ) == "1.23"
    @test format( 1.00, precision=3, stripzeros=true ) == "1"

    @test format( 1.0, conversion="e", stripzeros=true ) == "1e+00"
    @test format( 1.0, conversion="e", precision=4 ) == "1.0000e+00"

    # hex output
    @test format( 1118, conversion="x" ) == "45e"
    @test format( 1118, width=4, conversion="x" ) == " 45e"
    @test format( 1118, width=4, zeropadding=true, conversion="x" ) == "045e"
    @test format( 1118, alternative=true, conversion="x" ) == "0x45e"
    @test format( 1118, width=4, alternative=true, conversion="x" ) == "0x45e"
    @test format( 1118, width=6, alternative=true, conversion="x", zeropadding=true ) == "0x045e"

    # mixed fractions
    @test format( 3//2, mixedfraction=true ) == "1_1/2"
    @test format( -3//2, mixedfraction=true ) == "-1_1/2"
    @test format( 3//100, mixedfraction=true ) == "3/100"
    @test format( -3//100, mixedfraction=true ) == "-3/100"
    @test format( 307//100, mixedfraction=true ) == "3_7/100"
    @test format( -307//100, mixedfraction=true ) == "-3_7/100"
    @test format( 307//100, mixedfraction=true, fractionwidth=6 ) == "3_07/100"
    @test format( -307//100, mixedfraction=true, fractionwidth=6 ) == "-3_07/100"
    @test format( -302//100, mixedfraction=true ) == "-3_1/50"
    # try to make the denominator 100
    @test format( -302//100, mixedfraction=true,tryden = 100 ) == "-3_2/100"
    @test format( -302//30, mixedfraction=true,tryden = 100 ) == "-10_1/15" # lose precision otherwise
    @test format( -302//100, mixedfraction=true,tryden = 100,fractionwidth=6 ) == "-3_02/100" # lose precision otherwise

    #commas
    @test format( 12345678, width=10, commas=true ) == "12,345,678"
    # it would try to squeeze out the commas
    @test format( 12345678, width=9, commas=true ) == "12345,678"
    # until it can't anymore
    @test format( 12345678, width=8, commas=true ) == "12345678"
    @test format( 12345678, width=7, commas=true ) == "12345678"

    # only the numerator would have commas
    @test format( 1111//1000, commas=true ) == "1,111/1000"

    # this shows how, with enough space, parens line up with empty spaces
    @test format(  12345678, width=12, commas=true, parens=true )== " 12,345,678 "
    @test format( -12345678, width=12, commas=true, parens=true )== "(12,345,678)"
    # same with unspecified width
    @test format(  12345678, commas=true, parens=true )== " 12,345,678 "
    @test format( -12345678, commas=true, parens=true )== "(12,345,678)"
end

test_commas()
test_format()
println( "done")


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Two Body orbit propagator.
#
#   This algorithm considers a perfect Keplerian orbit. In other words, no
#   perturbation is considered during the propagation and the Earth is modeled
#   as a perfect sphere.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export TwoBody_Structure
export twobody_init, twobody!

################################################################################
#                                  Functions
################################################################################

"""
    twobody_init(epoch::Tepoch, a_0::Number, e_0::Number, i_0::Number, Ω_0::Number, ω_0::Number, f_0::Number; μ::T = m0) where {Tepoch, T}

Initialize the data structure of two body orbit propagator algorithm.

!!! note
    The type used in the propagation will be the same as used to define the
    gravitational constant `μ`.

# Args

- `epoch::Number`: Epoch of the initial mean orbital elements [Julian Day].
- `a_0::Number`: Initial mean semi-major axis [m].
- `e_0::Number`: Initial mean eccentricity.
- `i_0::Number`: Initial mean inclination [rad].
- `Ω_0::Number`: Initial mean right ascension of the ascending node [rad].
- `ω_0::Number`: Initial mean argument of perigee [rad].
- `f_0::Number`: Initial mean true anomaly [rad].

# Keywords

- `μ::T`: Standard gravitational parameter of the central body [m^3/s^2].
    (**Default** = `m0`)

# Returns

The structure [`TwoBody_Structure`](@ref) with the initialized parameters.
"""
function twobody_init(
    epoch::Tepoch,
    a_0::Number,
    e_0::Number,
    i_0::Number,
    Ω_0::Number,
    ω_0::Number,
    f_0::Number;
    μ::T = m0
) where {Tepoch, T}
    # The propagator is only defined for 0 <= e < 1.
    if !(0 <= e_0 < 1)
        throw(ArgumentError("The two body propagator only supports eccentricities in the interval [0,1)"))
    end

    # Compute the mean motion using the semi-major axis.
    n_0 = sqrt(μ / T(a_0)^3)

    # Compute the initial mean anomaly.
    M_0 = f_to_M(T(e_0), T(f_0))

    # Create and return the Two Body orbit propagator structure.
    TwoBody_Structure{Tepoch, T}(
        epoch = epoch,
        a_0   = a_0,
        n_0   = n_0,
        e_0   = e_0,
        i_0   = i_0,
        Ω_0   = Ω_0,
        ω_0   = ω_0,
        M_0   = M_0,
        f_0   = f_0,
        μ     = μ,
        Δt    = 0,
        f_k   = f_0,
        M_k   = M_0
    )
end

"""
    twobody!(tbd::TwoBody_Structure{Tepoch, T}, t::Number) where {Tepoch, T}

Propagate the orbit defined in `tbd` (see [`TwoBody_Structure`](@ref)) until the
time `t` [s].

!!! note
    The internal values in `tbd` will be modified.

# Returns

- The position vector represented in the inertial frame at time `t` [m].
- The velocity vector represented in the inertial frame at time `t` [m/s]

# Remarks

The inertial frame in which the output is represented depends on which frame it
was used to generate the orbit parameters.
"""
function twobody!(tbd::TwoBody_Structure{Tepoch, T}, t::Number) where {Tepoch, T}
    # Time elapsed since epoch.
    Δt     = T(t)
    tbd.Δt = Δt

    # Update the mean anomaly.
    tbd.M_k = tbd.M_0 + tbd.n_0 * Δt

    # Convert the mean anomaly to true anomaly.
    tbd.f_k = M_to_f(tbd.e_0, tbd.M_k)

    # Compute the position and velocity vectors given the orbital elements.
    r_i_k, v_i_k = kepler_to_rv(tbd.a_0, tbd.e_0, tbd.i_0, tbd.Ω_0, tbd.ω_0, tbd.f_k)

    # Return the position and velocity vector represented in the inertial
    # reference frame.
    return r_i_k, v_i_k
end


text:
# Bogumił Kamiński, 2022

# Codes for chapter 7

# Code for listing 7.1

using HTTP
using JSON3
query = "https://api.nbp.pl/api/exchangerates/rates/a/usd/" *
        "2020-06-01/?format=json"
response = HTTP.get(query)
json = JSON3.read(response.body)

# Code for the remainder of section 7.1.2

response.body

String(response.body)

response.body

json.table
json.currency
json.code
json.rates

json.rates[1].mid

only(json.rates).mid

only([])
only([1, 2])

# Code for listing 7.2

query = "https://api.nbp.pl/api/exchangerates/rates/a/usd/" *
        "2020-06-06/?format=json"
response = HTTP.get(query)

# Code for listing 7.3

query = "https://api.nbp.pl/api/exchangerates/rates/a/usd/" *
        "2020-06-01/?format=json"
try
    response = HTTP.get(query)
    json = JSON3.read(response.body)
    only(json.rates).mid
catch e
    if e isa HTTP.ExceptionRequest.StatusError
        missing
    else
        rethrow(e)
    end
end

query = "https://api.nbp.pl/api/exchangerates/rates/a/usd/" *
        "2020-06-06/?format=json"
try
    response = HTTP.get(query)
    json = JSON3.read(response.body)
    only(json.rates).mid
catch e
    if e isa HTTP.ExceptionRequest.StatusError
        missing
    else
        rethrow(e)
    end
end

# Code for section 7.2

ismissing(missing)
ismissing(1)

1 + missing
sin(missing)

1 == missing
1 > missing
1 < missing

if missing
    print("this is not printed")
end
missing && true

coalesce(missing, true)
coalesce(missing, false)

isequal(1, missing)
isequal(missing, missing)
isless(1, missing)
isless(missing, missing)

isless(Inf, missing)

a = [1]
b = [1]
isequal(a, b)
a === b

x = [1, missing, 3, 4, missing]

coalesce.(x, 0)

sum(x)

y = skipmissing(x)

sum(y)

sum(skipmissing(x))

fun(x::Int, y::Int) = x + y
fun(1, 2)
fun(1, missing)

using Missings
fun2 = passmissing(fun)
fun2(1, 2)
fun2(1, missing)

# Code for section 7.3

using Dates
d = Date("2020-06-01")

typeof(d)
year(d)
month(d)
day(d)

dayofweek(d)
dayname(d)

Date(2020, 6, 1)

dates = Date.(2020, 6, 1:30)

Day(1)

d
d + Day(1)

Date(2020, 5, 20):Day(1):Date(2020, 7, 5)

collect(Date(2020, 5, 20):Day(1):Date(2020, 7, 5))

# Code for listing 7.6

function get_rate(date::Date)
    query = "https://api.nbp.pl/api/exchangerates/rates/" *
            "a/usd/$date/?format=json"
    try
        response = HTTP.get(query)
        json = JSON3.read(response.body)
        return only(json.rates).mid
    catch e
        if e isa HTTP.ExceptionRequest.StatusError
            return missing
        else
            rethrow(e)
        end
    end
end

# Code for showing how string interpolation works

"https://api.nbp.pl/api/exchangerates/rates/" *
"a/usd/$(dates[1])/?format=json"

"https://api.nbp.pl/api/exchangerates/rates/" *
"a/usd/$dates[1]/?format=json"

# Code for listing 7.7

rates = get_rate.(dates)

# Code for section 7.4

using Statistics
mean(rates)
std(rates)

mean(skipmissing(rates))
std(skipmissing(rates))

# Code for listing 7.8

using FreqTables
proptable(dayname.(dates), ismissing.(rates); margins=1)

# Code showing how to specify a complex condition using broadcasting

dayname.(dates) .== "Thursday" .&& ismissing.(rates)

# Code for listing 7.9

dates[dayname.(dates) .== "Thursday" .&& ismissing.(rates)]

# Codes for plotting exchange rate data

using Plots
plot(dates, rates; xlabel="day", ylabel="PLN/USD", legend=false)

rates_ok = .!ismissing.(rates)

plot(dates[rates_ok], rates[rates_ok];
     xlabel="day", ylabel="PLN/USD", legend=false)

using Impute
rates_filled = Impute.interp(rates)

scatter!(dates, rates_filled)

max_stars_repo_name: Mo-Gul/JuliaForDataAnalysis

max_issues_repo_name: Mo-Gul/JuliaForDataAnalysis

max_forks_repo_name: Mo-Gul/JuliaForDataAnalysis


text:
using ColorTypes
using ColorTypes.FixedPointNumbers
using Test

@testset "compatibility tests for ARGB32/AGray32" begin
    # cf. PR #146
    @test ARGB32 <: AlphaColor{RGB24, N0f8, 4}
    @test ARGB32 <: AbstractARGB
    @test !(ARGB32 <: AbstractRGBA)
    @test AGray32 <: AlphaColor{Gray24, N0f8, 2}
    @test AGray32 <: AbstractAGray
    @test !(AGray32 <: AbstractGrayA)
    @test supertype(ARGB32) === AlphaColor{RGB24, N0f8, 4} # v0.9 or earlier
    @test supertype(ARGB32) === AbstractARGB{RGB24, N0f8}
    @test supertype(AGray32) === AlphaColor{Gray24, N0f8, 2} # v0.9 or earlier
    @test supertype(AGray32) === AbstractAGray{Gray24, N0f8}
end

### Prevent future commits from unexporting abstract types
@testset "abstract typealiases" begin
    dispatcher2(::Colorant) = error("Colorant")
    dispatcher2(::AbstractGray) = :AbstractGray
    dispatcher2(::Color3)       = :Color3
    dispatcher2(::TransparentGray) = :TransparentGray
    dispatcher2(::Transparent3)    = :Transparent3
    dispatcher2(::TransparentRGB)  = :TransparentRGB

    @test dispatcher2(Gray(0.2))   == :AbstractGray
    @test dispatcher2(Gray24(0.2)) == :AbstractGray
    @test dispatcher2(RGB(1.0,1.0,1.0)) == :Color3
    @test dispatcher2(HSV(100,0.6,0.4)) == :Color3
    @test dispatcher2(GrayA(0.2,0.5))   == :TransparentGray
    @test dispatcher2(AGray(0.2,0.5))   == :TransparentGray
    @test dispatcher2(AGray32(0.2,0.5)) == :TransparentGray
    @test dispatcher2(HSVA(100,0.6,0.4,0.5)) == :Transparent3
    @test dispatcher2(AHSV(100,0.6,0.4,0.5)) == :Transparent3
    @test dispatcher2(RGBA(1.0,1.0,1.0,0.5))   == :TransparentRGB
    @test dispatcher2(ARGB(1.0,1.0,1.0,0.5))   == :TransparentRGB
    @test dispatcher2(ARGB32(0.5,1.0,1.0,1.0)) == :TransparentRGB

    dispatcher3(::TransparentGray) = error("TransparentGray")
    dispatcher3(::Transparent3)    = error("Transparent3")
    dispatcher3(::TransparentRGB)  = error("TransparentRGB")
    dispatcher3(::AbstractAGray) = :AbstractAGray
    dispatcher3(::AbstractGrayA) = :AbstractGrayA
    dispatcher3(::AbstractARGB) = :AbstractARGB
    dispatcher3(::AbstractRGBA) = :AbstractRGBA

    @test dispatcher3(GrayA(0.2,0.5))   == :AbstractGrayA
    @test dispatcher3(AGray(0.2,0.5))   == :AbstractAGray
    @test dispatcher3(AGray32(0.2,0.5)) == :AbstractAGray
    @test dispatcher3(RGBA(1.0,1.0,1.0,0.5))   == :AbstractRGBA
    @test dispatcher3(ARGB(1.0,1.0,1.0,0.5))   == :AbstractARGB
    @test dispatcher3(ARGB32(0.5,1.0,1.0,1.0)) == :AbstractARGB

    normeddispatcher(::Colorant) = false
    normeddispatcher(::ColorantNormed) = true
    @test normeddispatcher(RGB(1.0,1.0,1.0)) == false
    @test normeddispatcher(Gray{Bool}(1)) == false
    @test normeddispatcher(RGB{N0f8}(1.0,1.0,1.0)) == true
    @test normeddispatcher(ARGB32(1.0,1.0,1.0)) == true
end

@testset "colorfields" begin
    @test ColorTypes.colorfields(RGB) == (:r, :g, :b)
    @test ColorTypes.colorfields(RGBA) == (:r, :g, :b, :alpha)
    @test ColorTypes.colorfields(ARGB) == (:r, :g, :b, :alpha)
    @test ColorTypes.colorfields(RGB{N0f8}) == (:r, :g, :b)
    @test ColorTypes.colorfields(RGBA{Float32}) == (:r, :g, :b, :alpha)
    @test ColorTypes.colorfields(ARGB{Float64}) == (:r, :g, :b, :alpha)
    @test ColorTypes.colorfields(RGB24) == (:color,)
    @test ColorTypes.colorfields(ARGB32) == (:color, :alpha)

    @test ColorTypes.colorfields(BGR) == (:r, :g, :b)
    @test ColorTypes.colorfields(XRGB) == (:r, :g, :b)
    @test ColorTypes.colorfields(RGBX) == (:r, :g, :b)

    @test ColorTypes.colorfields(Gray) == (:val,)
    @test ColorTypes.colorfields(GrayA) == (:val, :alpha)
    @test ColorTypes.colorfields(AGray) == (:val, :alpha)
    @test ColorTypes.colorfields(Gray24) == (:color,)
    @test ColorTypes.colorfields(AGray32) == (:color, :alpha)

    @test ColorTypes.colorfields(HSV) == (:h, :s, :v)
    @test ColorTypes.colorfields(HSVA) == (:h, :s, :v, :alpha)
    @test ColorTypes.colorfields(AHSV) == (:h, :s, :v, :alpha)

    @test_throws ArgumentError ColorTypes.colorfields(AbstractRGB)
    @test_throws MethodError ColorTypes.colorfields(N0f8)

    # for instances
    @test ColorTypes.colorfields(RGB{N0f8}(1,0,0)) == (:r, :g, :b)
    @test ColorTypes.colorfields(ARGB(1.0,0.8,0.6,0.4)) == (:r, :g, :b, :alpha)
    @test ColorTypes.colorfields(RGBA{Float32}(1.0,0.8,0.6,0.4)) == (:r, :g, :b, :alpha)
    @test_throws MethodError ColorTypes.colorfields(1N0f8)
end

@testset "rgb constructors" begin
    Crgb = filter(T -> T <: AbstractRGB, ColorTypes.parametric3)
    @testset "$C constructor" for C in Crgb
        for val1 in (0.2, 0.2f0, N0f8(0.2), N4f12(0.2), N0f16(0.2))
            for val2 in (0.6, 0.6f0, N0f8(0.6), N4f12(0.6), N0f16(0.6))
                c = C(val1,val2,val1)
                @test isa(c, C)
                @test C{N0f8}(val1,val2,val1) === C(0.2N0f8,0.6N0f8,0.2N0f8)
                @test C{N0f16}(val1,val2,val1) === C(0.2N0f16,0.6N0f16,0.2N0f16)
            end
            # 1-arg constructor
            @test C(val1) === C{typeof(val1)}(0.2,0.2,0.2)
            @test C{N0f8}(val1) === C{N0f8}(0.2,0.2,0.2)
        end

        @test_throws ArgumentError C(2,1,0) # integers

        for val in (1.2, 1.2f0, N4f12(1.2), -0.2)
            @test_throws ArgumentError C{N0f8}(val,val,val)
            @test_throws ArgumentError C{N0f16}(val,val,val)
            @test isa(C(val,val,val), C)
        end
    end
end

@testset "transparent rgb constructors" begin
    Crgb = filter(T -> T <: AbstractRGB, ColorTypes.parametric3)
    # all(T -> alphacolor(T) === ARGB, (RGB, XRGB, RGBX))
    Ctransparent = unique(vcat(coloralpha.(Crgb), alphacolor.(Crgb)))
    @testset "$C constructor" for C in Ctransparent
        C0 = color_type(C)

        for val1 in (0.2, 0.2f0, N0f8(0.2), N4f12(0.2), N0f16(0.2))
            for val2 in (0.6, 0.6f0, N0f8(0.6), N4f12(0.6), N0f16(0.6))
                c0 = C0(val1,val2,val1)
                et = eltype(c0)
                @test isa(C(val1,val2,val1), C)
                @test isa(C(val1,val2,val1,0.8), C)
                @test C(c0) === C(val1,val2,val1,convert(et, 1))
                @test C(c0) === C{et}(val1,val2,val1,1)
                @test C(c0,0.8) === C(val1,val2,val1,convert(et, 0.8))
                @test C{N0f8}(val1,val2,val1) === C(0.2N0f8,0.6N0f8,0.2N0f8,1N0f8)
                @test C{N0f16}(val1,val2,val1,0.8) === C(0.2N0f16,0.6N0f16,0.2N0f16,0.8N0f16)
            end
            # 1-arg constructor
            @test_broken C(val1) === C{typeof(val1)}(0.2,0.2,0.2,1)
            @test_broken C{N0f8}(val1) === C{N0f8}(0.2,0.2,0.2,1)
        end

        @test_throws ArgumentError C(2,1,0) # integers
        @test_throws ArgumentError C(2,1,0,1) # integers

        for val in (1.2, 1.2f0, N4f12(1.2), -0.2)
            c0 = C0(val,val,val)
            @test_throws ArgumentError C{N0f8}(val,val,val)
            @test_throws ArgumentError C{N0f16}(val,val,val,0.8)
            @test_throws ArgumentError C{N0f8}(c0)
            @test_throws ArgumentError C{N0f16}(c0,0.8)
            @test isa(C(val,val,val), C)
            @test isa(C(val,val,val,0.8), C)
            @test isa(C(val,val,val,val), C) # no exception thrown
        end
    end
end

@testset "RGB24/ARGB32 constructors" begin
    for val1 in (0.2, 0.2f0, N0f8(0.2), N4f12(0.2), N0f16(0.2))
        for val2 in (0.6, 0.6f0, N0f8(0.6), N4f12(0.6), N0f16(0.6))
            @test RGB24(val1,val2,val1) === RGB24(0.2,0.6,0.2)
            @test ARGB32(val1,val2,val1) === ARGB32(0.2,0.6,0.2,1)
            @test ARGB32(val1,val2,val1,0.8) === ARGB32(0.2,0.6,0.2,0.8)
            c0 = RGB24(val1,val2,val1)
            @test ARGB32(c0) === ARGB32(val1,val2,val1,1)
            @test ARGB32(c0,0.8) === ARGB32(val1,val2,val1,0.8)
        end
        # 1-arg constructor
        @test RGB24(val1) === RGB24(0.2,0.2,0.2)
        @test ARGB32(val1) === ARGB32(0.2,0.2,0.2,1)
    end

    @test_throws ArgumentError RGB24(2,1,0) # integers
    @test_throws ArgumentError ARGB32(2,1,0) # integers
    @test_throws ArgumentError ARGB32(2,1,0,1) # integers

    # https://github.com/JuliaGraphics/ColorTypes.jl/pull/183#issuecomment-616958191
    @test_broken ARGB32(-0.00196, 0.0, 1.00196) === ARGB32(0, 0, 1)

    for val in (1.2, 1.2f0, N4f12(1.2), -0.2)
        @test_throws ArgumentError RGB24(val,val,val)
        @test_throws ArgumentError ARGB32(val,val,val)
        @test_throws ArgumentError ARGB32(val,val,val,0.8)
        c0 = RGB(val,val,val)
        @test_throws ArgumentError ARGB32(c0)
        @test_throws ArgumentError ARGB32(c0,0.8)
    end
end

@testset "RGB construction with integer arguments" begin
    ret = @test_throws ArgumentError RGB(255, 17, 48)
    @test occursin("255, 17, 48", ret.value.msg)
    @test occursin("0-255", ret.value.msg)
    ret = @test_throws ArgumentError RGB(256, 17, 48)
    @test occursin("256, 17, 48", ret.value.msg)
    @test !occursin("0-255", ret.value.msg)
    # for RGB24 and ARGB32 with UInt32 arguments, see "bit pattern ambiguities"
    # in "test/conversion.jl"
end

@testset "color construction from grayscale" begin
    @test RGB(Gray(0.2), 0.3, 0.4) === RGB(0.2, 0.3, 0.4)
    @test RGB(0.2, Gray(0.3), 0.4) === RGB(0.2, 0.3, 0.4)
    @test RGB(0.2, 0.3, Gray(0.4)) === RGB(0.2, 0.3, 0.4)
    @test RGB(Gray(0.2), Gray(0.3), Gray(0.4)) === RGB(0.2, 0.3, 0.4)
    @test RGB24(Gray(0.2), 0.3, 0.4) === RGB24(0.2, 0.3, 0.4)
    @test ARGB32(0.2, 0.3, 0.4, Gray(0.5)) === ARGB32(0.2, 0.3, 0.4, 0.5)
    @test RGB(0.2, Gray24(0.3), 0.4) === RGB(0.2, 0.3N0f8, 0.4)
    @test_throws MethodError HSV(0.2, 0.3, Gray(0.4))
end

@testset "gray constructors" begin
    for val in (0.2, 0.2f0, N0f8(0.2), N4f12(0.2), N0f16(0.2),
                Gray{N0f8}(0.2), Gray{N4f12}(0.2), Gray24(0.2))
        @test isa(Gray(val), Gray)
        @test Gray{N0f8}(val) === Gray{N0f8}(0.2)
        @test Gray{N0f16}(val) === Gray{N0f16}(0.2)
        @test Gray24(val) === Gray24(0.2)
        @test GrayA{N0f8}(val) === GrayA{N0f8}(0.2, 1)
        @test AGray{N0f16}(val) === AGray{N0f16}(0.2, 1)
        @test AGray32(val) === AGray32(0.2, 1)
        @test AGray32(val, 0.8) === AGray32(0.2, 0.8)
    end
    for val in (1.2, 1.2f0, N4f12(1.2), Gray{N4f12}(1.2), 2, -0.2)
        !isa(val, Int) && @test isa(Gray(val), Gray)
        @test_throws ArgumentError Gray{N0f8}(val)
        @test_throws ArgumentError Gray{N0f16}(val)
        @test_throws ArgumentError Gray24(val)
        @test_throws ArgumentError GrayA{N0f8}(val)
        @test_throws ArgumentError AGray{N0f16}(val)
        @test_throws ArgumentError AGray32(val)
        @test_throws ArgumentError AGray32(val, 0.8)
    end
    @test Gray() === Gray{N0f8}(0)
    @test Gray(Gray()) === Gray()  # no StackOverflowError
    @test Gray(1) === Gray{N0f8}(1)
    @test Gray(true) === Gray{Bool}(1)
    @test Gray(GrayA(0.2,0.8)) === Gray{Float64}(0.2)
    @test Gray(AGray32(0.2,0.8)) === Gray{N0f8}(0.2)
    @test Gray24(AGray(0.2,0.8)) === Gray24(0.2)
    @test Gray24(AGray32(0.2,0.8)) === Gray24(0.2)
    @test eltype(broadcast(Gray, rand(5))) == Gray{Float64}
    @test eltype(broadcast(Gray, rand(Float32,5))) == Gray{Float32}
end

@testset "parametric3 constructors" begin
    @testset "$C constructor" for C in ColorTypes.parametric3
        et = (C <: AbstractRGB) ? N0f8 : Float32
        @test isa(C(1,0,0), C)
        @test isa(C(1,0,0), C{et})
        @test isa(C{Float32}(1, 0.5, 0), C{Float32})
        @test C(1N0f8, 0.6N0f8, 0N0f8) === C{et}(1, 0.6, 0) # issue #80
        @test C() === C{et}(0,0,0)
        @test C(C()) === C()  # no StackOverflowError
    end
end

@testset "transparent3 constructors" begin
    Cp3 = ColorTypes.parametric3
    Ctransparent = unique(vcat(coloralpha.(Cp3), alphacolor.(Cp3)))
    @testset "transparent3 constructors: $C" for C in Ctransparent
        et = (C <: TransparentRGB) ? N0f8 : Float32
        @test isa(C(1,0,0), C)
        @test isa(C(1,0,0), C{et})
        @test isa(C{Float32}(1, 0.5, 0), C{Float32})
        @test C(1,0,0,0.8) === C{Float64}(1,0,0,0.8)
        @test C(1,0,0) === C{et}(1,0,0,1)
        @test C(1,0,0,1) === C{et}(1,0,0,1)
        if C <: TransparentRGB
            @test C(1N0f8, 0.6N0f8, 0N0f8) === C{et}(1, 0.6, 0, 1)
        else
            @test_broken C(1N0f8, 0.6N0f8, 0N0f8) === C{et}(1, 0.6, 0, 1) # issue #156
        end
        @test C() === C{et}(0,0,0,1)
        @test C(C()) === C()  # no StackOverflowError
    end
end

@testset "coloralpha for types" begin
    @test @inferred(coloralpha(RGB)) === RGBA
    @test @inferred(coloralpha(RGBA)) === RGBA
    @test @inferred(coloralpha(ARGB)) === RGBA
    @test @inferred(coloralpha(RGB{N0f8})) === RGBA
    @test @inferred(coloralpha(RGBA{Float32})) === RGBA
    @test @inferred(coloralpha(ARGB{Float64})) === RGBA
    @test_throws MethodError coloralpha(RGB24)
    @test_throws MethodError coloralpha(ARGB32)

    @test @inferred(coloralpha(BGR)) === BGRA
    @test @inferred(coloralpha(BGR{N0f8})) === BGRA
    @test @inferred(coloralpha(XRGB)) === RGBA
    @test @inferred(coloralpha(RGBX)) === RGBA

    @test @inferred(coloralpha(Gray)) === GrayA
    @test @inferred(coloralpha(GrayA)) === GrayA
    @test @inferred(coloralpha(AGray)) === GrayA
    @test @inferred(coloralpha(Gray{N0f8})) === GrayA
    @test @inferred(coloralpha(GrayA{Float32})) === GrayA
    @test @inferred(coloralpha(AGray{Float64})) === GrayA
    @test_throws MethodError coloralpha(Gray24)
    @test_throws MethodError coloralpha(AGray32)

    @test @inferred(coloralpha(HSV)) === HSVA
    @test @inferred(coloralpha(HSVA)) === HSVA
    @test @inferred(coloralpha(AHSV)) === HSVA
    @test @inferred(coloralpha(HSV{Float16})) === HSVA
    @test @inferred(coloralpha(HSVA{Float32})) === HSVA
    @test @inferred(coloralpha(AHSV{Float64})) === HSVA

    @test_throws MethodError coloralpha(AbstractRGB)
    @test_throws MethodError coloralpha(N0f8)
end

@testset "alphacolor for types" begin
    @test @inferred(alphacolor(RGB)) === ARGB
    @test @inferred(alphacolor(RGBA)) === ARGB
    @test @inferred(alphacolor(ARGB)) === ARGB
    @test @inferred(alphacolor(RGB{N0f8})) === ARGB
    @test @inferred(alphacolor(RGBA{Float32})) === ARGB
    @test @inferred(alphacolor(ARGB{Float64})) === ARGB
    @test @inferred(alphacolor(RGB24)) === ARGB32
    @test @inferred(alphacolor(ARGB32)) === ARGB32

    @test @inferred(alphacolor(BGR)) === ABGR
    @test @inferred(alphacolor(BGR{N0f8})) === ABGR
    @test @inferred(alphacolor(XRGB)) === ARGB
    @test @inferred(alphacolor(RGBX)) === ARGB

    @test @inferred(alphacolor(Gray)) === AGray
    @test @inferred(alphacolor(GrayA)) === AGray
    @test @inferred(alphacolor(AGray)) === AGray
    @test @inferred(alphacolor(Gray{N0f8})) === AGray
    @test @inferred(alphacolor(GrayA{Float32})) === AGray
    @test @inferred(alphacolor(AGray{Float64})) === AGray
    @test @inferred(alphacolor(Gray24)) === AGray32
    @test @inferred(alphacolor(AGray32)) === AGray32

    @test @inferred(alphacolor(HSV)) === AHSV
    @test @inferred(alphacolor(HSVA)) === AHSV
    @test @inferred(alphacolor(AHSV)) === AHSV
    @test @inferred(alphacolor(HSV{Float16})) === AHSV
    @test @inferred(alphacolor(HSVA{Float32})) === AHSV
    @test @inferred(alphacolor(AHSV{Float64})) === AHSV

    @test_throws MethodError alphacolor(AbstractRGB)
    @test_throws MethodError alphacolor(N0f8)
end


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

@doc raw"""APIResourceList is a list of APIResource, it is used to expose the name of the resources supported in a specific group and version, and if the resource is namespaced.

    IoK8sApimachineryPkgApisMetaV1APIResourceList(;
        apiVersion=nothing,
        groupVersion=nothing,
        kind=nothing,
        resources=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - groupVersion::String : groupVersion is the group and version this APIResourceList is for.
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - resources::Vector{IoK8sApimachineryPkgApisMetaV1APIResource} : resources contains the name of the resources and if they are namespaced.
"""
mutable struct IoK8sApimachineryPkgApisMetaV1APIResourceList <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    groupVersion::Any # spec type: Union{ Nothing, String } # spec name: groupVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    resources::Any # spec type: Union{ Nothing, Vector{IoK8sApimachineryPkgApisMetaV1APIResource} } # spec name: resources

    function IoK8sApimachineryPkgApisMetaV1APIResourceList(;apiVersion=nothing, groupVersion=nothing, kind=nothing, resources=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResourceList, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResourceList, Symbol("groupVersion"), groupVersion)
        setfield!(o, Symbol("groupVersion"), groupVersion)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResourceList, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApimachineryPkgApisMetaV1APIResourceList, Symbol("resources"), resources)
        setfield!(o, Symbol("resources"), resources)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1APIResourceList

const _property_map_IoK8sApimachineryPkgApisMetaV1APIResourceList = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("groupVersion")=>Symbol("groupVersion"), Symbol("kind")=>Symbol("kind"), Symbol("resources")=>Symbol("resources"))
const _property_types_IoK8sApimachineryPkgApisMetaV1APIResourceList = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("groupVersion")=>"String", Symbol("kind")=>"String", Symbol("resources")=>"Vector{IoK8sApimachineryPkgApisMetaV1APIResource}")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1APIResourceList }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1APIResourceList))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1APIResourceList }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1APIResourceList[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1APIResourceList }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1APIResourceList[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1APIResourceList)
    (getproperty(o, Symbol("groupVersion")) === nothing) && (return false)
    (getproperty(o, Symbol("resources")) === nothing) && (return false)
    true
end

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


text:
function generate_sphere_system(n, m)
	x0 = zeros(n)

	Rs = rand(m) .+ 1
	centers = zeros(n, m)

	for i = 1:m
	    dir = randn(n)
	    normalize!(dir)
	    centers[:,i] = x0 + Rs[i]*dir
	end

	sphere(x, c, R) = dot(x, x) + dot(c, c) - 2*dot(x, c) - R^2

	function c!(cval, x)
	    for i = 1:m
	        cval[i] = sphere(x, view(centers, :, i), Rs[i])
	    end
	end

	function jac!(J, cval, x)
		fill!(J, 0.0)

	    for i = 1:m
	    	cval[i] = sphere(x, view(centers, :, i), Rs[i])
	    	J[i,:] .= 2.0.*(x .- view(centers, :, i))
	    end
	end

	return x0, c!, jac!
end

function generate_sin_system(n, m)
	x0 = zeros(n)

	function c!(cval, x)
	    for i = 1:m
	        cval[i] = x[2*i] - sin(x[2*i-1])
	    end
	end

	function jac!(J, cval, x)
		fill!(J, 0.0)

	    for i = 1:m
	    	cval[i] = x[2*i] - sin(x[2*i-1])
	    	J[i, 2*i] = 1.0
	    	J[i, 2*i-1] = -cos(x[2*i-1])
	    end
	end

	return x0, c!, jac!
end

@testset verbose=true "Retractions                 " begin

	n = 1000
	m = 100

	x0, c!, jac! = generate_sin_system(n, m)

	# do SVD
	cval = zeros(m)
	cval2 = zeros(m)
	J = zeros(m,n)
	jac!(J, cval, x0)

	U, S, Vt = svd(Array(J'))
	Vt = Array(Vt);

	nr = LFPSQP.NR(U, S, Vt, 1.0, 1000, LFPSQP.NRWork(m), false, LFPSQP.InequalityData())
	pp = LFPSQP.ProjPenalty(jac!, U, S, Vt, m, 0.01, 1.0, 100, 200, LFPSQP.ProjPenaltyWork(m, n, m, n), false, LFPSQP.InequalityDecomp(), LFPSQP.InequalityData())

	# generate random step and project onto tangent plane
	step = randn(n)
	step -= U*(U'*step)
	step .*= 5.0 / norm(step)

	xtilde = x0 + step
	xtilde_copy = copy(xtilde)
	xnew = zeros(n)

	c!(cval, xtilde)
	# @show norm(cval)

	@testset "Newton retraction" begin
		for tol in [1e-6, 1e-8]
			nr.tol = tol

			flag, i, _ = LFPSQP.retract!(cval, xnew, c!, xtilde, x0, nr)
			c!(cval2, xnew)

			# @show flag, i
			@test flag == 0 && norm(cval, Inf) < tol
			@test all(cval .== cval2)			# cval should be c! evaluated at xnew
			@test all(xtilde .== xtilde_copy)
			@test abs(dot(step, xnew - xtilde)) < 1e-6
		end
	end

	@testset "CG iteration" begin
		# solve (μI + J'J)x = b
		J = randn(m, n)
		p = zeros(n)
		z = zeros(n)
		tmp_m = zeros(m)
		tol = 1e-6
		maxiter = 100

		for μ in [1e-1, 1e-2, 1e-4]
			x = zeros(n) # initial guess
			b = randn(n)
			r = copy(b)

			M! = LFPSQP.no_precondition

			flag, i = LFPSQP.pcg!(μ, J, M!, x, r, p, z, tmp_m, tol, maxiter)

			@test flag == 0
			@test norm(r) < tol
			@test norm(μ.*x + J'*(J*x) - b) < tol

			# test exact preconditioning
			x = zeros(n) # initial guess
			b = randn(n)
			r = copy(b)

			exact_inverse = inv(μ*Matrix(I, n, n) + J'*J)
			M!(z, r) = z .= (μ*Matrix(I, n, n) + J'*J) \ r

			flag, i = LFPSQP.pcg!(μ, J, M!, x, r, p, z, tmp_m, tol, maxiter)

			@test flag == 0 && i == 1
			@test norm(r) < tol
			@test norm(μ.*x + J'*(J*x) - b) < tol
		end
	end

	@testset "Project penalty retraction" begin
		for tol in [1e-6, 1e-8, 1e-10]
			pp.tol = tol

			flag, i, pcg_i = LFPSQP.retract!(cval, xnew, c!, xtilde, x0, pp)
			c!(cval2, xnew)

			# @show flag, i, pcg_i
			@test flag == 0 && norm(cval, Inf) < tol
			@test all(cval .== cval2)			# cval should be c! evaluated at xnew
			@test all(xtilde .== xtilde_copy)
			@test norm(step) >= norm(xnew - x0) - tol
		end
	end
end


text:
"""
    anchor(key, [value = nothing])

Define a "manually steppable" point.
"""
anchor

use_anchors() = false

function enable_anchors()
    prev = use_anchors()
    @eval use_anchors() = true
    return prev
end

function disable_anchors()
    prev = use_anchors()
    @eval use_anchors() = false
    return prev
end

function anchor_impl(key, value)
    found = get(task_local_storage(), _ANCHOR_KEY, nothing)
    found === nothing && return
    task, on, enabled = found
    if enabled[] && task isa Task && on(key)
        yieldto(task, Some((key, value)))
    end
    return
end

function anchor(x, value = nothing)
    if use_anchors()
        anchor_impl(x, value)
    end
    return
end

struct AnchorError
    value::Any
end

struct AnchoredTask
    task::Task
    enabled::typeof(Ref(true))
end

always(_) = true

const _ANCHOR_KEY = :REAGENTS_ANCHOR_KEY

function withanchor(f, on = always)
    @nospecialize
    use_anchors() || error("anchors not enabled")
    enabled = Ref(true)
    parent = current_task()
    task = @task try
        task_local_storage(f, _ANCHOR_KEY, (parent, on, enabled))
    catch err
        enabled[] && yieldto(parent, AnchorError(err))
        rethrow()
    end
    return AnchoredTask(task, enabled)
end

function nextanchor!(a::AnchoredTask)
    a.enabled[] || error("`finish!(a)` already called")
    y = yieldto(a.task, nothing)
    if y isa Some
        return something(y)
    elseif y isa AnchorError
        schedule(a.task)
        wait(a.task)::Union{}
    else
        error("unexpected: $y")
    end
end

function finish!(a::AnchoredTask)
    a.enabled[] = false
    schedule(a.task)
    return fetch(a.task)
end


text:
module ModuleMathAliases
export √, ∑ , ⊗
include("math_aliases.jl")
end


text:
# simplified formula:	a = (500000-1000.b)/(1000-b), max b = 499
for b = 499:-1:1
    (a, r) = divrem(500000-1000*b, 1000-b)
    if r == 0
        println(a*b*(1000-a-b))
        break
    end
end

max_stars_repo_path: 009 - Special Pythagorean Triplet/009 - Special_Pythagorean_Triplet.jl

max_stars_repo_name: LanguageAdventures/Project-Euler-Julia

max_issues_repo_name: LanguageAdventures/Project-Euler-Julia

max_forks_repo_name: LanguageAdventures/Project-Euler-Julia


text:
## Dynamics
abstract type AbstractModel <: AbstractFunction end
output_dim(model::AbstractModel) = state_dim(model)

"""
    ContinuousDynamics

A dynamics model of the form:

```math
\\dot{x} = f(x,u)
```

where ``x`` is the `n`-dimensional state vector and ``u`` is the `m`-dimensional
control vector. The following methods should be defined on any model:

    state_dim(model)
    control_dim(model)

The output dimension `output_dim` is automatically defined to be equal to the
state dimension.

## Defining the dynamics
To define the dynamics using out-of-place (returning an SVector) evaluation,
the user must define one of the following:

    dynamics(model, z::AbstractKnotPoint)
    dynamics(model, x, u, t)
    dynamics(model, x, u)

The user shouldn't assume that the inputs are static arrays, although for best
computational performance these methods should be called by passing in static arrays.

To define the dynamics using in-place evaluation, the user must define one of the following:

    dynamics!(model, xdot, z::AbstractKnotPoint)
    dynamics!(model, xdot, x, u, t)
    dynamics!(model, xdot, x, u)

A user-defined Jacobian can be provided by defining one of the following methods

    jacobian!(model, J, xdot, z::AbstractKnotPoint)
    jacobian!(model, J, xdot, x, u, t)
    jacobian!(model, J, xdot, x, u)

## Non-Euclidean state vectors
By default, all elements of the state vector are assume to be in Euclidean space.
This assumption can be relaxed by defining a few extra functions. For non-Euclidean
state vectors, the Euclidean difference between two states, or the error state, is no
longer computed using pure subtraction. We can define a custom function for taking the
differences between 2 states vectors based on the type of the state vector.

We can query the type of the state vector using [`statevectortype(model)`](@ref),
which returns a [`StateVectorType`](@ref) trait (by default, [`EuclideanState`](@ref)).
After defining a new `StateVectorType` and the following methods (described in more
detail in the documentation for [`StateVectorType`](@ref)):

    state_diff!
    errstate_dim
    state_diff_jacobian!
    ∇state_diff_jacobian!

We can use dynamics with non-Euclidean state vectors. The most commonly encounter
case of this is when using 3D rotations. Support for states with 3D rotations is
provided by the [`RotationState`](@ref) `StateVectorType`.
"""
abstract type ContinuousDynamics <: AbstractModel end

"""
    dynamics(model, z::AbstractKnotPoint)
    dynamics(model, x, u, t)
    dynamics(model, x, u)

Evaluate the continuous time dynamics, returning the output ``\\dot{x}``. For best
performance, the output should usually be a `StaticArrays.SVector`. This method is
called when using the `StaticReturn` [`FunctionSignature`](@ref).

Calling `evaluate` on a [`ContinuousDynamics`](@ref) model will call this function.
"""
@inline dynamics(model::ContinuousDynamics, z::AbstractKnotPoint) =
    dynamics(model, state(z), control(z), time(z))
@inline dynamics(model::ContinuousDynamics, x, u, t) = dynamics(model, x, u)

"""
    dynamics!(model, xdot, z::AbstractKnotPoint)
    dynamics!(model, xdot, x, u, t)
    dynamics!(model, xdot, x, u)

Evaluate the continuous time dynamics, storing the output in `xdot`.
This method is called when using the `InPlace` [`FunctionSignature`](@ref).

Calling `evaluate!` on a [`ContinuousDynamics`](@ref) model will call this function.
"""
@inline dynamics!(model::ContinuousDynamics, xdot, z::AbstractKnotPoint) =
    dynamics!(model, xdot, state(z), control(z), time(z))
@inline dynamics!(model::ContinuousDynamics, ẋ, x, u, t) = dynamics!(model, ẋ, x, u)

"""
    dynamics!(sig, model, xdot, z::AbstractKnotPoint)

Evaluate the continuous time dynamics, storing the output in `xdot`, using the
[`FunctionSignature`](@ref) `sig` to determine which method to call.
"""
dynamics!(::InPlace, model::ContinuousDynamics, xdot, z::AbstractKnotPoint) =
    dynamics!(model, xdot, z)
dynamics!(::StaticReturn, model::ContinuousDynamics, xdot, z::AbstractKnotPoint) =
    xdot .= dynamics!(model, z)

@inline evaluate(model::ContinuousDynamics, x, u, p) =
    dynamics(model, x, u, p.t)
@inline evaluate!(model::ContinuousDynamics, ẋ, x, u, p) =
    dynamics!(model, ẋ, x, u, p.t)

@inline jacobian!(::FunctionSignature, ::UserDefined, model::ContinuousDynamics, J, ẋ, z) =
    jacobian!(model, J, ẋ, state(z), control(z), time(z))


text:
struct NNKolmogorov{C,O,S,E} <: NeuralNetDiffEqAlgorithm
    chain::C
    opt::O
    sdealg::S
    ensemblealg::E
end
NNKolmogorov(chain  ; opt=Flux.ADAM(0.1) , sdealg = EM() , ensemblealg = EnsembleThreads()) = NNKolmogorov(chain , opt , sdealg , ensemblealg)

function DiffEqBase.solve(
    prob::Union{KolmogorovPDEProblem,SDEProblem},
    alg::NNKolmogorov;
    abstol = 1f-6,
    verbose = false,
    maxiters = 300,
    trajectories = 1000,
    save_everystep = false,
    use_gpu = false,
    dt,
    dx,
    kwargs...
    )

    tspan = prob.tspan
    sigma = prob.g
    μ = prob.f
    noise_rate_prototype = prob.noise_rate_prototype
    if prob isa SDEProblem
        xspan = prob.kwargs.data.xspan
        d = prob.kwargs.data.d
        u0 = prob.u0
        phi(xi) = pdf(u0 , xi)
    else
        xspan = prob.xspan
        d = prob.d
        phi = prob.phi
    end
    ts = tspan[1]:dt:tspan[2]
    xs = xspan[1]:dx:xspan[2]
    N = size(ts)
    T = tspan[2]

    #hidden layer
    chain  = alg.chain
    opt    = alg.opt
    sdealg = alg.sdealg
    ensemblealg = alg.ensemblealg
    ps     = Flux.params(chain)
    xi     = rand(xs ,d ,trajectories)
    #Finding Solution to the SDE having initial condition xi. Y = Phi(S(X , T))
    sdeproblem = SDEProblem(μ,sigma,xi,tspan,noise_rate_prototype = noise_rate_prototype)
    function prob_func(prob,i,repeat)
      SDEProblem(prob.f , prob.g , xi[: , i] , prob.tspan ,noise_rate_prototype = prob.noise_rate_prototype)
    end
    output_func(sol,i) = (sol[end],false)
    ensembleprob = EnsembleProblem(sdeproblem , prob_func = prob_func , output_func = output_func)
    sim = solve(ensembleprob, sdealg, ensemblealg , dt=dt, trajectories=trajectories,adaptive=false)
    x_sde  = reshape([],d,0)
    # sol = solve(sdeproblem, sdealg ,dt=0.01 , save_everystep=false , kwargs...)
    # x_sde = sol[end]
    for u in sim.u
        x_sde = hcat(x_sde , u)
    end
    y = phi(x_sde)
    if use_gpu == true
        y = y |>gpu
        xi = xi |> gpu
    end
    data   = Iterators.repeated((xi , y), maxiters)
    if use_gpu == true
        data = data |>gpu
    end

    #MSE Loss Function
    loss(x , y) =Flux.mse(chain(x), y)

    cb = function ()
        l = loss(xi, y)
        verbose && println("Current loss is: $l")
        l < abstol && Flux.stop()
    end

    Flux.train!(loss, ps, data, opt; cb = cb)
    chainout = chain(xi)
    xi , chainout
 end #solve


text:
include("model.jl")
include(Knet.dir("data", "mnist.jl"))
using Debugger
using NPZ
using LinearAlgebra
Knet.seed!(1)

σ = .17
num_glimpses = 2
patch_size = 8
num_patches = 1
glimpse_scale = 2.
num_channels = 1
loc_hidden = 128
glimpse_hidden = 128
hidden_size = 256
num_classes = 10

ram = RAM(
    patch_size,
    num_patches,
    glimpse_scale,
    num_channels,
    loc_hidden,
    glimpse_hidden,
    σ,
    hidden_size,
    num_classes,
    num_glimpses,
)

dtrn, dtst = mnistdata(;batchsize=32)
x, y = first(dtrn)
B = length(y)
lt, ht = initstates(ram, B)
# ram.controller.h = ht
atype = artype(ram)
etype = eltype(ram)

fc_w = npzread("weights/baseliner.fc.weight.npy") .+ 0
fc_b = npzread("weights/baseliner.fc.bias.npy") .+ 0
ram.baseline_net.w = Param(atype{etype}(fc_w))
ram.baseline_net.b = Param(atype{etype}(fc_b))

fc_w = npzread("weights/classifier.fc.weight.npy") .+ 0
fc_b = npzread("weights/classifier.fc.bias.npy") .+ 0
ram.softmax_layer.w = Param(atype{etype}(fc_w))
ram.softmax_layer.b = Param(atype{etype}(fc_b))

fc_w = npzread("weights/locator.fc.weight.npy") .+ 0
fc_b = npzread("weights/locator.fc.bias.npy") .+ 0
ram.location_net.layer.w = Param(atype{etype}(fc_w))
ram.location_net.layer.b = Param(atype{etype}(fc_b))

fc_w = npzread("weights/sensor.fc1.weight.npy") .+ 0
fc_b = npzread("weights/sensor.fc1.bias.npy") .+ 0
ram.glimpse_net.fc1.w = Param(atype{etype}(fc_w))
ram.glimpse_net.fc1.b = Param(atype{etype}(fc_b))

fc_w = npzread("weights/sensor.fc2.weight.npy") .+ 0
fc_b = npzread("weights/sensor.fc2.bias.npy") .+ 0
ram.glimpse_net.fc2.w = Param(atype{etype}(fc_w))
ram.glimpse_net.fc2.b = Param(atype{etype}(fc_b))

fc_w = npzread("weights/sensor.fc3.weight.npy") .+ 0
fc_b = npzread("weights/sensor.fc3.bias.npy") .+ 0
ram.glimpse_net.fc3.w = Param(atype{etype}(fc_w))
ram.glimpse_net.fc3.b = Param(atype{etype}(fc_b))

fc_w = npzread("weights/sensor.fc4.weight.npy") .+ 0
fc_b = npzread("weights/sensor.fc4.bias.npy") .+ 0
ram.glimpse_net.fc4.w = Param(atype{etype}(fc_w))
ram.glimpse_net.fc4.b = Param(atype{etype}(fc_b))

fc_w = npzread("weights/rnn.h2h.weight.npy") .+ 0
fc_b = npzread("weights/rnn.h2h.bias.npy") .+ 0
ram.controller.h2h.w = Param(atype{etype}(fc_w))
ram.controller.h2h.b = Param(atype{etype}(fc_b))

fc_w = npzread("weights/rnn.i2h.weight.npy") .+ 0
fc_b = npzread("weights/rnn.i2h.bias.npy") .+ 0
ram.controller.i2h.w = Param(atype{etype}(fc_w))
ram.controller.i2h.b = Param(atype{etype}(fc_b))

x = permutedims(npzread("weights/x.npy"), (4,3,2,1))
y = npzread("weights/y.npy") .+1
g_noise = atype{etype}(npzread("weights/noise.npy")' .+ 0.0)
lt, ht = initstates(ram, 32)

loss(x, ygold) = sum(ram(x,ygold)[1:3])

max_stars_repo_name: ilkerkesen/RAM

max_issues_repo_name: ilkerkesen/RAM

max_forks_repo_name: ilkerkesen/RAM


text:
module uniquac

include("general.jl")
using Unitful
using NLsolve

∑(x)=sum(x)

τ(T::Quantity,Δu,i::Int,j::Int) = exp(-Δu[i,j]/(1.98721u"cal/(mol*K)"*T))
τ(T::Float64,Δu,i::Int,j::Int) = exp(-ustrip(u"cal/mol",Δu[i,j])/(1.98721*T))
V(x,r,i::Int,n::Int) = r[i]/∑(r[j]*x[j] for j=1:n)
F(x,q,i::Int,n::Int) = q[i]/∑(q[j]*x[j] for j=1:n)
γc(x,q,r,i::Int,n::Int) = 1-V(x,r,i,n)+log(V(x,r,i,n))-5*q[i]*(1-V(x,r,i,n)/F(x,q,i,n) +log(V(x,r,i,n)/F(x,q,i,n)))
γR(T,x,Δu,q,r,i::Int,n::Int) = q[i]*(1-log(∑(q[j]*x[j]*τ(T,Δu,j,i) for j=1:n)/∑(q[j]*x[j] for j=1:n))
-∑(q[j]*x[j]*τ(T,Δu,i,j)/∑(q[k]*x[k]*τ(T,Δu,k,j) for k=1:n) for j=1:n))
γ(T,x,Δu,q,r,i::Int,n::Int) = exp(γc(x,q,r,i,n) + γR(T,x,Δu,q,r,i,n))
P(T,x,CA,Δu,q,r,i::Int,n::Int) = x[i]*γ(T,x,Δu,q,r,i,n)*general.P_Antoine(CA,T,i)
Pₜₒₜ(T,x,CA,Δu,q,r,n::Int) = ∑(P(T,x,CA,Δu,q,r,i,n) for i=1:n)
y_calc(T,x,CA,Δu,q,r,i::Int,n::Int)=P(T,x,CA,Δu,q,r,i,n)/Pₜₒₜ(T,x,CA,Δu,q,r,n)

"""
    pxy(n,T,CA,Δu,q,r,xx;uni=u"Torr")

**Los campos de entrada son:**
- `n :: Int` Es el numero de componentes
- `T :: Float` Es la temperatura del sistema
- `CA :: n×3 Array{Float,2}` Constantes de Antoine para los `n` componentes
- `Δu :: n×n Array{Float,2}` Parametros de interacción entre los componentes en grados kelvin
- `q :: n-element Array{Float64,1}` área de Van der Waals
- `r :: n-element Array{Float64,1}` volumen de Van der Waals
- `xx :: StepRange ó Vector` Son los puntos en el liquido donde se buscara el equilibrio con el vapor
- `uni` Son las unidades en las que se desea el resultado, por defecto esta en Torr

**Salida:**

Regresa las fracciones de liquido `x`,
las fracciones del vapor `y` y la presión `P` en Torr (`m` es la longitud de `xx`).

- `x :: n×m Array{Float64,2}` liquido
- `y :: n×m Array{Float64,2}` vapor
- `P :: m-element Array{Float64,1}` presión
"""
function pxy(n::Int,T,CA::Array,Δu::Array,q::Vector,r::Vector,xx;uni=u"Torr")
    T = T|> u"K"
    xx = general.nc(xx...)
    m = size(xx,2)
    y = Array{Float64}(undef,n,m)
    P = Vector{Float64}(undef,m)*uni
    for i=1:m
        x = xx[:,i]
        for j=1:n
            y[j,i] = y_calc(T,x,CA,Δu,q,r,j,n)
        end
        P[i] = Pₜₒₜ(T,x,CA,Δu,q,r,n) |> uni
    end

    return xx,y,P

end

"""
    txy(n,Pobj,CA,Δu,q,r,xx;uni=u"K")

**Los campos de entrada son:**
- `n :: Int` Es el numero de componentes
- `Pobj :: Float` Es la presión del sistema
- `CA :: n×3 Array{Float,2}` Constantes de Antoine para los `n` componentes
- `Δu :: n×n Array{Float,2}` Parametros de interacción entre los componentes en grados kelvin
- `q :: n-element Array{Float64,1}` área de Van der Waals
- `r :: n-element Array{Float64,1}` volumen de Van der Waals
- `xx :: StepRange ó Vector` Son los puntos en el liquido donde se buscara el equilibrio con el vapor
- `uni` Son las unidades en las que se desea el resultado, por defecto esta en grados kelvin

**Salida:**

Regresa las fracciones de liquido `x`,
las fracciones del vapor `y` y la temperatura `T` en grados kelvin (`m` es la longitud de `xx`).

- `x :: n×m Array{Float64,2}` liquido
- `y :: n×m Array{Float64,2}` vapor
- `T :: m-element Array{Float64,1}` temperatura
"""
function txy(n::Int,Pobj,CA::Array,Δu::Array,q::Vector,r::Vector,xx;uni=u"K")
    Pobj = Pobj |> u"Torr"
    xx = general.nc(xx...)
    m = size(xx,2)
    y = Array{Float64}(undef,n,m)
    T = Vector{Float64}(undef,m)*u"K"
    Tᵢ = 300.0
    #temperatura inicial
    for i=1:n
        if xx[i,1]==1
            T[1] = general.T_Antoine(CA,Pobj,i)
            Tᵢ  = ustrip(u"K",T[1])
        end
    end
    for i = 1:m
        x = xx[:,i]
        P2ₜₒₜ(T) = Pₜₒₜ(T,x,CA,Δu,q,r,n) - ustrip(u"Torr",Pobj)
        #println(i)
        Pₛₒₗ = nlsolve(n_ary(P2ₜₒₜ),[Tᵢ])
        Tᵢ = Pₛₒₗ.zero[1]
        T[i] = Tᵢ*u"K"
        for j = 1:n
            y[j,i] = y_calc(T[i],x,CA,Δu,q,r,j,n)
        end
    end

    return xx,y,T .|> uni

end
end  # module

max_stars_repo_name: EmilioAlvizo/TermoQuimica

max_issues_repo_name: EmilioAlvizo/TermoQuimica

max_forks_repo_name: EmilioAlvizo/TermoQuimica


text:
function _draw_action_from_MvNormal!(a::Vector{Float64}, P::MvNormal)
    # this assumes it is 2x2

    μ::Vector{Float64} = P.μ
    Σ::Matrix{Float64} = P.Σ.mat

    a[1] = randn()*sqrt(Σ[1,1]) + μ[1]

    # compute conditional values
    μ₂ = μ[2] + Σ[1,2]*(a[1] - μ[1])/Σ[1,1]
    Σ₂ = Σ[2,2] - Σ[1,2]*Σ[1,2]/Σ[1,1]
    a[2] = randn()*sqrt(Σ₂) + μ₂

    a
end
function _get_e_and_denom(P::MvNormal, a::Vector{Float64})

    μ::Vector{Float64} = P.μ
    Σ::Matrix{Float64} = P.Σ.mat

    Δx = a[1] - μ[1]
    Δy = a[2] - μ[2]
    det = Σ[1,1]*Σ[2,2] - Σ[1,2]*Σ[1,2]

    e = exp(-0.5*((Σ[2,2]*Δx - Σ[1,2]*Δy)*Δx + (Σ[1,1]*Δy - Σ[1,2]*Δx)*Δy) / det)
    denom = sqrt(4*π*π*det)

    (e,denom)
end
function _pdf_of_MvNormal(P::MvNormal{PDMats.PDMat{Float64,Matrix{Float64}},Vector{Float64}}, a::Vector{Float64})
    # this assumes it is 2x2

    e, denom = _get_e_and_denom(P, a)
    e / denom
end
function _logpdf_of_MvNormal(P::MvNormal{PDMats.PDMat{Float64,Matrix{Float64}},Vector{Float64}}, a::Vector{Float64})
    # this assumes it is 2x2

    e, denom = _get_e_and_denom(P, a)
    log(e) - log(denom)
end


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

"""
    InpaintAlgo

An inpainting algorithm with given parameters.
"""
abstract type InpaintAlgo end

"""
    inpaint(img, mask, algo)

Inpaint `img` on pixels marked as `true` in `mask` using
algorithm `algo`.
"""
function inpaint(img::AbstractArray{T,N}, mask::AbstractArray{Bool,N},
                 algo::InpaintAlgo) where {N,T}
  # sanity checks
  @assert size(img) == size(mask) "image and mask must have same size"

  # dispatch appropriate implementation
  inpaint_impl(img, mask, algo)
end

#------------------
# IMPLEMENTATIONS
#------------------
include("criminisi.jl")


text:
module SimplePosetAlgorithms

using SimpleGraphs
using SimplePosets
using SimpleGraphAlgorithms

export max_chain, max_antichain, width

"""
`max_chain(P)` returns a largest set of pairwise *comparable* elements
in the `SimplePoset`.
"""
function max_chain(P::SimplePoset)
    return max_clique(ComparabilityGraph(P))
end

"""
`max_antichain(P)` returns a largest set of pairwise *incomparable*
elements in the `SimplePoset`.
"""
function max_antichain(P::SimplePoset)
    return max_indep_set(ComparabilityGraph(P))
end

"""
`width(P)` gives the size of a largest antichain in the poset `P`.
"""
width(P::SimplePoset) = length(max_antichain(P))

include("realizer.jl")

end # of module SimplePosetAlgorithms


text:
using ControlSystems, RobustAndOptimalControl
using Plots
using LinearAlgebra
"""
This is a simple SISO example with integrator dynamics corresponding to the
quad tank process in the lab.

The example can be set to visualize and save plots using the variables
  makeplots - true/false (true if plots are to be generated, false for testing)
  SavePlots - true/false (true if plots are to be saved, false for testing)
"""
makeplots = true

# Define the proces parameters
k1, k2, kc, g = 3.33, 3.35, 0.5, 981
A1, A3, A2, A4 = 28, 28, 32, 32
a1, a3, a2, a4= 0.071, 0.071, 0.057, 0.057
h01, h02, h03, h04 = 12.4, 12.7, 1.8, 1.4
T1, T2 = (A1/a1)*sqrt(2*h01/g), (A2/a2)*sqrt(2*h02/g)
T3, T4 = (A3/a3)*sqrt(2*h03/g), (A4/a4)*sqrt(2*h04/g)
c1, c2 = (T1*k1*kc/A1), (T2*k2*kc/A2)
γ1, γ2 = 0.7, 0.6

# Define the process dynamics
A = [-1/T1     0 A3/(A1*T3)          0;
     0     -1/T2          0 A4/(A2*T4);
     0         0      -1/T3          0;
     0         0          0      -1/T4];
B = [γ1*k1/A1     0;
     0                γ2*k2/A2;
     0                (1-γ2)*k2/A3;
     (1-γ1)*k1/A4 0              ];
C = [kc 0 0 0;
     0 kc 0 0];
D = zeros(2,2)
G = ss(A,B,C,D);

# Sensitivity weight function
WS = makeweight(100, (0.1, 1), 1/10) * I(2)

# Output sensitivity weight function
# WUelement = 5*tf([1,1],[0.1,1])
WU = tf(0.01) .* I(2)

# Complementary sensitivity weight function
WT = makeweight(1/10, (0.1, 1), 10) * I(2)

# Form augmented P dynamics in state-space
P = hinfpartition(G, WS, WU, WT)

# Check that the assumptions are satisfied
flag = hinfassumptions(P)

# Synthesize the H-infinity optimal controller
flag, C, γ = hinfsynthesize(P)

Pcl, S, CS, T = hinfsignals(P, G, C)

## Plot the specifications
# TODO figure out why I get segmentation errors when using ss instead of tf for
# the weighting functions, makes no sense at all
if makeplots
  specificationplot([S, CS, T], [WS[1,1], 0.01, WT[1,1]], γ)

## Plot the closed loop gain from w to z
# TODO figure out why the legends don't seem to work in this case

  specificationplot(Pcl, γ; s_labels=["\$\\sigma(P_{cl}(j\\omega))\$"], w_labels=["\$\\gamma\$"])

  times = [i for i in range(0, stop=300, length=10000)]
  plot(step(T, times))
end