text:
module ColorVectorSpace

using ColorTypes, FixedPointNumbers, SpecialFunctions
using TensorCore
import TensorCore: ⊙, ⊗

using FixedPointNumbers: ShorterThanInt

import Base: ==, +, -, *, /, ^, <, ~
import Base: abs, clamp, convert, copy, div, eps, float,
    isfinite, isinf, isnan, isless, length, mapreduce, oneunit,
    promote_op, promote_rule, zero, trunc, floor, round, ceil, bswap,
    mod, mod1, rem, atan, hypot, max, min, real, typemin, typemax
# More unaryOps (mostly math functions)
import Base:      conj, sin, cos, tan, sinh, cosh, tanh,
                  asin, acos, atan, asinh, acosh, atanh,
                  sec, csc, cot, asec, acsc, acot,
                  sech, csch, coth, asech, acsch, acoth,
                  sinc, cosc, cosd, cotd, cscd, secd,
                  sind, tand, acosd, acotd, acscd, asecd,
                  asind, atand, rad2deg, deg2rad,
                  log, log2, log10, log1p, exponent, exp,
                  exp2, exp10, expm1, cbrt, sqrt,
                  significand, frexp, modf
import LinearAlgebra: norm, ⋅, dot, promote_leaf_eltypes  # norm1, norm2, normInf
import SpecialFunctions: gamma, logabsgamma, lfact
using Statistics
import Statistics: middle # and `_mean_promote`

export RGBRGB, complement, nan, dotc, dot, ⋅, hadamard, ⊙, tensor, ⊗, norm, varmult, stdmult

MathTypes{T,C<:Union{AbstractGray{T},AbstractRGB{T}}} = Union{C,TransparentColor{C,T}}

## Version compatibility with ColorTypes
### TODO: Remove the definitons other than `one` when dropping ColorTypes v0.10 support

if !hasmethod(zero, (Type{TransparentGray},))
    zero(::Type{C}) where {C<:TransparentGray} = C(0,0)
    zero(::Type{C}) where {C<:AbstractRGB}     = C(0,0,0)
    zero(::Type{C}) where {C<:TransparentRGB}  = C(0,0,0,0)
    zero(p::Colorant) = zero(typeof(p))
end

if !hasmethod(one, (Type{Gray24},))
    Base.one(::Type{Gray24}) = Gray24(1)
end

if !hasmethod(one, (Type{TransparentGray},)) # specification change is planned for ColorTypes v0.12
    Base.one(::Type{C}) where {C<:TransparentGray} = C(1,1)
    Base.one(::Type{C}) where {C<:AbstractRGB}     = C(1,1,1)
    Base.one(::Type{C}) where {C<:TransparentRGB}  = C(1,1,1,1)
    Base.one(p::Colorant) = one(typeof(p))
end

if !hasmethod(isfinite, (Colorant,))
    isfinite(c::Colorant) = mapreducec(isfinite, &, true, c)
    isinf(c::Colorant) = mapreducec(isinf, |, false, c)
    isnan(c::Colorant) = mapreducec(isnan, |, false, c)
end

if which(<, Tuple{AbstractGray,AbstractGray}).module === Base
    (<)(g1::AbstractGray, g2::AbstractGray) = gray(g1) < gray(g2)
    (<)(c::AbstractGray, r::Real) = gray(c) < r
    (<)(r::Real, c::AbstractGray) = r < gray(c)
end
if !hasmethod(isless, Tuple{AbstractGray,Real})
    isless(c::AbstractGray, r::Real) = isless(gray(c), r)
    isless(r::Real, c::AbstractGray) = isless(r, gray(c))
end

if isdefined(ColorTypes, :nan)
    using ColorTypes: nan
else
    nan(::Type{T}) where {T<:AbstractFloat} = convert(T, NaN)
    nan(::Type{C}) where {T<:AbstractFloat, C<:MathTypes{T}} = mapc(_ -> nan(T), zero(C))
end

if which(real, (Type{<:AbstractGray},)).module === Base
    real(::Type{C}) where {C<:AbstractGray} = real(eltype(C))
end

# To help type inference
promote_rule(::Type{T}, ::Type{C}) where {T<:Real,C<:AbstractGray} = promote_type(T, eltype(C))

promote_leaf_eltypes(x::Union{AbstractArray{T},Tuple{T,Vararg{T}}}) where {T<:MathTypes} = eltype(T)

if isdefined(Statistics, :_mean_promote)
    Statistics._mean_promote(x::MathTypes, y::MathTypes) = mapc(FixedPointNumbers.Treduce, y)
end

## Traits and key utilities

# Return eltypes for arithmetic operations
multype(::Type{A}, ::Type{B}) where {A,B} = coltype(typeof(zero(A)*zero(B)))
sumtype(::Type{A}, ::Type{B}) where {A,B} = coltype(typeof(zero(A)+zero(B)))
divtype(::Type{A}, ::Type{B}) where {A,B} = coltype(typeof(zero(A)/oneunit(B)))
divtype(::Type{A}, ::Type{B}) where {A,B<:FixedPoint} = coltype(typeof(zero(A)/typemax(B)))
powtype(::Type{A}, ::Type{B}) where {A,B} = coltype(typeof(zero(A)^zero(B)))
sumtype(a::Colorant, b::Colorant) = coltype(sumtype(eltype(a),eltype(b)))

multype(::Type{Bool}, ::Type{B}) where {B} = B
multype(::Type{A}, ::Type{Bool}) where {A} = A <: Integer ? multype(A, N0f8) : A
multype(::Type{Bool}, ::Type{Bool}) = Bool
sumtype(::Type{Bool}, ::Type{B}) where {B} = B <: Integer ? N0f8 : sumtype(N0f8, B) # FIXME
sumtype(::Type{A}, ::Type{Bool}) where {A} = sumtype(Bool, A)
sumtype(::Type{Bool}, ::Type{Bool}) = N0f8 # FIXME
divtype(::Type{Bool}, ::Type{B}) where {B} = divtype(B <: Integer ? N0f8 : UInt8, B)
divtype(::Type{Bool}, ::Type{B}) where {B<:FixedPoint} = divtype(N0f8, B) # FIXME
divtype(::Type{A}, ::Type{Bool}) where {A} = divtype(A, A <: FixedPoint ? N0f8 : UInt8) # FIXME
divtype(::Type{A}, ::Type{Bool}) where {A<:Integer} = Float32 # FIXME
divtype(::Type{Bool}, ::Type{Bool}) = divtype(N0f8, UInt8)
powtype(::Type{Bool}, ::Type{B}) where {B} = B <: Integer ? Bool : B
powtype(::Type{A}, ::Type{Bool}) where {A} = A
powtype(::Type{Bool}, ::Type{Bool}) = Bool

coltype(::Type{T}) where {T<:Fractional} = T
coltype(::Type{T}) where {T<:Number}     = floattype(T)

acctype(::Type{T}) where {T<:FixedPoint}        = floattype(T)
acctype(::Type{T}) where {T<:ShorterThanInt}    = Int
acctype(::Type{Rational{T}}) where {T<:Integer} = typeof(zero(T)/oneunit(T))
acctype(::Type{T}) where {T<:Real}              = T

acctype(::Type{T1}, ::Type{T2}) where {T1,T2} = acctype(promote_type(T1, T2))

acc(x::Number) = convert(acctype(typeof(x)), x)

color_rettype(::Type{C1}, ::Type{C2}) where {C1, C2} = base_colorant_type(promote_type(C1, C2))
color_rettype(c1::Colorant, c2::Colorant) = color_rettype(typeof(c1), typeof(c2))

arith_colorant_type(::C) where {C<:Colorant} = arith_colorant_type(C)
function arith_colorant_type(::Type{C}) where {C<:Colorant}
    Cb = base_colorant_type(C)
    isconcretetype(C) && C === Cb && return _arith_colorant_type(C) # non-parametric
    return Cb
end
_arith_colorant_type(::Type{<:AbstractGray})    = Gray
_arith_colorant_type(::Type{<:TransparentGray}) = AGray
_arith_colorant_type(::Type{<:AbstractGrayA})   = GrayA
_arith_colorant_type(::Type{<:AbstractRGB})     = RGB
_arith_colorant_type(::Type{<:TransparentRGB})  = ARGB
_arith_colorant_type(::Type{<:AbstractRGBA})    = RGBA

parametric(::Type{C}, ::Type{T}) where {C,T} = C{T}
parametric(::Type{C}, ::Type{T}) where {T, C<:Colorant{T}} = C # e.g. parametric(RGB24, N0f8) == RGB24

rettype(::typeof(+), a::C, b::C) where {C <: Colorant} = C
rettype(::typeof(-), a::C, b::C) where {C <: Colorant} = C
rettype(::typeof(+), a, b) = parametric(color_rettype(a, b), sumtype(a, b))
rettype(::typeof(-), a, b) = parametric(color_rettype(a, b), sumtype(a, b))
rettype(::typeof(*), a, b) = parametric(color_rettype(a, b), multype(eltype(a), eltype(b))) # gray * gray
rettype(::typeof(*), a::Real, b) = arith_colorant_type(b){multype(typeof(a), eltype(b))}
rettype(::typeof(/), a, b::Real) = arith_colorant_type(a){divtype(eltype(a), typeof(b))}
rettype(::typeof(/), a, b) = arith_colorant_type(a){divtype(eltype(a), eltype(b))}
rettype(::typeof(^), a, b)          = arith_colorant_type(a){powtype(eltype(a), typeof(b))}
rettype(::typeof(^), a, b::Integer) = arith_colorant_type(a){powtype(eltype(a), Int)}

# Useful for leveraging iterator algorithms. Don't use this externally, as the implementation may change.
channels(c::AbstractGray)    = (gray(c),)
channels(c::TransparentGray) = (gray(c), alpha(c))
channels(c::AbstractRGB)     = (red(c), green(c), blue(c))
channels(c::TransparentRGB)  = (red(c), green(c), blue(c), alpha(c))

if reinterpret(N0f16, 0xBADb)^2 === N0f16(float(reinterpret(N0f16, 0xBADb))^2)
    function _mul(x::N0f8, y::N0f8)
        m = widemul(reinterpret(x), reinterpret(y))
        z = (((m + 0x80) >> 0x8) + m + 0x80) >> 0x8
        reinterpret(N0f8, z % UInt8)
    end
    function _mul(x::N0f16, y::N0f16)
        m = widemul(reinterpret(x), reinterpret(y))
        z = (((m + 0x8000) >> 0x10) + m + 0x8000) >> 0x10
        reinterpret(N0f16, z % UInt16)
    end
end
_mul(x::T, y::T) where {T} = x * y
_mul(x, y) = (T = multype(typeof(x), typeof(y)); _mul(T(x), T(y)))

function _div(x::T, y::T) where {T <: FixedPoint}
    F = floattype(T)
    # range check should be done in the color constructor
    F(reinterpret(x)) / F(reinterpret(y))
end
function _div(x::AbstractFloat, y::Normed)
    F = divtype(typeof(x), typeof(y))
    F(x) * F(reinterpret(oneunit(y))) / F(reinterpret(y))
end
_div(x::AbstractFloat, y::Integer) = _mul(x, oneunit(x) / y)
_div(x::T, y::T) where {T} = x / y
_div(x::Bool, y::Bool) = (T = divtype(typeof(x), typeof(y)); _div(T(x), T(y)))
_div(x::FixedPoint, y::Bool) = (T = divtype(typeof(x), typeof(y)); T(_div(T(x), T(y)))) # FIXME
_div(x::Bool, y::FixedPoint) = (T = divtype(typeof(x), typeof(y)); T(_div(T(x), T(y)))) # FIXME
_div(x, y) = (T = divtype(typeof(x), typeof(y)); _div(T(x), T(y)))

@inline _mapc(::Type{C}, f, c) where {C<:MathTypes} = C(f.(channels(c))...)
@inline _mapc(::Type{C}, f, a, b) where {C<:MathTypes} = C(f.(channels(a), channels(b))...)

## Generic algorithms
Base.add_sum(c1::MathTypes,c2::MathTypes) = mapc(Base.add_sum, c1, c2)
Base.add_sum(c1::MathTypes{Bool}, c2::MathTypes{Bool}) = mapc((x1, x2) -> FixedPointNumbers.Treduce(x1) + FixedPointNumbers.Treduce(x2), c1, c2)
Base.reduce_first(::typeof(Base.add_sum), c::MathTypes) = mapc(x->Base.reduce_first(Base.add_sum, x), c)
Base.reduce_first(::typeof(Base.add_sum), c::MathTypes{Bool}) = mapc(x->Base.reduce_first(Base.add_sum, N0f8(x)), c)
function Base.reduce_empty(::typeof(Base.add_sum), ::Type{T}) where {T<:MathTypes}
    z = Base.reduce_empty(Base.add_sum, eltype(T))
    return zero(base_colorant_type(T){typeof(z)})
end
Base.reduce_empty(::typeof(Base.add_sum), ::Type{T}) where {T<:MathTypes{Bool}} =
    Base.reduce_empty(Base.add_sum, base_colorant_type(T){N0f8})

## Rounding & mod
for f in (:trunc, :floor, :round, :ceil, :eps, :bswap)
    @eval $f(g::Gray{T}) where {T} = Gray{T}($f(gray(g)))
end
eps(::Type{Gray{T}}) where {T} = Gray(eps(T))

for f in (:trunc, :floor, :round, :ceil)
    @eval $f(::Type{T}, g::Gray) where {T<:Integer} = Gray{T}($f(T, gray(g)))
end

for f in (:mod, :rem, :mod1)
    @eval $f(x::Gray, m::Gray) = Gray($f(gray(x), gray(m)))
    @eval $f(x::Gray, m::Real) = Gray($f(gray(x), m))
end

dotc(x::T, y::T) where {T<:Real} = acc(x)*acc(y)
dotc(x::Real, y::Real) = dotc(promote(x, y)...)

"""
    y = complement(x)

Take the complement `1-x` of `x`.  If `x` is a color with an alpha channel,
the alpha channel is left untouched. Don't forget to add a dot when `x` is
an array: `complement.(x)`
"""
complement(x::Union{Number,Colorant}) = oneunit(x)-x
complement(x::TransparentColor) = typeof(x)(complement(color(x)), alpha(x))

## Math on Colors. These implementations encourage inlining and,
## for the case of Normed types, nearly halve the number of multiplications (for RGB)

# Common
copy(c::MathTypes) = c
(*)(f::Real, c::MathTypes) = _mapc(rettype(*, f, c), v -> _mul(f, v), c)
(*)(c::MathTypes, f::Real) = (*)(f, c)
(/)(c::MathTypes, f::Real) = _mapc(rettype(/, c, f), v -> _div(v, f), c)
(+)(c::MathTypes) = mapc(+, c)
(+)(c::MathTypes{Bool}) = c
(-)(c::MathTypes) = mapc(-, c)
(-)(c::MathTypes{Bool}) = c
abs(c::MathTypes) = mapc(abs, c)
norm(c::MathTypes, p::Real=2) = (cc = channels(c); norm(cc, p)/(p == 0 ? length(cc) : length(cc)^(1/p)))
(⊙)(a::C, b::C) where {C<:MathTypes} = _mapc(rettype(*, a, b), _mul, a, b)
(⋅)(a::C, b::C) where {C<:MathTypes} = throw(MethodError(dot, (a, b)))
(⊗)(a::C, b::C) where {C<:MathTypes} = throw(MethodError(tensor, (a, b)))

## Mixed types
(+)(a::MathTypes, b::MathTypes) = (+)(promote(a, b)...)
(-)(a::MathTypes, b::MathTypes) = (-)(promote(a, b)...)
(⊙)(a::MathTypes, b::MathTypes) = (⊙)(promote(a, b)...)
(⋅)(a::MathTypes, b::MathTypes) = (⋅)(promote(a, b)...) # not fully supported, but used for error hints
(⊗)(a::MathTypes, b::MathTypes) = (⊗)(promote(a, b)...) # not fully supported, but used for error hints

# Scalar RGB
function (/)(c::C, f::AbstractFloat) where {C<:Union{AbstractRGB, TransparentRGB}}
    r = oneunit(divtype(eltype(c), typeof(f))) / f
    _mapc(rettype(/, c, f), v -> v * r, c)
end
(+)(a::AbstractRGB, b::AbstractRGB) = _mapc(rettype(+, a, b), +, a, b)
(-)(a::AbstractRGB, b::AbstractRGB) = _mapc(rettype(-, a, b), -, a, b)
(+)(a::TransparentRGB, b::TransparentRGB) = _mapc(rettype(+, a, b), +, a, b)
(-)(a::TransparentRGB, b::TransparentRGB) = _mapc(rettype(-, a, b), -, a, b)

# New multiplication operators
function (⋅)(x::AbstractRGB, y::AbstractRGB)
    T = acctype(eltype(x), eltype(y))
    (T(red(x))*T(red(y)) + T(green(x))*T(green(y)) + T(blue(x))*T(blue(y)))/3
end
# ⊗ defined below

# These constants come from squaring the conversion to grayscale
# (rec601 luma), and normalizing
dotc(x::T, y::T) where {T<:AbstractRGB} = 0.200f0 * acc(red(x))*acc(red(y)) + 0.771f0 * acc(green(x))*acc(green(y)) + 0.029f0 * acc(blue(x))*acc(blue(y))
dotc(x::AbstractRGB, y::AbstractRGB) = dotc(promote(x, y)...)

# Scalar Gray
const unaryOps = (:~, :conj, :abs,
                  :sin, :cos, :tan, :sinh, :cosh, :tanh,
                  :asin, :acos, :atan, :asinh, :acosh, :atanh,
                  :sec, :csc, :cot, :asec, :acsc, :acot,
                  :sech, :csch, :coth, :asech, :acsch, :acoth,
                  :sinc, :cosc, :cosd, :cotd, :cscd, :secd,
                  :sind, :tand, :acosd, :acotd, :acscd, :asecd,
                  :asind, :atand, :rad2deg, :deg2rad,
                  :log, :log2, :log10, :log1p, :exponent, :exp,
                  :exp2, :exp10, :expm1, :cbrt, :sqrt,
                  :significand,
                  :gamma, :lfact, :frexp, :modf,
                  :(SpecialFunctions.erf), :(SpecialFunctions.erfc),
                  :(SpecialFunctions.erfcx), :(SpecialFunctions.erfi), :(SpecialFunctions.dawson),
                  :(SpecialFunctions.airyai),
                  :(SpecialFunctions.airyaiprime), :(SpecialFunctions.airybi), :(SpecialFunctions.airybiprime),
                  :(SpecialFunctions.besselj0), :(SpecialFunctions.besselj1), :(SpecialFunctions.bessely0), :(SpecialFunctions.bessely1),
                  :(SpecialFunctions.eta), :(SpecialFunctions.zeta), :(SpecialFunctions.digamma))
for op in unaryOps
    @eval ($op)(c::AbstractGray) = Gray($op(gray(c)))
end
function logabsgamma(c::AbstractGray)
    lagc, s = logabsgamma(gray(c))
    return Gray(lagc), s
end

middle(c::AbstractGray) = arith_colorant_type(c)(middle(gray(c)))
middle(x::C, y::C) where {C<:AbstractGray} = arith_colorant_type(C)(middle(gray(x), gray(y)))

(/)(n::Number, c::AbstractGray) = base_color_type(c)(_div(real(n), gray(c)))
(+)(a::AbstractGray,       b::AbstractGray)       = _mapc(rettype(+, a, b), +, a, b)
(+)(a::AbstractGray{Bool}, b::AbstractGray{Bool}) = _mapc(rettype(+, a, b), ⊻, a, b)
(+)(a::TransparentGray,    b::TransparentGray)    = _mapc(rettype(+, a, b), +, a, b)
(-)(a::AbstractGray,       b::AbstractGray)       = _mapc(rettype(-, a, b), -, a, b)
(-)(a::AbstractGray{Bool}, b::AbstractGray{Bool}) = _mapc(rettype(-, a, b), ⊻, a, b)
(-)(a::TransparentGray,    b::TransparentGray)    = _mapc(rettype(-, a, b), -, a, b)
(*)(a::AbstractGray, b::AbstractGray) = a ⊙ b
(^)(a::AbstractGray, b::Integer) = rettype(^, a, b)(gray(a)^convert(Int,b))
(^)(a::AbstractGray, b::Real)    = rettype(^, a, b)(gray(a)^b)
(/)(a::AbstractGray, b::AbstractGray) = rettype(/, a, b)(_div(gray(a), gray(b)))
(+)(a::AbstractGray, b::Number) = base_color_type(a)(gray(a)+b)
(+)(a::AbstractGray{Bool}, b::Number) = arith_colorant_type(a){sumtype(Bool, typeof(b))}(gray(a) + b)
(+)(a::AbstractGray{Bool}, b::Bool) = a + N0f8(b) # FIXME
(+)(a::Number, b::AbstractGray) = b+a
(-)(a::AbstractGray, b::Number) = base_color_type(a)(gray(a)-b)
(-)(a::Number, b::AbstractGray) = base_color_type(b)(a-gray(b))
(-)(a::AbstractGray{Bool}, b::Number) = arith_colorant_type(a){sumtype(Bool, typeof(b))}(gray(a) - b)
(-)(a::Number, b::AbstractGray{Bool}) = arith_colorant_type(b){sumtype(typeof(a), Bool)}(a - gray(b))
(-)(a::AbstractGray{Bool}, b::Bool) = a - N0f8(b) # FIXME
(-)(a::Bool, b::AbstractGray{Bool}) = N0f8(a) - b # FIXME

(⋅)(x::C, y::C) where {C<:AbstractGray} = _mul(gray(x), gray(y))
(⊗)(x::C, y::C) where {C<:AbstractGray} = x ⊙ y

max(a::T, b::T) where {T<:AbstractGray} = T(max(gray(a),gray(b)))
max(a::AbstractGray, b::AbstractGray) = max(promote(a,b)...)
max(a::Number, b::AbstractGray) = max(promote(a,b)...)
max(a::AbstractGray, b::Number) = max(promote(a,b)...)
min(a::T, b::T) where {T<:AbstractGray} = T(min(gray(a),gray(b)))
min(a::AbstractGray, b::AbstractGray) = min(promote(a,b)...)
min(a::Number, b::AbstractGray) = min(promote(a,b)...)
min(a::AbstractGray, b::Number) = min(promote(a,b)...)

atan(x::AbstractGray, y::AbstractGray)  = atan(gray(x), gray(y))
hypot(x::AbstractGray, y::AbstractGray) = hypot(gray(x), gray(y))

dotc(x::AbstractGray, y::AbstractGray) = _mul(acc(gray(x)), acc(gray(y)))

typemin(::Type{C}) where {C<:AbstractGray} = C(typemin(eltype(C)))
typemax(::Type{C}) where {C<:AbstractGray} = C(typemax(eltype(C)))
typemin(c::AbstractGray) = typemin(typeof(c))
typemax(c::AbstractGray) = typemax(typeof(c))

## RGB tensor products

"""
    RGBRGB(rr, gr, br, rg, gg, bg, rb, gb, bb)

Represent the [tensor product](https://en.wikipedia.org/wiki/Tensor_product) of two RGB values.

# Example

```jldoctest
julia> a, b = RGB(0.2f0, 0.3f0, 0.5f0), RGB(0.77f0, 0.11f0, 0.22f0)
(RGB{Float32}(0.2f0,0.3f0,0.5f0), RGB{Float32}(0.77f0,0.11f0,0.22f0))

julia> a ⊗ b
RGBRGB{Float32}:
 0.154  0.022  0.044
 0.231  0.033  0.066
 0.385  0.055  0.11
"""
struct RGBRGB{T}
    rr::T
    gr::T
    br::T
    rg::T
    gg::T
    bg::T
    rb::T
    gb::T
    bb::T
end
function RGBRGB{T}(mat::AbstractMatrix) where T
    size(mat) == (3, 3) || throw(DimensionMismatch("`mat` must be a 3×3 matrix."))
    @inbounds RGBRGB{T}(NTuple{9, T}(T(mat[I]) for I in CartesianIndices(mat))...)
end
RGBRGB(mat::AbstractMatrix{T}) where T = RGBRGB{T}(mat)

Base.eltype(::Type{RGBRGB{T}}) where T = T
Base.Matrix{T}(p::RGBRGB) where T = T[p.rr p.rg p.rb;
                                      p.gr p.gg p.gb;
                                      p.br p.bg p.bb]
Base.Matrix(p::RGBRGB{T}) where T = Matrix{T}(p)

function Base.show(io::IO, @nospecialize(p::RGBRGB))
    print(io, "RGBRGB{", eltype(p), "}([")
    ioc = IOContext(io, :typeinfo => eltype(p))
    print(ioc, p.rr, " ", p.rg, " ", p.rb, "; ")
    print(ioc, p.gr, " ", p.gg, " ", p.gb, "; ")
    print(ioc, p.br, " ", p.bg, " ", p.bb, "])")
end
function Base.show(io::IO, ::MIME"text/plain", @nospecialize(p::RGBRGB))
    println(io, "RGBRGB{", eltype(p), "}:")
    ioc = IOContext(io, :typeinfo => eltype(p), :compact => get(io, :compact, true))
    Base.print_matrix(ioc, Matrix(p))
end

Base.zero(::Type{RGBRGB{T}}) where T = (z = zero(T); RGBRGB(z, z, z, z, z, z, z, z, z))
Base.zero(a::RGBRGB) = zero(typeof(a))

+(a::RGBRGB) = a
-(a::RGBRGB) = RGBRGB(-a.rr, -a.gr, -a.br, -a.rg, -a.gg, -a.bg, -a.rb, -a.gb, -a.bb)
+(a::RGBRGB, b::RGBRGB) = RGBRGB(a.rr + b.rr, a.gr + b.gr, a.br + b.br,
                                 a.rg + b.rg, a.gg + b.gg, a.bg + b.bg,
                                 a.rb + b.rb, a.gb + b.gb, a.bb + b.bb)
-(a::RGBRGB, b::RGBRGB) = +(a, -b)
*(α::Real, a::RGBRGB) = RGBRGB(α*a.rr, α*a.gr, α*a.br, α*a.rg, α*a.gg, α*a.bg, α*a.rb, α*a.gb, α*a.bb)
*(a::RGBRGB, α::Real) = α*a
/(a::RGBRGB, α::Real) = (1/α)*a

function ⊗(a::AbstractRGB, b::AbstractRGB)
    ar, ag, ab = red(a), green(a), blue(a)
    br, bg, bb = red(b), green(b), blue(b)
    agbr, abbg, arbb, abbr, arbg, agbb = ag*br, ab*bg, ar*bb, ab*br, ar*bg, ag*bb
    return RGBRGB(ar*br, agbr, abbr, arbg, ag*bg, abbg, arbb, agbb, ab*bb)
end

"""
    varmult(op, itr; corrected::Bool=true, mean=Statistics.mean(itr), dims=:)

Compute the variance of elements of `itr`, using `op` as the multiplication operator.
The keyword arguments behave identically to those of `Statistics.var`.

# Example

```julia
julia> cs = [RGB(0.2, 0.3, 0.4), RGB(0.5, 0.3, 0.2)]
2-element Array{RGB{Float64},1} with eltype RGB{Float64}:
 RGB{Float64}(0.2,0.3,0.4)
 RGB{Float64}(0.5,0.3,0.2)

julia> varmult(⋅, cs)
0.021666666666666667

julia> varmult(⊙, cs)
RGB{Float64}(0.045,0.0,0.020000000000000004)

julia> varmult(⊗, cs)
RGBRGB{Float64}:
  0.045  0.0  -0.03
  0.0    0.0   0.0
 -0.03   0.0   0.02
```
"""
function varmult(op, itr; corrected::Bool=true, dims=:, mean=Statistics.mean(itr; dims=dims))
    if dims === (:)
        v = mapreduce(c->(Δc = c-mean; op(Δc, Δc)), +, itr; dims=dims)
        n = length(itr)
    else
        # TODO: avoid temporary creation
        v = mapreduce(Δc->op(Δc, Δc), +, itr .- mean; dims=dims)
        n = length(itr) // length(v)
    end
    return v / (corrected ? max(1, n-1) : max(1, n))
end
function stdmult(op, itr; kwargs...)
    _sqrt(x::Real)        = sqrt(x)
    _sqrt(c::Colorant)    = mapc(sqrt, c)
    _sqrt(rgbrgb::RGBRGB) = RGBRGB(sqrt.(Matrix(rgbrgb)))

    result = varmult(op, itr; kwargs...)
    return isa(result, Union{Real,Colorant,RGBRGB}) ? _sqrt(result) : _sqrt.(result)
end

Base.length(r::StepRange{<:AbstractGray,<:AbstractGray}) = length(StepRange(gray(r.start), gray(r.step), gray(r.stop)))
Base.length(r::StepRange{<:AbstractGray})                = length(StepRange(gray(r.start), r.step,       gray(r.stop)))

Base.abs2(g::AbstractGray) = abs2(gray(g))
function Base.abs2(c::AbstractRGB)
    Base.depwarn("""
    The return value of `abs2` will change to ensure that `abs2(g::Gray) ≈ abs2(RGB(g::Gray))`.
    For `RGB` colors, this results in dividing the previous output by 3.

    To avoid this warning, use `ColorVectorSpace.Future.abs2` instead of `abs2`; currently,
    `abs2` returns the old value (for compatibility), and `ColorVectorSpace.Future.abs2` returns the new value.
    When making this change, you may also need to adjust constants like color-difference thresholds
    to compensate for the change in the returned value.

    If you are getting this from `var`, use `varmult` instead.
    """, :abs2)
    return mapreducec(v->float(v)^2, +, zero(eltype(c)), c)
end

module Future
using ..ColorTypes
using ..ColorVectorSpace: ⋅, dot
"""
    ColorVectorSpace.Future.abs2(c)

Return a scalar "squared magnitude" for color types. For RGB and gray, this is just the mean-square
channelwise intensity.

Compatibility note: this gives a different result from `Base.abs2(c)`, but eventually `Base.abs2` will switch
to the definition used here. Using `ColorVectorSpace.Future.abs2` thus future-proofs your code.
For more information about the transition, see ColorVectorSpace's README.
"""
abs2(c::Union{Real,AbstractGray,AbstractRGB}) = c ⋅ c
end

function __init__()
    if isdefined(Base, :Experimental) && isdefined(Base.Experimental, :register_error_hint)
        Base.Experimental.register_error_hint(MethodError) do io, exc, argtypes, kwargs
            if exc.f === dot && length(argtypes) == 2
                if any(C -> C <: Union{TransparentGray, TransparentRGB}, argtypes)
                    print(io, "\n`dot` (or `⋅`) for transparent colors is not supported ")
                    print(io, "because the normalization is designed for opaque grays and RGBs.")
                    print(io, "\nYou can use `reducec(+, 0, mapc(float, a) ⊙ mapc(float, b))` ")
                    print(io, "to get the dot product without normalization.")
                end
            elseif exc.f === tensor && length(argtypes) == 2
                if any(C -> C <: Union{TransparentGray, TransparentRGB}, argtypes)
                    print(io, "\n`tensor` (or `⊗`) for transparent colors is not supported ")
                    print(io, "due to the non-high demand.")
                end
            elseif exc.f === abs2 && length(argtypes) == 1
                C = argtypes[1]
                if C <: MathTypes
                    print(io, "\nIt is ambiguous whether `abs2(::$C)` should return a real number ")
                    print(io, "or a color object.")
                    print(io, "\nYou can use `_abs2(c) = mapreducec(v->float(v)^2, +, 0, c)` ")
                    print(io, "to get the value compatible with ColorVectorSpace v0.8 or earlier.")
                end
            end
        end
    end
end

## Precompilation

if Base.VERSION >= v"1.4.2"
    include("precompile.jl")
    _precompile_()
end

end


text:
# Test cornerfeatures.jl

println("testing cornerfeatures")

using ImageProjectiveGeometry
using Base.Test

tol = 1;

# structuretensor, harris, noble, shi_tomasi, coherence
# hessianfeatures
# fastradial
# nonmaxsuppts

# Create an image with an isolated point
img = zeros(100,100)
r0 = 20;
c0 = 40;
img[r0,c0] = 10

(cimg, r, c) = harris(img, N=1,radius=5,subpixel=false)
@test abs(r[1]-r0) < tol && abs(c[1]-c0) < tol

(cimg, r, c) = noble(img, N=1,radius=5,subpixel=false)
@test abs(r[1]-r0) < tol && abs(c[1]-c0) < tol

(cimg, r, c) = shi_tomasi(img, N=1,radius=5,subpixel=false)
@test abs(r[1]-r0) < tol && abs(c[1]-c0) < tol

hdet = hessianfeatures(img, 1)
(r,c) = nonmaxsuppts(hdet, radius=5, N=1)
@test abs(r[1]-r0) < tol && abs(c[1]-c0) < tol

(S,So) = fastradial(img, [1, 3, 5])
(r,c) = nonmaxsuppts(S, radius=5, thresh=maximum(S)/2);
@test abs(r[1]-r0) < tol && abs(c[1]-c0) < tol

# Check nonmaxsuppts does subpixel precision. Create a 2x2 block of
# pixels in the image and locate the centre of that
img[r0:r0+1, c0:c0+1] = 10
(cimg, r, c) = shi_tomasi(img, N=1,radius=5,subpixel=true)
tol = 0.1
@test abs(r[1]-(r0+0.5)) < tol && abs(c[1]-(c0+0.5)) < tol

ext: jl


text:
abstract type AbstractQFunction end

struct QFunctionNone <: AbstractQFunction end
Base.getindex(::QFunctionNone) = C.CEED_QFUNCTION_NONE[]

"""
    QFunction

A libCEED `CeedQFunction` object, typically created using the [`@interior_qf`](@ref) macro.

A `QFunction` can also be created from the "Q-function gallery" using
[`create_interior_qfunction`](@ref). The identity Q-function can be created using
[`create_identity_qfunction`](@ref).
"""
mutable struct QFunction <: AbstractQFunction
    ref::RefValue{C.CeedQFunction}
    user_qf::Union{Nothing,UserQFunction}
    ctx::Union{Nothing,Context}
    function QFunction(ref, user_qf)
        obj = new(ref, user_qf, nothing)
        finalizer(obj) do x
            # ccall(:jl_safe_printf, Cvoid, (Cstring, Cstring), "Finalizing %s.\n", repr(x))
            destroy(x)
        end
        return obj
    end
end
QFunction(ref::Ref{C.CeedQFunction}) = QFunction(ref, nothing)
destroy(qf::QFunction) = C.CeedQFunctionDestroy(qf.ref) # COV_EXCL_LINE
Base.getindex(qf::QFunction) = qf.ref[]
Base.show(io::IO, ::MIME"text/plain", qf::QFunction) =
    ceed_show(io, qf, C.CeedQFunctionView)

function create_interior_qfunction(c::Ceed, f::UserQFunction; vlength=1)
    ref = Ref{C.CeedQFunction}()
    # Use empty string as source location to indicate to libCEED that there is
    # no C source for this Q-function
    C.CeedQFunctionCreateInterior(c[], vlength, f.fptr, "", ref)
    # COV_EXCL_START
    if !isnothing(f.cuf)
        C.CeedQFunctionSetCUDAUserFunction(ref[], f.cuf)
    elseif iscuda(c) && !isdefined(@__MODULE__, :CUDA)
        error(
            string(
                "In order to use user Q-functions with a CUDA backend, the CUDA.jl package ",
                "must be loaded",
            ),
        )
    end
    # COV_EXCL_STOP
    QFunction(ref, f)
end

"""
    create_interior_qfunction(ceed::Ceed, name::AbstractString)

Create a [`QFunction`](@ref) from the Q-function gallery, using the provided name.

# Examples

- Build and apply the 3D mass operator
```
build_mass_qf = create_interior_qfunction(c, "Mass3DBuild")
apply_mass_qf = create_interior_qfunction(c, "MassApply")
```
- Build and apply the 3D Poisson operator
```
build_poi_qf = create_interior_qfunction(c, "Poisson3DBuild")
apply_poi_qf = create_interior_qfunction(c, "Poisson3DApply")
```
"""
function create_interior_qfunction(c::Ceed, name::AbstractString)
    ref = Ref{C.CeedQFunction}()
    C.CeedQFunctionCreateInteriorByName(c.ref[], name, ref)
    QFunction(ref)
end

"""
    create_identity_qfunction(c::Ceed, size, inmode::EvalMode, outmode::EvalMode)

Create an identity [`QFunction`](@ref). Inputs are written into outputs in the order given.
This is useful for [`Operators`](@ref Operator) that can be represented with only the action
of a [`ElemRestriction`](@ref) and [`Basis`](@ref), such as restriction and prolongation
operators for p-multigrid. Backends may optimize `CeedOperators` with this Q-function to
avoid the copy of input data to output fields by using the same memory location for both.
"""
function create_identity_qfunction(c::Ceed, size, inmode::EvalMode, outmode::EvalMode)
    ref = Ref{C.CeedQFunction}()
    C.CeedQFunctionCreateIdentity(c[], size, inmode, outmode, ref)
    QFunction(ref)
end

function add_input!(qf::AbstractQFunction, name::AbstractString, size, emode)
    C.CeedQFunctionAddInput(qf[], name, size, emode)
end

function add_output!(qf::AbstractQFunction, name::AbstractString, size, emode)
    C.CeedQFunctionAddOutput(qf[], name, size, emode)
end

"""
    set_context!(qf::QFunction, ctx::Context)

Associate a [`Context`](@ref) object `ctx` with the given Q-function `qf`.
"""
function set_context!(qf::QFunction, ctx)
    # Preserve the context data from the GC by storing a reference
    qf.ctx = ctx
    C.CeedQFunctionSetContext(qf[], ctx[])
end

function get_field_sizes(qf::AbstractQFunction)
    ninputs = Ref{CeedInt}()
    noutputs = Ref{CeedInt}()

    C.CeedQFunctionGetNumArgs(qf[], ninputs, noutputs)

    inputs = Ref{Ptr{C.CeedQFunctionField}}()
    outputs = Ref{Ptr{C.CeedQFunctionField}}()
    C.CeedQFunctionGetFields(qf[], inputs, outputs)

    input_sizes = zeros(CeedInt, ninputs[])
    output_sizes = zeros(CeedInt, noutputs[])

    for i = 1:ninputs[]
        field = unsafe_load(inputs[], i)
        C.CeedQFunctionFieldGetSize(field, pointer(input_sizes, i))
    end

    for i = 1:noutputs[]
        field = unsafe_load(outputs[], i)
        C.CeedQFunctionFieldGetSize(field, pointer(output_sizes, i))
    end

    input_sizes, output_sizes
end

"""
    apply!(qf::QFunction, Q, vin, vout)

Apply the action of a [`QFunction`](@ref) to an array of input vectors, and store the result
in an array of output vectors.
"""
function apply!(qf::QFunction, Q, vin, vout)
    vins = map(x -> x[], vin)
    vouts = map(x -> x[], vout)
    GC.@preserve vin vout begin
        C.CeedQFunctionApply(qf[], Q, pointer(vins), pointer(vouts))
    end
end

max_stars_repo_name: jedbrown/libCEED

max_issues_repo_name: jedbrown/libCEED

max_forks_repo_name: jedbrown/libCEED


text:
"
vrr_ss(aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.

This is an auto generated routine specific to when the A shell
has s-type angular momentum, and the C shell has s-type
angular momentum.
"
function vrr_ss(aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax = 1
    vrrs = zeros(Float64,1,1,mmax)

    P = gaussian_product_center(aexpn,A,bexpn,B)
    Q = gaussian_product_center(cexpn,C,dexpn,D)
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    W = gaussian_product_center(zeta,P,eta,Q)
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2/ze
    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)
    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    vrrs[1,1,1] = Kab*Kcd*Fgamma(0,T)/sqrt(ze)

    return vrrs[:,:,1]
end
"
vrr_ps(aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.

This is an auto generated routine specific to when the A shell
has p-type angular momentum, and the C shell has s-type
angular momentum.
"
function vrr_ps(aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax = 2
    vrrs = zeros(Float64,4,1,mmax)

    P = gaussian_product_center(aexpn,A,bexpn,B)
    Q = gaussian_product_center(cexpn,C,dexpn,D)
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    W = gaussian_product_center(zeta,P,eta,Q)
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2/ze
    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)
    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    vrrs[1,1,1] = Kab*Kcd*Fgamma(0,T)/sqrt(ze)
    vrrs[1,1,2] = Kab*Kcd*Fgamma(1,T)/sqrt(ze)

    vrrs[2,1,1] = PA[1]*vrrs[1,1,1] + WP[1]*vrrs[1,1,2]
    vrrs[3,1,1] = PA[2]*vrrs[1,1,1] + WP[2]*vrrs[1,1,2]
    vrrs[4,1,1] = PA[3]*vrrs[1,1,1] + WP[3]*vrrs[1,1,2]

    return vrrs[:,:,1]
end
"
vrr_sp(aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.

This is an auto generated routine specific to when the A shell
has s-type angular momentum, and the C shell has p-type
angular momentum.
"
function vrr_sp(aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax = 2
    vrrs = zeros(Float64,1,4,mmax)

    P = gaussian_product_center(aexpn,A,bexpn,B)
    Q = gaussian_product_center(cexpn,C,dexpn,D)
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    W = gaussian_product_center(zeta,P,eta,Q)
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2/ze
    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)
    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    vrrs[1,1,1] = Kab*Kcd*Fgamma(0,T)/sqrt(ze)
    vrrs[1,1,2] = Kab*Kcd*Fgamma(1,T)/sqrt(ze)

    vrrs[1,2,1] = QC[1]*vrrs[1,1,1] + WQ[1]*vrrs[1,1,2]
    vrrs[1,3,1] = QC[2]*vrrs[1,1,1] + WQ[2]*vrrs[1,1,2]
    vrrs[1,4,1] = QC[3]*vrrs[1,1,1] + WQ[3]*vrrs[1,1,2]

    return vrrs[:,:,1]
end
"
vrr_pp(aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.

This is an auto generated routine specific to when the A shell
has p-type angular momentum, and the C shell has p-type
angular momentum.
"
function vrr_pp(aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax = 3
    vrrs = zeros(Float64,4,4,mmax)

    P = gaussian_product_center(aexpn,A,bexpn,B)
    Q = gaussian_product_center(cexpn,C,dexpn,D)
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    W = gaussian_product_center(zeta,P,eta,Q)
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2/ze
    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)
    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    vrrs[1,1,1] = Kab*Kcd*Fgamma(0,T)/sqrt(ze)
    vrrs[1,1,2] = Kab*Kcd*Fgamma(1,T)/sqrt(ze)
    vrrs[1,1,3] = Kab*Kcd*Fgamma(2,T)/sqrt(ze)

    vrrs[2,1,1] = PA[1]*vrrs[1,1,1] + WP[1]*vrrs[1,1,2]
    vrrs[2,1,2] = PA[1]*vrrs[1,1,2] + WP[1]*vrrs[1,1,3]
    vrrs[3,1,1] = PA[2]*vrrs[1,1,1] + WP[2]*vrrs[1,1,2]
    vrrs[3,1,2] = PA[2]*vrrs[1,1,2] + WP[2]*vrrs[1,1,3]
    vrrs[4,1,1] = PA[3]*vrrs[1,1,1] + WP[3]*vrrs[1,1,2]
    vrrs[4,1,2] = PA[3]*vrrs[1,1,2] + WP[3]*vrrs[1,1,3]
    vrrs[1,2,1] = QC[1]*vrrs[1,1,1] + WQ[1]*vrrs[1,1,2]
    vrrs[1,2,2] = QC[1]*vrrs[1,1,2] + WQ[1]*vrrs[1,1,3]
    vrrs[2,2,1] = QC[1]*vrrs[2,1,1] + WQ[1]*vrrs[2,1,2] + 1/2ze*vrrs[1,1,2]
    vrrs[3,2,1] = QC[1]*vrrs[3,1,1] + WQ[1]*vrrs[3,1,2]
    vrrs[4,2,1] = QC[1]*vrrs[4,1,1] + WQ[1]*vrrs[4,1,2]
    vrrs[1,3,1] = QC[2]*vrrs[1,1,1] + WQ[2]*vrrs[1,1,2]
    vrrs[1,3,2] = QC[2]*vrrs[1,1,2] + WQ[2]*vrrs[1,1,3]
    vrrs[2,3,1] = QC[2]*vrrs[2,1,1] + WQ[2]*vrrs[2,1,2]
    vrrs[3,3,1] = QC[2]*vrrs[3,1,1] + WQ[2]*vrrs[3,1,2] + 1/2ze*vrrs[1,1,2]
    vrrs[4,3,1] = QC[2]*vrrs[4,1,1] + WQ[2]*vrrs[4,1,2]
    vrrs[1,4,1] = QC[3]*vrrs[1,1,1] + WQ[3]*vrrs[1,1,2]
    vrrs[1,4,2] = QC[3]*vrrs[1,1,2] + WQ[3]*vrrs[1,1,3]
    vrrs[2,4,1] = QC[3]*vrrs[2,1,1] + WQ[3]*vrrs[2,1,2]
    vrrs[3,4,1] = QC[3]*vrrs[3,1,1] + WQ[3]*vrrs[3,1,2]
    vrrs[4,4,1] = QC[3]*vrrs[4,1,1] + WQ[3]*vrrs[4,1,2] + 1/2ze*vrrs[1,1,2]

    return vrrs[:,:,1]
end
"
vrr_ds(aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.

This is an auto generated routine specific to when the A shell
has d-type angular momentum, and the C shell has s-type
angular momentum.
"
function vrr_ds(aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax = 3
    vrrs = zeros(Float64,10,1,mmax)

    P = gaussian_product_center(aexpn,A,bexpn,B)
    Q = gaussian_product_center(cexpn,C,dexpn,D)
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    W = gaussian_product_center(zeta,P,eta,Q)
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2/ze
    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)
    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    vrrs[1,1,1] = Kab*Kcd*Fgamma(0,T)/sqrt(ze)
    vrrs[1,1,2] = Kab*Kcd*Fgamma(1,T)/sqrt(ze)
    vrrs[1,1,3] = Kab*Kcd*Fgamma(2,T)/sqrt(ze)

    vrrs[2,1,1] = PA[1]*vrrs[1,1,1] + WP[1]*vrrs[1,1,2]
    vrrs[2,1,2] = PA[1]*vrrs[1,1,2] + WP[1]*vrrs[1,1,3]
    vrrs[3,1,1] = PA[2]*vrrs[1,1,1] + WP[2]*vrrs[1,1,2]
    vrrs[3,1,2] = PA[2]*vrrs[1,1,2] + WP[2]*vrrs[1,1,3]
    vrrs[4,1,1] = PA[3]*vrrs[1,1,1] + WP[3]*vrrs[1,1,2]
    vrrs[4,1,2] = PA[3]*vrrs[1,1,2] + WP[3]*vrrs[1,1,3]
    vrrs[5,1,1] = PA[1]*vrrs[2,1,1] + WP[1]*vrrs[2,1,2] + 1/2zeta*(vrrs[1,1,1]-eta/ze*vrrs[1,1,2])
    vrrs[6,1,1] = PA[1]*vrrs[3,1,1] + WP[1]*vrrs[3,1,2]
    vrrs[7,1,1] = PA[1]*vrrs[4,1,1] + WP[1]*vrrs[4,1,2]
    vrrs[8,1,1] = PA[2]*vrrs[3,1,1] + WP[2]*vrrs[3,1,2] + 1/2zeta*(vrrs[1,1,1]-eta/ze*vrrs[1,1,2])
    vrrs[9,1,1] = PA[2]*vrrs[4,1,1] + WP[2]*vrrs[4,1,2]
    vrrs[10,1,1] = PA[3]*vrrs[4,1,1] + WP[3]*vrrs[4,1,2] + 1/2zeta*(vrrs[1,1,1]-eta/ze*vrrs[1,1,2])

    return vrrs[:,:,1]
end
"
vrr_sd(aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.

This is an auto generated routine specific to when the A shell
has s-type angular momentum, and the C shell has d-type
angular momentum.
"
function vrr_sd(aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax = 3
    vrrs = zeros(Float64,1,10,mmax)

    P = gaussian_product_center(aexpn,A,bexpn,B)
    Q = gaussian_product_center(cexpn,C,dexpn,D)
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    W = gaussian_product_center(zeta,P,eta,Q)
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2/ze
    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)
    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    vrrs[1,1,1] = Kab*Kcd*Fgamma(0,T)/sqrt(ze)
    vrrs[1,1,2] = Kab*Kcd*Fgamma(1,T)/sqrt(ze)
    vrrs[1,1,3] = Kab*Kcd*Fgamma(2,T)/sqrt(ze)

    vrrs[1,2,1] = QC[1]*vrrs[1,1,1] + WQ[1]*vrrs[1,1,2]
    vrrs[1,2,2] = QC[1]*vrrs[1,1,2] + WQ[1]*vrrs[1,1,3]
    vrrs[1,3,1] = QC[2]*vrrs[1,1,1] + WQ[2]*vrrs[1,1,2]
    vrrs[1,3,2] = QC[2]*vrrs[1,1,2] + WQ[2]*vrrs[1,1,3]
    vrrs[1,4,1] = QC[3]*vrrs[1,1,1] + WQ[3]*vrrs[1,1,2]
    vrrs[1,4,2] = QC[3]*vrrs[1,1,2] + WQ[3]*vrrs[1,1,3]
    vrrs[1,5,1] = QC[1]*vrrs[1,2,1] + WQ[1]*vrrs[1,2,2] + 1/2eta*(vrrs[1,1,1]-zeta/ze*vrrs[1,1,2])
    vrrs[1,6,1] = QC[1]*vrrs[1,3,1] + WQ[1]*vrrs[1,3,2]
    vrrs[1,7,1] = QC[1]*vrrs[1,4,1] + WQ[1]*vrrs[1,4,2]
    vrrs[1,8,1] = QC[2]*vrrs[1,3,1] + WQ[2]*vrrs[1,3,2] + 1/2eta*(vrrs[1,1,1]-zeta/ze*vrrs[1,1,2])
    vrrs[1,9,1] = QC[2]*vrrs[1,4,1] + WQ[2]*vrrs[1,4,2]
    vrrs[1,10,1] = QC[3]*vrrs[1,4,1] + WQ[3]*vrrs[1,4,2] + 1/2eta*(vrrs[1,1,1]-zeta/ze*vrrs[1,1,2])

    return vrrs[:,:,1]
end
"
vrr_dp(aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.

This is an auto generated routine specific to when the A shell
has d-type angular momentum, and the C shell has p-type
angular momentum.
"
function vrr_dp(aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax = 4
    vrrs = zeros(Float64,10,4,mmax)

    P = gaussian_product_center(aexpn,A,bexpn,B)
    Q = gaussian_product_center(cexpn,C,dexpn,D)
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    W = gaussian_product_center(zeta,P,eta,Q)
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2/ze
    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)
    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    vrrs[1,1,1] = Kab*Kcd*Fgamma(0,T)/sqrt(ze)
    vrrs[1,1,2] = Kab*Kcd*Fgamma(1,T)/sqrt(ze)
    vrrs[1,1,3] = Kab*Kcd*Fgamma(2,T)/sqrt(ze)
    vrrs[1,1,4] = Kab*Kcd*Fgamma(3,T)/sqrt(ze)

    vrrs[2,1,1] = PA[1]*vrrs[1,1,1] + WP[1]*vrrs[1,1,2]
    vrrs[2,1,2] = PA[1]*vrrs[1,1,2] + WP[1]*vrrs[1,1,3]
    vrrs[2,1,3] = PA[1]*vrrs[1,1,3] + WP[1]*vrrs[1,1,4]
    vrrs[3,1,1] = PA[2]*vrrs[1,1,1] + WP[2]*vrrs[1,1,2]
    vrrs[3,1,2] = PA[2]*vrrs[1,1,2] + WP[2]*vrrs[1,1,3]
    vrrs[3,1,3] = PA[2]*vrrs[1,1,3] + WP[2]*vrrs[1,1,4]
    vrrs[4,1,1] = PA[3]*vrrs[1,1,1] + WP[3]*vrrs[1,1,2]
    vrrs[4,1,2] = PA[3]*vrrs[1,1,2] + WP[3]*vrrs[1,1,3]
    vrrs[4,1,3] = PA[3]*vrrs[1,1,3] + WP[3]*vrrs[1,1,4]
    vrrs[5,1,1] = PA[1]*vrrs[2,1,1] + WP[1]*vrrs[2,1,2] + 1/2zeta*(vrrs[1,1,1]-eta/ze*vrrs[1,1,2])
    vrrs[5,1,2] = PA[1]*vrrs[2,1,2] + WP[1]*vrrs[2,1,3] + 1/2zeta*(vrrs[1,1,2]-eta/ze*vrrs[1,1,3])
    vrrs[6,1,1] = PA[1]*vrrs[3,1,1] + WP[1]*vrrs[3,1,2]
    vrrs[6,1,2] = PA[1]*vrrs[3,1,2] + WP[1]*vrrs[3,1,3]
    vrrs[7,1,1] = PA[1]*vrrs[4,1,1] + WP[1]*vrrs[4,1,2]
    vrrs[7,1,2] = PA[1]*vrrs[4,1,2] + WP[1]*vrrs[4,1,3]
    vrrs[8,1,1] = PA[2]*vrrs[3,1,1] + WP[2]*vrrs[3,1,2] + 1/2zeta*(vrrs[1,1,1]-eta/ze*vrrs[1,1,2])
    vrrs[8,1,2] = PA[2]*vrrs[3,1,2] + WP[2]*vrrs[3,1,3] + 1/2zeta*(vrrs[1,1,2]-eta/ze*vrrs[1,1,3])
    vrrs[9,1,1] = PA[2]*vrrs[4,1,1] + WP[2]*vrrs[4,1,2]
    vrrs[9,1,2] = PA[2]*vrrs[4,1,2] + WP[2]*vrrs[4,1,3]
    vrrs[10,1,1] = PA[3]*vrrs[4,1,1] + WP[3]*vrrs[4,1,2] + 1/2zeta*(vrrs[1,1,1]-eta/ze*vrrs[1,1,2])
    vrrs[10,1,2] = PA[3]*vrrs[4,1,2] + WP[3]*vrrs[4,1,3] + 1/2zeta*(vrrs[1,1,2]-eta/ze*vrrs[1,1,3])
    vrrs[1,2,1] = QC[1]*vrrs[1,1,1] + WQ[1]*vrrs[1,1,2]
    vrrs[1,2,2] = QC[1]*vrrs[1,1,2] + WQ[1]*vrrs[1,1,3]
    vrrs[1,2,3] = QC[1]*vrrs[1,1,3] + WQ[1]*vrrs[1,1,4]
    vrrs[2,2,1] = QC[1]*vrrs[2,1,1] + WQ[1]*vrrs[2,1,2] + 1/2ze*vrrs[1,1,2]
    vrrs[2,2,2] = QC[1]*vrrs[2,1,2] + WQ[1]*vrrs[2,1,3] + 1/2ze*vrrs[1,1,3]
    vrrs[3,2,1] = QC[1]*vrrs[3,1,1] + WQ[1]*vrrs[3,1,2]
    vrrs[3,2,2] = QC[1]*vrrs[3,1,2] + WQ[1]*vrrs[3,1,3]
    vrrs[4,2,1] = QC[1]*vrrs[4,1,1] + WQ[1]*vrrs[4,1,2]
    vrrs[4,2,2] = QC[1]*vrrs[4,1,2] + WQ[1]*vrrs[4,1,3]
    vrrs[5,2,1] = QC[1]*vrrs[5,1,1] + WQ[1]*vrrs[5,1,2] + 2/2ze*vrrs[2,1,2]
    vrrs[6,2,1] = QC[1]*vrrs[6,1,1] + WQ[1]*vrrs[6,1,2] + 1/2ze*vrrs[3,1,2]
    vrrs[7,2,1] = QC[1]*vrrs[7,1,1] + WQ[1]*vrrs[7,1,2] + 1/2ze*vrrs[4,1,2]
    vrrs[8,2,1] = QC[1]*vrrs[8,1,1] + WQ[1]*vrrs[8,1,2]
    vrrs[9,2,1] = QC[1]*vrrs[9,1,1] + WQ[1]*vrrs[9,1,2]
    vrrs[10,2,1] = QC[1]*vrrs[10,1,1] + WQ[1]*vrrs[10,1,2]
    vrrs[1,3,1] = QC[2]*vrrs[1,1,1] + WQ[2]*vrrs[1,1,2]
    vrrs[1,3,2] = QC[2]*vrrs[1,1,2] + WQ[2]*vrrs[1,1,3]
    vrrs[1,3,3] = QC[2]*vrrs[1,1,3] + WQ[2]*vrrs[1,1,4]
    vrrs[2,3,1] = QC[2]*vrrs[2,1,1] + WQ[2]*vrrs[2,1,2]
    vrrs[2,3,2] = QC[2]*vrrs[2,1,2] + WQ[2]*vrrs[2,1,3]
    vrrs[3,3,1] = QC[2]*vrrs[3,1,1] + WQ[2]*vrrs[3,1,2] + 1/2ze*vrrs[1,1,2]
    vrrs[3,3,2] = QC[2]*vrrs[3,1,2] + WQ[2]*vrrs[3,1,3] + 1/2ze*vrrs[1,1,3]
    vrrs[4,3,1] = QC[2]*vrrs[4,1,1] + WQ[2]*vrrs[4,1,2]
    vrrs[4,3,2] = QC[2]*vrrs[4,1,2] + WQ[2]*vrrs[4,1,3]
    vrrs[5,3,1] = QC[2]*vrrs[5,1,1] + WQ[2]*vrrs[5,1,2]
    vrrs[6,3,1] = QC[2]*vrrs[6,1,1] + WQ[2]*vrrs[6,1,2] + 1/2ze*vrrs[2,1,2]
    vrrs[7,3,1] = QC[2]*vrrs[7,1,1] + WQ[2]*vrrs[7,1,2]
    vrrs[8,3,1] = QC[2]*vrrs[8,1,1] + WQ[2]*vrrs[8,1,2] + 2/2ze*vrrs[3,1,2]
    vrrs[9,3,1] = QC[2]*vrrs[9,1,1] + WQ[2]*vrrs[9,1,2] + 1/2ze*vrrs[4,1,2]
    vrrs[10,3,1] = QC[2]*vrrs[10,1,1] + WQ[2]*vrrs[10,1,2]
    vrrs[1,4,1] = QC[3]*vrrs[1,1,1] + WQ[3]*vrrs[1,1,2]
    vrrs[1,4,2] = QC[3]*vrrs[1,1,2] + WQ[3]*vrrs[1,1,3]
    vrrs[1,4,3] = QC[3]*vrrs[1,1,3] + WQ[3]*vrrs[1,1,4]
    vrrs[2,4,1] = QC[3]*vrrs[2,1,1] + WQ[3]*vrrs[2,1,2]
    vrrs[2,4,2] = QC[3]*vrrs[2,1,2] + WQ[3]*vrrs[2,1,3]
    vrrs[3,4,1] = QC[3]*vrrs[3,1,1] + WQ[3]*vrrs[3,1,2]
    vrrs[3,4,2] = QC[3]*vrrs[3,1,2] + WQ[3]*vrrs[3,1,3]
    vrrs[4,4,1] = QC[3]*vrrs[4,1,1] + WQ[3]*vrrs[4,1,2] + 1/2ze*vrrs[1,1,2]
    vrrs[4,4,2] = QC[3]*vrrs[4,1,2] + WQ[3]*vrrs[4,1,3] + 1/2ze*vrrs[1,1,3]
    vrrs[5,4,1] = QC[3]*vrrs[5,1,1] + WQ[3]*vrrs[5,1,2]
    vrrs[6,4,1] = QC[3]*vrrs[6,1,1] + WQ[3]*vrrs[6,1,2]
    vrrs[7,4,1] = QC[3]*vrrs[7,1,1] + WQ[3]*vrrs[7,1,2] + 1/2ze*vrrs[2,1,2]
    vrrs[8,4,1] = QC[3]*vrrs[8,1,1] + WQ[3]*vrrs[8,1,2]
    vrrs[9,4,1] = QC[3]*vrrs[9,1,1] + WQ[3]*vrrs[9,1,2] + 1/2ze*vrrs[3,1,2]
    vrrs[10,4,1] = QC[3]*vrrs[10,1,1] + WQ[3]*vrrs[10,1,2] + 2/2ze*vrrs[4,1,2]

    return vrrs[:,:,1]
end
"
vrr_pd(aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.

This is an auto generated routine specific to when the A shell
has p-type angular momentum, and the C shell has d-type
angular momentum.
"
function vrr_pd(aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax = 4
    vrrs = zeros(Float64,4,10,mmax)

    P = gaussian_product_center(aexpn,A,bexpn,B)
    Q = gaussian_product_center(cexpn,C,dexpn,D)
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    W = gaussian_product_center(zeta,P,eta,Q)
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2/ze
    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)
    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    vrrs[1,1,1] = Kab*Kcd*Fgamma(0,T)/sqrt(ze)
    vrrs[1,1,2] = Kab*Kcd*Fgamma(1,T)/sqrt(ze)
    vrrs[1,1,3] = Kab*Kcd*Fgamma(2,T)/sqrt(ze)
    vrrs[1,1,4] = Kab*Kcd*Fgamma(3,T)/sqrt(ze)

    vrrs[2,1,1] = PA[1]*vrrs[1,1,1] + WP[1]*vrrs[1,1,2]
    vrrs[2,1,2] = PA[1]*vrrs[1,1,2] + WP[1]*vrrs[1,1,3]
    vrrs[2,1,3] = PA[1]*vrrs[1,1,3] + WP[1]*vrrs[1,1,4]
    vrrs[3,1,1] = PA[2]*vrrs[1,1,1] + WP[2]*vrrs[1,1,2]
    vrrs[3,1,2] = PA[2]*vrrs[1,1,2] + WP[2]*vrrs[1,1,3]
    vrrs[3,1,3] = PA[2]*vrrs[1,1,3] + WP[2]*vrrs[1,1,4]
    vrrs[4,1,1] = PA[3]*vrrs[1,1,1] + WP[3]*vrrs[1,1,2]
    vrrs[4,1,2] = PA[3]*vrrs[1,1,2] + WP[3]*vrrs[1,1,3]
    vrrs[4,1,3] = PA[3]*vrrs[1,1,3] + WP[3]*vrrs[1,1,4]
    vrrs[1,2,1] = QC[1]*vrrs[1,1,1] + WQ[1]*vrrs[1,1,2]
    vrrs[1,2,2] = QC[1]*vrrs[1,1,2] + WQ[1]*vrrs[1,1,3]
    vrrs[1,2,3] = QC[1]*vrrs[1,1,3] + WQ[1]*vrrs[1,1,4]
    vrrs[2,2,1] = QC[1]*vrrs[2,1,1] + WQ[1]*vrrs[2,1,2] + 1/2ze*vrrs[1,1,2]
    vrrs[2,2,2] = QC[1]*vrrs[2,1,2] + WQ[1]*vrrs[2,1,3] + 1/2ze*vrrs[1,1,3]
    vrrs[3,2,1] = QC[1]*vrrs[3,1,1] + WQ[1]*vrrs[3,1,2]
    vrrs[3,2,2] = QC[1]*vrrs[3,1,2] + WQ[1]*vrrs[3,1,3]
    vrrs[4,2,1] = QC[1]*vrrs[4,1,1] + WQ[1]*vrrs[4,1,2]
    vrrs[4,2,2] = QC[1]*vrrs[4,1,2] + WQ[1]*vrrs[4,1,3]
    vrrs[1,3,1] = QC[2]*vrrs[1,1,1] + WQ[2]*vrrs[1,1,2]
    vrrs[1,3,2] = QC[2]*vrrs[1,1,2] + WQ[2]*vrrs[1,1,3]
    vrrs[1,3,3] = QC[2]*vrrs[1,1,3] + WQ[2]*vrrs[1,1,4]
    vrrs[2,3,1] = QC[2]*vrrs[2,1,1] + WQ[2]*vrrs[2,1,2]
    vrrs[2,3,2] = QC[2]*vrrs[2,1,2] + WQ[2]*vrrs[2,1,3]
    vrrs[3,3,1] = QC[2]*vrrs[3,1,1] + WQ[2]*vrrs[3,1,2] + 1/2ze*vrrs[1,1,2]
    vrrs[3,3,2] = QC[2]*vrrs[3,1,2] + WQ[2]*vrrs[3,1,3] + 1/2ze*vrrs[1,1,3]
    vrrs[4,3,1] = QC[2]*vrrs[4,1,1] + WQ[2]*vrrs[4,1,2]
    vrrs[4,3,2] = QC[2]*vrrs[4,1,2] + WQ[2]*vrrs[4,1,3]
    vrrs[1,4,1] = QC[3]*vrrs[1,1,1] + WQ[3]*vrrs[1,1,2]
    vrrs[1,4,2] = QC[3]*vrrs[1,1,2] + WQ[3]*vrrs[1,1,3]
    vrrs[1,4,3] = QC[3]*vrrs[1,1,3] + WQ[3]*vrrs[1,1,4]
    vrrs[2,4,1] = QC[3]*vrrs[2,1,1] + WQ[3]*vrrs[2,1,2]
    vrrs[2,4,2] = QC[3]*vrrs[2,1,2] + WQ[3]*vrrs[2,1,3]
    vrrs[3,4,1] = QC[3]*vrrs[3,1,1] + WQ[3]*vrrs[3,1,2]
    vrrs[3,4,2] = QC[3]*vrrs[3,1,2] + WQ[3]*vrrs[3,1,3]
    vrrs[4,4,1] = QC[3]*vrrs[4,1,1] + WQ[3]*vrrs[4,1,2] + 1/2ze*vrrs[1,1,2]
    vrrs[4,4,2] = QC[3]*vrrs[4,1,2] + WQ[3]*vrrs[4,1,3] + 1/2ze*vrrs[1,1,3]
    vrrs[1,5,1] = QC[1]*vrrs[1,2,1] + WQ[1]*vrrs[1,2,2] + 1/2eta*(vrrs[1,1,1]-zeta/ze*vrrs[1,1,2])
    vrrs[1,5,2] = QC[1]*vrrs[1,2,2] + WQ[1]*vrrs[1,2,3] + 1/2eta*(vrrs[1,1,2]-zeta/ze*vrrs[1,1,3])
    vrrs[2,5,1] = QC[1]*vrrs[2,2,1] + WQ[1]*vrrs[2,2,2] + 1/2eta*(vrrs[2,1,1]-zeta/ze*vrrs[2,1,2]) + 1/2ze*vrrs[1,2,2]
    vrrs[3,5,1] = QC[1]*vrrs[3,2,1] + WQ[1]*vrrs[3,2,2] + 1/2eta*(vrrs[3,1,1]-zeta/ze*vrrs[3,1,2])
    vrrs[4,5,1] = QC[1]*vrrs[4,2,1] + WQ[1]*vrrs[4,2,2] + 1/2eta*(vrrs[4,1,1]-zeta/ze*vrrs[4,1,2])
    vrrs[1,6,1] = QC[1]*vrrs[1,3,1] + WQ[1]*vrrs[1,3,2]
    vrrs[1,6,2] = QC[1]*vrrs[1,3,2] + WQ[1]*vrrs[1,3,3]
    vrrs[2,6,1] = QC[1]*vrrs[2,3,1] + WQ[1]*vrrs[2,3,2] + 1/2ze*vrrs[1,3,2]
    vrrs[3,6,1] = QC[1]*vrrs[3,3,1] + WQ[1]*vrrs[3,3,2]
    vrrs[4,6,1] = QC[1]*vrrs[4,3,1] + WQ[1]*vrrs[4,3,2]
    vrrs[1,7,1] = QC[1]*vrrs[1,4,1] + WQ[1]*vrrs[1,4,2]
    vrrs[1,7,2] = QC[1]*vrrs[1,4,2] + WQ[1]*vrrs[1,4,3]
    vrrs[2,7,1] = QC[1]*vrrs[2,4,1] + WQ[1]*vrrs[2,4,2] + 1/2ze*vrrs[1,4,2]
    vrrs[3,7,1] = QC[1]*vrrs[3,4,1] + WQ[1]*vrrs[3,4,2]
    vrrs[4,7,1] = QC[1]*vrrs[4,4,1] + WQ[1]*vrrs[4,4,2]
    vrrs[1,8,1] = QC[2]*vrrs[1,3,1] + WQ[2]*vrrs[1,3,2] + 1/2eta*(vrrs[1,1,1]-zeta/ze*vrrs[1,1,2])
    vrrs[1,8,2] = QC[2]*vrrs[1,3,2] + WQ[2]*vrrs[1,3,3] + 1/2eta*(vrrs[1,1,2]-zeta/ze*vrrs[1,1,3])
    vrrs[2,8,1] = QC[2]*vrrs[2,3,1] + WQ[2]*vrrs[2,3,2] + 1/2eta*(vrrs[2,1,1]-zeta/ze*vrrs[2,1,2])
    vrrs[3,8,1] = QC[2]*vrrs[3,3,1] + WQ[2]*vrrs[3,3,2] + 1/2eta*(vrrs[3,1,1]-zeta/ze*vrrs[3,1,2]) + 1/2ze*vrrs[1,3,2]
    vrrs[4,8,1] = QC[2]*vrrs[4,3,1] + WQ[2]*vrrs[4,3,2] + 1/2eta*(vrrs[4,1,1]-zeta/ze*vrrs[4,1,2])
    vrrs[1,9,1] = QC[2]*vrrs[1,4,1] + WQ[2]*vrrs[1,4,2]
    vrrs[1,9,2] = QC[2]*vrrs[1,4,2] + WQ[2]*vrrs[1,4,3]
    vrrs[2,9,1] = QC[2]*vrrs[2,4,1] + WQ[2]*vrrs[2,4,2]
    vrrs[3,9,1] = QC[2]*vrrs[3,4,1] + WQ[2]*vrrs[3,4,2] + 1/2ze*vrrs[1,4,2]
    vrrs[4,9,1] = QC[2]*vrrs[4,4,1] + WQ[2]*vrrs[4,4,2]
    vrrs[1,10,1] = QC[3]*vrrs[1,4,1] + WQ[3]*vrrs[1,4,2] + 1/2eta*(vrrs[1,1,1]-zeta/ze*vrrs[1,1,2])
    vrrs[1,10,2] = QC[3]*vrrs[1,4,2] + WQ[3]*vrrs[1,4,3] + 1/2eta*(vrrs[1,1,2]-zeta/ze*vrrs[1,1,3])
    vrrs[2,10,1] = QC[3]*vrrs[2,4,1] + WQ[3]*vrrs[2,4,2] + 1/2eta*(vrrs[2,1,1]-zeta/ze*vrrs[2,1,2])
    vrrs[3,10,1] = QC[3]*vrrs[3,4,1] + WQ[3]*vrrs[3,4,2] + 1/2eta*(vrrs[3,1,1]-zeta/ze*vrrs[3,1,2])
    vrrs[4,10,1] = QC[3]*vrrs[4,4,1] + WQ[3]*vrrs[4,4,2] + 1/2eta*(vrrs[4,1,1]-zeta/ze*vrrs[4,1,2]) + 1/2ze*vrrs[1,4,2]

    return vrrs[:,:,1]
end
"
vrr_dd(aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.

This is an auto generated routine specific to when the A shell
has d-type angular momentum, and the C shell has d-type
angular momentum.
"
function vrr_dd(aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax = 5
    vrrs = zeros(Float64,10,10,mmax)

    P = gaussian_product_center(aexpn,A,bexpn,B)
    Q = gaussian_product_center(cexpn,C,dexpn,D)
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    W = gaussian_product_center(zeta,P,eta,Q)
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2/ze
    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)
    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    vrrs[1,1,1] = Kab*Kcd*Fgamma(0,T)/sqrt(ze)
    vrrs[1,1,2] = Kab*Kcd*Fgamma(1,T)/sqrt(ze)
    vrrs[1,1,3] = Kab*Kcd*Fgamma(2,T)/sqrt(ze)
    vrrs[1,1,4] = Kab*Kcd*Fgamma(3,T)/sqrt(ze)
    vrrs[1,1,5] = Kab*Kcd*Fgamma(4,T)/sqrt(ze)

    vrrs[2,1,1] = PA[1]*vrrs[1,1,1] + WP[1]*vrrs[1,1,2]
    vrrs[2,1,2] = PA[1]*vrrs[1,1,2] + WP[1]*vrrs[1,1,3]
    vrrs[2,1,3] = PA[1]*vrrs[1,1,3] + WP[1]*vrrs[1,1,4]
    vrrs[2,1,4] = PA[1]*vrrs[1,1,4] + WP[1]*vrrs[1,1,5]
    vrrs[3,1,1] = PA[2]*vrrs[1,1,1] + WP[2]*vrrs[1,1,2]
    vrrs[3,1,2] = PA[2]*vrrs[1,1,2] + WP[2]*vrrs[1,1,3]
    vrrs[3,1,3] = PA[2]*vrrs[1,1,3] + WP[2]*vrrs[1,1,4]
    vrrs[3,1,4] = PA[2]*vrrs[1,1,4] + WP[2]*vrrs[1,1,5]
    vrrs[4,1,1] = PA[3]*vrrs[1,1,1] + WP[3]*vrrs[1,1,2]
    vrrs[4,1,2] = PA[3]*vrrs[1,1,2] + WP[3]*vrrs[1,1,3]
    vrrs[4,1,3] = PA[3]*vrrs[1,1,3] + WP[3]*vrrs[1,1,4]
    vrrs[4,1,4] = PA[3]*vrrs[1,1,4] + WP[3]*vrrs[1,1,5]
    vrrs[5,1,1] = PA[1]*vrrs[2,1,1] + WP[1]*vrrs[2,1,2] + 1/2zeta*(vrrs[1,1,1]-eta/ze*vrrs[1,1,2])
    vrrs[5,1,2] = PA[1]*vrrs[2,1,2] + WP[1]*vrrs[2,1,3] + 1/2zeta*(vrrs[1,1,2]-eta/ze*vrrs[1,1,3])
    vrrs[5,1,3] = PA[1]*vrrs[2,1,3] + WP[1]*vrrs[2,1,4] + 1/2zeta*(vrrs[1,1,3]-eta/ze*vrrs[1,1,4])
    vrrs[6,1,1] = PA[1]*vrrs[3,1,1] + WP[1]*vrrs[3,1,2]
    vrrs[6,1,2] = PA[1]*vrrs[3,1,2] + WP[1]*vrrs[3,1,3]
    vrrs[6,1,3] = PA[1]*vrrs[3,1,3] + WP[1]*vrrs[3,1,4]
    vrrs[7,1,1] = PA[1]*vrrs[4,1,1] + WP[1]*vrrs[4,1,2]
    vrrs[7,1,2] = PA[1]*vrrs[4,1,2] + WP[1]*vrrs[4,1,3]
    vrrs[7,1,3] = PA[1]*vrrs[4,1,3] + WP[1]*vrrs[4,1,4]
    vrrs[8,1,1] = PA[2]*vrrs[3,1,1] + WP[2]*vrrs[3,1,2] + 1/2zeta*(vrrs[1,1,1]-eta/ze*vrrs[1,1,2])
    vrrs[8,1,2] = PA[2]*vrrs[3,1,2] + WP[2]*vrrs[3,1,3] + 1/2zeta*(vrrs[1,1,2]-eta/ze*vrrs[1,1,3])
    vrrs[8,1,3] = PA[2]*vrrs[3,1,3] + WP[2]*vrrs[3,1,4] + 1/2zeta*(vrrs[1,1,3]-eta/ze*vrrs[1,1,4])
    vrrs[9,1,1] = PA[2]*vrrs[4,1,1] + WP[2]*vrrs[4,1,2]
    vrrs[9,1,2] = PA[2]*vrrs[4,1,2] + WP[2]*vrrs[4,1,3]
    vrrs[9,1,3] = PA[2]*vrrs[4,1,3] + WP[2]*vrrs[4,1,4]
    vrrs[10,1,1] = PA[3]*vrrs[4,1,1] + WP[3]*vrrs[4,1,2] + 1/2zeta*(vrrs[1,1,1]-eta/ze*vrrs[1,1,2])
    vrrs[10,1,2] = PA[3]*vrrs[4,1,2] + WP[3]*vrrs[4,1,3] + 1/2zeta*(vrrs[1,1,2]-eta/ze*vrrs[1,1,3])
    vrrs[10,1,3] = PA[3]*vrrs[4,1,3] + WP[3]*vrrs[4,1,4] + 1/2zeta*(vrrs[1,1,3]-eta/ze*vrrs[1,1,4])
    vrrs[1,2,1] = QC[1]*vrrs[1,1,1] + WQ[1]*vrrs[1,1,2]
    vrrs[1,2,2] = QC[1]*vrrs[1,1,2] + WQ[1]*vrrs[1,1,3]
    vrrs[1,2,3] = QC[1]*vrrs[1,1,3] + WQ[1]*vrrs[1,1,4]
    vrrs[1,2,4] = QC[1]*vrrs[1,1,4] + WQ[1]*vrrs[1,1,5]
    vrrs[2,2,1] = QC[1]*vrrs[2,1,1] + WQ[1]*vrrs[2,1,2] + 1/2ze*vrrs[1,1,2]
    vrrs[2,2,2] = QC[1]*vrrs[2,1,2] + WQ[1]*vrrs[2,1,3] + 1/2ze*vrrs[1,1,3]
    vrrs[2,2,3] = QC[1]*vrrs[2,1,3] + WQ[1]*vrrs[2,1,4] + 1/2ze*vrrs[1,1,4]
    vrrs[3,2,1] = QC[1]*vrrs[3,1,1] + WQ[1]*vrrs[3,1,2]
    vrrs[3,2,2] = QC[1]*vrrs[3,1,2] + WQ[1]*vrrs[3,1,3]
    vrrs[3,2,3] = QC[1]*vrrs[3,1,3] + WQ[1]*vrrs[3,1,4]
    vrrs[4,2,1] = QC[1]*vrrs[4,1,1] + WQ[1]*vrrs[4,1,2]
    vrrs[4,2,2] = QC[1]*vrrs[4,1,2] + WQ[1]*vrrs[4,1,3]
    vrrs[4,2,3] = QC[1]*vrrs[4,1,3] + WQ[1]*vrrs[4,1,4]
    vrrs[5,2,1] = QC[1]*vrrs[5,1,1] + WQ[1]*vrrs[5,1,2] + 2/2ze*vrrs[2,1,2]
    vrrs[5,2,2] = QC[1]*vrrs[5,1,2] + WQ[1]*vrrs[5,1,3] + 2/2ze*vrrs[2,1,3]
    vrrs[6,2,1] = QC[1]*vrrs[6,1,1] + WQ[1]*vrrs[6,1,2] + 1/2ze*vrrs[3,1,2]
    vrrs[6,2,2] = QC[1]*vrrs[6,1,2] + WQ[1]*vrrs[6,1,3] + 1/2ze*vrrs[3,1,3]
    vrrs[7,2,1] = QC[1]*vrrs[7,1,1] + WQ[1]*vrrs[7,1,2] + 1/2ze*vrrs[4,1,2]
    vrrs[7,2,2] = QC[1]*vrrs[7,1,2] + WQ[1]*vrrs[7,1,3] + 1/2ze*vrrs[4,1,3]
    vrrs[8,2,1] = QC[1]*vrrs[8,1,1] + WQ[1]*vrrs[8,1,2]
    vrrs[8,2,2] = QC[1]*vrrs[8,1,2] + WQ[1]*vrrs[8,1,3]
    vrrs[9,2,1] = QC[1]*vrrs[9,1,1] + WQ[1]*vrrs[9,1,2]
    vrrs[9,2,2] = QC[1]*vrrs[9,1,2] + WQ[1]*vrrs[9,1,3]
    vrrs[10,2,1] = QC[1]*vrrs[10,1,1] + WQ[1]*vrrs[10,1,2]
    vrrs[10,2,2] = QC[1]*vrrs[10,1,2] + WQ[1]*vrrs[10,1,3]
    vrrs[1,3,1] = QC[2]*vrrs[1,1,1] + WQ[2]*vrrs[1,1,2]
    vrrs[1,3,2] = QC[2]*vrrs[1,1,2] + WQ[2]*vrrs[1,1,3]
    vrrs[1,3,3] = QC[2]*vrrs[1,1,3] + WQ[2]*vrrs[1,1,4]
    vrrs[1,3,4] = QC[2]*vrrs[1,1,4] + WQ[2]*vrrs[1,1,5]
    vrrs[2,3,1] = QC[2]*vrrs[2,1,1] + WQ[2]*vrrs[2,1,2]
    vrrs[2,3,2] = QC[2]*vrrs[2,1,2] + WQ[2]*vrrs[2,1,3]
    vrrs[2,3,3] = QC[2]*vrrs[2,1,3] + WQ[2]*vrrs[2,1,4]
    vrrs[3,3,1] = QC[2]*vrrs[3,1,1] + WQ[2]*vrrs[3,1,2] + 1/2ze*vrrs[1,1,2]
    vrrs[3,3,2] = QC[2]*vrrs[3,1,2] + WQ[2]*vrrs[3,1,3] + 1/2ze*vrrs[1,1,3]
    vrrs[3,3,3] = QC[2]*vrrs[3,1,3] + WQ[2]*vrrs[3,1,4] + 1/2ze*vrrs[1,1,4]
    vrrs[4,3,1] = QC[2]*vrrs[4,1,1] + WQ[2]*vrrs[4,1,2]
    vrrs[4,3,2] = QC[2]*vrrs[4,1,2] + WQ[2]*vrrs[4,1,3]
    vrrs[4,3,3] = QC[2]*vrrs[4,1,3] + WQ[2]*vrrs[4,1,4]
    vrrs[5,3,1] = QC[2]*vrrs[5,1,1] + WQ[2]*vrrs[5,1,2]
    vrrs[5,3,2] = QC[2]*vrrs[5,1,2] + WQ[2]*vrrs[5,1,3]
    vrrs[6,3,1] = QC[2]*vrrs[6,1,1] + WQ[2]*vrrs[6,1,2] + 1/2ze*vrrs[2,1,2]
    vrrs[6,3,2] = QC[2]*vrrs[6,1,2] + WQ[2]*vrrs[6,1,3] + 1/2ze*vrrs[2,1,3]
    vrrs[7,3,1] = QC[2]*vrrs[7,1,1] + WQ[2]*vrrs[7,1,2]
    vrrs[7,3,2] = QC[2]*vrrs[7,1,2] + WQ[2]*vrrs[7,1,3]
    vrrs[8,3,1] = QC[2]*vrrs[8,1,1] + WQ[2]*vrrs[8,1,2] + 2/2ze*vrrs[3,1,2]
    vrrs[8,3,2] = QC[2]*vrrs[8,1,2] + WQ[2]*vrrs[8,1,3] + 2/2ze*vrrs[3,1,3]
    vrrs[9,3,1] = QC[2]*vrrs[9,1,1] + WQ[2]*vrrs[9,1,2] + 1/2ze*vrrs[4,1,2]
    vrrs[9,3,2] = QC[2]*vrrs[9,1,2] + WQ[2]*vrrs[9,1,3] + 1/2ze*vrrs[4,1,3]
    vrrs[10,3,1] = QC[2]*vrrs[10,1,1] + WQ[2]*vrrs[10,1,2]
    vrrs[10,3,2] = QC[2]*vrrs[10,1,2] + WQ[2]*vrrs[10,1,3]
    vrrs[1,4,1] = QC[3]*vrrs[1,1,1] + WQ[3]*vrrs[1,1,2]
    vrrs[1,4,2] = QC[3]*vrrs[1,1,2] + WQ[3]*vrrs[1,1,3]
    vrrs[1,4,3] = QC[3]*vrrs[1,1,3] + WQ[3]*vrrs[1,1,4]
    vrrs[1,4,4] = QC[3]*vrrs[1,1,4] + WQ[3]*vrrs[1,1,5]
    vrrs[2,4,1] = QC[3]*vrrs[2,1,1] + WQ[3]*vrrs[2,1,2]
    vrrs[2,4,2] = QC[3]*vrrs[2,1,2] + WQ[3]*vrrs[2,1,3]
    vrrs[2,4,3] = QC[3]*vrrs[2,1,3] + WQ[3]*vrrs[2,1,4]
    vrrs[3,4,1] = QC[3]*vrrs[3,1,1] + WQ[3]*vrrs[3,1,2]
    vrrs[3,4,2] = QC[3]*vrrs[3,1,2] + WQ[3]*vrrs[3,1,3]
    vrrs[3,4,3] = QC[3]*vrrs[3,1,3] + WQ[3]*vrrs[3,1,4]
    vrrs[4,4,1] = QC[3]*vrrs[4,1,1] + WQ[3]*vrrs[4,1,2] + 1/2ze*vrrs[1,1,2]
    vrrs[4,4,2] = QC[3]*vrrs[4,1,2] + WQ[3]*vrrs[4,1,3] + 1/2ze*vrrs[1,1,3]
    vrrs[4,4,3] = QC[3]*vrrs[4,1,3] + WQ[3]*vrrs[4,1,4] + 1/2ze*vrrs[1,1,4]
    vrrs[5,4,1] = QC[3]*vrrs[5,1,1] + WQ[3]*vrrs[5,1,2]
    vrrs[5,4,2] = QC[3]*vrrs[5,1,2] + WQ[3]*vrrs[5,1,3]
    vrrs[6,4,1] = QC[3]*vrrs[6,1,1] + WQ[3]*vrrs[6,1,2]
    vrrs[6,4,2] = QC[3]*vrrs[6,1,2] + WQ[3]*vrrs[6,1,3]
    vrrs[7,4,1] = QC[3]*vrrs[7,1,1] + WQ[3]*vrrs[7,1,2] + 1/2ze*vrrs[2,1,2]
    vrrs[7,4,2] = QC[3]*vrrs[7,1,2] + WQ[3]*vrrs[7,1,3] + 1/2ze*vrrs[2,1,3]
    vrrs[8,4,1] = QC[3]*vrrs[8,1,1] + WQ[3]*vrrs[8,1,2]
    vrrs[8,4,2] = QC[3]*vrrs[8,1,2] + WQ[3]*vrrs[8,1,3]
    vrrs[9,4,1] = QC[3]*vrrs[9,1,1] + WQ[3]*vrrs[9,1,2] + 1/2ze*vrrs[3,1,2]
    vrrs[9,4,2] = QC[3]*vrrs[9,1,2] + WQ[3]*vrrs[9,1,3] + 1/2ze*vrrs[3,1,3]
    vrrs[10,4,1] = QC[3]*vrrs[10,1,1] + WQ[3]*vrrs[10,1,2] + 2/2ze*vrrs[4,1,2]
    vrrs[10,4,2] = QC[3]*vrrs[10,1,2] + WQ[3]*vrrs[10,1,3] + 2/2ze*vrrs[4,1,3]
    vrrs[1,5,1] = QC[1]*vrrs[1,2,1] + WQ[1]*vrrs[1,2,2] + 1/2eta*(vrrs[1,1,1]-zeta/ze*vrrs[1,1,2])
    vrrs[1,5,2] = QC[1]*vrrs[1,2,2] + WQ[1]*vrrs[1,2,3] + 1/2eta*(vrrs[1,1,2]-zeta/ze*vrrs[1,1,3])
    vrrs[1,5,3] = QC[1]*vrrs[1,2,3] + WQ[1]*vrrs[1,2,4] + 1/2eta*(vrrs[1,1,3]-zeta/ze*vrrs[1,1,4])
    vrrs[2,5,1] = QC[1]*vrrs[2,2,1] + WQ[1]*vrrs[2,2,2] + 1/2eta*(vrrs[2,1,1]-zeta/ze*vrrs[2,1,2]) + 1/2ze*vrrs[1,2,2]
    vrrs[2,5,2] = QC[1]*vrrs[2,2,2] + WQ[1]*vrrs[2,2,3] + 1/2eta*(vrrs[2,1,2]-zeta/ze*vrrs[2,1,3]) + 1/2ze*vrrs[1,2,3]
    vrrs[3,5,1] = QC[1]*vrrs[3,2,1] + WQ[1]*vrrs[3,2,2] + 1/2eta*(vrrs[3,1,1]-zeta/ze*vrrs[3,1,2])
    vrrs[3,5,2] = QC[1]*vrrs[3,2,2] + WQ[1]*vrrs[3,2,3] + 1/2eta*(vrrs[3,1,2]-zeta/ze*vrrs[3,1,3])
    vrrs[4,5,1] = QC[1]*vrrs[4,2,1] + WQ[1]*vrrs[4,2,2] + 1/2eta*(vrrs[4,1,1]-zeta/ze*vrrs[4,1,2])
    vrrs[4,5,2] = QC[1]*vrrs[4,2,2] + WQ[1]*vrrs[4,2,3] + 1/2eta*(vrrs[4,1,2]-zeta/ze*vrrs[4,1,3])
    vrrs[5,5,1] = QC[1]*vrrs[5,2,1] + WQ[1]*vrrs[5,2,2] + 1/2eta*(vrrs[5,1,1]-zeta/ze*vrrs[5,1,2]) + 2/2ze*vrrs[2,2,2]
    vrrs[6,5,1] = QC[1]*vrrs[6,2,1] + WQ[1]*vrrs[6,2,2] + 1/2eta*(vrrs[6,1,1]-zeta/ze*vrrs[6,1,2]) + 1/2ze*vrrs[3,2,2]
    vrrs[7,5,1] = QC[1]*vrrs[7,2,1] + WQ[1]*vrrs[7,2,2] + 1/2eta*(vrrs[7,1,1]-zeta/ze*vrrs[7,1,2]) + 1/2ze*vrrs[4,2,2]
    vrrs[8,5,1] = QC[1]*vrrs[8,2,1] + WQ[1]*vrrs[8,2,2] + 1/2eta*(vrrs[8,1,1]-zeta/ze*vrrs[8,1,2])
    vrrs[9,5,1] = QC[1]*vrrs[9,2,1] + WQ[1]*vrrs[9,2,2] + 1/2eta*(vrrs[9,1,1]-zeta/ze*vrrs[9,1,2])
    vrrs[10,5,1] = QC[1]*vrrs[10,2,1] + WQ[1]*vrrs[10,2,2] + 1/2eta*(vrrs[10,1,1]-zeta/ze*vrrs[10,1,2])
    vrrs[1,6,1] = QC[1]*vrrs[1,3,1] + WQ[1]*vrrs[1,3,2]
    vrrs[1,6,2] = QC[1]*vrrs[1,3,2] + WQ[1]*vrrs[1,3,3]
    vrrs[1,6,3] = QC[1]*vrrs[1,3,3] + WQ[1]*vrrs[1,3,4]
    vrrs[2,6,1] = QC[1]*vrrs[2,3,1] + WQ[1]*vrrs[2,3,2] + 1/2ze*vrrs[1,3,2]
    vrrs[2,6,2] = QC[1]*vrrs[2,3,2] + WQ[1]*vrrs[2,3,3] + 1/2ze*vrrs[1,3,3]
    vrrs[3,6,1] = QC[1]*vrrs[3,3,1] + WQ[1]*vrrs[3,3,2]
    vrrs[3,6,2] = QC[1]*vrrs[3,3,2] + WQ[1]*vrrs[3,3,3]
    vrrs[4,6,1] = QC[1]*vrrs[4,3,1] + WQ[1]*vrrs[4,3,2]
    vrrs[4,6,2] = QC[1]*vrrs[4,3,2] + WQ[1]*vrrs[4,3,3]
    vrrs[5,6,1] = QC[1]*vrrs[5,3,1] + WQ[1]*vrrs[5,3,2] + 2/2ze*vrrs[2,3,2]
    vrrs[6,6,1] = QC[1]*vrrs[6,3,1] + WQ[1]*vrrs[6,3,2] + 1/2ze*vrrs[3,3,2]
    vrrs[7,6,1] = QC[1]*vrrs[7,3,1] + WQ[1]*vrrs[7,3,2] + 1/2ze*vrrs[4,3,2]
    vrrs[8,6,1] = QC[1]*vrrs[8,3,1] + WQ[1]*vrrs[8,3,2]
    vrrs[9,6,1] = QC[1]*vrrs[9,3,1] + WQ[1]*vrrs[9,3,2]
    vrrs[10,6,1] = QC[1]*vrrs[10,3,1] + WQ[1]*vrrs[10,3,2]
    vrrs[1,7,1] = QC[1]*vrrs[1,4,1] + WQ[1]*vrrs[1,4,2]
    vrrs[1,7,2] = QC[1]*vrrs[1,4,2] + WQ[1]*vrrs[1,4,3]
    vrrs[1,7,3] = QC[1]*vrrs[1,4,3] + WQ[1]*vrrs[1,4,4]
    vrrs[2,7,1] = QC[1]*vrrs[2,4,1] + WQ[1]*vrrs[2,4,2] + 1/2ze*vrrs[1,4,2]
    vrrs[2,7,2] = QC[1]*vrrs[2,4,2] + WQ[1]*vrrs[2,4,3] + 1/2ze*vrrs[1,4,3]
    vrrs[3,7,1] = QC[1]*vrrs[3,4,1] + WQ[1]*vrrs[3,4,2]
    vrrs[3,7,2] = QC[1]*vrrs[3,4,2] + WQ[1]*vrrs[3,4,3]
    vrrs[4,7,1] = QC[1]*vrrs[4,4,1] + WQ[1]*vrrs[4,4,2]
    vrrs[4,7,2] = QC[1]*vrrs[4,4,2] + WQ[1]*vrrs[4,4,3]
    vrrs[5,7,1] = QC[1]*vrrs[5,4,1] + WQ[1]*vrrs[5,4,2] + 2/2ze*vrrs[2,4,2]
    vrrs[6,7,1] = QC[1]*vrrs[6,4,1] + WQ[1]*vrrs[6,4,2] + 1/2ze*vrrs[3,4,2]
    vrrs[7,7,1] = QC[1]*vrrs[7,4,1] + WQ[1]*vrrs[7,4,2] + 1/2ze*vrrs[4,4,2]
    vrrs[8,7,1] = QC[1]*vrrs[8,4,1] + WQ[1]*vrrs[8,4,2]
    vrrs[9,7,1] = QC[1]*vrrs[9,4,1] + WQ[1]*vrrs[9,4,2]
    vrrs[10,7,1] = QC[1]*vrrs[10,4,1] + WQ[1]*vrrs[10,4,2]
    vrrs[1,8,1] = QC[2]*vrrs[1,3,1] + WQ[2]*vrrs[1,3,2] + 1/2eta*(vrrs[1,1,1]-zeta/ze*vrrs[1,1,2])
    vrrs[1,8,2] = QC[2]*vrrs[1,3,2] + WQ[2]*vrrs[1,3,3] + 1/2eta*(vrrs[1,1,2]-zeta/ze*vrrs[1,1,3])
    vrrs[1,8,3] = QC[2]*vrrs[1,3,3] + WQ[2]*vrrs[1,3,4] + 1/2eta*(vrrs[1,1,3]-zeta/ze*vrrs[1,1,4])
    vrrs[2,8,1] = QC[2]*vrrs[2,3,1] + WQ[2]*vrrs[2,3,2] + 1/2eta*(vrrs[2,1,1]-zeta/ze*vrrs[2,1,2])
    vrrs[2,8,2] = QC[2]*vrrs[2,3,2] + WQ[2]*vrrs[2,3,3] + 1/2eta*(vrrs[2,1,2]-zeta/ze*vrrs[2,1,3])
    vrrs[3,8,1] = QC[2]*vrrs[3,3,1] + WQ[2]*vrrs[3,3,2] + 1/2eta*(vrrs[3,1,1]-zeta/ze*vrrs[3,1,2]) + 1/2ze*vrrs[1,3,2]
    vrrs[3,8,2] = QC[2]*vrrs[3,3,2] + WQ[2]*vrrs[3,3,3] + 1/2eta*(vrrs[3,1,2]-zeta/ze*vrrs[3,1,3]) + 1/2ze*vrrs[1,3,3]
    vrrs[4,8,1] = QC[2]*vrrs[4,3,1] + WQ[2]*vrrs[4,3,2] + 1/2eta*(vrrs[4,1,1]-zeta/ze*vrrs[4,1,2])
    vrrs[4,8,2] = QC[2]*vrrs[4,3,2] + WQ[2]*vrrs[4,3,3] + 1/2eta*(vrrs[4,1,2]-zeta/ze*vrrs[4,1,3])
    vrrs[5,8,1] = QC[2]*vrrs[5,3,1] + WQ[2]*vrrs[5,3,2] + 1/2eta*(vrrs[5,1,1]-zeta/ze*vrrs[5,1,2])
    vrrs[6,8,1] = QC[2]*vrrs[6,3,1] + WQ[2]*vrrs[6,3,2] + 1/2eta*(vrrs[6,1,1]-zeta/ze*vrrs[6,1,2]) + 1/2ze*vrrs[2,3,2]
    vrrs[7,8,1] = QC[2]*vrrs[7,3,1] + WQ[2]*vrrs[7,3,2] + 1/2eta*(vrrs[7,1,1]-zeta/ze*vrrs[7,1,2])
    vrrs[8,8,1] = QC[2]*vrrs[8,3,1] + WQ[2]*vrrs[8,3,2] + 1/2eta*(vrrs[8,1,1]-zeta/ze*vrrs[8,1,2]) + 2/2ze*vrrs[3,3,2]
    vrrs[9,8,1] = QC[2]*vrrs[9,3,1] + WQ[2]*vrrs[9,3,2] + 1/2eta*(vrrs[9,1,1]-zeta/ze*vrrs[9,1,2]) + 1/2ze*vrrs[4,3,2]
    vrrs[10,8,1] = QC[2]*vrrs[10,3,1] + WQ[2]*vrrs[10,3,2] + 1/2eta*(vrrs[10,1,1]-zeta/ze*vrrs[10,1,2])
    vrrs[1,9,1] = QC[2]*vrrs[1,4,1] + WQ[2]*vrrs[1,4,2]
    vrrs[1,9,2] = QC[2]*vrrs[1,4,2] + WQ[2]*vrrs[1,4,3]
    vrrs[1,9,3] = QC[2]*vrrs[1,4,3] + WQ[2]*vrrs[1,4,4]
    vrrs[2,9,1] = QC[2]*vrrs[2,4,1] + WQ[2]*vrrs[2,4,2]
    vrrs[2,9,2] = QC[2]*vrrs[2,4,2] + WQ[2]*vrrs[2,4,3]
    vrrs[3,9,1] = QC[2]*vrrs[3,4,1] + WQ[2]*vrrs[3,4,2] + 1/2ze*vrrs[1,4,2]
    vrrs[3,9,2] = QC[2]*vrrs[3,4,2] + WQ[2]*vrrs[3,4,3] + 1/2ze*vrrs[1,4,3]
    vrrs[4,9,1] = QC[2]*vrrs[4,4,1] + WQ[2]*vrrs[4,4,2]
    vrrs[4,9,2] = QC[2]*vrrs[4,4,2] + WQ[2]*vrrs[4,4,3]
    vrrs[5,9,1] = QC[2]*vrrs[5,4,1] + WQ[2]*vrrs[5,4,2]
    vrrs[6,9,1] = QC[2]*vrrs[6,4,1] + WQ[2]*vrrs[6,4,2] + 1/2ze*vrrs[2,4,2]
    vrrs[7,9,1] = QC[2]*vrrs[7,4,1] + WQ[2]*vrrs[7,4,2]
    vrrs[8,9,1] = QC[2]*vrrs[8,4,1] + WQ[2]*vrrs[8,4,2] + 2/2ze*vrrs[3,4,2]
    vrrs[9,9,1] = QC[2]*vrrs[9,4,1] + WQ[2]*vrrs[9,4,2] + 1/2ze*vrrs[4,4,2]
    vrrs[10,9,1] = QC[2]*vrrs[10,4,1] + WQ[2]*vrrs[10,4,2]
    vrrs[1,10,1] = QC[3]*vrrs[1,4,1] + WQ[3]*vrrs[1,4,2] + 1/2eta*(vrrs[1,1,1]-zeta/ze*vrrs[1,1,2])
    vrrs[1,10,2] = QC[3]*vrrs[1,4,2] + WQ[3]*vrrs[1,4,3] + 1/2eta*(vrrs[1,1,2]-zeta/ze*vrrs[1,1,3])
    vrrs[1,10,3] = QC[3]*vrrs[1,4,3] + WQ[3]*vrrs[1,4,4] + 1/2eta*(vrrs[1,1,3]-zeta/ze*vrrs[1,1,4])
    vrrs[2,10,1] = QC[3]*vrrs[2,4,1] + WQ[3]*vrrs[2,4,2] + 1/2eta*(vrrs[2,1,1]-zeta/ze*vrrs[2,1,2])
    vrrs[2,10,2] = QC[3]*vrrs[2,4,2] + WQ[3]*vrrs[2,4,3] + 1/2eta*(vrrs[2,1,2]-zeta/ze*vrrs[2,1,3])
    vrrs[3,10,1] = QC[3]*vrrs[3,4,1] + WQ[3]*vrrs[3,4,2] + 1/2eta*(vrrs[3,1,1]-zeta/ze*vrrs[3,1,2])
    vrrs[3,10,2] = QC[3]*vrrs[3,4,2] + WQ[3]*vrrs[3,4,3] + 1/2eta*(vrrs[3,1,2]-zeta/ze*vrrs[3,1,3])
    vrrs[4,10,1] = QC[3]*vrrs[4,4,1] + WQ[3]*vrrs[4,4,2] + 1/2eta*(vrrs[4,1,1]-zeta/ze*vrrs[4,1,2]) + 1/2ze*vrrs[1,4,2]
    vrrs[4,10,2] = QC[3]*vrrs[4,4,2] + WQ[3]*vrrs[4,4,3] + 1/2eta*(vrrs[4,1,2]-zeta/ze*vrrs[4,1,3]) + 1/2ze*vrrs[1,4,3]
    vrrs[5,10,1] = QC[3]*vrrs[5,4,1] + WQ[3]*vrrs[5,4,2] + 1/2eta*(vrrs[5,1,1]-zeta/ze*vrrs[5,1,2])
    vrrs[6,10,1] = QC[3]*vrrs[6,4,1] + WQ[3]*vrrs[6,4,2] + 1/2eta*(vrrs[6,1,1]-zeta/ze*vrrs[6,1,2])
    vrrs[7,10,1] = QC[3]*vrrs[7,4,1] + WQ[3]*vrrs[7,4,2] + 1/2eta*(vrrs[7,1,1]-zeta/ze*vrrs[7,1,2]) + 1/2ze*vrrs[2,4,2]
    vrrs[8,10,1] = QC[3]*vrrs[8,4,1] + WQ[3]*vrrs[8,4,2] + 1/2eta*(vrrs[8,1,1]-zeta/ze*vrrs[8,1,2])
    vrrs[9,10,1] = QC[3]*vrrs[9,4,1] + WQ[3]*vrrs[9,4,2] + 1/2eta*(vrrs[9,1,1]-zeta/ze*vrrs[9,1,2]) + 1/2ze*vrrs[3,4,2]
    vrrs[10,10,1] = QC[3]*vrrs[10,4,1] + WQ[3]*vrrs[10,4,2] + 1/2eta*(vrrs[10,1,1]-zeta/ze*vrrs[10,1,2]) + 2/2ze*vrrs[4,4,2]

    return vrrs[:,:,1]
end


text:
using DispatchOps
using Test

@testset "DispatchOps.jl" begin
    # Write your tests here.
end


text:
# Record the current counter state.
#
# Useful for recording what the state of the CPU counters is before running any experiments
# so we can restore the state when we are dont.
struct CoreMSRProgramming
    fixed::Int
    programmable::Vector{Int}
end

struct CounterState
    # Maps a CPU to the raw bits of the counter state.
    state::Dict{IndexZero{Int}, CoreMSRProgramming}
end

function CounterState(; cpus = 1:numcpus())
    state = Dict{IndexZero{Int},CoreMSRProgramming}()

    # Collect the counter states for all cpus
    for cpu in cpus
        # Open the MSR forlder for this CPU
        handle = Handle(msrpath(cpu))

        # Read the fixed-counter state
        fixed = read(handle, Int, IA32_FIXED_CTR_CTRL_MSR)

        # Read the programmable counter states
        programmable = Int64[]
        for i in 1:numcounters()
            register = EVENT_SELECT_MSRS[i]
            push!(programmable, read(handle, Int, register))
        end
        close(handle)
        state[indexzero(cpu)] = CoreMSRProgramming(fixed, programmable)
    end

    return CounterState(state)
end

"""
    restore!(cs::CounterState)

Restore the current CPU to the state recorded in `cs`.
"""
function restore!(CS::CounterState)
    state = CS.state
    registers = EVENT_SELECT_MSRS[1:numcounters()]

    for (cpu, programming) in state
        writemsr(cpu, IA32_FIXED_CTR_CTRL_MSR, programming.fixed)
        writemsr.(Ref(cpu), registers, programming.programmable)
    end
    return nothing
end


text:
include("1.jl")
include("5.jl")
function smallsquare(r::Robot)
    angleofsquare(r)
    i = mod(snake(r,West),2)
    if (i == 1)
        i = 1
    else i = 3
    end

    for side in (HorizonSide(mod(i, 4)),Nord,HorizonSide(mod(i+2, 4)),Sud)
        putmarker!(r)
        putmarkhole1(r,side,mod(i+2,4))
        move!(r,HorizonSide(mod(Int(side) + (i+2), 4)))
    end

end
function snake(r::Robot,side1::HorizonSide)
    i = 0
    while !isborder(r,Nord)
        if isborder(r,Ost)
            movehole1(r,side1)
        else
            movehole1(r,inverse(side1))
        end
        if isborder(r,Ost) || isborder(r,West)
            move!(r,Nord)
        end
        i+=1
    end
    return i
end
function movehole1(r::Robot,side::HorizonSide)
    while (!isborder(r,side))
        move!(r,side)
        if (isborder(r,Nord))
            break
        end
    end

end
function putmarkhole1(r::Robot,side::HorizonSide,i::Int64)

    while isborder(r,HorizonSide(mod(Int(side) + i, 4)))
        move!(r,side)
        putmarker!(r)
    end
end

max_stars_repo_name: bubasooo/proggggggggg

max_issues_repo_name: bubasooo/proggggggggg

max_forks_repo_name: bubasooo/proggggggggg


text:
using MPIMeasurements
using Base.Test
using Unitful
using Compat
using HDF5

# define Grid
shape = [5,5,1]
fov = [200.0,200.0,1.0]Unitful.mm
center = [0.0,0.0,0.0]Unitful.mm
positions = MeanderingGridPositions( RegularGridPositions(shape,fov,center) )

# create Scanner
#scanner = MPIScanner("HeadScanner.toml")
scanner = MPIScanner("DummyScanner.toml")
robot = getRobot(scanner)
safety = getSafety(scanner)
gaussmeter = getGaussMeter(scanner)

#TODO mT sollte man hier nicht angeben müssen. Das sollte im Gaussmeter gekapselt sein
mfMeasObj = MagneticFieldMeas(gaussmeter, Unitful.mT, positions)

res = performTour!(robot, safety, positions, mfMeasObj)

filenameField = joinpath(homedir(),"TestBackground.h5")
saveMagneticFieldAsHDF5(mfMeasObj, filenameField)

pos, field = loadMagneticField(filenameField)

MPISimulations.plotMagneticField(field, pos, 3, 1)


text:
function intMnli(Mnli, r, n, l, i; q="best")
    if all(Mnli[n,l,i] .== 0)
        return [0.0],[0.0]
    end
    Mfunc = Spline1D(r,Mnli[n,l,i])
    if q == "best"
        q = best_q(Mfunc, i, 0, rmini, rmaxi)
    end
    #@show q, n, l, i
    k, intM = xicalc(Mfunc, i, 0, N=N, kmin=rmini, kmax=rmaxi, r0=r0i, q=q)
    return k, intM
end

function P22lggs(Pk, l; BiasOperatorDict22, f, Mnli, r)
    P = zeros(Float64,N)
    k,_ = intMnli(Mnli, r, 0, 0, 0)
    for i = 0:8
        for n = 0:4
            try
                mint = k.^n .* intMnli(Mnli, r, n, l, i)[2]
                P .+= (mint .- mint[1])#k.^n .* intMnli(Mnli, r, n, l, i)[2]
                #@show P[1], n, l, i
            catch y
                if isa(y,UndefRefError)

                else
                    @show y, n, l, i
                end
            end
        end
    end
    return k, 2*(2*pi)^3 .* P
end

function P22ggsmultipole(Pk, μs; BiasOperatorDict22, f)
    L(n,x) = Legendre(n,x)
    r, Mnli = generateMnliFast(Pk, BiasOperatorDict22, f)
    Pl = l -> P22lggs(Pk, l, BiasOperatorDict22=BiasOperatorDict22, f=f, Mnli=Mnli, r=r)
    k, P0 = Pl(0)
    _, P2 = Pl(2)
    _, P4 = Pl(4)
    _, P6 = Pl(6)
    _, P8 = Pl(8)
    p22mp = Array{Float64}(undef, length(k), length(μs))
    for i = 1:length(μs)
        @. p22mp[:,i] = L(0,μs[i])*P0 + L(2,μs[i])*P2 + L(4,μs[i])*P4 + L(6,μs[i])*P6 + L(8,μs[i])*P8
    end
    return k, p22mp
end

max_stars_repo_name: JosephTomlinson/GeneralBiasPk

max_issues_repo_name: JosephTomlinson/GeneralBiasPk

max_forks_repo_name: JosephTomlinson/GeneralBiasPk


text:
"""
    WeightedMarginLoss{L,W} <: MarginLoss

Can an be used to represent a re-weighted version of some type of
binary loss `L`. The weight-factor `W`, which must be in `[0, 1]`,
denotes the relative weight of the positive class, while the
relative weight of the negative class will be `1 - W`.
"""
struct WeightedMarginLoss{L<:MarginLoss,W} <: SupervisedLoss
    loss::L
end

_werror() = throw(ArgumentError("The given \"weight\" has to be a number in the interval [0, 1]"))

@generated function WeightedMarginLoss(loss::L, ::Val{W}) where {L<:MarginLoss,W}
    typeof(W) <: Number && 0 <= W <= 1 || _werror()
    :(WeightedMarginLoss{L,W}(loss))
end

function WeightedMarginLoss(loss::MarginLoss, w::Number)
    WeightedMarginLoss(loss, Val(w))
end

@generated function WeightedMarginLoss(s::WeightedMarginLoss{T,W1}, ::Val{W2}) where {T,W1,W2}
    :(WeightedMarginLoss(s.loss, Val($(W1*W2))))
end

function WeightedMarginLoss(s::WeightedMarginLoss{L,W}, w::Number) where {L<:MarginLoss,W}
    WeightedMarginLoss(s.loss, Val(W*w))
end

for fun in (:value, :deriv, :deriv2)
    @eval function ($fun)(l::WeightedMarginLoss{L,W}, target::Number, output::Number) where {L,W}
        # We interpret the W to be the weight of the positive class
        if target > 0
            W * ($fun)(l.loss, target, output)
        else
            (1-W) * ($fun)(l.loss, target, output)
        end
    end
end

# An α-weighted version of a classification callibrated margin loss is
# itself classification callibrated if and only if α == 1/2
isclasscalibrated(l::WeightedMarginLoss{T,W}) where {T,W} = W == 0.5 && isclasscalibrated(l.loss)

# TODO: Think about this semantic
issymmetric(::WeightedMarginLoss) = false

for prop in [:isminimizable, :isdifferentiable,
             :istwicedifferentiable,
             :isconvex, :isstrictlyconvex,
             :isstronglyconvex, :isnemitski,
             :isunivfishercons, :isfishercons,
             :islipschitzcont, :islocallylipschitzcont,
             :isclipable, :ismarginbased,
             :isdistancebased]
    @eval ($prop)(l::WeightedMarginLoss) = ($prop)(l.loss)
end

for prop_param in (:isdifferentiable, :istwicedifferentiable)
    @eval ($prop_param)(l::WeightedMarginLoss, at) = ($prop_param)(l.loss, at)
end


text:
using Test

include("../src/FiniteGroups.jl")

@testset "coset" begin
    @test coset(1, Set([1, 2, 3]), *) == Set([1, 2, 3])
    @test coset(0, Set([-1, 3, 4, 5]), *) == Set([0])
    # @test coset(Set([1]), Set([Set([1 2]), Set([1, 3])]), ((x, y)->coset(x, y, +))) == Set([Set([2, 4]), Set([2, 3])])

end

@testset "set composition" begin
    @test set_composition(Set([1]), Set([1, 2]), (+)) == Set([2, 3])
end

@testset  "validate group" begin
    G = Group(Set([1, 2]), (x, y)-> x * y % 3)
    @test_throws AssertionError Group(Set([-1, 1]), (+))
end

@testset  "coset{group}" begin
    G = Group(Set([1, 2]), (x, y)-> x * y % 3)
    @test coset(0, G) == Set([0])
    @test coset(2, G) == Set([1, 2])
end

@testset "naive_isequal" begin
    add_mod_2(x, y) = (x + y) % 2
    mod_3(x, y) = mod(x*y, 3)

    G = Group(Set([0, 1]), add_mod_2)
    @test naive_isequal(G, Group(Set([0]), add_mod_2)) == false
    @test naive_isequal(G, G)
    @test_throws MethodError naive_isequal(G, SubGroup(G, Set([1]), mod_3))
end

@testset  "subgroup definition" begin
    mod_3(x, y) = mod(x*y, 3)
    G = Group(Set([1, 2]), mod_3)
    @test SubGroup(G, Set([1])) == SubGroup(G, Set([1]), mod_3)
    @test G == SubGroup(G, G.set)
    @test SubGroup(G, G.set) == G
    @test SubGroup(G, Set([1])) != G
end

@testset "normal subgroup" begin
    # In a cyclic group, all subgroups are normal
    ×₇(x, y) = mod(x * y, 7)
    G = Group(Set([1 2 3 4 5 6]), ×₇)
    N = SubGroup(G, Set([1 6]))
    @test isnormal(N)
    @test isnormal(SubGroup(G, G.set))
    @test isnormal(SubGroup(G, Set([1])))
    @test inv(G, 2) == 4
    @test conjugate(G, 2, 4) == 2
    # Symmetry group of triangle
    mapping = Dict( # Non-commutative for `:r`, `:s`, `:t`
        :e => Dict(:e => :e, :a => :a, :b => :b, :t => :t, :r => :r, :s => :s),
        :a => Dict(:e => :a, :a => :b, :b => :e, :t => :s, :r => :t, :s => :r),
        :b => Dict(:e => :b, :a => :e, :b => :a, :t => :r, :r => :s, :s => :t),
        :r => Dict(:e => :r, :a => :s, :b => :t, :t => :b, :r => :e, :s => :a),
        :s => Dict(:e => :s, :a => :t, :b => :r, :t => :a, :r => :b, :s => :e),
        :t => Dict(:e => :t, :a => :r, :b => :s,:t => :e, :r => :a, :s => :b),

    )

    f(x, y) = mapping[x][y]
    G = Group(Set([:e :a :b :r :s :t]), f)
    S = SubGroup(G, Set([:e :r]))
    @test isnormal(S) == false
end

@testset "conjugacy classes" begin

    mapping = Dict( # Non-commutative for `:r`, `:s`, `:t`
    :e => Dict(:e => :e, :a => :a, :b => :b, :t => :t, :r => :r, :s => :s),
    :a => Dict(:e => :a, :a => :b, :b => :e, :t => :s, :r => :t, :s => :r),
    :b => Dict(:e => :b, :a => :e, :b => :a, :t => :r, :r => :s, :s => :t),
    :r => Dict(:e => :r, :a => :s, :b => :t, :t => :b, :r => :e, :s => :a),
    :s => Dict(:e => :s, :a => :t, :b => :r, :t => :a, :r => :b, :s => :e),
    :t => Dict(:e => :t, :a => :r, :b => :s,:t => :e, :r => :a, :s => :b),

    )
    f(x, y) = mapping[x][y]
    G = Group(Set([:e :a :b :r :s :t]), f)

    @test conjugacy_class(G, :t) == conjugacy_class(G, :s) == conjugacy_class(G, :r) == Set([:s :r :t])
    @test conjugacy_class(G, :a) == Set([:a :b])
    @test conjugacy_classes(G) == Set([Set([:a :b]), Set([:s :r :t]), Set([:e])])

end

@testset "order" begin
    ×₇(x, y) = mod(x * y, 7)
    G = Group(Set([1 2 3 4 5 6]), ×₇)
    @test order(G) == 6
    @test order(G, 1) == 1
    @test order(G, 2) == 3
end

@testset "generators" begin
    ×₇(x, y) = mod(x * y, 7)
    G = Group(Set([1 2 3 4 5 6]), ×₇)
    N = SubGroup(G, Set([1 6]))
    @test generate_subgroup(G, 6) == N
    @test generate_subgroup(G, 1) == SubGroup(G, Set([1]))
    @test generate_subgroup(G, 3) == G
    @test find_generators(G) == Set([3 5])
    @test iscyclic(G)

    mapping = Dict( # Non-commutative for `:r`, `:s`, `:t`
    :e => Dict(:e => :e, :a => :a, :b => :b, :t => :t, :r => :r, :s => :s),
    :a => Dict(:e => :a, :a => :b, :b => :e, :t => :s, :r => :t, :s => :r),
    :b => Dict(:e => :b, :a => :e, :b => :a, :t => :r, :r => :s, :s => :t),
    :r => Dict(:e => :r, :a => :s, :b => :t, :t => :b, :r => :e, :s => :a),
    :s => Dict(:e => :s, :a => :t, :b => :r, :t => :a, :r => :b, :s => :e),
    :t => Dict(:e => :t, :a => :r, :b => :s,:t => :e, :r => :a, :s => :b),

)
    f(x, y) = mapping[x][y]
    G = Group(Set([:e :a :b :r :s :t]), f)
    @test find_generators(G) == Set([])
    @test iscyclic(G) == false
end

@testset "ismorphism" begin
    U12 = Group(Set([1 5 7 11]), (x, y)->mod(x*y, 12))
    mapping = Dict(
        :e => Dict(:e=>:e, :a=>:a, :r=>:r, :s=>:s),
        :a => Dict(:e=>:a, :a=>:e, :r=>:s, :s=>:r),
        :r => Dict(:e=>:r, :a=>:s, :r=>:e, :s=>:a),
        :s => Dict(:e =>:s, :a=>:r, :r=>:a, :s=>:e)
    )
    f(x, y) = mapping[x][y]
    S_rectangle = Group(keys(mapping) |> Set, f)

    trivial_homomorphism = Dict(k => :e for k in U12)
    _assert_homomorphism_property(trivial_homomorphism, U12, S_rectangle)

    # Automorphism
    @test get_isomorphism(S_rectangle, S_rectangle) == Dict(:e=>:e, :a=>:a, :r=>:r, :s=>:s)

    # Kernel, image
    @test image(trivial_homomorphism, U12) == Set([:e])
end

@testset "group actions" begin
    U12 = Group(Set([1 5 7 11]), (x, y)->mod(x*y, 12))
    GA = GroupAction(U12, Set(1:3), (x,y) -> mod(x*y, 4))  # 4 is a divisor of 12, so will be preserving composites
    @test orbit(GA, 1) == Set([1 3])
    @test orbits(GA) == Set([Set([1, 3]), Set([2])])
    @test stabilizer(GA, 1) == Set([1 5])
    @test fixed_set(GA, 1) == Set(1:3)
    @test fixed_set(GA, 11) == Set([2])
end


text:
"""
# LocalPackageServer

Implementation of a Julia Storage Server for a single registry, where
the registry and the registered packages are dynamically retrieved
from their git repositories. It can also be configured as a
low-feature Julia Package Server with a built in storage server.
See the package's `README.md` for more information.
"""
module LocalPackageServer

using Pkg
using Dates
using HTTP
using Random
using Tar

include("config.jl")
include("resource.jl")
include("meta.jl")
include("gitserver.jl")

function __init__()
    # Set default HTTP useragent
    HTTP.setuseragent!("LocalPackageServer (HTTP.jl)")
end

start(config) = start(Config(config))

function start(config::Config)
    host = config.host
    port = config.port
    mkpath(config.cache_dir)
    @info "server listening on $(host):$(port)" Dates.now()
    HTTP.listen(host, port) do http
        resource = http.message.target
        @info "resource requested" resource Dates.now()
        # If the user is asking for `/meta`, generate the
        # requisite JSON object and send it back.
        if occursin(meta_re, resource)
            serve_meta(config, http)
            return
        end

        # If the user asked for something that is an actual
        # resource, send it directly.
        if occursin(resource_re, resource)
            path = fetch(config, resource)
            if path !== nothing
                serve_file(http, path)
                return
            end
        end
        HTTP.setstatus(http, 404)
        startwrite(http)
    end
end

end # module


text:
## methods to output the GASS results, mainly in CASA format
##

### save an array in CASA format using cfg mapping
###
function save_CASAarr(fileprefix::String, cfg::cfg, bestsub)
    headerCasa= "# observatory=ALMA\n# coordsys=LOC (local tangent plane)\n# x y z diam pad#\n"

    for i in 1:length(bestsub)
        arrStr= string(headerCasa)
        subname= cfg.sub.Subarray_Name[i]

        for j in 1:length(bestsub[i])
            arrStr *= @sprintf("%f %f %f %2.1f %s\n",cfg.arr.X[bestsub[i][j]], cfg.arr.Y[bestsub[i][j]] ,
                cfg.arr.Z[bestsub[i][j]] , cfg.arr.diam[bestsub[i][j]] , cfg.arr.name[bestsub[i][j]])
        end

        filename= fileprefix*"_"*subname*".cfg"
        f= open(filename, "w")
        print(f,arrStr)
        close(f)
        println("## $filename written...")
    end
end

## sort population wrt. score
##
function sort_population(p::population)
    s= sortperm(p.score, rev=true)
    psorted= population(p.age, p.subarr[s,:], p.fitness[s,:], p.score[s], p.param[s,:])
    return(psorted)
end

## full run with JLD output ..
##
function full_run(fileprefix::String, cfg::cfg, jld= false)
    res= gass_optimization(cfg)

    lastPopulation= cfg.ga.Number_Iterations
    psort= sort_population(res[lastPopulation])
    save_CASAarr(fileprefix, cfg, psort.subarr[1,:])

    if jld
       JLD.save(fileprefix*".jld", "gass", res)
       println("## $fileprefix.jld written..")
    end

    return(res)
end

max_stars_repo_name: bosscha/GASS

max_issues_repo_name: bosscha/GASS

max_forks_repo_name: bosscha/GASS


text:
module MakieRipserer

using AbstractPlotting
using IterTools
using PersistenceDiagrams
using PlotUtils
using ProgressMeter
using Ripserer

using AbstractPlotting: Plot, PointBased, Triangle
using Ripserer:
    AbstractSimplex, AbstractFiltration, AbstractRipsFiltration, AbstractChainElement
using GeometryBasics: GLTriangleFace

import GeometryBasics

export plot_barcode, FilteredChain, ObservableChain

const DEFAULT_PALETTE = :default

# This allows us to provide either color names or integers for colors.
function get_color(p, name)
    lift(name, p[:palette]) do color, scheme
        if color isa Integer
            colors = PlotUtils.get_colorscheme(scheme)
            colors[mod1(color, length(colors))]
        else
            color
        end
    end
end

include("conversion.jl")
include("filteredchain.jl")
include("observablechain.jl")
include("diagrams.jl")
include("observablediagrams.jl")
include("app.jl")
include("ripserer.jl")

end


text:
####################################################
#                                                  #
# Compute the no policy / no intervention outcome  #
#                                                  #
####################################################

"""
    nopolicy(; TT=52)
Simulate the model in Chari, Kirpalani, and Phelan (2021) for `TT` periods under no government response.
The output is a [`Dict`](https://docs.julialang.org/en/v1/base/collections/#Dictionaries) containing
the necessary variables and time-series for generating the figures.
"""
function nopolicy(; TT=52)

    # Parameters
    ISO = 0;
    TEST = 1;
    OPT = 0; #0 for fmincon, 1 for global, 2 for interp
    beta = .99; # Discount factor
    gamma = 7/18; # Exit rate
    delta = .005; # Death rate
    sigma = 2; # CES parameter on production function
    B1 = 1; # Productivity of work sector
    B2 = .1; # Productivity of home sector
    nu = 1; # Testing cost parameter
    c = (1/.01)*50*52/63000 .* YY_noep(B1, B2, sigma, 1); # Testing cost parameter
    z = 7500*52/63000 .* YY_noep(B1, B2, sigma, 1); # Treatment cost
    d = 15*3*52*(1/.96^(1/52)-1); # Death cost
    xi = .8; # Infected productivity los
    tautilde = 13/18*1/2; # Prob. of becoming symptomatic
    pi2 = .01; # Infection rate at home
    pi1 = 1.5; # Infection rate at work
    thets = .0044; #Signal probability from susceptible person
    theti = .38; #Signal probability from infected person
    ut = 1;

    # Combine into (immutable) NamedTuple
    Par = (ISO=ISO, TEST=TEST, OPT=OPT, beta=beta, gamma=gamma, delta=delta, sigma=sigma, B1=B1, B2=B2,
    nu=nu, c=c, z=z, d=d, xi=xi, tautilde=tautilde, pi2=pi2, pi1=pi1, thets=thets, theti=theti, ut=ut);

    # Empty vectors for simulation results
    I0 = 0.02;
    S0 = 1-I0/(1-Par.tautilde);
    R0  = 0;
    IT = zeros(TT);
    ST = zeros(TT);
    RT = zeros(TT);
    qT = zeros(TT);
    lamT = zeros(TT);
    muT = zeros(TT);
    #Del = zeros(TT);
    VT = zeros(TT);
    YT = zeros(TT);
    L1 = zeros(TT);
    L2 = zeros(TT);
    Rknot = zeros(TT);
    Test = zeros(TT);
    IN = zeros(TT);
    C1T = zeros(TT);
    C2T = zeros(TT);
    IT[1] = I0;
    ST[1] = S0;
    RT[1] = R0;
    VV = 0.0;
    VVC = zeros(TT);
    masstest = zeros(TT);
    mctest = zeros(TT);

    # Simulate optimal path
    for i in 1:TT
        S = ST[i];
        I = IT[i];
        R = RT[i];
        IN[i] = 1-S-R;
        qT[i] = 0;
        lamT[i] = .9901;
        muT[i] = .9901;
        Del = 0;
        LI1 = lamT[i].*(S+(1-Del).*I)+muT[i].*(1-Par.delta).*R;
        LI2 = (1-lamT[i]).*(S+(1-Del).*I)+(1-muT[i]).*(1-Par.delta).*R;
        if LI1>0
            Ipr1 = (lamT[i].^2 .* Par.pi1.*(1-Del).*I)./LI1;
            RIpr1 = (lamT[i].^2 .* Par.pi1)./LI1;
        else
            Ipr1 = 0;
            RIpr1 = 0;
        end
        if LI2>0
            Ipr2 = ((1-lamT[i]).^2 .* Par.pi2.*(1-Del).*I)./LI2;
            RIpr2 = ((1-lamT[i]).^2 .* Par.pi2)./LI2;
        else
            Ipr2 = 0;
            RIpr2 = 0;
        end
        Ipr = Ipr1 + Ipr2;
        Ipr = min(Ipr,1);
        if i<TT
            IT[i+1] = min(max((1-Del).*(1-Par.tautilde).*(1-Par.gamma).*I+S.*Ipr,0),1);
            ST[i+1] = min(max((1-Ipr).*S,0),1);
            RT[i+1] = max(min(R + Par.gamma.*IN[i],1),0);
        end

        YT[i] = YY(Par,S,I,R,qT[i],lamT[i],muT[i]);
        C1T[i],C2T[i] = CC(Par,S,I,R,qT[i],lamT[i],muT[i]);
        L1[i] = (lamT[i].*(S+Par.xi.*(1-Del).*I)+muT[i].*(1-Par.delta).*R);
        L2[i] = ((1-lamT[i]).*(S+Par.xi.*(1-Del).*I)+(1-muT[i]).*(1-Par.delta).*R);

        if I>0
            Rknot[i]=(1 ./ (Par.gamma)) .* Ipr ./ (1-S-R);
        else
            Rknot[i] = 0;
        end
        VV += Par.beta.^(i-1) * (log(YT[i]-C1T[i])-C2T[i]);
        VVC[i] = VV;
    end

    # Calculate output variables
    VV = VV + Par.beta^TT/(1-Par.beta)*VTerm(Par,ST[TT],IT[TT],RT[TT]);
    cumdeaths = cumsum(Par.gamma.*Par.delta.*max.(1 .- ST .- RT,0));
    Totaldeath =cumdeaths[TT];
    Ycumsum = YT[1]-YY_noep(Par,1);
    for i = 2:TT
        Ycumsum = Ycumsum + Par.beta.^(i-1) .* (YT[i]-YY_noep(Par,1));
    end
    Ycumsum = (1-Par.beta).*Ycumsum;
    compvar2 = 1/YY_noep(Par,1)*exp((1-Par.beta)*VV);
    compvar2 = 1-compvar2;
    #suscepss = ST[TT];
    #Iss = 1-ST[TT]-RT[TT];
    #maxI = maximum(1 .- ST .- RT);
    #R01 = Rknot[1];
    #R026 = Rknot[26];
    #R052 = Rknot[52];

    return Dict("Par" => Par, "TT" => TT, "ST" => ST, "IT" => IT, "RT" => RT,
    "YT" => YT, "C1T" => C1T, "Rknot" => Rknot, "masstest" => masstest, "mctest" => mctest)

end


text:
using AutomotiveDrivingModels
using NGSIM
using HDF5, JLD

ΔT = NGSIM_TIMESTEP # [s]

immutable Interval
    id::Int
    frame_lo::Int
    frame_hi::Int
end
type IntervalInProgress
    id::Int
    frame_lo::Int
    updated::Bool
end

# include("tim2dfeatures.jl")
# include("radar_features.jl")
include("multifeatureset.jl")

extract_core = true
extract_temporal = true
extract_well_behaved = false
extract_neighbor_features = false
extract_carlidar_rangerate = true
carlidar_nbeams = 20
roadlidar_nbeams = 0
roadlidar_nlanes = 0
carlidar_max_range = 100.0 # [m]
roadlidar_max_range = 50.0 # [m]

extractor = MultiFeatureExtractor(
    extract_core,
    extract_temporal,
    extract_well_behaved,
    extract_neighbor_features,
    extract_carlidar_rangerate,
    carlidar_nbeams,
    roadlidar_nbeams,
    roadlidar_nlanes,
    carlidar_max_range=carlidar_max_range,
    roadlidar_max_range=roadlidar_max_range,
    )

N_FEATURES = length(extractor)

# const ROADWAY = open(io->read(io, Roadway), "roadway_passive_aggressive.txt", "r")
# for trajdata_path in ["trajdata_passive_aggressive1.txt"]
# println("FILE: ", trajdata_path); tic()
#     trajdata = open(io->read(io, Trajdata, ROADWAY), trajdata_path, "r")

culling_fidelity = 5.0 # [m]
lane_portion_max_range = roadlidar_max_range + culling_fidelity/2 + 5.0 # [m] note extra margin of error

err = 0.0

for trajdata_path in TRAJDATA_PATHS
    println("FILE: ", trajdata_path); tic()
    trajdata = load_trajdata(trajdata_path)

    extractor.road_lidar_culling = RoadwayLidarCulling(trajdata.roadway, lane_portion_max_range, culling_fidelity)
    @assert extractor.road_lidar_culling.is_leaf == false

    segments = pull_continuous_segments(trajdata, 1)

    # extract the dataset for each segment
    scene = Scene()
    rec = SceneRecord(2, ΔT)
    n_states = sum(i->length(i.frame_lo:i.frame_hi)-1, segments)
    features = Array(Float64, N_FEATURES)
    featureset = Array(Float64, length(features), n_states)
    targetset = Array(Float64, 2, n_states)
    interval_starts = Array(Int, length(segments))

    i = 0
    for (segind, seg) in enumerate(segments)
        empty!(rec)
        interval_starts[segind] = i+1
        tic()
        for frame in seg.frame_lo : seg.frame_hi
            get!(scene, trajdata, frame)
            update!(rec, scene)
            vehicle_index = get_index_of_first_vehicle_with_id(rec, seg.egoid)
            @assert(vehicle_index != 0)
            pull_features!(extractor, features, rec, trajdata.roadway, vehicle_index)

            if frame > seg.frame_lo
                targetset[1, i] = get(ACC, rec, trajdata.roadway, vehicle_index)
                targetset[2, i] = get(TURNRATEG, rec, trajdata.roadway, vehicle_index)
            end

            if frame < seg.frame_hi
                i += 1
                featureset[:,i] = features
            end
        end
        toc()
    end

    # Save it all to a file
    println("saving...")
    outfile = ("./../2d_drive_data/data_" * splitdir(splitext(trajdata_path)[1])[2]
    * "_clb" * string(carlidar_nbeams)
    * "_rlb" * string(roadlidar_nbeams)
    * "_rll" * string(roadlidar_nlanes)
    * "_clmr" * string(Int(carlidar_max_range))
    * "_rlmr" * string(Int(roadlidar_max_range))
    * ".jld")

    #outfile = "./data_" * splitdir(splitext(trajdata_path)[1])[2] * ".jld"
    JLD.save(outfile, "features", featureset, "targets", targetset, "intervals", interval_starts, "timestep", ΔT)
    toc()
end

max_stars_repo_name: Sawato/gail-driver

max_issues_repo_name: fcdtc/gail-driver

max_forks_repo_name: fcdtc/gail-driver


text:
include(joinpath("structures", "gauging.jl"))
include(joinpath("structures", "ratingcurve.jl"))
include(joinpath("structures", "compoundratingcurve.jl"))