text:
"""
    meanprop(p::AbstractParticles, prop::Symbol; massprop=:mass, massweighted::Bool=true)

Returns the mean of the specified property over all particles.

The mean can be weighted by the mass.
"""
function meanprop(p::AbstractParticles, prop::Symbol; massprop=:mass, massweighted=true)
    if massweighted
        return _meanprop(p[prop], p[massprop])
    else
        return _meanprop(p[prop])
    end
end

_meanprop(a::Number) = a
_meanprop(a::Number, ::Any) = a
_meanprop(a::Fill) = a[1]
_meanprop(a::Fill, ::AbstractVector) = a[1]
_meanprop(a::Fill, ::Fill) = a[1]
_meanprop(a::AbstractVector) = mean(a)
_meanprop(a::AbstractMatrix) = mean(a; dims=2)
_meanprop(a::AbstractVector, m::AbstractVector) = mean(a, weights(ustrip_lazy(m)))
_meanprop(a::AbstractMatrix, m::AbstractVector) = mean(a, weights(ustrip_lazy(m)); dims=2)
_meanprop(a::AbstractArray, ::Number) = _meanprop(a)
_meanprop(a::AbstractVector, ::Fill) = _meanprop(a)
_meanprop(a::AbstractMatrix, ::Fill) = _meanprop(a)

"""
    sumprop(p::AbstractParticles, prop::Symbol)

Returns the sum of the specified property over all particles.
"""
sumprop(p::AbstractParticles, prop::Symbol) = _sumprop(p, p[prop])

_sumprop(p::AbstractParticles, a::Number) = particle_number(p) * a
_sumprop(::AbstractParticles, a::AbstractVecOrMat) = _sumprop(a)
_sumprop(a::AbstractVector) = sum(a)
_sumprop(a::AbstractMatrix) = sum(a; dims=2)

"""
    meanpos(p::AbstractParticles, prop::Symbol=:pos; massprop=:mass, massweighted=true)

Returns the mean of the particles' positions.

The mean can be weighted by the mass.
"""
meanpos(p::AbstractParticles, prop::Symbol=:pos; massprop=:mass, massweighted=true) =
    meanprop(p, prop; massprop, massweighted)

"""
    meanvel(p::AbstractParticles, prop::Symbol=:pos; massprop=:mass, massweighted=true)

Returns the mean of the particles' velocities.

The mean can be weighted by the mass.
"""
meanvel(p::AbstractParticles, prop::Symbol=:vel; massprop=:mass, massweighted=true) =
    meanprop(p, prop; massprop, massweighted)

"""
    angmom(p::AbstractParticles)

Returns the individual angular momenta of the particles.

If `angmomprop` is a `Symbol`, the already computed angular momenta are returned directly if they exist.

# Keyword Arguments
- `origin=Nothing`: origin of the coordinate system for the computation of the angular momentum
- `velorigin=Nothing`: relative velocity zero-point, oftentimes the velocity of the center of mass
- `posprop=:pos`: property name of the position
- `velprop=:vel`: property name of the velocity
- `massprop=:mass`: property name of the mass
- `angmomprop=nothing`: property name of the angular momentum - the method checks whether it already exists or not. Warning: The method assumes the already saved property to have used the same `origin` and `velorigin`.
"""
function angmom(
    p::AbstractParticles;
    origin=nothing,
    velorigin=nothing,
    posprop=:pos,
    velprop=:vel,
    massprop=:mass,
    angmomprop=nothing,
)
    !isnothing(angmomprop) && haskey(p, angmomprop) && return p[angmomprop]
    return angmom(p[posprop], p[velprop], p[massprop]; origin, velorigin)
end

"""
    angmomtot(p::AbstractParticles)

Returns the total angular momentum of the particles.

If `angmomprop` is a `Symbol`, the total angular momentum is computed from the already existing angular momenta.

# Keyword Arguments
- `origin=Nothing`: origin of the coordinate system for the computation of the angular momentum
- `velorigin=Nothing`: relative velocity zero-point, oftentimes the velocity of the center of mass
- `posprop=:pos`: property name of the position
- `velprop=:vel`: property name of the velocity
- `massprop=:mass`: property name of the mass
- `angmomprop=nothing`: property name of the angular momentum - the method checks whether it already exists or not. Warning: The method assumes the already saved property to have used the same `origin` and `velorigin`.
"""
function angmomtot(
    p::AbstractParticles;
    origin=nothing,
    velorigin=nothing,
    posprop=:pos,
    velprop=:vel,
    massprop=:mass,
    angmomprop=nothing,
)
    !isnothing(angmomprop) && haskey(p, angmomprop) && return dropdims(sumprop(p, angmomprop); dims=2)
    return angmomtot(p[posprop], p[velprop], p[massprop]; origin, velorigin)
end

"""
    angmomtot_stable(p::AbstractParticles)

Returns the total angular momentum of the particles, like [`angmomtot`](@ref).

Uses the stable summation algorithm `Base.sum` for the summation, but is almost three times slower than
[`angmomtot`](@ref).

See [`angmomtot`](@ref) for the available keyword arguments.
"""
function angmomtot_stable(
    p::AbstractParticles;
    origin=nothing,
    velorigin=nothing,
    posprop=:pos,
    velprop=:vel,
    massprop=:mass,
    angmomprop=nothing,
)
    !isnothing(angmomprop) && haskey(p, angmomprop) && return dropdims(sumprop(p, angmomprop); dims=2)
    return angmomtot_stable(p[posprop], p[velprop], p[massprop]; origin, velorigin)
end

"""
    colnorm(a::AbstractMatrix[, origin::AbstractVector])

Returns a new `Vector` with the columnwise norms of `a` around `origin`.

This method does not check for overflow of the squared norm.
"""
function colnorm(a::AbstractMatrix{T}) where {T}
    dst = Vector{T}(undef, size(a, 2))
    return _colnorm_unsafe!(dst, a)
end

function colnorm(a::AbstractMatrix{T}, origin::AbstractVector{T0}) where {T,T0}
    @assert size(a, 1) == length(origin)
    T === T0 || (origin = convert(Vector{T}, origin))
    dst = Vector{T}(undef, size(a, 2))
    return _colnorm_unsafe!(dst, a, origin)
end

"""
    colnorm2(a::AbstractMatrix[, origin:AbstractVector])

Returns a new `Vector` with the columnwise squared norms of `a` around `origin`.
"""
function colnorm2(a::AbstractMatrix{T}) where {T}
    dst = Vector{T}(undef, size(a, 2))
    return _colnorm2_unsafe!(dst, a)
end

function colnorm2(a::AbstractMatrix{Q}) where {T,Q<:Quantity{T}}
    u = unit(Q)^2
    Q² = Quantity{T,dimension(u),typeof(u)}
    dst = Vector{Q²}(undef, size(a, 2))
    return _colnorm2_unsafe!(dst, a)
end

function colnorm2(a::AbstractMatrix{Union{Missing,Q}}) where {T,Q<:Quantity{T}}
    u = unit(Q)^2
    Q² = Quantity{T,dimension(u),typeof(u)}
    dst = Vector{Union{Missing,Q²}}(undef, size(a, 2))
    return _colnorm2_unsafe!(dst, a)
end

function colnorm2(a::AbstractMatrix{T}, origin::AbstractVector{T0}) where {T,T0}
    @assert size(a, 1) == length(origin)
    T === T0 || (origin = convert(Vector{T}, origin))
    dst = Vector{T}(undef, size(a, 2))
    return _colnorm2_unsafe!(dst, a, origin)
end

function colnorm2(a::AbstractMatrix{Q}, origin::AbstractVector{Q0}) where {T,Q<:Quantity{T},Q0}
    @assert size(a, 1) == length(origin)
    Q === Q0 || (origin = convert(Vector{Q}, origin))
    u = unit(Q)^2
    Q² = Quantity{T,dimension(u),typeof(u)}
    dst = Vector{Q²}(undef, size(a, 2))
    return _colnorm2_unsafe!(dst, a, origin)
end

function colnorm2(a::AbstractMatrix{Union{Missing,Q}}, origin::AbstractVector{Q0}) where {T,Q<:Quantity{T},Q0}
    @assert size(a, 1) == length(origin)
    Q === Q0 || (origin = convert(Vector{Q}, origin))
    u = unit(Q)^2
    Q² = Quantity{T,dimension(u),typeof(u)}
    dst = Vector{Union{Missing,Q²}}(undef, size(a, 2))
    return _colnorm2_unsafe!(dst, a, origin)
end

function _colnorm_unsafe! end
function _colnorm2_unsafe! end

for (name, zeroT2) in zip(["_colnorm_unsafe", "_colnorm2_unsafe"], [:(zero(T)^2), :(zero(T))])
    quote
        function $(Symbol(name, "!"))(dst::AbstractVector{T}, a::AbstractMatrix) where {T}
            dims = size(a, 1)
            if dims == 2
                return $(Symbol(name, "_2D!"))(dst, a)
            elseif dims == 3
                return $(Symbol(name, "_3D!"))(dst, a)
            else
                return $(Symbol(name, "_ND!"))(dst, a; zeroT2=$zeroT2)
            end
        end

        function $(Symbol(name, "!"))(
            dst::AbstractVector{T},
            a::AbstractMatrix,
            origin::AbstractVector,
        ) where {T}
            dims = size(a, 1)
            if dims == 2
                return $(Symbol(name, "_2D!"))(dst, a, origin)
            elseif dims == 3
                return $(Symbol(name, "_3D!"))(dst, a, origin)
            else
                return $(Symbol(name, "_ND!"))(dst, a, origin; zeroT2=$zeroT2)
            end
        end
    end |> eval
end

function _colnorm2_unsafe_ND!(dst::AbstractVector{T}, a::AbstractMatrix; zeroT2=zero(T)^2) where {T}
    @inbounds @simd for i in axes(a, 2)
        dsti = zeroT2
        for j in axes(a, 1)
            dsti += a[j, i]^2
        end
        dst[i] = dsti
    end
    return dst
end

function _colnorm2_unsafe_2D!(dst::AbstractVector, a::AbstractMatrix)
    @inbounds @simd for i in axes(a, 2)
        dst[i] = a[1, i]^2 + a[2, i]^2
    end
    return dst
end

function _colnorm2_unsafe_3D!(dst::AbstractVector, a::AbstractMatrix)
    @inbounds @simd for i in axes(a, 2)
        dst[i] = a[1, i]^2 + a[2, i]^2 + a[3, i]^2
    end
    return dst
end

function _colnorm2_unsafe_ND!(
    dst::AbstractVector{T},
    a::AbstractMatrix,
    origin::AbstractVector;
    zeroT2=zero(T)^2,
) where {T}
    @inbounds @simd for i in axes(a, 2)
        dsti = zeroT2
        for j in axes(a, 1)
            dsti += (a[j, i] - origin[j])^2
        end
        dst[i] = dsti
    end
    return dst
end

function _colnorm2_unsafe_2D!(dst::AbstractVector, a::AbstractMatrix, origin::AbstractVector)
    @inbounds @simd for i in axes(a, 2)
        dst[i] = (a[1, i] - origin[1])^2 + (a[2, i] - origin[2])^2
    end
    return dst
end

function _colnorm2_unsafe_3D!(dst::AbstractVector, a::AbstractMatrix, origin::AbstractVector)
    @inbounds @simd for i in axes(a, 2)
        dst[i] = (a[1, i] - origin[1])^2 + (a[2, i] - origin[2])^2 + (a[3, i] - origin[3])^2
    end
    return dst
end

function _colnorm_unsafe_ND!(dst::AbstractVector{T}, a::AbstractMatrix; zeroT2=zero(T)^2) where {T}
    @inbounds @simd for i in axes(a, 2)
        dsti = zeroT2
        for j in axes(a, 1)
            dsti += a[j, i]^2
        end
        dst[i] = sqrt(dsti)
    end
    return dst
end

function _colnorm_unsafe_2D!(dst::AbstractVector, a::AbstractMatrix)
    @inbounds @simd for i in axes(a, 2)
        dst[i] = sqrt(a[1, i]^2 + a[2, i]^2)
    end
    return dst
end

function _colnorm_unsafe_3D!(dst::AbstractVector, a::AbstractMatrix)
    @inbounds @simd for i in axes(a, 2)
        dst[i] = sqrt(a[1, i]^2 + a[2, i]^2 + a[3, i]^2)
    end
    return dst
end

function _colnorm_unsafe_ND!(
    dst::AbstractVector{T},
    a::AbstractMatrix,
    origin::AbstractVector;
    zeroT2=zero(T)^2,
) where {T}
    @inbounds @simd for i in axes(a, 2)
        dsti = zeroT2
        for j in axes(a, 1)
            dsti += (a[j, i] - origin[j])^2
        end
        dst[i] = sqrt(dsti)
    end
    return dst
end

function _colnorm_unsafe_2D!(dst::AbstractVector, a::AbstractMatrix, origin::AbstractVector)
    @inbounds @simd for i in axes(a, 2)
        dst[i] = sqrt((a[1, i] - origin[1])^2 + (a[2, i] - origin[2])^2)
    end
    return dst
end

function _colnorm_unsafe_3D!(dst::AbstractVector, a::AbstractMatrix, origin::AbstractVector)
    @inbounds @simd for i in axes(a, 2)
        dst[i] = sqrt((a[1, i] - origin[1])^2 + (a[2, i] - origin[2])^2 + (a[3, i] - origin[3])^2)
    end
    return dst
end

"""
    coldot(a::AbstractMatrix, b::AbstractMatrix)

Returns the dot products of the matrix columns as a new matrix.

Both input matrices need to have the same dimensions ``d × N``, representing ``N`` ``d``-dimensional vectors.
"""
function coldot(a::AbstractMatrix, b::AbstractMatrix)
    @assert size(a) == size(b)
    return _coldot_unsafe(a, b)
end

function _coldot_unsafe(a::AbstractMatrix, b::AbstractMatrix)
    dst = similar(a, Any, size(a, 2))
    return _coldot_unsafe!(dst, a, b)
end

function _coldot_unsafe(a::AbstractMatrix{T1}, b::AbstractMatrix{T2}) where {T1<:Number,T2<:Number}
    T = typeof(zero(T1) * zero(T2))
    dst = similar(a, T, size(a, 2))
    return _coldot_unsafe!(dst, a, b)
end

function _coldot_unsafe(
    a::AbstractMatrix{Union{T1,Missing}},
    b::AbstractMatrix{T2},
) where {T1<:Number,T2<:Number}
    T = typeof(zero(T1) * zero(T2))
    dst = similar(a, Union{T,Missing}, size(a, 2))
    return _coldot_unsafe!(dst, a, b)
end

function _coldot_unsafe(
    a::AbstractMatrix{T1},
    b::AbstractMatrix{Union{T2,Missing}},
) where {T1<:Number,T2<:Number}
    T = typeof(zero(T1) * zero(T2))
    dst = similar(a, Union{T,Missing}, size(a, 2))
    return _coldot_unsafe!(dst, a, b)
end

function _coldot_unsafe(
    a::AbstractMatrix{Union{T1,Missing}},
    b::AbstractMatrix{Union{T2,Missing}},
) where {T1<:Number,T2<:Number}
    T = typeof(zero(T1) * zero(T2))
    dst = similar(a, Union{T,Missing}, size(a, 2))
    return _coldot_unsafe!(dst, a, b)
end

function _coldot_unsafe!(dst::AbstractVector, a::AbstractMatrix, b::AbstractMatrix)
    dims = size(a, 1)
    if dims == 2
        return _coldot_unsafe_2D!(dst, a, b)
    elseif dims == 3
        return _coldot_unsafe_3D!(dst, a, b)
    else
        return _coldot_unsafe_ND!(dst, a, b)
    end
end

function _coldot_unsafe_2D!(dst::AbstractVector, a::AbstractMatrix, b::AbstractMatrix)
    @inbounds @simd for i in axes(a, 2)
        dst[i] = a[1, i] * b[1, i] + a[2, i] * b[2, i]
    end
    return dst
end

function _coldot_unsafe_3D!(dst::AbstractVector, a::AbstractMatrix, b::AbstractMatrix)
    @inbounds @simd for i in axes(a, 2)
        dst[i] = a[1, i] * b[1, i] + a[2, i] * b[2, i] + a[3, i] * b[3, i]
    end
    return dst
end

function _coldot_unsafe_ND!(dst::AbstractVector{T}, a::AbstractMatrix, b::AbstractMatrix) where {T}
    @inbounds for i in axes(a, 2)
        dsti = a[1, i] * b[1, i]
        for j in Iterators.drop(axes(a, 1), 1)
            dsti += a[j, i] * b[j, i]
        end
        dst[i] = dsti
    end
    return dst
end

"""
    colcross(a::AbstractMatrix, b::AbstractMatrix)

Returns the cross products of the matrix columns as a new matrix.

Both input matrices need to have the same dimensions ``3 × N``.
"""
function colcross(a::AbstractMatrix, b::AbstractMatrix)
    @assert size(a) == size(b)
    @assert size(a, 1) == 3
    return _colcross_unsafe(a, b)
end

function _colcross_unsafe(a::AbstractMatrix, b::AbstractMatrix)
    dst = similar(a, Any)
    return _colcross_unsafe!(dst, a, b)
end

function _colcross_unsafe(a::AbstractMatrix{T1}, b::AbstractMatrix{T2}) where {T1<:Number,T2<:Number}
    T = typeof(zero(T1) * zero(T2))
    dst = similar(a, T)
    return _colcross_unsafe!(dst, a, b)
end

function _colcross_unsafe(
    a::AbstractMatrix{Union{T1,Missing}},
    b::AbstractMatrix{T2},
) where {T1<:Number,T2<:Number}
    T = typeof(zero(T1) * zero(T2))
    dst = similar(a, Union{T,Missing})
    return _colcross_unsafe!(dst, a, b)
end

function _colcross_unsafe(
    a::AbstractMatrix{T1},
    b::AbstractMatrix{Union{T2,Missing}},
) where {T1<:Number,T2<:Number}
    T = typeof(zero(T1) * zero(T2))
    dst = similar(a, Union{T,Missing})
    return _colcross_unsafe!(dst, a, b)
end

function _colcross_unsafe(
    a::AbstractMatrix{Union{T1,Missing}},
    b::AbstractMatrix{Union{T2,Missing}},
) where {T1<:Number,T2<:Number}
    T = typeof(zero(T1) * zero(T2))
    dst = similar(a, Union{T,Missing})
    return _colcross_unsafe!(dst, a, b)
end

function _colcross_unsafe!(dst::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix)
    @inbounds @simd for i in axes(dst, 2)
        dst[1, i] = a[2, i] * b[3, i] - a[3, i] * b[2, i]
        dst[2, i] = a[3, i] * b[1, i] - a[1, i] * b[3, i]
        dst[3, i] = a[1, i] * b[2, i] - a[2, i] * b[1, i]
    end
    return dst
end

"""
    angmom(pos::AbstractMatrix, vel::AbstractMatrix, mass; origin=nothing, velorigin=nothing)

Returns the angular momentum based on the positions, velocities, and masses.

Both input matrices need to have the same dimensions ``3 × N``.
The mass can be a vector or a scalar value (useful if all particles have the same mass).

The coordinate system's origin can be passed as a vector of length 3, as well as the velocity origin, which
determines the relative velocity at which the angular momentum should be determined. Typically, this will be the
velocity of the center of mass.
"""
function angmom(
    pos::AbstractMatrix,
    vel::AbstractMatrix,
    mass::AbstractVector;
    origin=nothing,
    velorigin=nothing,
)
    @assert size(pos, 2) == size(vel, 2) == length(mass)
    @assert size(pos, 1) == size(vel, 1) == 3
    return _angmom_unsafe(pos, vel, mass, origin, velorigin)
end

function angmom(pos::AbstractMatrix, vel::AbstractMatrix, mass::Number; origin=nothing, velorigin=nothing)
    @assert size(pos, 2) == size(vel, 2)
    @assert size(pos, 1) == size(vel, 1) == 3
    return _angmom_unsafe(pos, vel, mass, origin, velorigin)
end

function _angmom_unsafe(pos::AbstractMatrix, vel::AbstractMatrix, mass, origin, velorigin)
    dst = similar(pos, Any)
    return _angmom_unsafe!(dst, pos, vel, mass, origin, velorigin)
end

_help_convert_origin(::Type{T}, origin::AbstractVector) where {T} = convert(Vector{T}, origin)
_help_convert_origin(::Type{T}, origin::Nothing) where {T} = nothing

function _angmom_unsafe(
    pos::AbstractMatrix{T1},
    vel::AbstractMatrix{T2},
    mass,
    origin,
    velorigin,
) where {T1<:Number,T2<:Number}
    T = typeof(zero(T1) * zero(T2) * _zero_of(mass))
    dst = similar(pos, T)
    return _angmom_unsafe!(
        dst,
        pos,
        vel,
        mass,
        _help_convert_origin(T1, origin),
        _help_convert_origin(T2, velorigin),
    )
end

function _angmom_unsafe(
    pos::AbstractMatrix{Union{T1,Missing}},
    vel::AbstractMatrix{T2},
    mass,
    origin,
    velorigin,
) where {T1<:Number,T2<:Number}
    T = typeof(zero(T1) * zero(T2) * _zero_of(mass))
    dst = similar(pos, Union{T,Missing})
    return _angmom_unsafe!(dst, pos, vel, mass, origin, velorigin)
end

function _angmom_unsafe(
    pos::AbstractMatrix{T1},
    vel::AbstractMatrix{Union{T2,Missing}},
    mass,
    origin,
    velorigin,
) where {T1<:Number,T2<:Number}
    T = typeof(zero(T1) * zero(T2) * _zero_of(mass))
    dst = similar(pos, Union{T,Missing})
    return _angmom_unsafe!(dst, pos, vel, mass, origin, velorigin)
end

function _angmom_unsafe(
    pos::AbstractMatrix{Union{T1,Missing}},
    vel::AbstractMatrix{Union{T2,Missing}},
    mass,
    origin,
    velorigin,
) where {T1<:Number,T2<:Number}
    T = typeof(zero(T1) * zero(T2) * _zero_of(mass))
    dst = similar(pos, Union{T,Missing})
    return _angmom_unsafe!(dst, pos, vel, mass, origin, velorigin)
end

for (T1, T2) in zip(
    [:Nothing, :AbstractVector, :Nothing, :AbstractVector],
    [:Nothing, :Nothing, :AbstractVector, :AbstractVector],
)
    for T3 in [:Number, :AbstractVector]
        if T1 === :Nothing
            p1 = :(p[1, i])
            p2 = :(p[2, i])
            p3 = :(p[3, i])
        else
            p1 = :((p[1, i] - origin[1]))
            p2 = :((p[2, i] - origin[2]))
            p3 = :((p[3, i] - origin[3]))
        end
        if T2 === :Nothing
            v1 = :(v[1, i])
            v2 = :(v[2, i])
            v3 = :(v[3, i])
        else
            v1 = :((v[1, i] - velorigin[1]))
            v2 = :((v[2, i] - velorigin[2]))
            v3 = :((v[3, i] - velorigin[3]))
        end
        if T3 === :Number
            mi = :m
        else
            mi = :(m[i])
        end

        quote
            function _angmom_unsafe!(
                dst::AbstractMatrix,
                p::AbstractMatrix,
                v::AbstractMatrix,
                m::$T3,
                origin::$T1,
                velorigin::$T2,
            )
                @inbounds @simd for i in axes(dst, 2)
                    dst[1, i] = $mi * ($p2 * $v3 - $p3 * $v2)
                    dst[2, i] = $mi * ($p3 * $v1 - $p1 * $v3)
                    dst[3, i] = $mi * ($p1 * $v2 - $p2 * $v1)
                end
                return dst
            end
        end |> eval
    end
end

_zero_of(::T) where {T<:Number} = zero(T)
_zero_of(::AbstractArray{T}) where {T<:Number} = zero(T)

"""
    angmomtot(pos::AbstractMatrix, vel::AbstractMatrix, mass; origin=nothing, velorigin=nothing)

Returns the total angular momentum based on the positions, velocities, and masses.

Both input matrices need to have the same dimensions ``3 × N``.
The mass can be a vector or a scalar value (useful if all particles have the same mass).

The coordinate system's origin can be passed as a vector of length 3, as well as the velocity origin, which
determines the relative velocity at which the angular momentum should be determined. Typically, this will be the
velocity of the center of mass.

This method never allocates the full angular momentum matrix for all particles like [`angmom`](@ref) does.
Use this when the individual particle angular momenta are not needed.
"""
function angmomtot(
    pos::AbstractMatrix,
    vel::AbstractMatrix,
    mass::AbstractVector;
    origin=nothing,
    velorigin=nothing,
)
    @assert size(pos, 2) == size(vel, 2) == length(mass)
    @assert size(pos, 1) == size(vel, 1) == 3
    return _angmomtot_unsafe(pos, vel, mass, origin, velorigin)
end

function angmomtot(pos::AbstractMatrix, vel::AbstractMatrix, mass::Number; origin=nothing, velorigin=nothing)
    @assert size(pos, 2) == size(vel, 2)
    @assert size(pos, 1) == size(vel, 1) == 3
    return _angmomtot_unsafe(pos, vel, mass, origin, velorigin)
end

function _angmomtot_unsafe(pos::AbstractMatrix, vel::AbstractMatrix, mass, origin, velorigin)
    T, zeroT = _type_and_zero_of_summed_angmomtot(pos, vel, mass)
    dst = Vector{T}(undef, 3)
    _angmomtot_unsafe!(dst, pos, vel, mass, origin, velorigin; zeroT)
end

function _type_and_zero_of_summed_angmomtot(
    p::AbstractMatrix{T1},
    v::AbstractMatrix{T2},
    m::T3,
) where {T1,T2,T3<:Number}
    z = zero(T1) * zero(T2)
    return typeof(z * zero(T3)), z
end

function _type_and_zero_of_summed_angmomtot(
    p::AbstractMatrix{T1},
    v::AbstractMatrix{T2},
    m::AbstractVector{T3},
) where {T1,T2,T3}
    z = zero(T1) * zero(T2) * zero(T3)
    return typeof(z), z
end

"""
    angmomtot_stable(pos::AbstractMatrix, vel::AbstractMatrix, mass; origin=nothing, velorigin=nothing)

Returns the total angular momentum based on the positions, velocities, and masses, like [`angmomtot`](@ref).

Uses the stable summation algorithm `Base.sum` for the summation, but is almost three times slower than
[`angmomtot`](@ref).

See [`angmomtot`](@ref) for the available keyword arguments.
"""
function angmomtot_stable end

function angmomtot_stable(
    pos::AbstractMatrix,
    vel::AbstractMatrix,
    mass::AbstractVector;
    origin=nothing,
    velorigin=nothing,
)
    @assert size(pos, 2) == size(vel, 2) == length(mass)
    @assert size(pos, 1) == size(vel, 1) == 3
    return _angmomtot_stable_unsafe(pos, vel, mass, origin, velorigin)
end

function angmomtot_stable(
    pos::AbstractMatrix,
    vel::AbstractMatrix,
    mass::Number;
    origin=nothing,
    velorigin=nothing,
)
    @assert size(pos, 2) == size(vel, 2)
    @assert size(pos, 1) == size(vel, 1) == 3
    return _angmomtot_stable_unsafe(pos, vel, mass, origin, velorigin)
end

for (T1, T2) in zip(
    [:Nothing, :AbstractVector, :Nothing, :AbstractVector],
    [:Nothing, :Nothing, :AbstractVector, :AbstractVector],
)
    if T1 === :Nothing
        p1 = :(p[1, i])
        p2 = :(p[2, i])
        p3 = :(p[3, i])
    else
        p1 = :((p[1, i] - origin[1]))
        p2 = :((p[2, i] - origin[2]))
        p3 = :((p[3, i] - origin[3]))
    end
    if T2 === :Nothing
        v1 = :(v[1, i])
        v2 = :(v[2, i])
        v3 = :(v[3, i])
    else
        v1 = :((v[1, i] - velorigin[1]))
        v2 = :((v[2, i] - velorigin[2]))
        v3 = :((v[3, i] - velorigin[3]))
    end

    quote
        function _angmomtot_unsafe!(
            dst::AbstractVector,
            p::AbstractMatrix,
            v::AbstractMatrix,
            m::Number,
            origin::$T1,
            velorigin::$T2;
            zeroT,
        )
            d1 = d2 = d3 = zeroT
            @inbounds @simd for i in axes(p, 2)
                d1 += $p2 * $v3 - $p3 * $v2
                d2 += $p3 * $v1 - $p1 * $v3
                d3 += $p1 * $v2 - $p2 * $v1
            end
            dst .= m * d1, m * d2, m * d3
            return dst
        end

        function _angmomtot_unsafe!(
            dst::AbstractVector,
            p::AbstractMatrix,
            v::AbstractMatrix,
            m::AbstractVector,
            origin::$T1,
            velorigin::$T2;
            zeroT,
        )
            dst .= zeroT
            @inbounds @simd for i in axes(p, 2)
                dst[1] += m[i] * ($p2 * $v3 - $p3 * $v2)
                dst[2] += m[i] * ($p3 * $v1 - $p1 * $v3)
                dst[3] += m[i] * ($p1 * $v2 - $p2 * $v1)
            end
            return dst
        end

        function _angmomtot_stable_unsafe(
            p::AbstractMatrix,
            v::AbstractMatrix,
            m::Number,
            origin::$T1,
            velorigin::$T2,
        )
            return [
                m * sum($p2 * $v3 - $p3 * $v2 for i in axes(p, 2))
                m * sum($p3 * $v1 - $p1 * $v3 for i in axes(p, 2))
                m * sum($p1 * $v2 - $p2 * $v1 for i in axes(p, 2))
            ]
        end

        function _angmomtot_stable_unsafe(
            p::AbstractMatrix,
            v::AbstractMatrix,
            m::AbstractVector,
            origin::$T1,
            velorigin::$T2,
        )
            return [
                sum(m[i] * ($p2 * $v3 - $p3 * $v2) for i in axes(p, 2))
                sum(m[i] * ($p3 * $v1 - $p1 * $v3) for i in axes(p, 2))
                sum(m[i] * ($p1 * $v2 - $p2 * $v1) for i in axes(p, 2))
            ]
        end
    end |> eval
end


text:
geometric_mean(a::Number, b::Number) = √(a*b)

function find_under_and_over(a::Number, possible_values::Array)
    possible_series_sorted = sort(possible_values)
    i, j = 1, 2
        while !(possible_series_sorted[i] ≤ a ≤ possible_series_sorted[j])
            i+=1; j+=1
        end
    return possible_series_sorted[i], possible_series_sorted[j]
end

function norm_to_between_100_and_1000(val::Number; return_OOM=true)
    power_of_10 = 0.0
    while val > 1000
        val/=10; power_of_10 += 1
    end
    while val < 100
        val*=10; power_of_10 += -1
    end
    if return_OOM
        return val, power_of_10
    else
        return val
    end
end

ext: jl


text:
function pressureFunc( cpfile::String, Vinf::Float64, ρinf::Float64,
                       wrat::Float64, c::Float64 )

    data = readdlm( cpfile )

    s = data[:,1]
    inds = find( s[1:end-1] - s[2:end] .> 0.5 )
    ind = inds[1]

    xup = data[     1:ind,1 ]
    xlo = data[ ind+1:end,1 ]
    cpup = data[     1:ind, 5 ]
    cplo = data[ ind+1:end, 5 ]

    if xup[1] > xup[2]
        inds2 = find( xup .== minimum(xup) )
        ind2 = inds2[1]
        xlo  = vcat( xup[ind2:-1:2],  xlo )
        cplo = vcat( cpup[ind2:-1:2], cplo )
        deleteat!( xup,  1:(ind2-1) )
        deleteat!( cpup, 1:(ind2-1) )
    end

    xup = vcat( -1e-4, xup )
    xlo = vcat( -1e-4, xlo )

    cpup = vcat( cpup[1], cpup )
    cplo = vcat( cplo[1], cplo )

    spl_up = Dierckx.Spline1D( xup, cpup )
    spl_lo = Dierckx.Spline1D( xlo, cplo )

    pu( p::Matrix{Float64} ) = hcat( wrat * 0.5 * ρinf * Vinf^2 * spl_up( p[:,1]/c ), 0.0*p[:,2] )
    pl( p::Matrix{Float64} ) = hcat( wrat * 0.5 * ρinf * Vinf^2 * spl_lo( p[:,1]/c ), 0.0*p[:,2] )

    return (pu, pl)

end


text:
# Unit tests for functions in utility.jl

M_c = [complex(1,0)       complex(0,-1)
       complex(-0.5,0.5)  complex(0,0)]

M_r = [1.0   0.0   0.0  -1.0
       0.0   1.0   1.0  0.0
      -0.5  0.5    0.0  0.0
      -0.5  -0.5   0.0  0.0]

# Variable bounds to test auxiliary_variable_bounds function
l = [-2, -1, 0, 2]
u = [-1,  2, 3, 2]

@testset "Tests: complex to real matrix function" begin
    test_M_r = QCO.complex_to_real_gate(M_c)
    @test test_M_r == M_r
end

@testset "Tests: real to complex matrix function" begin
    test_M_c = QCO.real_to_complex_gate(M_r)
    @test test_M_c == M_c
end

@testset "Tests: auxiliary variable product bounds" begin
    m = Model()
    @variable(m, l[i] <= x[i=1:length(l)] <= u[i])
    # x1*x2
    @test QCO.auxiliary_variable_bounds([x[1], x[2]]) == (-4, 2)
    # x1*x3
    @test QCO.auxiliary_variable_bounds([x[1], x[3]]) == (-6, 0)
    # x1*x4
    @test QCO.auxiliary_variable_bounds([x[1], x[4]]) == (-4, -2)
    # x4*x4
    @test QCO.auxiliary_variable_bounds([x[4], x[4]]) == (4, 4)
    # x1*x2*x3
    @test QCO.auxiliary_variable_bounds([x[1], x[2], x[3]]) == (-12, 6)
    @test QCO.auxiliary_variable_bounds([x[3], x[1], x[2]]) == (-12, 6)
    @test QCO.auxiliary_variable_bounds([x[2], x[1], x[3]]) == (-12, 6)
    # x1*x2*x3*x4
    @test QCO.auxiliary_variable_bounds([x[1], x[2], x[3], x[4]]) == (-24, 12)
    @test QCO.auxiliary_variable_bounds([x[4], x[1], x[3], x[2]]) == (-24, 12)
    @test QCO.auxiliary_variable_bounds([x[4], x[3], x[2], x[1]]) == (-24, 12)
end

@testset "Tests: unique_matrices and unique_idx" begin

    D = zeros(3,3,4)
    D[:,:,1] = Matrix(LA.I,3,3)
    D[:,:,2] = rand(3,3)
    D[:,:,3] = Matrix(LA.I,3,3)
    D[:,:,4] = rand(3,3)

    D[1,1,3] = 1 + 2*tol_0
    D[3,3,3] = 1 - tol_0
    D[1,2,3] = 0 - tol_0
    D[1,3,3] = 0 + tol_0

    D[isapprox.(D, 0, atol=tol_0)] .= 0

    D_unique, D_unique_idx = QCO.unique_matrices(D)

    @test D_unique_idx == [1,2,4]
    @test D_unique[:,:,1] == D[:,:,1]
    @test D_unique[:,:,2] == D[:,:,2]
    @test D_unique[:,:,3] == D[:,:,4]

end

@testset "Tests: commuting matrices" begin
    params = Dict{String, Any}(
       "num_qubits" => 2,
       "maximum_depth" => 5,

       "elementary_gates" => ["H_1", "H_2", "CNot_1_2", "Identity"],
       "target_gate" => QCO.CNotRevGate()
       )

    data = QCO.get_data(params)
    C2,C2_Identity = QCO.get_commutative_gate_pairs(data["gates_dict"])
    @test length(C2) == 4
    @test length(C2_Identity) == 0
    @test isapprox(data["gates_real"][:,:,C2[1][1]] * data["gates_real"][:,:,C2[1][2]], data["gates_real"][:,:,C2[1][2]] * data["gates_real"][:,:,C2[1][1]], atol=tol_0)
end

@testset "Tests: kron_single_qubit_gate" begin
    I1 = QCO.kron_single_qubit_gate(4, QCO.IGate(1), "q1")
    I2 = QCO.kron_single_qubit_gate(4, QCO.IGate(1), "q2")
    I3 = QCO.kron_single_qubit_gate(4, QCO.IGate(1), "q3")
    I4 = QCO.kron_single_qubit_gate(4, QCO.IGate(1), "q4")
    @test I1 == I2 == I3 == I4
end

@testset "Tests: get_redundant_gate_product_pairs" begin

    params = Dict{String, Any}(
        "num_qubits" => 3,
        "maximum_depth" => 15,
        "elementary_gates" => ["T_1", "T_2", "T_3", "H_3", "CNot_1_2", "CNot_1_3", "CNot_2_3", "Tdagger_2", "Tdagger_3", "Identity"],
        "target_gate" => QCO.ToffoliGate()
        )
    data = QCO.get_data(params)
    redundant_pairs = QCO.get_redundant_gate_product_pairs(data["gates_dict"])
    @test length(redundant_pairs) == 0

end

@testset "Tests: get_idempotent_gates" begin

    params = Dict{String, Any}(
        "num_qubits" => 2,
        "maximum_depth" => 14,
        "elementary_gates" => ["Y_1", "Y_2", "Z_1", "Z_2", "T_2", "Tdagger_1", "Sdagger_1", "SX_1", "SXdagger_2", "CNot_2_1", "CNot_1_2", "Identity"],
        "target_gate" => -QCO.HCoinGate())

    data = QCO.get_data(params)
    idempotent_gates = QCO.get_idempotent_gates(data["gates_dict"])
    @test idempotent_gates[1] == 6
    @test idempotent_gates[2] == 7
    @test "Tdagger_1" in data["gates_dict"]["$(idempotent_gates[1])"]["type"]
    @test "Sdagger_1" in data["gates_dict"]["$(idempotent_gates[2])"]["type"]
    @test "Z_1" in data["gates_dict"]["3"]["type"]

    M1 = data["gates_dict"]["$(idempotent_gates[1])"]["matrix"]
    M2 = data["gates_dict"]["$(idempotent_gates[2])"]["matrix"]

    @test isapprox(M1^2, data["gates_dict"]["7"]["matrix"], atol = tol_0)
    @test isapprox(M2^2, data["gates_dict"]["3"]["matrix"], atol = tol_0)
end

@testset "Tests: _parse_gate_string" begin
    v1 = QCO._parse_gate_string("RX_1", qubits=true)
    @test v1[1] == 1

    v2 = QCO._parse_gate_string("CU3_51", qubits=true)
    @test v2[1] == 51

    v3 = QCO._parse_gate_string("CRZ_9_1", qubits=true)
    @test length(v3) == 2
    @test v3[1] == 9
    @test v3[2] == 1
end

@testset "Tests: is_gate_real" begin
    M_test = [0 1E-7im -1E-7im; 1E-10im 0 0; -1E-8im -1E-10im 0]
    @test QCO.is_gate_real(M_test)

    M_test = [0 1E-7im -1E-7im; 1E-10im 0 0; -1E-8im -1E-5im 0]
    @test !QCO.is_gate_real(M_test)
end

@testset "Tests: _get_constraint_slope_intercept" begin
    v1 = (0.0, 0.0)
    v2 = (1.0, 1.0)
    m,c = QCO._get_constraint_slope_intercept(v1, v2)
    @test m == 1
    @test c == 0

    m,c = QCO._get_constraint_slope_intercept(v2, v1)
    @test m == 1
    @test c == 0

    v2 = (0, 0.5)
    m,c = QCO._get_constraint_slope_intercept(v1, v2)
    @test !isfinite(m)
    @test !isfinite(c)

    v1 = (1.0, 2.0)
    v2 = (-3.0, 4.0)
    m,c = QCO._get_constraint_slope_intercept(v1, v2)
    @test m == -0.5
    @test c == 2.5

    QCO._get_constraint_slope_intercept(v1, v1)
end


text:
# Mini Hallway problem defined in http://cs.brown.edu/research/ai/pomdp/examples/mini-hall2.POMDP.
# Original idea published in Littman, Cassandra and Kaelbling's ML-95 paper.

# Basic parameters are:
# discount: 0.950000
# values: reward
# states: 13
# actions: 3
# observations: 9
# The rest is available at link

struct MiniHallway <: POMDP{Int, Int, Int}
    T::Array{Deterministic, 1}
end

function MiniHallway()
    T = Array{Deterministic, 1}(undef, 13)

    # Transitions for action 1 (and all actions in state 13) as I did not find a function for it
    T[1] = Deterministic(1); T[2] = Deterministic(2); T[3] = Deterministic(7);
    T[4] = Deterministic(4); T[5] = Deterministic(1); T[6] = Deterministic(10);
    T[8] = Deterministic(8); T[7] = Deterministic(1); T[9] = Deterministic(9);
    T[10] = Deterministic(10); T[11] = Deterministic(13); T[12] = Deterministic(8);
    T[13] = Deterministic(13)

    return MiniHallway(T)
end

##################
# mdps interface #
##################
POMDPs.states(m::MiniHallway) = 1:13
POMDPs.stateindex(m::MiniHallway, ss::Int)::Int = ss
POMDPs.isterminal(m::MiniHallway, ss::Int)::Bool = ss == 13

function POMDPs.transition(m::MiniHallway, ss::Int, a::Int)
    if a == 1 || ss == 13
        return m.T[ss]
    else
        if a == 2
            return ss % 4 == 0 ? Deterministic(ss - 3) : Deterministic(ss + 1)
        elseif a == 3
            return (ss - 1) % 4 == 0 ? Deterministic(ss + 3) : Deterministic(ss - 1)
        end
    end
end

POMDPs.actions(m::MiniHallway) = 1:3
POMDPs.actionindex(m::MiniHallway, a::Int)::Int = a

POMDPs.reward(m::MiniHallway, ss::Int, a::Int, sp::Int) = float(ss != sp && sp == 13)
POMDPs.reward(m::MiniHallway, ss::Int, a::Int) = mean_reward(m, ss, a)
POMDPs.discount(m::MiniHallway)::Float64 = 0.95

####################
# pomdps interface #
####################
POMDPs.initialstate(m::MiniHallway) = SparseCat(1:13, append!(fill(1/12, 12), 0.))
POMDPs.observations(m::MiniHallway) = 1:9

function POMDPs.observation(m::MiniHallway, a::Int, sp::Int)::Deterministic
    if sp <= 8
        return Deterministic(sp)
    elseif sp <= 10
        return Deterministic(sp - 2)
    elseif sp <= 12
        return Deterministic(sp - 6)
    else
        return Deterministic(9)
    end

    return m.O[sp]
end
POMDPs.observation(m::MiniHallway, s::Int, a::Int, sp::Int)::Deterministic = observation(m, a, sp)
POMDPs.obsindex(m::MiniHallway, o::Int)::Int = o

Base.broadcastable(m::MiniHallway) = Ref(m)


text:
using Comrade, ComradeGalactic, GalacticOptimJL, Distributions
using Test

load_ehtim()
include(joinpath(@__DIR__, "../../../test/test_util.jl"))

@testset "ComradeGalactic.jl" begin
    m, vis, amp, lcamp, cphase = load_data()
    prior = test_prior()
    lklhd = RadioLikelihood(lcamp, cphase)
    post = Posterior(lklhd, prior, test_model)

    tpost = ascube(post)
    f = OptimizationFunction(tpost, GalacticOptim.AutoForwardDiff{4}())
    x0 = [ 0.1,
           0.4,
           0.5,
           0.1,
           0.3,
           0.1,
           0.4,
           0.3,
           0.8,
           0.8]
    prob = GalacticOptim.OptimizationProblem(f, x0, nothing, lb=fill(0.001, 10), ub=fill(0.999,10))
    sol = solve(prob, LBFGS(); maxiters=100_000)

    xopt = transform(tpost, sol)
    @test isapprox(xopt.f1/xopt.f2, 2.0, atol=1e-3)
    @test isapprox(xopt.σ1*2*sqrt(2*log(2)), μas2rad(40.0), rtol=1e-3)
    @test isapprox(xopt.σ1*xopt.τ1*2*sqrt(2*log(2)), μas2rad(20.0), rtol=1e-3)
    @test isapprox(-xopt.ξ1, π/3, atol=1e-3)
    @test isapprox(xopt.σ2*2*sqrt(2*log(2)), μas2rad(20.0), atol=1e-3)
    @test isapprox(xopt.σ2*xopt.τ2*2*sqrt(2*log(2)), μas2rad(10.0), rtol=1e-3)
    @test isapprox(-xopt.ξ2, π/6, atol=1e-3)
    @test isapprox(xopt.x, μas2rad(30.0), rtol=1e-3)
    @test isapprox(xopt.y, μas2rad(30.0), rtol=1e-3)

end


text:
using Base.Test
using RobustStats

x=[1.672064, 0.7876588, 0.317322, 0.9721646, 0.4004206, 1.665123, 3.059971, 0.09459603, 1.27424, 3.522148,
   0.8211308, 1.328767, 2.825956, 0.1102891, 0.06314285, 2.59152, 8.624108, 0.6516885, 5.770285, 0.5154299]
y = sort(x)

test_approx(a,b, atol) = @test isapprox(a,b, atol=atol)
# Tests are based on examples from README

@test mean(x) ≈ 1.853401259
@test tmean(x) ≈ 1.2921802666666669
@test tmean(x, tr=0) ≈ 1.853401259
@test tmean(x, tr=0.3) ≈ 1.1466045875000002
@test tmean(x, tr=0.5) ≈ 1.1232023

# Winsorized data
winval(x)      # just verify that it doesn't crash
@test winmean(x) ≈ 1.4205834800000001
@test winvar(x) ≈ 0.998659015947531
@test wincov(x, x) ≈ 0.998659015947531
@test wincov(x, x.^2) ≈ 3.2819238397424004
@test RobustStats.winstd(x) ≈ 0.9993292830431474
@test trimse(x) ≈ 0.3724280347984342
out = trimci(x)
@test out.estimate ≈ 1.2921802666666669
@test out.p ≈ 0.005243565819244678

q1,q3 = idealf(x)
@assert length(x) == 20
@test q1 ≈ y[5]*7/12+y[6]*5/12
@test q3 ≈ y[15]*5/12+y[16]*7/12

@test pbvar(x) ≈ 2.0009575278957623
@test bivar(x) ≈ 1.5885279811329132

@test tauloc(x) ≈ 1.26966525675108
@test tauvar(x) ≈ 1.53008203090696

ob = outbox(x)
@test ob.nout == 1
@test ob.outid == [17]
@test ob.outval[1] ≈ 8.624108

@test msmedse(x) ≈ 0.4708261134886094
bci = binomci(2, 10)
@test bci.p_hat == 2/10
@test bci.confint[1] ≈ 0.02735907527075918
@test bci.confint[2] ≈ 0.5562210572204956
bci = acbinomci(2, 10)
@test bci.p_hat == 2/10
@test bci.confint[1] ≈ 0.04588727048388419
@test bci.confint[2] ≈ 0.5206324094338493

@test onestep(x) ≈ 1.3058109021286803
@test mom(x) ≈ 1.2596462322222222
mci = momci(x, seed=2, nboot=2000, nullvalue=0.6)
@test mci.ci[1] ≈ 0.5042229275
@test mci.ci[2] ≈ 2.120978892222222
@test mci.p ≈ 0.11099999999999999

ci = bootstrapci(x, est=onestep, nullvalue=0.6)
@test ci.estimate ≈ 1.3058109021286803
@test ci.ci[1] ≈ 0.6877225172104009
@test ci.ci[2] ≈ 2.259070870017077
@test ci.p ≈ 0.026000000000000023
se = bootstrapse(x, est=onestep)
@test se ≈ 0.41956761772722817

s = sint(x)
@test s.ci[1] ≈ 0.5474825849287168
@test s.ci[2] ≈ 2.3752324887983707
s = sint(x, 0.6)
@test s.ci[1] ≈ 0.5474825849287168
@test s.ci[2] ≈ 2.3752324887983707
@test s.p ≈ 0.07186070461307993

tm =trimpb(x, nullvalue=0.75)
@test tm.estimate ≈ 1.2921802666666669
@test tm.ci[1] ≈ 0.6905386583333333
@test tm.ci[2] ≈ 2.1963814608333334
@test tm.p ≈ 0.08600000000000008

cv = pcorb(x, x.^5)
@test cv.estimate ≈ 0.8026388958896737
@test cv.ci[1] ≈ 0.6836999351727973
@test cv.ci[2] ≈ 0.9634782953906714

srand(3)
y2 = randn(20)+3;
res = yuend(x, y2)
test_approx(res.estimate, -1.547776, 1e-6)
test_approx(res.se, 0.460304, 1e-6)
test_approx(res.p, 0.006336, 1e-4)

# Basic tests of the weighted means
a = collect(-2:2)
test_approx(bisquareWM(a,3,.1,1e-5), 0.0, 1e-4)
test_approx(huberWM(a,3,.1,1e-5), 0.0, 1e-4)
push!(a, 97)
test_approx(bisquareWM(a,3,.1,1e-5), 0.0, 1e-4)
push!(a, 98)
test_approx(bisquareWM(a,3,.1,1e-5), 0.0, 1e-4)
push!(a, 99)
test_approx(bisquareWM(a,3,.1,1e-5), 0.0, 1e-4)
push!(a, 98)
test_approx(bisquareWM(a,3,.1,1e-5), 0.0, 1e-4)
append!(a, [98,98,98])
test_approx(bisquareWM(a,3,97,1e-4), 98.0, 1e-3)

@test_throws ErrorException bisquareWM(collect(0:9), 0.1)

# Test 10 sets of 100 N(0,1) random numbers
# Add up to 8 values of much larger tailiness
# Expected mean is within [-.2, +.2] (2-sigma), so test for
# being in the [-1, +1] range, even with added large values.
for i = 1:10
    r = randn(100)
    for j = 1:5
        b = bisquareWM(r, 4, 0, 0.01)
        h = huberWM(r, 1.3, 0, 0.01)
        test_approx(b, 0, 1)
        test_approx(h, 0, 1)
        test_approx(b, h, 1.5)
        append!(r, randn(2)*100)
    end
end

# Trimean: exact quartiles
@test trimean(0:12) ≈ 6.0
# Trimean: exact median, inexact 1st and 3rd quartiles
@test trimean(0:10) ≈ 5.0
# Trimean: inexact quartiles
@test trimean(0:9) ≈ 4.5
@test trimean(x) ≈ 1.35575501875

# Shortest half-range
a = [10,3,5,6,6.5,7,8,0,13]  # Odd # of values. 9 values, with "half" containing 5
@test shorthrange(a) ≈ 3
@test a[1] ≈ 10 # Be sure it didn't rearrange a
@test shorthrange!(a) ≈ 3
@test a[1] ≈ 0  # Should have sorted a
a = [0,4.5,6,7,8,10,14,99]  # Even # of values. 8 values, with each "half" containing 5
@test shorthrange(a) ≈ 10-4.5

# Weighted high median
whm = RobustStats._weightedhighmedian
whm! = RobustStats._weightedhighmedian!

# Check that whm and whm! give correct and equal answers.
function _verify_whm{T<:Real,U<:Integer}(a::Vector{T}, wts::Vector{U})
	answer = whm(a,wts)
	wlow = whigh = wexact =  0
	for i=1:length(a)
		if a[i] < answer
			wlow += wts[i]
		elseif a[i] > answer
			whigh += wts[i]
		else
			wexact += wts[i]
		end
	end
	wtotal = wlow + wexact + whigh
	@test 2*wlow <= wtotal && 2*whigh < wtotal

	answermangled = whm!(a,wts)
	@test answer == answermangled
end

_verify_whm(collect(1:5), [1,1,1,1,1])
_verify_whm(collect(1:5), [1,2,3,4,5])
_verify_whm(collect(1:5), [1,2,3,4,9])
_verify_whm(collect(1:5), [1,2,3,4,10])
_verify_whm(collect(1:5), [2,1,1,1,1])
_verify_whm(collect(1:5), [1,1,1,2,1])

_verify_whm(collect(5:-1:1), [1,1,1,1,1])
_verify_whm(collect(5:-1:1), [1,2,3,4,5])
_verify_whm(collect(5:-1:1), [1,2,3,4,9])
_verify_whm(collect(5:-1:1), [1,2,3,4,10])
_verify_whm(collect(5:-1:1), [1,2,3,4,11])
_verify_whm(collect(5:-1:1), [2,1,1,1,1])
_verify_whm(collect(5:-1:1), [1,1,1,2,1])
_verify_whm([1,4,2,5,3,6], [1,4,2,5,3,6])
_verify_whm([1,4,2,5,3,6], [1,4,2,5,3,5])
_verify_whm([1,4,2,5,3,6], [1,4,2,5,3,4])

datalengths = [10,11,1000,1111]
wttypes = [Int8, Int16, Int32, Int64, UInt8, UInt16, UInt32, UInt64]
for N in datalengths
	for wttype in wttypes
		a = randn(N)
        w = rand(wttype, N)
		# Careful! Can't have sum of weights overflow, and can't have negative weights
		for k=1:N
			if w[k]<0
				w[k] = 1
			end
			w[k] %= 8192
		end
		result = whm(a,w)
		_verify_whm(a, w)
	end
end
@test_throws ArgumentError whm(collect(1:5),collect(1:4))

@test RobustStats._slow_scaleQ([1,2,3,4,5,10]) ≈ scaleQ!([1,2,3,4,5,10])
@test RobustStats._slow_scaleS!([1,2,3,4,5,10]) ≈ scaleS!([1,2,3,4,5,10])
@test RobustStats._slow_scaleQ([1,2,3,4,5,10.5]) ≈ scaleQ!([1,2,3,4,5,10.5])
@test RobustStats._slow_scaleS!([1,2,3,4,5,10.5]) ≈ scaleS!([1,2,3,4,5,10.5])

NTESTS = 10
for N in [2,3,4,5,6,7,8,9,10,11,12,15,20,25,50,100,151,200,225,250,299,350,399,500]
    for _=1:NTESTS
        a = randn(N)
        Q = RobustStats._slow_scaleQ(a)
        S = RobustStats._slow_scaleS!(a)
        @test Q ≈ scaleQ(a)
        @test Q ≈ scaleQ!(copy(a))
        @test S ≈ scaleS(a)
        @test S ≈ scaleS!(a)
        @test S ≈ scaleS!(a) # scaleS! sorts the array only, so result is unchanged
        @test Q ≈ scaleQ!(a)
        @test Q ≈ scaleQ!(a) # scaleQ! sorts the array only, so result is unchanged
    end
    # println("Success on $NTESTS scaleQ and scaleS tests with size $(N)")
end


text:
using LightGraphs
using DataStructures

struct SpanningTreeInstance{T} <: CombinatorialInstance
  graph::AbstractGraph{T}
  rewards::Dict{Edge{T}, Float64}
end

graph(i::SpanningTreeInstance{T}) where T = i.graph
# rewards(i::SpanningTreeInstance{T}) where T = i.rewards
function reward(i::SpanningTreeInstance{T}, e::Edge{T}) where T
  if e in keys(i.rewards)
    return i.rewards[e]
  end
  return i.rewards[reverse(e)]
end

struct SpanningTreeSolution{T} <: CombinatorialSolution
  instance::SpanningTreeInstance{T}
  tree::Vector{Edge{T}}
end


text:
"""
    aa2int(char)

Amino-acid one-letter code to integer (1 to 21, with gap = 21).
"""
function aa2int(c::Union{Char,UInt8})
    # from: https://github.com/carlobaldassi/GaussDCA.jl
    alphabet = Int8.((1,21, 2, 3, 4, 5, 6, 7, 8,21, 9,10,11,12,21,13,14,15,16,17,21,18,19,21,20))
              # A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R, S, T, U, V, W, X, Y
    i = UInt8(c) - 0x40
    1 <= i <= 25 && return alphabet[i]
    return Int8(21)
end

"""
    aa2int(seq::String)

Convert a string amino-acid sequence to integers.
"""
aa2int(s::String) = aa2int.(collect(s))

"""
    aa2int(seq::String[])

Given a vector of string amino-acid sequences, converts it
to an integer matrix.
"""
aa2int(ss::AbstractVector{String}) = hcat(aa2int.(ss)...)


text:
"""
    ZAFCorrection

Pulls together the ZA, F and coating corrections into a single structure.
"""
struct ZAFCorrection
    za::MatrixCorrection
    f::FluorescenceCorrection
    coating::CoatingCorrection
    """
        ZAFCorrection(za::MatrixCorrection, f::FluorescenceCorrection, coating::CoatingCorrection)
    """
    ZAFCorrection(za::MatrixCorrection, f::FluorescenceCorrection, coating::CoatingCorrection) = new(za, f, coating)
end

NeXLCore.material(zaf::ZAFCorrection) = material(zaf.za)

"""
    Z(unk::ZAFCorrection, std::ZAFCorrection)

Computes the atomic number correction.
"""
Z(unk::ZAFCorrection, std::ZAFCorrection) = Z(unk.za, std.za)

"""
    A(unk::ZAFCorrection, std::ZAFCorrection, cxr::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat)

Computes the absorption correction.
"""
A(unk::ZAFCorrection, std::ZAFCorrection, cxr::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat) =
    A(unk.za, std.za, cxr, θunk, θstd)

"""
    ZA(unk::ZAFCorrection, std::ZAFCorrection, cxr::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat)

Computes the combined atomic number and fluorescence correction.
"""
ZA(unk::ZAFCorrection, std::ZAFCorrection, cxr::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat) =
    ZA(unk.za, std.za, cxr, θunk, θstd)

"""
    coating(unk::ZAFCorrection, std::ZAFCorrection, cxr::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat)

Computes the coating correction.
"""
coating(unk::ZAFCorrection, std::ZAFCorrection, cxr::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat) =
    transmission(unk.coating, cxr, θunk) / transmission(std.coating, cxr, θstd)

"""
    generation(unk::ZAFCorrection, std::ZAFCorrection, ass::AtomicSubShell)

Computes a correction factor for differences X-ray generation due to differences in beam energy.
"""
generation(unk::ZAFCorrection, std::ZAFCorrection, ass::AtomicSubShell) =
    ionizationcrosssection(ass, beamEnergy(unk)) / ionizationcrosssection(ass, beamEnergy(std))

"""
    F(unk::ZAFCorrection, std::ZAFCorrection, cxr::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat)

Computes the secondary fluorescence correction.
"""
F(unk::ZAFCorrection, std::ZAFCorrection, cxr::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat) =
    F(unk.f, cxr, θunk) / F(std.f, cxr, θstd)

"""
    ZAFc(unk::ZAFCorrection, std::ZAFCorrection, cxr::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat)

Computes the combined correction for atomic number, absorption, secondary fluorescence and generation.
"""
ZAFc(unk::ZAFCorrection, std::ZAFCorrection, cxr::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat) =
    ZA(unk, std, cxr, θunk, θstd) * F(unk, std, cxr, θunk, θstd) * coating(unk, std, cxr, θunk, θstd)

"""
    k(unk::MultiZAF, std::MultiZAF, θunk::AbstractFloat, θstd::AbstractFloat)

The computed k-ratio for the unknown relative to standard.
"""
function k(unk::ZAFCorrection, std::ZAFCorrection, cxr::CharXRay, θunk::AbstractFloat, θstd::AbstractFloat)
    elm = element(cxr)
    return (nonneg(material(unk), elm) / nonneg(material(std), elm)) * ZAFc(unk, std, cxr, θunk, θstd)
end

beamEnergy(zaf::ZAFCorrection) = beamEnergy(zaf.za)

Base.show(io::IO, cc::ZAFCorrection) = print(io, "ZAF[", cc.za, ", ", cc.f, ", ", cc.coating, "]")

"""
    NeXLUncertainties.asa(::Type{DataFrame}, unk::ZAFCorrection, std::ZAFCorrection, trans::AbstractVector{Transition},
    θunk::AbstractFloat, θstd::AbstractFloat)::DataFrame

Tabulate a matrix correction relative to the specified unknown and standard for the iterable of Transition, trans.
"""
function NeXLUncertainties.asa(#
    ::Type{DataFrame},
    unk::ZAFCorrection,
    std::ZAFCorrection,
    trans::NTuple{N,Transition},
    θunk::AbstractFloat,
    θstd::AbstractFloat,
)::DataFrame where {N}
    @assert isequal(atomicsubshell(unk.za), atomicsubshell(std.za))
    "The atomic sub-shell for the standard and unknown don't match."
    cxrs = characteristic(
        element(atomicsubshell(unk.za)),
        trans,
        1.0e-9,
        0.999 * min(beamEnergy(unk.za), beamEnergy(std.za)),
    )
    stds, celm, stdE0, unks, unkE0, xray = String[], Float64[], Float64[], String[], Float64[], CharXRay[]
    z, a, f, c, zaf, k, unkToa, stdToa = #
        Float64[], Float64[], Float64[], Float64[], Float64[], Float64[], Float64[], Float64[]
    for cxr in cxrs
        if isequal(inner(cxr), atomicsubshell(std.za))
            elm = element(cxr)
            push!(unks, name(material(unk.za)))
            push!(celm, material(unk.za)[elm])
            push!(unkE0, beamEnergy(unk.za))
            push!(unkToa, θunk)
            push!(stds, name(material(std.za)))
            push!(stdE0, beamEnergy(std.za))
            push!(stdToa, θstd)
            push!(xray, cxr)
            push!(z, Z(unk, std))
            push!(a, A(unk, std, cxr, θunk, θstd))
            push!(f, F(unk, std, cxr, θunk, θstd))
            push!(c, coating(unk, std, cxr, θunk, θstd))
            tot = ZAFc(unk, std, cxr, θunk, θstd)
            push!(zaf, tot)
            push!(k, tot * material(unk.za)[elm] / material(std.za)[elm])
        end
    end
    return DataFrame(
        Unknown = unks,
        E0unk = unkE0,
        TOAunk = unkToa,
        Standard = stds,
        E0std = stdE0,
        TOAstd = stdToa,
        Xray = xray,
        Z = z,
        A = a,
        F = f,
        c = c,
        ZAF = zaf,
        C = celm,
        k = k,
    )
end

NeXLUncertainties.asa( #
    ::Type{DataFrame},
    unk::ZAFCorrection,
    std::ZAFCorrection,
    θunk::AbstractFloat,
    θstd::AbstractFloat,
)::DataFrame = asa(DataFrame, unk, std, alltransitions, θunk, θstd)

function NeXLUncertainties.asa( #
    ::Type{DataFrame},
    zafs::Dict{ZAFCorrection,ZAFCorrection},
    θunk::AbstractFloat,
    θstd::AbstractFloat,
)::DataFrame
    df = DataFrame()
    for (unk, std) in zafs
        append!(df, asa(DataFrame, unk, std, θunk, θstd))
    end
    return df
end

"""
    zafcorrection(
      mctype::Type{<:MatrixCorrection},
      fctype::Type{<:FluorescenceCorrection},
      cctype::Type{<:CoatingCorrection},
      mat::Material,
      ashell::AtomicSubShell,
      e0::Real,
      coating::Union{Film,AbstractVector{Film},Missing}
    )

Constructs an ZAFCorrection object using the mctype correction model with
the fluorescence model for the specified parameters.
"""
function zafcorrection(
    mctype::Type{<:MatrixCorrection},
    fctype::Type{<:FluorescenceCorrection},
    cctype::Type{<:CoatingCorrection},
    mat::Material,
    ashell::AtomicSubShell,
    e0::Real,
    coating::Union{Film,AbstractVector{Film},Missing},
)
    norm = asnormalized(mat)  # Ensures convergence of the interation algorithms...
    return ZAFCorrection(
        matrixcorrection(mctype, norm, ashell, e0),
        fluorescencecorrection(fctype, norm, ashell, e0),
        coatingcorrection(cctype, coating),
    )
end

"""
    zafcorrection(
       mctype::Type{<:MatrixCorrection},
       fctype::Type{<:FluorescenceCorrection},
       cctype::Type{<:CoatingCorrection},
       unk::Material,
       std::Material,
       ashell::AtomicSubShell,
       e0::Real;
       unkCoating::Union{Film,AbstractVector{Film},Missing} = missing,
       stdCoating::Union{Film,AbstractVector{Film},Missing} = missing,
    )

Creates a matched pair of ZAFCorrection objects using the matrix correction algorithm
for the specified unknown and standard.
"""
function zafcorrection(
    mctype::Type{<:MatrixCorrection},
    fctype::Type{<:FluorescenceCorrection},
    cctype::Type{<:CoatingCorrection},
    unk::Material,
    std::Material,
    ashell::AtomicSubShell,
    e0::Real;
    unkCoating::Union{Film,AbstractVector{Film},Missing} = missing,
    stdCoating::Union{Film,AbstractVector{Film},Missing} = missing,
)
    @assert energy(ashell) < e0 "The shell energy must be less than the beam energy."
    return (
        zafcorrection(mctype, fctype, cctype, unk, ashell, e0, unkCoating),
        zafcorrection(mctype, fctype, cctype, std, ashell, e0, stdCoating),
    )
end

function NeXLUncertainties.asa( #
    ::Type{DataFrame},
    unk::Material,
    stds::Dict{Element,Material},
    e0::Real,
    θunk::AbstractFloat,
    θstd::AbstractFloat;
    mctype::Type{<:MatrixCorrection} = XPP,
    fctype::Type{<:FluorescenceCorrection} = ReedFluorescence,
    cctype::Type{<:CoatingCorrection} = Coating,
)
    df = DataFrame()
    for (elm, std) in stds
        for ashell in atomicsubshells(elm, e0)
            append!(
                df,
                asa( #
                    DataFrame,
                    zafcorrection(mctype, fctype, cctype, unk, std, ashell, e0)...,
                    alltransitions,
                    θunk,
                    θstd,
                ),
            )
        end
    end
    return df
end

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
#! /usr/bin/julia

# Rosetta Code, Sorting algorithms/Selection sort

function selectionsort!{T<:Real}(a::Array{T,1})
    len = length(a)
    if len < 2
        return nothing
    end
    for i in 1:len-1
        (lmin, j) = findmin(a[i+1:end])
        if lmin < a[i]
            a[i+j] = a[i]
            a[i] = lmin
        end
    end
    return nothing
end

a = [rand(-100:100) for i in 1:20]
println("Before Sort:")
println(a)
selectionsort!(a)
println("\nAfter Sort:")
println(a)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
## Expression parsing (not part of MOI API)

function _parse(::Optimizer, expr)
    error("_parse not yet implemented for type: $(typeof(expr))")
    return
end

function _parse(model::Optimizer, f::MOI.VariableIndex)
    # A Concert Variable is already an expression for CPLEX.
    return _info(model, f).variable
end

function _parse(model::Optimizer, f::MOI.VectorOfVariables)
    return [_parse(model, v) for v in f.variables]
end

function _parse(
    model::Optimizer,
    terms::Vector{MOI.ScalarAffineTerm{T}},
) where {T <: Integer}
    cp = model.inner
    coeffs = Int32[t.coefficient for t in terms] # Type is forced by JNI.
    vars = IloIntVar[_parse(model, t.variable) for t in terms]
    return cpo_java_scalprod(cp, coeffs, vars)
end

function _parse(
    model::Optimizer,
    terms::Vector{MOI.ScalarAffineTerm{T}},
) where {T <: Real}
    cp = model.inner
    coeffs = T[t.coefficient for t in terms]
    vars = NumVar[_parse(model, t.variable) for t in terms]
    return cpo_java_scalprod(cp, coeffs, vars)
end

function _parse(
    model::Optimizer,
    terms::Vector{MOI.ScalarQuadraticTerm{T}},
) where {T <: Real}
    # IloCP.scalprod only works for unweighted vectors of variables, which is not always the case here.
    cp = model.inner

    function prod_vars(term::MOI.ScalarQuadraticTerm{T})
        return cpo_java_prod(
            cp,
            _parse(model, term.variable_1),
            _parse(model, term.variable_2),
        )
    end
    function prod(term::MOI.ScalarQuadraticTerm{T})
        return cpo_java_prod(cp, term.coefficient, prod_vars(term))
    end

    e = prod(terms[1])
    for t in terms[2:end]
        e = cpo_java_sum(cp, e, prod(t))
    end
    return e
end

function _parse(
    model::Optimizer,
    f::MOI.ScalarAffineFunction{T},
) where {T <: Real}
    f = MOI.Utilities.canonical(f)
    e = _parse(model, f.terms)
    if !iszero(f.constant)
        cp = model.inner
        e = cpo_java_sum(cp, e, cpo_java_constant(cp, f.constant))
    end
    return e
end

function _parse(
    model::Optimizer,
    f::MOI.VectorAffineFunction{T},
) where {T <: Real}
    f = MOI.Utilities.canonical(f)
    cp = model.inner
    es = [
        cpo_java_constant(cp, f.constants[idx]) for
        idx in 1:MOI.output_dimension(f)
    ]
    for i in 1:MOI.output_dimension(f)
        if i <= length(f.terms)
            idx = f.terms[i].output_index
            term = f.terms[i].scalar_term
            es[idx] = cpo_java_sum(
                cp,
                [
                    es[idx],
                    cpo_java_prod(
                        cp,
                        IntExpr[cpo_java_constant(cp, term.coefficient)],
                        IntExpr[_parse(model, term.variable)],
                    ),
                ],
            )
        end
    end
    return es
end

function _parse(
    model::Optimizer,
    f::MOI.ScalarQuadraticFunction{T},
) where {T <: Real}
    f = MOI.Utilities.canonical(f)
    cp = model.inner
    e = _parse(model, f.quadratic_terms)
    if length(f.affine_terms) > 0
        e = cpo_java_sum(cp, e, _parse(model, f.affine_terms))
    end
    if !iszero(f.constant)
        e = cpo_java_sum(cp, e, cpo_java_constant(cp, f.constant))
    end
    return e
end

function _parse(
    model::Optimizer,
    f::MOI.VectorOfVariables,
    ::MOI.SecondOrderCone,
)
    # SOC is the cone: t ≥ ||x||₂ ≥ 0. In quadratic form, this is
    # t² - Σᵢ xᵢ² ≥ 0 and t ≥ 0.
    # This function returns the expression t² - Σᵢ xᵢ².
    vars = [_info(model, v) for v in f.variables]
    e = cpo_java_square(model, vars[1].variable)
    for i in 2:length(vars)
        e = cpo_java_diff(model, e, vars[i].variable)
    end
    return e
end


text:
# wrapper to convert Nonlinear solver into LinearQuadratic solver

# To enable LPQP support from a Nonlinear solver, define, e.g.,
# LinearQuadraticModel(s::IpoptSolver) = NonlinearToLPQPBridge(NonlinearModel(s))

type QuadConstr
    linearidx::Vector{Int}
    linearval::Vector{Float64}
    quadrowidx::Vector{Int}
    quadcolidx::Vector{Int}
    quadval::Vector{Float64}
    sense::Char
    rhs::Float64
end
getdata(q::QuadConstr) = q.linearidx, q.linearval, q.quadrowidx, q.quadcolidx, q.quadval

type NonlinearToLPQPBridge <: AbstractLinearQuadraticModel
    nlpmodel::AbstractNonlinearModel
    LPdata
    Qobj::Tuple{Vector{Int},Vector{Int},Vector{Float64}}
    Qconstr::Vector{QuadConstr}
    numvar::Int
    numconstr::Int
    warmstart::Vector{Float64}
    vartypes::Vector{Symbol}
end

function NonlinearToLPQPBridge(m::AbstractNonlinearModel)
    return NonlinearToLPQPBridge(m,nothing,(Int[],Int[],Float64[]),QuadConstr[],0,0,Float64[],Symbol[])
end

export NonlinearToLPQPBridge

function loadproblem!(model::NonlinearToLPQPBridge, A, l, u, c, lb, ub, sense)
    model.LPdata = (A,l,u,c,lb,ub,sense)
    model.numvar = size(A,2)
    model.numconstr = size(A,1)
end

function setquadobj!(model::NonlinearToLPQPBridge, rowidx, colidx, quadval)
    model.Qobj = (rowidx, colidx, quadval)
end

function addquadconstr!(model::NonlinearToLPQPBridge, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
    push!(model.Qconstr, QuadConstr(linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs))
    model.numconstr += 1
end

function optimize!(model::NonlinearToLPQPBridge)
    A,l,u,c,lb,ub,sense = model.LPdata
    lb = copy(lb)
    ub = copy(ub)
    for q in model.Qconstr
        if q.sense == '<'
            push!(lb,-Inf)
            push!(ub,q.rhs)
        elseif q.sense == '>'
            push!(lb,q.rhs)
            push!(ub,Inf)
        else
            @assert q.sense == '='
            push!(lb,q.rhs)
            push!(ub,q.rhs)
        end
    end

    loadproblem!(model.nlpmodel, model.numvar, model.numconstr, l, u, lb, ub, sense, LPQPEvaluator(model))
    if any(model.vartypes .!= :Cont)
        setvartype!(model.nlpmodel, model.vartypes)
    end
    optimize!(model.nlpmodel)
end

for f in [:status, :getsolution, :getobjval, :getobjbound, :getreducedcosts]
    @eval $f(model::NonlinearToLPQPBridge) = $f(model.nlpmodel)
end

getsense(model::NonlinearToLPQPBridge) = model.LPdata[7]
numvar(model::NonlinearToLPQPBridge) = model.numvar
numconstr(model::NonlinearToLPQPBridge) = model.numconstr

numlinconstr(model::NonlinearToLPQPBridge) = size(model.LPdata[1],1)
numquadconstr(model::NonlinearToLPQPBridge) = length(model.Qconstr)

setvartype!(model::NonlinearToLPQPBridge,vartypes::Vector{Symbol}) = (model.vartypes = vartypes)

function getconstrsolution(model::NonlinearToLPQPBridge)
    x = getsolution(model)
    A = model.LPdata[1]
    return A*x
end

function getconstrduals(model::NonlinearToLPQPBridge)
    du = getconstrduals(model.nlpmodel)
    A = model.LPdata[1]
    return du[1:size(A,1)]
end

function getquadconstrduals(model::NonlinearToLPQPBridge)
    du = getconstrduals(model.nlpmodel)
    return du[(numlinconstr(model)+1):end]
end

type LPQPEvaluator <: AbstractNLPEvaluator
    A::SparseMatrixCSC{Float64,Int}
    c::Vector{Float64}
    Qi::Vector{Int} # quadratic objective terms
    Qj::Vector{Int}
    Qv::Vector{Float64}
    Qconstr::Vector{QuadConstr}
end

function LPQPEvaluator(model::NonlinearToLPQPBridge)
    A = model.LPdata[1]
    c = model.LPdata[4]
    Qi,Qj,Qv = model.Qobj
    return LPQPEvaluator(A,c,Qi,Qj,Qv,model.Qconstr)
end

features_available(d::LPQPEvaluator) = [:Grad, :Jac, :Hess]

function initialize(d::LPQPEvaluator, requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in features_available(d))
            error("Unsupported feature $feat")
            # TODO: implement Jac-vec and Hess-vec products
            # for solvers that need them
        end
    end
end

function eval_f(d::LPQPEvaluator, x)
    obj = dot(d.c,x)
    Qi = d.Qi
    Qj = d.Qj
    Qv = d.Qv
    for k in 1:length(Qi)
        if Qi[k] == Qj[k]
            obj += Qv[k]*x[Qi[k]]*x[Qi[k]]/2
        else
            obj += Qv[k]*x[Qi[k]]*x[Qj[k]]
        end
    end
    return obj
end

function eval_grad_f(d::LPQPEvaluator, grad_f, x)
    Qi = d.Qi
    Qj = d.Qj
    Qv = d.Qv
    for j = 1:length(d.c)
        grad_f[j] = d.c[j]
    end
    for k in 1:length(Qi)
        if Qi[k] == Qj[k]
            grad_f[Qi[k]] += Qv[k]*x[Qj[k]]
        else
            grad_f[Qi[k]] += Qv[k]*x[Qj[k]]
            grad_f[Qj[k]] += Qv[k]*x[Qi[k]]
        end
    end
end

function eval_g(d::LPQPEvaluator, g, x)
    g_val = d.A*x
    g[1:size(d.A,1)] = g_val
    for i in 1:length(d.Qconstr)
        linearidx, linearval, quadrowidx, quadcolidx, quadval = getdata(d.Qconstr[i])
        val = 0.0
        for k in 1:length(linearidx)
            val += linearval[k]*x[linearidx[k]]
        end
        for k in 1:length(quadrowidx)
            val += quadval[k]*x[quadrowidx[k]]*x[quadcolidx[k]]
        end
        g[i+size(d.A,1)] = val
    end
end

function jac_structure(d::LPQPEvaluator)
    jac_nnz = nnz(d.A)
    for i in 1:length(d.Qconstr)
        jac_nnz += length(d.Qconstr[i].linearidx)
        jac_nnz += 2*length(d.Qconstr[i].quadrowidx)
    end

    I = Vector{Int}(jac_nnz)
    J = Vector{Int}(jac_nnz)
    idx = 1
    for col = 1:size(d.A,2)
        for pos = d.A.colptr[col]:(d.A.colptr[col+1]-1)
            I[idx] = d.A.rowval[pos]
            J[idx] = col
            idx += 1
        end
    end
    for i in 1:length(d.Qconstr)
        linearidx, linearval, quadrowidx, quadcolidx, quadval = getdata(d.Qconstr[i])
        for k in 1:length(linearidx)
            I[idx] = i+size(d.A,1)
            J[idx] = linearidx[k]
            idx += 1
        end
        for k in 1:length(quadrowidx)
            I[idx] = i+size(d.A,1)
            I[idx+1] = i+size(d.A,1)
            J[idx] = quadrowidx[k]
            J[idx+1] = quadcolidx[k]
            idx += 2
        end
    end
    return I,J
end

function eval_jac_g(d::LPQPEvaluator, J, x)
    idx = 1
    for col = 1:size(d.A,2)
        for pos = d.A.colptr[col]:(d.A.colptr[col+1]-1)
            J[idx] = d.A.nzval[pos]
            idx += 1
        end
    end
    for i in 1:length(d.Qconstr)
        linearidx, linearval, quadrowidx, quadcolidx, quadval = getdata(d.Qconstr[i])
        for k in 1:length(linearidx)
            J[idx] = linearval[k]
            idx += 1
        end
        for k in 1:length(quadrowidx)
            J[idx] = quadval[k]*x[quadcolidx[k]]
            J[idx+1] = quadval[k]*x[quadrowidx[k]]
            idx += 2
        end
    end
end

function hesslag_structure(d::LPQPEvaluator)
    hess_nnz = length(d.Qi)
    for i in 1:length(d.Qconstr)
        hess_nnz += length(d.Qconstr[i].quadrowidx)
    end

    I = Vector{Int}(hess_nnz)
    J = Vector{Int}(hess_nnz)

    for k in 1:length(d.Qi)
        I[k] = d.Qi[k]
        J[k] = d.Qj[k]
    end
    idx = length(d.Qi)+1
    for i in 1:length(d.Qconstr)
        linearidx, linearval, quadrowidx, quadcolidx, quadval = getdata(d.Qconstr[i])
        for k in 1:length(quadrowidx)
            qidx1 = quadrowidx[k]
            qidx2 = quadcolidx[k]
            if qidx2 > qidx1
                qidx1, qidx2 = qidx2, qidx1
            end
            I[idx] = qidx1
            J[idx] = qidx2
            idx += 1
        end
    end
    return I,J
end

function eval_hesslag(d::LPQPEvaluator, H, x, σ, μ)
    for k in 1:length(d.Qi)
        H[k] = σ*d.Qv[k]
    end
    idx = length(d.Qi) + 1
    for i in 1:length(d.Qconstr)
        linearidx, linearval, quadrowidx, quadcolidx, quadval = getdata(d.Qconstr[i])
        for k in 1:length(quadrowidx)
            l = μ[size(d.A,1)+i]
            if quadrowidx[k] == quadcolidx[k]
                H[idx] = l*2*quadval[k]
            else
                H[idx] = l*quadval[k]
            end
            idx += 1
        end
    end
end

isobjlinear(d::LPQPEvaluator) = (length(d.Qi) == 0)
isobjquadratic(d::LPQPEvaluator) = true
isconstrlinear(d::LPQPEvaluator,i::Integer) = (i <= size(d.A,1))


text:
abstract type PadéExponentiator{T} <: AbstractExponentiator{T} end

# This implementation only supports time-independent operators, since
# the implicit step is pre-factorized.
struct CrankNicolson{T,R,F,B,V} <: PadéExponentiator{T}
    fwd::F
    bwd::B
    Δt::R
    tmp::V
end

CrankNicolson(A, Δt::R, μ::Number=one(R)) where R =
    CrankNicolson(I + μ*Δt*A/2, factorize(I - μ*Δt*A/2),
                  Vector{typeof(μ)}(undef, size(A,1)))

# Since the implicit step is pre-factorized, it is not possible to
# integrate over the time-step. Returning `nothing` signifies this.
operator(::CrankNicolson) = nothing

inplace(::CrankNicolson) = true
order(::CrankNicolson) = 2

function expv!(_, CN::CrankNicolson{T,R}, Δt::R, u)
    @assert Δt == CN.Δt
    mul!(CN.tmp, CN.fwd, u)
    ldiv!(u, CN.bwd, CN.tmp)
end


text:
using DataDeps

register(DataDep(
    "HAR",
    """
    Dataset: HAR
    Website: https://archive.ics.uci.edu/ml/datasets/human+activity+recognition+using+smartphones

    Observations: 10299 (7352 + 2947)
    Features:     561
    Classes:      6

    Jorge L. Reyes-Ortiz(1,2), Davide Anguita(1), Alessandro Ghio(1), Luca Oneto(1) and Xavier Parra(2)

        1 - Smartlab - Non-Linear Complex Systems Laboratory
        DITEN - Università degli Studi di Genova, Genoa (I-16145), Italy.

        2 - CETpD - Technical Research Centre for Dependency Care and Autonomous Living
        Universitat Politècnica de Catalunya (BarcelonaTech). Vilanova i la Geltrú (08800), Spain

    activityrecognition '@' smartlab.ws

    The experiments have been carried out with a group of 30 volunteers within an age bracket of 19-48 years. Each person performed six activities (WALKING, WALKING_UPSTAIRS, WALKING_DOWNSTAIRS, SITTING, STANDING, LAYING) wearing a smartphone (Samsung Galaxy S II) on the waist. Using its embedded accelerometer and gyroscope, we captured 3-axial linear acceleration and 3-axial angular velocity at a constant rate of 50Hz. The experiments have been video-recorded to label the data manually. The obtained dataset has been randomly partitioned into two sets, where 70% of the volunteers was selected for generating the training data and 30% the test data.

    The sensor signals (accelerometer and gyroscope) were pre-processed by applying noise filters and then sampled in fixed-width sliding windows of 2.56 sec and 50% overlap (128 readings/window). The sensor acceleration signal, which has gravitational and body motion components, was separated using a Butterworth low-pass filter into body acceleration and gravity. The gravitational force is assumed to have only low frequency components, therefore a filter with 0.3 Hz cutoff frequency was used. From each window, a vector of features was obtained by calculating variables from the time and frequency domain.

    Check the README.txt file for further details about this dataset.

    A video of the experiment including an example of the 6 recorded activities with one of the participants can be seen in the following link: http://www.youtube.com/watch?v=XOEN9W05_4A

    An updated version of this dataset can be found at

    http://archive.ics.uci.edu/ml/datasets/Smartphone-Based+Recognition+of+Human+Activities+and+Postural+Transitions.

    It includes labels of postural transitions between activities and also the full raw inertial signals instead of the ones pre-processed into windows.
    """,
    "https://archive.ics.uci.edu/ml/machine-learning-databases/00240/UCI%20HAR%20Dataset.zip",
    "2045e435c955214b38145fb5fa00776c72814f01b203fec405152dac7d5bfeb0",
    # customize post fetch method to combine files and save in compressed format
    post_fetch_method = (path -> begin
        DataDeps.unpack(path, keep_originals=false)

        # Read the data and labels files. X data are in fixed width format.
        tmpdir = "UCI HAR Dataset"
        train_dir = joinpath(tmpdir, "train")
        test_dir = joinpath(tmpdir, "test")
        data = vcat(
            CSV.read(joinpath(train_dir, "X_train.txt"), DataFrame, header=false, delim=' ', ignorerepeated=true),
            CSV.read(joinpath(test_dir, "X_test.txt"), DataFrame, header=false, delim=' ', ignorerepeated=true)
        )
        labels = vcat(
            CSV.read(joinpath(train_dir, "y_train.txt"), DataFrame, header=false),
            CSV.read(joinpath(test_dir, "y_test.txt"), DataFrame, header=false)
        )
        tmpdf = hcat(labels, data, makeunique=true)

        # Process the file as usual.
        MVDA.process_dataset(tmpdf,
            target_index=1,
            feature_indices=2:ncol(tmpdf),
            ext=".csv.gz"
        )

        # Clean up.
        rm("__MACOSX", recursive=true)
        rm(tmpdir, recursive=true)
    end),
))


text:
"""
Rational Krylov helper module. Contains functionality needed for methods
like NLEIGS and CORK.
"""
module RKHelper

include("inpolygon.jl")
include("discretizepolygon.jl")
include("rk_nep.jl")
include("rk_utils.jl")
include("linsolvercache.jl")

end


text:
using Gadfly, RDatasets, Dates

set_default_plot_size(6inch, 3inch)

economics = dataset("ggplot2", "economics")
dates = Date[Date(d) for d in economics[:Date]]
economics[:Date] = dates

p = plot(economics, x=:Date, y=:Unemploy, Geom.line)


text:
# __precompile__()
# module YAMLParameters

using YAML
# using Parameters

export paramyaml, add_paramyaml_specializer, load_yaml_config

"""
Convert string to symbol using conversionmap dictionary
If not found in the dictionary - intern symbol from the string
"""
function paramid{T <: AbstractString}(
        s :: AbstractString,
        conversionmap :: Dict{T, Symbol} = Dict{AbstractString, Symbol}())
    get(conversionmap, s, symbol(s))
end

"""
Gets name, id and description from the dictionary
"""
function nameiddesc(dict :: Dict)
    name = get(dict, "name", :notfound)
    id = get(dict, "id",:notfound)
    if name === :notfound && id === :notfound
        error("Either name of id must be provided!")
    elseif name === :notfound
        name = id
    elseif id === :notfound
        id = name
    end
    id = symbol(id)
    (name, id, get(dict, "description", ""))
end

"""
Converts dictionary to Parameter.
Dictionary must contain "value" key
"""
function parameter_yaml(m :: Dict)
    (name, id, description) = nameiddesc(m)
    value = m["value"]
    units = get(m, "units", "-")
    Parameter(
        name = name,
        id = id,
        value = value,
        description = description,
        units = units)
end

"""
Converts dictionary (with "children" key) to ParameterContainer
"""
function container_yaml(m :: Dict)
    (name, id, desc) = nameiddesc(m)
    children = m["children"]
    ParameterContainer(
        name = name,
        id = id,
        description = desc,
        children = map(paramyaml, children)
    )
end

"""
Converts dictionary (with "options" key) to ParameterOptions
"""
function options_yaml(m :: Dict)
    (name, id, desc) = nameiddesc(m)
    options = m["options"]
    selection = get(m, "selection", 1)
    ParameterOptions(
        name = name,
        id = id,
        description = desc,
        options = map(paramyaml,options),
        selection = selection
    )
end

"""
Converts dictionary (with "value" and "perturbation" keys) to
PerturbedParameter
"""
function perturbed_yaml(m :: Dict)
    (name, id, desc) = nameiddesc(m)
    value = m["value"]
    perturbation = m["perturbation"]
    units = get(m, "units", "-")
    PerturbedParameter(
        name = name,
        id = id,
        value = value,
        units = units,
        perturbation = perturbation
    )
end

"""
Dispatch dictionary for conversion
Each entry must be (default, Dict())
where Dict may contain further specializing arguments
"""
const paramyaml_conversions = Dict()

"""
Add new specializer for converting dictionary to <: AbstractParameter
"""
function add_paramyaml_specializer(
        f :: Function,
        names,
        dict = paramyaml_conversions)
    if length(names) == 1
        dict[names[1]] = (f, Dict())
    else
        if !haskey(dict, names[1])
            error("Overspecializing")
        else
            default, subdict = dict[names[1]]
            add_paramyaml_specializer(f, names[2:end], subdict)
        end
    end
    :done
end

begin
    add_paramyaml_specializer(parameter_yaml, ["value"])
    add_paramyaml_specializer(container_yaml, ["children"])
    add_paramyaml_specializer(options_yaml, ["options"])
    add_paramyaml_specializer(perturbed_yaml, ["value", "perturbation"])
end

"""
Construct parameter from dictionary `m` obtained from loading YAML file.

YAML files generally cannot contain parameter transformers. In which case,
a default transformer is used (usually `identity`). Otherwise, transformers
can be supplied by `trnsformers` parameter.

Optional dict is used to find the correct correspondence of `m` to a
particular parameter type. See `add_paramyaml_specializer`
"""
function paramyaml(m :: Dict, dict = paramyaml_conversions, default = identity)
    for key in keys(dict)
        if haskey(m, key)
            subdefault, subdict = dict[key]
            return paramyaml(m, subdict, subdefault)
        end
    end
    # didn't find a key: call default
    default(m)
end

function paramyaml(m :: Dict, transformers :: Dict{Symbol, Function})
    p = paramyaml(m)
    traverse(p, p ->
        begin
            if haskey(transformers, id(p))
                settransformer!(p, transformers[id(p)])
            end
            return true
        end)
    p
end

"""
Construct parameter from YAML config file
"""
function load_yaml_config(filename,
            transformers :: Dict{Symbol, Function} = Dict{Symbol, Function}())
    data = YAML.load_file(filename)
    paramyaml(data, transformers)
end

function update_param_from_config(p :: AbstractParameter, filename)
    data = YAML.load_file(filename)

end


text:
module MyModule

include("funcA.jl")
include("funcB.jl")

export myfuncA

end

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
module TestUserApi

using ..MyPetStore
using Swagger
using Test
using Random

const TEST_USER = "jlswag"
const TEST_USER1 = "jlswag1"
const TEST_USER2 = "jlswag2"

function test(uri)
    println("testing UserApi...")
    client = Swagger.Client(uri)
    api = UserApi(client)

    println("   - loginUser")
    login_result = loginUser(api, TEST_USER, "testpassword")
    @test !isempty(login_result)

    #println("   - createUser")
    #user1 = User(; id=100, username=TEST_USER1, firstName="test1", lastName="user1", email="jlswag1@example.com", password="testpass1", phone="1000000001", userStatus=0)
    #@test createUser(api, user1) === nothing

    #println("   - createUsersWithArrayInput")
    #user2 = User(; id=200, username=TEST_USER2, firstName="test2", lastName="user2", email="jlswag2@example.com", password="testpass2", phone="1000000002", userStatus=0)
    #@test createUsersWithArrayInput(api, [user1, user2]) === nothing

    #println("   - createUsersWithListInput")
    #@test createUsersWithListInput(api, [user1, user2]) === nothing

    println("   - getUserByName")
    @test_throws Swagger.ApiException getUserByName(api, randstring())
    @test_throws Swagger.ApiException getUserByName(api, TEST_USER)
    #getuser_result = getUserByName(api, TEST_USER)
    #@test isa(getuser_result, User)

    #println("   - updateUser")
    #@test updateUser(api, TEST_USER2, getuser_result) === nothing
    #println("   - deleteUser")
    #@test deleteUser(api, TEST_USER2) === nothing

    println("   - logoutUser")
    logout_result = logoutUser(api)
    @test logout_result === nothing

    nothing
end

end # module TestUserApi


text:
module ThreePhotons
using SpecialFunctions
using Distributions
using Optim
using ProgressMeter
using CurveFit
using PyPlot
using JLD

export
	#structures
	atomfactor,
	atomfactorlist,
	atomType,
	loadPDB,
	centerStructure,
	getDensity_cartesian,
	getIntensity_cartesian,
	getFourier_cartesian,
	createSphericalHarmonicsStructure,
	createCubicStructure,
	cubeToSphericalHarmonics,

	#cubic volumes
	CubeVolume,
	dr,
	qmax,
	getDensityCube,
	average,
	calculate_center_of_mass,
	translate_cube,
	center_cube,
	displacement_fitting,
	forward,
	absoluteSquare,
	mirrorCube,
	similarity,
	saveCube,
	loadCube,
	getVolumeInterpolated,
	rotateStructure,
	shell_correlation,
	shell_correlation_ISC,
	shell_correlation_FSC,
	ISC,
	FSC,

	#spherical harmonics
	Volume,
	SphericalVolume,
	SphericalHarmonicsVolume,
	SurfaceVolume,
	wigner,
	sphPlm,
	Plm,
	Ylm,
	get_size,
	num_val,
	num_coeff,
	dr,
	dq,
	qmax,
	average,
	disable_radial_interpolation,
	transform,
	backtransform,
	getSurfaceVolume,
	getSphericalHarmonicsVolume,
	seanindex,
	getc,
	setc,
	cvec_get,
	cvec_set,
	i_to_ang,
	ang_to_i,
	getSurface,
	getDensity,
	densityShell,
	fourierShell,
	intensityShell,
	deleteTerms,
	sphBesselTransform,
	forward,
	backward,
	absoluteSquare,
	shell_correlation,
	shell_correlation_ISC,
	shell_correlation_FSC,
	ISC,
	FSC,
	similarity,
	Umat,
	comp_to_real,
	real_to_comp,
	getVolumeInterpolated,
	getCube,
	saveCube,
	d_work,
	Dlms,
	rotationMatrix,
	rotateStructure,
	get_value,
	negativityCheck,

	#structure determination
	randomStartStructure,
	randomPertubation,
	rotation_search,
	get_MC_step,
	calculate_temperature,
	saveState,
	regular_action,
	rotate_all_at_once,
	rotate_hierarchical,
	increase_stepsize,
	decrease_stepsize,
	evaluate_step,
	checkRotationSearch,

	#phasing
	sRAAR,

	#plotting
	plot_random_2photon_slices,
	plot_random_3photon_slices,
	plotPoints2D,
	plotPointCloud,
	plot_scattering_image,
	compare_c2_grid,
	compare_histogram_with_theory,
	initialize_CUDA

import Base.real, Base.sum, Base.abs, Base.+, Base.-, Base.*, Base./, Base.length, Base.maximum

include("utilities.jl")
include("spherical_harmonics.jl")
include("cubic.jl")
include("structure.jl")
include("correlations.jl")
include("datagen.jl")
include("phases.jl")
include("data_processing.jl")
include("determination.jl")
include("exp_data.jl")
try
	include("cuda.jl")
catch e
	println(e)
	println("!!! Init of CUDA failed. Fallback to CPU.")
end

"""
A Julia package for the structure determination from single molecule X-ray scattering experiments at very low photon counts.
"""

end #module


text:
# This file includes the network

# --------------------------------------- ODENetwork construction and modificatio ------------------------------- #

"""
    $TYPEDEF

# Fields

    $TYPEDFIELDS
"""
mutable struct ODENetwork{T1<:AbstractNode,T2<:AbstractSimpleWeightedGraph,T3<:AbstractMatrix} <: AbstractNetwork
    "Dynamical nodes"
    nodes::Vector{T1}
    "Coupling strength"
    ϵ::Float64
    "Graph explaining topology"
    graph::T2
    "Inner coupling matrix"
    coupling::T3
end
ODENetwork(nodes = Lorenz[], ϵ = 1.0, graph = zeros(length(nodes), length(nodes)), coupling = diagm(ones(dimension(eltype(nodes))))) =
    ODENetwork(nodes, ϵ, SimpleWeightedGraph(graph), coupling)
ODENetwork(::Type{T}, ϵ = 1.0, graph = fill(0.0, 0, 0), coupling = diagm(ones(dimension(T)))) where {T<:AbstractNode} =
    ODENetwork([T() for i = 1:size(graph)[1]], ϵ, graph, coupling)

"""
    $TYPEDEF

# Fields

    $TYPEDFIELDS
"""
mutable struct SDENetwork{T1<:AbstractNode,T2<:AbstractSimpleWeightedGraph,T3<:AbstractMatrix,T4<:AbstractMatrix} <: AbstractNetwork
    "Dynamical nodes"
    nodes::Vector{T1}
    "Coupling strength"
    ϵ::Float64
    "Graph explaining topology"
    graph::T2
    "Noise coupling matrix"
    noisemat::T3
    "Inner coupling matrix"
    coupling::T4
end
function SDENetwork(nodes = Lorenz[], ϵ = 1.0, graph = zeros(length(nodes), length(nodes)), noisemat = nothing, coupling = fullstatecoupling(nodes))
    weightedgraph = SimpleWeightedGraph(graph)
    if noisemat === nothing
        noisemat = getnoisemat(weightedgraph)
    end
    SDENetwork(nodes, ϵ, weightedgraph, noisemat, coupling)
end
SDENetwork(::Type{T}, ϵ = 1.0, graph = fill(0.0, 0, 0), noisemat = nothing, coupling = diagm(ones(dimension(T)))) where {T<:AbstractNode} =
    SDENetwork([T() for i = 1:size(graph)[1]], ϵ, graph, noisemat, coupling)

"""
    $TYPEDSIGNATURES

Returns full noise coupling matrix corresponding to graph.
"""
function getnoisemat(graph)
    numnodes = nv(graph)
    numedges = ne(graph)
    mat = zeros(numnodes, numedges)
    for (k, edge) in enumerate(edges(graph))
        mat[edge.src, k] = 1
        mat[edge.dst, k] = -1
    end
    mat
end

# Returns the innner coupling matrix corresponding to the fully-state coupled network.
fullstatecoupling(nodes) = diagm(ones(dimension(eltype(nodes))))

"""
    $TYPEDSIGNATURES

Adds `node` to `net`
"""
function addnode!(net::AbstractNetwork, node::AbstractNode = eltype(net.nodes)())
    add_vertex!(net.graph)
    push!(net.nodes, node)
    net
end

"""
    $TYPEDSIGNATURES

Remoes the node with index `idx` from `net`.
"""
function remnode!(net::AbstractNetwork, idx::Int = lastindex(net.nodes))
    deleteat!(net.nodes, idx)
    rem_vertex!(net.graph, idx)
    net
end

"""
    $TYPEDSIGNATURES

Adds and edge in `net` between the nodes `src` and `dst` with weight `weight`
"""
function addedge!(net::AbstractNetwork, src::Int, dst::Int, weight::Real = 1.0, noisestrength::Real = 0.0)
    add_edge!(net.graph, src, dst, weight)
    is_directed(net.graph) && add_edge!(net.graph, dst, src, weight)
    # Update the noise matrix of the network. Tehe
    newcol = zeros(size(net.noisemat, 1), 1)
    newcol[src] = noisestrength
    newcol[dst] = -noisestrength
    net.noisemat = [net.noisemat newcol]
    net
end

"""
    $TYPEDSIGNATURES

Removes the edge between `src` and `dst` from `net`
"""
function remedge!(net::AbstractNetwork, src::Int, dst::Int)
    rem_edge!(net.graph, src, dst)
    is_directed(net.graph) && rem_edge!(net.graph, dst, src)
    # Find and delete the column of the noise matrix of the network.
    rows = net.noisemat[[src, dst], :]
    idx = [!iszero(col[1]) && !iszero(col[2]) for col in eachcol(rows)]
    net.noisemat = net.noisemat[:, .!idx]
    net
end


text:
infile = "./src/Day19/input.txt"
input = readlines(infile)

sep = findfirst(x -> x == "", input)
rules = input[1:(sep-1)]
messages = input[(sep+1):end]

ruledict = Dict{Int, String}()
for r in rules
    ids, spec = split(r, ": ")
    id = parse(Int, ids)
    ruledict[id] = string(spec)
end

function parserule(rid)
    spec = ruledict[rid]
    parts = string.(split(spec, " "))
    fresh = Vector{String}()
    for (i,p) in enumerate(parts)
        if contains(p,'"')
            # string literal
            push!(fresh, p)
        elseif contains(p, '|')
            # or function
            push!(fresh, p)
        else
            subid = parse(Int, p)
            if subid == rid
                if subid == 8
                    override = "$(fresh[i-1])+"
                    fresh = override
                    break
                elseif subid == 11
                    left = parserule(parse(Int, parts[i-1]))
                    right = parserule(parse(Int, parts[i+1]))
                    override = ["($(left)){$j}($(right)){$j}" for j in 1:5]
                    fresh = override
                    break
                end
            else
                subrule = parserule(subid)
                wrap = "($subrule)"
                push!(fresh, wrap)
            end
        end
    end
    updated = join(fresh, " ")
    ruledict[rid] = updated
end

backup = deepcopy(ruledict)
# ruledict = deepcopy(backup)
parserule(0)
function cleanrule(rid)
    rule = ruledict[rid]
    clean = replace(rule, "\"" => "")
    short = replace(clean, " " => "")
    # trimmed = replace(short, r"\(a\)\(a\)" => "(aa)")
    # trimmed = replace(trimmed, r"\(a\)\(b\)" => "(ab)")
    # trimmed = replace(trimmed, r"\(b\)\(a\)" => "(ba)")
    # trimmed = replace(trimmed, r"\(b\)\(b\)" => "(bb)")
    # trimmed = replace(trimmed, r"\(\(([ab]+)\)\)" => s"(\1)")
    ruledict[rid] = short
end
cleanrule(0)

reg = Regex("^$(ruledict[0])\$")
matches = match.(reg, messages)
matchey = filter(x -> x !== nothing, matches)
length(matchey)
backup[0]

# part2
# resets

ruledict[0] = backup[0]
ruledict[8] = "42 | 42 8"
ruledict[11] = "42 31 | 42 11 31"
parserule(0)
cleanrule(0)
cleanrule(8)
cleanrule(11)
cleanrule(42)
cleanrule(31)

customrules = ["($(ruledict[42])){1,}(($(ruledict[42])){$n}($(ruledict[31])){$n})" for n in 1:5]

regc = [Regex("^$(c)\$") for c in customrules]

matches = [length(filter(!isnothing, match.(re, messages))) for re in regc]
sum(matches)
mat = filter(!isnothing, matches);
length(mat)

reg42 = Regex("^($(ruledict[42])){2,}")
reg31 = Regex("($(ruledict[8])){1,}\$")
match42 = match.(reg42, messages);
match31 = match.(reg31, messages);

valid = 0
for m in messages
    left = match(reg42, m)
    if left === nothing
        continue
    end
    v = replace(m, reg42 => "_")
    println(v)
    right = match(reg31, v)
    if right === nothing
        continue
    else
        valid += 1
    end
end

valid