text:
import Base: +, -, *, /, \, exp, log, ≈, ==, inv, conj

"""
    UnitQuaternion{T} <: Rotation

4-parameter attitute representation that is singularity-free. Quaternions with unit norm
represent a double-cover of SO(3). The `UnitQuaternion` does NOT strictly enforce the unit
norm constraint, but certain methods will assume you have a unit quaternion.
Follows the Hamilton convention for quaternions.

# Constructors
```julia
UnitQuaternion(w,x,y,z)
UnitQuaternion(q::AbstractVector)
```
where `w` is the scalar (real) part, `x`,`y`, and `z` are the vector (imaginary) part,
and `q = [w,x,y,z]`.
"""
struct UnitQuaternion{T} <: Rotation{3,T}
    w::T
    x::T
    y::T
    z::T

    @inline function UnitQuaternion{T}(w, x, y, z, normalize::Bool = true) where T
        if normalize
            inorm = inv(sqrt(w*w + x*x + y*y + z*z))
            new{T}(w*inorm, x*inorm, y*inorm, z*inorm)
        else
            new{T}(w, x, y, z)
        end
    end

    UnitQuaternion{T}(q::UnitQuaternion) where T = new{T}(q.w, q.x, q.y, q.z)
end

# ~~~~~~~~~~~~~~~ Constructors ~~~~~~~~~~~~~~~ #
# Use default map
function UnitQuaternion(w,x,y,z, normalize::Bool = true)
    types = promote(w,x,y,z)
    UnitQuaternion{eltype(types)}(w,x,y,z, normalize)
end

# Pass in Vectors
@inline function (::Type{Q})(q::AbstractVector, normalize::Bool = true) where Q <: UnitQuaternion
    check_length(q, 4)
    Q(q[1], q[2], q[3], q[4], normalize)
end
@inline (::Type{Q})(q::StaticVector{4}, normalize::Bool = true) where Q <: UnitQuaternion =
    Q(q[1], q[2], q[3], q[4], normalize)

# Copy constructors
UnitQuaternion(q::UnitQuaternion) = q

# UnitQuaternion <=> Quat
# (::Type{Q})(q::Quat) where Q <: UnitQuaternion = Q(q.w, q.x, q.y, q.z, false)
# (::Type{Q})(q::UnitQuaternion) where Q <: Quat = Q(q.w, q.x, q.y, q.z, false)
# const AllQuats{T} = Union{<:Quat{T}, <:UnitQuaternion{T}}

# ~~~~~~~~~~~~~~~ StaticArrays Interface ~~~~~~~~~~~~~~~ #
function (::Type{Q})(t::NTuple{9}) where Q<:UnitQuaternion
    #=
    This function solves the system of equations in Section 3.1
    of https://arxiv.org/pdf/math/0701759.pdf. This cheap method
    only works for matrices that are already orthonormal (orthogonal
    and unit length columns). The nearest orthonormal matrix can
    be found by solving Wahba's problem:
    https://en.wikipedia.org/wiki/Wahba%27s_problem as shown below.

    not_orthogonal = randn(3,3)
    u,s,v = svd(not_orthogonal)
    is_orthogoral = u * diagm([1, 1, sign(det(u * transpose(v)))]) * transpose(v)
    =#

    a = 1 + t[1] + t[5] + t[9]
    b = 1 + t[1] - t[5] - t[9]
    c = 1 - t[1] + t[5] - t[9]
    d = 1 - t[1] - t[5] + t[9]
    max_abcd = max(a, b, c, d)
    if a == max_abcd
        b = t[6] - t[8]
        c = t[7] - t[3]
        d = t[2] - t[4]
    elseif b == max_abcd
        a = t[6] - t[8]
        c = t[2] + t[4]
        d = t[7] + t[3]
    elseif c == max_abcd
        a = t[7] - t[3]
        b = t[2] + t[4]
        d = t[6] + t[8]
    else
        a = t[2] - t[4]
        b = t[7] + t[3]
        c = t[6] + t[8]
    end
    return Q(a, b, c, d)
end

function Base.getindex(q::UnitQuaternion, i::Int)
    if i == 1
        ww = (q.w * q.w)
        xx = (q.x * q.x)
        yy = (q.y * q.y)
        zz = (q.z * q.z)

        ww + xx - yy - zz
    elseif i == 2
        xy = (q.x * q.y)
        zw = (q.w * q.z)

        2 * (xy + zw)
    elseif i == 3
        xz = (q.x * q.z)
        yw = (q.y * q.w)

        2 * (xz - yw)
    elseif i == 4
        xy = (q.x * q.y)
        zw = (q.w * q.z)

        2 * (xy - zw)
    elseif i == 5
        ww = (q.w * q.w)
        xx = (q.x * q.x)
        yy = (q.y * q.y)
        zz = (q.z * q.z)

        ww - xx + yy - zz
    elseif i == 6
        yz = (q.y * q.z)
        xw = (q.w * q.x)

        2 * (yz + xw)
    elseif i == 7
        xz = (q.x * q.z)
        yw = (q.y * q.w)

        2 * (xz + yw)
    elseif i == 8
        yz = (q.y * q.z)
        xw = (q.w * q.x)

        2 * (yz - xw)
    elseif i == 9
        ww = (q.w * q.w)
        xx = (q.x * q.x)
        yy = (q.y * q.y)
        zz = (q.z * q.z)

        ww - xx - yy + zz
    else
        throw(BoundsError(r,i))
    end
end

function Base.Tuple(q::UnitQuaternion)
    ww = (q.w * q.w)
    xx = (q.x * q.x)
    yy = (q.y * q.y)
    zz = (q.z * q.z)
    xy = (q.x * q.y)
    zw = (q.w * q.z)
    xz = (q.x * q.z)
    yw = (q.y * q.w)
    yz = (q.y * q.z)
    xw = (q.w * q.x)

    # initialize rotation part
    return (ww + xx - yy - zz,
            2 * (xy + zw),
            2 * (xz - yw),
            2 * (xy - zw),
            ww - xx + yy - zz,
            2 * (yz + xw),
            2 * (xz + yw),
            2 * (yz - xw),
            ww - xx - yy + zz)
end

# ~~~~~~~~~~~~~~~ Getters ~~~~~~~~~~~~~~~ #
@inline scalar(q::UnitQuaternion) = q.w
@inline vector(q::UnitQuaternion) = SVector{3}(q.x, q.y, q.z)

"""
    params(R::Rotation)

Return an `SVector` of the underlying parameters used by the rotation representation.

# Example
```julia
p = MRP(1.0, 2.0, 3.0)
Rotations.params(p) == @SVector [1.0, 2.0, 3.0]  # true
```
"""
@inline params(q::UnitQuaternion) = SVector{4}(q.w, q.x, q.y, q.z)

# ~~~~~~~~~~~~~~~ Initializers ~~~~~~~~~~~~~~~ #
Base.rand(::Type{<:UnitQuaternion{T}}) where T =
    normalize(UnitQuaternion{T}(randn(T), randn(T), randn(T), randn(T)))
Base.rand(::Type{UnitQuaternion}) = Base.rand(UnitQuaternion{Float64})
@inline Base.zero(::Type{Q}) where Q <: UnitQuaternion = Q(1.0, 0.0, 0.0, 0.0)
@inline Base.one(::Type{Q}) where Q <: UnitQuaternion = Q(1.0, 0.0, 0.0, 0.0)

# ~~~~~~~~~~~~~~~ Math Operations ~~~~~~~~~~~~~~~ #

# Inverses
conj(q::Q) where Q <: UnitQuaternion = Q(q.w, -q.x, -q.y, -q.z, false)
inv(q::UnitQuaternion) = conj(q)
(-)(q::Q) where Q <: UnitQuaternion = Q(-q.w, -q.x, -q.y, -q.z, false)

# Norms
LinearAlgebra.norm(q::UnitQuaternion) = sqrt(q.w^2 + q.x^2 + q.y^2 + q.z^2)
vecnorm(q::UnitQuaternion) = sqrt(q.x^2 + q.y^2 + q.z^2)

function LinearAlgebra.normalize(q::Q) where Q <: UnitQuaternion
    n = inv(norm(q))
    Q(q.w*n, q.x*n, q.y*n, q.z*n)
end

# Identity
(::Type{Q})(I::UniformScaling) where Q <: UnitQuaternion = one(Q)

# Exponentials and Logarithms
"""
    pure_quaternion(v::AbstractVector)
    pure_quaternion(x, y, z)

Create a `UnitQuaternion` with zero scalar part (i.e. `q.w == 0`).
"""
function pure_quaternion(v::AbstractVector)
    check_length(v, 3)
    UnitQuaternion(zero(eltype(v)), v[1], v[2], v[3], false)
end

@inline pure_quaternion(x::Real, y::Real, z::Real) =
    UnitQuaternion(zero(x), x, y, z, false)

function exp(q::Q) where Q <: UnitQuaternion
    θ = vecnorm(q)
    sθ,cθ = sincos(θ)
    es = exp(q.w)
    M = es*sθ/θ
    Q(es*cθ, q.x*M, q.y*M, q.z*M, false)
end

function expm(ϕ::AbstractVector)
    check_length(ϕ, 3)
    θ = norm(ϕ)
    sθ,cθ = sincos(θ/2)
    M = 1//2 *sinc(θ/π/2)
    UnitQuaternion(cθ, ϕ[1]*M, ϕ[2]*M, ϕ[3]*M, false)
end

function log(q::Q, eps=1e-6) where Q <: UnitQuaternion
    # Assumes unit quaternion
    θ = vecnorm(q)
    if θ > eps
        M = atan(θ, q.w)/θ
    else
        M = (1-(θ^2/(3q.w^2)))/q.w
    end
    pure_quaternion(M*vector(q))
end

function logm(q::UnitQuaternion)
    # Assumes unit quaternion
    2*vector(log(q))
end

# Composition
"""
    (*)(q::UnitQuaternion, w::UnitQuaternion)

Quternion Composition

Equivalent to
```julia
lmult(q) * SVector(w)
rmult(w) * SVector(q)
```

Sets the output mapping equal to the mapping of `w`
"""
function (*)(q::UnitQuaternion, w::UnitQuaternion)
    UnitQuaternion(q.w * w.w - q.x * w.x - q.y * w.y - q.z * w.z,
                   q.w * w.x + q.x * w.w + q.y * w.z - q.z * w.y,
                   q.w * w.y - q.x * w.z + q.y * w.w + q.z * w.x,
                   q.w * w.z + q.x * w.y - q.y * w.x + q.z * w.w, false)
end

"""
    (*)(q::UnitQuaternion, r::StaticVector)

Rotate a vector

Equivalent to `hmat()' lmult(q) * rmult(q)' hmat() * r`
"""
function Base.:*(q::UnitQuaternion, r::StaticVector)  # must be StaticVector to avoid ambiguity
    check_length(r, 3)
    w = q.w
    v = vector(q)
    (w^2 - v'v)*r + 2*v*(v'r) + 2*w*cross(v,r)
end

"""
    (*)(q::UnitQuaternion, w::Real)

Scalar multiplication of a quaternion. Breaks unit norm.
"""
function (*)(q::Q, w::Real) where Q<:UnitQuaternion
    return Q(q.w*w, q.x*w, q.y*w, q.z*w, false)
end
(*)(w::Real, q::UnitQuaternion) = q*w

"""
    (/)(q::UnitQuaternion, w::Real)

Scalar divison of a quaternion. Breaks unit norm.
"""
function (/)(q::Q, w::Real) where Q<:UnitQuaternion
    return Q(q.w/w, q.x/w, q.y/w, q.z/w, false)
end
(\)(w::Real, q::UnitQuaternion) = q/w
(/)(w::Real, q::UnitQuaternion) = w*conj(q)/norm(q)^2      # Equivalent to w*conj(q1) if q1 has unit length
(\)(q::UnitQuaternion, w::Real) = w/q

(\)(q1::UnitQuaternion, q2::UnitQuaternion) = conj(q1)*q2  # Equivalent to inv(q1)*q2
(/)(q1::UnitQuaternion, q2::UnitQuaternion) = q1*conj(q2)  # Equivalent to q1*inv(q2)

(\)(q::UnitQuaternion, r::SVector{3}) = conj(q)*r          # Equivalent to inv(q)*r

"""
    rotation_between(from, to)

Compute the quaternion that rotates vector `from` so that it aligns with vector
`to`, along the geodesic (shortest path).
"""
rotation_between(from::AbstractVector, to::AbstractVector) = rotation_between(SVector{3}(from), SVector{3}(to))
function rotation_between(from::SVector{3}, to::SVector{3})
    # Robustified version of implementation from https://www.gamedev.net/topic/429507-finding-the-quaternion-betwee-two-vectors/#entry3856228
    normprod = sqrt(dot(from, from) * dot(to, to))
    T = typeof(normprod)
    normprod < eps(T) && throw(ArgumentError("Input vectors must be nonzero."))
    w = normprod + dot(from, to)
    v = abs(w) < 100 * eps(T) ? perpendicular_vector(from) : cross(from, to)
    @inbounds return UnitQuaternion(w, v[1], v[2], v[3]) # relies on normalization in constructor
end

# ~~~~~~~~~~~~~~~ Kinematics ~~~~~~~~~~~~~~~ $
"""
    kinematics(R::Rotation{3}, ω::AbstractVector)

The time derivative of the rotation R, according to the definition

``Ṙ = \\lim_{Δt → 0} \\frac{R(t + Δt) - R(t)}{Δt}``

where `ω` is the angular velocity. This is equivalent to

``Ṙ = \\lim_{Δt → 0} \\frac{R δR - R}{Δt}``

where ``δR`` is some small rotation, parameterized by a small rotation ``δθ`` about
an axis ``r``, such that ``lim_{Δt → 0} \\frac{δθ r}{Δt} = ω``

The kinematics are extremely useful when computing the dynamics of rigid bodies, since
`Ṙ = kinematics(R,ω)` is the first-order ODE for the evolution of the attitude dynamics.

See "Fundamentals of Spacecraft Attitude Determination and Control" by Markley and Crassidis
Sections 3.1-3.2 for more details.
"""
function kinematics(q::Q, ω::AbstractVector) where Q <: UnitQuaternion
    1//2 * params(q*Q(0.0, ω[1], ω[2], ω[3], false))
end

# ~~~~~~~~~~~~~~~ Linear Algebraic Conversions ~~~~~~~~~~~~~~~ #
"""
    lmult(q::UnitQuaternion)
    lmult(q::StaticVector{4})

`lmult(q2)*params(q1)` returns a vector equivalent to `q2*q1` (quaternion composition)
"""
function lmult(q::UnitQuaternion)
    SA[
        q.w -q.x -q.y -q.z;
        q.x  q.w -q.z  q.y;
        q.y  q.z  q.w -q.x;
        q.z -q.y  q.x  q.w;
    ]
end
lmult(q::StaticVector{4}) = lmult(UnitQuaternion(q, false))

"""
    rmult(q::UnitQuaternion)
    rmult(q::StaticVector{4})

`rmult(q1)*params(q2)` return a vector equivalent to `q2*q1` (quaternion composition)
"""
function rmult(q::UnitQuaternion)
    SA[
        q.w -q.x -q.y -q.z;
        q.x  q.w  q.z -q.y;
        q.y -q.z  q.w  q.x;
        q.z  q.y -q.x  q.w;
    ]
end
rmult(q::SVector{4}) = rmult(UnitQuaternion(q, false))

"""
    tmat()

`tmat()*params(q)`return a vector equivalent to `inv(q)`, where `q` is a `UnitQuaternion`
"""
function tmat(::Type{T}=Float64) where T
    SA{T}[
        1  0  0  0;
        0 -1  0  0;
        0  0 -1  0;
        0  0  0 -1;
    ]
end

"""
    vmat()

`vmat()*params(q)`` returns the imaginary
    (vector) part of the quaternion `q` (equivalent to `vector(q)``)
"""
function vmat(::Type{T}=Float64) where T
    SA{T}[
        0 1 0 0;
        0 0 1 0;
        0 0 0 1
    ]
end

"""
    hmat()
    hmat(r::AbstractVector)

`hmat()*r` or `hmat(r)` converts `r` into a pure quaternion, where `r` is 3-dimensional.

`hmat() == vmat()'`
"""
function hmat(::Type{T}=Float64) where T
    SA{T}[
        0 0 0;
        1 0 0;
        0 1 0;
        0 0 1.;
    ]
end

function hmat(r)
    @assert length(r) == 3
    SA[0, r[1], r[2], r[3]]
end

# ~~~~~~~~~~~~~~~ Useful Jacobians ~~~~~~~~~~~~~~~ #
"""
    ∇differential(q::UnitQuaternion)

Jacobian of `lmult(q) QuatMap(ϕ)`, when ϕ is near zero.

Useful for converting Jacobians from R⁴ to R³ and
    correctly account for unit norm constraint. Jacobians for different
    differential quaternion parameterization are the same up to a constant.
"""
function ∇differential(q::UnitQuaternion)
    SA[
        -q.x -q.y -q.z;
         q.w -q.z  q.y;
         q.z  q.w -q.x;
        -q.y  q.x  q.w;
    ]
end

"""
    ∇²differential(q::UnitQuaternion, b::AbstractVector)

Jacobian of `(∂/∂ϕ lmult(q) QuatMap(ϕ))`b, evaluated at ϕ=0, and `b` has length 4.
"""
function ∇²differential(q::UnitQuaternion, b::AbstractVector)
    check_length(b, 4)
    b1 = -params(q)'b
    Diagonal(@SVector fill(b1,3))
end

"""
    ∇rotate(R::Rotation{3}, r::AbstractVector)

Jacobian of `R*r` with respect to the rotation
"""
function ∇rotate(q::UnitQuaternion, r::AbstractVector)
    check_length(r, 3)
    rhat = UnitQuaternion(zero(eltype(r)), r[1], r[2], r[3], false)
    2vmat()*rmult(q)'rmult(rhat)
end

"""
    ∇composition1(R2::Rotation{3}, R1::Rotation{3})

Jacobian of `R2*R1` with respect to `R1`
"""
function ∇composition1(q2::UnitQuaternion, q1::UnitQuaternion)
    lmult(q2)
end

"""
    ∇composition2(R2::Rotation{3}, R1::Rotation{3})

Jacobian of `R2*R1` with respect to `R2`
"""
function ∇composition2(q2::UnitQuaternion, q1::UnitQuaternion)
    rmult(q1)
end


text:
module Frameworks

using Printf: @printf
using TimerOutputs: TimerOutputs, TimerOutput, @timeit
using ..Terms: Parameters
using ...Prerequisites: decimaltostr
using ...Prerequisites.CompositeStructures: NamedContainer
using ...Prerequisites.TypeTraits: efficientoperations

import ...Interfaces: id, update!, prepare!, register!, run!, add!

export App, Engine
export Assignment, Algorithm
export dependences, rundependences!

"""
    App

Abstract type for all apps.
"""
abstract type App end

"""
    ==(app1::App, app2::App) -> Bool

Judge whether two apps are equivalent to each other.
"""
Base.:(==)(app1::App, app2::App) = ==(efficientoperations, app1, app2)

"""
    isequal(app1::App, app2::App) -> Bool

Judge whether two apps are equivalent to each other.
"""
Base.isequal(app1::App, app2::App) = isequal(efficientoperations, app1, app2)

"""
    update!(app::App; kwargs...) -> App

Update the status of an app.
"""
update!(app::App; kwargs...) = app

"""
    Engine

Abstract type for all engines.
"""
abstract type Engine end

"""
    ==(engine1::Engine, engine2::Engine) -> Bool

Judge whether two engines are equivalent to each other.
"""
Base.:(==)(engine1::Engine, engine2::Engine) = ==(efficientoperations, engine1, engine2)

"""
    isequal(engine1::Engine, engine2::Engine) -> Bool

Judge whether two engines are equivalent to each other.
"""
Base.isequal(engine1::Engine, engine2::Engine) = isequal(efficientoperations, engine1, engine2)

"""
    update!(engine::Engine; kwargs...) -> Engine

Update the status of an engine.
"""
update!(engine::Engine; kwargs...) = engine

"""
    Assignment( id::Symbol, app::App, parameters::Parameters;
                map::Function=identity,
                dependences::Tuple{Vararg{Symbol}}=(),
                data::Any=nothing,
                savedata::Bool=true,
                virgin::Bool=true,
                kwargs...
                )

An assignment associated with an app.
"""
mutable struct Assignment{A<:App, P<:Parameters, M<:Function, D<:Tuple{Vararg{Symbol}}, R<:Any, I}
    app::A
    parameters::P
    map::M
    dependences::D
    data::R
    savedata::Bool
    virgin::Bool
end
function Assignment(id::Symbol, app::App, parameters::Parameters;
                    map::Function=identity,
                    dependences::Tuple{Vararg{Symbol}}=(),
                    data::Any=nothing,
                    savedata::Bool=true,
                    virgin::Bool=true,
                    kwargs...)
    A, P, M, D = app|>typeof, parameters|>typeof, map|>typeof, dependences|>typeof
    R = (data === nothing) ? Any : data|>typeof
    return Assignment{A, P, M, D, R, id}(app, parameters, map, dependences, data, savedata, virgin)
end

"""
    ==(assign1::Assignment, assign2::Assignment) -> Bool

Judge whether two assignments are equivalent to each other.
"""
Base.:(==)(assign1::Assignment, assign2::Assignment) = ==(efficientoperations, assign1, assign2)

"""
    isequal(assign1::Assignment, assign2::Assignment) -> Bool

Judge whether two assignments are equivalent to each other.
"""
Base.isequal(assign1::Assignment, assign2::Assignment) = isequal(efficientoperations, assign1, assign2)

"""
    valtype(assign::Assignment)
    valtype(::Type{<:Assignment{<:App, <:Parameters, <:Function, <:Tuple{Vararg{Symbol}}, R}}) where R

The type of the data(result) of an assignment.
"""
Base.valtype(assign::Assignment) = assign |> typeof |> valtype
Base.valtype(::Type{<:Assignment{<:App, <:Parameters, <:Function, <:Tuple{Vararg{Symbol}}, R}}) where {R} = R

"""
    id(assign::Assignment) -> Symbol
    id(::Type{<:Assignment{<:App, <:Parameters, <:Function, <:Tuple{Vararg{Symbol}}, <:Any, I}}) where I -> Symbol

The id of an assignment.
"""
id(assign::Assignment) = assign |> typeof |> id
id(::Type{<:Assignment{<:App, <:Parameters, <:Function, <:Tuple{Vararg{Symbol}}, <:Any, I}}) where {I} = I

"""
    update!(assign::Assignment; kwargs...) -> Assignment

Update the parameters of an assignment and the status of its associated app.
"""
@generated function update!(assign::Assignment; kwargs...)
    exprs = []
    names = fieldnames(fieldtype(assign, :parameters))
    for (i, name) in enumerate(names)
        name = QuoteNode(name)
        push!(exprs, :(get(kwargs, $name, getfield(assign.parameters, $i))))
    end
    return quote
        assign.parameters = Parameters{$names}($(exprs...))
        update!(assign.app; assign.map(assign.parameters)...)
        return assign
    end
end

"""
    Algorithm(  name::String, engine::Engine;
                din::String=".",
                dout::String=".",
                parameters::Union{Parameters, Nothing}=nothing,
                map::Function=identity,
                assignments::Tuple{Vararg{Assignment}}=(),
                kwargs...
                )

An algorithm associated with an engine.
"""
mutable struct Algorithm{E<:Engine, P<:Parameters, M<:Function, S<:NamedContainer{Assignment}}
    name::String
    engine::E
    din::String
    dout::String
    parameters::P
    map::M
    sassignments::S
    dassignments::Dict{Symbol, Assignment}
    timer::TimerOutput
end
function Algorithm( name::String, engine::Engine;
                    din::String=".",
                    dout::String=".",
                    parameters::Union{Parameters, Nothing}=nothing,
                    map::Function=identity,
                    assignments::Tuple{Vararg{Assignment}}=(),
                    kwargs...
                    )
    (parameters === nothing) && (parameters = Parameters(engine))
    assignments = namedassignments(assignments)
    return Algorithm(name, engine, din, dout, parameters, map, assignments, Dict{Symbol, Assignment}(), TimerOutput())
end
@generated function namedassignments(assignments::Tuple{Vararg{Assignment}})
    names, values = [], []
    for i = 1:fieldcount(assignments)
        push!(names, fieldtype(assignments, 1)|>id)
        push!(values, :(assignments[$i]))
    end
    names = NTuple{fieldcount(assignments), Symbol}(names)
    values = Expr(:tuple, values...)
    return :(NamedContainer{$names}($values))
end

"""
    update!(alg::Algorithm; kwargs...) -> Algorithm

Update the parameters of an algorithm and its associated engine.
"""
@generated function update!(alg::Algorithm; kwargs...)
    exprs = []
    names = fieldnames(fieldtype(alg, :parameters))
    for (i, name) in enumerate(names)
        name = QuoteNode(name)
        push!(exprs, :(get(kwargs, $name, getfield(alg.parameters, $i))))
    end
    return quote
        alg.parameters = Parameters{$names}($(exprs...))
        update!(alg.engine; alg.map(alg.parameters)...)
        return alg
    end
end

"""
    repr(alg::Algorithm, mask::Tuple{Vararg{Symbol}}=(); ndecimal::Int=10) -> String

Get the repr representation of an algorithm.

Optionally, some parameters of the algorithm can be masked. Besides, the maximum number of decimals of the parameters can also be specified.
"""
function Base.repr(alg::Algorithm, mask::Tuple{Vararg{Symbol}}=(); ndecimal::Int=10)
    result = [string(alg.name), repr(alg.engine)]
    for (name, value) in pairs(alg.parameters)
        (name ∉ mask) && push!(result, decimaltostr(value, ndecimal))
    end
    return join(result, "_")
end

"""
    show(io::IO, alg::Algorithm)

Show an algorithm.
"""
function Base.show(io::IO, alg::Algorithm)
    @printf io "%s_%s" alg.name alg.engine
    for (name, value) in pairs(alg.parameters)
        @printf io "_%s" decimaltostr(value, 10)
    end
end

"""
    get(alg::Algorithm, id::Symbol) -> Assignment
    get(alg::Algorithm, ::Val{id}) where id -> Assignment

Find the assignment registered on a algorithm by its id.
"""
Base.get(alg::Algorithm, id::Symbol) = (id ∈ keys(alg.sassignments)) ? getfield(alg.sassignments, id) : alg.dassignments[id]
@generated Base.get(alg::Algorithm, ::Val{id}) where {id} = (id ∈ fieldnames(fieldtype(alg, :sassignments))) ? :(getfield(alg.sassignments, id)) : :(alg.dassignments[id])

"""
    summary(alg::Algorithm)

Provide a summary of an algorithm.
"""
function Base.summary(alg::Algorithm)
    @info "Summary of $(alg.name)($(nameof(typeof(alg.engine)))):"
    @info string(alg.timer)
end

"""
    prepare!(alg::Algorithm, assign::Assignment) -> Nothing

Prepare an assignment registered on a algorithm.
"""
prepare!(alg::Algorithm, assign::Assignment) = nothing

"""
    run!(alg::Algorithm, assign::Assignment) -> Nothing

Run an assignment registered on a algorithm.
"""
run!(alg::Algorithm, assign::Assignment) = nothing

"""
    register!(alg::Algorithm, id::Symbol, app::App; kwargs...) -> Algorithm

Add an assignment on a algorithm by providing the contents of the assignment, and run this assignment.
"""
function register!(alg::Algorithm, id::Symbol, app::App; kwargs...)
    add!(alg, id, app; kwargs...)
    run!(alg, id, true)
end

"""
    add!(alg::Algorithm, id::Symbol, app::App; kwargs...) -> Algorithm

Add an assignment on a algorithm by providing the contents of the assignment.

The difference between `add!` and `register!` is that the `add!` function does not run the newly added assignment but the `register!` function does.
"""
function add!(alg::Algorithm, id::Symbol, app::App; kwargs...)
    @assert id ∉ keys(alg.sassignments) "add! error: id($id) conflict."
    alg.dassignments[id] = Assignment(id, app, merge(alg.parameters, get(kwargs, :parameters, Parameters{()}())); kwargs...)
    return alg
end

"""
    run!(alg::Algorithm, id::Symbol, timing::Bool=true) -> Algorithm
    run!(alg::Algorithm, ::Val{id}, timing::Bool=true) where id -> Algorithm

Run an assignment with the given id registered on an algorithm. Optionally, the run process can be timed by setting the `timing` argument to be `true`.
"""
run!(alg::Algorithm, id::Symbol, timing::Bool = true) = run!(alg, Val(id), timing)
function run!(alg::Algorithm, ::Val{id}, timing::Bool = true) where id
    assign = get(alg, Val(id))
    if timing
        @timeit alg.timer string(id) algrunassignment!(alg, assign)
        @info "App $id($(nameof(assign.app|>typeof))): time consumed $(TimerOutputs.time(alg.timer[string(id)]) / 10^9)s."
    else
        algrunassignment!(alg, assign)
    end
    return alg
end
function algrunassignment!(alg::Algorithm, assign::Assignment)
    ismatched = match(assign.parameters, alg.parameters)
    !assign.virgin && ismatched && return
    !ismatched && update!(alg; assign.parameters...)
    prepare!(alg, assign)
    run!(alg, assign)
    assign.virgin = false
end

"""
    dependences(alg::Algorithm, assign::Assignment, ::Tuple{}=()) -> Tuple{Vararg{Symbol}}
    dependences(alg::Algorithm, assign::Assignment, mask::Tuple{Vararg{Symbol}}) -> Tuple{Vararg{Symbol}}

Get the dependences of an assignment and return their ids.
"""
dependences(alg::Algorithm, assign::Assignment, ::Tuple{}=()) = assign.dependences
dependences(alg::Algorithm, assign::Assignment, mask::Tuple{Vararg{Symbol}}) = Tuple(filter(x->(x ∉ mask), collect(assign.dependences)))

"""
    rundependences!(alg::Algorithm, assign::Assignment, mask::Tuple{Vararg{Symbol}}=()) -> Algorithm

Run the dependences of an assignment. Optionally, some dependences can be jumped by specifying the `mask` argument.
"""
function rundependences!(alg::Algorithm, assign::Assignment, mask::Tuple{Vararg{Symbol}}=())
    for id in dependences(alg, assign, mask)
        assign = get(alg, id)
        algrundependence!(alg, assign)
    end
    return alg
end
function algrundependence!(alg::Algorithm, assign::Assignment)
    ismatched = match(alg.parameters, assign.parameters)
    !assign.virgin && ismatched && return
    !ismatched && update!(assign; alg.parameters...)
    prepare!(alg, assign)
    run!(alg, assign)
    assign.virgin = false
end

end  #module

ext: jl


text:
include("toolbox.jl")
include("lbfgsb-wrapper.jl")

function cvxreg_qp_ineq_APG(Xflat::Array{T,1}, Y::Array{T,1}, rho::T, I::Array{Int64,1}, J::Array{Int64,1}, Wlen::Int64,
        lamb::LambdaActiveSet, n::Int64,d::Int64, phi_tmp::Array{T,1}, xi_tmp::Array{T,1}; max_steps=1000, TOL=1e-4) where T<:AbstractFloat
    θ::Float64 = 1.0;
    θ_new::Float64 = 1.0;
    L = power_iter(Xflat,I,J,Wlen,rho,n,d);
    λ_x = copy(lamb.nzval);
    λ_y = copy(lamb.nzval);
    λ_new = zeros(eltype(Y),Wlen);
    s_tmp = zeros(eltype(Y),Wlen);
    grad = zeros(eltype(Y),Wlen);
    fnew = 0.;
    f = get_objective_value(Xflat,Y,rho,I,J, Wlen, λ_x,n, d, phi_tmp, xi_tmp, s_tmp);
    k = 0;
    for k in 1:max_steps
        get_gradient_evaluation!(grad, Xflat,Y,rho,I,J, Wlen, λ_y,n,d, phi_tmp, xi_tmp);
        λ_new .= min.(λ_y - grad/L,0);
        fnew = get_objective_value(Xflat,Y,rho,I,J, Wlen, λ_new,n, d, phi_tmp, xi_tmp, s_tmp);
        if (fnew > f)
                θ = 1;
                λ_new .= λ_x;
                λ_y .= λ_x;
        end
        if (fnew <= f && fnew - f >-TOL && fnew>f*(1+TOL))
                break
        end
        θ_new = (1+sqrt(1+4*θ^2))/2;
        γ = (θ-1)/θ_new;
        λ_y .= λ_new .+ γ*(λ_new .- λ_x);
        λ_x .= λ_new;
        θ = θ_new;
        f = fnew;
    end
    lamb.nzval .= λ_new
    phi_tmp .+= Y
    phi_tmp .*= -1
    xi_tmp .*= -1

    return fnew
end

function cvxreg_qp_ineq_APG_profiling(Xflat::Array{T,1}, Y::Array{T,1}, rho::T, I::Array{Int64,1}, J::Array{Int64,1}, Wlen::Int64,
        lamb::LambdaActiveSet, n::Int64,d::Int64, phi_tmp::Array{T,1}, xi_tmp::Array{T,1}; max_steps=1000, TOL=1e-4, start_time = 0., end_time = Inf) where T<:AbstractFloat
    all_times = zeros(Float64,max_steps+1);
    cur_time = @elapsed begin
        θ::Float64 = 1.0;
        θ_new::Float64 = 1.0;
        pow_time = @elapsed L = power_iter(Xflat,I,J,Wlen,rho,n,d);
        λ_x = copy(lamb.nzval);
        λ_y = copy(lamb.nzval);
        λ_new = zeros(eltype(Y),Wlen);
        s_tmp = zeros(eltype(Y),Wlen);
        grad = zeros(eltype(Y),Wlen);
        fnew = 0.;
        objs = zeros(Float64,max_steps+1);
        f = get_objective_value(Xflat,Y,rho,I,J, Wlen, λ_x,n, d, phi_tmp, xi_tmp, s_tmp);
        k = 0;
        iter = 0;
        objs[1] = f;
        status= 0;
    end
    all_times[1] = cur_time;
    for k in 1:max_steps
        cur_time += @elapsed begin
            get_gradient_evaluation!(grad, Xflat,Y,rho,I,J, Wlen, λ_y,n,d, phi_tmp, xi_tmp);
            λ_new .= min.(λ_y - grad/L,0);
            fnew = get_objective_value(Xflat,Y,rho,I,J, Wlen, λ_new,n, d, phi_tmp, xi_tmp, s_tmp);
        end
        objs[k+1] = fnew;
        all_times[k+1] = cur_time;
        iter += 1;
        cur_time += @elapsed begin
            if (fnew > f)
                    θ = 1;
                    λ_new .= λ_x;
                    λ_y .= λ_x;
            end
            if (fnew <= f && fnew - f >-TOL && fnew>f*(1+TOL))
                status = 1
                    break
            end
            θ_new = (1+sqrt(1+4*θ^2))/2;
            γ = (θ-1)/θ_new;
            λ_y .= λ_new .+ γ*(λ_new .- λ_x);
            λ_x .= λ_new;
            θ = θ_new;
            f = fnew;
        end
        if start_time + cur_time >= end_time
            println("early stopping during optimization due to time limit")
            status = -1
            break
        end
    end
    lamb.nzval .= λ_new
    phi_tmp .+= Y
    phi_tmp .*= -1
    xi_tmp .*= -1
    return status, fnew, all_times[1:iter+1], objs[1:iter+1], pow_time, L, iter
end

function cvxreg_qp_ineq_APG_ls(Xflat::Array{T,1}, Y::Array{T,1}, rho::T, I::Array{Int64,1}, J::Array{Int64,1}, Wlen::Int64,
        lamb::LambdaActiveSet, L::T, n::Int64, d::Int64, phi_tmp::Array{T,1}, xi_tmp::Array{T,1};
        β_1 = 0.9, β_2 = 0.98, max_steps=1000::Int64, TOL=1e-4::T, maxls = 20::Int64) where T<:AbstractFloat
    θ::Float64 = 1.0;
    θ_new::Float64 = 1.0;
    cur_ls = 0;
    if L <= 0
        L = power_iter(Xflat,I,J,Wlen,rho,n,d);
        cur_ls = maxls;
    end
    λ_x = copy(lamb.nzval);
    λ_y = copy(lamb.nzval);
    λ_new = zeros(eltype(Y),Wlen);
    s_tmp = zeros(eltype(Y),Wlen);
    grad = zeros(eltype(Y),Wlen);
    fnew = 0.;
    f = get_objective_value(Xflat,Y,rho,I,J, Wlen, λ_x,n, d, phi_tmp, xi_tmp, s_tmp);
    k = 0;
    for k in 1:max_steps
        get_gradient_evaluation!(grad, Xflat,Y,rho,I,J, Wlen, λ_y,n,d, phi_tmp, xi_tmp);
        while (true)
            λ_new .= min.(λ_y - grad/L,0);
            fnew = get_objective_value(Xflat,Y,rho,I,J, Wlen, λ_new,n, d, phi_tmp, xi_tmp, s_tmp);
            if cur_ls >= maxls || fnew <= f + dot(grad, λ_new - λ_y) + L/2 * norm(λ_new - λ_y,2)^2
                cur_ls = 0
                break
            else
                L /= ((k==1) ? β_1 : β_2);
                cur_ls += 1
                if cur_ls >= maxls
                    L = power_iter(Xflat,I,J,Wlen,rho,n,d,maxiter=50);
                    cur_ls = 0;
                    break
                end
            end
        end
        if (fnew > f)
            θ = 1;
            λ_new .= λ_x;
            λ_y .= λ_x;
        end
        if (fnew <= f && fnew - f >-TOL && fnew>f*(1+TOL))
            break
        end
        θ_new = (1+sqrt(1+4*θ^2))/2;
        γ = (θ-1)/θ_new;
        λ_y .= λ_new .+ γ*(λ_new .- λ_x);
        λ_x .= λ_new;
        θ = θ_new;
        f = fnew;
    end
    lamb.nzval .= λ_new
    phi_tmp .+= Y
    phi_tmp .*= -1
    xi_tmp .*= -1

    return L, fnew
end

function cvxreg_qp_ineq_APG_ls_profiling(Xflat::Array{T,1}, Y::Array{T,1}, rho::T, I::Array{Int64,1}, J::Array{Int64,1}, Wlen::Int64,
        lamb::LambdaActiveSet, L::T, n::Int64, d::Int64, phi_tmp::Array{T,1}, xi_tmp::Array{T,1};
        β_1 = 0.9, β_2 = 0.98, max_steps=1000::Int64, TOL=1e-4::T, maxls = 20::Int64, start_time = 0., end_time = Inf) where T<:AbstractFloat
    all_times = zeros(Float64,max_steps+1);
    cur_time = @elapsed begin
        θ::Float64 = 1.0;
        θ_new::Float64 = 1.0;
        cur_ls = 0;
        if L <= 0
            L = power_iter(Xflat,I,J,Wlen,rho,n,d);
            cur_ls = maxls;
        end
        λ_x = copy(lamb.nzval);
        λ_y = copy(lamb.nzval);
        λ_new = zeros(eltype(Y),Wlen);
        s_tmp = zeros(eltype(Y),Wlen);
        grad = zeros(eltype(Y),Wlen);
        fnew = 0.;
        objs = zeros(Float64,max_steps+1);
        f_evals = zeros(Float64,max_steps);
        Ls = zeros(Float64,max_steps);
        f = get_objective_value(Xflat,Y,rho,I,J, Wlen, λ_x,n, d, phi_tmp, xi_tmp, s_tmp);
        k = 0;
        iter = 0;
        objs[1] = f;
        status = 0;
    end
    all_times[1] = cur_time;
    for k in 1:max_steps
        cur_time += @elapsed begin
            get_gradient_evaluation!(grad, Xflat,Y,rho,I,J, Wlen, λ_y,n,d, phi_tmp, xi_tmp);
            while (true)
                λ_new .= min.(λ_y - grad/L,0);
                fnew = get_objective_value(Xflat,Y,rho,I,J, Wlen, λ_new,n, d, phi_tmp, xi_tmp, s_tmp);
                if cur_ls >= maxls || fnew <= f + dot(grad, λ_new - λ_y) + L/2 * norm(λ_new - λ_y,2)^2
                    f_evals[k] = cur_ls+1;
                    cur_ls = 0
                    break
                else
                    L /= ((k==1) ? β_1 : β_2);
                    cur_ls += 1
                    if cur_ls >= maxls
                        L = power_iter(Xflat,I,J,Wlen,rho,n,d,maxiter=50);
                        f_evals[k] = cur_ls+1;
                        cur_ls = 0;
                        break
                    end
                end
            end
        end
        Ls[k] = L;
        objs[k+1] = fnew;
        all_times[k+1] = cur_time;
        iter += 1;
        cur_time += @elapsed begin
            if (fnew > f)
                θ = 1;
                λ_new .= λ_x;
                λ_y .= λ_x;
            end
            if (fnew <= f && fnew - f >-TOL && fnew>f*(1+TOL))
                status = 1;
                break
            end
            θ_new = (1+sqrt(1+4*θ^2))/2;
            γ = (θ-1)/θ_new;
            λ_y .= λ_new .+ γ*(λ_new .- λ_x);
            λ_x .= λ_new;
            θ = θ_new;
            f = fnew;
        end
        if start_time + cur_time >= end_time
            println("early stopping during optimization due to time limit")
            status = -1
            break
        end
    end
    lamb.nzval .= λ_new
    phi_tmp .+= Y
    phi_tmp .*= -1
    xi_tmp .*= -1

    return status, L, fnew, all_times[1:iter+1], objs[1:iter+1], f_evals[1:iter], Ls[1:iter], iter
end

function cvxreg_qp_ineq_LBFGS(Xflat::Array{T,1}, Y::Array{T,1}, rho::T, I::Array{Int64,1}, J::Array{Int64,1}, Wlen::Int64,
        lamb::LambdaActiveSet, n::Int64, d::Int64, phi_tmp::Array{T,1}, xi_tmp::Array{T,1};
        max_steps=1000::Int64, TOL=1e-4::T) where T<:AbstractFloat
    s_tmp = zeros(eltype(Y),Wlen);
    f(G,λ) = get_objective_value_gradient!(G,Xflat,Y,rho,I,J, Wlen, λ,n, d, phi_tmp, xi_tmp, s_tmp)
    obj, minimizer = lbfgsb(f,lamb.nzval; ub=zeros(Wlen), ftol=TOL, maxiter=max_steps, profiling=false);
    lamb.nzval .= minimizer
    phi_tmp .+= Y
    phi_tmp .*= -1
    xi_tmp .*= -1
    return obj
end

function cvxreg_qp_ineq_LBFGS_profiling(Xflat::Array{T,1}, Y::Array{T,1}, rho::T, I::Array{Int64,1}, J::Array{Int64,1}, Wlen::Int64,
        lamb::LambdaActiveSet, n::Int64, d::Int64, phi_tmp::Array{T,1}, xi_tmp::Array{T,1};
        max_steps=1000::Int64, TOL=1e-4::T) where T<:AbstractFloat
    all_times = zeros(Float64, 1+max_steps)
    cur_time = @elapsed begin
        s_tmp = zeros(eltype(Y),Wlen);
        f(G,λ) = get_objective_value_gradient!(G,Xflat,Y,rho,I,J, Wlen, λ,n, d, phi_tmp, xi_tmp, s_tmp)
    end
    all_times[1] = cur_time
    total_time = @elapsed obj,minimizer,objs,f_evals,iter = lbfgsb(f,lamb.nzval; ub=zeros(Wlen), ftol=TOL, maxiter=max_steps, profiling = true);
    all_times[2:iter+1] = cumsum(f_evals)/sum(f_evals)*total_time
    lamb.nzval .= minimizer
    phi_tmp .+= Y
    phi_tmp .*= -1
    xi_tmp .*= -1
    return obj, all_times[1:iter+1], objs, f_evals, iter
end


text:
module Experimental

using NBodyIPs: NBodyFunction, BondLengthDesc, ClusterBLDesc,
                descriptor, bodyorder, NBodyIP, transform, cutoff,
                fast
using NBodyIPs.Polys: NBPoly, StNBPoly
using NBodyIPs.EnvIPs: EnvIP, EnvPoly
using JuLIP: AbstractCalculator

function faster_blpot(IP)
   components = AbstractCalculator[]
   for V in IP.components
      modify = false
      if (V isa NBPoly)
         if (bodyorder(V) >= 2) && (descriptor(V) isa BondLengthDesc)
            modify=true
         end
      end
      if modify
         desc_V = descriptor(V)
         desc = ClusterBLDesc(transform(desc_V), desc_V.cutoff)
         Vcl = NBPoly(V.t, V.c * bodyorder(V), desc)
         Vff = fast(Vcl)
      elseif V isa NBPoly # NBPoly but not a 3B or higher BL
         Vff = fast(V)
      else
         Vff = V
      end
      push!(components, Vff)
   end
   return NBodyIP(components)
end

function faster_envpot(IP)
   @warn("""Warning - this is an experimental code; I'm not checking
            compatibility of the individual components!""")
   components = AbstractCalculator[]
   envs = EnvIP[]
   for V in IP.components
      if !(V isa EnvIP)
         push!(components, V)
      else
         push!(envs, V)
      end
   end
   envdegs = [V.t for V in envs]           # env degrees
   @assert maximum(envdegs) <= 4    # TODO: why?!?!?
   bos = [bodyorder(V.Vr) for V in envs]   #
   for bo = 2:4  # loop over body-orders
      ibo = findall(bos .== bo)  # find the relevant subset
      isempty(ibo) && continue
      # check that _all_ envdegs are present
      @assert sort(envdegs[ibo]) == 0:(length(ibo)-1)
      Vrs = Vector{Any}(undef, length(ibo))
      for i in ibo
         Vrs[envs[i].t+1] = envs[i].Vr
      end
      # check that all Vrs have the same type
      # and also make them fast!
      Vrs = [fast(V) for V in Vrs]
      @assert isconcretetype(eltype(Vrs))
      # construct the combined EnvPoly
      Vn = envs[ibo[1]].Vn
      str_Vn = envs[ibo[1]].str_Vn
      push!(components, EnvPoly(Vrs, Vn, str_Vn))
   end
   return NBodyIP(components)
end

end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


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

name = "LibVPX"
version = v"1.9.0"

# Collection of sources required to build LibVPX
sources = [
    ArchiveSource("https://github.com/webmproject/libvpx/archive/v$(version).tar.gz",
                  "d279c10e4b9316bf11a570ba16c3d55791e1ad6faa4404c67422eb631782c80a"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/libvpx-*/
sed -i 's/cp -p/cp/' build/make/Makefile

mkdir vpx_build && cd vpx_build
apk add diffutils yasm

if [[ "${target}" == i686-linux-* ]]; then
    export TARGET=x86-linux-gcc
elif [[ "${target}" == x86_64-linux-* ]]; then
    export TARGET=x86_64-linux-gcc
elif [[ "${target}" == arm-linux-* ]]; then
    export TARGET=armv7-linux-gcc
elif [[ "${target}" == aarch64-linux-* ]]; then
    export TARGET=arm64-linux-gcc
elif [[ "${target}" == powerpc64le-linux-* ]]; then
    export TARGET=ppc64le-linux-gcc
elif [[ "${target}" == x86_64-apple-* ]]; then
    export TARGET=x86_64-darwin14-gcc
elif [[ "${target}" == i686-w64-mingw32 ]]; then
    export TARGET=x86-win32-gcc
elif [[ "${target}" == x86_64-w64-mingw32 ]]; then
    export TARGET=x86_64-win64-gcc
    export CFLAGS="${CFLAGS} -fno-asynchronous-unwind-tables"
elif [[ "${target}" == *freebsd* ]]; then
    export TARGET=generic-gnu
    export CONFIG_OPTS="--disable-multithread"
fi

../configure --prefix=$prefix --target=${TARGET} \
    --as=yasm \
    --enable-postproc \
    --enable-pic \
    --enable-vp8 \
    --enable-vp9 \
    --enable-vp9-highbitdepth \
    --enable-runtime-cpu-detect \
    ${CONFIG_OPTS}
echo "SRC_PATH_BARE=.." >> config.mk
echo "target=libs" >> config.mk
make -j${nproc}
make install

# pkgconfig file on Windows is installed to ${libdir}/pkgconfig,
# we have to move it to ${prefix}/lib/pkgconfig/
if [[ "${target}" == *-mingw* ]] && [[ -d "${libdir}/pkgconfig" ]] ; then
    mkdir -p "${prefix}/lib"
    mv "${libdir}/pkgconfig" "${prefix}/lib/pkgconfig"
fi
"""

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

# The products that we will ensure are always built
products = [
    # While we want this, it's not available on windows.  :/
    #LibraryProduct("libvpx", :libvpx),
    FileProduct("\${libdir}/libvpx.a", :libvpx),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
using NeXLMatrixCorrection
using Test

@testset "XPP" begin
    @testset "K240" begin # Test XPP against DTSA-II using K240

        toa = deg2rad(40.0)
        k240 = NeXLCore.material(
            "K240",
            n"O" => 0.340023,
            n"Mg" => 0.030154,
            n"Si" => 0.186986,
            n"Ti" => 0.059950,
            n"Zn" => 0.040168,
            n"Zr" => 0.074030,
            n"Ba" => 0.268689,
        )

        cxr = n"O K-L3"
        xpp = NeXLMatrixCorrection.XPP(k240, inner(cxr), 20.0e3)
        @test isapprox(xpp.A, 393.01, atol = 0.01)
        @test isapprox(xpp.a, 6754, atol = 1.0)
        @test isapprox(xpp.B, -378725, atol = 1.0)
        @test isapprox(xpp.b, 7786, atol = 1.0)
        @test isapprox(xpp.ϕ0, 1.63011, atol = 0.0001)
        @test isapprox(xpp.F, 0.0016754, atol = 0.000001)

        xppO = NeXLMatrixCorrection.XPP(pure(element(cxr)), inner(cxr), 20.0e3)
        ZA = Fχ(xpp, cxr, toa) / Fχ(xppO, cxr, toa)
        Z = F(xpp) / F(xppO)
        @test isapprox(Z, 1.1759, atol = 0.0001)
        @test isapprox(ZA / Z, 0.2951, atol = 0.0001)

        cxr = n"Zr L3-M5"
        xpp = NeXLMatrixCorrection.XPP(k240, inner(cxr), 20.0e3)
        @test isapprox(xpp.A, 129.53, atol = 0.01)
        @test isapprox(xpp.a, 6668, atol = 1.0)
        @test isapprox(xpp.B, -177201, atol = 1.0)
        @test isapprox(xpp.b, 8239, atol = 1.0)
        @test isapprox(xpp.ϕ0, 1.60571, atol = 0.0001)
        @test isapprox(xpp.F, 0.0012887, atol = 0.00001)

        xppO = NeXLMatrixCorrection.XPP(pure(element(cxr)), inner(cxr), 20.0e3)
        ZA = Fχ(xpp, cxr, toa) / Fχ(xppO, cxr, toa)
        Z = F(xpp) / F(xppO)
        @test isapprox(Z, 0.8502, atol = 0.0001)
        @test isapprox(ZA / Z, 0.7463, atol = 0.0001)

        cxr = n"Mg K-L3"
        xpp = NeXLMatrixCorrection.XPP(k240, inner(cxr), 20.0e3)
        @test isapprox(xpp.A, 215.96, atol = 0.01)
        @test isapprox(xpp.a, 6681, atol = 1.0)
        @test isapprox(xpp.B, -252614, atol = 1.0)
        @test isapprox(xpp.b, 7973, atol = 1.0)
        @test isapprox(xpp.ϕ0, 1.620603, atol = 0.0001)
        @test isapprox(xpp.F, 0.001468275, atol = 0.00001)

        xppO = NeXLMatrixCorrection.XPP(pure(element(cxr)), inner(cxr), 20.0e3)
        ZA = Fχ(xpp, cxr, toa) / Fχ(xppO, cxr, toa)
        Z = F(xpp) / F(xppO)
        @test isapprox(Z, 1.0890, atol = 0.0001)
        @test isapprox(ZA / Z, 0.3889, atol = 0.0001)
    end

    @testset "K240 ZAF" begin

        k240 = NeXLCore.material(
            "K240",
            n"O" => 0.340023,
            n"Mg" => 0.030154,
            n"Si" => 0.186986,
            n"Ti" => 0.059950,
            n"Zn" => 0.040168,
            n"Zr" => 0.074030,
            n"Ba" => 0.268689,
        )

        sio2 = atomicfraction("Quartz", n"Si" => 1, n"O" => 2)
        mgo = atomicfraction("MgO", n"Mg" => 1, n"O" => 1)
        baf2 = atomicfraction("Barium Fluoride", n"Ba" => 1, n"F" => 2)
        ti, zn, zr = pure(n"Ti"), pure(n"Zn"), pure(n"Zr")

        e0, θ = 17.0e3, deg2rad(40.0)

        zafSi = zafcorrection(XPP, ReedFluorescence, Coating, k240, sio2, n"Si K", e0)
        zafMg = zafcorrection(XPP, ReedFluorescence, Coating, k240, mgo, n"Mg K", e0)
        zafBa = zafcorrection(XPP, ReedFluorescence, Coating, k240, baf2, n"Ba L3", e0)
        zafTi = zafcorrection(XPP, ReedFluorescence, Coating, k240, ti, n"Ti K", e0)
        zafZn = zafcorrection(XPP, ReedFluorescence, Coating, k240, zn, n"Zn K", e0)
        zafZr = zafcorrection(XPP, ReedFluorescence, Coating, k240, zr, n"Zr L3", e0)
        zafO = zafcorrection(XPP, ReedFluorescence, Coating, k240, sio2, n"O K", e0)

        @test isapprox(Z(zafSi...), 1.1345, atol = 0.001) # 1.1411
        @test isapprox(Z(zafMg...), 1.1280, atol = 0.001) # 1.1279
        @test isapprox(Z(zafBa...), 0.8259, atol = 0.001) # 0.8261
        @test isapprox(Z(zafTi...), 0.9446, atol = 0.001) # 0.9439
        @test isapprox(Z(zafZn...), 0.8976, atol = 0.001) # 0.8978
        @test isapprox(Z(zafZr...), 0.8442, atol = 0.001) # 0.8443
        @test isapprox(Z(zafO...), 1.1316, atol = 0.001) # 1.1382

        @test isapprox(A(zafSi..., n"Si K-L3", θ, θ), 0.7652, atol = 0.001) # 0.7407
        @test isapprox(A(zafMg..., n"Mg K-L3", θ, θ), 0.5996, atol = 0.001) # 0.6021
        @test isapprox(A(zafBa..., n"Ba L3-M5", θ, θ), 1.0124, atol = 0.001) # 1.0128
        @test isapprox(A(zafTi..., n"Ti K-L3", θ, θ), 0.9607, atol = 0.001) # 0.9567
        @test isapprox(A(zafZn..., n"Zn K-L3", θ, θ), 0.9842, atol = 0.001) # 0.9837
        @test isapprox(A(zafZr..., n"Zr L3-M5", θ, θ), 0.7929, atol = 0.001) # 0.6826
        @test isapprox(A(zafO..., n"O K-L3", θ, θ), 0.7750, atol = 0.001) # 0.8310

        @test isapprox(ZA(zafSi..., n"Si K-L3", θ, θ), 1.1345 * 0.7652, atol = 0.001)
        @test isapprox(ZA(zafMg..., n"Mg K-L3", θ, θ), 1.1280 * 0.5996, atol = 0.001)
        @test isapprox(ZA(zafBa..., n"Ba L3-M5", θ, θ), 0.8259 * 1.0124, atol = 0.001)
        @test isapprox(ZA(zafTi..., n"Ti K-L3", θ, θ), 0.9446 * 0.9607, atol = 0.001)
        @test isapprox(ZA(zafZn..., n"Zn K-L3", θ, θ), 0.8973 * 0.9842, atol = 0.001)
        @test isapprox(ZA(zafZr..., n"Zr L3-M5", θ, θ), 0.8442 * 0.7929, atol = 0.001)
        @test isapprox(ZA(zafO..., n"O K-L3", θ, θ), 1.1316 * 0.7750, atol = 0.001)

        @test isapprox(F(zafSi..., n"Si K-L3", θ, θ), 1.0015, atol = 0.001) # 1.0029
        @test isapprox(F(zafMg..., n"Mg K-L3", θ, θ), 1.0041, atol = 0.001) # 1.0038
        @test isapprox(F(zafBa..., n"Ba L3-M5", θ, θ), 0.9998, atol = 0.001) # 0.9997
        @test isapprox(F(zafTi..., n"Ti K-L3", θ, θ), 1.0042, atol = 0.001) # 1.0007
        @test isapprox(F(zafZn..., n"Zn K-L3", θ, θ), 1.000, atol = 0.001) # 1.0000
        @test isapprox(F(zafZr..., n"Zr L3-M5", θ, θ), 1.0020, atol = 0.001) # 1.0106
        @test isapprox(F(zafO..., n"O K-L3", θ, θ), 0.9996, atol = 0.001) # 0.9997
        # TryZAF results
        # K240	C	    Z	    A	    F	    k		        C`	    Z`	    A`	    F`	    k`
        # O	    0.3400	1.1829	0.3373	1.0003	0.1357	SiO2	0.6348	1.1382	0.8310	0.9997	0.6002
        # Mg	0.0302	1.0927	0.4579	1.0038	0.0152	MgO	    0.0501	1.1279	0.6021	1.0038	0.0341
        # Si	0.1870	1.0806	0.6562	1.0029	0.1330	SiO2	0.3984	1.1411	0.7407	1.0029	0.3377
        # Ti	0.0600	0.9439	0.9567	1.0007	0.0542	Ti
        # Zn	0.4020	0.8978	0.9837	1.0000	0.3550	Zn
        # Zr	0.0740	0.8443	0.6826	1.0106	0.0431	Zr
        # Ba	0.2687	0.7457	1.0361	1.0038	0.2084	BaF2	0.3430	0.8261	1.0128	1.0038	0.2881
        # SiO2
        # Si	0.4694	0.9470	0.8859	1.0000	0.3938
        # O	    0.5356	1.0393	0.4059	1.0006	0.2261
        # MgO
        # Mg	0.6030	0.9688	0.7605	1.0000	0.4443
        # O	    0.3970	1.0500	0.5357	1.0018	0.2237
        # BaF2
        # Ba	0.7833	0.9027	1.0230	1.0000	0.7233
        # F	    0.2167	1.3080	0.6009	1.0001	0.1703
        #    print(asa(DataFrame, Dict( [ zafSi, zafMg, zafBa, zafTi, zafZn, zafZr, zafO ]), θ))
    end

    @testset "U3O8 at 25 keV" begin

        e0, toa = 25.0e3, deg2rad(40.0)
        u3o8 = atomicfraction("U3O8", n"U" => 3, n"O" => 8)
        k227 = material("K227", n"O" => 0.1639, n"Si" => 0.0935, n"Pb" => 0.7427)
        u = pure(n"U")

        zafO = zafcorrection(XPP, ReedFluorescence, Coating, u3o8, k227, n"O K", e0)
        zafU_L = zafcorrection(XPP, ReedFluorescence, Coating, u3o8, u, n"U L3", e0)
        zafU_M = zafcorrection(XPP, ReedFluorescence, Coating, u3o8, u, n"U M5", e0)

        @test isapprox(Z(zafO...), 1.0735, atol = 0.001)
        @test isapprox(Z(zafU_L...), 0.8965, atol = 0.001)
        @test isapprox(Z(zafU_M...), 0.9197, atol = 0.001)

        @test isapprox(A(zafO..., n"O K-L3", toa, toa), 1.5772, atol = 0.001)
        @test isapprox(A(zafU_L..., n"U L3-M5", toa, toa), 1.0089, atol = 0.001)
        @test isapprox(A(zafU_M..., n"U M5-N7", toa, toa), 1.0533, atol = 0.001)

        @test isapprox(F(zafO..., n"O K-L3", toa, toa), 0.9999, atol = 0.001)
        @test isapprox(F(zafU_L..., n"U L3-M5", toa, toa), 0.9999, atol = 0.001)
        @test isapprox(F(zafU_M..., n"U M5-N7", toa, toa), 0.9999, atol = 0.001)

    end

    @testset "K240 multiZAF" begin

        k240 = NeXLCore.material(
            "K240",
            n"O" => 0.340023,
            n"Mg" => 0.030154,
            n"Si" => 0.186986,
            n"Ti" => 0.059950,
            n"Zn" => 0.040168,
            n"Zr" => 0.074030,
            n"Ba" => 0.268689,
        )

        sio2 = atomicfraction("Quartz", n"Si" => 1, n"O" => 2)
        mgo = atomicfraction("MgO", n"Mg" => 1, n"O" => 1)
        baf2 = atomicfraction("Barium Fluoride", n"Ba" => 1, n"F" => 2)
        ti, zn, zr = pure(n"Ti"), pure(n"Zn"), pure(n"Zr")

        e0, θ = 17.0e3, deg2rad(40.0)

        cxrSi = characteristic(n"Si", ktransitions, 0.001, e0)
        zafSi = zafcorrection(XPP, ReedFluorescence, Coating, k240, sio2, cxrSi, e0)
        cxrMg = characteristic(n"Mg", ktransitions, 0.001, e0)
        zafMg = zafcorrection(XPP, ReedFluorescence, Coating, k240, mgo, cxrMg, e0)
        cxrBa = characteristic(n"Ba", ltransitions, 0.001, e0)
        zafBa = zafcorrection(XPP, ReedFluorescence, Coating, k240, baf2, cxrBa, e0)
        cxrTi = characteristic(n"Ti", ktransitions, 0.001, e0)
        zafTi = zafcorrection(XPP, ReedFluorescence, Coating, k240, ti, cxrTi, e0)
        cxrZn = characteristic(n"Zn", kalpha, 0.001, e0)
        zafZn = zafcorrection(XPP, ReedFluorescence, Coating, k240, zn, cxrZn, e0)
        cxrZr = characteristic(n"Zr", ltransitions, 0.001, e0)
        zafZr = zafcorrection(XPP, ReedFluorescence, Coating, k240, zr, cxrZr, e0)
        cxrO = characteristic(n"O", ktransitions, 0.001, e0)
        zafO = zafcorrection(XPP, ReedFluorescence, Coating, k240, sio2, cxrO, e0)

        @test isapprox(Z(zafSi...), 1.1345, atol = 0.001)
        @test isapprox(Z(zafMg...), 1.1280, atol = 0.001)
        @test_broken isapprox(Z(zafBa...), 0.8330, atol = 0.001)
        @test isapprox(Z(zafTi...), 0.9446, atol = 0.001)
        @test isapprox(Z(zafZn...), 0.8973, atol = 0.001)
        @test isapprox(Z(zafZr...), 0.8442, atol = 0.001)
        @test isapprox(Z(zafO...), 1.1316, atol = 0.001)

        @test isapprox(A(zafMg..., θ, θ), 0.5996, atol = 0.001)
        @test_broken isapprox(A(zafBa..., θ, θ), 1.0172, atol = 0.001)
        @test isapprox(A(zafTi..., θ, θ), 0.9615, atol = 0.001)
        @test isapprox(A(zafZn..., θ, θ), 0.9841, atol = 0.001)
        @test_broken isapprox(A(zafZr..., θ, θ), 0.7596, atol = 0.001)
        @test isapprox(A(zafO..., θ, θ), 0.7750, atol = 0.001)

        @test isapprox(F(zafSi..., θ, θ), 1.0015, atol = 0.001) # 1.0026
        @test isapprox(F(zafMg..., θ, θ), 1.0041, atol = 0.001) # 1.0033
        @test isapprox(F(zafBa..., θ, θ), 0.9999, atol = 0.001) # 1.0043
        @test isapprox(F(zafTi..., θ, θ), 1.0042, atol = 0.002) # 1.0009
        @test isapprox(F(zafZn..., θ, θ), 1.000, atol = 0.001)  # 1.000
        @test isapprox(F(zafZr..., θ, θ), 1.0015, atol = 0.001) # 1.0106
        @test isapprox(F(zafO..., θ, θ), 0.9996, atol = 0.001)  # 1.0006
        #    print(asa(DataFrame, Dict( [ zafSi, zafMg, zafBa, zafTi, zafZn, zafZr, zafO ]), θ))
    end

    @testset "Quantify K458" begin
        # Material	K458 = [O(0.3193 mass frac),Si(0.2308 mass frac),Zn(0.0307 mass frac),Ba(0.4192 mass frac),Σ=1.0000,3.5 g/cc]
        # Detector	Bruker 5 eV/ch - FWHM[Mn Kα]=128.6 eV - 2019-02-27 00:00
        # Algorithm	XPP - Pouchou & Pichoir Simplified
        # MAC	NIST-Chantler 2005
        # E0	15 keV
        # Take-off	40°

        # IUPAC	Seigbahn	Standard	Energy	 ZAF	  Z	  A	  F	k-ratio
        # O K-L3	O Kα1	SiO2	0.5249	1.2446	1.1591	1.0743	0.9996	0.746227
        # Si K-L3	Si Kα1	SiO2	1.7397	0.8937	1.1635	0.7672	1.0012	0.441263
        # Si K-M3	Si Kβ1	SiO2	1.8290	0.9150	1.1635	0.7854	1.0012	0.451796
        # Zn K-L3	Zn Kα1	Pure Zn	8.6389	0.9048	0.9201	0.9833	1.0000	0.027776
        # Zn K-M3	Zn Kβ1	Pure Zn	9.5720	0.9083	0.9201	0.9872	1.0000	0.027884
        # Zn L3-M5	Zn Lα1	Pure Zn	1.0116	0.4506	0.9261	0.4860	1.0011	0.013833
        # Ba L3-M5	Ba Lα1	BaCl	4.4663	0.8490	0.8175	1.0384	1.0001	0.447794
        lbl = label("K458")
        unkProps = Dict(:BeamEnergy=>15.0e3, :TakeOffAngle=>deg2rad(40.0))
        stdProps = unkProps # Same for both (in this case...)
        krs = [
            KRatio([n"O K-L3"], unkProps, stdProps, mat"SiO2", 0.746227 ),
            KRatio([n"Si K-L3"], unkProps, stdProps, mat"SiO2", 0.441263 ),
            KRatio([n"Zn K-L3"], unkProps, stdProps, mat"Zn", 0.027776 ),
            KRatio([n"Ba L3-M5"], unkProps, stdProps, mat"BaCl", 0.447794 )
        ]
        res = quantify(lbl, krs)
        @test res.converged
        @test isapprox(value(res.comp[n"Si"]), 0.2308, atol=0.0004)
        @test isapprox(value(res.comp[n"Ba"]), 0.4192, atol=0.0021)
        @test isapprox(value(res.comp[n"O"]), 0.3193, atol=0.0001)
        @test isapprox(value(res.comp[n"Zn"]), 0.0307, atol=0.0001)
    end

end

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
using Random

@testset "Orbitals" begin
    @testset "kappa" begin
        import AtomicLevels: assert_orbital_ℓj
        @test assert_orbital_ℓj(0, HalfInteger(1//2)) === nothing
        @test assert_orbital_ℓj(1, HalfInteger(1//2)) === nothing
        @test assert_orbital_ℓj(1, 3//2) === nothing
        @test assert_orbital_ℓj(2, 2.5) === nothing
        @test_throws ArgumentError assert_orbital_ℓj(0, 0)
        @test_throws ArgumentError assert_orbital_ℓj(0, 1)
        @test_throws ArgumentError assert_orbital_ℓj(0, 3//2)
        @test_throws ArgumentError assert_orbital_ℓj(5, 1//2)
        @test_throws MethodError assert_orbital_ℓj(HalfInteger(1), HalfInteger(1//2))

        import AtomicLevels: κ2ℓ
        @test_throws ArgumentError κ2ℓ(0)
        @test κ2ℓ(-1) == 0
        @test κ2ℓ( 1) == 1
        @test κ2ℓ(-2) == 1
        @test κ2ℓ( 2) == 2
        @test κ2ℓ(-3) == 2
        @test κ2ℓ( 3) == 3

        import AtomicLevels: κ2j
        @test_throws ArgumentError κ2j(0)
        @test κ2j(-1) === half(1)
        @test κ2j( 1) === half(1)
        @test κ2j(-2) === half(3)
        @test κ2j( 2) === half(3)
        @test κ2j(-3) === half(5)
        @test κ2j( 3) === half(5)

        import AtomicLevels: ℓj2κ
        @test ℓj2κ(0, half(1)) == -1
        @test str2κ("s") == -1
        @test κ"s" == -1
        @test ℓj2κ(1, half(1)) == 1
        @test str2κ("p-") == 1
        @test κ"p-" == 1
        @test ℓj2κ(1, 3//2) == -2
        @test str2κ("p") == -2
        @test κ"p" == -2
        @test ℓj2κ(2, 3//2) == 2
        @test str2κ("d-") == 2
        @test κ"d-" == 2
        @test ℓj2κ(2, 5//2) == -3
        @test str2κ("d") == -3
        @test κ"d" == -3
        @test ℓj2κ(3, 5//2) == 3
        @test str2κ("f-") == 3
        @test κ"f-" == 3
        @test ℓj2κ(3, 7//2) == -4
        @test str2κ("f") == -4
        @test κ"f" == -4
        @test_throws ArgumentError ℓj2κ(0, half(3))
        @test_throws ArgumentError ℓj2κ(0, 0)
        @test_throws ArgumentError ℓj2κ(6, 1//2)
    end

    @testset "Construction" begin
        @test o"1s"   == Orbital(1, 0)
        @test o"2p"   == Orbital(2, 1)
        @test o"2[1]" == Orbital(2, 1)

        @test parse(Orbital, "1s") == Orbital(1, 0)
        @test parse(Orbital{Int}, "1s") == Orbital(1, 0)
        @test parse(Orbital{Symbol}, "ks") == Orbital(:k, 0)

        @test ro"1s"   == RelativisticOrbital(1, -1) # κ=-1 => s orbital
        @test ro"2p-"  == RelativisticOrbital(2,  1, half(1))
        @test ro"2p-"  == RelativisticOrbital(2,  1, 1//2)
        @test ro"2p"   == RelativisticOrbital(2,  1, 3//2)
        @test ro"2[1]" == RelativisticOrbital(2,  1, 3//2)

        @test o"kp" == Orbital(:k,1)
        @test o"ϵd" == Orbital(:ϵ,2)

        @test ro"kp"  == RelativisticOrbital(:k, 1, 3//2)
        @test ro"ϵd-" == RelativisticOrbital(:ϵ, 2, 3//2)

        @test_throws ArgumentError parse(Orbital, "2p-")
        @test_throws ArgumentError parse(Orbital, "sdkfl")

        @test_throws ArgumentError Orbital(0, 0)
        @test_throws ArgumentError Orbital(1, 1)
        @test_throws ArgumentError RelativisticOrbital(0, 0)
        @test_throws ArgumentError RelativisticOrbital(1, 1)
        @test_throws ArgumentError RelativisticOrbital(1, 0, 3//2)

        @test nonrelorbital(o"2p") == o"2p"
        @test nonrelorbital(ro"2p") == o"2p"
        @test nonrelorbital(ro"2p-") == o"2p"
    end

    @testset "Properties" begin
        let o = o"3d"
            @test o.n == 3
            @test o.ℓ == 2
        end
        let o = o"ks"
            @test o.n == :k
            @test o.ℓ == 0
        end

        @test propertynames(o"1s") == (:n , :ℓ)
        @test propertynames(ro"1s") == (:n, :κ, :j, :ℓ)

        let o = ro"3d"
            @test o.n == 3
            @test o.κ == -3
            @test o.j == 5//2
            @test o.ℓ == 2
        end
        let o = ro"ks"
            @test o.n == :k
            @test o.κ == -1
            @test o.j == 1//2
            @test o.ℓ == 0
        end
        let o = ro"2p-"
            @test o.n == 2
            @test o.κ == 1
            @test o.j == 1//2
            @test o.ℓ == 1
        end
    end

    @testset "Order" begin
        @test sort(shuffle([o"1s", o"2s", o"2p", o"3s", o"3p"])) ==
            [o"1s", o"2s", o"2p", o"3s", o"3p"]
        @test sort([o"ls", o"kp", o"2p", o"1s"]) ==
            [o"1s", o"2p", o"kp", o"ls"]

        @test sort(shuffle([ro"1s", ro"2s", ro"2p-", ro"2p", ro"3s", ro"3p-", ro"3p"])) ==
            [ro"1s", ro"2s", ro"2p-", ro"2p", ro"3s", ro"3p-", ro"3p"]
        @test sort([ro"ls", ro"kp", ro"2p", ro"1s"]) ==
            [ro"1s", ro"2p", ro"kp", ro"ls"]
    end

    @testset "Range of orbitals" begin
        @test os"6[s-d] 5[d]" == [o"5d", o"6s", o"6p", o"6d"]
        @test os"k[s-g] l[s-g]" == [o"ks", o"kp", o"kd", o"kf", o"kg",
                                    o"ls", o"lp", o"ld", o"lf", o"lg"]

        @test ros"6[s-d] 5[d]" == [ro"5d-", ro"5d", ro"6s", ro"6p-", ro"6p", ro"6d-", ro"6d"]
        @test ros"k[s-g] l[s-g]" == [ro"ks", ro"kp-", ro"kp", ro"kd-", ro"kd", ro"kf-", ro"kf", ro"kg-", ro"kg",
                                     ro"ls", ro"lp-", ro"lp", ro"ld-", ro"ld", ro"lf-", ro"lf", ro"lg-", ro"lg"]
    end

    @testset "Flip j" begin
        @test AtomicLevels.flip_j(ro"1s") == ro"1s"
        @test AtomicLevels.flip_j(ro"2p-") == ro"2p"
        @test AtomicLevels.flip_j(ro"2p") == ro"2p-"
        @test AtomicLevels.flip_j(ro"3d-") == ro"3d"
        @test AtomicLevels.flip_j(ro"3d") == ro"3d-"

        @test AtomicLevels.flip_j(ro"kd-") == ro"kd"
        @test AtomicLevels.flip_j(ro"kd") == ro"kd-"
    end

    @testset "Degeneracy" begin
        @test degeneracy(o"1s") == 2
        @test degeneracy(o"2p") == 6
        @test degeneracy(o"3d") == 10

        @test degeneracy(o"kp") == 6

        @test degeneracy(ro"1s") == 2
        @test degeneracy(ro"2p-") == 2
        @test degeneracy(ro"2p") == 4
        @test degeneracy(ro"3d-") == 4
        @test degeneracy(ro"3d") == 6

        @test degeneracy(ro"kp-") == 2
        @test degeneracy(ro"kp") == 4
    end

    @testset "Parity" begin
        @test iseven(parity(o"1s"))
        @test isodd(parity(o"2p"))
        @test iseven(parity(o"3s"))
        @test isodd(parity(o"3p"))
        @test iseven(parity(o"3d"))

        @test isodd(parity(o"kp"))

        @test iseven(parity(ro"1s"))
        @test isodd(parity(ro"2p"))
        @test iseven(parity(ro"3s"))
        @test isodd(parity(ro"3p"))
        @test iseven(parity(ro"3d"))

        @test isodd(parity(ro"kp"))
    end

    @testset "Symmetry" begin
        @test symmetry(o"1s") == 0
        @test symmetry(o"2s") == symmetry(o"1s")
        @test symmetry(o"2s") != symmetry(o"2p")

        @test symmetry(ro"1s") == symmetry(ro"2s")
        @test symmetry(ro"2p") == symmetry(ro"3p")
        @test symmetry(ro"2p") != symmetry(ro"2p-")
    end

    @testset "Bound" begin
        @test isbound(o"1s")
        @test isbound(o"3d")
        @test isbound(ro"1s")
        @test isbound(ro"3d")
        @test !isbound(o"ks")
        @test !isbound(o"ϵd")
        @test !isbound(ro"ks")
        @test !isbound(ro"ϵd")
    end

    @testset "Angular momenta" begin
        @test angular_momenta(o"2s") == (0,half(1))
        @test angular_momenta(o"2p") == (1,half(1))
        @test angular_momenta(o"4f") == (3,half(1))
        @test angular_momentum_ranges(o"4f") == (-3:3,-half(1):half(1))

        @test angular_momenta(ro"1s") == (half(1),)
        @test angular_momenta(ro"2p-") == (half(1),)
        @test angular_momenta(ro"2p") == (half(3),)
        @test angular_momenta(ro"3d-") == (half(3),)
        @test angular_momenta(ro"3d") == (half(5),)
        @test angular_momentum_ranges(ro"3d") == (-half(5):half(5),)
    end

    @testset "Spin orbitals" begin
        up, down = half(1),-half(1)

        @test_throws ArgumentError SpinOrbital(o"1s", 1, up)
        @test_throws ArgumentError SpinOrbital(o"ks", 1, up)
        @test_throws ArgumentError SpinOrbital(o"2p", -3, up)
        @test_throws ArgumentError SpinOrbital(o"2p", 1, half(3))
        @test_throws ArgumentError SpinOrbital(ro"2p-", half(3))

        soα = SpinOrbital(o"1s", 0, up)
        soβ = SpinOrbital(o"1s", 0, down)

        po₋α = SpinOrbital(o"2p", -1, up)
        po₀α = SpinOrbital(o"2p", 0, up)
        po₊α = SpinOrbital(o"2p", 1, up)
        po₋β = SpinOrbital(o"2p", -1, down)
        po₀β = SpinOrbital(o"2p", 0, down)
        po₊β = SpinOrbital(o"2p", 1, down)

        @test degeneracy(soα) == 1
        @test soα < soβ

        @test parity(soα) == parity(soβ) == p"even"

        @test parity(po₋α) == p"odd"
        @test parity(po₀α) == p"odd"
        @test parity(po₊α) == p"odd"
        @test parity(po₋β) == p"odd"
        @test parity(po₀β) == p"odd"
        @test parity(po₊β) == p"odd"

        @test symmetry(soα) != symmetry(soβ)

        @test symmetry(po₋α) != symmetry(po₀α)
        @test symmetry(po₋α) != symmetry(po₊α)
        @test symmetry(po₊α) != symmetry(po₀α)

        @test symmetry(po₋β) != symmetry(po₀β)
        @test symmetry(po₋β) != symmetry(po₊β)
        @test symmetry(po₊β) != symmetry(po₀β)

        @test symmetry(po₋α) != symmetry(po₋β)
        @test symmetry(po₀α) != symmetry(po₀β)
        @test symmetry(po₊α) != symmetry(po₊β)

        @test symmetry(po₋α) == symmetry(SpinOrbital(o"3p", -1, up))

        @test isbound(soα)
        @test !isbound(SpinOrbital(o"ks", 0, up))

        @test spin_orbitals(o"1s") == [soα, soβ]
        @test spin_orbitals(o"2p") == [po₋α, po₋β, po₀α, po₀β, po₊α, po₊β]
        for orb in [o"1s", o"2p", o"3d", o"4f", o"5g"]
            @test length(spin_orbitals(orb)) == degeneracy(orb)
        end

        @test sos"1[s] 2[p]" == [soα, soβ, po₋α, po₋β, po₀α, po₀β, po₊α, po₊β]

        @test "$(soα)" == "1s₀α"
        @test "$(po₊β)" == "2p₁β"

        @testset "Construction" begin
            @test so"1s(0,α)" == SpinOrbital(o"1s", (0,1/2))
            @test so"1s(0,β)" == SpinOrbital(o"1s", (0,-1/2))
            @test so"2p(1,β)" == SpinOrbital(o"2p", (1,-1/2))
            @test so"1s(0,-1/2)" == SpinOrbital(o"1s", (0,-1/2))
            @test so"1s(0,-0.5)" == SpinOrbital(o"1s", (0,-1/2))
            @test so"1s(0,1/2)" == SpinOrbital(o"1s", (0,1/2))
            @test_throws ArgumentError parse(SpinOrbital{RelativisticOrbital}, "1s(0,1/2)")
            @test rso"1s(1/2)" == SpinOrbital(ro"1s", (1/2))
            @test rso"1s(α)" == SpinOrbital(ro"1s", (1/2)) # This should not work, but hey...
            @test_throws ArgumentError parse(SpinOrbital{RelativisticOrbital}, "1s(3/2)")
            @test rso"2p(3/2)" == SpinOrbital(ro"2p", (3/2))
            @test rso"3p(-3/2)" == SpinOrbital(ro"3p", (-3/2))
            @test rso"3p(-1/2)" == SpinOrbital(ro"3p", (-1/2))
            @test_throws ArgumentError parse(SpinOrbital{RelativisticOrbital}, "3p(5/2)")
            @test rso"3p(1/2)" == SpinOrbital(ro"3p", (1/2))
            @test rso"3p-(1/2)" == SpinOrbital(ro"3p-", (1/2))
            @test_throws ArgumentError parse(SpinOrbital{RelativisticOrbital}, "3p-(3/2)")
            @test_throws ArgumentError parse(SpinOrbital{RelativisticOrbital}, "3p-(-3/2)")
            @test rso"3p-(-1/2)" == SpinOrbital(ro"3p-", (-1/2))
            @test_throws ArgumentError parse(SpinOrbital{RelativisticOrbital}, "3d-(-5/2)")
            @test_throws ArgumentError parse(SpinOrbital{RelativisticOrbital}, "3d-(..)")
            @test rso"3d-(-3/2)" == SpinOrbital(ro"3d-", (-3/2))
            @test rso"3d(5/2)" == SpinOrbital(ro"3d", (5/2))

            @test parse(SpinOrbital{Orbital}, "1s(0,α)") == SpinOrbital(Orbital(1, 0), (0, up))
            @test parse(SpinOrbital{Orbital{Int}}, "1s(0,α)") == SpinOrbital(Orbital(1, 0), (0, up))
            @test parse(SpinOrbital{Orbital{Symbol}}, "ks(0,α)") == SpinOrbital(Orbital(:k, 0), (0, up))

            @test so"1s₀β" == so"1s(0,β)"
            @test so"2p₋₁α" == so"2p(-1,α)"
            @test so"k[31]₋₁₃α" == so"k[31](-13,α)"

            for o in [SpinOrbital(o"1s", (0,-1/2)),
                      SpinOrbital(o"1s", (0,1/2)),
                      SpinOrbital(o"2p", (1,-1/2)),
                      SpinOrbital(ro"1s", (1/2)),
                      SpinOrbital(ro"2p", (3/2)),
                      SpinOrbital(ro"3d", (5/2)),
                      SpinOrbital(ro"3d-", (-3/2)),
                      SpinOrbital(ro"3p", (-1/2)),
                      SpinOrbital(ro"3p", (-3/2)),
                      SpinOrbital(ro"3p", (1/2)),
                      SpinOrbital(ro"3p-", (-1/2)),
                      SpinOrbital(ro"3p-", (1/2))]
                O = typeof(o.orb)
                @test parse(SpinOrbital{O}, string(o)) == o
            end
        end
    end

    @testset "Internal methods" begin
        @test AtomicLevels.mqtype(o"2s") == Int
        @test AtomicLevels.mqtype(o"ks") == Symbol
        @test AtomicLevels.mqtype(ro"2s") == Int
        @test AtomicLevels.mqtype(ro"ks") == Symbol
    end

    @testset "String representation" begin
        @test string(o"1s") == "1s"
        @test string(ro"2p-") == "2p-"
        @test ascii(ro"2p-") == "2p-"
    end

    @testset "Hashing" begin
        @test hash(o"3s") == hash(o"3s")
        @test hash(ro"3p-") == hash(ro"3p-")
    end

    @testset "Serialization" begin
        @testset "Orbitals" begin
            o = o"1s"
            p = o"kg"
            q = Orbital(:k̃, 14)
            oo = SpinOrbital(o, (0, 1/2))
            r = Orbital(Symbol("[$(oo)]"), 14)

            no,np,nq,nr = let io = IOBuffer()
                foreach(Base.Fix1(write, io), (o,p,q,r))

                seekstart(io)
                [read(io, Orbital)
                 for i = 1:4]
            end

            @test no == o
            @test np == p
            @test nq == q
            @test nr == r
        end

        @testset "Relativistic orbitals" begin
            o = ro"1s"
            p = ro"kg"
            q = RelativisticOrbital(:k̃, 14)
            oo = SpinOrbital(o, (1/2))
            r = RelativisticOrbital(Symbol("[$(oo)]"), 14)

            no,np,nq,nr = let io = IOBuffer()
                foreach(Base.Fix1(write, io), (o,p,q,r))

                seekstart(io)
                [read(io, RelativisticOrbital)
                 for i = 1:4]
            end

            @test no == o
            @test np == p
            @test nq == q
            @test nr == r
        end

        @testset "Spin-orbitals" begin
            a = so"1s(0,α)"
            b = rso"2p-(1/2)"
            c = rso"kd-(-1.5)"
            d = so"3d(-2,-0.5)"
            e = SpinOrbital(Orbital(Symbol("[$(d)]"), 14), (-13, -0.5))

            na,nb,nc,nd,ne = let io = IOBuffer()
                foreach(Base.Fix1(write, io), (a,b,c,d,e))

                seekstart(io)
                [read(io, SpinOrbital)
                 for i = 1:5]
            end

            @test na == a
            @test nb == b
            @test nc == c
            @test nd == d
            @test ne == e
        end
    end

    @testset "Broadcasting" begin
        @test ([o"1s", o"2p"] .== o"1s") == [true, false]
        @test ([ro"1s", ro"2p-"] .== ro"1s") == [true, false]
        @test ([so"1s(0,α)", so"2p(0,α)"] .== so"1s(0,α)") == [true, false]
        @test ([rso"1s(1/2)", rso"2p-(1/2)"] .== rso"1s(1/2)") == [true, false]
    end
end


text:
using LinearAlgebra
using Plots

# generate the data for the example
function generateSamples(n, k, mu::Float64)
  # generate theta
  theta = rand(Float64, k)
  @. theta = (theta -0.5) * 10

  # generate y with noise
  y = Vector{Float64}(undef, n+k -1)
  for j in 1:n
    y[j] = 1 * (rand(Float64) -0.5)
  end
  for j in n+1:n+k-1
    y[j] = 0
  end

  # generate X, update y
  supportPoints = rand(Float64, n)
  @. supportPoints = (supportPoints - 0.5)
  supportPoints = sort(supportPoints)
  X = Matrix{Float64}(undef, k, n+k-1)

  for j in 1:n
    X[1,j] = 1
    y[j] = theta[1]
    for i in 2:k
      X[i,j] = X[i-1, j] * supportPoints[j]
      y[j] += X[i,j] * theta[i]
    end
  end

  # identidy matrix multiplied with vector of weights, first row
  for j in n+1:n+k-1
    X[1,j] = 0
    for i in 1:k
      if (j - n + 1 == i)
        X[i,j] = sqrt(mu * (j-n))
      end
    end
  end

  return X, y, theta
end

function evalPolynomial(grid, coeff)
  k = length(coeff)
  len = length(grid)
  f = zeros(Float64, len)
  for i in 1:k
    for j in 1:len
      f[j] *= grid[j]
      f[j] += coeff[k + 1 - i]
    end
  end
  return f
end

function plotAll(X::Matrix{Float64}, y::Vector{Float64}, theta::Vector{Float64}, n, k)
  x = Vector{Float64}(undef, n)
  val = Vector{Float64}(undef, n)
  for j in 1:n
    x[j] = X[2,j]
    val[j] = y[j]
  end
  plot(x, val ,label="", seriestype=:scatter) #xlims=[x[1], x[n]]

  grid = collect(range(x[1], step=0.01, stop=x[n]))
  f = evalPolynomial(grid, theta)
  plot!(grid, f, label="theta")

  theta2 = X' \ y # solve in least squares sense
  f2 = evalPolynomial(grid, theta2)
  plot!(grid, f2, label="regression theta")
end

function wrap(n, k, mu::Float64)
  X, y, theta = generateSamples(n, k, mu)
  plotAll(X, y, theta, n, k)
end

max_stars_repo_name: Mathemalsky/BenignOverfitting

max_issues_repo_name: Mathemalsky/BenignOverfitting

max_forks_repo_name: Mathemalsky/BenignOverfitting


text:
using DataDeps

register(DataDep(
  "beer-aroma",
  "http://www.stat.ufl.edu/~winner/data/aroma_beer.txt",
  "http://www.stat.ufl.edu/~winner/data/aroma_beer.dat",
  "eaf147c291030154be364fe408dc0d01bac0407482eb91e965aca111abc0b977",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=9,
        feature_indices=2:8,
        id_indices=1,
        delim=' ',
        ignorerepeated=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
include("../utils.jl")
include("group_utils.jl")

using Manifolds: invariant_metric_dispatch, default_metric_dispatch

@testset "Circle group" begin
    G = CircleGroup()
    @test repr(G) == "CircleGroup()"

    @test base_manifold(G) === Circle{ℂ}()

    @test (@inferred invariant_metric_dispatch(G, LeftAction())) === Val(true)
    @test (@inferred invariant_metric_dispatch(G, RightAction())) === Val(true)
    @test (@inferred Manifolds.biinvariant_metric_dispatch(G)) === Val(true)
    @test (@inferred default_metric_dispatch(MetricManifold(G, EuclideanMetric()))) ===
          Val(true)
    @test has_invariant_metric(G, LeftAction())
    @test has_invariant_metric(G, RightAction())
    @test has_biinvariant_metric(G)
    @test is_default_metric(MetricManifold(G, EuclideanMetric()))

    @testset "identity overloads" begin
        ig = Identity(G, [Complex(1.0)])
        @test identity(G, ig) === ig
        @test inv(G, ig) === ig
        @test allocate_result(G, get_coordinates, ig, Complex(1.0), DefaultBasis()) isa
              Array{Complex{Float64},1}
        y = [Complex(0.0)]
        @test identity!(G, y, [Complex(1.0)]) === y
        @test y == [Complex(1.0)]
        y = [1.0 * im]
        v = [Complex(0.5)]
        @test translate_diff(G, ig, y, v) === v
    end

    @testset "scalar points" begin
        pts = [1.0 + 0.0im, 0.0 + 1.0im, (1.0 + 1.0im) / √2]
        vpts = [0.0 + 0.5im]
        @test compose(G, pts[2], pts[1]) ≈ pts[2] * pts[1]
        @test translate_diff(G, pts[2], pts[1], vpts[1]) ≈ pts[2] * vpts[1]
        test_group(
            G,
            pts,
            vpts,
            vpts;
            test_diff=true,
            test_mutating=false,
            test_invariance=true,
        )
    end

    @testset "vector points" begin
        pts = [[1.0 + 0.0im], [0.0 + 1.0im], [(1.0 + 1.0im) / √2]]
        vpts = [[0.0 + 0.5im]]
        @test compose(G, pts[2], pts[1]) ≈ pts[2] .* pts[1]
        @test translate_diff(G, pts[2], pts[1], vpts[1]) ≈ pts[2] .* vpts[1]
        test_group(
            G,
            pts,
            vpts,
            vpts;
            test_diff=true,
            test_mutating=true,
            test_invariance=true,
        )
    end

    @testset "Group forwards to decorated" begin
        pts = [[1.0 + 0.0im], [0.0 + 1.0im], [(1.0 + 1.0im) / √2]]
        test_manifold(
            G,
            pts,
            test_forward_diff=false,
            test_reverse_diff=false,
            test_vector_spaces=false,
            test_project_tangent=true,
            test_musical_isomorphisms=false,
            test_default_vector_transport=true,
            is_mutating=true,
            exp_log_atol_multiplier=2.0,
            is_tangent_atol_multiplier=2.0,
        )
    end
end


text:
# Perform Sanity Tests for Agents, Arm models and Experiments

# @testset "Sanity Checks" begin
#
# end


text:
using PyCall
using Conda

if PyCall.conda
	Conda.add("pip")
	pip = joinpath(Conda.SCRIPTDIR, "pip")
    run(`$pip install python-keystoneclient`)
	run(`$pip install python-swiftclient`)
else
	try
		pyimport("keystoneclient")
		pyimport("swiftclient")
	catch ee
		typeof(ee) <: PyCall.PyError || rethrow(ee)
		warn("""
Python Dependancies not installed
Please either:

 - Rebuild PyCall to use Conda, by running in the julia REPL:
    - `ENV["PYTHON"]=""; Pkg.build("PyCall"); Pkg.build("SwiftObjectStores")`
 - Or install the depencences yourself, eg by running pip
	- `pip install python-keystoneclient`
	- `pip install python-swiftclient`
	"""
		)
	end
end


text:
module Rifraf

using BioSequences
using Distributions
using Parameters
using StatsBase

import Base.length
import Base.reverse

export rifraf,
       RifrafParams,
       RifrafResult,
       ErrorModel,
       Scores,
       normalize,
       sample_sequences,
       write_samples,
       read_samples,
       sample_mixture,
       RifrafSequence,
       length,
       reverse,
       DNASeq,
       Phred,
       Prob,
       LogProb,
       Score,
       BandedArray,
       sparsecol,
       data_row,
       row_range,
       data_row_range,
       inband,
       flip!,
       newbandwidth!,
       align,
       calibrate_phreds,
       read_fastq_records,
       read_fastq,
       write_fastq,
       read_fasta_records,
       read_fasta,
       write_fasta

include("bandedarrays.jl")
include("types.jl")
include("phred.jl")
include("util.jl")
include("fastxio.jl")
include("proposals.jl")
include("errormodel.jl")
include("rifrafsequences.jl")
include("align.jl")
include("model.jl")
include("sample.jl")

end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module AbstractSimulatorImpl

export AbstractSimulator

abstract AbstractSimulator

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


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

name = "Elfutils"
version = v"0.177"

# Collection of sources required to build Elfutils
sources = [
    ArchiveSource("https://sourceware.org/elfutils/ftp/$(version.major).$(version.minor)/elfutils-$(version.major).$(version.minor).tar.bz2",
                  "fa489deccbcae7d8c920f60d85906124c1989c591196d90e0fd668e3dc05042e"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/elfutils-*/
export CC=gcc
export CXX=g++
CFLAGS="-Wno-error=unused-result" CPPFLAGS="-I${prefix}/include" ./configure --prefix=${prefix} --host=${target}
make -j${nproc}
make install
"""

# Only build for Linux, and disable musl/FreeBSD until we have a custom libargp, perhaps
# the https://github.com/xhebox/libuargp which is known to be musl-friendly.
platforms = [p for p in supported_platforms() if p isa Linux && libc(p) != :musl]

# The products that we will ensure are always built
products = [
    LibraryProduct("libasm", :libasm),
    LibraryProduct("libdw", :libdw),
    LibraryProduct("libelf", :libelf),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "Zlib_jll",
    "Bzip2_jll",
    "XZ_jll",
]

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


text:
module IEEEFloats

export IEEEFloat, ieeefloat,

       sign_mask, significand_mask, exponent_mask,
       sign_bits, significand_bits, exponent_bits,

       exponent_bias, exponent_biasing, exponent_unbiasing,
       unbiased_exponent_max, unbiased_exponent_min,
       biased_exponent_max, biased_exponent_min,

       get_sign, get_exponent, get_signficand,
       set_sign, set_exponent, set_signficand,
       isolate_sign, isolate_exponent, isolate_significand,
       clear_sign, clear_exponent, clear_significand

using Base: IEEEFloat, precision, significand_mask, exponent_mask

# compat for Julia v1.0.5 Base: exponent_bits, signnificand_bits
include("v105compatible.jl")

# `f(xT)` denotes a function `f` that accepts either a value `x::T` or a type `T

sign_bits(::Type{T}) where T<:IEEEFloat = 1

include("support_more_types.jl")  # Float128, Float8
include("corresponding_sizes.jl") # ieeefloat(xT), unsigned(xT), signed(xT)

bitwidth(::Type{T}) where T = sizeof(T) * 8

exponent_bias(n_exponent_bits) = 2^(n_exponent_bits - 1) - 1
exponent_bias(::Type{T}) where T = exponent_bias(exponent_bits(T))

# unbiased_exponent_max(T) == exponent_bias(T)
unbiased_exponent_max(::Type{T}) where T = 2^(exponent_bits(T) - 1) - 1  # for normal values
unbiased_exponent_min(::Type{T}) where T = 1 - exponent_max(T)           # for normal values

biased_exponent_max(::Type{T}) where T = 2^(exponent_bits(T)) - 1
biased_exponent_min(::Type{T}) where T = 0

exponent_biasing(::Type{T}, unbiased) where T = unbiased + exponent_bias(T)
exponent_unbiasing(::Type{T}, biased) where T = biased - exponent_bias(T)

include("use_fields.jl") # get_x_field, set_x_field, isolate_x_field, clear_x_field

# ============== #

decode_normal_value(nexpbits, nsigbits, significand, unbiasedexponent=0, isneg::Bool=false) =
    (isneg ? -1 : 1) *
    (1 + (significand / 2^nsigbits)) *
    (2.0^(unbiasedexponent - exponent_bias(nexpbits)))
decode_subnormal_value(nexpbits, nsigbits, significand, unbiasedexponent=0, isneg::Bool=false) =
    (isneg ? -1 : 1) *
    (significand) *
    (2.0^(1 - exponent_bias(nexpbits) - nsigbits))

 #=
const bias128 = 16_383 # 15 bits
const bias64 = exponent_bias(exponent_bits(Float64)) # 1023
const bias32 = exponent_bias(exponent_bits(Float32)) # 127
const bias16 = exponent_bias(exponent_bits(Float16)) # 15
const bias8 = exponent_bias(3) # 3

exponent_bias(::Type{T}) where T<:IEEEFloat = exponent_max(T)
exponent_max(::Type{T}) where T =   2^(exponent_bits(T) - 1) - 1  # normal values
exponent_min(::Type{T}) where T =  1 - exponent_max(T)             # normal values

biased_exponentmax(::Type{T}) where T = 2^(exponent_bits(T)) - 1
biased_exponentmin(::Type{T}) where T = 0
=#

end # module IEEEFloats


text:
module StochasticArithmetic

using Statistics

export UP, DWN, RND, SFloat64, SFloat32, value, @reliable_digits

include("EFT.jl")
include("determinism.jl")
include("rounding.jl")
include("sfloat.jl")

end # module


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

mutable struct IoK8sApiAppsV1RollingUpdateStatefulSetStrategy <: SwaggerModel
    partition::Any # spec type: Union{ Nothing, Int32 } # spec name: partition

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

const _property_map_IoK8sApiAppsV1RollingUpdateStatefulSetStrategy = Dict{Symbol,Symbol}(Symbol("partition")=>Symbol("partition"))
const _property_types_IoK8sApiAppsV1RollingUpdateStatefulSetStrategy = Dict{Symbol,String}(Symbol("partition")=>"Int32")
Base.propertynames(::Type{ IoK8sApiAppsV1RollingUpdateStatefulSetStrategy }) = collect(keys(_property_map_IoK8sApiAppsV1RollingUpdateStatefulSetStrategy))
Swagger.property_type(::Type{ IoK8sApiAppsV1RollingUpdateStatefulSetStrategy }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1RollingUpdateStatefulSetStrategy[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1RollingUpdateStatefulSetStrategy }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1RollingUpdateStatefulSetStrategy[property_name]

function check_required(o::IoK8sApiAppsV1RollingUpdateStatefulSetStrategy)
    true
end

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


text:
using GenTeX
using Test

@testset "LaTeX" begin
    maths = [raw"$a$", raw"\(b\)", raw"$$ c = d $$", raw"\[ e = f \]"]
    s = join(maths, ' ')
    matches = map(x -> string(x.match), eachmatch(GenTeX.latex_regex, s))
    @test maths == matches

    tmpdir = tempname() * '/'
    mkdir(tmpdir)
    @test startswith(GenTeX.substitute_latex(raw"$$x_y$$", 1, tmpdir), "<center")
    @test startswith(GenTeX.substitute_latex(raw"$x$", 1, tmpdir), "<img")
    @test startswith(GenTeX.substitute_latex(raw"$$\text{mean} = 1$$", 1, tmpdir), "<center")
    rm(tmpdir, recursive=true)
end


text:
#!/usr/bin/env julia

# run with 5 ranks

using Gasp

dl = Dlog("./dlog"; ranks_per=2, space_per=128)

if grank() == 1
    for i = 1:3
        run(`cat dlog/messages-$i.dlog`)
    end
end

sync()

message(dl, "foobar")

dl = nothing
gc()

if grank() == 1
    for i = 1:3
        run(`cat ./dlog/messages-$i.dlog`)
    end
end


text:
@testset "EarthEnv HabitatHeterogeneity" begin
    using RasterDataSources: rastername, rasterurl, rasterpath

    @test rastername(EarthEnv{HabitatHeterogeneity}, :cv; res="1km") == "cv_1km.tif"
    hh_path = joinpath(ENV["RASTERDATASOURCES_PATH"], "EarthEnv", "habitat_heterogeneity")
    @test rasterpath(EarthEnv{HabitatHeterogeneity}) == hh_path
    @test rasterpath(EarthEnv{HabitatHeterogeneity}, :cv; res="1km") == joinpath(hh_path, "1km", "cv_1km.tif")

    @test rasterurl(EarthEnv{HabitatHeterogeneity}, :cv; res="1km") ==
        URI(scheme="https", host="data.earthenv.org", path="/habitat_heterogeneity/1km/cv_01_05_1km_uint16.tif")
    getraster(EarthEnv{HabitatHeterogeneity}, :cv)
    @test isfile(joinpath(hh_path, "25km", "cv_25km.tif"))
end


text:
module BroadcastableStructs

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end BroadcastableStructs

export BroadcastableStruct, BroadcastableCallable

import Setfield
using ZygoteRules: @adjoint

const constructorof = try
    Setfield.constructorof
catch
    Setfield.constructor_of
end

@inline foldlargs(op, x) = x
@inline foldlargs(op, x1, x2, xs...) = foldlargs(op, op(x1, x2), xs...)
# Unroll by hand (optimization):
@inline foldlargs(op, x1, x2) = op(x1, x2)
@inline foldlargs(op, x1, x2, x3) = op(op(x1, x2), x3)
@inline foldlargs(op, x1, x2, x3, x4) = op(op(op(x1, x2), x3), x4)
@inline foldlargs(op, x1, x2, x3, x4, x5) = op(op(op(op(x1, x2), x3), x4), x5)
@inline foldlargs(op, x1, x2, x3, x4, x5, x6) = op(op(op(op(op(x1, x2), x3), x4), x5), x6)
@inline foldlargs(op, x1, x2, x3, x4, x5, x6, x7) =
    op(op(op(op(op(op(x1, x2), x3), x4), x5), x6), x7)
@inline foldlargs(op, x1, x2, x3, x4, x5, x6, x7, x8) =
    op(op(op(op(op(op(op(x1, x2), x3), x4), x5), x6), x7), x8)
@inline foldlargs(op, x1, x2, x3, x4, x5, x6, x7, x8, x9) =
    op(op(op(op(op(op(op(op(x1, x2), x3), x4), x5), x6), x7), x8), x9)
@inline foldlargs(op, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10) =
    op(op(op(op(op(op(op(op(op(x1, x2), x3), x4), x5), x6), x7), x8), x9), x10)

abstract type BroadcastableStruct end

fieldvalues(obj) = ntuple(i -> getfield(obj, i), nfields(obj))

Broadcast.broadcastable(obj::BroadcastableStruct) =
    Broadcast.broadcasted(constructorof(typeof(obj)), fieldvalues(obj)...)

Base.ndims(T::Type{<:BroadcastableStruct}) =
    mapreduce(ndims, max, fieldtypes(T); init=0)

#=
Base.axes(obj::BroadcastableStruct) = axes(Broadcast.broadcastable(obj))
Base.length(obj::BroadcastableStruct) = prod(length.(axes(obj)))

Base.getindex(obj::BroadcastableStruct, i::Int...) =
    Broadcast.broadcastable(obj)[i...]
=#

abstract type BroadcastableCallable <: BroadcastableStruct end

@inline Broadcast.broadcasted(c::BroadcastableCallable, args...) =
    Broadcast.broadcasted(calling(c), deconstruct(c)..., args...)

@inline deconstruct(obj::T) where T =
    foldlargs((), fieldvalues(obj)...) do fields, x
        if x isa BroadcastableStruct
            (fields..., deconstruct(x)...)
        else
            (fields..., x)
        end
    end

struct StructSchema{Constructor, Fields} end
StructSchema(::Type{Constructor}, Fields::Tuple) where {Constructor} =
    StructSchema{Constructor, Fields}()

# Build a type-level tree reflecting the field locations of
# `BroadcastableStruct`s.
@inline function _structschema(obj::Constructor) where {Constructor}
    fields = foldlargs((), fieldvalues(obj)...) do fields, x
        if x isa BroadcastableStruct
            (fields..., _structschema(x))
        else
            (fields..., nothing)
        end
    end
    return StructSchema(Constructor, fields)
end

@inline _reconstruct(::Type{T}, fields) where T = constructorof(T)(fields...)

@inline reconstruct(f, obj::BroadcastableStruct, allargs...) =
    reconstruct(f, _structschema(obj), allargs...)

@inline function reconstruct(f::F, ::StructSchema{Constructor, Fields}, allargs...) where {F, Constructor, Fields}
    fields, args = foldlargs(((), allargs), Fields...) do (fields, allargs), x
        if x isa StructSchema
            y, rest = reconstruct(f, x, allargs...)
            ((fields..., y), rest)
        else
            ((fields..., allargs[1]), Base.tail(allargs))
        end
    end
    return f(Constructor, fields), args
end

# Manually flatten broadcast to avoid unbroadcast MethodError:
# https://github.com/FluxML/Zygote.jl/issues/313
calling(obj::T) where T =
    let schema = _structschema(obj)
        # Make sure to _not_ reference `obj` inside the closure.  It
        # must not have non-isbit types (like arrays) as its members.
        @inline function(allargs...)
            f, args = reconstruct(_reconstruct, schema, allargs...)
            return f(args...)
        end
    end

@adjoint fieldvalues(obj::T) where T = fieldvalues(obj), function(v)
    (NamedTuple{fieldnames(T)}(v),)
end

end # module


text:
function onehot(y::Vector)
    n_samples = size(y, 1)
    n_features = length(unique(y))
    Y = zeros(n_features, n_samples)

    # all elements of y must be > 0 because they are treated
    # as indices
    assert(all(0 .< y))
    for i in 1:n_samples
        Y[y[i], i] = 1
    end
    return Y
end

maxindices{T}(X::AA{T, 2}) = [indmax(X[:, i]) for i in 1:size(X, 2)]

function accuracy{T}(y::AA{T, 1}, ŷ::AA{T, 1})
    assert(size(y) == size(ŷ))
    sum(y .== ŷ) / length(y)
end

to_tuple_if_int(x::Real) = (x, x)

max_stars_repo_name: IshitaTakeshi/NeuralNetwork

max_issues_repo_name: IshitaTakeshi/NeuralNetwork

max_forks_repo_name: IshitaTakeshi/NeuralNetwork