text:
################################################################################
### Helpers
################################################################################

function vmuladd!(C::Matrix{T}, A::Matrix{T}, B::Matrix{T}, factor::T = T(1)) where {T <: Real}
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[m,k] * B[k,n]
        end
        C[m,n] += factor * Cmn
    end
end
function vmuladd!(C::Matrix{T}, A::Matrix{T}, B::Diagonal{T}, factor::T = T(1)) where {T <: Real}
    @turbo for m in 1:size(A, 1), n in 1:size(A, 2)
        C[m,n] += factor * A[m,n] * B.diag[n]
    end
end
function vmuladd!(C::Matrix{T}, A::Matrix{T}, X::Adjoint{T}, factor::T = T(1)) where {T <: Real}
    B = X.parent
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[m,k] * B[n, k]
        end
        C[m,n] += factor * Cmn
    end
end
function vmuladd!(C::Matrix{T}, X::Adjoint{T}, B::Matrix{T}, factor::T = T(1)) where {T <: Real}
    A = X.parent
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[k,m] * B[k,n]
        end
        C[m,n] += factor * Cmn
    end
end
function vmul!(C::Matrix{T}, A::Matrix{T}, X::Adjoint{T}, factor::T) where {T <: Real}
    B = X.parent
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[m,k] * B[n, k]
        end
        C[m,n] = factor * Cmn
    end
end
function vmul!(C::Matrix{T}, X::Adjoint{T}, B::Matrix{T}, factor::T) where {T <: Real}
    A = X.parent
    @turbo for m in 1:size(A, 1), n in 1:size(B, 2)
        Cmn = zero(eltype(C))
        for k in 1:size(A, 2)
            Cmn += A[k,m] * B[k,n]
        end
        C[m,n] = factor * Cmn
    end
end

################################################################################
### Overloads for complex StructArrays
################################################################################

const CMat64 = StructArray{Complex{Float64},2,NamedTuple{(:re, :im), Tuple{AT, AT}}, I} where {AT <: AbstractArray{Float64, 2}, I}
const CVec64 = StructArray{Complex{Float64},1,NamedTuple{(:re, :im), Tuple{AT, AT}}, I} where {AT <: AbstractArray{Float64, 1}, I}

function vmul!(C::CMat64, A::CMat64, B::CMat64)
    @warn "Complex StructArrays are untested not really optimized." maxlog=10
    vmul!(   C.re, A.re, B.re)     # C.re = A.re * B.re
    vmuladd!(C.re, A.im, B.im, -1.0) # C.re = C.re - A.im * B.im
    vmul!(   C.im, A.re, B.im)     # C.im = A.re * B.im
    vmuladd!(C.im, A.im, B.re)     # C.im = C.im + A.im * B.re
end

function vmul!(C::CMat64, A::CMat64, B::Diagonal{T}) where {T <: Real}
    @warn "Complex StructArrays are untested not really optimized." maxlog=10
    vmul!(C.re, A.re, B)
    vmul!(C.im, A.im, B)
end
function vmul!(C::CMat64, A::CMat64, B::Diagonal{ComplexF64, CVec64})
    @warn "Complex StructArrays are untested not really optimized." maxlog=10
    vmul!(   C.re, A.re, B.re)
    vmuladd!(C.re, A.im, B.im, -1.0)
    vmul!(   C.im, A.re, B.im)
    vmuladd!(C.im, A.im, B.re)
end

function vmul!(C::CMat64, A::CMat64, X::Adjoint{T, CMat64}) where {T <: ComplexF64}
    @warn "Complex StructArrays are untested not really optimized." maxlog=10
    B = X.parent
    vmul!(   C.re, A.re, adjoint(B.re))
    vmuladd!(C.re, A.im, adjoint(B.im), 1.0)
    vmul!(   C.im, A.re, adjoint(B.im), -1.0)
    vmuladd!(C.im, A.im, adjoint(B.re))
end
function vmul!(C::CMat64, X::Adjoint{T}, B::CMat64) where {T <: Real}
    @warn "Complex StructArrays are untested not really optimized." maxlog=10
    A = X.parent
    vmul!(   C.re, adjoint(A.re), B.re)
    vmuladd!(C.re, adjoint(A.im), B.im, 1.0)
    vmul!(   C.im, adjoint(A.re), B.im)
    vmuladd!(C.im, adjoint(A.im), B.re, -1.0)
end

function rvmul!(A::CMat64, B::Diagonal{T}) where {T <: Real}
    @warn "Complex StructArrays are untested not really optimized." maxlog=10
    rvmul!(A.re, B)
    rvmul!(A.im, B)
end
function rvmul!(A::CMat64, B::Diagonal{ComplexF64})
    @warn "Complex StructArrays are untested not really optimized." maxlog=10
    @inbounds for m in 1:size(A, 1), n in 1:size(A, 2)
        A[m,n] = A[m, n] * B.diag[n]
    end
end

function lvmul!(A::Diagonal{T}, B::CMat64) where {T <: Real}
    @warn "Complex StructArrays are untested not really optimized." maxlog=10
    lvmul!(A, B.re)
    lvmul!(A, B.im)
end
function lvmul!(A::Diagonal{ComplexF64}, B::CMat64)
    @warn "Complex StructArrays are untested not really optimized." maxlog=10
    @inbounds for m in 1:size(B, 1), n in 1:size(B, 2)
        B[m,n] = A.diag[m] * B[m, n]
    end
end

rvadd!(A::CMat64, D::Diagonal{T}) where {T <: Real} = rvadd!(A.re, D)

##################################################
### TODO
##################################################

function rdivp!(A::CMat64, T::CMat64, O::CMat64, pivot)
    # assume Diagonal is ±1!
    @inbounds begin
        N = size(A, 1)

        # Apply pivot
        for (j, p) in enumerate(pivot)
            @turbo for i in 1:N
                O.re[i, j] = A.re[i, p]
            end
            @turbo for i in 1:N
                O.im[i, j] = A.im[i, p]
            end
        end

        # do the rdiv
        # @turbo will segfault on `k in 1:0`, so pull out first loop
        invT = conj(T[1, 1]) / abs2(T[1, 1])
        re = real(invT); im = imag(invT)
        @turbo for i in 1:N
            A.re[i, 1] = O.re[i, 1] * re
        end
        @turbo for i in 1:N
            A.re[i, 1] -= O.im[i, 1] * im
        end
        @turbo for i in 1:N
            A.im[i, 1] = O.im[i, 1] * re
        end
        @turbo for i in 1:N
            A.im[i, 1] += O.re[i, 1] * im
        end

        # TODO Is this optimal?
        for j in 2:N
            invT = conj(T[j, j]) / abs2(T[j, j])
            re = real(invT); im = imag(invT)

            @turbo for i in 1:N
                x = O.re[i, j]
                for k in 1:j-1
                    x -= A.re[i, k] * T.re[k, j]
                end
                for k in 1:j-1
                    x += A.im[i, k] * T.im[k, j]
                end
                A.re[i, j] = x * re
                A.im[i, j] = x * im
            end

            @turbo for i in 1:N
                x = O.im[i, j]
                for k in 1:j-1
                    x -= A.im[i, k] * T.re[k, j]
                end
                for k in 1:j-1
                    x -= A.re[i, k] * T.im[k, j]
                end
                A.re[i, j] -= x * im
                A.im[i, j] += x * re
            end
        end
    end
    A
end

################################################################################
### UDT
################################################################################

@inline function reflector!(x::CMat64, normu, j=1, n=size(x, 1))
    @inbounds begin
        ξ1 = x[j, j]
        if iszero(normu)
            return zero(ξ1) #zero(ξ1/normu)
        end
        normu = sqrt(normu)
        ν = LinearAlgebra.copysign(normu, real(ξ1))
        ξ1 += ν
        invξ1 = 1.0 / ξ1
        x.re[j, j] = -ν
        @turbo for i = j+1:n
            x.re[i, j] = x.re[i, j] * real(invξ1)
        end
        @turbo for i = j+1:n
            x.re[i, j] = -x.im[i, j] * imag(invξ1)
        end
        @turbo for i = j+1:n
            x.im[i, j] = x.im[i, j] * real(invξ1)
        end
        @turbo for i = j+1:n
            x.im[i, j] = x.re[i, j] * imag(invξ1)
        end
    end
    ξ1/ν
end

@inline function reflectorApply!(x::CVec64, τ::Number, A::CMat64)
    m, n = size(A)
    @inbounds for j = 1:n
        # dot
        vAj_re = A.re[1, j]
        @turbo for i = 2:m
            vAj_re += x.re[i] * A.re[i, j]
        end
        @turbo for i = 2:m
            vAj_re += x.im[i] * A.im[i, j]
        end

        vAj_im = A.im[1, j]
        @turbo for i = 2:m
            vAj_im -= x.im[i] * A.re[i, j]
        end
        @turbo for i = 2:m
            vAj_im += x.re[i] * A.im[i, j]
        end

        temp = real(τ) * vAj_re + imag(τ) * vAj_im
        vAj_im = real(τ) * vAj_im - imag(τ) * vAj_re
        vAj_re = temp

        # ger
        A.re[1, j] -= vAj_re
        @turbo for i = 2:m
            A.re[i, j] -= x.re[i] * vAj_re
        end
        @turbo for i = 2:m
            A.re[i, j] += x.im[i] * vAj_im
        end

        A.im[1, j] -= vAj_im
        @turbo for i = 2:m
            A.im[i, j] -= x.re[i] * vAj_im
        end
        @turbo for i = 2:m
            A.im[i, j] -= x.im[i] * vAj_re
        end
    end
    return A
end

function udt_AVX_pivot!(
        U::CMat64,
        D::AbstractArray{Float64, 1},
        input::CMat64,
        pivot::AbstractArray{Int64, 1} = Vector(UnitRange(1:size(input, 1))),
        temp::AbstractArray{C, 1} = Vector{C}(undef, length(D)),
        apply_pivot::Val = Val(true)
    ) where {C <: Complex}
    # Assumptions:
    # - all matrices same size
    # - input can be changed (input becomes T)
    @warn "Complex StructArrays are untested not really optimized." maxlog=10

    n = size(input, 1)
    @inbounds for i in 1:n
        pivot[i] = i
    end

    # TODO optimize
    @inbounds for j = 1:n
        # Find column with maximum norm in trailing submatrix
        jm, maxval = indmaxcolumn(input, j, n)

        if jm != j
            # Flip elements in pivoting vector
            tmpp = pivot[jm]
            pivot[jm] = pivot[j]
            pivot[j] = tmpp

            # Update matrix with
            @turbo for i = 1:n
                tmp = input.re[i,jm]
                input.re[i,jm] = input.re[i,j]
                input.re[i,j] = tmp
            end
            @turbo for i = 1:n
                tmp = input.im[i,jm]
                input.im[i,jm] = input.im[i,j]
                input.im[i,j] = tmp
            end
        end

        # Compute reflector of columns j
        τj = reflector!(input, maxval, j, n)
        temp[j] = τj

        # Update trailing submatrix with reflector
        x = LinearAlgebra.view(input, j:n, j)
        MonteCarlo.reflectorApply!(x, τj, LinearAlgebra.view(input, j:n, j+1:n))
    end

    copyto!(U.re, I)
    copyto!(U.im, 0.0)
    @inbounds begin
        U.re[n, n] -= temp.re[n]
        U.im[n, n] -= temp.im[n]
        for k = n-1:-1:1
            for j = k:n
                vBj_re = U.re[k,j]
                @turbo for i = k+1:n
                    vBj_re += input.re[i,k] * U.re[i,j]
                end
                @turbo for i = k+1:n
                    vBj_re += input.im[i,k] * U.im[i,j]
                end
                vBj_im = U.im[k,j]
                @turbo for i = k+1:n
                    vBj_im += input.re[i,k] * U.im[i,j]
                end
                @turbo for i = k+1:n
                    vBj_im -= input.im[i,k] * U.re[i,j]
                end

                re = temp.re[k] * vBj_re - temp.im[k] * vBj_im
                vBj_im = temp.im[k] * vBj_re + temp.re[k] * vBj_im
                vBj_re = re

                U.re[k,j] -= vBj_re
                U.im[k,j] -= vBj_im
                @turbo for i = k+1:n
                    U.re[i,j] -= input.re[i,k] * vBj_re
                end
                @turbo for i = k+1:n
                    U.re[i,j] += input.im[i,k] * vBj_im
                end
                @turbo for i = k+1:n
                    U.im[i,j] -= input.im[i,k] * vBj_re
                end
                @turbo for i = k+1:n
                    U.im[i,j] -= input.re[i,k] * vBj_im
                end

            end
        end
    end

    @turbo for i in 1:n
        D[i] = abs(input.re[i, i])
    end

    _apply_pivot!(input, D, temp, pivot, apply_pivot)

    nothing
end

function _apply_pivot!(input::CMat64, D, temp, pivot, ::Val{true})
    # TODO: optimize
    n = size(input, 1)
    @inbounds for i in 1:n
        d = 1.0 / D[i]
        for j in 1:i-1
            temp[pivot[j]] = zero(C)
        end
        for j in i:n
            temp[pivot[j]] = d * input[i, j]
        end
        for j in 1:n
            input[i, j] = temp[j]
        end
    end
end
function _apply_pivot!(input::CMat64, D, temp, pivot, ::Val{false})
    n = size(input, 1)
    @inbounds for i in 1:n
        d = 1.0 / D[i]
        @turbo for j in i:n
            input.re[i, j] = d * input.re[i, j]
        end
    end
    @inbounds for i in 1:n
        d = 1.0 / D[i]
        @turbo for j in i:n
            input.im[i, j] = d * input.im[i, j]
        end
    end
end


text:
is_merged_pull_request(event) = get(get(event, "pull_request", Dict()), "merged", false)

function handle_merged_pull_request(event)
    number = event["pull_request"]["number"]
    @info "Processing pull request $number"
    repo, version = repo_and_version_of_pull_request_body(event["pull_request"]["body"])
    if repo === nothing
        @info "Failed to parse GitHub repository from pull request"
        return
    end
    maybe_notify(event, repo, version)
end

ext: jl


text:
# HGPold contains older implementations of the HGP recurrence relations.
# These are no longer tested or linked to, and are kept here for reference
# purposes only.
#
# These functions implement recursive versions of the integral code,
# hence the trailing "_r" in the names. These versions are not as fast
# than the versions in HGP.jl.

"vdiff(a,i,n) - Move vector a by n unit vectors in the i direction"
function vdiff(a::MVector{3,Int},i::Int,n::Int)
    b = copy(a)
    b[i] += n
    return b
end

function coulomb_hgp_r(a::PGBF,b::PGBF,c::PGBF,d::PGBF)
    return a.norm*b.norm*c.norm*d.norm*hrr_r(a.expn,a.xyz...,a.I,a.J,a.K,
        b.expn,b.xyz...,b.I,b.J,b.K,
        c.expn,c.xyz...,c.I,c.J,c.K,
        d.expn,d.xyz...,d.I,d.J,d.K)
end
coulomb_hgp_r(a::CGBF,b::CGBF,c::CGBF,d::CGBF) = contract(coulomb_hgp_r,a,b,c,d)

function hrr_r(aexpn,ax,ay,az,aI,aJ,aK,
    bexpn,bx,by,bz,bI,bJ,bK,
    cexpn,cx,cy,cz,cI,cJ,cK,
    dexpn,dx,dy,dz,dI,dJ,dK)
    if bI > 0
        return hrr_r(aexpn,ax,ay,az,aI+1,aJ,aK,bexpn,bx,by,bz,bI-1,bJ,bK,
            cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK) +
        (ax-bx)*hrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI-1,bJ,bK,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK)
    elseif bJ > 0
        return hrr_r(aexpn,ax,ay,az,aI,aJ+1,aK,bexpn,bx,by,bz,bI,bJ-1,bK,
            cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK) +
            (ay-by)*hrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ-1,bK,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK)
    elseif bK > 0
        return hrr_r(aexpn,ax,ay,az,aI,aJ,aK+1,bexpn,bx,by,bz,bI,bJ,bK-1,
            cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK) +
            (az-bz)*hrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK-1,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK)
    elseif dI > 0
        return hrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK,
            cexpn,cx,cy,cz,cI+1,cJ,cK,dexpn,dx,dy,dz,dI-1,dJ,dK) +
            (cx-dx)*hrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI-1,dJ,dK)
    elseif dJ > 0
        return hrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK,
            cexpn,cx,cy,cz,cI,cJ+1,cK,dexpn,dx,dy,dz,dI,dJ-1,dK) +
            (cy-dy)*hrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ-1,dK)
    elseif dK > 0
        return hrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK,
            cexpn,cx,cy,cz,cI,cJ,cK+1,dexpn,dx,dy,dz,dI,dJ,dK-1) +
            (cz-dz)*hrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,bI,bJ,bK,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,dI,dJ,dK-1)
    end
    return vrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
               cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,0)
end

function vrr_r(aexpn,ax,ay,az,aI,aJ,aK,
        bexpn,bx,by,bz,
        cexpn,cx,cy,cz,cI,cJ,cK,
        dexpn,dx,dy,dz,m)
    px,py,pz = gaussian_product_center(aexpn,ax,ay,az,bexpn,bx,by,bz)
    qx,qy,qz = gaussian_product_center(cexpn,cx,cy,cz,dexpn,dx,dy,dz)
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    wx,wy,wz = gaussian_product_center(zeta,px,py,pz,eta,qx,qy,qz)

    val = 0
    if cK>0
        val = (qz-cz)*vrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
            cexpn,cx,cy,cz,cI,cJ,cK-1,dexpn,dx,dy,dz,m) +
            (wz-qz)*vrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK-1,dexpn,dx,dy,dz,m+1)
        if cK>1
            val += 0.5*(cK-1)/eta*(
                vrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                    cexpn,cx,cy,cz,cI,cJ,cK-2,dexpn,dx,dy,dz,m) -
            zeta/(zeta+eta)*
                vrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                    cexpn,cx,cy,cz,cI,cJ,cK-2,dexpn,dx,dy,dz,m+1) )
        end
        if aK>0
            val += 0.5*aK/(zeta+eta)*
                vrr_r(aexpn,ax,ay,az,aI,aJ,aK-1,bexpn,bx,by,bz,
                    cexpn,cx,cy,cz,cI,cJ,cK-1,dexpn,dx,dy,dz,m+1)
        end
        return val
    elseif cJ>0
        val = (qy-cy)*vrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
            cexpn,cx,cy,cz,cI,cJ-1,cK,dexpn,dx,dy,dz,m) +
        (wy-qy)*vrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ-1,cK-1,dexpn,dx,dy,dz,m+1)
        #println("val4=$val")
        if cJ>1
            val += 0.5*(cJ-1)/eta*(
            vrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ-2,cK,dexpn,dx,dy,dz,m) -
            zeta/(zeta+eta)*
            vrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ-2,cK,dexpn,dx,dy,dz,m+1)
            )
        #println("val5=$val")
        end
        if aJ>0
            val += 0.5*aJ/(zeta+eta)*
            vrr_r(aexpn,ax,ay,az,aI,aJ-1,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ-1,cK,dexpn,dx,dy,dz,m+1)
        end
        #println("val6=$val")
        return val
    elseif cI>0
        val = (qx-cx)*vrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
            cexpn,cx,cy,cz,cI-1,cJ,cK,dexpn,dx,dy,dz,m) +
        (wx-qx)*vrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI-1,cJ,cK-1,dexpn,dx,dy,dz,m+1)
        #println("val7=$val")
        if cI>1
            val += 0.5*(cI-1)/eta*(
            vrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI-2,cJ,cK,dexpn,dx,dy,dz,m) -
            zeta/(zeta+eta)*
            vrr_r(aexpn,ax,ay,az,aI,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI-2,cJ,cK,dexpn,dx,dy,dz,m+1)
            )
        end
        #println("val8=$val")
        if aI>0
            val += 0.5*aI/(zeta+eta)*
            vrr_r(aexpn,ax,ay,az,aI-1,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI-1,cJ,cK,dexpn,dx,dy,dz,m+1)
        end
        #println("val9=$val")
        return val
    elseif aK>0
        val = (pz-az)*vrr_r(aexpn,ax,ay,az,aI,aJ,aK-1,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m) +
        (wz-pz)*vrr_r(aexpn,ax,ay,az,aI,aJ,aK-1,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)
        #println("val10=$val")
        if aK>1
            val += 0.5*(aK-1)/zeta*(
            vrr_r(aexpn,ax,ay,az,aI,aJ,aK-2,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI-1,cJ,cK,dexpn,dx,dy,dz,m) -
            eta/(zeta+eta)*
            vrr_r(aexpn,ax,ay,az,aI,aJ,aK-2,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI-1,cJ,cK,dexpn,dx,dy,dz,m+1)
            )
        end
        #println("val11=$val")
        return val
    elseif aJ>0
        val = (py-ay)*vrr_r(aexpn,ax,ay,az,aI,aJ-1,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m)+
        (wy-py)*vrr_r(aexpn,ax,ay,az,aI,aJ-1,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)
        #println("val12=$val")
        if aJ>1
            val += 0.5*(aJ-1)/zeta*(
            vrr_r(aexpn,ax,ay,az,aI,aJ-2,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m) -
            eta/(zeta+eta)*
            vrr_r(aexpn,ax,ay,az,aI,aJ-2,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)
            )
        end
        #println("val13=$val")
        return val
    elseif aI>0
        val = (px-ax)*vrr_r(aexpn,ax,ay,az,aI-1,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m) +
        (wx-px)*vrr_r(aexpn,ax,ay,az,aI-1,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)
        #println("val14=$val")
        if aI>1
            val += 0.5*(aI-1)/zeta*(
            vrr_r(aexpn,ax,ay,az,aI-2,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m) -
            eta/(zeta+eta)*
            vrr_r(aexpn,ax,ay,az,aI-2,aJ,aK,bexpn,bx,by,bz,
                cexpn,cx,cy,cz,cI,cJ,cK,dexpn,dx,dy,dz,m+1)
            )
        end
        #println("val15=$val")
        return val
    end

    rab2 = dist2(ax-bx,ay-by,az-bz)
    rcd2 = dist2(cx-dx,cy-dy,cz-dz)
    rpq2 = dist2(px-qx,py-qy,pz-qz)
    T = zeta*eta/(zeta+eta)*rpq2
    Kab = sqrt(2)pi^1.25/zeta*exp(-aexpn*bexpn*rab2/zeta)
    Kcd = sqrt(2)pi^1.25/eta*exp(-cexpn*dexpn*rcd2/eta)
    #println("rab2=$rab2,rcd2=$rcd2,rpq2=$rpq2,T=$T,Kab=$Kab,Kcd=$Kcd")
    return Kab*Kcd/sqrt(zeta+eta)*Fgamma(m,T)
end

# HGPother contains other method I implemented and tested, but these
# are either slower or less efficient with space, and should not be
# used. The fast versions are the ones in HGP2.jl

"prunem - Keep only the dictionary keys with m (last index) = 0"
prunem(d::Dict) = Dict(k[1:end-1] => v for (k,v) in d if k[end] == 0)

"vrr2 - iterative version of HGP vertical recurrence relations"
function vrr2(amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    values = Dict()
    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)
    mmax=amax+cmax

    # HGP equation 6, with b=d=0:
    #   [a+1,c]m = (Pi-Ai)[a,c]m + (Wi-Pi)[a,c]m+1
    #        + a_i/2zeta ([a-1,c]m - eta/zeta+eta[a-1,c]m+1)        # eq 6a
    #        + ci/2(zeta+eta)[a,c-1]m+1

    # First generate (0,0,0, 0,0,0, m) using eq 12
    for m in 0:mmax
        values[(0,0,0, 0,0,0, m)] = Kab*Kcd*Fgamma(m,T)/sqrt(ze)
    end

    # Now generate (ax,ay,az,0,0,0,m)
    # Eq 6a, with c=0 also:
    #   [a+1,0]m = (Pi-Ai)[a,0]m + (Wi-Pi)[a,0]m+1
    #        + a_i/2zeta ([a-1,0]m - eta/zeta+eta[a-1,0]m+1)        # eq 6b
    for ashell in 1:amax
        for ap in shell_indices[ashell]
            apx,apy,apz = ap
            i = argmax(ap) # Choose argmax(ap) as the direction to use for building new terms
            a = vdiff(ap,i,-1)
            am = vdiff(ap,i,-2)
            ax,ay,az = a
            amx,amy,amz = am
            for m in 0:(mmax-ashell)
                values[(apx,apy,apz,0,0,0,m)] = (P[i]-A[i])*values[(ax,ay,az,0,0,0,m)]+(W[i]-P[i])*values[(ax,ay,az,0,0,0,m+1)]
                if am[i] >= 0
                    values[(apx,apy,apz,0,0,0,m)] += a[i]/(2*zeta)*(values[(amx,amy,amz,0,0,0,m)]-eta/ze*values[(amx,amy,amz,0,0,0,m+1)])
                end
            end
        end
    end

    # Next build (0,0,0,cx,cy,cz,m)
    # The c-based version of 6a is:
    #   [0,c+1]m = (Qi-Bi)[0,c]m + (Wi-Qi)[0,c]m+1
    #       + ci/2eta ([0,c-1]m - zeta/zeta+eta[0,c-1]m+1)         # eq 6c
    #
    for cshell in 1:cmax
        for cp in shell_indices[cshell]
            cpx,cpy,cpz = cp
            i = argmax(cp)  # Choose argmax(cp) as the direction to use for building new terms
            c = vdiff(cp,i,-1)
            cm = vdiff(cp,i,-2)
            cx,cy,cz = c
            cmx,cmy,cmz = cm
            for m in 0:(mmax-cshell)
                values[(0,0,0,cpx,cpy,cpz,m)] = (Q[i]-C[i])*values[(0,0,0,cx,cy,cz,m)]+(W[i]-Q[i])*values[(0,0,0,cx,cy,cz,m+1)]
                if cm[i] >= 0
                    values[(0,0,0,cpx,cpy,cpz,m)] += c[i]/(2*eta)*(values[(0,0,0,cmx,cmy,cmz,m)]-zeta/ze*values[(0,0,0,cmx,cmy,cmz,m+1)])
                end
            end
        end
    end

    # Now build (ax,ay,az,cx,cy,cz,m)
    # The c-based version of 6a is:
    #   [a,c+1]m = (Qj-Bi)[a,c]m + (Wj-Qj)[a,c]m+1
    #       + c_j/2eta ([a,c-1]m - zeta/zeta+eta[a,c-1]m+1)         # eq 6d
    #       + a_j/2(zeta+eta)[a-1,c]m+1
    for ashell in 1:amax
        for a in shell_indices[ashell]
            ax,ay,az = a
            for cshell in 1:cmax
                for cp in shell_indices[cshell]
                    cpx,cpy,cpz = cp
                    j = argmax(cp)  # Choose argmax(cp) as the direction to use for building new terms
                    c = vdiff(cp,j,-1)
                    cm = vdiff(cp,j,-2)
                    cx,cy,cz = c
                    cmx,cmy,cmz = cm

                    am = vdiff(a,j,-1)
                    amx,amy,amz = am
                    for m in 0:(mmax-ashell-cshell)
                        values[(ax,ay,az,cpx,cpy,cpz,m)] = (Q[j]-C[j])*values[(ax,ay,az,cx,cy,cz,m)]+(W[j]-Q[j])*values[(ax,ay,az,cx,cy,cz,m+1)]
                        if cm[j] >= 0
                            values[(ax,ay,az,cpx,cpy,cpz,m)] += c[j]/(2*eta)*(values[(ax,ay,az,cmx,cmy,cmz,m)]-zeta/ze*values[(ax,ay,az,cmx,cmy,cmz,m+1)])
                        end
                        if am[j] >= 0
                            values[(ax,ay,az,cpx,cpy,cpz,m)] += a[j]/(2*ze)*(values[(amx,amy,amz,cx,cy,cz,m+1)])
                        end
                    end
                end
            end
        end
    end
    return prunem(values)
end

"""
hrr - compute two electron integrals using the horizontal recurrence relations

The HRRs are given in HGP eq 18:
    (a(b+1),cd) = ((a+1)b,cd) + (Ai-Bi)(ab,cd)

Note that the HRRs apply to either contracted or primitive integrals.
"""
function hrr2(ashell,bshell,cshell,dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D)

    # Get the relevant vrr terms.
    # Interesting that the Gaussian exponents are simply a pass-through to vrr.
    vrrs = vrr2(ashell+bshell,cshell+dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D)

    values = Dict()
    # Put vrr values into the hrr values dictionary
    for (ax,ay,az,cx,cy,cz) in keys(vrrs)
        values[(ax,ay,az,0,0,0,cx,cy,cz,0,0,0)] = vrrs[(ax,ay,az,cx,cy,cz)]
    end

    # First build (ab,c0) from (a0,c0)
    for bs in 1:bshell
        for bp in shell_indices[bs]
            bpx,bpy,bpz = bp
            j = argmax(bp)
            bx,by,bz = vdiff(bp,j,-1)
            for as in ashell:(ashell+bshell-bs) # -1 guess
                for a in shell_indices[as]
                    ax,ay,az = a
                    apx,apy,apz = vdiff(a,j,1)
                    for cs in 0:(cshell+dshell)
                        for c in shell_indices[cs]
                            cx,cy,cz = c
                            values[(ax,ay,az,bpx,bpy,bpz,cx,cy,cz,0,0,0)] = values[(apx,apy,apz,bx,by,bz,cx,cy,cz,0,0,0)] +
                                (A[j]-B[j])*values[(ax,ay,az,bx,by,bz,cx,cy,cz,0,0,0)]
                        end
                    end
                end
            end
        end
    end
    # now build (ab,cd) from (ab,c0)
    for ds in 1:dshell
        for dp in shell_indices[ds]
            dpx,dpy,dpz = dp
            j = argmax(dp)
            dx,dy,dz = vdiff(dp,j,-1)
            for cs in cshell:(cshell+dshell-ds)
                for c in shell_indices[cs]
                    cx,cy,cz = c
                    cpx,cpy,cpz = vdiff(c,j,1)
                    for a in shell_indices[ashell]
                        ax,ay,az = a
                        for b in shell_indices[bshell]
                            bx,by,bz = b
                            values[(ax,ay,az,bx,by,bz,cx,cy,cz,dpx,dpy,dpz)] = values[(ax,ay,az,bx,by,bz,cpx,cpy,cpz,dx,dy,dz)] +
                                (C[j]-D[j])*values[(ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz)]
                        end
                    end
                end
            end
        end
    end
    # We could also just return the subset of values where (b) = bshell and (d) = dshell.
    # But we've already calculated them
    return values
end

"hrr3 - hrr using arrays rather than dicts.
This method appears to give the right answer, but it's a very wasteful method of
storing integrals, because of the 12-d array.
"
function hrr3(ashell,bshell,cshell,dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D)

    # Get the relevant vrr terms.
    vrrs = vrr_widearray(ashell+bshell,cshell+dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D)

    alim = ashell+bshell
    blim = bshell
    clim = cshell+dshell
    dlim = dshell
    values = OffsetArray(zeros(Float64,alim+1,alim+1,alim+1,blim+1,blim+1,blim+1,clim+1,clim+1,clim+1,dlim+1,dlim+1,dlim+1),
                            0:alim,0:alim,0:alim,0:blim,0:blim,0:blim,0:clim,0:clim,0:clim,0:dlim,0:dlim,0:dlim)
    # Put vrr values into the hrr values dictionary
    values[:,:,:,0,0,0,:,:,:,0,0,0] = vrrs[:,:,:,:,:,:]

    # First build (ab,c0) from (a0,c0)
    for bs in 1:bshell
        for bp in shell_indices[bs]
            bpx,bpy,bpz = bp
            j = argmax(bp)
            bx,by,bz = vdiff(bp,j,-1)
            for as in ashell:(ashell+bshell-bs) # -1 guess
                for a in shell_indices[as]
                    ax,ay,az = a
                    apx,apy,apz = vdiff(a,j,1)
                    for cs in 0:(cshell+dshell)
                        for c in shell_indices[cs]
                            cx,cy,cz = c
                            values[ax,ay,az,bpx,bpy,bpz,cx,cy,cz,0,0,0] = values[apx,apy,apz,bx,by,bz,cx,cy,cz,0,0,0] +
                                (A[j]-B[j])*values[ax,ay,az,bx,by,bz,cx,cy,cz,0,0,0]
                        end
                    end
                end
            end
        end
    end
    # now build (ab,cd) from (ab,c0)
    for ds in 1:dshell
        for dp in shell_indices[ds]
            dpx,dpy,dpz = dp
            j = argmax(dp)
            dx,dy,dz = vdiff(dp,j,-1)
            for cs in cshell:(cshell+dshell-ds)
                for c in shell_indices[cs]
                    cx,cy,cz = c
                    cpx,cpy,cpz = vdiff(c,j,1)
                    for a in shell_indices[ashell]
                        ax,ay,az = a
                        for b in shell_indices[bshell]
                            bx,by,bz = b
                            values[ax,ay,az,bx,by,bz,cx,cy,cz,dpx,dpy,dpz] = values[ax,ay,az,bx,by,bz,cpx,cpy,cpz,dx,dy,dz] +
                                (C[j]-D[j])*values[ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz]
                        end
                    end
                end
            end
        end
    end
    # We could also just return the subset of values where (b) = bshell and (d) = dshell.
    # But we've already calculated them
    return values
end
"""
vrr_widearray(amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)

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

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.
`amax` and `cmax` are related to the sum of the shell angular
momenta for the `a+b`, and `c+d` shells, respectively.

The function returns a six-dimensional array over the possible
powers of the `a+b` and `c+d` shell functions.
"""
function vrr_widearray(amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax=amax+cmax
    vrrs = OffsetArray(zeros(Float64,amax+1,amax+1,amax+1,cmax+1,cmax+1,cmax+1,mmax+1),
         0:amax,0:amax,0:amax, 0:cmax,0:cmax,0:cmax,0: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)

    # HGP equation 6, with b=d=0:
    #   [a+1,c]m = (Pi-Ai)[a,c]m + (Wi-Pi)[a,c]m+1
    #        + a_i/2zeta ([a-1,c]m - eta/zeta+eta[a-1,c]m+1)        # eq 6a
    #        + ci/2(zeta+eta)[a,c-1]m+1

    # First generate (0,0,0, 0,0,0, m) using eq 12
    #vrrs[0,0,0, 0,0,0,0:mmax] = Kab*Kcd/sqrt(ze)*Fgamma.(0:mmax,T)
    for m in 0:mmax
        vrrs[0,0,0, 0,0,0, m] = Kab*Kcd*Fgamma(m,T)/sqrt(ze)
    end

    # Now generate (ax,ay,az,0,0,0,m)
    # Eq 6a, with c=0 also:
    #   [a+1,0]m = (Pi-Ai)[a,0]m + (Wi-Pi)[a,0]m+1
    #        + a_i/2zeta ([a-1,0]m - eta/zeta+eta[a-1,0]m+1)        # eq 6b
    for ashell in 1:amax
        for ap in shell_indices[ashell]
            apx,apy,apz = ap
            i = argmax(ap) # Choose argmax(ap) as the direction to use for building new terms
            a = vdiff(ap,i,-1)
            am = vdiff(ap,i,-2)
            ax,ay,az = a
            amx,amy,amz = am
            for m in 0:(mmax-ashell)
                vrrs[apx,apy,apz,0,0,0,m] = (P[i]-A[i])*vrrs[ax,ay,az,0,0,0,m]+(W[i]-P[i])*vrrs[ax,ay,az,0,0,0,m+1]
            end
            if am[i] >= 0
                for m in 0:(mmax-ashell)
                    vrrs[apx,apy,apz,0,0,0,m] += a[i]/(2*zeta)*(vrrs[amx,amy,amz,0,0,0,m]-eta/ze*vrrs[amx,amy,amz,0,0,0,m+1])
                end
            end
        end
    end

    # Next build (0,0,0,cx,cy,cz,m)
    # The c-based version of 6a is:
    #   [0,c+1]m = (Qi-Bi)[0,c]m + (Wi-Qi)[0,c]m+1
    #       + ci/2eta ([0,c-1]m - zeta/zeta+eta[0,c-1]m+1)         # eq 6c
    #
    for cshell in 1:cmax
        for cp in shell_indices[cshell]
            cpx,cpy,cpz = cp
            i = argmax(cp)  # Choose argmax(cp) as the direction to use for building new terms
            c = vdiff(cp,i,-1)
            cm = vdiff(cp,i,-2)
            cx,cy,cz = c
            cmx,cmy,cmz = cm
            for m in 0:(mmax-cshell)
                vrrs[0,0,0,cpx,cpy,cpz,m] = (Q[i]-C[i])*vrrs[0,0,0,cx,cy,cz,m]+(W[i]-Q[i])*vrrs[0,0,0,cx,cy,cz,m+1]
            end
            if cm[i] >= 0
                for m in 0:(mmax-cshell)
                    vrrs[0,0,0,cpx,cpy,cpz,m] += c[i]/(2*eta)*(vrrs[0,0,0,cmx,cmy,cmz,m]-zeta/ze*vrrs[0,0,0,cmx,cmy,cmz,m+1])
                end
            end
        end
    end

    # Now build (ax,ay,az,cx,cy,cz,m)
    # The c-based version of 6a is:
    #   [a,c+1]m = (Qj-Bi)[a,c]m + (Wj-Qj)[a,c]m+1
    #       + c_j/2eta ([a,c-1]m - zeta/zeta+eta[a,c-1]m+1)         # eq 6d
    #       + a_j/2(zeta+eta)[a-1,c]m+1
    for ashell in 1:amax
        for a in shell_indices[ashell]
            ax,ay,az = a
            for cshell in 1:cmax
                for cp in shell_indices[cshell]
                    cpx,cpy,cpz = cp
                    j = argmax(cp)  # Choose argmax(cp) as the direction to use for building new terms
                    c = vdiff(cp,j,-1)
                    cm = vdiff(cp,j,-2)
                    cx,cy,cz = c
                    cmx,cmy,cmz = cm

                    am = vdiff(a,j,-1)
                    amx,amy,amz = am
                    for m in 0:(mmax-ashell-cshell)
                        vrrs[ax,ay,az,cpx,cpy,cpz,m] = (Q[j]-C[j])*vrrs[ax,ay,az,cx,cy,cz,m]+(W[j]-Q[j])*vrrs[ax,ay,az,cx,cy,cz,m+1]
                    end
                    if cm[j] >= 0
                        for m in 0:(mmax-ashell-cshell)
                            vrrs[ax,ay,az,cpx,cpy,cpz,m] += c[j]/(2*eta)*(vrrs[ax,ay,az,cmx,cmy,cmz,m]-zeta/ze*vrrs[ax,ay,az,cmx,cmy,cmz,m+1])
                        end
                    end
                    if am[j] >= 0
                        for m in 0:(mmax-ashell-cshell)
                            vrrs[ax,ay,az,cpx,cpy,cpz,m] += a[j]/(2*ze)*(vrrs[amx,amy,amz,cx,cy,cz,m+1])
                        end
                    end
                end
            end
        end
    end
    return vrrs[:,:,:,:,:,:,0]
end
function hrr_shells(ashell,bshell,cshell,dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    # Get the relevant vrr terms.
     vrrs = vrr(ashell+bshell,cshell+dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
     hrrs = zeros(Float64,nao[ashell+bshell],nao[bshell],nao[cshell+dshell],nao[dshell])
     hrrs[:,1,:,1] = vrrs[:,:]

     # First build (ab,c0) from (a0,c0)
     for bs in 1:bshell
         for bp in shell_indices[bs]
             bpindex = m2ao[bp]
             j = argmax(bp)
             b = vdiff(bp,j,-1)
             bindex = m2ao[b]
             for as in 0:(ashell+bshell-bs)
                 for a in shell_indices[as]
                     aindex = m2ao[a]
                     ap = vdiff(a,j,1)
                     apindex = m2ao[ap]
                     for cs in 0:(cshell+dshell)
                         for c in shell_indices[cs]
                             cindex = m2ao[c]
                             hrrs[aindex,bpindex,cindex,1] = hrrs[apindex,bindex,cindex,1] +
                                 (A[j]-B[j])*hrrs[aindex,bindex,cindex,1]
                         end
                     end
                 end
             end
         end
     end
     # now build (ab,cd) from (ab,c0)
     for ds in 1:dshell
         for dp in shell_indices[ds]
             dpindex = m2ao[dp]
             j = argmax(dp)
             d = vdiff(dp,j,-1)
             dindex = m2ao[d]
             for cs in 0:(cshell+dshell-ds)
                 for c in shell_indices[cs]
                     cindex = m2ao[c]
                     cp = vdiff(c,j,1)
                     cpindex = m2ao[cp]
                     for as in 0:ashell
                         for a in shell_indices[as]
                             aindex = m2ao[a]
                             for bs in 0:bshell
                                 for b in shell_indices[bs]
                                     bindex = m2ao[b]
                                     hrrs[aindex,bindex,cindex,dpindex] = hrrs[aindex,bindex,cpindex,dindex] +
                                         (C[j]-D[j])*hrrs[aindex,bindex,cindex,dindex]
                                 end
                             end
                         end
                     end
                 end
             end
         end
     end
     return hrrs[1:nao[ashell],:,1:nao[cshell],:]
 end

"""
hrr_dict(ashell,bshell,cshell,dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's horizontal 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.
`ashell`, `bshell`, `cshell` and `dshell` are the shell angular
momenta for the `a`, `b`, `c`, and `d` shells, respectively.

The function returns a dictionary containing entries for the
relevant integrals. E.g., `hrrs[ax,ay,az,bx,by,bz,cx,cy,cz,dpx,dpy,dpz]`
contains the integral corresponding to the bfs with powers `ax,ay,az`,
`bx,by,bz`, `cx,cy,cz`, `dx,dy,dz`.

`hrr_dict` is slower than `hrr`, but is kept for convenience.
"""
function hrr_dict(ashell,bshell,cshell,dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    vrrs = vrr_widearray(ashell+bshell,cshell+dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    hrrs = Dict{NTuple{12,Int},Float64}()

    for as in 0:(ashell+bshell)
        for (ax,ay,az) in shell_indices[as]
            for cs in 0:(cshell+dshell)
                for (cx,cy,cz) in shell_indices[cs]
                    hrrs[ax,ay,az,0,0,0,cx,cy,cz,0,0,0] = vrrs[ax,ay,az,cx,cy,cz]
                end
            end
        end
    end

    # First build (ab,c0) from (a0,c0)
    for bs in 1:bshell
        for bp in shell_indices[bs]
            bpx,bpy,bpz = bp
            j = argmax(bp)
            bx,by,bz = vdiff(bp,j,-1)
            for as in 0:(ashell+bshell-bs)
                for a in shell_indices[as]
                    ax,ay,az = a
                    apx,apy,apz = vdiff(a,j,1)
                    for cs in 0:(cshell+dshell)
                        for c in shell_indices[cs]
                            cx,cy,cz = c
                            hrrs[ax,ay,az,bpx,bpy,bpz,cx,cy,cz,0,0,0] = hrrs[apx,apy,apz,bx,by,bz,cx,cy,cz,0,0,0] +
                                (A[j]-B[j])*hrrs[ax,ay,az,bx,by,bz,cx,cy,cz,0,0,0]
                        end
                    end
                end
            end
        end
    end

    # now build (ab,cd) from (ab,c0)
    for ds in 1:dshell
        for dp in shell_indices[ds]
            dpx,dpy,dpz = dp
            j = argmax(dp)
            dx,dy,dz = vdiff(dp,j,-1)
            for cs in 0:(cshell+dshell-ds)
                for c in shell_indices[cs]
                    cx,cy,cz = c
                    cpx,cpy,cpz = vdiff(c,j,1)
                    for as in 0:ashell
                        for a in shell_indices[as]
                            ax,ay,az = a
                            for bs in 0:bshell
                                for b in shell_indices[bs]
                                    bx,by,bz = b
                                    hrrs[ax,ay,az,bx,by,bz,cx,cy,cz,dpx,dpy,dpz] = hrrs[ax,ay,az,bx,by,bz,cpx,cpy,cpz,dx,dy,dz] +
                                        (C[j]-D[j])*hrrs[ax,ay,az,bx,by,bz,cx,cy,cz,dx,dy,dz]
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return hrrs
end
function vrr_shells(amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax=amax+cmax+1
    vrrs = zeros(Float64,nao[amax],nao[cmax],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)

    # HGP equation 6, with b=d=0:
    #   [a+1,c]m = (Pi-Ai)[a,c]m + (Wi-Pi)[a,c]m+1
    #        + a_i/2zeta ([a-1,c]m - eta/zeta+eta[a-1,c]m+1)        # eq 6a
    #        + ci/2(zeta+eta)[a,c-1]m+1

    # First generate (1,1,m) using eq 12
    for m in 1:(mmax)
        vrrs[1,1, m] = Kab*Kcd*Fgamma(m-1,T)/sqrt(ze)
    end

    # Generate (A,1,m)
    # Eq 6a, with c=0 also:
    #   [a+1,0]m = (Pi-Ai)[a,1]m + (Wi-Pi)[a,1]m+1
    #        + a_i/2zeta ([a-1,0]m - eta/zeta+eta[a-1,0]m+1)        # eq 6b

    for ashell in 1:amax
        for ap in shell_indices[ashell]
            apindex = m2ao[ap]
            i = argmax(ap) # Choose argmax(ap) as the direction to use for building new terms
            a = vdiff(ap,i,-1)
            aindex = m2ao[a]
            am = vdiff(ap,i,-2)
            for m in 1:(mmax-ashell)
                vrrs[apindex,1,m] = (P[i]-A[i])*vrrs[aindex,1,m]+(W[i]-P[i])*vrrs[aindex,1,m+1]
            end
            if am[i] >= 0
                amindex = m2ao[am]
                for m in 1:(mmax-ashell)
                    vrrs[apindex,1,m] += a[i]/(2*zeta)*(vrrs[amindex,1,m]-eta/ze*vrrs[amindex,1,m+1])
                end
            end
        end
    end

    # Next build (1,C,m)
    # The c-based version of 6a is:
    #   [0,c+1]m = (Qi-Bi)[1,c]m + (Wi-Qi)[1,c]m+1
    #       + ci/2eta ([1,c-1]m - zeta/zeta+eta[1,c-1]m+1)         # eq 6c
    #
    for cshell in 1:cmax
        for cp in shell_indices[cshell]
            cpindex = m2ao[cp]
            i = argmax(cp)  # Choose argmax(cp) as the direction to use for building new terms
            c = vdiff(cp,i,-1)
            cindex = m2ao[c]
            cm = vdiff(cp,i,-2)
           for m in 1:(mmax-cshell)
                vrrs[1,cpindex,m] = (Q[i]-C[i])*vrrs[1,cindex,m]+(W[i]-Q[i])*vrrs[1,cindex,m+1]
            end
            if cm[i] >= 0
                cmindex = m2ao[cm]
                for m in 1:(mmax-cshell)
                    vrrs[1,cpindex,m] += c[i]/(2*eta)*(vrrs[1,cmindex,m]-zeta/ze*vrrs[1,cmindex,m+1])
                end
            end
        end
    end

    # Now build (A,C,m)
    # The c-based version of 6a is:
    #   [a,c+1]m = (Qj-Bi)[a,c]m + (Wj-Qj)[a,c]m+1
    #       + c_j/2eta ([a,c-1]m - zeta/zeta+eta[a,c-1]m+1)         # eq 6d
    #       + a_j/2(zeta+eta)[a-1,c]m+1
    for ashell in 1:amax
        for a in shell_indices[ashell]
            aindex = m2ao[a]
            for cshell in 1:cmax
                for cp in shell_indices[cshell]
                    cpindex = m2ao[cp]
                    j = argmax(cp)  # Choose argmax(cp) as the direction to use for building new terms
                    c = vdiff(cp,j,-1)
                    cindex = m2ao[c]
                    cm = vdiff(cp,j,-2)
                    am = vdiff(a,j,-1)
                    for m in 1:(mmax-ashell-cshell)
                        vrrs[aindex,cpindex,m] = (Q[j]-C[j])*vrrs[aindex,cindex,m]+(W[j]-Q[j])*vrrs[aindex,cindex,m+1]
                    end
                    if cm[j] >= 0
                        cmindex = m2ao[cm]
                        for m in 1:(mmax-ashell-cshell)
                            vrrs[aindex,cpindex,m] += c[j]/(2*eta)*(vrrs[aindex,cmindex,m]-zeta/ze*vrrs[aindex,cmindex,m+1])
                        end
                    end
                    if am[j] >= 0
                        amindex = m2ao[am]
                        for m in 1:(mmax-ashell-cshell)
                            vrrs[aindex,cpindex,m] += a[j]/(2*ze)*(vrrs[amindex,cindex,m+1])
                        end
                    end
                end
            end
        end
    end
    return vrrs[:,:,1]
end

function chrr_shell(ash::Shell,bsh::Shell,csh::Shell,dsh::Shell)
    ashell,bshell,cshell,dshell = ash.L,bsh.L,csh.L,dsh.L
    A,B,C,D = ash.xyz,bsh.xyz,csh.xyz,dsh.xyz
    vrrs = cvrr(ash,bsh,csh,dsh)
    hrrs = zeros(Float64,nao[ashell+bshell],nao[bshell],nao[cshell+dshell],nao[dshell])
    hrrs[:,1,:,1] = vrrs[:,:]

    # First build (ab,c0) from (a0,c0)
    for bs in 1:bshell
        for bp in shell_indices[bs]
            bpindex = m2ao[bp]
            j = argmax(bp)
            b = vdiff(bp,j,-1)
            bindex = m2ao[b]
            for as in 0:(ashell+bshell-bs)
                for a in shell_indices[as]
                    aindex = m2ao[a]
                    ap = vdiff(a,j,1)
                    apindex = m2ao[ap]
                    for cs in 0:(cshell+dshell)
                        for c in shell_indices[cs]
                            cindex = m2ao[c]
                            hrrs[aindex,bpindex,cindex,1] = hrrs[apindex,bindex,cindex,1] +
                                (A[j]-B[j])*hrrs[aindex,bindex,cindex,1]
                        end
                    end
                end
            end
        end
    end
    # now build (ab,cd) from (ab,c0)
    for ds in 1:dshell
        for dp in shell_indices[ds]
            dpindex = m2ao[dp]
            j = argmax(dp)
            d = vdiff(dp,j,-1)
            dindex = m2ao[d]
            for cs in 0:(cshell+dshell-ds)
                for c in shell_indices[cs]
                    cindex = m2ao[c]
                    cp = vdiff(c,j,1)
                    cpindex = m2ao[cp]
                    for as in 0:ashell
                        for a in shell_indices[as]
                            aindex = m2ao[a]
                            for bs in 0:bshell
                                for b in shell_indices[bs]
                                    bindex = m2ao[b]
                                    hrrs[aindex,bindex,cindex,dpindex] = hrrs[aindex,bindex,cpindex,dindex] +
                                        (C[j]-D[j])*hrrs[aindex,bindex,cindex,dindex]
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    return hrrs#[1:nao[ashell],:,1:nao[cshell],:]
end

#= Works in progress:
function vrr_shell_recursive(amax,cmax,mmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
   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)
   if amax+cmax == 0
       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)
       vrrs = zeros(Float64,1,1,mmax)
       for m in 1:mmax
           vrrs[1,1,m] = Kab*Kcd/sqrt(ze)*Fgamma(m-1,T)
       end
       return vrrs
   end
   if cmax > 0
       vrrsm = vrr_shell_recursive(amax,cmax-1,mmax+1, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
       # (a,c+dir) = shift(dir)(a,c)
       # if exists(a-dir): (a,c+dir) += (a-dir,c)
       # if exists(c-dir): (a,c+dir) += (a,c-dir)
   return vrrs
end
function vrr_shell_dict(amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
   mmax = amax+cmax+1
   vrrs = Dict{NTuple{Int,2},Matrix{Float64}}
   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

   # HGP equation 6, with b=d=0:
   #   [a+1,c]m = (Pi-Ai)[a,c]m + (Wi-Pi)[a,c]m+1
   #        + a_i/2zeta ([a-1,c]m - eta/zeta+eta[a-1,c]m+1)        # eq 6a
   #        + ci/2(zeta+eta)[a,c-1]m+1
   # First generate (1,1,m) using eq 12
   vrrs[0,0] = zeros(Float64,1,1,mmax)
   for m in 1:mmax
       vrrs[0,0][1,1,m] = Kab*Kcd/sqrt(ze)*Fgamma(m-1,T)
   end
   for ashell in 1:amax
       na = length(shell_indices[ashell])
       nc = 1
       ma = mmax-ashell
       vrrs[1,0] = zeros(Float64,na,nc,ma)
       ic = nc
       for iap = 1:na
           ap = shell_indices[iap]
           dir = argmax(ap)
           ia = shell_minus(iap,dir)
           for m in 1:ma
               vrrs[ashell,1][iap,ic,m] = PA[dir]*vrrs[ashell-1,0][ia,ic,m]+ WP[dir]vrrs[ashell-1,0][ia,ic,m+1]
           end
           iam = shell_minus(ia,dir)
           if exists(iam)
               for m in 1:ma
                   vrrs[ashell,1][iap,ic,m] += index_values2(ia,ashell-1,dir)*(vrrs[ashell-2,0][iam,ic,m] -eta/ze*vrrs[ashell-2,0][iam,ic,m+1])
               end
           end
       end
   end
   # Next build (0,0,0,cx,cy,cz,m)
   # The c-based version of 6a is:
   #   [0,c+1]m = (Qi-Bi)[0,c]m + (Wi-Qi)[0,c]m+1
   #       + ci/2eta ([0,c-1]m - zeta/zeta+eta[0,c-1]m+1)         # eq 6c
   #
   for cshell in 1:cmax
       for cp in shell_indices[cshell]
           cpx,cpy,cpz = cp
           i = argmax(cp)  # Choose argmax(cp) as the direction to use for building new terms
           c = vdiff(cp,i,-1)
           cm = vdiff(cp,i,-2)
           cx,cy,cz = c
           cmx,cmy,cmz = cm
           for m in 0:(mmax-cshell)
               values[(0,0,0,cpx,cpy,cpz,m)] = QC[i]*values[(0,0,0,cx,cy,cz,m)]+(W[i]-Q[i])*values[(0,0,0,cx,cy,cz,m+1)]
               if cm[i] >= 0
                   values[(0,0,0,cpx,cpy,cpz,m)] += c[i]/(2*eta)*(values[(0,0,0,cmx,cmy,cmz,m)]-zeta/ze*values[(0,0,0,cmx,cmy,cmz,m+1)])
               end
           end
       end
   end
   # Now build (ax,ay,az,cx,cy,cz,m)
   # The c-based version of 6a is:
   #   [a,c+1]m = (Qj-Bi)[a,c]m + (Wj-Qj)[a,c]m+1
   #       + c_j/2eta ([a,c-1]m - zeta/zeta+eta[a,c-1]m+1)         # eq 6d
   #       + a_j/2(zeta+eta)[a-1,c]m+1
   for ashell in 1:amax
       for a in shell_indices[ashell]
           ax,ay,az = a
           for cshell in 1:cmax
               for cp in shell_indices[cshell]
                   cpx,cpy,cpz = cp
                   j = argmax(cp)  # Choose argmax(cp) as the direction to use for building new terms
                   c = vdiff(cp,j,-1)
                   cm = vdiff(cp,j,-2)
                   cx,cy,cz = c
                   cmx,cmy,cmz = cm
                   am = vdiff(a,j,-1)
                   amx,amy,amz = am
                   for m in 0:(mmax-ashell-cshell)
                       values[(ax,ay,az,cpx,cpy,cpz,m)] = (Q[j]-C[j])*values[(ax,ay,az,cx,cy,cz,m)]+(W[j]-Q[j])*values[(ax,ay,az,cx,cy,cz,m+1)]
                       if cm[j] >= 0
                           values[(ax,ay,az,cpx,cpy,cpz,m)] += c[j]/(2*eta)*(values[(ax,ay,az,cmx,cmy,cmz,m)]-zeta/ze*values[(ax,ay,az,cmx,cmy,cmz,m+1)])
                       end
                       if am[j] >= 0
                           values[(ax,ay,az,cpx,cpy,cpz,m)] += a[j]/(2*ze)*(values[(amx,amy,amz,cx,cy,cz,m+1)])
                       end
                   end
               end
           end
       end
   end
   return prunem(values)
end
=#


text:
# Sample code to test code
using Random
using CSV
using DataFrames
using StatsBase

Random.seed!(28)
include("../cdfJulia\\CDFJulia.jl")
include("utilities.jl")
sp= CSV.read("spyData.csv", DataFrame)
rs = getrnd(sp.R, length(sp.R))

nlag = 50
n2lag = 5*nlag
aR = autocor(sp.R, 1:nlag)
bR = autocor(abs.(sp.R), 1:n2lag)
cR = crosscor(sp.R, abs.(sp.R), -nlag:nlag)
bR2 = autocor((sp.R).^2, 1:n2lag)

cmin = efun6(aR, bR, cR, bR2, rs, nlag, n2lag)
cgoal = 1e-3

function wombo(aR,bR,cR,bR2,cmin,rs0,nlag,n2lag,cgoal)
  temp = 1e-4
  maxT = 20000
  maxG = 2000
  maxS = 200
  wc   = 0.9
  afac = 0.5
  maxA = 120e6
  rs   = copy(rs0)
  cost = copy(cmin)
  nTot  = 0
  nGood = 0
  nAll  = 0
  ra    = rs
  tini  = temp
  E     = zeros(convert(Int64,ceil(maxA/1e6)),1)
  tE    = copy(cmin)
  for tt=1:maxA
    nTot = nTot+1
    nAll = nAll+1
    if (nAll%100000) == 0
      println(nAll/1e6)
      E[convert(Int64,ceil(nAll/1e6))] = tE
    end
    cmax = cost-temp*log(rand())
    trs  = pfun(rs)
    tE   = efun6(aR,bR,cR,bR2,trs,nlag,n2lag)
    if tE<=cmax
      nGood = nGood+1
      rs    = trs
      cost  = tE
    end
    if cost<=cgoal
      ra = rs
      break
    end
  if ((nTot<maxT) && (nGood<maxG))
    if (cost<(cmin*wc))
      cmin = cost
      ra   = rs
    end
    continue
  end
  if ((temp==tini) && (nTot>1.5*nGood))
    tini = 10*temp
    temp = tini
  else
    if (nGood<=maxS)
      afac=sqrt(afac)
      maxT=maxT*sqrt(2)
      temp=tini
    else
      temp=temp*afac
    end
  end
  nTot  = 0
  nGood = 0
end
return ra,E
end

ra,E = wombo(aR,bR,cR,bR2,cmin,rs,nlag,n2lag,cgoal)

max_stars_repo_name: silvaac/SMCJulia

max_issues_repo_name: silvaac/SMCJulia

max_forks_repo_name: silvaac/SMCJulia


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# RUN: julia --startup-file=no %s %t && llvm-link -S %t/* -o %t/module.ll
# RUN: cat %t/module.ll | FileCheck %s
# RUN: cat %t/module.ll | opt -load libjulia%shlibext -LowerSIMDLoop -S - | FileCheck %s -check-prefix=LOWER
# RUN: julia --startup-file=no %s %t -O && llvm-link -S %t/* -o %t/module.ll
# RUN: cat %t/module.ll | FileCheck %s -check-prefix=FINAL

## Notes:
# This script uses the `emit` function (defined llvmpasses.jl) to emit either
# optimized or unoptimized LLVM IR. Each function is emitted individually and
# `llvm-link` is used to create a single module that can be passed to opt.
# The order in which files are emitted and linked is important since `lit` will
# process the test cases in order.
#
# There are three different test prefixes defined:
# - `CHECK`: Checks the result of codegen
# - `LOWER`: Checks the result of -LowerSIMDLoop
# - `FINAL`: Checks the result of running the entire pipeline
include(joinpath("..", "testhelpers", "llvmpasses.jl"))

# CHECK-LABEL: @julia_simdf_
# LOWER-LABEL: @julia_simdf_
# FINAL-LABEL: @julia_simdf_
function simdf(X)
    acc = zero(eltype(X))
    @simd for x in X
        acc += x
# CHECK: call void @julia.loopinfo_marker(), {{.*}}, !julia.loopinfo [[LOOPINFO:![0-9]+]]
# LOWER-NOT: llvm.mem.parallel_loop_access
# LOWER: fadd fast double
# LOWER-NOT: call void @julia.loopinfo_marker()
# LOWER: br {{.*}}, !llvm.loop [[LOOPID:![0-9]+]]
# FINAL: fadd fast <{{[0-9]+}} x double>
    end
    acc
end

# CHECK-LABEL: @julia_simdf2_
# LOWER-LABEL: @julia_simdf2_
function simdf2(X)
    acc = zero(eltype(X))
    @simd ivdep for x in X
        acc += x
# CHECK: call void @julia.loopinfo_marker(), {{.*}}, !julia.loopinfo [[LOOPINFO2:![0-9]+]]
# LOWER: llvm.mem.parallel_loop_access
# LOWER-NOT: call void @julia.loopinfo_marker()
# LOWER: fadd fast double
# LOWER: br {{.*}}, !llvm.loop [[LOOPID2:![0-9]+]]
    end
    acc
end

@noinline iteration(i) = (@show(i); return nothing)

# CHECK-LABEL: @julia_loop_unroll
# LOWER-LABEL: @julia_loop_unroll
# FINAL-LABEL: @julia_loop_unroll
@eval function loop_unroll(N)
    for i in 1:N
        iteration(i)
        $(Expr(:loopinfo, (Symbol("llvm.loop.unroll.count"), 3)))
# CHECK: call void @julia.loopinfo_marker(), {{.*}}, !julia.loopinfo [[LOOPINFO3:![0-9]+]]
# LOWER-NOT: call void @julia.loopinfo_marker()
# LOWER: br {{.*}}, !llvm.loop [[LOOPID3:![0-9]+]]
# FINAL: call void @j_iteration
# FINAL: call void @j_iteration
# FINAL: call void @j_iteration
# FINAL-NOT: call void @j_iteration
# FINAL: br
    end
end

# Example from a GPU kernel where we want to unroll the outer loop
# and the inner loop is a boundschecked single iteration loop.
# The `@show` is used to bloat the loop and `X-COUNT-10:` seems
# not to be working so we duplicate the checks. FIXME LLVM8
# CHECK-LABEL: @julia_loop_unroll2
# LOWER-LABEL: @julia_loop_unroll2
# FINAL-LABEL: @julia_loop_unroll2
@eval function loop_unroll2(J, I)
    for i in 1:10
        for j in J
            1 <= j <= I && continue
            @show (i,j)
            iteration(i)
        end
        $(Expr(:loopinfo, (Symbol("llvm.loop.unroll.full"),)))
# CHECK: call void @julia.loopinfo_marker(), {{.*}}, !julia.loopinfo [[LOOPINFO4:![0-9]+]]
# LOWER-NOT: call void @julia.loopinfo_marker()
# LOWER: br {{.*}}, !llvm.loop [[LOOPID4:![0-9]+]]
# FINAL: call void @j_iteration
# FINAL: call void @j_iteration
# FINAL: call void @j_iteration
# FINAL: call void @j_iteration
# FINAL: call void @j_iteration
# FINAL: call void @j_iteration
# FINAL: call void @j_iteration
# FINAL: call void @j_iteration
# FINAL: call void @j_iteration
# FINAL: call void @j_iteration
# FINAL-NOT: call void @j_iteration
    end
end

# FINAL-LABEL: @julia_notunroll
function notunroll(J, I)
    for i in 1:10
        for j in J
            1 <= j <= I && continue
            @show (i,j)
            iteration(i)
# FINAL: call void @j_iteration
# FINAL-NOT: call void @j_iteration
        end
    end
end

## Check all the MD nodes
# CHECK: [[LOOPINFO]] = !{!"julia.simdloop"}
# CHECK: [[LOOPINFO2]] = !{!"julia.simdloop", !"julia.ivdep"}
# CHECK: [[LOOPINFO3]] = !{[[LOOPUNROLL:![0-9]+]]}
# CHECK: [[LOOPUNROLL]] = !{!"llvm.loop.unroll.count", i64 3}
# CHECK: [[LOOPINFO4]] = !{[[LOOPUNROLL2:![0-9]+]]}
# CHECK: [[LOOPUNROLL2]] = !{!"llvm.loop.unroll.full"}
# LOWER: [[LOOPID]] = distinct !{[[LOOPID]]}
# LOWER: [[LOOPID2]] = distinct !{[[LOOPID2]]}
# LOWER: [[LOOPID3]] = distinct !{[[LOOPID3]], [[LOOPUNROLL:![0-9]+]]}
# LOWER: [[LOOPUNROLL]] = !{!"llvm.loop.unroll.count", i64 3}
# LOWER: [[LOOPID4]] = distinct !{[[LOOPID4]], [[LOOPUNROLL2:![0-9]+]]}
# LOWER: [[LOOPUNROLL2]] = !{!"llvm.loop.unroll.full"}

# Maintaining the order is important
emit(simdf, Vector{Float64})
emit(simdf2, Vector{Float64})
emit(loop_unroll, Int64)
emit(loop_unroll2, Int64, Int64)
emit(notunroll, Int64, Int64)

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using ArtifactUtils, Pkg.Artifacts
using Test
using TOML
import Pkg

function _artifact(name, loc)
    eval(Expr(:macrocall, Symbol("@artifact_str"), LineNumberNode(1, Symbol(loc)), name))
end

@testset "ArtifactUtils.jl" begin
    mktempdir() do tempdir
        artifact_file = joinpath(tempdir, "Artifacts.toml")
        add_artifact!(
            artifact_file,
            "JuliaMono",
            "https://github.com/cormullion/juliamono/releases/download/v0.007/JuliaMono.tar.gz",
            force=true,
        )
        artifacts = TOML.parsefile(artifact_file)
        @test artifacts == Dict{String,Any}(
            "JuliaMono" => Dict{String,Any}(
                "git-tree-sha1" => "65279f9c8a3dd1e2cb654fdedbe8cd58889ae1bc",
                "download" => Any[
                    Dict{String,Any}(
                        "sha256" => "f1ab65231cda7981531398644a58fd5fde8f367b681e1b8e9c35d9b2aacfcb1c",
                        "url" => "https://github.com/cormullion/juliamono/releases/download/v0.007/JuliaMono.tar.gz",
                    ),
                ],
            ),
        )
        ensure_artifact_installed("JuliaMono", artifact_file)
        @test ispath(_artifact("JuliaMono", tempdir))
        hash = ArtifactUtils.sha256sum(joinpath(_artifact("JuliaMono", tempdir), "JuliaMono-Regular.ttf"))
        @test hash == "dc0af40e8bc944a5d38c049b9a5b33e80b1e7a9621ac30fb440e5a2b6a4192d7"
    end
end


text:
module Typography

import Genie

include(joinpath("typography", "header.jl"))

end


text:
# Multinomial Logistic regression

using SGDOptim
using ArrayViews

function error_rate(W::Matrix{Float64}, X::Matrix{Float64}, y::Vector{Int})
    u = W * X
    r = Int[indmax(view(u,:,i)) for i = 1:size(X,2)]
    countnz(r .!= y) / length(y)
end

function mnlogireg_sgd(Wg::Matrix{Float64}, n::Int, σ::Float64)

    # prepare experimental data
    k, d = size(Wg)
    X = randn(d, n)
    u = Wg * X + σ * randn(k, n)
    y = Int[indmax(view(u,:,i)) for i = 1:n]

    # construct the risk model
    rmodel = riskmodel(MvLinearPred(d, k), MultiLogisticLoss())

    # initialize solution
    W0 = randn(k, d)

    # optimize
    W = sgd(rmodel, W0,
        minibatch_seq(X, y, 10),          # configure the way data are supplied
        reg = SqrL2Reg(1.0e-4),           # regularization
        lrate = t->1.0 / (100.0 + t),     # learing rate policy
        cbinterval = 100,                 # how frequently callback is invoked
        callback = simple_trace)          # the callback function

    # compare solution with initial guess
    println()
    @printf("Initial  :  error.rate = %5.2f%%\n", error_rate(W0, X, y) * 100.0)
    @printf("Solution :  error.rate = %5.2f%%\n", error_rate(W,  X, y) * 100.0)
    @printf("gTruth   :  error.rate = %5.2f%%\n", error_rate(Wg, X, y) * 100.0)
end

mnlogireg_sgd(randn(3, 5), 10000, 0.2)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using .NestedSamplers               # used in this file, because this file ist only load if it is required

include("ens_bounds.jl")
include("ens_proposals.jl")

"""
    struct EllipsoidalNestedSampling <: AbstractSamplingAlgorithm

*Experimental feature, not part of stable public API.*

Uses the julia package
[NestedSamplers.jl](https://github.com/TuringLang/NestedSamplers.jl) to use nested sampling algorithm.

Constructors:

* ```$(FUNCTIONNAME)(; fields...)```

Fields:

$(TYPEDFIELDS)

!!! note

    This functionality is only available when the
    [NestedSamplers.jl](https://github.com/TuringLang/NestedSamplers.jl) package
    is loaded (e.g. via
    `import`).
"""
@with_kw struct EllipsoidalNestedSampling{TR<:AbstractDensityTransformTarget} <: AbstractSamplingAlgorithm
    trafo::TR = PriorToUniform()

    "Number of live-points."
    num_live_points::Int = 1000

    "Volume around the live-points."
    bound::ENSBound = ENSEllipsoidBound()

    "Algorithm used to choose new live-points."
    proposal::ENSProposal = ENSAutoProposal()

    "Scale factor for the volume."
    enlarge::Float64 = 1.25

    # "Not sure about how this works yet."
    # update_interval::Float64 =

    "Number of iterations before the first bound will be fit."
    min_ncall::Int64 = 2*num_live_points

    "Efficiency before fitting the first bound."
    min_eff::Float64 = 0.1

    # "The following four are the possible convergence criteria to end the algorithm."
    dlogz::Float64 = 0.01
    max_iters = Inf
    max_ncalls = 10^7
    maxlogl = Inf
end
export EllipsoidalNestedSampling

function bat_sample_impl(rng::AbstractRNG, target::AnyDensityLike, algorithm::EllipsoidalNestedSampling)

    density_notrafo = convert(AbstractDensity, target)
    shaped_density, trafo = bat_transform(algorithm.trafo, density_notrafo)                 # BAT prior transformation
    vs = varshape(shaped_density)
    density = unshaped(shaped_density)
    dims = totalndof(vs)

    model = NestedModel(logdensityof(density), identity);                                   # identity, because ahead the BAT prior transformation is used instead
    bounding = ENSBounding(algorithm.bound)
    prop = ENSprop(algorithm.proposal)
    sampler = Nested(dims, algorithm.num_live_points;
                        bounds=bounding, proposal=prop,
                        enlarge=algorithm.enlarge, min_ncall=algorithm.min_ncall, min_eff=algorithm.min_eff
    )

    samples_w, state = sample(model, sampler;                                               # returns samples with weights as one vector and the actual state
        dlogz = algorithm.dlogz, maxiter = algorithm.max_iters,
        maxcall = algorithm.max_ncalls, maxlogl = algorithm.maxlogl, chain_type=Array
    )

    weights = samples_w[:, end]                                                             # the last elements of the vectors are the weights
    nsamples = size(samples_w,1)
    samples = [samples_w[i, 1:end-1] for i in 1:nsamples]                                   # the other ones (between 1 and end-1) are the samples
    logvals = map(logdensityof(density), samples)                                           # posterior values of the samples
    samples_trafo = vs.(BAT.DensitySampleVector(samples, logvals, weight = weights))
    samples_notrafo = inverse(trafo).(samples_trafo)                                            # Here the samples are retransformed

    logintegral = Measurements.measurement(state.logz, state.logzerr)
    ess = bat_eff_sample_size(samples_notrafo, KishESS()).result

    return (
        result = samples_notrafo, result_trafo = samples_trafo, trafo = trafo,
        logintegral = logintegral, ess = ess,
        info = state
    )
end


text:
# auto-generated tests from julia-repl docstrings
using Test, CyclotomicNumbers
function mytest(f::String,a::String,b::String)
  println(f," ",a)
  omit=a[end]==';'
  a=replace(a,"\\\\"=>"\\")
  a=repr(MIME("text/plain"),eval(Meta.parse(a)),context=:limit=>true)
  if omit a="nothing" end
  a=replace(a,r" *(\n|$)"s=>s"\1")
  a=replace(a,r"\n$"s=>"")
  b=replace(b,r" *(\n|$)"s=>s"\1")
  b=replace(b,r"\n$"s=>"")
  i=1
  while i<=lastindex(a) && i<=lastindex(b) && a[i]==b[i]
    i=nextind(a,i)
  end
  if a!=b print("exec=$(repr(a[i:end]))\nmanl=$(repr(b[i:end]))\n") end
  a==b
end
@testset verbose = true "Gapjm" begin
@testset "CyclotomicNumbers.jl" begin
@test mytest("CyclotomicNumbers.jl","E(3,2)","Root1: ζ₃²")
@test mytest("CyclotomicNumbers.jl","E(3)+E(4)","Cyc{Int64}: ζ₁₂⁴-ζ₁₂⁷-ζ₁₂¹¹")
@test mytest("CyclotomicNumbers.jl","E(12,11)-E(12,7)","Cyc{Int64}: √3")
@test mytest("CyclotomicNumbers.jl","repr(E(12,11)-E(12,7),context=(:limit=>true,:quadratic=>false))","\"-ζ₁₂⁷+ζ₁₂¹¹\"")
@test mytest("CyclotomicNumbers.jl","a=E(3)+E(3,2)","Cyc{Int64}: -1")
@test mytest("CyclotomicNumbers.jl","conductor(a)","1")
@test mytest("CyclotomicNumbers.jl","typeof(Int(a))","Int64")
@test mytest("CyclotomicNumbers.jl","inv(1+E(4))","Cyc{Float64}: 0.5-0.5ζ₄")
@test mytest("CyclotomicNumbers.jl","1//(1+E(4))","Cyc{Rational{Int64}}: (1-ζ₄)/2")
@test mytest("CyclotomicNumbers.jl","Cyc(1//2+im)","Cyc{Rational{Int64}}: (1+2ζ₄)/2")
@test mytest("CyclotomicNumbers.jl","conj(1+E(4))","Cyc{Int64}: 1-ζ₄")
@test mytest("CyclotomicNumbers.jl","real(E(3))","Cyc{Rational{Int64}}: -1/2")
@test mytest("CyclotomicNumbers.jl","Rational{Int}(real(E(3)))","-1//2")
@test mytest("CyclotomicNumbers.jl","imag(E(3))","Cyc{Rational{Int64}}: √3/2")
@test mytest("CyclotomicNumbers.jl","c=Cyc(E(9))","Cyc{Int64}: -ζ₉⁴-ζ₉⁷")
@test mytest("CyclotomicNumbers.jl","Root1(c)","Root1: ζ₉")
@test mytest("CyclotomicNumbers.jl","Root1(1+E(4))","nothing")
@test mytest("CyclotomicNumbers.jl","Root1(;r=1//4)","Root1: ζ₄")
@test mytest("CyclotomicNumbers.jl","c=E(4)*E(3)","Root1: ζ₁₂⁷")
@test mytest("CyclotomicNumbers.jl","c=Complex{Float64}(E(3))","-0.4999999999999999 + 0.8660254037844387im")
@test mytest("CyclotomicNumbers.jl","0.0+E(3)","Cyc{Float64}: 1.0ζ₃")
@test mytest("CyclotomicNumbers.jl","E(3)+1//2","Cyc{Rational{Int64}}: √-3/2")
@test mytest("CyclotomicNumbers.jl","E(3)+im","Cyc{Int64}: ζ₁₂⁴-ζ₁₂⁷-ζ₁₂¹¹")
@test mytest("CyclotomicNumbers.jl","complex(E(4))","0 + 1im")
@test mytest("CyclotomicNumbers.jl","complex(E(3))","-0.4999999999999999 + 0.8660254037844387im")
@test mytest("CyclotomicNumbers.jl","-1<Cyc(0)<1","true")
@test mytest("CyclotomicNumbers.jl","a=E(3)+E(4)","Cyc{Int64}: ζ₁₂⁴-ζ₁₂⁷-ζ₁₂¹¹")
@test mytest("CyclotomicNumbers.jl","collect(pairs(a))","3-element Vector{Pair{Int64, Int64}}:\n  4 => 1\n  7 => -1\n 11 => -1")
@test mytest("CyclotomicNumbers.jl","a[6],a[7]","(0, -1)")
@test mytest("CyclotomicNumbers.jl","coefficients(a)","12-element Vector{Int64}:\n  0\n  0\n  0\n  0\n  1\n  0\n  0\n -1\n  0\n  0\n  0\n -1")
@test mytest("CyclotomicNumbers.jl","valtype(a)","Int64")
@test mytest("CyclotomicNumbers.jl","conductor(E(9))","9")
@test mytest("CyclotomicNumbers.jl","conductor([E(3),1//2,E(4)])","12")
@test mytest("CyclotomicNumbers.jl","coefficients(Cyc(E(9)))","9-element Vector{Int64}:\n  0\n  0\n  0\n  0\n -1\n  0\n  0\n -1\n  0")
@test mytest("CyclotomicNumbers.jl","galois(1+E(4),-1)","Cyc{Int64}: 1-ζ₄")
@test mytest("CyclotomicNumbers.jl","galois(root(5),2)==-root(5)","true")
@test mytest("CyclotomicNumbers.jl","conjugates(1+root(5))","2-element Vector{Cyc{Int64}}:\n 1+√5\n 1-√5")
@test mytest("CyclotomicNumbers.jl","r=Root1(-E(9,2)-E(9,5))","Root1: ζ₉⁸")
@test mytest("CyclotomicNumbers.jl","order(r)","9")
@test mytest("CyclotomicNumbers.jl","exponent(r)","8")
@test mytest("CyclotomicNumbers.jl","Cyc(r)","Cyc{Int64}: -ζ₉²-ζ₉⁵")
@test mytest("CyclotomicNumbers.jl","Root1(-E(9,4)-E(9,5))","nothing")
@test mytest("CyclotomicNumbers.jl","Quadratic(E(3,2)-2E(3))","(1-3√-3)/2")
@test mytest("CyclotomicNumbers.jl","Quadratic(1+E(5))","nothing")
@test mytest("CyclotomicNumbers.jl","root(-1)","Cyc{Int64}: ζ₄")
@test mytest("CyclotomicNumbers.jl","root(E(4))","Root1: ζ₈")
@test mytest("CyclotomicNumbers.jl","root(27,6)","Cyc{Int64}: √3")
end
end


text:
using Plots
gr(fmt=:png);
using QuantEcon

phi_1, phi_2, phi_3, phi_4 = 0.5, -0.2, 0, 0.5
sigma = 0.1

A = [phi_1 phi_2 phi_3 phi_4
     1.0 0.0 0.0 0.0
     0.0 1.0 0.0 0.0
     0.0 0.0 1.0 0.0]

C = [sigma 0.0 0.0 0.0]'
G = [1.0 0.0 0.0 0.0]

T = 30
ar = LSS(A, C, G, mu_0=ones(4))

ymin, ymax = -0.8, 1.25

ys = []
yTs = []

for i=1:20
    x, y = simulate(ar, T+15)
    y = vec(y)
    push!(ys, y)
    push!(yTs, y[T])
end
p1 = plot(ys, linewidth=1, alpha=0.5)
scatter!(fill(T, length(yTs)), yTs, color=:black)
plot!(ylims=(ymin, ymax), xticks=[], ylabel="y_t")
vline!([T], color=:black, legend=:none)
annotate!(T+1, -0.8, "T")

p2 = histogram(yTs, bins=16, normed=true, orientation=:h, alpha=0.5)
plot!(ylims=(ymin, ymax))
plot(p1, p2, layout=2, legend=:none)

max_stars_repo_name: Shunsuke-Hori/lecture-source-jl

max_issues_repo_name: Shunsuke-Hori/lecture-source-jl

max_forks_repo_name: Shunsuke-Hori/lecture-source-jl


text:
using MPI
import MUMPS # needs to be intialized before Ipopt don't understand why

type linear_solver_MUMPS <: abstract_linear_system_solver
	_SparseMatrix::SparseMatrixCSC{Float64,Int64}
	_factor::MUMPS.Mumps{Float64} # TO DO, give type

  function linear_solver_MUMPS()
      return new();
  end
end

# intialize and finalize function???
# MPI ???
function mumps_sym(sym::Symbol)
	if sym == :symmetric
			return MUMPS.mumps_symmetric
	elseif sym == :unsymmetric
			return MUMPS.mumps_unsymmetric
	elseif sym == :definite
			return MUMPS.mumps_definite
	else
			error("this symmetry symbol is not understood")
	end
end

function initialize!(solver::linear_solver_MUMPS, A::SparseMatrixCSC{Float64,Int64}, sym::Symbol)
		if ~MPI.Initialized()
			MPI.Init()
		end

		icntl = MUMPS.default_icntl[:]; # copy
    icntl[4] = 1;
		icntl[14] = 100.0;

    solver._factor = MUMPS.Mumps{Float64}(mumps_sym(sym), icntl, MUMPS.default_cntl64	);  # Real, general unsymmetric
		solver._SparseMatrix = A;

    return solver;
end

function finalize!(solver::linear_solver_MUMPS)
    MUMPS.finalize(solver._factor);
end

function ls_factor!(solver::linear_solver_MUMPS, n::Int64, m::Int64)
		start_advanced_timer("MUMPS/associate_matrix")
    MUMPS.associate_matrix(solver._factor, solver._SparseMatrix);
		pause_advanced_timer("MUMPS/associate_matrix")

		MUMPS.factorize(solver._factor);

		return solver._factor.infog[12] == m; # inertia
end

function ls_solve(solver::linear_solver_MUMPS, my_rhs::AbstractArray)
    MUMPS.associate_rhs(solver._factor, my_rhs);
    MUMPS.solve(solver._factor);
    return MUMPS.get_solution(solver._factor);
end

max_forks_repo_name: ohinder/Everest


text:
export a_is_b_min_c_MT!

"""
in-place multithreaded a = b - c
"""
function a_is_b_min_c_MT!(a::Vector{TF},b::Vector{TF},c::Vector{TF}) where {TF<:Real}
    Threads.@threads for j=1:length(a)
      @inbounds a[j] = b[j] - c[j]
    end
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
using DataStructure101
const DS = DataStructure101

export dfs

function dfs!(g::DS.AbstractGraph, v::Int64, traversed::Vector{Int64}=Int[])
    push!(traversed, v)

    neighbors = DS.neighbor(g, v)
    for n in neighbors
        !(n in traversed) && dfs!(g, n, traversed)
    end

    return traversed
end

dfs(g::DS.AbstractGraph, v::Int64) = dfs!(g, v)


text:
using Pope, HDF5
using Pope.LJH
using Test
using ReferenceMicrocalFiles

@testset "BufferedHDF5Dataset2D" begin
    for name in ["ds","1/ds"]
        h5 = h5open(tempname(),"w")
        chunksize = 100
        # create an extendible dataset to hold reduced pulses with 6 element basis
        buffer = Pope.BufferedHDF5Dataset2D{Float32}(h5,name,6,1000)
        write(buffer, Vector{Float32}(1:6))
        write(buffer, Vector{Float32}(7:12))
        write(buffer, [Vector{Float32}(13:18), Vector{Float32}(19:24)])
        Pope.write_to_hdf5(buffer)
        @test read(buffer.ds) == reshape(Vector{Float32}(1:24),(6,4))
        close(h5)
        rm(h5.filename)
    end
end

@testset "BasisAnalyzer" begin
    FrameTime, PretrigNSamples, NRow, row = 9.6e-6, 100, 30, 1
    # rowcount of LJH2 files = framecount*NRow+row, framecount=frame1index
    r = LJH.LJHRecord{FrameTime, PretrigNSamples, NRow}(1:1000,10*NRow+row,20)
    r3 = LJH.LJH3Record{FrameTime}(1:1000,PretrigNSamples,10,20)
    projectors = rand(6,1000)
    analyzer = Pope.BasisAnalyzer(projectors, projectors')
    dp = analyzer(r)
    dp3 = analyzer(r3)
    @test Pope.nbases(analyzer) == 6
    @test Pope.nsamples(analyzer) == 1000
    @test dp.reduced == dp3.reduced
    @test dp.residual_std == dp3.residual_std
    @test dp.timestamp_usec == dp3.timestamp_usec == LJH.timestamp_usec(r)
    @test dp.first_rising_sample == dp.first_rising_sample == PretrigNSamples
    @test dp.nsamples == dp3.nsamples == length(r)
end

@testset "BasisBufferedWriter" begin
    h5 = h5open(tempname(),"w")
    g = HDF5.g_create(h5,"1")
    b = Pope.BasisBufferedWriter(g, 6, 1000, 0.001, start=true)
    dataproduct = Pope.BasisDataProduct(collect(1:6),1,2,3,4,5)
    write(b,dataproduct)
    Pope.stop(b)
    fetch(b)
    @test read(h5["1/reduced"]) == reshape(1:6,(6,1))
    @test read(h5["1/residual_std"]) == [dataproduct.residual_std]
    @test read(h5["1/frame1index"]) == [dataproduct.frame1index]
    @test read(h5["1/timestamp_usec"]) == [dataproduct.timestamp_usec]
    @test read(h5["1/first_rising_sample"]) == [dataproduct.first_rising_sample]
    @test read(h5["1/nsamples"]) == [dataproduct.nsamples]
    close(h5)
    rm(h5.filename)
end


text:
@testset "dfs" begin
    A = load_matrix_network("dfs_example")
    dfs(MatrixNetwork(A),1)

    n = 10
    O = ones(Int64,n-1)
    Z = zeros(Int64,n)
    A_ref = Tridiagonal(O,Z,O)
    B = sparse(Matrix(A_ref))
    @testset "SparseMatrixCSC" begin
        dfs(A,1)

        (d,dt,ft,pred) = dfs(B,1)
        @test d == collect(0:n-1)
        @test dt == collect(0:n-1)
        @test ft == collect(2*n-1:-1:n)
        @test pred == collect(0:n-1)
    end
    @testset "CSR" begin
        dfs(sparse_to_csr(A)..., 1)

        (d,dt,ft,pred) = dfs(sparse_to_csr(B)..., 1)
        @test d == collect(0:n-1)
        @test dt == collect(0:n-1)
        @test ft == collect(2*n-1:-1:n)
        @test pred == collect(0:n-1)
    end
    @testset "triplet" begin
        dfs(findnz(A)[1], findnz(A)[2], 1)

        (d,dt,ft,pred) = dfs(findnz(B)[1], findnz(B)[2], 1)
        @test d == collect(0:n-1)
        @test dt == collect(0:n-1)
        @test ft == collect(2*n-1:-1:n)
        @test pred == collect(0:n-1)
    end
end


text:
# GPU - apply along the features axis
function hist_kernel!(h::CuDeviceArray{T,3}, x::CuDeviceMatrix{T}, id, 𝑖, 𝑗, K) where {T<:AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y
    if i <= length(𝑖) && j <= length(𝑗)
        for k in 1:K
            @inbounds pt = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], k, 𝑗[j])
            @inbounds CUDA.atomic_add!(pointer(h, pt), x[𝑖[i],k])
        end
    end
    return
end

# for 2D input: 𝑤
function hist_kernel!(h::CuDeviceMatrix{T}, x::CuDeviceVector{T}, id, 𝑖, 𝑗) where {T<:AbstractFloat}
    i = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    j = threadIdx().y + (blockIdx().y - 1) * blockDim().y
    if i <= length(𝑖) && j <= length(𝑗)
        @inbounds pt = Base._to_linear_index(h, id[𝑖[i], 𝑗[j]], 𝑗[j])
        @inbounds CUDA.atomic_add!(pointer(h, pt), x[𝑖[i]])
    end
    return
end

# base approach - block built along the cols first, the rows (limit collisions)
function update_hist_gpu!(hδ::CuArray{T,3}, hδ²::CuArray{T,3}, h𝑤::CuMatrix{T},
    δ::CuMatrix{T}, δ²::CuMatrix{T}, 𝑤::CuVector{T},
    X_bin::CuMatrix{Int}, 𝑖::CuVector{Int}, 𝑗::CuVector{Int}, K; MAX_THREADS=1024) where {T<:AbstractFloat}

    hδ .*= 0.0
    hδ² .*= 0.0
    h𝑤 .*= 0.0

    thread_j = min(MAX_THREADS, length(𝑗))
    thread_i = min(MAX_THREADS ÷ thread_j, length(𝑖))
    threads = (thread_i, thread_j)
    blocks = ceil.(Int, (length(𝑖), length(𝑗)) ./ threads)
    @cuda blocks=blocks threads=threads hist_kernel!(hδ, δ, X_bin, 𝑖, 𝑗, K)
    @cuda blocks=blocks threads=threads hist_kernel!(hδ², δ², X_bin, 𝑖, 𝑗, K)
    @cuda blocks=blocks threads=threads hist_kernel!(h𝑤, 𝑤, X_bin, 𝑖, 𝑗)

    return
end

function find_split_gpu!(hist_δ::AbstractMatrix{T}, hist_δ²::AbstractMatrix{T}, hist_𝑤::AbstractVector{T},
    params::EvoTypes, node::TrainNode_gpu{T,S}, info::SplitInfo_gpu{T,S}, edges::Vector{T}) where {T,S}

    # initialize tracking
    ∑δL = copy(node.∑δ) .* 0
    ∑δ²L = copy(node.∑δ²) .* 0
    ∑𝑤L = node.∑𝑤 * 0
    ∑δR = copy(node.∑δ)
    ∑δ²R = copy(node.∑δ²)
    ∑𝑤R = node.∑𝑤

    # println("∑δ²L: ", ∑δ²L, " ∑δ²R:", ∑δ²R)
    # println("hist_δ²: ", hist_δ²)

    @inbounds for bin in 1:(length(hist_δ)-1)
        @views ∑δL .+= hist_δ[bin,:]
        @views ∑δ²L .+= hist_δ²[bin,:]
        ∑𝑤L += hist_𝑤[bin]
        @views ∑δR .-= hist_δ[bin,:]
        @views ∑δ²R .-= hist_δ²[bin,:]
        ∑𝑤R -= hist_𝑤[bin]

        # println("∑δ²L: ", ∑δ²L, " | ∑δ²R:", ∑δ²R, " | hist_δ²[bin,:]: ", hist_δ²[bin,:])

        gainL, gainR = get_gain(params.loss, ∑δL, ∑δ²L, ∑𝑤L, params.λ), get_gain(params.loss, ∑δR, ∑δ²R, ∑𝑤R, params.λ)
        gain = gainL + gainR

        if gain > info.gain && ∑𝑤L >= params.min_weight + 1e-12 && ∑𝑤R >= params.min_weight + 1e-12
            info.gain = gain
            info.gainL = gainL
            info.gainR = gainR
            @views info.∑δL .= ∑δL
            @views info.∑δ²L .= ∑δ²L
            info.∑𝑤L = ∑𝑤L
            @views info.∑δR .= ∑δR
            @views info.∑δ²R .= ∑δ²R
            info.∑𝑤R = ∑𝑤R
            info.cond = edges[bin]
            info.𝑖 = bin
        end # info update if gain
    end # loop on bins
end


text:
# Mandatory
export polyhedron, getinequalities, getgenerators, eliminate, detecthlinearities!, detectvlinearities!, removeredundantinequalities!, removeredundantgenerators!, isredundantinequality, isredundantgenerator, isstronglyredundantinequality, isstronglyredundantgenerator, inequalitiesarecomputed, generatorsarecomputed, loadpolyhedron!

if VERSION < v"0.5-"
  export normalize
  normalize(v,p=2) = v / norm(v,p)
end

polyhedron(repr::Representation) = polyhedron(repr, getlibraryfor(eltype(repr)))
Base.push!{N}(p::Polyhedron{N}, ine::HRepresentation{N})                             = error("push! not implemented for $(typeof(p)) for HRepresentation")
Base.push!{N}(p::Polyhedron{N}, ext::VRepresentation{N})                             = error("push! not implemented for $(typeof(p)) for VRepresentation")
hrepiscomputed(p::Polyhedron)                                                        = error("hrepiscomputed not implemented for $(typeof(p))")
gethrep(p::Polyhedron)                                                               = error("gethrep not implemented for $(typeof(p))")
vrepiscomputed(p::Polyhedron)                                                        = error("vrepiscomputed not implemented for $(typeof(p))")
getvrep(p::Polyhedron)                                                               = error("getvrep not implemented for $(typeof(p))")
implementseliminationmethod(p::Polyhedron, ::Type{Val{:FourierMotzkin}})             = false
eliminate(p::Polyhedron, delset::IntSet, ::Type{Val{:FourierMotzkin}})               = error("Fourier-Motzkin elimination not implemented for $(typeof(p))")
implementseliminationmethod(p::Polyhedron, ::Type{Val{:BlockElimination}})           = false
eliminate(p::Polyhedron, delset::IntSet, ::Type{Val{:BlockElimination}})             = error("Block elimination not implemented for $(typeof(p))")
#loadpolyhedron!(p::Polyhedron, filename::AbstractString, extension::Type{Val{:ine}}) = error("not implemented")
#loadpolyhedron!(p::Polyhedron, filename::AbstractString, extension::Type{Val{:ext}}) = error("not implemented") # FIXME ExtFileVRepresentation or just ExtFile

# These can optionally be reimplemented for speed by a library
export numberofinequalities, numberofgenerators, dim, affinehull, getredundantinequalities, getstronglyredundantinequalities, getredundantgenerators, getstronglyredundantgenerators, transforminequalities, transformgenerators, project, radialprojectoncut

loadpolyhedron!(p::Polyhedron, filename::AbstractString, extension::Symbol) = loadpolyhedron!(p, filename, Val{extension})

function loadpolyhedron!(p::Polyhedron, filename::AbstractString, extension::AbstractString)
  s = findfirst(["ext", "ine"], filename)
  if s == 0
    error("Invalid extension $extension, please give 'ext' for V-representation or 'ine' for H-representation")
  end
  loadpolyhedron!(p, filename, [:ext, :ine][s])
end

eliminate(p::Polyhedron, method::Symbol) = eliminate(p, Val{method})
eliminate(p::Polyhedron, delset::IntSet, method::Symbol) = eliminate(p, delset::IntSet, Val{method})

eliminate{N}(p::Polyhedron{N}, method::Type{Val{:ProjectGenerators}}) = eliminate(p, IntSet(N), method)

function eliminate{N}(p::Polyhedron{N}, delset::IntSet=IntSet(N))
  fm = implementseliminationmethod(p, Val{:FourierMotzkin})
  be = implementseliminationmethod(p, Val{:BlockElimination})
  if (!fm && !be) || generatorsarecomputed(p)
    method = :ProjectGenerators
  elseif fm && (!be || delset == IntSet(N))
    method = :FourierMotzkin
  else
    method = :BlockElimination
  end
  eliminate(p, delset, Val{method})
end

function eliminate{N}(p::Polyhedron{N}, delset::IntSet, ::Type{Val{:ProjectGenerators}})
  ext = getgenerators(p)
  I = eye(Int, N)
  polyhedron(I[setdiff(IntSet(1:N), collect(delset)),:] * ext, getlibrary(p))
end

function Base.convert{N, S, T}(::Type{Polyhedron{N, S}}, p::Polyhedron{N, T})
  if !inequalitiesarecomputed(p) && generatorsarecomputed(p)
    f = (i, x) -> changeeltype(typeof(x), S)(x)
    if decomposedvfast(p)
      polyhedron(PointIterator(p, f), RayIterator(p, f), getlibraryfor(p, S))
    else
      polyhedron(VRepIterator(p, f), getlibraryfor(p, S))
    end
  else
    if decomposedvfast(p)
      polyhedron(IneqIterator(p, f), EqIterator(p, f), getlibraryfor(p, S))
    else
      polyhedron(HRepIterator(p, f), getlibraryfor(p, S))
    end
  end
end

# eliminate the last dimension by default
eliminate{N,T}(p::Polyhedron{N,T})  = eliminate(p::Polyhedron, IntSet([N]))

# function transformgenerators{N}(p::Polyhedron{N}, P::AbstractMatrix)
#   # Each generator x is transformed to P * x
#   # If P is orthogonal, the new axis are the rows of P.
#   if size(P, 2) != N
#     error("The number of columns of P must match the dimension of the polyhedron")
#   end
#   ext = P * getgenerators(p)
#   polyhedron(ext, getlibraryfor(p, eltype(ext)))
# end
#
# function transforminequalities(p::Polyhedron, P::AbstractMatrix)
#   # The new axis are the column of P.
#   # Let y be the coordinates of a point x in these new axis.
#   # We have x = P * y so y = P \ x.
#   # We have
#   # b = Ax = A * P * (P \ x) = (A * P) * y
#   ine = getinequalities(p) * P
#   polyhedron(ine, getlibraryfor(p, eltype(ine)))
# end

# function (*){N,S}(A::AbstractMatrix{S}, p::Polyhedron{N})
#   if size(A, 2) != N
#     error("Incompatible dimension")
#   end
#   if generatorsarecomputed(p)
#     transformgenerators(p, A)
#   else # FIXME not wokring
#     ine = SimpleHRepresentation(getinequalities(p))
#     nnew = size(A, 1)
#     if false
#       # CDD works with delset not at the end ?
#       newA = [ine.A spzeros(S, size(ine.A, 1), nnew);
#                   A  -speye(S, nnew, nnew)]
#       delset = IntSet(nnew+(1:N))
#     else
#       newA = [spzeros(S, size(ine.A, 1), nnew) ine.A;
#                -speye(S, nnew, nnew) A]
#       delset = IntSet(1:N)
#     end
#     newb = [ine.b; spzeros(S, nnew)]
#     newlinset = ine.linset ∪ IntSet(N+(1:nnew))
#     newine = SimpleHRepresentation(newA, newb, newlinset)
#     newpoly = polyhedron(newine, getlibraryfor(p, eltype(newine)))
#     eliminate(newpoly, IntSet(nnew+(1:N)))
#   end
# end

function project{N,T}(p::Polyhedron{N,T}, P::Array)
  # Function to make x orthogonal to an orthonormal basis in Q
  # We first make the columns of P orthonormal
  if size(P, 1) != N
    error("The columns of P should have the same dimension than the polyhedron")
  end
  m = size(P, 2)
  if m > N
    error("P should have more columns than rows")
  end
  Q = Array{Float64}(P) # normalize will make it nonrational
  Proj = zeros(eltype(P), N, N)
  for i = 1:m
    Q[:,i] = normalize(Q[:,i] - Proj * Q[:,i])
    Proj += Q[:,i] * Q[:,i]'
  end
  if m == N
    basis = Q
  else
    # For the rest, we take the canonical basis and we look at
    # I - Proj * I
    I = eye(Float64, N)
    R = I - Proj
    # We take the n-m that have highest norm
    order = sortperm([dot(R[:,i], R[:,i]) for i in 1:N])
    R = I[:,order[m+1:N]]
    for i in 1:N-m
      R[:,i] = normalize(R[:,i] - Proj * R[:,i])
      Proj += R[:,i] * R[:,i]'
    end
    basis = [Q R]
  end
  eliminate(p * basis, IntSet(m+1:N))
end

# TODO rewrite, it is just cutting a cone with a half-space, nothing more
# function radialprojectoncut{N}(p::Polyhedron{N}, cut::Vector, at)
#   if myeqzero(at)
#     error("at is zero")
#   end
#   if length(cut) != N
#     error("The dimensions of the cut and of the polyhedron do not match")
#   end
#   ext = SimpleVRepresentation(getgenerators(p))
#   V = copy(ext.V)
#   R = copy(ext.R)
#   for i in 1:size(V, 1)
#     v = vec(ext.V[i,:])
#     if !myeq(dot(cut, v), at)
#       error("The nonhomogeneous part should be in the cut")
#     end
#   end
#   for i in 1:size(R, 1)
#     v = vec(ext.R[i,:])
#     if myeqzero(v)
#       # It can happen since I do not necessarily have removed redundancy
#       v = zeros(eltype(v), length(v))
#     elseif !myeq(dot(cut, v), at)
#       if myeqzero(dot(cut, v))
#         error("A ray is parallel to the cut") # FIXME is ok if some vertices are on the cut ? (i.e. at == 0, cut is not needed)
#       end
#       v = v * at / dot(cut, v)
#     end
#     R[i,:] = v
#   end
#   # no more rays nor linearity since at != 0
#   ext2 = SimpleVRepresentation([V; R])
#   polyhedron(ext2, getlibraryfor(p, eltype(ext2)))
# end

#function fulldim{N,T}(p::Polyhedron{N,T})
#  N
#end

function dim(p::Polyhedron)
  detecthlinearities!(p)
  fulldim(p) - neqs(p)
end

# function affinehull(p::Polyhedron)
#   detecthlinearities!(p)
#   typeof(p)(affinehull(getinequalities(p)))
# end

function isvredundant{N,T}(p::Polyhedron{N,T}, v::VRepElement)
  for h in hrep(p)
    if vertex in h
      return cert ? (false, Nullable{HRepElement{N,T}}(h)) : false
    end
  end
  cert ? (true, Nullable{HRepElement{N,T}}(nothing)) : true
end


text:
using DelimitedFiles

no_of_depth_increases = 0
depths = readdlm("input.txt", '\n', Int)
  for i = 2:length(depths)
    if depths[i] > depths[i-1]
      global no_of_depth_increases += 1
    end
  end
return no_of_depth_increases

print("Total Number of depth increases is: ", no_of_depth_increases)


text:
module PetriNets

using Catlab, Catlab.CategoricalAlgebra, Catlab.Graphics, Catlab.Programs
using CombinatorialSpaces.ExteriorCalculus
using Catlab.Programs.DiagrammaticPrograms: NamedGraph
using ..Simulations
using ..AlgebraicPetri
using LinearAlgebra
using Unicode

export pn2dec, expand_pres!, gen_functions

i2sub = Dict(
  '0'=>'₀', '1'=>'₁', '2'=>'₂', '3'=>'₃', '4'=>'₄', '5'=>'₅',
  '6'=>'₆', '7'=>'₇', '8'=>'₈', '9'=>'₉', '-'=>'₋'
)
i2sup = Dict(
  '0'=>'⁰', '1'=>'¹', '2'=>'²', '3'=>'³', '4'=>'⁴', '5'=>'⁵',
  '6'=>'⁶', '7'=>'⁷', '8'=>'⁸', '9'=>'⁹', '-'=>'⁻'
)

function neg_term(graph, obs, homs, prs, val::Int64)
  res = add_part!(graph, :V, vname=Symbol("$(gensym())"))
  push!(obs, prs.syntax.Form0(prs[:X]))
  add_part!(graph, :E, ename=nothing, src=val, tgt=res)
  push!(homs, prs[Symbol("neg₀")])
  return res
end

function scale_term(graph, obs, homs, prs, val::Int64, sc::Int64)
  res = add_part!(graph, :V, vname="$(gensym())")
  push!(obs, prs.syntax.Form0(prs[:X]))
  add_part!(graph, :E, ename=nothing, src=val, tgt=res)
  push!(homs, prs[Symbol("mult_$(sc)")])
  return res
end

function exp_term(graph, obs, homs, prs, val::Int64, sc::Int64)
  res = add_part!(graph, :V, vname=Symbol("$(gensym())"))
  push!(obs, prs.syntax.Form0(prs[:X]))
  add_part!(graph, :E, ename=nothing, src=val, tgt=res)
  push!(homs, prs[Symbol("^$(sc)")])
  return res
end

function rate_term(graph, obs, homs, prs, val::Int64, tr_name::Symbol, tgt::Int64)
  add_part!(graph, :E, ename=nothing, src=val, tgt=tgt)
  push!(homs, prs[Symbol("k_", tr_name)])
  return tgt
end

function rate_term(graph, obs, homs, prs, val::Int64, tr_name::Symbol)
  res = add_part!(graph, :V, vname=Symbol("$(gensym())"))
  push!(obs, prs.syntax.Form0(prs[:X]))
  rate_term(graph, obs, homs, prs, val, tr_name, res)
end

function sum_dec!(graph, obs, homs, prs, vals::Int64...)
  bin_exp!(graph, obs, homs, prs, :sum₀, vals...)
end

function prod_dec!(graph, obs, homs, prs, vals::Int64...)
  bin_exp!(graph, obs, homs, prs, :prod₀, vals...)
end

function bin_exp!(graph, obs, homs, prs, exp, vals::Int64...)
  if length(vals) == 1
    first(vals)
  else
    form0 = prs.syntax.Form0(prs[:X])
    form0x0 = prs.syntax.otimes(form0, form0)
    cur_val = first(vals)
    bin_int = add_parts!(graph, :V, length(vals)-1,
              vname=[Symbol("$(gensym())") for i in 1:(length(vals)-1)])
    bin_res = add_parts!(graph, :V, length(vals)-1,
              vname=[Symbol("$(gensym())") for i in 1:(length(vals)-1)])
    append!(obs, fill(form0x0, length(vals)-1))
    append!(obs, fill(form0, length(vals)-1))
    for i in 1:(length(vals)-1)
      add_parts!(graph, :E, 3,
                ename =[nothing, nothing, nothing],
                tgt=[vals[i+1], cur_val, bin_res[i]],
                src=[bin_int[i], bin_int[i], bin_int[i]])
      append!(homs, [prs[:proj₁_⁰⁰₀], prs[:proj₂_⁰⁰₀], prs[exp]])
      cur_val = bin_res[i]
    end
    cur_val
  end
end

function pn2dec(prs, pn::LabelledReactionNet)
  synt = prs.syntax

  # Initial state variables
  graph = NamedGraph{Symbol,Union{Symbol,Nothing}}()
  svars = add_parts!(graph, :V, ns(pn), vname=snames(pn))
  tvars = add_parts!(graph, :V, nt(pn), vname=tnames(pn))
  obs  = Vector{Any}(fill(synt.Form0(prs[:X]), ns(pn) + nt(pn))) #"($(join(vcat(snames(pn), tnames(pn)), ", ")))::Form0{X}"
  homs = Vector{Any}()
  map(1:nt(pn)) do t
    inputs = pn[incident(pn, t, :it), :is]
    res = prod_dec!(graph, obs, homs, prs, svars[inputs]...)
    rate_term(graph, obs, homs, prs, res, tname(pn, t), tvars[t])
  end

  res_s = map(1:ns(pn)) do s
    inps = tvars[pn[incident(pn, s, :os), :ot]]
    otps = tvars[pn[incident(pn, s, :is), :it]]
    otp_scaled = [neg_term(graph, obs, homs, prs, o) for o in otps]
    if isempty(vcat(inps, otp_scaled))
      res = add_part!(graph, :V, vname=Symbol("$(gensym())"))
      push!(obs, prs.syntax.Form0(prs[:X]))
      res
    else
      sum_dec!(graph, obs, homs, prs, vcat(inps, otp_scaled)...)
    end
  end

  set_subparts!(graph, res_s,
    vname=[length("$name") > 1 ?
            Symbol("∂ₜ", "$name") : Symbol(Unicode.normalize("$name"* '̇'))
            for name in snames(pn)])
  FinFunctor(obs, homs, FinCat(graph), FinCat(prs))
end

function expand_pres!(pres, pn::LabelledReactionNet)
  synt = pres.syntax
  form0 = synt.Form0(pres[:X])
  form0x0 = synt.otimes(form0, form0)
  gens = Dict(:neg₀ => synt.Hom(:neg₀, form0, form0),
              :sum₀ => synt.Hom(:sum₀, form0x0, form0),
              :prod₀ => synt.Hom(:prod₀, form0x0, form0),
              :proj₁_⁰⁰₀ => synt.Hom(:proj₁_⁰⁰₀, form0x0, form0),
              :proj₂_⁰⁰₀ => synt.Hom(:proj₂_⁰⁰₀, form0x0, form0))
  for t in tnames(pn)
    gens[t] = synt.Hom(Symbol("k_$t"), form0, form0)
  end

  for (k,v) in gens
    if !has_generator(pres, k)
      add_generator!(pres, v)
    end
  end
end

function gen_functions(pn::LabelledReactionNet, s)
  funcs = Dict{Symbol, Any}()
  funcs[:sum₀] = Dict(:operator => (s,a,b) -> (s .= a .+ b), :type => InPlaceFunc())
  funcs[:prod₀] = Dict(:operator => (s,a,b) -> (s .= a .* b), :type => InPlaceFunc())
  funcs[:neg₀] = Dict(:operator => (s,a) -> (s .= -1 * a), :type => InPlaceFunc())
  for t in 1:nt(pn)
    funcs[Symbol("k_", tname(pn, t))] = Dict(:operator => rate(pn, t) * I(nv(s)), :type => MatrixFunc())
  end
  funcs
end

end


text:
using MomentClosure, Symbolics
using Test

@variables x, y, z
@variables u, v, w
@parameters p, q, r

ps = [p]
eq = p*x
subs = Dict(x => u)
sub_eq = MomentClosure.poly_subs(eq, subs, ps)

@test isequal(sub_eq, p*u)

ps = [p, q, r]
eq = ((p*x + q*y)^2 + r*z)^2
subs = Dict(y^4 => u, x*y*z => v, z*x^2 => w)
sub_eq = MomentClosure.poly_subs(eq, subs, ps)

@test isequal(sub_eq, (p^4)*(x^4) + (q^4)*u +
                      (r^2)*(z^2) + 4p*x*(q^3)*(y^3) +
                      4q*y*(p^3)*(x^3) + 2r*(p^2)*w + 2r*z*(q^2)*(y^2) +
                      6(p^2)*(q^2)*(x^2)*(y^2) + 4p*q*r*v)

@variables x, y, z
@variables u, v, w

f = x * y
n = 2
@test isequal(0, MomentClosure.nth_differential(f, x, n))
@test isequal(u*v + v * (x - u) + u * (y - v) + (x - u)*(y - v), MomentClosure.taylor_expand(f, [x, y], [u, v], n))
f = exp(x)
n = 10
@test isequal(exp(x), MomentClosure.nth_differential(f, x, 10))
@test isequal(sum(x^i/factorial(i) for i in 0:n), MomentClosure.taylor_expand(f, [x], [0], n))

f = 1/x
for n in 1:9 # for higher values of n, numerical inaccuracies preclude comparison ...
    @test isequal(sum((-1)^i*(x-1)^i for i in 0:n), MomentClosure.taylor_expand(f, [x], [1.0], n))
end

f = log(x)
for n in 1:10
    @test isequal((-1)^(n-1)*factorial(n-1)*x^(-n), MomentClosure.nth_differential(f, x, n))
end


text:
# TODO: add more target distributions and make them iteratable

# Dimension of testing distribution
const D = 5

# Deterministic tolerance
const DETATOL = 1e-3 * D
# Random tolerance
const RNDATOL = 5e-2 * D

using Distributions: logpdf, MvNormal, InverseGamma, Normal
using DiffResults: GradientResult, value, gradient
using ForwardDiff: gradient!

ℓπ(θ) = logpdf(MvNormal(zeros(D), ones(D)), θ)

function ∂ℓπ∂θ(θ)
    res = GradientResult(θ)
    gradient!(res, ℓπ, θ)
    return (value(res), gradient(res))
end

function ℓπ_gdemo(θ)
    s = exp(θ[1])
    m = θ[2]
    logprior = logpdf(InverseGamma(2, 3), s) + log(s) + logpdf(Normal(0, sqrt(s)), m)
    loglikelihood = logpdf(Normal(m, sqrt(s)), 1.5) + logpdf(Normal(m, sqrt(s)), 2.0)
    return logprior + loglikelihood
end

function ∂ℓπ∂θ_gdemo(θ)
    res = GradientResult(θ)
    gradient!(res, ℓπ_gdemo, θ)
    return (value(res), gradient(res))
end


text:
#! /usr/bin/julia

# Rosetta Code, Count in factors

function factor_print{T<:Integer}(n::T)
    const SEP = " \u00d7 "
    -2 < n || return "-1"*SEP*factor_print(-n)
    if isprime(n) || n < 2
        return string(n)
    end
    a = T[]
    for (k, v) in factor(n)
        append!(a, k*ones(T, v))
    end
    sort!(a)
    join(a, SEP)
end

lo = -4
hi = 40
println("Factor print ", lo, " to ", hi)
for i in lo:hi
    println(@sprintf("%5d = ", i), factor_print(i))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code