text:
using Mosek
using LinearAlgebra

function ind2sub_replacement(a, b)
    indices = CartesianIndices(a)[b]
    subk = Array{Int32}(undef, length(indices))
    subl = Array{Int32}(undef, length(indices))
    for (i, ind) in enumerate(indices)
        subk[i], subl[i] = Tuple(ind)
    end
    return (subk, subl)
end

# The 'prob' structure specifies the moment problem
#
# minimize    prob.obj'*y
# subject to  prob.mom[k]*y is PSD,  k=1,...,length(prob.mom)
#             prob.eq[k]*y = 0,      k=1,...,length(prob.eq)
#             y[1] = 1
#
# We formulate the dual problem for MOSEK
#
# maximize     t
# subject to   sum_j dot(prob.mom[j][:,1], Xj) + sum_k dot(prb.eq[k][:,1], Zk)) = prob.obj[1] - t
#              sum_j dot(prob.mom[j][:,i], Xj) + sum_k dot(prb.eq[k][:,1], Zk)) = prob.obj[i], i=2,...,length(prob.obj)
#              Xj is PSD, j=1,...,length(prob.mom)
#              Zk is symmetric but free,  k=1,...,length(prob.eq)
#
function solve_mosek(prob::MomentProb; tolrelgap=1e-10, showlog=true)

    printstream(msg::AbstractString) = print(msg)

    # Create a task object and attach log stream printer
    task = maketask()
    if showlog  putstreamfunc(task,MSK_STREAM_LOG,printstream)  end

    # The momemt problem is in dual form, so we dualize it
    numcon = length(prob.obj)
    numbarvar = length(prob.mom)
    barvardim = Int[ sqrt(size(prob.mom[k],1)) for k=1:numbarvar ]

    eqdim = Int[ sqrt(size(prob.eq[k],1)) for k=1:length(prob.eq) ]
    eqidx = Array{Int}(undef, length(prob.eq)+1)
    eqidx[1] = 0
    for k=1:length(prob.eq)
        eqidx[k+1] = eqidx[k] + eqdim[k]*(eqdim[k]+1)/2
    end
    numvar = 1 + eqidx[ end ]

    # add free variables from equality constraints
    appendvars(task, Int32(numvar))

    putvarboundslice(task, 1, numvar+1,
                     [ MSK_BK_FR  for i in 1:numvar ],
                     [ -Inf       for i in 1:numvar ],
                     [ +Inf       for i in 1:numvar ])

    putcj(task, 1, 1.0)

    # Append matrix variables of sizes in 'BARVARDIM'.
    appendbarvars(task, barvardim)

    # Add constraints
    numconst = 1
    appendcons(task, 1)
    putaij(task, 1, 1, 1.0)

    I = Array{Int}(undef,0)
    for i=1:numcon

        added_const = ( i == 1);

        for j=1:numbarvar
            nj = Int64(barvardim[j])
            k1, k2 = prob.mom[j].colptr[i], prob.mom[j].colptr[i+1]-1
            if k2 >= k1

                if !added_const
                    appendcons(task, 1)
                    added_const = true
                end

                subk, subl = ind2sub_replacement((nj,nj), prob.mom[j].rowval[k1:k2])

		        trilidx = subk .>= subl
                aij = appendsparsesymmat(task, nj, subk[trilidx], subl[trilidx], prob.mom[j].nzval[k1:k2][trilidx])
                putbaraij(task, numconst, j, [aij], [1.0])
            end
        end

        for j=1:length(prob.eq)
            k1, k2 = prob.eq[j].colptr[i], prob.eq[j].colptr[i+1]-1
            if k2 >= k1
                if !added_const
                    appendcons(task, 1)
                    added_const = true
                end

                subk, subl = ind2sub_replacement((eqdim[j], eqdim[j]), prob.eq[j].rowval[k1:k2])
                trilidx = subk .>= subl
                subj = trilind( prob.eq[j].rowval[k1:k2][trilidx], eqdim[j] ) .+ (eqidx[j] + 1)

                putaijlist(task, numconst*ones(Int, length(subj)), subj, prob.eq[j].nzval[k1:k2][trilidx])
#                 subj = trilind( prob.eq[j].rowval[k1:k2], eqdim[j] ) + eqidx[j] + 1
#                 putaijlist(task, numconst*ones(Int, length(subj)), subj, prob.eq[j].nzval[k1:k2])
            end
        end

        if added_const
            push!(I, i)
            putconbound(task, numconst, MSK_BK_FX, prob.obj[i], prob.obj[i])
            numconst = numconst + 1
        end
    end

    # Input the objective sense (minimize/maximize)
    putobjsense(task,MSK_OBJECTIVE_SENSE_MAXIMIZE)

    putparam(task, "MSK_IPAR_NUM_THREADS", "4")
    putparam(task, "MSK_DPAR_INTPNT_CO_TOL_REL_GAP", string(tolrelgap))

    # Write .task file
    writetask(task, "polyopt.txt")

    # Solve the problem and print summary
    optimize(task)
    solutionsummary(task,MSK_STREAM_MSG)

    # Get status information about the solution
    solsta = getsolsta(task,MSK_SOL_ITR)

    X = [ symm(getbarxj(task, MSK_SOL_ITR, j), Int(sqrt(size(prob.mom[j],1)))) for j=1:length(prob.mom) ]
    t = getxxslice(task, MSK_SOL_ITR, 1, 2)[1]
    Z = [ symm(getxxslice(task, MSK_SOL_ITR, 2+eqidx[k], 2+eqidx[k+1]), eqdim[k]) for k=1:length(prob.eq) ]
    Z = [ 0.5*(Zk + diagm(diag(Zk))) for Zk = Z ]

    y = gety(task, MSK_SOL_ITR)
    if length(y) < length(prob.obj) y = sparsevec(I, y) end

    if solsta == MSK_SOL_STA_OPTIMAL
        return (X, Z, t, y, "Optimal")
#    elseif solsta == MSK_SOL_STA_NEAR_OPTIMAL
#        return (X, Z, t, y, "Near optimal")
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        return (X, Z, t, y, "Dual infeasibility")
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        return (X, Z, t, y, "Primal infeasibility")
#    elseif solsta == MSK_SOL_STA_NEAR_DUAL_INFEAS_CER
#        return (X, Z, t, y, "Near dual infeasibility")
#    elseif solsta == MSK_SOL_STA_NEAR_PRIM_INFEAS_CER
#        return (X, Z, t, y, "Near primal infeasibility")
    elseif solsta == MSK_SOL_STA_UNKNOWN
        return (X, Z, t, y, "Unknown")
    else
        error("Other solution status")
    end
end

trilind(k::Vector{Int}, n::Int) = Int[i + (j-1)*(n-1) - (j-1)*(j-2)/2 for (i,j) = zip(ind2sub_replacement((n, n),k)...) ]

function solve_mosek_no_elim(prob::BSOSProb; tolrelgap=1e-10, showlog=true)

    printstream(msg::AbstractString) = print(msg)

    # Create a task object and attach log stream printer
    task = maketask()
    if showlog  putstreamfunc(task,MSK_STREAM_LOG,printstream)  end

    f, Al, El, As = prob.obj, prob.Al, prob.El, prob.As

    m = length(Al)
    diml = [ size(a,1) for a in Al ]
    dimf = [ size(a,2) for a in Al ]
    dimX = Int[ round(Int, sqrt(size(as,1))) for as in As ]

    # variables are indexed as (l1,...,lm,f1,...,fm,t)
    idx_l = zeros(Int, length(diml)+1)
    idx_l[1] = 1
    for j=1:length(diml)
        idx_l[j+1] = idx_l[j] + diml[j]
    end

    idx_f = zeros(Int, length(dimf)+1)
    idx_f[1] = idx_l[end]
    for j=1:length(dimf)
        idx_f[j+1] = idx_f[j] + dimf[j]
    end

    l = 1
    for k=1:length(El)
        l += nnz(El[k])
    end
    const_idx = Array{Int,1}(undef, l)
    const_idx[1] = 1
    l = 1
    for k=1:length(El)
        const_idx[l + (1:nnz(El[k]))] = El[k].rowval
        l += nnz(El[k])
    end
    const_idx = unique(const_idx)
#
#     const_idx = Int[ 1 ]
#     for i=1:length(Al)
#         push!(const_idx, prob.El[i].rowval...)
#     end
#     const_idx = unique(const_idx)
    const_map = Dict(zip(const_idx, collect(1:length(const_idx))))

    numvar = sum(diml) + sum(dimf) + 1
    numcon = sum(dimf) + length(const_idx)
    numbarvar = m

    appendvars(task, numvar)

    # bounds on lj
    offs = 0
    for l=1:length(diml)
        for j=1:diml[l]
            if prob.lb[l][j] == -Inf
                putvarbound(task, offs+j, MSK_BK_FR, -Inf, Inf)
                #println("VARBOUND($(offs+j)): FREE")
            else
                putvarbound(task, offs+j, MSK_BK_LO, 0.0, Inf)
                #println("VARBOUND($(offs+j)): LOWER 0.0")
            end
        end
        offs += diml[l]
    end

    # fj and t free
    for j=sum(diml)+(1:sum(dimf)+1)
        putvarbound(task, j, MSK_BK_FR, -Inf, Inf)
        #println("VARBOUND($(j)): FREE")
    end

    putcj(task, numvar, 1.0)

    # Append matrix variables
    appendbarvars(task, dimX)

    # Add constraints
    appendcons(task, numcon)

    idx_const = 1

    for j=1:m
        #println("Block $(j)")
        for i=1:dimf[j]
            k1, k2 = Al[j].colptr[i], Al[j].colptr[i+1]-1
            subj = [idx_l[j]-1 + Al[j].rowval[k1:k2]; idx_f[j]-1 + i]
            val  = [Al[j].nzval[k1:k2]; -1.0]
            putarow(task, idx_const, subj, val)
            #println("CONSTRAINT A($(idx_const)): $(subj), $(val)")

            k1, k2 = As[j].colptr[i], As[j].colptr[i+1]-1
            if k2>=k1
                subk, subl = ind2sub_replacement( (dimX[j], dimX[j]), As[j].rowval[k1:k2] )
                I = subk .>= subl
                aij = appendsparsesymmat(task, dimX[j], subk[I], subl[I], As[j].nzval[k1:k2][I])
                #println("BARAIJ($(idx_const+i-1),$(j)): $(subk[I]), $(subl[I]), $(As[j].nzval[k1:k2][I])")
                putbaraij(task, idx_const, j, [aij], [1.0])
            end

            putconbound(task, idx_const, MSK_BK_FX, 0.0, 0.0)
            idx_const += 1
        end

    end

    putaij(task, idx_const, numvar, 1.0)
    for k=1:length(El)
        for (j,i) in enumerate(prob.El[k].rowval)
            #println("putaij: $(idx_const-1+const_map[i]), $(idx_f[k]-1+j), $(prob.El[k].nzval[j])")
            putaij(task, idx_const-1+const_map[i], idx_f[k]-1+j, prob.El[k].nzval[j])
        end
    end

    for j=idx_const:numcon
        putconbound(task, j, MSK_BK_FX, 0.0, 0.0)
    end

    for (j,i) in enumerate(f.nzind)
        #println("putconbound: $(idx_const-1+const_map[i]), $(f.nzval[j])")
        putconbound(task, idx_const-1+const_map[i], MSK_BK_FX, f.nzval[j], f.nzval[j])
    end

    # Input the objective sense (minimize/maximize)
    putobjsense(task,MSK_OBJECTIVE_SENSE_MAXIMIZE)

    putparam(task, "MSK_IPAR_NUM_THREADS", "4")
    putparam(task, "MSK_DPAR_INTPNT_CO_TOL_REL_GAP", string(tolrelgap))

    # Write .task file
    writetask(task, "polyopt.txt")

    # Solve the problem and print summary
    optimize(task)
    solutionsummary(task,MSK_STREAM_MSG)

    # Get status information about the solution
    solsta = getsolsta(task,MSK_SOL_ITR)

    X = [ symm(getbarxj(task, MSK_SOL_ITR, j), Int(sqrt(size(As[j],1)))) for j=1:m ]
    l = [ getxxslice(task, MSK_SOL_ITR, idx_l[j], idx_l[j+1]) for j=1:m ]
    f = [ getxxslice(task, MSK_SOL_ITR, idx_f[j], idx_f[j+1]) for j=1:m ]
    t = getxxslice(task, MSK_SOL_ITR, numvar, numvar+1)[1]
    y = gety(task, MSK_SOL_ITR)
    if solsta == MSK_SOL_STA_OPTIMAL
        return (X, t, l, f, y, "Optimal")
#    elseif solsta == MSK_SOL_STA_NEAR_OPTIMAL
#        return (X, t, l, f, y, "Near optimal")
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        return (X, t, l, f, y, "Dual infeasibility")
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        return (X, t, l, f, y, "Primal infeasibility")
#    elseif solsta == MSK_SOL_STA_NEAR_DUAL_INFEAS_CER
#        return (X, t, l, f, y, "Near dual infeasibility")
#    elseif solsta == MSK_SOL_STA_NEAR_PRIM_INFEAS_CER
#        return (X, t, l, f, y, "Near primal infeasibility")
    elseif solsta == MSK_SOL_STA_UNKNOWN
        return (X, t, l, f, y, "Unknown")
    else
        error("Other solution status")
    end
end

function solve_mosek(prob::BSOSProb; tolrelgap=1e-10, showlog=true)

    function getrow(Ak, Ek, j)
        sub, val = Int[], Float64[]

        # i = findfirst(Ek.rowval, j)
		i = findfirst(isequal(j),Ek.rowval)
        if i>0
            k1, k2 = Ak.colptr[i], Ak.colptr[i+1]-1
            sub = view(Ak.rowval, k1:k2)
            val = view(Ak.nzval, k1:k2)
        else
            sub, val = Int[], Float64[]
        end
        sub, val
    end

    function constraint_size(A, E, j)
        sz = (j == 1 ? 1 : 0)
        for k=1:length(A)
            sz += length(getrow(A[k], E[k], j)[1])
        end
        sz
    end

    printstream(msg::AbstractString) = print(msg)

    # Create a task object and attach log stream printer
    task = maketask()
    if showlog  putstreamfunc(task,MSK_STREAM_LOG,printstream)  end

    Al, El, As = prob.Al, prob.El, prob.As

    m = length(Al)
    diml = [ size(a,1) for a in Al ]
    dimX = Int[ round(Int, sqrt(size(as,1))) for as in As ]

    # variables are indexed as (l1,...,lm,t)
    idx_l = zeros(Int, length(diml)+1)
    idx_l[1] = 1
    for j=1:length(diml)
        idx_l[j+1] = idx_l[j] + diml[j]
    end

    l = 1
    for k=1:length(El)
        l += nnz(El[k])
    end
    const_idx = Array{Int,1}(undef, l)
    const_idx[1] = 1
    l = 1
    for k=1:length(El)
        const_idx[l .+ (1:nnz(El[k]))] = El[k].rowval
        l += nnz(El[k])
    end

    const_idx = sort(unique(const_idx))

    nonzero_markers = falses(length(const_idx))
    for (i,k)=enumerate(const_idx)

        if (nonzero_markers[i]) continue end

        if i==1
            nonzero_markers[i] = true
            continue
        end

        for j=1:m
            subj, valj = getrow(Al[j], El[j], k)
            if length(subj) > 0
                nonzero_markers[i] = true
                break
            end
        end

        if (nonzero_markers[i]) continue end

        for j=1:m
            subj, valj = getrow(As[j], El[j], k)
            if length(subj) >0
                nonzero_markers[i] = true
                break
            end
        end
    end
    const_idx = const_idx[ nonzero_markers ]

    numvar = sum(diml) + 1
    numcon = length(const_idx)
    numbarvar = m

    appendvars(task, numvar)

    # bounds on lj
    offs = 0
    for l=1:m
        for j=1:diml[l]
            if prob.lb[l][j] == -Inf
                putvarbound(task, offs+j, MSK_BK_FR, -Inf, Inf)
                #println("VARBOUND($(offs+j)): FREE")
            else
                putvarbound(task, offs+j, MSK_BK_LO, 0.0, Inf)
                #println("VARBOUND($(offs+j)): LOWER 0.0")
            end
        end
        offs += diml[l]
    end

    # t free
    putvarbound(task, numvar, MSK_BK_FR, -Inf, Inf)

    putcj(task, numvar, 1.0)

    # Append matrix variables
    appendbarvars(task, dimX)

    # Add constraints
    appendcons(task, numcon)

    f = prob.obj[const_idx]

    for (i,k)=enumerate(const_idx)
        offs = 0
        for j=1:m
            subj, valj = getrow(Al[j], El[j], k)
            if length(subj) > 0
                for r=1:length(subj)
                    putaij(task, i, subj[r]+offs, valj[r])
                end
            end
            offs += size(Al[j],1)
        end
        if i==1
            putaij(task, i, numvar, 1.0)
        end
        putconbound(task, i, MSK_BK_FX, f[i], f[i])

        for j=1:m
            subj, valj = getrow(As[j], El[j], k)
            if length(subj) >0
                subk, subl = ind2sub_replacement( (dimX[j], dimX[j]), subj )
                I = subk .>= subl
                aij = appendsparsesymmat(task, dimX[j], subk[I], subl[I], valj[I])
                putbaraij(task, i, j, [aij], [1.0])
            end
        end
    end

    # Input the objective sense (minimize/maximize)
    putobjsense(task,MSK_OBJECTIVE_SENSE_MAXIMIZE)

    putparam(task, "MSK_IPAR_NUM_THREADS", "4")
    putparam(task, "MSK_DPAR_INTPNT_CO_TOL_REL_GAP", string(tolrelgap))

    # Write .task file
    writetask(task, "polyopt.txt")

    # Solve the problem and print summary
    optimize(task)
    solutionsummary(task,MSK_STREAM_MSG)

    # Get status information about the solution
    solsta = getsolsta(task,MSK_SOL_ITR)

    X = [ symm(getbarxj(task, MSK_SOL_ITR, j), Int(sqrt(size(As[j],1)))) for j=1:m ]
    l = [ getxxslice(task, MSK_SOL_ITR, idx_l[j], idx_l[j+1]) for j=1:m ]
    t = getxxslice(task, MSK_SOL_ITR, numvar, numvar+1)[1]
    y = sparsevec(const_idx, gety(task, MSK_SOL_ITR), length(prob.obj))
    if solsta == MSK_SOL_STA_OPTIMAL
        return (X, t, l, y, "Optimal")
    elseif solsta == MSK_SOL_STA_NEAR_OPTIMAL
        return (X, t, l,  y, "Near optimal")
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        return (X, t, l,  y, "Dual infeasibility")
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        return (X, t, l,  y, "Primal infeasibility")
    elseif solsta == MSK_SOL_STA_NEAR_DUAL_INFEAS_CER
        return (X, t, l,  y, "Near dual infeasibility")
    elseif solsta == MSK_SOL_STA_NEAR_PRIM_INFEAS_CER
        return (X, t, l,  y, "Near primal infeasibility")
    elseif solsta == MSK_SOL_STA_UNKNOWN
        return (X, t, l,  y, "Unknown")
    else
        error("Other solution status")
    end
end

function symm(x::Array{T,1}, n::Int) where {T<:Number}
    X = zeros(n,n)
    k = 0
    for j=1:n
        X[j:n,j] = x[k .+ (1:n-j+1)]
        k += n-j+1
    end

    Matrix(Symmetric(X,:L))
end

max_stars_repo_name: shudianzhao/Polyopt

max_issues_repo_name: shudianzhao/Polyopt

max_forks_repo_name: shudianzhao/Polyopt


text:
using StructTypes

export PCResult

mutable struct PCAtomInt
    aid::Int64
    value::Int64

    PCAtomInt() = new()
end
StructTypes.StructType(::Type{PCAtomInt}) = StructTypes.Mutable()

mutable struct PCAtomString
    aid::Int64
    value::String

    PCAtomString() = new()
end
StructTypes.StructType(::Type{PCAtomString}) = StructTypes.Mutable()

@enum(PCUrnDataType,
    STRING = 1,      # String                             [maps to a VisibleString]
    STRINGLIST,    # List of Strings                    [maps to VisibleString list]
    INT,           # 32-Bit Signed Integer              [maps to an INTEGER]
    INTVEC,        # Vector of 32-Bit Signed Integer    [maps to INTEGER vector]
    UINT,          # 32-Bit Unsigned Integer            [maps to an INTEGER]
    UINTVEC,       # Vector of 32-Bit Unsigned Integer  [maps to INTEGER vector]
    DOUBLE,        # 64-Bit Float                       [maps to a REAL]
    DOUBLEVEC,     # Vector of Double                   [maps to REAL vector]
    BOOL,          # Boolean or Binary value            [maps to a BOOLEAN]
    BOOLVEC,       # Boolean Vector                     [maps to BOOLEAN vector]
    UINT64,        # 64-Bit Unsigned Integer (Hex form) [maps to a VisibleString]
    BINARY,        # Binary Data Blob                   [maps to an OCTET STRING]
    URL,           # URL                                [maps to a VisibleString]
    UNICODE,       # UniCode String                     [maps to a VisibleString]
    DATE,          # ISO8601 Date                       [maps to a Date]
    FINGERPRINT,   # Binary Fingerprint (Gzip'ped bit   [maps to an OCTET STRING]
                   #   list w/ 4-Byte prefix denoting bit list length)
    UNKNOWN = 255    # Unknown Data Type               [maps to a set of VisibleString]
)
StructTypes.StructType(::Type{PCUrnDataType}) = StructTypes.NumberType()
StructTypes.numbertype(::Type{PCUrnDataType}) = UInt8

mutable struct PCUrn
    label::String                          # Generic Name or Label for Display  [e.g., "Log P"]
    name::Union{Nothing,String}            # Qualified Name  [e.g., "XlogP"]
    datatype::Union{Nothing,PCUrnDataType} # Specific Data Type of Value  [e.g., binary]
    parameters::Union{Nothing,String}      # Implementation Parameter  [e.g., "metal=0"]
    implementation::Union{Nothing,String}  # Implementation Name  [e.g., "E_XlogP"]
    version::Union{Nothing,String}         # Implementation Version  [e.g., "3.317"]
    software::Union{Nothing,String}        # Implementation Software  [e.g., "Cactvs"]
    source::Union{Nothing,String}          # Implementation Organization  [e.g., "xemistry.com"]
    release::Union{Nothing,String}         # NCBI Implementation Release  [e.g., "10.25.2005"]

    PCUrn() = new()
end
StructTypes.StructType(::Type{PCUrn}) = StructTypes.Mutable()

mutable struct PCDateStd
    year::UInt32                   # full year (including 1900)
    month::Union{Nothing,UInt8}    # month (1-12)
    day::Union{Nothing,UInt8}      # day of month (1-31)
    season::Union{Nothing,String}  # for "spring", "may-june", etc
    hour::Union{Nothing,UInt8}     # hour of day (0-23)
    minute::Union{Nothing,UInt8}   # minute of hour (0-59)
    second::Union{Nothing,UInt8}   # second of minute (0-59)

    PCDateStd() = new()
end
StructTypes.StructType(::Type{PCDateStd}) = StructTypes.Mutable()

mutable struct PCDate
    str::Union{Nothing,String}    # for those unparsed dates
    std::Union{Nothing,PCDateStd} # use this if you can

    PCDate() = new()
end
StructTypes.StructType(::Type{PCDate}) = StructTypes.Mutable()

mutable struct PCInfoDataValue
    bval::Bool            # Boolean or Binary
    bvec::Vector{Bool}    # Boolean Vector
    ival::Int64           # Integer (signed or unsigned)
    ivec::Vector{Int64}   # Integer Vector
    fval::Float64         # Float or Double
    fvec::Vector{Float64} # Double Vector
    sval::String          # String
    slist::Vector{String} # List of Strings
    date::PCDate          # Date
    binary::Vector{UInt8} # Binary Data
    bitlist::Vector{Bool} # Bit List (specialized version of Boolean vector)

    PCInfoDataValue() = new()
end
StructTypes.StructType(::Type{PCInfoDataValue}) = StructTypes.Mutable()

mutable struct PCInfoData
    urn::PCUrn #   Universal Resource Name  [for Value Qualification]
    value::PCInfoDataValue

    PCInfoData() = new()
end
StructTypes.StructType(::Type{PCInfoData}) = StructTypes.Mutable()

mutable struct PCCount
    heavy_atom::UInt32        # Total count of non-Hydrogen (Heavy) Atoms

    # StereoChemistry Counts
    atom_chiral::UInt32       # Total count of (SP3) Chiral Atoms
    atom_chiral_def::UInt32   # Total count of Defined (SP3) Chiral Atoms
    atom_chiral_undef::UInt32 # Total count of Undefined (SP3) Chiral Atoms
    bond_chiral::UInt32       # Total count of (SP2) Chiral Bonds
    bond_chiral_def::UInt32   # Total count of (SP2) Defined Chiral Bonds
    bond_chiral_undef::UInt32 # Total count of (SP2) Undefined Chiral Bonds

    # Isotopic Counts
    isotope_atom::UInt32      # Total count of Atoms with Isotopic Information

   # Discrete Structure Counts
    covalent_unit::UInt32     # Total count of covalently_bonded units in the record
    tautomers::UInt32         # Number of possible tautomers (Max. 999)

    PCCount() = new()
end
StructTypes.StructType(::Type{PCCount}) = StructTypes.Mutable()

# we need to map the json names, which contain -, to _-using symbols
count_fieldnames = []
StructTypes.foreachfield((i, name, FT, v) -> push!(count_fieldnames, (name, Symbol(replace(String(name), "_" => "-")))), PCCount())
StructTypes.names(::Type{PCCount}) = Tuple(count_fieldnames)

@enum(PCBondAnnotation,
    BA_CROSSED = 1,    # Double Bond that can be both Cis/Trans
    BA_DASHED,       # Hydrogen-Bond (3D Only?)
    BA_WAVY,         # Unknown Stereochemistry
    BA_DOTTED,       # Complex/Fractional
    BA_WEDGE_UP,     # Above-Plane
    BA_WEDGE_DOWN,   # Below-Plane
    BA_ARROW,        # Dative
    BA_AROMATIC,     # Aromatic
    BA_RESONANCE,    # Resonance
    BA_BOLD,         # Fat Bond (Non-Specific User Interpreted Information)
    BA_FISCHER,      # Interpret Bond Stereo using Fischer Conventions
    BA_CLOSECONTACT, # Identification of Atom-Atom Close Contacts (3D Only)
    BA_UNKNOWN = 255   # Unspecified or Unknown Atom-Atom Annotation
)

StructTypes.StructType(::Type{PCBondAnnotation}) = StructTypes.NumberType()
StructTypes.numbertype(::Type{PCBondAnnotation}) = UInt8

mutable struct PCDrawAnnotations
    annotation::Vector{PCBondAnnotation}
    aid1::Vector{Int64}
    aid2::Vector{Int64}

    PCDrawAnnotations() = new()
end
StructTypes.StructType(::Type{PCDrawAnnotations}) = StructTypes.Mutable()

mutable struct PCConformer
    x::Vector{Float32}
    y::Vector{Float32}
    z::Union{Nothing,Vector{Float32}}

    style::Union{Nothing,PCDrawAnnotations}
    data::Union{Nothing,Vector{PCInfoData}} # Data Associated with this Conformer

    PCConformer() = new()
end
StructTypes.StructType(::Type{PCConformer}) = StructTypes.Mutable()

@enum(PCCoordinateType,
        CT_2D = 1,             # 2D Coordinates
        CT_3D,               # 3D Coordinates (should also indicate units, below)
        CT_SUBMITTED,        # Depositor Provided Coordinates
        CT_EXPERIMENTAL,     # Experimentally Determined Coordinates
        CT_COMPUTED,         # Computed Coordinates
        CT_STANDARDIZED,     # Standardized Coordinates
        CT_AUGMENTED,        # Hybrid Original with Computed Coordinates (e.g., explicit H)
        CT_ALIGNED,          # Template used to align drawing
        CT_COMPACT,          # Drawing uses shorthand forms (e.g., COOH, OCH3, Et, etc.)
        CT_UNITS_ANGSTROMS,  # (3D) Coordinate units are Angstroms
        CT_UNITS_NANOMETERS, # (3D) Coordinate units are nanometers
        CT_UNITS_PIXEL,      # (2D) Coordinate units are pixels
        CT_UNITS_POINTS,     # (2D) Coordinate units are points
        CT_UNITS_STDBONDS,   # (2D) Coordinate units are standard bond lengths (1.0)
        CT_UNITS_UNKNOWN = 255 # Coordinate units are unknown or unspecified
)
StructTypes.StructType(::Type{PCCoordinateType}) = StructTypes.NumberType()
StructTypes.numbertype(::Type{PCCoordinateType}) = UInt8

mutable struct PCCoordinates
    type::Vector{PCCoordinateType} # Coordinate Type Information (vector)
    aid::Vector{Int64}             # Conformer Atom IDs (vector) -- to be kept synchronized with Conformers

    conformers::Union{Nothing,Vector{PCConformer}}  # Conformers for this Coordinate Set
    atomlabels::Union{Nothing,Vector{PCAtomString}} # Atom labels for Conformer Set
    data::Union{Nothing,Vector{PCInfoData}}         # Data Associated with these Coordinates

    PCCoordinates() = new()
end
StructTypes.StructType(::Type{PCCoordinates}) = StructTypes.Mutable()

@enum(PCStereoGroupType,
  SGT_ABSOLUTE = 1,    # Absolute configuration is known
  SGT_OR,            # Relative configuration is known (absolute configuration is unknown)
  SGT_AND,           # Mixture of stereoisomers
  SGT_UNKNOWN = 255    # Unknown configuration type
)
StructTypes.StructType(::Type{PCStereoGroupType}) = StructTypes.NumberType()
StructTypes.numbertype(::Type{PCStereoGroupType}) = UInt8

mutable struct PCStereoGroup
    type::PCStereoGroupType
    aid::Vector{Int64}       # Atom Identifiers of atoms in this group

    PCStereoGroup() = new()
end
StructTypes.StructType(::Type{PCStereoGroup}) = StructTypes.Mutable()

mutable struct PCStereoPentagonalBiPyramid
    center::Int64   # Atom ID of Atom Center
    top::Int64      # Atom ID of Atom In-Plane and at the Top
    bottom::Int64   # Atom ID of Atom In-Plane and at the Bottom
    left::Int64     # Atom ID of Atom In-Plane and at the Left
    labove::Int64   # Atom ID of Atom Above the Plane on the Left
    lbelow::Int64   # Atom ID of Atom Below the Plane on the Left
    rabove::Int64   # Atom ID of Atom Above the Plane on the Right
    rbelow::Int64   # Atom ID of Atom Below the Plane on the Right

    PCStereoPentagonalBiPyramid() = new()
end
StructTypes.StructType(::Type{PCStereoPentagonalBiPyramid}) = StructTypes.Mutable()

mutable struct PCStereoTShape
    center::Int64   # Atom ID of Atom Center
    top::Int64      # Atom ID of Atom In-Plane and at the Top
    bottom::Int64   # Atom ID of Atom In-Plane and at the Bottom
    above::Int64    # Atom ID of Atom Above the Plane

    PCStereoTShape() = new()
end
StructTypes.StructType(::Type{PCStereoTShape}) = StructTypes.Mutable()

mutable struct PCStereoTrigonalBiPyramid
    center::Int64   # Atom ID of Atom Center
    above::Int64    # Atom ID of Atom Above the Plane
    below::Int64    # Atom ID of Atom Below the Plane
    top::Int64      # Atom ID of Atom In-Plane and at the Top
    bottom::Int64   # Atom ID of Atom In-Plane and at the Bottom
    right::Int64    # Atom ID of Atom In-Plane and to the Right

    PCStereoTrigonalBiPyramid() = new()
end
StructTypes.StructType(::Type{PCStereoTrigonalBiPyramid}) = StructTypes.Mutable()

mutable struct PCStereoOctahedral
    center::Int64   # Atom ID of Atom Center
    top::Int64      # Atom ID of Atom In-Plane and at the Top
    bottom::Int64   # Atom ID of Atom In-Plane and at the Bottom
    labove::Int64   # Atom ID of Atom Above the Plane on the Left
    lbelow::Int64   # Atom ID of Atom Below the Plane on the Left
    rabove::Int64   # Atom ID of Atom Above the Plane on the Right
    rbelow::Int64   # Atom ID of Atom Below the Plane on the Right

    PCStereoOctahedral() = new()
end
StructTypes.StructType(::Type{PCStereoOctahedral}) = StructTypes.Mutable()

@enum PCStereoCenterSquarePlanarType CSP_USHAPE = 1 CSP_ZSHAPE CSP_XSHAPE CSP_ANY CSP_UNKNOWN = 255
StructTypes.StructType(::Type{PCStereoCenterSquarePlanarType}) = StructTypes.NumberType()
StructTypes.numbertype(::Type{PCStereoCenterSquarePlanarType}) = UInt8

mutable struct PCStereoSquarePlanar
    center::Int64       # Atom ID of Atom Center
    lbelow::Int64       # Atom ID of Left Below Plane Atom
    rbelow::Int64       # Atom ID of Right Below Plane Atom
    labove::Int64       # Atom ID of Left Above Plane Atom
    rabove::Int64       # Atom ID of Right Above Plane Atom

    parity::Union{Nothing,PCStereoCenterSquarePlanarType} # StereoCenter Type

    PCStereoSquarePlanar() = new()
end
StructTypes.StructType(::Type{PCStereoSquarePlanar}) = StructTypes.Mutable()

@enum PCStereoCenterPlanarDesignation CP_SAME = 1 CP_OPPOSITE CP_ANY CP_UNKNOWN = 255
StructTypes.StructType(::Type{PCStereoCenterPlanarDesignation}) = StructTypes.NumberType()
StructTypes.numbertype(::Type{PCStereoCenterPlanarDesignation}) = UInt8

@enum PCStereoCenterPlanarType CP_PLANAR = 1 CP_CUMULENIC
StructTypes.StructType(::Type{PCStereoCenterPlanarType}) = StructTypes.NumberType()
StructTypes.numbertype(::Type{PCStereoCenterPlanarType}) = UInt8

mutable struct PCStereoPlanar
    left::Int64       # Atom ID of Left Double Bond Atom
    ltop::Int64       # Atom ID of Top Atom attached to the Left Double Bond Atom
    lbottom::Int64    # Atom ID of Bottom Atom attached to the Left Double Bond Atom
    right::Int64      # Atom ID of Top Atom attached to the Right Double Bond Atom
    rtop::Int64       # Atom Identifier of Atom Below the Plane
    rbottom::Int64    # Atom ID of Bottom Atom attached to the Right Double Bond Atom

    parity::Union{Nothing,PCStereoCenterPlanarDesignation} # StereoCenter Designation
    type::Union{Nothing,PCStereoCenterPlanarType}          # Type of StereoCenter, SP2 Planar, if not specified

    PCStereoPlanar() = new()
end
StructTypes.StructType(::Type{PCStereoPlanar}) = StructTypes.Mutable()

@enum PCStereoCenterTetrahedralDesignation TP_CLOCKWISE = 1 TP_COUNTERCLOCKWISE TP_ANY TP_UNKNOWN = 255
StructTypes.StructType(::Type{PCStereoCenterTetrahedralDesignation}) = StructTypes.NumberType()
StructTypes.numbertype(::Type{PCStereoCenterTetrahedralDesignation}) = UInt8

@enum PCStereoCenterTetrahedralType TP_TETRAHEDRAL = 1 TP_CUMULENIC TP_BIARYL
StructTypes.StructType(::Type{PCStereoCenterTetrahedralType}) = StructTypes.NumberType()
StructTypes.numbertype(::Type{PCStereoCenterTetrahedralType}) = UInt8

mutable struct PCStereoTetrahedral
    center::Int64 # Atom Identifier of Atom Center
    above::Int64  # Atom Identifier of Atom Above the Plane
    top::Int64    # Atom Identifier of Atom In-Plane and at the Top
    bottom::Int64 # Atom Identifier of Atom In-Plane and at the Bottom
    below::Int64  # Atom Identifier of Atom Below the Plane

    parity::Union{Nothing,PCStereoCenterTetrahedralDesignation} # StereoCenter Designation
    type::Union{Nothing,PCStereoCenterTetrahedralType}          # Type of StereoCenter, Tetrahedral, if not specified

    PCStereoTetrahedral() = new()
end
StructTypes.StructType(::Type{PCStereoTetrahedral}) = StructTypes.Mutable()

mutable struct PCStereoCenter
    tetrahedral::Union{Nothing,PCStereoTetrahedral}        # Tetrahedral (SP3) StereoCenter
    planar::Union{Nothing,PCStereoPlanar}                  # Planar (SP2) StereoCenter
    squareplanar::Union{Nothing,PCStereoSquarePlanar}      # Square Planar (SP4) StereoCenter
    octahedral::Union{Nothing,PCStereoOctahedral}          # Octahedral (OC-6) / Square Pyramid (SPY-5) StereoCenters
    bipyramid::Union{Nothing,PCStereoTrigonalBiPyramid}    # Trigonal BiPyramid (TBPY-4 and TBPY-5) StereoCenters
    tshape::Union{Nothing,PCStereoTShape}                  # T-Shaped (TS-3) StereoCenters
    pentagonal::Union{Nothing,PCStereoPentagonalBiPyramid} # Pentagonal BiPyramid (PBPY-7) StereoCenters

    PCStereoCenter() = new()
end
StructTypes.StructType(::Type{PCStereoCenter}) = StructTypes.Mutable()

@enum PCBondType PCBOND_SINGLE = 1 PCBOND_DOUBLE PCBOND_TRIPLE PCBOND_QUADRUPLE PCBOND_DATIVE PCBOND_COMPLEX PCBOND_IONIC PCBOND_BOND_UNKNOWN = 255
StructTypes.StructType(::Type{PCBondType}) = StructTypes.NumberType()
StructTypes.numbertype(::Type{PCBondType}) = UInt8

mutable struct PCBonds
    aid1::Vector{Int64}
    aid2::Vector{Int64}
    order::Vector{PCBondType}

    PCBonds() = new()
end
StructTypes.StructType(::Type{PCBonds}) = StructTypes.Mutable()

mutable struct PCMMDBSource
    mmdb_id::Int64
    molecule_id::Int64
    molecule_name::Vector{String}
    residue_id::Union{Nothing,Int64}
    residue_name::Union{Nothing,String}
    atom_id::Union{Nothing,Int64}
    atom_name::Union{Nothing,String}

    PCMMDBSource() = new(0, 0, [], nothing, nothing, nothing, nothing)
end
StructTypes.StructType(::Type{PCMMDBSource}) = StructTypes.Mutable()

# we need to map the json names, which contain -, to _-using symbols
fieldnames = []
StructTypes.foreachfield((i, name, FT, v) -> push!(fieldnames, (name, Symbol(replace(String(name), "_" => "-")))), PCMMDBSource())
StructTypes.names(::Type{PCMMDBSource}) = Tuple(fieldnames)

mutable struct PCAtomSource
    aid::Int64
    source::PCMMDBSource

    PCAtomSource() = new()
end
StructTypes.StructType(::Type{PCAtomSource}) = StructTypes.Mutable()

@enum PCAtomRadicalType SINGLET = 1 DOUBLET TRIPLET QUARTET QUINTET HEXTET HEPTET OCTET NONE = 255
StructTypes.StructType(::Type{PCAtomRadicalType}) = StructTypes.NumberType()
StructTypes.numbertype(::Type{PCAtomRadicalType}) = UInt8

mutable struct PCAtomRadical
    aid::Int64
    value::PCAtomRadicalType

    PCAtomRadical() = new()
end
StructTypes.StructType(::Type{PCAtomRadical}) = StructTypes.Mutable()

# Type definitions used by JSON3, derived from https://www.ncbi.nlm.nih.gov/IEB/ToolBox/CPP_DOC/asn_spec/pcsubstance.asn.html
@enum(PCElement,
      H = 1, HE, LI, BE, B, C, N, O, F, NE, NA, MG, AL, SI, P, S, CL, AR, K, CA, SC,
      TI, V, CR, MN, FE, CO, NI, CU, ZN, GA, GE, AS, SE, BR, KR, RB, SR, Y, ZR, NB,
      MO, TC, RU, RH, PD, AG, CD, IN, SN, SB, TE, I, XE, CS, BA, LA, CE, PR, ND, PM,
      SM, EU, GD, TB, DY, HO, ER, TM, YB, LU, HF, TA, W, RE, OS, IR, PT, AU, HG, TL,
      PB, BI, PO, AT, RN, FR, RA, AC, TH, PA, U, NP, PU, AM, CM, BK, CF, ES, FM, MD,
      NO, LR, RF, DB, SG, BH, HS, MT, DS, RG, CN, NH, FL, MC, LV, TS, OG,
      LP = 252, R = 253, D = 254, A = 255)
StructTypes.StructType(::Type{PCElement}) = StructTypes.NumberType()
StructTypes.numbertype(::Type{PCElement}) = UInt8

mutable struct PCAtoms
    aid::Vector{Int64}
    element::Vector{PCElement}

    label::Union{Nothing,Vector{PCAtomString}}
    isotope::Union{Nothing,Vector{PCAtomInt}}
    charge::Union{Nothing,Vector{PCAtomInt}}
    radical::Union{Nothing,Vector{PCAtomRadical}}
    source::Union{Nothing,Vector{PCAtomSource}}
    comment::Union{Nothing,Vector{PCAtomString}}

    PCAtoms() = new([], [], nothing, nothing, nothing, nothing, nothing, nothing)
end
StructTypes.StructType(::Type{PCAtoms}) = StructTypes.Mutable()

mutable struct PCCid
    cid::Int64

    PCCid() = new()
end
StructTypes.StructType(::Type{PCCid}) = StructTypes.Mutable()

mutable struct PCSid
    sid::Int64

    PCSid() = new()
end
StructTypes.StructType(::Type{PCSid}) = StructTypes.Mutable()

mutable struct PCXid
    xid::Int64

    PCXid() = new()
end
StructTypes.StructType(::Type{PCXid}) = StructTypes.Mutable()

@enum PCCompoundTypeID CTI_DEPOSITED = 0 CTI_STANDARDIZED CTI_COMPONENT CTI_NEUTRALIZED CTI_MIXTURE CTI_TAUTOMER CTI_PKASTATE CTI_UNKNOWN = 255
StructTypes.StructType(::Type{PCCompoundTypeID}) = StructTypes.NumberType()
StructTypes.numbertype(::Type{PCCompoundTypeID}) = UInt8

mutable struct PCCompoundType
    type::Union{Nothing,PCCompoundTypeID}
    id::Union{Nothing,PCCid,PCSid,PCXid}

    PCCompoundType() = new(CTI_UNKNOWN, nothing)
end
StructTypes.StructType(::Type{PCCompoundType}) = StructTypes.Mutable()

mutable struct PCCompound
    id::PCCompoundType
    atoms::Union{Nothing,PCAtoms}
    bonds::Union{Nothing,PCBonds}
    stereo::Union{Nothing,Vector{PCStereoCenter}}      # StereoCenter Descriptions
    coords::Union{Nothing,Vector{PCCoordinates}}       # 2D/3D Coordinate Sets of Compound
    charge::Int32                                      # Provided Total Formal Charge  (Signed Integer)
    props::Union{Nothing,Vector{PCInfoData}}           # Derived (computed) Properties
    stereogroups::Union{Nothing,Vector{PCStereoGroup}} # Relative stereochemistry groups
    count::Union{Nothing,PCCount}                      # Counts of various properties
    vbalt::Union{Nothing,Vector{PCCompound}}           # Alternate Valence-Bond Forms

    PCCompound() = new()
end
StructTypes.StructType(::Type{PCCompound}) = StructTypes.Mutable()

mutable struct PCResult
    PC_Compounds::Vector{PCCompound}

    PCResult() = new()
end
StructTypes.StructType(::Type{PCResult}) = StructTypes.Mutable()

ext: jl


text:
ACTIVATIONS = [
    :line, :σ, :sigmoid, :hardσ, :hardsigmoid, :hardtanh, :relu,
    :leakyrelu, :relu6, :rrelu, :prelu, :elu, :gelu, :swish, :selu,
    :celu, :softplus, :softsign, :logσ, :logcosh,
    :mish, :tanhshrink, :softshrink, :trelu, :lisht,
    :gaussian, :GCU, :SQU, :NCU, :SSU, :DSU
]
for f in ACTIVATIONS
    @eval export $(f)
end

"""
    line(x)
Linear function. This is the expression:
```math
line(x) = x
```
"""
line(x) = x

@doc raw"""
    σ(x)
Standard sigmoid activation function. Also, this function can be called with `σ`. This is the expression:
```math
\sigma(x) = \frac{1}{1+e^{-x}}
```
"""
function sigmoid(x)
    t = exp(-abs(x))
    return (x>=0) ? inv(1 + t) : t / (t + 1)
end
const σ = sigmoid

@doc raw"""
    hardsigmoid(x) = max(0, min(1, (x + 2.5) / 6))
Piecewise linear approximation of sigmoid. Also, this function can be called with `hardσ`. This is the expression:
```math
hardsigmoid(x) = \left\{
\begin{array}{ll}
1 & (x \geq \frac{1}{4}) \\
\frac{1}{5} x & (- \frac{1}{4} \lt x \lt \frac{1}{4}) \\
0 & (x \leq - \frac{1}{4})
\end{array}
\right.
```
"""
hardsigmoid(x) = max(0, min(1, (x + 2.5) / 6))
const hardσ = hardsigmoid

@doc raw"""
    lisht(x)
This is the expression:
```math
lisht(x) = x\tanh(x)
```
"""
lisht(x) = x*tanh(x)

@doc raw"""
    relu(x) = max(0, x)
`relu` is `Rectified Linear Unit`. This is the expression:
```math
relu(x) = \left\{
\begin{array}{ll}
x & (x \geq 0) \\
0 & (x \lt 0)
\end{array}
\right.
```
"""
relu(x) = max(0, x)

@doc raw"""
    leakyrelu(x; α=0.01) = (x>0) ? x : α*x
Leaky Rectified Linear Unit. This is the expression:
```math
leakyrelu(x) = \left\{
\begin{array}{ll}
\alpha x & (x \lt 0) \\
x & (x \geq 0)
\end{array}
\right.
```
"""
leakyrelu(x; α=0.01f0) = max(x, α*x)

@doc raw"""
    rrelu(min, max)
Randomized Rectified Linear Unit. The expression is the as [`leakyrelu`](@ref), but `α` is a random number between `min` and `max`.
Also, since this function is defined as a structure, use it as follows:
```
Dense(10=>5, rrelu(0.001, 0.1))
```
"""
struct rrelu
    α::Union{Float64, Float32}
    function rrelu(min::Float64, max::Float64)
        new(round(rand()*(max-min+1))+min)
    end
end

(rrelu::rrelu)(x) = max(x, rrelu.α*x)

@doc raw"""
    prelu(; α=0.01)
Parametric Ractified LinearUnit. The expression is the as [`leakyrelu`](@ref), but `α` is determined by learning. Also, when using this function, use [`Denseσ`](@ref) instead of [`Dense`](@ref).
"""
struct prelu
    α::Union{Float64, Float32}
    prelu(; α=0.01) = new(α)
end

(prelu::prelu)(x) = max(x, prelu.α*x)

@doc raw"""
    relu6(x)
Relu function with an upper limit of 6. This is the expression:
```math
relu6(x) = \left\{
\begin{array}{ll}
6 & (x \gt 6) \\
x & (x \geq 0) \\
0 & (x \lt 0)
\end{array}
\right.
```
"""
relu6(x) = min(6, max(0, x))

@doc raw"""
    elu(x, α=1)
Exponential Linear Unit activation function. You can also specify the coefficient explicitly, e.g. elu(x, 1). This is the expression:
```math
elu(x, α) = \left\{
\begin{array}{ll}
x & (x \geq 0) \\
\alpha(e^x-1) & (x \lt 0)
\end{array}
\right.
```
"""
elu(x; α=1) = (x<0) ? (exp(x)-1)α : x

@doc raw"""
    celu(x; α=1)
Continuously Differentiable Exponential Linear Unit. This is the expression:
```math
\alpha = 1 \\
celu(x) = \left\{
\begin{array}{ll}
x & (x \geq 0) \\
\alpha(e^\frac{x}{\alpha}-1) & (x \lt 0)
\end{array}
\right.
```
"""
celu(x; α=1) = (x>=0) ? x : α*(exp(x/α)-1)

@doc raw"""
    gelu(x)
Gaussian Error Linear Unit. This is the expression(``\phi`` is a distribution function of standard normal distribution.):
```math
gelu(x) = x\phi(x)
```
However, in the implementation, it is calculated with the following expression.
```math
\sigma(x) = \frac{1}{1+e^{-x}} \\
gelu(x) = x\sigma(1.702x)
```
"""
gelu(x) = σ(1.702*x)x

@doc raw"""
    selu(x)
Scaled exponential linear units. This is the expression
```math
\lambda = 1.0507009873554804934193349852946 \\
\alpha = 1.6732632423543772848170429916717 \\
selu(x) = \lambda \left\{
\begin{array}{ll}
x & (x \geq 0) \\
\alpha(e^x-1) & (x \lt 0)
\end{array}
\right.
```
"""
function selu(x)
    λ = oftype(float(x), selu_λ)
    α = oftype(float(x), selu_α)
    return ((x > 0) ? x : (exp(x)-1)α)*λ
end

const selu_λ = 1.0507009873554804934193349852946
const selu_α = 1.6732632423543772848170429916717

@doc raw"""
    trelu(x; θ=1)
Threshold gated Rectified Linear Unit. This is the expression:
```math
\theta = 1 \\
trelu(x) = \left\{
\begin{array}{ll}
x & (x \gt 0) \\
0 & (x \leq 0)
\end{array}
\right.
```
"""
trelu(x; θ=1) = (x<=θ) ? 0 : x

@doc raw"""
    logσ(x)
logarithmic sigmoid function. This is the expression:
```math
\sigma(x) = \frac{1}{1+e^{-x}} \\
logsigmoid(x) = \log(\sigma(x))
```
"""
logsigmoid(x) = log(σ(x))
const logσ = logsigmoid

@doc raw"""
    logcosh(x)
Log-Cosh function. This is the expression:
```math
logcosh(x) = \log(\cosh(x))
```
"""
logcosh(x) = log(cosh(x))

@doc raw"""
    hardtanh(x)
Linear tanh function. This is the expression:
```math
hardtanh(x) = \left\{
\begin{array}{ll}
1 & (x \geq 1) \\
x & (-1 \lt x \lt 1) \\
-1 & (x \leq -1)
\end{array}
\right.
```
"""
hardtanh(x) = min(1, max(-1, x))

@doc raw"""
    tanhshrink(x)
Shrink tanh function. This is the expression:
```math
tanhshrink(x) = 1-\tanh(x)
```
"""
tanhshrink(x) = 1-tanh(x)

@doc raw"""
    softshrink(x; λ=0.5)
This is the expression:
```math
\lambda=0.5 \\
softshrink(x) = \left\{
\begin{array}{ll}
x-\lambda & (x \gt \lambda) \\
0 & (-\lambda \leq x \leq \lambda) \\
x+\lambda & (x \lt -\lambda) \\
\end{array}
\right.
```
"""
softshrink(x; λ=0.5) = min(x+λ, max(x-λ, 0))

@doc raw"""
    softsign(x) = x / (1+abs(x))
The softsign activation function. This is the expression:
```math
softsign(x) = \frac{x}{1+|x|}
```
"""
softsign(x) = x / (1 + abs(x))

@doc raw"""
    softplus(x) = log(1 + exp(x))
the softplus activation function. This is the expression:
```math
softplus(x) = \ln(1+e^x)
```
"""
softplus(x) = log(1 + exp(x))

@doc raw"""
    mish(x) = x * tanh(softplus(x))
The mish function. This is the expression:
```math
softplus(x) = \ln(1+e^x) \\
mish(x) = x\tanh(softplus(x))
```
"""
mish(x) = x * tanh(softplus(x))

@doc raw"""
    swish(x; β=1)
The swish function. This is the expression:
```math
\sigma(x) = \frac{1}{1+e^{-x}} \\
swish(x) = x\sigma(\beta x)
```
"""
swish(x; β=1) = x * σ(β*x)

@doc raw"""
    gaussian(x)
The Gauss Function. This is the expression:
```math
Gaussian(x) = e^{-x^{2}}
```
"""
gaussian(x) = exp(-x^2)

@doc raw"""
    GCU(x)
Growing Cosine Unit. This is the expression:
```math
GCU(x) = x\cos(x)
```
"""
GCU(x) = x*cos(x)

@doc raw"""
    SQU(x)
Shifted Quadratic Unit. SQU is a biologically inspired activation that enables single neurons to learn the XOR function. This is the expression:
```math
SQU(x) = x^{2}+x
```
"""
SQU(x) = x^2+x

@doc raw"""
    NCU(x)
Non-Monotonic Cubic Unit. This is the expression:
```math
NCU(x) = x-x^{3}
```
"""
NCU(x) = x-x^3

@doc raw"""
    SSU(x)
Shifted Sinc Unit. This is the expression:
```math
SSU(x) = \pi sinc(x - \pi)
```
"""
SSU(x) = pi*sinc(x-pi)

@doc raw"""
    DSU(x)
Decaying Sine Unit. This is the expression:
```math
DSU(x) = \frac{\pi}{2}(sinc(x-\pi)-sinc(x+\pi))
```
"""
DSU(x) = (pi/2)*(sinc(x-pi)-sinc(x+pi))


text:
# ==================
# Inhomogeneous case
# ==================

# no invariant
function reach_inhomog_GLGM06!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
                               Ω0::Zonotope{N, VN, MN},
                               Φ::AbstractMatrix,
                               NSTEPS::Integer,
                               δ::Float64,
                               max_order::Integer,
                               X::Universe,
                               U::LazySet,
                               reduction_method::AbstractReductionMethod,
                               Δt0::TimeInterval,
                               disjointness_method::AbstractDisjointnessMethod) where {N, VN, MN}

    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    Wk₊ = U
    Φ_power_k = similar(Φ)
    copyto!(Φ_power_k, Φ)
    Φ_power_k_cache = similar(Φ)

    k = 2
    @inbounds while k <= NSTEPS
        Rₖ = minkowski_sum(linear_map(Φ_power_k, Ω0), Wk₊)
        Rₖ = _reduce_order(Rₖ, max_order, reduction_method)

        Δt += δ
        F[k] = ReachSet(Rₖ, Δt)

        Wk₊ = minkowski_sum(Wk₊, linear_map(Φ_power_k, U))
        Wk₊ = _reduce_order(Wk₊, max_order, reduction_method)

        mul!(Φ_power_k_cache, Φ_power_k, Φ)
        copyto!(Φ_power_k, Φ_power_k_cache)
        k += 1
    end
    return F
end

# with invariant
function reach_inhomog_GLGM06!(F::Vector{ReachSet{N, Zonotope{N, VN, MN}}},
                               Ω0::Zonotope{N, VN, MN},
                               Φ::AbstractMatrix,
                               NSTEPS::Integer,
                               δ::Float64,
                               max_order::Integer,
                               X::LazySet,
                               U::LazySet,
                               reduction_method::AbstractReductionMethod,
                               Δt0::TimeInterval,
                               disjointness_method::AbstractDisjointnessMethod) where {N, VN, MN}

    # initial reach set
    Δt = (zero(N) .. δ) + Δt0
    @inbounds F[1] = ReachSet(Ω0, Δt)

    Wk₊ = U
    Φ_power_k = similar(Φ)
    copyto!(Φ_power_k, Φ)
    Φ_power_k_cache = similar(Φ)

    k = 2
    @inbounds while k <= NSTEPS
        Rₖ = minkowski_sum(linear_map(Φ_power_k, Ω0), Wk₊)
        Rₖ = _reduce_order(Rₖ, max_order, reduction_method)
        _is_intersection_empty(X, Rₖ, disjointness_method) && break
        Δt += δ
        F[k] = ReachSet(Rₖ, Δt)

        Wk₊ = minkowski_sum(Wk₊, linear_map(Φ_power_k, U))
        Wk₊ = _reduce_order(Wk₊, max_order, reduction_method)

        mul!(Φ_power_k_cache, Φ_power_k, Φ)
        copyto!(Φ_power_k, Φ_power_k_cache)
        k += 1
    end
    if k < NSTEPS + 1
        resize!(F, k-1)
    end
    return F
end


text:
### Computer Vision

# Export methods

export  AFFeatures,
        orb,
        sift,
        gloh,
        diffOfGaussians,
        fast,
        harris,
        susan,
        hammingMatcher,
        nearestNeighbour,
        matchTemplate

# Export constants

export  AF_SAD,
        AF_ZSAD,
        AF_LSAD,
        AF_ZSSD,
        AF_LSSD,
        AF_NCC,
        AF_ZNCC,
        AF_SHD

# Feature Type

immutable AFFeatures
    ptr::Ptr{Void}
end

# Constants

AF_SAD = 0
AF_ZSAD = 1
AF_LSAD = 2
AF_ZSSD = 3
AF_LSSD = 4
AF_NCC = 5
AF_ZNCC = 6
AF_SHD = 7

# Feature Descriptors

function orb(a::AFArray; fast_thr = 20., max_feat = 400, scl_fctr = 1.5, levels = 4, blur_img = false)
    feat = new_ptr()
    desc = new_ptr()
    af_orb(feat, desc, a, fast_thr, Cuint(max_feat), scl_fctr, Cuint(levels), blur_img)
    AFFeatures(feat[]),
    AFArray{backend_eltype(desc[])}(desc[])
end

for (op,fn) in ((:sift, :af_sift), (:gloh, :af_gloh))

    @eval function ($op)(a::AFArray; n_layers = 3, constant_thr = 0.04, edge_thr = 0.04,
                    init_sigma = 1.6, double_input = true, intensity_scale = 0.00390625,
                    feature_ratio = 0.05)
        feat = new_ptr()
        desc = new_ptr()
        eval($fn)(feat, desc, a, Cuint(n_layers), contrast_thr, edge_thr, init_sigma,
                double_input, intensity_scale, feature_ratio)
        AFFeatures(feat[]),
        AFArray{backend_eltype(desc[])}(desc[])
    end

end

# Feature Detectors

function diffOfGaussians(a::AFArray, radius1::Int, radius2::Int)
    out = new_ptr()
    af_dog(out, a, radius1, radius)
    AFArray{backend_eltype(out[])}(out[])
end

function fast(a::AFArray; thr = 20., arc_length = 9, non_max = true,
                feature_ratio = 0.05, edge = 3)
    out = new_ptr()
    af_fast(out, a, thr, Cuint(arc_length), non_max, feature_ratio, Cuint(edge))
    AFFeatures(out[])
end

function harris(a::AFArray; max_corners = 500, min_response = 1e-5,
                sigma = 1., block_size = 0, k_thr = 0.04)
    out = new_ptr()
    af_harris(out, a, Cuint(max_corners), min_response, sigma, Cuint(block_size), k_thr)
    AFFeatures(out[])
end

function susan(a::AFArray; radius = 3, diff_thr = 32.0,
                geom_thr = 10.0, feature_ratio = 0.05, edge = 3)
    out = new_ptr()
    af_susan(out, a, Cuint(radius), diff_thr, geom_thr, feature_ratio, Cuint(edge))
    AFFeatures(out[])
end

# Feature Matchers

function hammingMatcher(a::AFArray, train::AFArray; dist_dim = 0, n_dist = 1)
    idx = out_ptr()
    dist = out_ptr()
    af_hamming_matcher(idx, dist, a, train, dist_dim, Cuint(n_dist))
    AFArray{backend_eltype(idx[])}(idx[]) + 1,
    AFArray{backend_eltype(dist[])}(dist[])
end

function nearestNeighbour(a::AFArray, train::AFArray; dist_dim = 9,
                            n_dist = 1, dist_type = AF_SSD)
    idx = new_ptr()
    dist = new_ptr()
    af_nearest_neighbour(idx, dist, query, train, dist_dim, n_dist, dist_type)
    AFArray{backend_eltype(idx[])}(idx[]) + 1,
    AFArray{backend_eltype(dist[])}(dist[])
end

function matchTemplate(search_img::AFArray, template_img::AFArray; m_typ = AF_SAD)
    out = new_ptr()
    af_match_template(out, search_img, template_img, m_type)
    AFArray{backend_eltype(out[])}(out[])
end


text:
using GeometricIntegrators
using ChargedParticleDynamics.GuidingCenter4d.TokamakFastSurface
using GeometricExamples

const Δt     = 1.
const ntime  = 25000
const nx     = 200
const ny     = 200
const nsave  = 50
const nplot  = 5
const run_id = "poincare_2nd_fast_vprk_dt1"

include("guiding_center_4d_settings_vprk.jl")

tableau_list = get_tableau_list_vprk_projection()

pinv = guiding_center_4d_iode_poincare_invariant_2nd(Δt, nx, ny, ntime, nsave)

include("guiding_center_4d_fast_poincare_invariant_2nd.jl")


text:
module StrangQuadraturesTests

using Test
using Gridap.Fields
using Gridap.ReferenceFEs

for k in keys(ReferenceFEs._strang_tri_k2n)
  reffe = ReferenceFE(TRI,lagrangian,Float64,k)
  quad1 = Quadrature(TRI,strang,k)
  quad2 = Quadrature(TRI,duffy,k)
  f = get_shapefuns(reffe)
  f1 = integrate(f,get_coordinates(quad1),get_weights(quad1))
  f2 = integrate(f,get_coordinates(quad2),get_weights(quad2))
  err = maximum(abs.(f1 .- f2))/maximum(abs.(f1))
  @test err < 1.0e-7
  #@show (2,k,err)
end

for k in keys(ReferenceFEs._strang_tet_k2n)
  reffe = ReferenceFE(TET,lagrangian,Float64,k)
  quad1 = Quadrature(TET,strang,k)
  quad2 = Quadrature(TET,duffy,k)
  f = get_shapefuns(reffe)
  f1 = integrate(f,get_coordinates(quad1),get_weights(quad1))
  f2 = integrate(f,get_coordinates(quad2),get_weights(quad2))
  err = maximum(abs.(f1 .- f2))/maximum(abs.(f1))
  @test err < 1.0e-13
  #@show (3,k,err)
end

end # module


text:
using Test
using SMACOF
using Distances
using LinearAlgebra
using Random
using Statistics
Random.seed!(1993)
ϵ = 2e-6

@testset "vanilla smacof" begin
    include("test_smacof.jl")
end

# @testset "wda smacof" begin
#     include("test_wda_smacof.jl")
# end

# @testset "classical scaling" begin
#     X = rand(50, 2)
#     Δ = dist(X)
#     X̂ = classical_scaling(Δ)
#     Δ̂ = dist(X̂)
#     @test norm(Δ - Δ̂) < ϵ
# end

# @testset "issymmetric" begin
#     X = randn(5, 2)
#     @test issymmetric(X) == false

#     Δ = dist(X)
#     @test issymmetric(Δ) == true

#     B = randn(5, 5)
#     B = B + tril(ones(5, 5))
#     @test issymmetric(B) == false
# end

# @testset "smacof float32" begin
#     X = rand(Float32, 50, 2)
#     Δ = dist(X)
#     # params = SmacofParams(Δ, anchors = X)
#     Y = smacof(Δ, anchors = X)
#     @test SMACOF.mse(Y, X) < ϵ
# end

# @testset "(wda) smacof with n=5" begin
#     X = rand(5, 2)
#     Δ = dist(X)
#     # params = SmacofParams(Δ, anchors = X)
#     Y_smacof = smacof(Δ, anchors = X, Xinit = rand(5, 2), verbose = false)
#     @test SMACOF.mse(Y_smacof, X) < ϵ

#     Y_wda = wda_smacof(Δ, anchors = X, DA_itmax = 300, verbose = false)
#     @test SMACOF.mse(Y_wda, X) < ϵ
# end

# @testset "higher dimensions" begin
#     # X = rand(100, 5)
#     # X[:, 3:end] .= 0.0
#     # X = SMACOF.rotate(X)

#     # Δ = dist(X)
#     # Y_smacof = smacof(Δ, anchors = X)
#     # @test SMACOF.mse(Y_smacof, X) < ϵ

#     function higher_dim_test(n, p)
#         X = zeros(n, p)
#         X[:, 1:2] = rand(n, 2)
#         X = SMACOF.rotate(X)

#         Δ = dist(X)
#         Y_smacof = smacof(Δ, anchors = X)
#         @test SMACOF.mse(Y_smacof, X) < ϵ
#     end

#     higher_dim_test(100, 5)
#     higher_dim_test(500, 10)
# end

# @testset "recover noiseless configuration" begin
#     nlst = [10, 20, 30]
#     X = Dict(n => rand(n, 2) for n in nlst)
#     Δ = Dict(n => dist(X[n]) for n in nlst)
#     ϵ = 1e-6

#     @testset "smacof" begin
#     for n in nlst
#             Xhat = smacof(Δ[n], anchors = X[n])
#             @test SMACOF.mse(X[n], Xhat) < ϵ
#         end
#     end

#     @testset "wda_smacof" begin
#         for n in nlst
#             Xhat = wda_smacof(Δ[n], anchors = X[n])
#             @test SMACOF.mse(X[n], Xhat) < ϵ
#         end
#     end

#     # With varying sizes
#     nlst = [10, 20, 30]
#     X = Dict(n => (n / 2) * rand(n, 2) for n in nlst)
#     Δ = Dict(n => dist(X[n]) for n in nlst)
#     ϵ = 1e-6

#     @testset "smacof" begin
#     for n in nlst
#             Xhat = smacof(Δ[n], anchors = X[n])
#             @test SMACOF.mse(X[n], Xhat) < ϵ
#         end
#     end

#     @testset "wda_smacof" begin
#         for n in nlst
#             Xhat = wda_smacof(Δ[n], anchors = X[n])
#             @test SMACOF.mse(X[n], Xhat) < ϵ
#         end
#     end
# end

# # # @test 1 == 0

# # @testset "Procrustes" begin
# # A = [0.0 0.0;
# #          0   1;
# #          1   0]
#             #     # TRANSLATION
# #     B = A .+ 0.1
#         #     Ahat =  SMACOF.align(B, A)
# #     for i in eachindex(Ahat)
# #     @test Ahat[i] ≈ A[i] atol = 1e-8
# #     end

# #     # ROTATION
# #     B = SMACOF.rotate(A)
# #         Ahat =  SMACOF.align(B, A)
# #     for i in eachindex(Ahat)
# #     @test Ahat[i] ≈ A[i] atol = 1e-8
# #     end

# #     p = SMACOF.Procrustes(B, A)
# #     Ahat = SMACOF.transform(B, p)
# #     @test norm(A - Ahat) ≈ 0.0 atol = 1e-8
# # end

# # @testset "noiseless configuration with some anchors" begin
# #     nlst = [10, 20, 30]
# #     n_anchors = 5
# #     X = Dict(n => rand(n, 2) for n in nlst)
# #     anchors = Dict(n => (idx = 1:n_anchors, pos = X[n][1:n_anchors, :]) for n in nlst)
# #     Δ = Dict(n => dist(X[n]) for n in nlst)
# #     ϵ = 1e-6

# #     @testset "smacof" begin
# #     for n in nlst
# #         Xhat = smacof(SmacofParams(Δ[n], anchors = anchors[n]))
# #         # Xhat = fit(Smacof(Δ[n]), anchors = anchors[n])
# #         @test SMACOF.mse(X[n], Xhat) < ϵ
# #         end
# #     end

# #     @testset "da_smacof" begin
# #             for n in nlst
# #         Xhat = da_smacof(Δ[n], anchors = anchors[n])
# #             @test SMACOF.mse(X[n], Xhat) < ϵ
# #         end
# #     end

# #     @testset "wda_smacof" begin
# #         for n in nlst
# #             Xhat = wda_smacof(Δ[n], anchors = anchors[n])
# #             @test SMACOF.mse(X[n], Xhat) < ϵ
# #         end
# #     end
# #             end

# #             # TODO: @test make_gif()
# # #       makegif returns true if successful


text:
aws.region = "us-east-1"
bucket_name = "ocaws.jl.test." * lowercase(Dates.format(now(Dates.UTC), "yyyymmddTHHMMSSZ"))

@testset "Create Bucket" begin
    s3_create_bucket(aws, bucket_name)
    @test bucket_name in s3_list_buckets(aws)
    s3_enable_versioning(aws, bucket_name)
    sleep(1)
end

@testset "Bucket Tagging" begin
    @test isempty(s3_get_tags(aws, bucket_name))
    tags = Dict("A" => "1", "B" => "2", "C" => "3")
    s3_put_tags(aws, bucket_name, tags)
    @test s3_get_tags(aws, bucket_name) == tags
    s3_delete_tags(aws, bucket_name)
    @test isempty(s3_get_tags(aws, bucket_name))
end

@testset "Create Objects" begin
    s3_put(aws, bucket_name, "key1", "data1.v1")
    s3_put(bucket_name, "key2", "data2.v1", tags = Dict("Key" => "Value"))
    s3_put(aws, bucket_name, "key3", "data3.v1")
    s3_put(aws, bucket_name, "key3", "data3.v2")
    s3_put(aws, bucket_name, "key3", "data3.v3"; metadata = Dict("foo" => "bar"))
    s3_put(aws, bucket_name, "key4", "data3.v4"; acl="bucket-owner-full-control")
    s3_put_tags(aws, bucket_name, "key3", Dict("Left" => "Right"))

    @test isempty(s3_get_tags(aws, bucket_name, "key1"))
    @test s3_get_tags(aws, bucket_name, "key2")["Key"] == "Value"
    @test s3_get_tags(aws, bucket_name, "key3")["Left"] == "Right"
    s3_delete_tags(aws, bucket_name, "key2")
    @test isempty(s3_get_tags(aws, bucket_name, "key2"))

    @test s3_get(aws, bucket_name, "key1") == b"data1.v1"
    @test s3_get(aws, bucket_name, "key2") == b"data2.v1"
    @test s3_get(bucket_name, "key3") == b"data3.v3"
    @test s3_get(bucket_name, "key4") == b"data3.v4"
    @test s3_get_meta(bucket_name, "key3")["x-amz-meta-foo"] == "bar"
end

@testset "ASync Get" begin
    @sync begin
        for i in 1:2
            @async begin
                @test s3_get(bucket_name, "key3") == b"data3.v3"
            end
        end
    end
end

@testset "Raw Return - XML" begin
    xml = "<?xml version='1.0'?><Doc><Text>Hello</Text></Doc>"
    s3_put(aws, bucket_name, "file.xml", xml, "text/xml")
    @test String(s3_get(aws, bucket_name, "file.xml", raw=true)) == xml
    @test s3_get(aws, bucket_name, "file.xml")["Text"] == "Hello"
end

@testset "Get byte range" begin
    teststr = "123456789"
    s3_put(aws, bucket_name, "byte_range", teststr)
    range = 3:6
    @test String(s3_get(aws, bucket_name, "byte_range"; byte_range=range)) == teststr[range]
end

@testset "Object Copy" begin
    s3_copy(bucket_name, "key1"; to_bucket=bucket_name, to_path="key1.copy")
    @test s3_get(aws, bucket_name, "key1.copy") == b"data1.v1"
end

@testset "Sign URL" begin
    for v in ["v2", "v4"]
        url = s3_sign_url(aws, bucket_name, "key1"; signature_version=v)
        curl_output = ""

        @repeat 3 try
            curl_output = read(`curl -s -o - $url`, String)
        catch e
            @delay_retry if true end
        end

        @test curl_output == "data1.v1"

        fn = "/tmp/jl_qws_test_key1"
        if isfile(fn)
            rm(fn)
        end

        @repeat 3 try
            s3_get_file(aws, bucket_name, "key1", fn)
        catch e
            sleep(1)
            @retry if true end
        end

        @test read(fn, String) == "data1.v1"
        rm(fn)
    end
end

@testset "Object exists" begin
    for key in ["key1", "key2", "key3", "key1.copy"]
        @test s3_exists(bucket_name, key)
    end
end

@testset "List Objects" begin
    for key in ["key1", "key2", "key3", "key1.copy"]
        @test key in [o["Key"] for o in s3_list_objects(aws, bucket_name)]
    end
end

@testset "Object Delete" begin
    s3_delete(aws, bucket_name, "key1.copy")
    @test !("key1.copy" in [o["Key"] for o in s3_list_objects(aws, bucket_name)])
end

@testset "Check Metadata" begin
    meta = s3_get_meta(aws, bucket_name, "key1")
    @test meta["ETag"] == "\"68bc8898af64159b72f349b391a7ae35\""
end

@testset "Check Object Versions" begin
    versions = s3_list_versions(aws, bucket_name, "key3")
    @test length(versions) == 3
    @test (s3_get(aws, bucket_name, "key3"; version=versions[3]["VersionId"]) == b"data3.v1")
    @test (s3_get(aws, bucket_name, "key3"; version=versions[2]["VersionId"]) == b"data3.v2")
    @test (s3_get(aws, bucket_name, "key3"; version=versions[1]["VersionId"]) == b"data3.v3")
end

@testset "Purge Versions" begin
    s3_purge_versions(aws, bucket_name, "key3")
    versions = s3_list_versions(aws, bucket_name, "key3")
    @test length(versions) == 1
    @test s3_get(aws, bucket_name, "key3") == b"data3.v3"
end

@testset "default Content-Type" begin
    # https://github.com/samoconnor/AWSS3.jl/issues/24
    ctype(key) = s3_get_meta(bucket_name, key)["Content-Type"]

    for k in [
        "file.foo",
        "file",
        "file_html",
        "file/html",
        "foobar.html/file.htm"]

        s3_put(aws, bucket_name, k, "x")
        @test ctype(k) == "application/octet-stream"
    end

    for (k, t) in [
        ("foo/bar/file.html",  "text/html"),
        ("x.y.z.js",           "application/javascript"),
        ("downalods/foo.pdf",  "application/pdf"),
        ("data/foo.csv",       "text/csv"),
        ("this.is.a.file.txt", "text/plain"),
        ("my.log",             "text/plain"),
        ("big.dat",            "application/octet-stream"),
        ("some.tar.gz",        "application/octet-stream"),
        ("data.bz2",           "application/octet-stream")]

        s3_put(aws, bucket_name, k, "x")
        @test ctype(k) == t
    end
end

@testset "Multi-Part Upload" begin
    MIN_S3_CHUNK_SIZE = 5 * 1024 * 1024 # 5 MB
    key_name = "multi-part-key"
    upload = s3_begin_multipart_upload(aws, bucket_name, key_name)
    tags = Vector{String}()

    for part_number in 1:5
        push!(tags, s3_upload_part(aws, upload, part_number, rand(UInt8, MIN_S3_CHUNK_SIZE)))
    end

    s3_complete_multipart_upload(aws, upload, tags)
    @test s3_exists(bucket_name, key_name)
end

@testset "Empty and Delete Bucket" begin
    AWSS3.s3_nuke_bucket(aws, bucket_name)
    @test !in(bucket_name, s3_list_buckets(aws))
end

@testset "Delete Non-Existant Bucket" begin
    @test_throws AWS.AWSException s3_delete_bucket(aws, bucket_name)
end


text:
using HDF5, JLD
using Distributions
using Random

export runStochasticGradientHamiltonianMonteCarloSampler,
       doDPPBayesianLearningSparseVectorData

# Computes the gradient of the log-prior density
function computePriorGradient(paramsMatrix, numItems, numItemTraits,
                              gaussianPriorPrecisions)
  paramsMatrixGradient = zeros(numItems, numItemTraits)

  for paramsMatrixColIndex = 1:numItemTraits
    for paramsMatrixRowIndex = 1:numItems
      paramsMatrixGradient[paramsMatrixRowIndex, paramsMatrixColIndex] =
        gaussianPriorPrecisions[paramsMatrixRowIndex] * paramsMatrix[paramsMatrixRowIndex, paramsMatrixColIndex]
    end
  end

  return paramsMatrixGradient
end

# Computes the gradient of the unnormalized log-posterior density
function computePosteriorGradient(paramsMatrix, trainingInstances, numTrainingInstances,
                                  numItems, numItemTraits, gaussianPriorPrecisions)
  posteriorGradient = zeros(numItems, numItemTraits)

  # Compute the gradient for the prior density, and the gradient for the likelihood function
  priorGradient = computePriorGradient(paramsMatrix, numItems, numItemTraits,
                                       gaussianPriorPrecisions)
  likelihoodGradient = computeGradient(paramsMatrix, trainingInstances,
                                       numTrainingInstances, numItems, numItemTraits)

  posteriorGradient = priorGradient + likelihoodGradient

  return posteriorGradient
end

# Runs Stochastic Gradient Hamiltonian Monte Carlo (SGHMC) to generate samples from the posterior
function runStochasticGradientHamiltonianMonteCarloSampler(trainingInstances,
  numTrainingInstances, numItems, numItemTraits, testInstances, numTestInstances,
  learnedModelOutputDirName, initialItemTraitMatrix = -1; minibatchSize = 1000)

  # Initialize itemTraitMatrix (V matrix)
  itemTraitMatrix = rand(numItems, numItemTraits)
  itemTraitMatrixPrev = zeros(numItems, numItemTraits)
  if initialItemTraitMatrix != -1
    itemTraitMatrix = initialItemTraitMatrix
  end

  posteriorGradient = zeros(numItems, numItemTraits)

  itemTraitMatrixSamplesSum = zeros(numItems, numItemTraits)

  # Max number of iterations (samples)
  maxNumIterations = 1000

  # Number of "burn-in" samples, which are discarded
  numBurnInSamples = 100

  # Discard all but every sampleLag samples.  Used for thinning of samples.
  sampleLag = 10

  # Step size (learning rate)
  # stepSizeLarger = 1.0e-5
  stepSizeLarger = 1.0e-4
  stepSizeIntermediate = 1.0e-6
  stepSizeSmaller = 3.0e-8
  stepSize = stepSizeLarger

  alphaMomentum = 0.01
  alphaMomentumIntermediate = 0.01
  alphaMomentumSmaller = 0.01

  numIterationsLargerStepSize = 17
  numIterationsIntermediateStepSize = 95

  # Gradient noise estimate (beta hat)
  gradientNoiseEstimate = 0.0

  # Gradient noise distribution
  gradientNoiseGaussian = Normal(0, 2 * (alphaMomentum - gradientNoiseEstimate) * stepSize)

  # Number of "leapfrog" steps
  numLeapfrogSteps = 10

  # Number of training instances to process per minibatch
  # minibatchSize = 1000

  delta = zeros(numItems, numItemTraits)

  # Hyperparameters on the gamma hyperprior for each item
  gammaShapeHyperParam = sqrt(numItemTraits)
  gammaRateHyperParam = 1.0

  # Shape parameter for the gamma hyperprior for each item
  gammaShape = gammaShapeHyperParam + (numItems * numItemTraits) / 2.0

  # Number of samples collected from the posterior
  numSamplesCollected = 1

  # Create directory for storing collected samples, if it doesn't already exist
  collectedSamplesDir =
    "$learnedModelOutputDirName/learnedDPPParams-k$numItemTraits-SGHMC-collectedSamples"
  if !isdir(collectedSamplesDir)
    mkpath(collectedSamplesDir)
  end

  trainingLogLike = computeLogLikelihood(itemTraitMatrix,
    trainingInstances, numTrainingInstances, numItems, 0)

  currTrainingInstanceIndex = 1
  shuffle!(trainingInstances)
  for iterCounter = 1:maxNumIterations
    # Get current minibatch
    numTrainingInstancesInMinibatch = minibatchSize
    if currTrainingInstanceIndex + minibatchSize > numTrainingInstances
      numTrainingInstancesInMinibatch = numTrainingInstances - currTrainingInstanceIndex
    end
    minibatchTrainingInstances =
      trainingInstances[currTrainingInstanceIndex:(currTrainingInstanceIndex + numTrainingInstancesInMinibatch)]

    # Rate parameter for the gamma hyperprior for each item
    sumMagnitudeItemVectors = 0.0
    for i = 1:numItems
      sumMagnitudeItemVectors += vecnorm(itemTraitMatrix[i, :]) ^ 2
    end
    gammaRate = gammaRateHyperParam + sumMagnitudeItemVectors / 2

    # Sample new precision for the Gaussian prior on each item from a
    # gamma hyperprior on each precision.  All items share the same sampled
    # precision parameter.
    gammaHyperprior = Gamma(gammaShape, 1 / gammaRate)
    gaussianPriorPrecisions = zeros(numItems)
    gaussianPriorPrecisionShared = rand(gammaHyperprior)
    for itemIndex = 1:numItems
      gaussianPriorPrecisions[itemIndex] = gaussianPriorPrecisionShared
    end

    # Compute new itemTraitMatrix and delta, using the SGHMC updates
    itemTraitMatrixPrev = itemTraitMatrix
    @time for j = 1:numLeapfrogSteps
      itemTraitMatrix = itemTraitMatrix + delta

      posteriorGradient =
        computePosteriorGradient(itemTraitMatrix, minibatchTrainingInstances, numTrainingInstancesInMinibatch,
                                 numItems, numItemTraits, gaussianPriorPrecisions)

      gradientNoiseSample = rand(gradientNoiseGaussian)

      delta = stepSize * posteriorGradient - alphaMomentum * delta + gradientNoiseSample
    end

    if iterCounter % 1 == 0
      println("$iterCounter iterations completed; stepSize: $stepSize; alphaMomentum: $alphaMomentum")

      avgTrainingLogLikelihood =
        computeLogLikelihood(itemTraitMatrix, trainingInstances, numTrainingInstances,
                             numItems, numItemTraits) / numTrainingInstances
      avgTestLogLikelihood =
        computeLogLikelihood(itemTraitMatrix, testInstances, numTestInstances,
                             numItems, numItemTraits) / numTestInstances
      println("\t avgTrainingLogLikelihood: $avgTrainingLogLikelihood")
      println("\t avgTestLogLikelihood: $avgTestLogLikelihood")
    end

    if iterCounter >= numIterationsLargerStepSize && iterCounter < numIterationsIntermediateStepSize
      stepSize = stepSizeIntermediate
      alphaMomentum = alphaMomentumIntermediate

      if alphaMomentum > 0
        gradientNoiseGaussian = Normal(0, 2 * (alphaMomentum - gradientNoiseEstimate) * stepSize)
      end
    end

    if iterCounter >= numIterationsIntermediateStepSize
      stepSize = stepSizeSmaller
      alphaMomentum = alphaMomentumSmaller

      if alphaMomentum > 0
        gradientNoiseGaussian = Normal(0, 2 * (alphaMomentum - gradientNoiseEstimate) * stepSize)
      end
    end

    if iterCounter > 2
      prevTrainingLogLike = trainingLogLike
      trainingLogLike = computeLogLikelihood(itemTraitMatrix,
        trainingInstances, numTrainingInstances, numItems, 0)

      # Tracking training log likelihood is not a "true" measure of convergence
      # here, but provides a useful heuristic
      isConvergedLogLikelihood(prevTrainingLogLike, trainingLogLike,
        1.0e-7, "training", verbose = true)
    end

    # After burn in, collect samples
    if iterCounter > numBurnInSamples
      if iterCounter % sampleLag == 0
        @async save("$collectedSamplesDir/learnedDPPParams-k$numItemTraits-SGHMC-sample-$numSamplesCollected.jld",
          "itemTraitMatrixSample", itemTraitMatrix)

        println("Collected $numSamplesCollected samples")
        numSamplesCollected += 1
      end
    end

    currTrainingInstanceIndex += numTrainingInstancesInMinibatch + 1
    if currTrainingInstanceIndex > numTrainingInstances
      # We've processed all of the training instances, so start processing the
      # training instances again from the beginning
      currTrainingInstanceIndex = 1

      shuffle!(trainingInstances)
    end
  end
end

# Performs learning for the Bayesian low-rank DPP model, using an MCMC sampler
# function, on a dataset in sparse vector format.
function doDPPBayesianLearningSparseVectorData(trainingBasketsDictFileName,
  trainingBasketsDictObjectName, testBasketsDictFileName,
  testBasketsDictObjectName, learnedModelOutputDirName, numItemTraits,
  runSamplerFunction::Function, initialItemTraitMatrxFileName = -1,
  initialItemTraitMatrxObjectName = -1)
  # srand(1234)
  Random.seed!(1234)

  # Load training data
  trainingBasketsDict = load(trainingBasketsDictFileName,
                                  trainingBasketsDictObjectName)
  println("Loaded $trainingBasketsDictFileName")

  # Load test data
  testBasketsDict = load(testBasketsDictFileName,
                              testBasketsDictObjectName)
  println("Loaded $testBasketsDictFileName")

  # Build set of training instances
  numTrainingInstances = length(collect(keys(trainingBasketsDict)))
  trainingInstances = fill(Array{Int}(1), numTrainingInstances)
  trainingInstanceIndex = 1
  numItems = 0
  for trainingInstanceBasketId in collect(keys(trainingBasketsDict))
    trainingInstanceItems = deepcopy(trainingBasketsDict[trainingInstanceBasketId].basketItems)
    heldOutItem = trainingBasketsDict[trainingInstanceBasketId].heldOutItem
    push!(trainingInstanceItems, heldOutItem)

    trainingInstances[trainingInstanceIndex] = trainingInstanceItems

    trainingInstanceIndex += 1

    numItems = trainingBasketsDict[trainingInstanceBasketId].numItemsInCatalog
  end

  # Build set of test instances
  numTestInstances = length(collect(keys(testBasketsDict)))
  testInstances = fill(Array{Int}(1), numTestInstances)
  testInstanceIndex = 1
  for testInstanceBasketId in collect(keys(testBasketsDict))
    testInstanceItems = deepcopy(testBasketsDict[testInstanceBasketId].basketItems)
    heldOutItem = testBasketsDict[testInstanceBasketId].heldOutItem
    push!(testInstanceItems, heldOutItem)

    testInstances[testInstanceIndex] = testInstanceItems

    testInstanceIndex += 1
  end

  # Initialize itemTraitMatrix, if provided
  itemTraitMatrixInit = -1
  if initialItemTraitMatrxFileName != -1
    itemTraitMatrixInit = load(initialItemTraitMatrxFileName, initialItemTraitMatrxObjectName)
    println("Loaded $initialItemTraitMatrxFileName")
  end

  # Run sampler to learn the posterior on the low-rank DPP parameters (V matrix)
  runSamplerFunction(trainingInstances, numTrainingInstances, numItems,
    numItemTraits, testInstances, numTestInstances, learnedModelOutputDirName,
    itemTraitMatrixInit)
end


text:
using Documenter, SciMLBenchmarksOutput

makedocs(
    sitename="The SciML Benchmarks: Differential Equations, Inverse Problems, Physics-Informed ML, Science-Guided AI",
    authors="Chris Rackauckas",
    modules=[SciMLBenchmarksOutput],
    clean=true,doctest=false,
    format = Documenter.HTML(#analytics = "UA-90474609-3",
                             assets = ["assets/favicon.ico"],
                             canonical="https://benchmarks.sciml.ai/stable/"),
    pages=[
        "Home" => "README.md",
        "Non-Stiff Ordinary Differential Equation (ODE) Solver Benchmarks" => Any[
            "markdown/MultiLanguage/wrapper_packages.md"
        ],
    ]
)

deploydocs(;
    repo="github.com/SciML/SciMLBenchmarksOutput",
    devbranch="main",
    branch = "main"
)


text:
import Base: require

export @require

isprecompiling() = ccall(:jl_generating_output, Cint, ()) == 1

type Hook end

@init @guard begin
  methods(require).mt.cache.sig = Tuple{typeof(require),Symbol,Hook}
  function Base.require(mod::Symbol)
    eval(:(Base.require))(mod, Main.Requires.Hook())
    Main.Requires.loadmod(string(mod))
  end
end

loaded(mod) = getthing(Main, mod) != nothing

const modlisteners = Dict{AbstractString,Vector{Function}}()

listenmod(f, mod) =
  loaded(mod) ? f() :
    modlisteners[mod] = push!(get(modlisteners, mod, Function[]), f)

loadmod(mod) =
  map(f->f(), get(modlisteners, mod, []))

importexpr(mod::Symbol) = Expr(:import, mod)
importexpr(mod::Expr) = Expr(:import, map(Symbol, split(string(mod), "."))...)

function withpath(f, path)
  tls = task_local_storage()
  hassource = haskey(tls, :SOURCE_PATH)
  hassource && (path′ = tls[:SOURCE_PATH])
  tls[:SOURCE_PATH] = path
  try
    return f()
  finally
    hassource ?
      (tls[:SOURCE_PATH] = path′) :
      delete!(tls, :SOURCE_PATH)
  end
end

function err(f, listener, mod)
  try
    f()
  catch e
    warn("Error requiring $mod from $listener:")
    showerror(STDERR, e, catch_backtrace())
    println(STDERR)
  end
end

macro require(mod, expr)
  ex = quote
    listenmod($(string(mod))) do
      withpath(@__FILE__) do
        err($(current_module()), $(string(mod))) do
          $(esc(:(eval($(Expr(:quote, Expr(:block,
                                           importexpr(mod),
                                           expr)))))))
        end
      end
    end
  end
  quote
    if isprecompiling()
      @init @guard $(ex)
    else
      $(ex)
    end
  end
end


text:
module MethTest

# package code goes here

export test_greater_population
include("hyper-geometric.jl")

export get_beta_parameters
export getDCpGStatistic
export CpGMethCount
export CpGMethEstimate
include("beta-binomial-z.jl")

end # module


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ 83a2c2cc-296b-11eb-38bf-fd97aa3c6fda
using PlutoUI, NBodyLeapFrog, Plots, NiLang

# ╔═╡ fdafd8a6-29d0-11eb-0d05-73cb4a836a20
using NiLang.AD

# ╔═╡ c08d80cc-29d0-11eb-002d-d1e08da6a7c1
md"# Simulating Celestial Bodies"

# ╔═╡ 814c93d0-296c-11eb-333d-216ce303ce34
plotly();

# ╔═╡ bd337dc6-29d0-11eb-1a58-b354615e11e1
md"## Computing the orbitals with time-reversible Forest-Ruth integrator"

# ╔═╡ 8e419bc2-296b-11eb-0e0f-afdab78b1392
function showres(vs; plt=plot([], []))
    a = map(x->x.x, vs)
    b = map(x->x.y, vs)
	labels = ["sun", "mercury", "venus", "earth", "mars",
		"jupyter", "saturn", "uranus", "neptune", "pluto"]
    for i in 1:size(a, 1)
        plot!(plt,a[i,:],b[i,:]; aspect_ratio=:equal, label=labels[i])
    end
    plt
end

# ╔═╡ 4afeee16-29d4-11eb-0a02-ff725633872c
planets = Bodies.chunit_day2year.(Bodies.set)

# ╔═╡ a4110f0e-29d0-11eb-1cf8-a950830d0cfa
md"## Differentiate the simulation"

# ╔═╡ bbe14adc-296b-11eb-1cb2-c149ac6e6321
@i function loss!(y::T, r0::AbstractVector{V3{T}}, v0::AbstractVector{V3{T}}, m; n, dt) where T
	for i=1:n
		NBodyLeapFrog.i_fr_step!(r0, v0, m; dt=dt)
	end

	# compute the squared distance between Mercury (2) and earth (4)
	@routine begin
		sqy ← zero(T)
		sqy += NBodyLeapFrog.sqdistance(r0[9], r0[10])
	end
	y += sqrt(sqy)

	~@routine # uncompute sqy
end

# ╔═╡ de67315c-296b-11eb-0ae2-93757c2e2331
r0, v0, m = let
	planets = Bodies.chunit_day2year.(Bodies.set)
	getfield.(planets, :r), getfield.(planets, :v), getfield.(planets, :m)
end;

# ╔═╡ 36acacfe-29ce-11eb-18b3-0fc163ea0dd5
"""simulate and keep history"""
@i function simulate!(rmat::AbstractMatrix{V3{T}}, vmat::AbstractMatrix{V3{T}}, planets; dt) where T
	@safe @assert size(rmat) == size(vmat) && size(rmat, 1) == length(planets)

	# store the mass into a vector
	@routine begin
		m ← zeros(T, length(planets))
		for j=1:length(planets)
			m[j] += planets[j].m
		end
	end
	# setup initial values
	for j=1:length(planets)
		rmat[j,1] += planets[j].r
		vmat[j,1] += planets[j].v
	end

	# simulate!
	for i=1:size(rmat, 2)-1
		for j=1:length(planets)
			rmat[j,i+1] += rmat[j,i]
			vmat[j,i+1] += vmat[j,i]
		end
		# single step Forest-Ruth algorithm
		@inbounds NBodyLeapFrog.i_fr_step!(view(rmat,:,i+1), view(vmat,:,i+1), m; dt=dt)
	end

	~@routine #uncompute m
end

# ╔═╡ 157f8d8e-29cf-11eb-3926-01d791d6d703
rmat, vmat = let
	rmat = zeros(V3{Float64}, length(planets), 4001)
	vmat = zeros(V3{Float64}, length(planets), 4001)
	simulate!(rmat, vmat, planets; dt=0.01)
end;

# ╔═╡ 6529ba3a-29d4-11eb-1123-59951753b134
showres(rmat)

# ╔═╡ 6b217206-296c-11eb-0e56-a742db172508
y, r_, v_, m_ = loss!(0.0, copy(r0), copy(v0), copy(m); dt=0.01, n=4000);

# ╔═╡ ba521b86-296b-11eb-2bd4-ebf29f78b87a
_, _, gr, gv, gm = Grad(loss!)(Val(1), 0.0, copy(r0), copy(v0), copy(m); dt=0.01, n=4000); # Val(1) means, the the first argument stores the loss

# ╔═╡ 3863e6ea-29d1-11eb-3d9c-8baf25f581ca
gr

# ╔═╡ 0fa5ad5a-29d2-11eb-258f-590d950d565b
function showgrad(r0, gr)
    x0 = map(x->x.x, r0)
    y0 = map(x->x.y, r0)
    gx0 = map(x->x.x.g, gr)
    gy0 = map(x->x.y.g, gr)
	labels = ["sun", "mercury", "venus", "earth", "mars",
		"jupyter", "saturn", "uranus", "neptune", "pluto"]
    quiver(x0, y0; quiver=collect(zip(gx0, gy0)))
end

# ╔═╡ 88441af8-29d2-11eb-3e77-5bace9a3a1cf
let
	plt = showgrad(r0, gr)
	showres(rmat, plt=plt)
end

# ╔═╡ Cell order:
# ╟─c08d80cc-29d0-11eb-002d-d1e08da6a7c1
# ╠═83a2c2cc-296b-11eb-38bf-fd97aa3c6fda
# ╠═814c93d0-296c-11eb-333d-216ce303ce34
# ╟─bd337dc6-29d0-11eb-1a58-b354615e11e1
# ╠═8e419bc2-296b-11eb-0e0f-afdab78b1392
# ╠═36acacfe-29ce-11eb-18b3-0fc163ea0dd5
# ╠═4afeee16-29d4-11eb-0a02-ff725633872c
# ╠═157f8d8e-29cf-11eb-3926-01d791d6d703
# ╠═6529ba3a-29d4-11eb-1123-59951753b134
# ╟─a4110f0e-29d0-11eb-1cf8-a950830d0cfa
# ╠═fdafd8a6-29d0-11eb-0d05-73cb4a836a20
# ╠═bbe14adc-296b-11eb-1cb2-c149ac6e6321
# ╠═de67315c-296b-11eb-0ae2-93757c2e2331
# ╠═6b217206-296c-11eb-0e56-a742db172508
# ╠═ba521b86-296b-11eb-2bd4-ebf29f78b87a
# ╠═3863e6ea-29d1-11eb-3d9c-8baf25f581ca
# ╠═0fa5ad5a-29d2-11eb-258f-590d950d565b
# ╠═88441af8-29d2-11eb-3e77-5bace9a3a1cf


text:
# Data taken from
# https://en.wikipedia.org/wiki/Electron_configurations_of_the_elements_(data_page)
const table_of_elements = [:H => ("hydrogen", c"1s1"),
                           :He => ("helium", c"1s2"),
                           :Li => ("lithium", c"[He]* 2s1"),
                           :Be => ("beryllium", c"[He]* 2s2"),
                           :B => ("boron", c"[He]* 2s2 2p1"),
                           :C => ("carbon", c"[He]* 2s2 2p2"),
                           :N => ("nitrogen", c"[He]* 2s2 2p3"),
                           :O => ("oxygen", c"[He]* 2s2 2p4"),
                           :F => ("fluorine", c"[He]* 2s2 2p5"),
                           :Ne => ("neon", c"[He]* 2s2 2p6"),
                           :Na => ("sodium", c"[Ne]* 3s1"),
                           :Mg => ("magnesium", c"[Ne]* 3s2"),
                           :Al => ("aluminium", c"[Ne]* 3s2 3p1"),
                           :Si => ("silicon", c"[Ne]* 3s2 3p2"),
                           :P => ("phosphorus", c"[Ne]* 3s2 3p3"),
                           :S => ("sulfur", c"[Ne]* 3s2 3p4"),
                           :Cl => ("chlorine", c"[Ne]* 3s2 3p5"),
                           :Ar => ("argon", c"[Ne]* 3s2 3p6"),
                           :K => ("potassium", c"[Ar]* 4s1"),
                           :Ca => ("calcium", c"[Ar]* 4s2"),
                           :Sc => ("scandium", c"[Ar]* 3d1 4s2"),
                           :Ti => ("titanium", c"[Ar]* 3d2 4s2"),
                           :V => ("vanadium", c"[Ar]* 3d3 4s2"),
                           :Cr => ("chromium", c"[Ar]* 3d5 4s1"),
                           :Mn => ("manganese", c"[Ar]* 3d5 4s2"),
                           :Fe => ("iron", c"[Ar]* 3d6 4s2"),
                           :Co => ("cobalt", c"[Ar]* 3d7 4s2"),
                           :Ni => ("nickel", c"[Ar]* 3d8 4s2"), # or [Ar]* 3d9 4s1(disputed - see Nickel article)
                           :Cu => ("copper", c"[Ar]* 3d10 4s1"),
                           :Zn => ("zinc", c"[Ar]* 3d10 4s2"),
                           :Ga => ("gallium", c"[Ar]* 3d10 4s2 4p1"),
                           :Ge => ("germanium", c"[Ar]* 3d10 4s2 4p2"),
                           :As => ("arsenic", c"[Ar]* 3d10 4s2 4p3"),
                           :Se => ("selenium", c"[Ar]* 3d10 4s2 4p4"),
                           :Br => ("bromine", c"[Ar]* 3d10 4s2 4p5"),
                           :Kr => ("krypton", c"[Ar]* 3d10 4s2 4p6"),
                           :Rb => ("rubidium", c"[Kr]* 5s1"),
                           :Sr => ("strontium", c"[Kr]* 5s2"),
                           :Y => ("yttrium", c"[Kr]* 4d1 5s2"),
                           :Zr => ("zirconium", c"[Kr]* 4d2 5s2"),
                           :Nb => ("niobium", c"[Kr]* 4d4 5s1"),
                           :Mo => ("molybdenum", c"[Kr]* 4d5 5s1"),
                           :Tc => ("technetium", c"[Kr]* 4d5 5s2"),
                           :Ru => ("ruthenium", c"[Kr]* 4d7 5s1"),
                           :Rh => ("rhodium", c"[Kr]* 4d8 5s1"),
                           :Pd => ("palladium", c"[Kr]* 4d10"),
                           :Ag => ("silver", c"[Kr]* 4d10 5s1"),
                           :Cd => ("cadmium", c"[Kr]* 4d10 5s2"),
                           :In => ("indium", c"[Kr]* 4d10 5s2 5p1"),
                           :Sn => ("tin", c"[Kr]* 4d10 5s2 5p2"),
                           :Sb => ("antimony", c"[Kr]* 4d10 5s2 5p3"),
                           :Te => ("tellurium", c"[Kr]* 4d10 5s2 5p4"),
                           :I => ("iodine", c"[Kr]* 4d10 5s2 5p5"),
                           :Xe => ("xenon", c"[Kr]* 4d10 5s2 5p6"),
                           :Cs => ("caesium", c"[Xe]* 6s1"),
                           :Ba => ("barium", c"[Xe]* 6s2"),
                           :La => ("lanthanum", c"[Xe]* 5d1 6s2"),
                           :Ce => ("cerium", c"[Xe]* 4f1 5d1 6s2"),
                           :Pr => ("praseodymium", c"[Xe]* 4f3 6s2"),
                           :Nd => ("neodymium", c"[Xe]* 4f4 6s2"),
                           :Pm => ("promethium", c"[Xe]* 4f5 6s2"),
                           :Sm => ("samarium", c"[Xe]* 4f6 6s2"),
                           :Eu => ("europium", c"[Xe]* 4f7 6s2"),
                           :Gd => ("gadolinium", c"[Xe]* 4f7 5d1 6s2"),
                           :Tb => ("terbium", c"[Xe]* 4f9 6s2"),
                           :Dy => ("dysprosium", c"[Xe]* 4f10 6s2"),
                           :Ho => ("holmium", c"[Xe]* 4f11 6s2"),
                           :Er => ("erbium", c"[Xe]* 4f12 6s2"),
                           :Tm => ("thulium", c"[Xe]* 4f13 6s2"),
                           :Yb => ("ytterbium", c"[Xe]* 4f14 6s2"),
                           :Lu => ("lutetium", c"[Xe]* 4f14 5d1 6s2"),
                           :Hf => ("hafnium", c"[Xe]* 4f14 5d2 6s2"),
                           :Ta => ("tantalum", c"[Xe]* 4f14 5d3 6s2"),
                           :W => ("tungsten", c"[Xe]* 4f14 5d4 6s2"),
                           :Re => ("rhenium", c"[Xe]* 4f14 5d5 6s2"),
                           :Os => ("osmium", c"[Xe]* 4f14 5d6 6s2"),
                           :Ir => ("iridium", c"[Xe]* 4f14 5d7 6s2"),
                           :Pt => ("platinum", c"[Xe]* 4f14 5d9 6s1"),
                           :Au => ("gold", c"[Xe]* 4f14 5d10 6s1"),
                           :Hg => ("mercury", c"[Xe]* 4f14 5d10 6s2"),
                           :Tl => ("thallium", c"[Xe]* 4f14 5d10 6s2 6p1"),
                           :Pb => ("lead", c"[Xe]* 4f14 5d10 6s2 6p2"),
                           :Bi => ("bismuth", c"[Xe]* 4f14 5d10 6s2 6p3"),
                           :Po => ("polonium", c"[Xe]* 4f14 5d10 6s2 6p4"),
                           :At => ("astatine", c"[Xe]* 4f14 5d10 6s2 6p5"),
                           :Rn => ("radon", c"[Xe]* 4f14 5d10 6s2 6p6"),
                           :Fr => ("francium", c"[Rn]* 7s1"),
                           :Ra => ("radium", c"[Rn]* 7s2"),
                           :Ac => ("actinium", c"[Rn]* 6d1 7s2"),
                           :Th => ("thorium", c"[Rn]* 6d2 7s2"),
                           :Pa => ("protactinium", c"[Rn]* 5f2 6d1 7s2"),
                           :U => ("uranium", c"[Rn]* 5f3 6d1 7s2"),
                           :Np => ("neptunium", c"[Rn]* 5f4 6d1 7s2"),
                           :Pu => ("plutonium", c"[Rn]* 5f6 7s2"),
                           :Am => ("americium", c"[Rn]* 5f7 7s2"),
                           :Cm => ("curium", c"[Rn]* 5f7 6d1 7s2"),
                           :Bk => ("berkelium", c"[Rn]* 5f9 7s2"),
                           :Cf => ("californium", c"[Rn]* 5f10 7s2"),
                           :Es => ("einsteinium", c"[Rn]* 5f11 7s2"),
                           :Fm => ("fermium", c"[Rn]* 5f12 7s2"),
                           :Md => ("mendelevium", c"[Rn]* 5f13 7s2"),
                           :No => ("nobelium", c"[Rn]* 5f14 7s2"),
                           :Lr => ("lawrencium", c"[Rn]* 5f14 7s2 7p1"),
                           :Rf => ("rutherfordium", c"[Rn]* 5f14 6d2 7s2"),
                           :Db => ("dubnium", c"[Rn]* 5f14 6d3 7s2"),
                           :Sg => ("seaborgium", c"[Rn]* 5f14 6d4 7s2"),
                           :Bh => ("bohrium", c"[Rn]* 5f14 6d5 7s2"),
                           :Hs => ("hassium", c"[Rn]* 5f14 6d6 7s2"),
                           :Mt => ("meitnerium", c"[Rn]* 5f14 6d7 7s2"),
                           :Ds => ("darmstadtium", c"[Rn]* 5f14 6d8 7s2"),
                           :Rg => ("roentgenium", c"[Rn]* 5f14 6d9 7s2"),
                           :Cn => ("copernicium", c"[Rn]* 5f14 6d10 7s2"),
                           :Nh => ("nihonium", c"[Rn]* 5f14 6d10 7s2 7p1"),
                           :Fl => ("flerovium", c"[Rn]* 5f14 6d10 7s2 7p2"),
                           :Mc => ("moscovium", c"[Rn]* 5f14 6d10 7s2 7p3"),
                           :Lv => ("livermorium", c"[Rn]* 5f14 6d10 7s2 7p4"),
                           :Ts => ("tennessine", c"[Rn]* 5f14 6d10 7s2 7p5"),
                           :Og => ("oganesson", c"[Rn]* 5f14 6d10 7s2 7p6")]


text:
## primitive experiment of

#read prop axiom
cp1=readcore("data/prop001.cnf")
#add template
alltemplateof(cp1)
#show it
printcore(cp1)

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
using FFTW

using Oceananigans.BoundaryConditions: PBC, ZFBC

#=
These functions return the transforms required to solve Poisson's equation with
periodic boundary conditions or staggered Neumann boundary.

Fast Fourier transforms (FFTs) are used in the periodic dimensions and
real-to-real discrete cosine transforms are used in the wall-bounded dimensions.
Note that the DCT-II is used for the DCT and the DCT-III for the IDCT
which correspond to REDFT10 and REDFT01 in FFTW.

They operatore on an array with the shape of `A`, which is needed to plan
efficient transforms. `A` will be mutated.
=#

plan_forward_transform(A::Array, ::PBC, dims, planner_flag=FFTW.PATIENT) =
    FFTW.plan_fft!(A, dims, flags=planner_flag)

plan_forward_transform(A::Array, ::ZFBC, dims, planner_flag=FFTW.PATIENT) =
    FFTW.plan_r2r!(A, FFTW.REDFT10, dims, flags=planner_flag)

plan_backward_transform(A::Array, ::PBC, dims, planner_flag=FFTW.PATIENT) =
    FFTW.plan_ifft!(A, dims, flags=planner_flag)

plan_backward_transform(A::Array, ::ZFBC, dims, planner_flag=FFTW.PATIENT) =
    FFTW.plan_r2r!(A, FFTW.REDFT01, dims, flags=planner_flag)

@hascuda begin
     plan_forward_transform(A::CuArray, ::Union{PBC, ZFBC}, dims) = plan_fft!(A, dims)
    plan_backward_transform(A::CuArray, ::Union{PBC, ZFBC}, dims) = plan_ifft!(A, dims)
end


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

mutable struct IoK8sApiCoreV1PodReadinessGate <: SwaggerModel
    conditionType::Any # spec type: Union{ Nothing, String } # spec name: conditionType

    function IoK8sApiCoreV1PodReadinessGate(;conditionType=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodReadinessGate, Symbol("conditionType"), conditionType)
        setfield!(o, Symbol("conditionType"), conditionType)
        o
    end
end # type IoK8sApiCoreV1PodReadinessGate

const _property_map_IoK8sApiCoreV1PodReadinessGate = Dict{Symbol,Symbol}(Symbol("conditionType")=>Symbol("conditionType"))
const _property_types_IoK8sApiCoreV1PodReadinessGate = Dict{Symbol,String}(Symbol("conditionType")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1PodReadinessGate }) = collect(keys(_property_map_IoK8sApiCoreV1PodReadinessGate))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodReadinessGate }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1PodReadinessGate[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodReadinessGate }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodReadinessGate[property_name]

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

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


text:
using Revise
using FEniCS
using PyCall
using PyPlot
using Infiltrator
using Bem2d

# Evaluation solution at a new point.  T. Ben Thomposon wrote this!
function interior_eval(u, x, y)
    py"""
    def G(u, x, y):
        return u(x, y)
    """
    py"G"(u.pyobject, x, y)
end

# Calculate strain
function strain(u)
    0.5 * (nabla_grad(u) + Transpose(nabla_grad(u)))
end

# Calculate
function stress(u, lambda, mu)
     lambda * nabla_div(u) * Identity(2) + 2 * mu * strain(u)
end

function femvsbemgravitybox()
    mu = 3e10
    lambda = mu
    rho = 2700
    g = 9.81
    width = 10e3

    ###
    ### FEM solution
    ###
    # Create mesh and define function space
    mesh = RectangleMesh(Point((-width, 0)), Point((width, 2 * width)), 100, 100)
    V = VectorFunctionSpace(mesh, "P", 1)
    bc = DirichletBC(V, Constant((0, 0)), "on_boundary && x[1]<1E-14") # What BCs are bing set???

    # Solve
    u = TrialFunction(V)
    d = geometric_dimension(u) # space dimension
    v = TestFunction(V)
    f = Constant((0, -rho * g)) # Vector of uniform body force
    T = Constant((0, 0))
    a = inner(stress(u, lambda, mu), strain(v)) * dx
    L = dot(f, v) * dx + dot(T, v) * ds
    u = FeFunction(V)
    lvsolve(a, L, u, bc)

    # Extract plottable quantities
    V = FunctionSpace(mesh, "P", 1)
    u_magnitude = sqrt(dot(u, u))
    u_magnitude = project(u_magnitude, V)
    close("all")
    # figure(figsize=(10, 5))
    # subplot(1, 3, 1)
    # plothandle = FEniCS.Plot(u_magnitude)
    # colorbar(plothandle)
    # title("FEM solution")

    ###
    ### BEM solution
    ###
    alpha = 7e-8 # scalar preconditioner for traction terms
    nu = 0.25
    nels = 20
    npts = 100
    offset = 10
    x, y = obsgrid(-width+offset, -2*width+offset, width-offset, 0-offset, npts)

    # Define BEM geometry
    els = Elements(Int(1e5))
    x1, y1, x2, y2 = discretizedline(-width, -2*width, width, -2*width, nels) # Bottom
    addelsez!(els, x1, y1, x2, y2, "B")
    x1, y1, x2, y2 = discretizedline(width, -2*width, width, 0, nels) # Right hand side
    addelsez!(els, x1, y1, x2, y2, "R")
    x1, y1, x2, y2 = discretizedline(width, 0, -width, 0, nels) # Top
    addelsez!(els, x1, y1, x2, y2, "T")
    x1, y1, x2, y2 = discretizedline(-width, 0, -width, -2*width, nels) # Left hand side
    addelsez!(els, x1, y1, x2, y2, "L")

    # Common indexing
    idx = getidxdict(els) # Should this return "all" - YES TODO
    bcidxU = idx["B"] # Boundaries with *displacement* BCs
    bcidxT = [idx["R"] ; idx["T"]; idx["L"]] # Boundaries with *traction* BCs
    bcidxall = collect(1:1:els.endidx) # All boundaries

    # Gravity square problem with quadratic elements
    T_pU_qall, _ = PUTQ(slip2dispstress, els, bcidxU, bcidxall, mu, nu)
    _, H_pT_qall = PUTQ(slip2dispstress, els, bcidxT, bcidxall, mu, nu)
    TH = [T_pU_qall ; alpha .* H_pT_qall] # Assemble combined linear operator

    # Particular solution and effective boundary conditions
    xnodes = transpose(els.xnodes[idx["B"], :])[:]
    ynodes = transpose(els.ynodes[idx["B"], :])[:]
    UB, _ = gravityparticularfunctions(xnodes, ynodes, g, rho, lambda, mu)
    bcs = zeros(6 * els.endidx)
    bcs[1:2:6*nels] = UB[:, 1] # Bottom boundary (x-component)
    bcs[2:2:6*nels] = UB[:, 2] # Bottom boundary (y-component)
    bcs *= -1 # This is neccesary for the right answer and is consistent with derivation

    # BEM solve to get particular solution
    Ueffparticular = inv(TH) * bcs

    Uinteriorcomplementary, Sinteriorcomplementary = quaddispstress(slip2dispstress, x, y, els, bcidxall, quadstack(Ueffparticular[1:2:end]), quadstack(Ueffparticular[2:2:end]), mu, nu)
    Uinteriorparticular, Sinteriorparticular = gravityparticularfunctions(x, y, g, rho, lambda, mu)
    U = @. Uinteriorcomplementary + Uinteriorparticular
    S = @. Sinteriorcomplementary + Sinteriorparticular
    Umag = sqrt.(U[:, 1].^2 + U[:, 2].^2)

    figure(figsize=(16, 4))
    subplot(1, 3, 1)
    contourf(reshape(x, npts, npts), reshape(y, npts, npts).+2*width, reshape(Umag, npts, npts), 50)
    colorbar()
    gca().set_xlim([-width, width])
    gca().set_ylim([0, 2*width])
    gca().set_aspect("equal")
    title("BEM")

    # Evaluate FEM solution at BEM observation coordinates
    fematbemux = zeros(length(x))
    fematbemuy = zeros(length(x))
    for i in 1:length(x)
        fematbemux[i], fematbemuy[i] = interior_eval(u, x[i], y[i] + 2 * width)
    end
    fematbemumag = sqrt.(fematbemux.^2 + fematbemuy.^2)
    fematbemumagresid = sqrt.((fematbemux.-U[:, 1]).^2 + (fematbemuy.-U[:, 2]).^2)

    subplot(1, 3, 2)
    contourf(reshape(x, npts, npts), reshape(y, npts, npts).+2*width,
             reshape(fematbemumag, npts, npts), 50)
    colorbar()
    gca().set_xlim([-width, width])
    gca().set_ylim([0, 2*width])
    gca().set_aspect("equal")
    title("FEM")

    subplot(1, 3, 3)
    contourf(reshape(x, npts, npts), reshape(y, npts, npts).+2*width,
             reshape(log10.(abs.(fematbemumagresid)), npts, npts), 50)
    colorbar()
    gca().set_xlim([-width, width])
    gca().set_ylim([0, 2*width])
    gca().set_aspect("equal")
    title("log10(BEM - FEM)")

end
femvsbemgravitybox()


text:
module Day22
__revise_mode__ = :eval

using Advent: save_input, eachlinegroup

function read_input(filename)
    Vector{Int}[
        [parse(Int, x) for x in Iterators.drop(group, 1)]
        for group in eachlinegroup(filename)
    ]
end

deck_score(deck) = sum(prod, deck |> Iterators.reverse |> enumerate)

function solve_a(filename)
    groups = read_input(filename)
    while !any(isempty, groups)
        nums = popfirst!.(groups)
        low, high = extrema(nums)
        push!(groups[argmax(nums)], high, low)
    end

    Iterators.filter(!isempty, groups) |> only |> deck_score
end

function play_game!(groups)
    seen = Set{UInt64}()
    while !any(isempty, groups)
        h = reduce(hash, Iterators.map(hash, groups))
        h in seen && return 1
        push!(seen, h)

        nums = popfirst!.(groups)

        if all(length(groups[i]) >= n for (i, n) in enumerate(nums))
            winner = play_game!([groups[i][1:n] for (i, n) in enumerate(nums)])
            push!(groups[winner], popat!(nums, winner), only(nums))
        else
            low, high = extrema(nums)
            push!(groups[argmax(nums)], high, low)
        end
    end
    findmax(groups)[2]
end

function solve_b(filename)
    groups = read_input(filename)
    winner = play_game!(groups)
    deck_score(groups[winner])
end

solve(filename = save_input(2020, 22)) = solve_a(filename), solve_b(filename)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive