text:
########### Defines differentiable step relaxations
@inline function cv_step(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  if (xU<=zero(T))
    return zero(T),zero(T)
  elseif (xL>=zero(xL))
    return one(T),zero(T)
  elseif (x>=zero(x))
    return (x/xU)^2,2*x/xU^2
  else
    return zero(T),zero(T)
  end
end
@inline function cc_step(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  if (xU<=zero(T))
     return zero(T),zero(T)
  elseif (xL>=zero(T))
     return one(T),zero(T)
  elseif (x>=zero(T))
     return one(T),zero(T)
  else
    return one(T)-(x/xL)^2,-2*x/xL^2
  end
end
@inline function cv_step_NS(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  if (xU<=zero(T))
    return zero(T),zero(T)
  elseif (xL>=zero(T))
    return one(T),zero(T)
  elseif (x>=zero(T))
    return line_seg(x,zero(T),zero(T),xU,one(T)),dline_seg(x,zero(T),zero(T),xU,one(T),1/x)
  else
    return zero(T),zero(T)
  end
end
@inline function cc_step_NS(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  if (xU<=zero(T))
     return zero(T),zero(T)
  elseif (xL>=zero(T))
     return one(T),zero(T)
  elseif (x>=zero(x))
     return one(T),zero(T)
  else
    return line_seg(x,xL,zero(T),zero(T),one(T)),dline_seg(x,xL,zero(T),zero(T),one(T),1/x)
  end
end
@inline function step(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  if (MC_param.mu >= 1)
    cc::T,dcc::T = cc_step(midcc,x.Intv.lo,x.Intv.hi)
    cv::T,dcv::T = cv_step(midcv,x.Intv.lo,x.Intv.hi)
    gcc1::T,gdcc1::T = cc_step(x.cv,x.Intv.lo,x.Intv.hi)
	  gcv1::T,gdcv1::T = cv_step(x.cv,x.Intv.lo,x.Intv.hi)
	  gcc2::T,gdcc2::T = cc_step(x.cc,x.Intv.lo,x.Intv.hi)
	  gcv2::T,gdcv2::T = cv_step(x.cc,x.Intv.lo,x.Intv.hi)
	  cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
	  cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc,dcc = cc_step_NS(midcc,x.Intv.lo,x.Intv.hi)
    cv,dcv = cv_step_NS(midcv,x.Intv.lo,x.Intv.hi)
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, step(x.Intv),x.cnst,x.IntvBox, x.xref)
end

@inline function step(x::Interval{T}) where {T<:AbstractFloat}
           isempty(x) && return emptyinterval(x)
           xmin::T = ((x.lo)<zero(T)) ? zero(T) : one(T)
           xmax::T = ((x.hi)>=zero(T)) ? one(T) : zero(T)
           return Interval{T}(xmin,xmax)
end

########### Defines sign
@inline sign(x::SMCg{N,T}) where {N,T<:AbstractFloat} = -step(x)+2*step(x)


text:
module StaffJoy

# module dependencies
importall JuMP
using Base.Test: @test
importall JSON
import Logging

# Functions available elsewhere
export schedule, unit_test, functional_test, is_unassigned_shift, anneal_availability

if "ENV" in keys(ENV) && (
    ENV["ENV"] == "stage" || ENV["ENV"] == "prod"
    )

    env = ENV["ENV"]
    Logging.configure(level=Logging.INFO)
else
    env = "dev"
    Logging.configure(level=Logging.DEBUG)
end

if "Gurobi" in keys(Pkg.installed())
    Logging.info("Gurobi detected")
    using Gurobi
    GUROBI = true
else
    GUROBI = false
    using Cbc
    Logging.info("Using CBC solver")
end

UNASSIGNED_PREFIX = "unassigned_shift_" # Just basically substring matching for when we generate unassigned shifts
UNASSIGNED_ID = 0 # user id for sending back to scheduler

# what's available elsewhere

LYFT_THROTTLE = .999
MIN_LYFT = 1.1 # If it's less than this, don't even try

# 0 generates all fixed shifts at min_shift_length
# 1 generates all fixed shifts at max_shift_length
# lower means: more shifts, easier feasibility, lower theoretical optimality
UNASSIGNED_RATIO = .4
BIFURCATE_THRESHOLD = 350 # labor sum per week

# Amount of time to search for at least one feasible solutions before adding unassigned
if "ITERATION_TIMEOUT" in keys(ENV)
    ITERATION_TIMEOUT = env["ITERATION_TIMEOUT"]
else
    ITERATION_TIMEOUT = 6 * 60 # 6 minutes
end

# Max amount of time to search for feasible solutions in a given lift
# IF AND ONLY IF at least one existing solution has been found.
# This is like "We have one solution, let's find lots for optimality"
if "SEARCH_TIMEOUT" in keys(ENV)
    SEARCH_TIMEOUT = env["SEARCH_TIMEOUT"]
else
    # Should be longer than ITERATION_TIMEOUT
    SEARCH_TIMEOUT = 20 * 60 # 20 minutes
end

# Longest time to run an individual calculation
if "CALCULATION_TIMEOUT" in keys(ENV)
    CALCULATION_TIMEOUT = env["ITERATION_TIMEOUT"]
else
    CALCULATION_TIMEOUT = 3 * 60 # 3 minutes
end

#
# Includes
#

# Module
include("StaffJoy/serial_schedulers.jl")
include("StaffJoy/unassigned.jl")
include("StaffJoy/week.jl")
include("StaffJoy/day.jl")
include("StaffJoy/helpers.jl")
include("StaffJoy/windowing.jl")
include("StaffJoy/bifurcate.jl")

# Unit Tests
include("StaffJoy/unit_tests/test.jl")
include("StaffJoy/unit_tests/helpers_test.jl")
include("StaffJoy/unit_tests/windowing_test.jl")
include("StaffJoy/unit_tests/unassigned_test.jl")
include("StaffJoy/unit_tests/week_test.jl")

# Functional Tests
TEST_THRESHOLD = 1.08  # Minimum efficiency
include("StaffJoy/functional_tests/test.jl")
include("StaffJoy/functional_tests/day_test.jl")
include("StaffJoy/functional_tests/windowing_test.jl")
include("StaffJoy/functional_tests/week_test.jl")
include("StaffJoy/functional_tests/unassigned_test.jl")
include("StaffJoy/functional_tests/bifurcate_test.jl")

# Source data for tests
include("StaffJoy/test_data/courier-bike-1.jl")
include("StaffJoy/test_data/courier-car-1.jl")
include("StaffJoy/test_data/courier-bike-monday.jl")
include("StaffJoy/test_data/courier-incorrect-availability-length.jl")
include("StaffJoy/test_data/day-prior-failure.jl")
include("StaffJoy/test_data/consecutive-week.jl")
include("StaffJoy/test_data/consecutive-week-with-preceding.jl")
include("StaffJoy/test_data/week-with-day-off.jl")
include("StaffJoy/test_data/delivery-unassigned.jl")

end#module

ext: jl

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
# This file defines the Subpopulation type. This type defines an arbitrary
#   subpopulation in the system.

export Subpopulation
"""
This type defines a subpopulation in the system. Subpopulations here are based on a node, base or compound, and a number of conditions based on time, the personnel member's history, and their attributes.

The type contains the following fields:
* `name::String`: the name of the subpopulation.
* `sourceNode`: the name of the node (base or compound) from which to generate the subpopulation.
* `timeConds::Vector{MPcondition}`: the list of conditions relating to time. These are conditions on the age, the tenure, and the time since the last transition of the personnel member.
* `historyConds::Vector{MPcondition}`: the list of conditions relating to the history of personnel members. These are conditions on the node they entered the system in, on nodes they were in during their career, and transitions they did.
* `attributeConds::Vector{MPcondition}`: the list of conditions relating to the attributes of the personnel members.

Constructor:
```
Subpopulation( name::String )
```
This constructor creates a `Subpopulation` object with name `name`, the entire population as its source node, and no extra conditions.
"""
mutable struct Subpopulation

    name::String
    sourceNode::String
    timeConds::Vector{MPcondition}
    historyConds::Vector{MPcondition}
    attributeConds::Vector{MPcondition}

    function Subpopulation( name::String )::Subpopulation

        newSubpop = new()
        newSubpop.name = name
        newSubpop.sourceNode = "active"
        newSubpop.timeConds = Vector{MPcondition}()
        newSubpop.historyConds = Vector{MPcondition}()
        newSubpop.attributeConds = Vector{MPcondition}()
        return newSubpop

    end  # Subpopulation( name )

end  # mutable struct Subpopulation


text:
using StaticArrays

abstract type AbstrLTriMatrix end

"""
    LTriMatrix(dim, x)

A lower-triangular matrix of dimension `dim X dim`, whose entries have
the same type as `x`. The parameter `x` is used only for its type. Its
value is not used. A matrix is constructed with a call such as

    a = LTriMatrix(4, BigInt(1))

which generates a 4X4 lower triangular matrix whose entries are of
type BigInt. Then, a[4,3] returns the entry at row 4, column 3, and
a[4,3] throws a BoundsError. The entries of the matrix are initialized
with zeroes.

For values such that `isbits(x)` is true, use BitsLTriMatrix, which is
faster.

Entries are stored in column-major order.

    x        offset of col 1 = 0
    x x      offset of col 2 = 4
    x x x    offset of col 3 = 7
    x x x x  offset of col 4 = 9

    A[i, j] = data[i-j+1 + offset[j]]
    A[4, 4] = data[1 + 9] = data[10]
"""
mutable struct LTriMatrix{N,M,T} <: AbstrLTriMatrix
    data :: Vector{T}        # elements in matrix
    offset :: SVector{N,Int} # num elements before column j

    function LTriMatrix(dim, x)
        t = typeof(x) # element type
        n = Int(dim)
        m = (n*(n+1)) ÷ 2
        mat = new{n, m, t}()

        mat.data = zeros(t, m)
        offset = MVector{n,Int}(undef)

        offset[1] = 0
        for i in 2:n
            offset[i] = offset[i-1] + n - i + 2
        end
        mat.offset = SVector(offset)
        mat
    end
end

"""
    BitsLTriMatrix(dim, x)

Like LTriMatrix, but stores data in a static array for speed. The
elements of the matrix must be a bits type.  The entries of the matrix
are not initialized.

A lower-triangular matrix of dimension `dim X dim`, whose entries have
the same type as `x`. The parameter `x` is used only for its type. Its
value is not used. The type of x must be a bits type, so
that `isbits(x)` returns `true`.

A matrix is constructed with a call such as

    a = LTriMatrix(4, 1.0)

which generates a 4X4 lower triangular matrix whose entries are of
type Float64. Then, a[4,3] returns the entry at row 4, column 3, and
a[4,3] throws a BoundsError.
"""
mutable struct BitsLTriMatrix{N,M,T} <: AbstrLTriMatrix
    data :: MVector{M,T}     # elements in matrix
    offset :: SVector{N,Int} # num elements before column j

    function BitsLTriMatrix(dim, x)
        t = typeof(x) # element type
        !isbitstype(t) && throw(DomainError(t,
            "BitsLTriMatrix can only hold bitstype values."))
        n = Int(dim)
        m = (n*(n+1)) ÷ 2
        mat = new{n, m, t}()

        mat.data = MVector{m, t}(undef)
        offset = MVector{n,Int}(undef)

        offset[1] = 0
        for i in 2:n
            offset[i] = offset[i-1] + n - i + 2
        end
        mat.offset = SVector(offset)
        mat
    end
end

@inline function Base.size(A::AbstrLTriMatrix) where {N,M,T}
    dim = length(A.offset)
    (dim, dim)
end

@inline Base.length(A::AbstrLTriMatrix) where {N,M,T} = length(A.data)

@inline Base.IndexStyle(::AbstrLTriMatrix) where {N,M,T} = IndexCartesian()

@inline function Base.checkbounds(A::AbstrLTriMatrix, i::Int,
                          j::Int) where {N,M,T}
    i < j && throw(BoundsError(A, (i,j)))
end

@inline function Base.getindex(A::AbstrLTriMatrix, i::Int,
                               j::Int) where {N,M,T}
    @boundscheck checkbounds(A, i, j)
    A.data[i-j+1 + A.offset[j]]
end

# Fails if value cannot be converted to type T.
@inline function Base.setindex!(A::AbstrLTriMatrix, value, i::Int,
                                j::Int) where {N,M,T}
    @boundscheck checkbounds(A, i, j)
    A.data[i-j+1 + A.offset[j]] = value
end

"""
Write a representation of an LTriMatrix to io.
"""
function Base.show(io::IO, A::AbstrLTriMatrix) where {N,M,T}
    for i in 1:length(A.offset)
        for j in 1:i
            print(io, j>1 ? " " : "", A[i,j])
        end
        println(io)
    end
end


text:
using BinDeps

include("config_init.jl")

if Sys.which("cif2cell") === nothing
    relpath(args...) = joinpath(@__DIR__, args...)

    pythonpath = relpath("python2")
    prevpythfound = false
    package_basepath = joinpath(@__DIR__, "../../")
    for d in readdir(package_basepath)
        test_pythonpath = joinpath(package_basepath, d, "deps/python2/")
        if ispath(test_pythonpath) && !ispath(pythonpath) && !isempty(readdir(test_pythonpath))
            @info "Previous python2 found at $test_pythonpath."
            cp(test_pythonpath, pythonpath)
            prevpythonfound = true
            break
        end
    end
    if !ispath(pythonpath)
        mkdir(pythonpath)
        dlpath = relpath("downloads")
        if !ispath(dlpath)
            mkdir(dlpath)
        end
        #From Conda installation
        url = "https://repo.continuum.io/miniconda/Miniconda2-latest-"
        if Sys.isapple()
            url *= "MacOSX"
        elseif Sys.islinux()
            url *= "Linux"
        elseif Sys.iswindows()
            url = "https://repo.continuum.io/miniconda/Miniconda2-4.5.4-"
            url *= "Windows"
        else
            error("Unsuported OS.")
        end
        url *= Sys.WORD_SIZE == 64 ? "-x86_64" : "-x86"
        url *= Sys.iswindows() ? ".exe" : ".sh"

        if Sys.isunix()
            installer = joinpath(dlpath, "installer.sh")
        end
        if Sys.iswindows()
            installer = joinpath(dlpath, "installer.exe")
        end
        download(url, installer)
        if Sys.isunix()
            chmod(installer, 33261)  # 33261 corresponds to 755 mode of the 'chmod' program
            run(`$installer -b -f -p $pythonpath`)
        end
        if Sys.iswindows()
            run(Cmd(`$installer /S /AddToPath=0 /RegisterPython=0 /D=$pythonpath`;
                    windows_verbatim = true))
        end

        tarpath = relpath("external", "cif2cell.tar.gz")
        run(unpack_cmd(tarpath, @__DIR__, ".gz", ".tar"))
        cif2celldir = relpath("cif2cell-1.2.10")
        cd(cif2celldir)
        pyex = Sys.iswindows() ? joinpath(pythonpath, "python") :
               joinpath(pythonpath, "bin", "python2")
        run(`$pyex setup.py install --prefix=$pythonpath`)
        cd("..")
        #stupid urlopen
        starfile = Sys.iswindows() ?
                   joinpath(pythonpath, "lib", "site-packages", "StarFile.py") :
                   joinpath(pythonpath, "lib", "python2.7", "site-packages", "StarFile.py")
        starsource = read(starfile, String)
        starsource = replace(starsource,
                             "filestream = urlopen(filename)" => "filestream = urlopen('file:' + filename)")
        write(starfile, starsource)
        #
        rm(relpath("cif2cell-1.2.10"); recursive = true)
        rm(relpath("downloads"); recursive = true)
    end
end

include("asset_init.jl")


text:
@assert isequal(DataVector[1, 2, NA],
	            DataArray(PooledDataVector[1, 2, NA]))

# Test vector() and matrix() conversion tools
dv = dataones(5)
@assert isa(vector(dv), Vector{Float64})
@assert isa(convert(Vector{Float64}, dv), Vector{Float64})
dv[1] = NA
# Should raise errors:
# vector(dv)
# convert(Vector{Float64}, dv)
@assert isa(vector(dv, Any), Vector{Any})
@assert isnan(vector(dv, Float64, NaN)[1])

dm = dataones(3, 3)
@assert isa(matrix(dm), Matrix{Float64})
@assert isa(convert(Matrix{Float64}, dm), Matrix{Float64})
dm[1, 1] = NA
# Should raise errors:
# matrix(dm)
# convert(Matrix{Float64}, dm)
@assert isa(matrix(dm, Any), Matrix{Any})
@assert isnan(matrix(dm, Float64, NaN)[1, 1])

df = DataFrame()
df["A"] = 1:5
df["B"] = ["A", "B", "C", "D", "E"]
@assert isa(matrix(df), Matrix{Any})
@assert isa(matrix(df, Any), Matrix{Any})

df = DataFrame()
df["A"] = 1:5
df["B"] = 1.0:5.0
@assert isa(matrix(df), Matrix{Real})
@assert isa(matrix(df, Any), Matrix{Any})
@assert isa(matrix(df, Float64), Matrix{Float64})

df = DataFrame()
df["A"] = 1.0:5.0
df["B"] = 1.0:5.0
@assert isa(matrix(df), Matrix{Float64})
@assert isa(matrix(df, Any), Matrix{Any})
@assert isa(matrix(df, Int64), Matrix{Int64})


text:
module RingArrays

using VirtualArrays

import Base: size, getindex, checkbounds, display, RefValue, showerror
export RingArray, size, checkbounds, display
export getindex, load_block, RingArrayBoundsError
export showerror, OverwriteError, RingArrayFullError, RingArrayBoundsError

type RingArray{T, N} <: AbstractArray{T, N}
    next_write::Int
    max_blocks::Int
    blocks::Array{AbstractArray{T, N}, 1}
    num_users::Array{Int, 1}
    block_size::Tuple
    block_length::Int
    range::UnitRange{Int}
    data_length::Int

    function RingArray(;max_blocks::Int=10, block_size::Tuple=(10,), data_length::Int=100)
        if length(block_size) != N
            throw(DimensionMismatch("block size $(block_size) does not match dimension of RingArray $(N)"))
        end
        return new(1, max_blocks,
            Array{AbstractArray{T, N}, 1}(max_blocks),
            zeros(Int, max_blocks), block_size,
            block_size[1], 1:0, data_length)
    end
end

type OverwriteError <: Exception
    ring::RingArray
end

type RingArrayBoundsError <: Exception
    ring::RingArray
    i
end

type RingArrayFullError <: Exception
    ring::RingArray
end

function showerror(io::IO, err::OverwriteError)
    next_write = err.ring.next_write
    num_users = err.ring.num_users[err.ring.next_write]
    print(io, "OverwriteError: Cannot overwrite block $(next_write) since it has $(num_users) views")
end

function showerror(io::IO, err::RingArrayBoundsError)
    range = tuple(err.ring.range, err.ring.block_size[2:end]...)
    print(io, "RingArrayBoundsError: Cannot index $(err.i), outside of range $range")
end

function showerror(io::IO, err::RingArrayFullError)
    print(io, "RingArrayFullError: Cannot load another block, max data length is $(err.ring.data_length)")
end

function display(ring::RingArray)
    display(ring.blocks)
end

function size{T, N}(ring::RingArray{T, N})
    if N <= 1
        return tuple(ring.data_length,)
    else
        return tuple(ring.data_length, ring.block_size[2:end]...)
    end
end

checkbounds(ring::RingArray) = true # because of warnings

function checkbounds{T, N}(ring::RingArray{T, N}, indexes::UnitRange{Int}...)
    num_ranges = length(indexes)
    if num_ranges >= N
        if ring.range.start > indexes[1].start || ring.range.stop < indexes[1].stop
            throw(RingArrayBoundsError(ring, indexes))
        end
        for i in 2:N
            if 1 > indexes[i].start || ring.block_size[i] < indexes[i].stop
                throw(RingArrayBoundsError(ring, indexes))
            end
        end
    elseif num_ranges == 1
        index = indexes[1]
        if ring.range.start > fix_zero_index(index.start, ring.data_length)  ||
            ring.range.stop < fix_zero_index(index.stop, ring.data_length)
            throw(RingArrayBoundsError(ring, indexes))
        end
        if divide(index.start, ring.data_length) != divide(index.stop, ring.data_length)
            throw(RingArrayBoundsError(ring, indexes))
        end
        if get_last_index(ring) < index.stop
            throw(RingArrayBoundsError(ring, indexes))
        end
    else
        if ring.range.start > indexes[1].start || ring.range.stop < indexes[1].stop
            throw(RingArrayBoundsError(ring, indexes))
        end
        for i in 2:num_ranges - 1
            if 1 > indexes[i].start || ring.block_size[i] < indexes[i].stop
                throw(RingArrayBoundsError(ring, indexes))
            end
        end
        largest_index_allowed = 1
        for i in num_ranges:N
            largest_index_allowed *= size(ring)[i]
        end
        if 1 > indexes[num_ranges].start || largest_index_allowed < indexes[num_ranges].stop
            throw(RingArrayBoundsError(ring, indexes))
        end
    end
    return true
end

function checkbounds{T, N}(ring::RingArray{T, N}, indexes::Int...)

    if ring.range.start > indexes[1] || ring.range.stop < indexes[1]
        throw(RingArrayBoundsError(ring, indexes))
    else
        for i in 2:N
            if 1 > indexes[i] || ring.block_size[i] < indexes[i]
                throw(RingArrayBoundsError(ring, indexes))
            end
        end
    end
    return true
end

getindex(ring::RingArray) = nothing # Warnings told me to make this

function getindex(ring::RingArray, i::Int...)
    i = expand_index(ring, i...)
    checkbounds(ring, i...)
    result = 0

    block_index = fix_zero_index(divide(i[1], ring.block_length), ring.max_blocks)
    index = collect([fix_zero_index(i[1], ring.block_length), i[2:end]...])

    return ring.blocks[block_index][index...]
end

function getindex(ring::RingArray, i::UnitRange...)
    add_users(ring, i...)
    return get_view(ring, i...)
end

function load_block{T, N}(ring::RingArray{T ,N}, block::AbstractArray{T, N})
    can_load_block!(ring, block)

    ring.blocks[ring.next_write] = block
    ring.next_write = fix_zero_index(ring.next_write + 1, ring.max_blocks)

    if ring.range.stop < ring.max_blocks * ring.block_length
        ring.range = ring.range.start:ring.range.stop + ring.block_length
    else
        ring.range = ring.range.start + ring.block_length:ring.range.stop + ring.block_length
    end
    nothing
end

# Util functions

function can_load_block!(ring::RingArray, block::AbstractArray)
    if ring.num_users[ring.next_write] > 0
        gc()
        if ring.num_users[ring.next_write] > 0
            throw(OverwriteError(ring))
        end
    end

    if ring.range.stop >= ring.data_length
        throw(RingArrayFullError(ring))
    end

    check_dimensions(ring, block)
end

function expand_index{T, N}(ring::RingArray{T, N}, i::Int...)
    result = collect(i)
    len_needed = N
    len_have = length(i)

    for at in len_have:len_needed - 1
        last_value = result[end]
        if at == 1
            result[end] = fix_zero_index(last_value, ring.data_length)
            push!(result, divide(last_value, ring.data_length))
        else
            result[end] = fix_zero_index(last_value, ring.block_size[at])
            push!(result, divide(last_value, ring.block_size[at]))
        end
    end

    return result
end

function fix_zero_index(value::Int, s::Int)
    value = value % s
    if value == 0
        return s
    else
        return value
    end
end

function divide(value::Int, s::Int)
    if value % s == 0
        return div(value, s)
    else
        return div(value, s) + 1
    end
end

function add_users(ring::RingArray, i::UnitRange...)

    first_block = divide(i[1].start, ring.block_length)
    last_block = divide(i[1].stop, ring.block_length)

    for i in first_block:last_block
        block_index = fix_zero_index(i, ring.max_blocks)
        ring.num_users[block_index] += 1
    end

end

function remove_users(ring::RingArray, i::UnitRange...)

    first_block = divide(i[1].start, ring.block_length)
    last_block = divide(i[1].stop, ring.block_length)

    for i in first_block:last_block
        block_index = fix_zero_index(i, ring.max_blocks)
        ring.num_users[block_index] -= 1
    end

end

function get_view(ring::RingArray, i::UnitRange...)

    view = virtual_vcat(sub(ring, i...))

    first_block = divide(i[1].start, ring.block_length)
    last_block = divide(i[1].stop, ring.block_length)

    function when_done(view)
        remove_users(ring, i...)
    end

    finalizer(view, when_done)

    return view
end

function check_dimensions(ring::RingArray, block::AbstractArray)
    if ring.block_size != size(block)
        throw(DimensionMismatch("block size $(size(block)) does not match what RingArray expects $(ring.block_size)"))
    end
end

function get_last_index(ring::RingArray)
    last_index = 1
    for i in size(ring)
        last_index *= i
    end
    return last_index
end

end


text:
"""
    normalize(x)

Transforms a univariate array into a 0 mean and unit standard deviation array.

# Examples

```jldoctest
julia> GrammianAngularFields.normalize(collect(1:5))
5-element Array{Float64,1}:
 -1.2649110640673518
 -0.6324555320336759
  0.0
  0.6324555320336759
  1.2649110640673518
```
"""
function normalize(x::AbstractVector)
    m = StatsBase.mean(x)
    (x .- m) ./ StatsBase.stdm(x, m)
end

"""
    standardize(x)

Transforms a univariate array into a 0, 1 range.
"""
function standardize(x::AbstractVector)
    (x .- minimum(x)) ./ (maximum(x)-minimum(x))
end

"""
    standardize(x, a, b)

Transforms a univariate array into a a, b range.
"""
function standardize(x::AbstractVector, a::Number, b::Number)
    standardize(x) .* (b-a) .+ a
end


text:
"""
`module JAC.InteractionStrength`
    ... a submodel of JAC that contains all methods for evaluating the interaction strength (reduced matrix elements)
        for various atomic interactions.
"""
module InteractionStrength

    using  GSL, JAC, ..AngularMomentum, ..Basics, ..Bsplines, ..Defaults, ..ManyElectron, ..Nuclear, ..Radial, ..RadialIntegrals

    """
    `struct  InteractionStrength.XLCoefficient`  ... defines a type for coefficients of the two-electron (Breit) interaction

        + mu        ::Int64      ... Kind of integral, presently fixed to mu = 5 to keep similarity with RATIP during development.
        + nu        ::Int64      ... Rank of the integral.
        + a         ::Orbital    ... Orbitals a, b, c, d.
        + b         ::Orbital
        + c         ::Orbital
        + d         ::Orbital
        + coeff     ::Float64    ... corresponding coefficient.
    """
    struct XLCoefficient
        mu          ::Int64
        nu          ::Int64
        a           ::Orbital
        b           ::Orbital
        c           ::Orbital
        d           ::Orbital
        coeff       ::Float64
    end

    """
    `InteractionStrength.dipole(a::Orbital, b::Orbital, grid::Radial.Grid)`
        ... computes the  <a|| d ||b>  reduced matrix element of the dipole operator for orbital functions a, b.
            A value::Float64 is returned.
    """
    function dipole(a::Orbital, b::Orbital, grid::Radial.Grid)
        wa = AngularMomentum.CL_reduced_me(a.subshell, 1, b.subshell) * RadialIntegrals.rkDiagonal(1, a, b, grid)
        return( wa )
    end

    """
    `InteractionStrength.bosonShift(a::Orbital, b::Orbital, potential::Array{Float64,1}, grid::Radial.Grid)`
        ... computes the  <a|| h^(boson-field) ||b>  reduced matrix element of the boson-field shift Hamiltonian for orbital
            functions a, b. This boson-field shift Hamiltonian just refers to the effective potential of the given
            isotope due to the (assumed) boson mass. A value::Float64 is returned.
    """
    function bosonShift(a::Orbital, b::Orbital, potential::Array{Float64,1}, grid::Radial.Grid)
        wa = RadialIntegrals.isotope_boson(a, b, potential, grid)
        ### wa = RadialIntegrals.isotope_boson(a, b, potential, grid)
        ##x println("**  <$(a.subshell) || h^(boson-field shift) || $(b.subshell)>  = $wa" )
        return( wa )
    end

    """
    `InteractionStrength.fieldShift(a::Orbital, b::Orbital, deltaPotential::Array{Float64,1}, grid::Radial.Grid)`
        ... computes the  <a|| h^(field-shift) ||b>  reduced matrix element of the field-shift Hamiltonian for orbital
            functions a, b. This field-shift Hamiltonian just refers to the difference of the nuclear potential
            deltaPotential of two isotopes, and which is already divided by the difference of the mean-square radii.
            A value::Float64 is returned.
    """
    function fieldShift(a::Orbital, b::Orbital, deltaPotential::Array{Float64,1}, grid::Radial.Grid)
        wa = RadialIntegrals.isotope_field(a, b, deltaPotential, grid)
        ### wa = RadialIntegrals.isotope_field(a, b, deltaPotential, grid)
        ##x println("**  <$(a.subshell) || h^(field-shift) || $(b.subshell)>  = $wa" )
        return( wa )
    end

    """
    `InteractionStrength.hamiltonian_nms(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid)`
        ... computes the  <a|| h_nms ||b>  reduced matrix element of the normal-mass-shift Hamiltonian for orbital
            functions a, b. A value::Float64 is returned.  For details, see Naze et al., CPC 184 (2013) 2187, Eq. (37).
    """
    function hamiltonian_nms(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid)
        if  a.subshell.kappa != b.subshell.kappa   return( 0. )   end
        wa = RadialIntegrals.isotope_nms(a, b, nm.Z, grid) / 2
        ##x println("**  <$(a.subshell) || h^nms || $(b.subshell)>  = $wa" )
        return( wa )
    end

    """
    `InteractionStrength.hfs_t1(a::Orbital, b::Orbital, grid::Radial.Grid)`
        ... computes the <a|| t^(1) ||b> reduced matrix element for the HFS coupling to the magnetic-dipole moment
            of the nucleus for orbital functions a, b. A value::Float64 is returned.
    """
    function hfs_t1(a::Orbital, b::Orbital, grid::Radial.Grid)
        # Use Andersson, Jönson (2008), CPC, Eq. (49) ... test for the proper definition of the C^L tensors.
        minusa = Subshell(1, -a.subshell.kappa)
        wb =   - (a.subshell.kappa + b.subshell.kappa) * AngularMomentum.CL_reduced_me_rb(minusa, 1, b.subshell)
        wc =   RadialIntegrals.rkNonDiagonal(-2, a, b, grid)
        wa =   wb * wc
        #
        ##x println("**  <$(a.subshell) || t1 || $(b.subshell)>  = $wa   = $wb * $wc" )
        return( wa )
    end

    """
    `InteractionStrength.hfs_t2(a::Orbital, b::Orbital, grid::Radial.Grid)`
        ... computes the <a|| t^(2) ||b> reduced matrix element for the HFS coupling to the electric-quadrupole moment of
            the nucleus for orbital functions a, b. A value::Float64 is returned.
     """
    function hfs_t2(a::Orbital, b::Orbital, grid::Radial.Grid)
        # Use Andersson, Jönson (2008), CPC, Eq. (49) ... test for the proper definition of the C^L tensors.
        wb = - AngularMomentum.CL_reduced_me_rb(a.subshell, 2, b.subshell)
        ## wc =   RadialIntegrals.rkDiagonal(-3, a, b, grid)
        wc =   RadialIntegrals.rkDiagonal(-3, a, b, grid)
        wa =   wb * wc
        #
        ##x println("**  <$(a.subshell) || t2 || $(b.subshell)>  = $wa   = $wb * $wc" )
        return( wa )
    end

    """
    `InteractionStrength.MbaAbsorptionCheng(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)`
        ... to compute the (single-electron reduced matrix element) interaction strength <b || O^(Mp, absorption) || a>
            for the interaction with the Mp multipole component of the radiation field and the transition frequency omega, and
            within the given gauge. A value::Float64 is returned.
    """
    function MbaAbsorptionCheng(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)
        wa = MbaEmissionCheng(mp, gauge, omega, b, a, grid)
        wa = conj(wa)
        return( wa )
    end

    """
    `InteractionStrength.MbaEmissionCheng(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)`
        ... to compute the (single-electron reduced matrix element) interaction strength <b || O^(Mp,emission) || a>
            for the interaction with the Mp multipole component of the radiation field and the transition frequency omega, and
            within the given gauge. A value::Float64 is returned.
            This procedure has been first worked out with Andrey; in this case, however, the phases are not under good control,
            and this gives rise to wrong amplitudes and rates. The procedure is currently not in use.
    """
    function MbaEmissionChengOld(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)
        kapa = a.subshell.kappa;   kapb = b.subshell.kappa;    q = omega / Defaults.getDefaults("speed of light: c")
        #
        if       gauge == Basics.Magnetic
            ChengI = AngularMomentum.ChengI(-kapa, kapb, AngularJ64(mp.L))
            wa     = -1.0im / sqrt(mp.L*(mp.L+1)) * ChengI * (kapb + kapa) * RadialIntegrals.GrantILplus(mp.L, q, a, b, grid::Radial.Grid)
        #
        elseif   gauge == Basics.Babushkin
        ##x elseif   gauge in [Basics.Babushkin, Basics.Coulomb ]
            ChengI = AngularMomentum.ChengI(kapa, kapb, AngularJ64(mp.L))
            wr     = (mp.L+1) * RadialIntegrals.GrantJL(mp.L, q, a, b, grid::Radial.Grid)
            wr     = wr +  (kapb-kapa-mp.L-1) * RadialIntegrals.GrantIL0(mp.L+1, q, a, b, grid::Radial.Grid)
            wr     = wr +  (kapb-kapa+mp.L+1) * RadialIntegrals.GrantIL0(mp.L+1, q, b, a, grid::Radial.Grid)
            wa     = 1.0im / sqrt(mp.L*(mp.L+1)) * ChengI * wr
        #
        elseif   gauge == Basics.Coulomb
            ChengI = AngularMomentum.ChengI(kapa, kapb, AngularJ64(mp.L))
            wr     = (mp.L+1) * (kapb-kapa+mp.L) / (2mp.L+1) * RadialIntegrals.GrantIL0(mp.L-1, q, a, b, grid::Radial.Grid)
            wr     = wr +  ((mp.L+1) * (kapb-kapa+mp.L) - (2mp.L+1) * (kapb-kapa)) / (2mp.L+1) *
                                                               RadialIntegrals.GrantIL0(mp.L+1, q, a, b, grid::Radial.Grid)
            wr     = wr + (mp.L+1) * (kapb-kapa-mp.L) / (2mp.L+1) * RadialIntegrals.GrantIL0(mp.L-1, q, b, a, grid::Radial.Grid)
            wr     = wr +  ((mp.L+1) * (kapb-kapa-mp.L) - (2mp.L+1) * (kapb-kapa)) / (2mp.L+1) *
                                                               RadialIntegrals.GrantIL0(mp.L+1, q, b, a, grid::Radial.Grid)
            wa     = 1.0im / sqrt(mp.L*(mp.L+1)) * ChengI * wr
        else     error("stop a")
        end
        wa = conj(wa)

        return( wa )
    end

    """
    `InteractionStrength.MbaEmissionJohnsonx(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)`
        ... to compute the (single-electron reduced matrix element) interaction strength <b || O^(Mp,emission) || a>
            for the interaction with the Mp multipole component of the radiation field and the transition frequency omega, and
            within the given gauge. A value::Float64 is returned.
            This procedure has been adapted from Jiri's work but modified for the Coulomb gauge which was apparently wrong following
            some former implementation with RATIP.
    """
    function MbaEmissionJohnsonx(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)
        kapa = a.subshell.kappa;   kapb = b.subshell.kappa;    q = omega / Defaults.getDefaults("speed of light: c")
        #
        if       gauge == Basics.Magnetic
            JohnsonI = AngularMomentum.JohnsonI(-kapb, kapa, AngularJ64(mp.L))
            wa     = JohnsonI * (kapa + kapb)/(mp.L+1) * RadialIntegrals.GrantILplus(mp.L, q, a, b, grid::Radial.Grid)
        #
        elseif   gauge == Basics.Babushkin
            JohnsonI = AngularMomentum.JohnsonI(kapb, kapa, AngularJ64(mp.L))
            wr     = - RadialIntegrals.GrantJL(mp.L, q, a, b, grid::Radial.Grid)
            wr     = wr +  (kapb-kapa)/(mp.L+1) * RadialIntegrals.GrantILplus(mp.L+1, q, a, b, grid::Radial.Grid)
            wr     = wr +  RadialIntegrals.GrantILminus(mp.L+1, q, a, b, grid::Radial.Grid)
            wa     = JohnsonI * wr
        #
        elseif   gauge == Basics.Coulomb
            # This reduced matrix element still follows RATIP apart from JohnsonI
            ChengI = AngularMomentum.JohnsonI(kapb, kapa, AngularJ64(mp.L))
            wr     = (mp.L+1) * (kapb-kapa+mp.L) / (2mp.L+1) * RadialIntegrals.GrantIL0(mp.L-1, q, a, b, grid::Radial.Grid)
            wr     = wr +  ((mp.L+1) * (kapb-kapa+mp.L) - (2mp.L+1) * (kapb-kapa)) / (2mp.L+1) *
                                                               RadialIntegrals.GrantIL0(mp.L+1, q, a, b, grid::Radial.Grid)
            wr     = wr + (mp.L+1) * (kapb-kapa-mp.L) / (2mp.L+1) * RadialIntegrals.GrantIL0(mp.L-1, q, b, a, grid::Radial.Grid)
            wr     = wr +  ((mp.L+1) * (kapb-kapa-mp.L) - (2mp.L+1) * (kapb-kapa)) / (2mp.L+1) *
                                                               RadialIntegrals.GrantIL0(mp.L+1, q, b, a, grid::Radial.Grid)
            wa     = 1.0 / (mp.L+1) * ChengI * wr
        else     error("stop a")
        end
        ##x Multiply with the multipolarity factors to keep different multipoles on the same footings
        ##x wa = wa * sqrt( (2mp.L+1)*(mp.L+1)/mp.L )
        ##x wa = conj(wa)

        return( wa )
    end

    """
    `InteractionStrength.MabEmissionJohnsony(mp::EmMultipole, gauge::EmGauge, omega::Float64, a::Orbital, b::Orbital, grid::Radial.Grid)`
        ... to compute the (single-electron reduced matrix element) interaction strength <a || O^(Mp,emission) || b>
            for the interaction with the Mp multipole component of the radiation field and the transition frequency omega, and
            within the given gauge. A value::Float64 is returned. This procedure has been re-worked due to the book by Johnson (2007).
    """
    function MabEmissionJohnsony(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)
        kapa = a.subshell.kappa;   kapb = b.subshell.kappa;    q = omega / Defaults.getDefaults("speed of light: c")
        #
        if       gauge == Basics.Magnetic
            JohnsonI = AngularMomentum.JohnsonI(-kapa, kapb, AngularJ64(mp.L))
            wa     = JohnsonI * (kapa + kapb)/(mp.L+1) * RadialIntegrals.GrantILplus(mp.L, q, a, b, grid::Radial.Grid)
            #
        elseif   gauge == Basics.Babushkin
            JohnsonI = AngularMomentum.JohnsonI(kapa, kapb, AngularJ64(mp.L))
            wr     = RadialIntegrals.GrantJL(mp.L, q, a, b, grid::Radial.Grid)
            wr     = wr +  (kapa-kapb)/(mp.L+1) * RadialIntegrals.GrantILplus(mp.L+1, q, a, b, grid::Radial.Grid)
            wr     = wr +  RadialIntegrals.GrantILminus(mp.L+1, q, a, b, grid::Radial.Grid)
            wa     = JohnsonI * wr
            #
        elseif   gauge == Basics.Coulomb
            # This reduced matrix element has been re-worked due to Johnson (2007)
            JohnsonI = AngularMomentum.JohnsonI(kapa, kapb, AngularJ64(mp.L))
            wr       = (1 - mp.L/(2mp.L+1)) * RadialIntegrals.GrantILplus(mp.L-1, q, a, b, grid::Radial.Grid)  -
                       mp.L/(2mp.L+1) * RadialIntegrals.GrantILplus(mp.L+1, q, a, b, grid::Radial.Grid)
            wr       = -(kapa-kapb) / (mp.L+1) * wr
            wr       = wr  +  mp.L/(2mp.L+1) * RadialIntegrals.GrantILminus(mp.L-1, q, a, b, grid::Radial.Grid)
            wr       = wr  +  mp.L/(2mp.L+1) * RadialIntegrals.GrantILminus(mp.L+1, q, a, b, grid::Radial.Grid)
            wa       = JohnsonI * wr
        else     error("stop a")
        end
        ##x wa = wa * sqrt( (2mp.L+1)*(mp.L+1)/mp.L )
        ##x wa = conj(wa)

        return( wa )
    end

    """
    `InteractionStrength.MbaEmissionAndrey(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)`
        ... to compute the (single-electron reduced matrix element) interaction strength <b || O^(Mp, emission) || a>
            for the interaction with the Mp multipole component of the radiation field and the transition frequency omega,
            and within the given gauge. A value::Float64 is returned. At present, only the magnetic matrix elements are
            implemented.
            This procedure has been worked out with Andrey but is currently not in use.
    """
    function MbaEmissionAndreyOld(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)
        kapa = a.subshell.kappa;   kapb = b.subshell.kappa;    q = omega / Defaults.getDefaults("speed of light: c")
        ja   = Basics.subshell_j(a.subshell);   jb   = Basics.subshell_j(b.subshell);
        #
        if       gauge == Basics.Magnetic
            wa = -1.0im * sqrt( (2mp.L+1) * (Basics.subshell_2j(a.subshell)+1) ) / sqrt(4pi) / sqrt(mp.L * (mp.L+1)) * (-1)^mp.L * (kapb + kapa)
            wa = wa * AngularMomentum.ClebschGordan(ja, AngularM64(1//2), AngularJ64(mp.L), AngularM64(0), jb, AngularM64(1//2))
            wa = wa * RadialIntegrals.GrantILplus(mp.L, q, a, b, grid::Radial.Grid)
            wa = conj(wa)
        #
        elseif   gauge == Basics.Babushkin
            wa = 0.
        #
        elseif   gauge == Basics.Coulomb
            wa = 0.
        else     error("stop a")
        end

        return( wa )
    end

    """
    `InteractionStrength.MbaEmissionMigdalek(cp::CorePolarization, b::Orbital, a::Orbital, grid::Radial.Grid)`
        ... to compute the (single-electron reduced matrix element) interaction strength <b || O^(E1, emission with core-polarization) || a>
            in length gauge. A value::Float64 is returned.
    """
    function MbaEmissionMigdalek(cp::CorePolarization, b::Orbital, a::Orbital, grid::Radial.Grid)
        mp   = E1;    kapa = a.subshell.kappa;   kapb = b.subshell.kappa
        ja   = Basics.subshell_j(a.subshell);   jb    = Basics.subshell_j(b.subshell);
        #
        JohnsonI = AngularMomentum.JohnsonI(kapa, kapb, AngularJ64(mp.L))
        wr       = RadialIntegrals.GrantJL_cp(mp.L, 0., a, b, grid::Radial.Grid, cp::CorePolarization)
        @show "MbaEmissionMigdalek", wr
        wa       = JohnsonI * wr

        return( wa )
    end

    """
    `InteractionStrength.multipoleTransition(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)`
        ... to compute the (single-electron reduced matrix element) multipole-transition interaction strength
            <b || T^(Mp, absorption) || a> due to Johnson (2007) for the interaction with the Mp multipole component of the
            radiation field and the transition frequency omega, and within the given gauge. A value::Float64 is returned.
    """
    function multipoleTransition(mp::EmMultipole, gauge::EmGauge, omega::Float64, b::Orbital, a::Orbital, grid::Radial.Grid)
        function besselPrime_jl(L::Int64, x::Float64)    return( GSL.sf_bessel_jl(L-1, x) - (L+1)/x * GSL.sf_bessel_jl(L, x) )       end

        kapa = a.subshell.kappa;   kapb = b.subshell.kappa;    q = omega / Defaults.getDefaults("speed of light: c")
        mtp  = min(size(a.P, 1), size(b.P, 1))
        !(grid.meshType == Radial.MeshGL())  &&  error("Only for Radial.MeshGL() implemented so far.")
        #
        if       gauge == Basics.Magnetic
            ChengI = AngularMomentum.ChengI(-kapa, kapb, AngularJ64(mp.L));   if  abs(ChengI) < 1.0e-10  return( 0. )   end
            wa = Complex(0.)
            for  i = 2:mtp
                wa = wa + (kapa+kapb) / (mp.L+1) * GSL.sf_bessel_jl(mp.L, q * grid.r[i]) * (a.P[i] * b.Q[i] + a.Q[i] * b.P[i]) * grid.wr[i]
            end
            wa = ChengI * wa
            #
        elseif   gauge == Basics.Velocity
            ChengI = AngularMomentum.ChengI(kapa, kapb, AngularJ64(mp.L));    if  abs(ChengI) < 1.0e-10  return( 0. )   end
            wa = Complex(0.)
            for  i = 2:mtp
                wa = wa - (kapa-kapb) / (mp.L+1) *
                          ( besselPrime_jl(mp.L, q * grid.r[i]) + GSL.sf_bessel_jl(mp.L, q * grid.r[i]) / (q * grid.r[i]) ) *
                          (a.P[i] * b.Q[i] + a.Q[i] * b.P[i]) * grid.wr[i]
                wa = wa + mp.L * GSL.sf_bessel_jl(mp.L, q * grid.r[i]) / (q * grid.r[i]) * (a.P[i] * b.Q[i] - a.Q[i] * b.P[i]) * grid.wr[i]
            end
            wa = ChengI * wa
            #
        elseif   gauge == Basics.Length
            ChengI = AngularMomentum.ChengI(kapa, kapb, AngularJ64(mp.L));    if  abs(ChengI) < 1.0e-10  return( 0. )   end
            wa = Complex(0.)
            for  i = 2:mtp
                wa = wa + GSL.sf_bessel_jl(mp.L, q * grid.r[i]) * (a.P[i] * b.P[i] + a.Q[i] * b.Q[i]) * grid.wr[i]
                        + GSL.sf_bessel_jl(mp.L+1, q * grid.r[i]) *
                          ( (kapa-kapb) / (mp.L+1) * (a.P[i] * b.Q[i] + a.Q[i] * b.P[i]) +
                            (a.P[i] * b.Q[i] - a.Q[i] * b.P[i]) ) * grid.wr[i]
            end
            wa = ChengI * wa
            #
        else     error("stop a")
        end

        return( wa )
    end

    """
    `InteractionStrength.schiffMoment(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid)`
        ... computes the  <a|| h^(Schiff-moment) ||b>  reduced matrix element of the Schiff-moment Hamiltonian for orbital
            functions a, b and for the nuclear density as given by the nuclear model. A value::Float64 is returned.
    """
    function schiffMoment(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid)
        printstyled("\nWarning -- InteractionStrength.schiffMoment():: Not yet implemented.", color=:cyan)
        wb = 1.0 + 2.0im
        return( wb )
    end

    """
    `InteractionStrength.weakCharge(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid)`
        ... computes the  <a|| h^(weak-charge) ||b>  reduced matrix element of the weak-charge Hamiltonian for orbital functions
            a, b and for the nuclear density as given by the nuclear model. A value::Float64 is returned.
    """
    function weakCharge(a::Orbital, b::Orbital, nm::Nuclear.Model, grid::Radial.Grid)
        printstyled("\nWarning -- InteractionStrength.weakCharge():: Not yet implemented.", color=:cyan)
        wb = 1.0 + 2.0im
        return( wb )
    end

    """
    `InteractionStrength.XL_Breit_WO(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid)`
        ... computes the the effective Breit interaction strengths X^L_Breit (abcd) for given rank L and orbital functions
            a, b, c and d  at the given grid but without optimization. A value::Float64 is returned. At present, only the zero-frequency
            Breit interaction is taken into account.
    """
    function XL_Breit_WO(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid)
        ja2 = Basics.subshell_2j(a.subshell)
        jb2 = Basics.subshell_2j(b.subshell)
        jc2 = Basics.subshell_2j(c.subshell)
        jd2 = Basics.subshell_2j(d.subshell)
        if  AngularMomentum.triangularDelta(ja2+1,jc2+1,L+L+1) * AngularMomentum.triangularDelta(jb2+1,jd2+1,L+L+1) == 0   ||  L == 0
            return( 0. )
        end
        #
        xcList   = XL_Breit0_coefficients(L,a,b,c,d)
        ##x println("\nBreit strength: L=$L, nx = $(length(xcList)) ")
        XL_Breit = XL_Breit0_densities(xcList, grid)
        #
        return( XL_Breit )
    end

    """
    `InteractionStrength.XL_Breit_reset_storage(keep::Bool; printout::Bool=false)`
        ... resets the global storage of XL_Breit interaction strength; nothing is returned.
    """
    function XL_Breit_reset_storage(keep::Bool; printout::Bool=false)
        if  keep
            if  printout   println("  reset GBL_Storage_XL_Breit storage ...")    end
            global GBL_Storage_XL_Breit = Dict{String, Float64}()
        else
        end
        return( nothing )
    end

    """
    `InteractionStrength.XL_Breit(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid; keep::Bool=false)`
        ... computes the the effective Breit interaction strengths X^L_Breit (abcd) for given rank L and orbital functions
            a, b, c and d  at the given grid. For keep=true, the procedure looks up the (global) directory GBL_Storage_XL_Coulomb
            and returns the corresponding value without re-calculation of the interaction strength; it also 'stores' the calculated
            value if not yet included. For keep=false, the interaction strength is always computed on-fly. A value::Float64 is returned.
            At present, only the zero-frequency Breit interaction is taken into account.
    """
    function XL_Breit(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid; keep::Bool=false)
        global GBL_Storage_XL_Breit
        ja2 = Basics.subshell_2j(a.subshell)
        jb2 = Basics.subshell_2j(b.subshell)
        jc2 = Basics.subshell_2j(c.subshell)
        jd2 = Basics.subshell_2j(d.subshell)
        if  AngularMomentum.triangularDelta(ja2+1,jc2+1,L+L+1) * AngularMomentum.triangularDelta(jb2+1,jd2+1,L+L+1) == 0   ||  L == 0
            return( 0. )
        end

        # Now distiguish due to the optional argument keep
        if  keep
            sa = "XL" * string(L) * " " * string(a.subshell) * string(b.subshell) * string(c.subshell) * string(d.subshell)
            if haskey(GBL_Storage_XL_Breit, sa )
                XL_Breit = GBL_Storage_XL_Breit[sa]
            else
                xcList   = XL_Breit0_coefficients(L,a,b,c,d)
                XL_Breit = XL_Breit0_densities(xcList, grid)
                global GBL_Storage_XL_Breit = Base.merge(GBL_Storage_XL_Breit, Dict( sa => XL_Breit))
            end
        else
            xcList   = XL_Breit0_coefficients(L,a,b,c,d)
            XL_Breit = XL_Breit0_densities(xcList, grid)
        end
        #

        return( XL_Breit )
    end

    """
    `InteractionStrength.XL_BreitDamped(tau::Float64, L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid)`
        ... computes the the effective Breit interaction strengths X^L_Breit (abcd) for given rank L and orbital functions
            a, b, c and d  at the given grid. A value::Float64 is returned. At present, only the zero-frequency Breit
            interaction is taken into account.
    """
    function XL_BreitDamped(tau::Float64, L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid)
        error("stop a")
    end

    """
    `InteractionStrength.XL_Breit0_coefficients(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital)`
        ... evaluates the combinations and pre-coefficients for the zero-frequency Breit interaction
            X^L_Breit (omega=0.; abcd) for given rank L and orbital functions a, b, c and d. A list of coefficients
            xcList::Array{XLCoefficient,1} is returned.
    """
    function XL_Breit0_coefficients(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital)
        xcList = XLCoefficient[]

        la = Basics.subshell_l(a.subshell);    ja2 = Basics.subshell_2j(a.subshell)
        lb = Basics.subshell_l(b.subshell);    jb2 = Basics.subshell_2j(b.subshell)
        lc = Basics.subshell_l(c.subshell);    jc2 = Basics.subshell_2j(c.subshell)
        ld = Basics.subshell_l(d.subshell);    jd2 = Basics.subshell_2j(d.subshell)

        xc = AngularMomentum.CL_reduced_me(a.subshell, L, c.subshell) * AngularMomentum.CL_reduced_me(b.subshell, L, d.subshell)
        if   rem(L,2) == 1    xc = - xc                end
        if   abs(xc)  <  1.0e-10    return( xcList )   end

        # Consider the individual contributions from sum_nu and sum_mu. First, take T^(nu,L)_mu = R^(nu,L)_mu
        nu = L - 1

        if  rem(la+lc+nu,2) == 1   &&   rem(lb+ld+nu,2) == 1   &&   L != 0
            wa = (L+1) / ( L*(L+L-1)*(L+L+1) )
            # mu = 1
            xcc = xc * wa * (c.subshell.kappa-a.subshell.kappa+L) * (d.subshell.kappa-b.subshell.kappa+L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, a, b, c, d, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, b, a, d, c, xcc) )   end
            # mu = 2
            xcc = xc * wa * (c.subshell.kappa-a.subshell.kappa-L) * (d.subshell.kappa-b.subshell.kappa-L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, c, d, a, b, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, d, c, b, a, xcc) )   end
            # mu = 3
            xcc = xc * wa * (c.subshell.kappa-a.subshell.kappa+L) * (d.subshell.kappa-b.subshell.kappa-L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, a, d, c, b, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, d, a, b, c, xcc) )   end
            # mu = 4
            xcc = xc * wa * (c.subshell.kappa-a.subshell.kappa-L) * (d.subshell.kappa-b.subshell.kappa+L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, c, b, a, d, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, b, c, d, a, xcc) )   end
        end
        #
        nu = L

        if  rem(la+lc+nu,2) == 1   &&   rem(lb+ld+nu,2) == 1   &&   L != 0
            wa = - (a.subshell.kappa + c.subshell.kappa) * (b.subshell.kappa + d.subshell.kappa) / (L*(L+1))
            # mu = 1
            xcc = xc * wa
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, a, b, c, d, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, b, a, d, c, xcc) )   end
            # mu = 2
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, c, d, a, b, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, d, c, b, a, xcc) )   end
            # mu = 3
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, a, d, c, b, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, d, a, b, c, xcc) )   end
            # mu = 4
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, c, b, a, d, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, b, c, d, a, xcc) )   end
        end
        #
        nu = L + 1

        if  rem(la+lc+nu,2) == 1   &&   rem(lb+ld+nu,2) == 1   &&   L != 0
            wa =  L / ( (L+1)*(L+L+1)*(L+L+3) )
            # mu = 1
            xcc = xc * wa * (c.subshell.kappa - a.subshell.kappa - L - 1) * (d.subshell.kappa - b.subshell.kappa - L - 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, a, b, c, d, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, b, a, d, c, xcc) )   end
            # mu = 2
            xcc = xc * wa * (c.subshell.kappa - a.subshell.kappa + L + 1) * (d.subshell.kappa - b.subshell.kappa + L + 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, c, d, a, b, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, d, c, b, a, xcc) )   end
            # mu = 3
            xcc = xc * wa * (c.subshell.kappa - a.subshell.kappa - L - 1) * (d.subshell.kappa - b.subshell.kappa + L + 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, a, d, c, b, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, d, a, b, c, xcc) )   end
            # mu = 4
            xcc = xc * wa * (c.subshell.kappa - a.subshell.kappa + L + 1) * (d.subshell.kappa - b.subshell.kappa - L - 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, c, b, a, d, xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, nu, b, c, d, a, xcc) )   end
        end

        # Add contributions of the S^k_mu integrals
        if  rem(la+lc+L-1,2) == 1   &&   rem(lb+ld+L+1,2) == 1
            # mu = 1
            wb =  1 / ( (L+L+1)*(L+L+1) )
            xcc = xc * wb * (c.subshell.kappa - a.subshell.kappa + L) * (d.subshell.kappa - b.subshell.kappa - L - 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, b, a, d, c,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, b, a, d, c,  -(L+L+1) / 2 * xcc) )   end
            # mu = 2
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa + L) * (c.subshell.kappa - a.subshell.kappa - L - 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, a, b, c, d,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, a, b, c, d,  -(L+L+1) / 2 * xcc) )   end
            # mu = 3
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa + L + 1) * (c.subshell.kappa - a.subshell.kappa - L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, d, c, b, a,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, d, c, b, a,  -(L+L+1) / 2 * xcc) )   end
            # mu = 4
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa - L) * (c.subshell.kappa - a.subshell.kappa + L + 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, c, d, a, b,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, c, d, a, b,  -(L+L+1) / 2 * xcc) )   end
            # mu = 5
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa + L + 1) * (c.subshell.kappa - a.subshell.kappa + L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, d, a, b, c,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, d, a, b, c,  -(L+L+1) / 2 * xcc) )   end
            # mu = 6
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa - L) * (c.subshell.kappa - a.subshell.kappa - L - 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, a, d, c, b,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, a, d, c, b,  -(L+L+1) / 2 * xcc) )   end
            # mu = 7
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa - L - 1) * (c.subshell.kappa - a.subshell.kappa - L)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, b, c, d, a,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, b, c, d, a,  -(L+L+1) / 2 * xcc) )   end
            # mu = 8
            xcc = xc * wb * (d.subshell.kappa - b.subshell.kappa + L) * (c.subshell.kappa - a.subshell.kappa + L + 1)
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L+1, c, b, a, d,   (L+L+1) / 2 * xcc) )   end
            if  abs(xcc) > 1.0e-10   push!( xcList, XLCoefficient(5, L-1, c, b, a, d,  -(L+L+1) / 2 * xcc) )   end
        end

        return( xcList )
    end

    """
    `InteractionStrength.XL_Breit0_densities(xcList::Array{XLCoefficient,1}, grid::Radial.Grid)`
        ... computes the the effective Breit interaction strengths X^L,0_Breit (abcd) for given rank L and a list of
            orbital functions a, b, c, d and angular coefficients at the given grid. A value::Float64 is returned.
            At present, only the zero-frequency Breit interaction is taken into account.
    """
    function XL_Breit0_densities(xcList::Array{XLCoefficient,1}, grid::Radial.Grid)

        if  grid.meshType == Radial.MeshGL()
            wa = 0.
            for  xc  in  xcList  ## [end:end]
                # Use the minimal extent of any involved orbitals; this need to be improved
                mtp_ac = min(size(xc.a.P, 1), size(xc.c.P, 1));   mtp_bd = min(size(xc.b.P, 1), size(xc.d.P, 1))
                for  r = 2:mtp_ac
                    for  s = 2:mtp_bd
                        wc = xc.coeff * grid.wr[r] * grid.wr[s]
                        if      s > r   continue
                        elseif  s == r
                            wa = wa + wc * (xc.a.P[r] * xc.c.Q[r]) * (grid.r[s]^xc.nu) / (grid.r[r]^(xc.nu+1)) * (xc.b.P[s] * xc.d.Q[s]) / 2.0
                        else
                            wa = wa + wc * (xc.a.P[r] * xc.c.Q[r]) * (grid.r[s]^xc.nu) / (grid.r[r]^(xc.nu+1)) * (xc.b.P[s] * xc.d.Q[s])
                        end
                    end
                end
            end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `InteractionStrength.matrixL_Coulomb(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, primitives::Bsplines.Primitives)`
        ... computes the partly-contracted (effective) Coulomb interaction matrices M^L_Coulomb (abcd) for given rank L and orbital functions
            a, b, c and d at the given grid. The matrix M^L is defined for the primitives and contracted over the two orbitals
            b, d (for a=c) or  b, c (for a=d).  An error message is issued if a != c && a != d. A matrix::Array{Float64,2} is returned.
    """
    function matrixL_Coulomb(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, primitives::Bsplines.Primitives)
        grid = primitives.grid;   nsL = primitives.grid.nsL;    nsS = primitives.grid.nsS
        wm = zeros( nsL+nsS, nsL+nsS )
        # Test for the triangular-delta conditions and calculate the reduced matrix elements of the C^L tensors
        la = Basics.subshell_l(a.subshell);    ja2 = Basics.subshell_2j(a.subshell)
        lb = Basics.subshell_l(b.subshell);    jb2 = Basics.subshell_2j(b.subshell)
        lc = Basics.subshell_l(c.subshell);    jc2 = Basics.subshell_2j(c.subshell)
        ld = Basics.subshell_l(d.subshell);    jd2 = Basics.subshell_2j(d.subshell)

        if  AngularMomentum.triangularDelta(ja2+1,jc2+1,L+L+1) * AngularMomentum.triangularDelta(jb2+1,jd2+1,L+L+1) == 0   ||
            rem(la+lc+L,2) == 1   ||   rem(lb+ld+L,2) == 1
            return( wm )
        end
        xc = AngularMomentum.CL_reduced_me(a.subshell, L, c.subshell) * AngularMomentum.CL_reduced_me(b.subshell, L, d.subshell)
        if   rem(L,2) == 1    xc = - xc    end

        if  a.subshell == c.subshell
            # Direct interaction; contract the full interaction array over the orbitals b and d
            for  i = 1:nsL
                for  k = 1:nsL
                    Ba = primitives.bsplinesL[i].bs;    Bc = primitives.bsplinesL[k].bs
                    wm[i,k] = RadialIntegrals.SlaterRkComponent_2dim(L, Ba, b.P, Bc, d.P, grid) +
                              RadialIntegrals.SlaterRkComponent_2dim(L, Ba, b.Q, Bc, d.Q, grid)
                end
            end
            for  i = 1:nsS
                for  k = 1:nsS
                    Ba = primitives.bsplinesS[i].bs;    Bc = primitives.bsplinesS[k].bs
                    wm[nsL+i,nsL+k] = RadialIntegrals.SlaterRkComponent_2dim(L, Ba, b.P, Bc, d.P, grid) +
                                      RadialIntegrals.SlaterRkComponent_2dim(L, Ba, b.Q, Bc, d.Q, grid)
                end
            end
        elseif true
            println("Skip exchange integrals")
            return( wm )
        elseif  a.subshell == d.subshell
            # Exchange interaction; contract the full interaction array over the orbitals b and c
            for  i = 1:nsL
                for  k = 1:nsL
                    Ba = primitives.bsplinesL[i].bs;    Bd = primitives.bsplinesL[k].bs
                    wm[i,k] = RadialIntegrals.SlaterRkComponent_2dim(L, Ba, b.P, c.P, Bd, grid)
                end
                for  k = 1:nsS
                    Ba = primitives.bsplinesL[i].bs;    Bd = primitives.bsplinesS[k].bs
                    wm[i,nsL+k] = RadialIntegrals.SlaterRkComponent_2dim(L, Ba, b.P, c.Q, Bd, grid)
                end
            end
            for  i = 1:nsS
                for  k = 1:nsL
                    Ba = primitives.bsplinesS[i].bs;    Bd = primitives.bsplinesL[k].bs
                    wm[nsL+i,k] = RadialIntegrals.SlaterRkComponent_2dim(L, Ba, b.Q, c.P, Bd, grid)
                end
                for  k = 1:nsS
                    Ba = primitives.bsplinesS[i].bs;    Bd = primitives.bsplinesS[k].bs
                    wm[nsL+i,nsL+k] = RadialIntegrals.SlaterRkComponent_2dim(L, Ba, b.Q, c.Q, Bd, grid)
                end
            end
        else    error("stop d")
        end

        return( wm )
    end

    #==
    """
    `InteractionStrength.XS_Coulomb(large::Bool, largep::Bool, coeffs::Array{JAC.AngularCoefficientsRatip2013.AngularVcoeff,1},
                                    orbitals::Dict{Subshell, Orbital}, grid::Radial.Grid; exchange::Bool=false)`
        ... computes the effective XS_Coulomb interaction function for given components (large, largep) and a set of (Coulomb-like) angular coefficients.
            The procedures assumes that all coefficients V^K (abcd) fulfill  a==c || a == d, and a error is issued if this is not the case.
            For the remaining two orbitals b, d (or c), the correct components are added due to boolean selectors.
            The XS_Coulomb interaction function is defined by

                S(r) = sum_k  Int ds  U(r,s)  V^k(abcd) b.P(large) cd.P(largep)

            where b.P(large) refers to b.P of orbital b if large=true and b.Q if large= false, and analogue for cd.P. A Sfunc::Float64[] is returned.
    """
    function XS_Coulomb(large::Bool, largep::Bool, coeffs::Array{JAC.AngularCoefficientsRatip2013.AngularVcoeff,1},
                        orbitals::Dict{Subshell, Orbital}, grid::Radial.Grid; exchange::Bool=false)
        function ul(k::Int64, r::Float64, s::Float64) :: Float64
            if     r <= s    return( r^k/s^(k+1) )
            elseif r > s     return( s^k/r^(k+1) )
            end
        end
        #
        Sfunc = zeros( grid.NoPoints );   na = 1
        for  r = 2:grid.NoPoints
            # Determine integrant over s
            ws = zeros( grid.NoPoints )
            for  coeff in  coeffs
                a = orbitals[coeff.a];   b = orbitals[coeff.b];   c = orbitals[coeff.c];   d = orbitals[coeff.d];   L = coeff.nu
                # Determine the weight xc of this coefficient
                la = Basics.subshell_l(a.subshell);    ja2 = Basics.subshell_2j(a.subshell)
                lb = Basics.subshell_l(b.subshell);    jb2 = Basics.subshell_2j(b.subshell)
                lc = Basics.subshell_l(c.subshell);    jc2 = Basics.subshell_2j(c.subshell)
                ld = Basics.subshell_l(d.subshell);    jd2 = Basics.subshell_2j(d.subshell)

                if  AngularMomentum.triangularDelta(ja2+1,jc2+1,L+L+1) * AngularMomentum.triangularDelta(jb2+1,jd2+1,L+L+1) == 0   ||
                    rem(la+lc+L,2) == 1   ||   rem(lb+ld+L,2) == 1      error("stop a")
                end
                xc = AngularMomentum.CL_reduced_me(a.subshell, L, c.subshell) * AngularMomentum.CL_reduced_me(b.subshell, L, d.subshell)
                if   rem(L,2) == 1    xc = - xc    end
                # Divide total coefficient by occupation w_a of the considered shell (DF equations vs. energy functional)
                xc = xc / (ja2 + 1)
                #
                ##x if na < 6   na = na + 1;    @show L, coeff.a, coeff.b, coeff.c, coeff.d, coeff.V * xc  end
                #
                # Now add the contributions for this coefficient
                if      coeff.a == coeff.c  &&  large  &&  largep
                    mtp_bd = min(size(b.P, 1), size(d.P, 1));    for  s = 2:mtp_bd   ws[s] = ws[s] + coeff.V * xc * b.P[s] * d.P[s]* ul(L, grid.r[r], grid.r[s])    end
                elseif  coeff.a == coeff.c  &&  !(large)  &&  !(largep)
                    mtp_bd = min(size(b.Q, 1), size(d.Q, 1));    for  s = 2:mtp_bd   ws[s] = ws[s] + coeff.V * xc * b.Q[s] * d.Q[s]* ul(L, grid.r[r], grid.r[s])    end
                elseif  coeff.a == coeff.d  &&  exchange  &&  large  &&  largep
                    mtp_bc = min(size(b.P, 1), size(c.P, 1));    for  s = 2:mtp_bc   ws[s] = ws[s] + coeff.V * xc * b.P[s] * c.P[s]* ul(L, grid.r[r], grid.r[s])    end
                elseif  coeff.a == coeff.d  &&  exchange  &&  !(large)  &&  largep
                    mtp_bc = min(size(b.Q, 1), size(c.P, 1));    for  s = 2:mtp_bc   ws[s] = ws[s] + coeff.V * xc * b.Q[s] * c.P[s]* ul(L, grid.r[r], grid.r[s])    end
                elseif  coeff.a == coeff.d  &&  exchange  &&  large     &&  !(largep)
                     mtp_bc = min(size(b.P, 1), size(c.Q, 1));    for  s = 2:mtp_bc   ws[s] = ws[s] + coeff.V * xc * b.P[s] * c.Q[s]* ul(L, grid.r[r], grid.r[s])    end
                elseif  coeff.a == coeff.d  &&  exchange  &&  !(large)  &&  !(largep)
                    mtp_bc = min(size(b.Q, 1), size(c.Q, 1));    for  s = 2:mtp_bc   ws[s] = ws[s] + coeff.V * xc * b.Q[s] * c.Q[s]* ul(L, grid.r[r], grid.r[s])    end
                elseif  coeff.a == coeff.d  &&  !(exchange)
                elseif  coeff.a == coeff.c  &&  (large     &&  !(largep)   ||   !(large)  &&  largep)
                else    @show L, coeff.a, coeff.b, coeff.c, coeff.d, coeff.V * xc
                end
            end
            #
            for  s = 2:grid.NoPoints   Sfunc[r] = Sfunc[r]  +  ws[s] * grid.wr[s]  end
        end

        return( Sfunc )
    end  ==#

    """
    `InteractionStrength.XL_Coulomb_WO(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid)`
        ... computes the the effective Coulomb interaction strengths X^L_Coulomb (abcd) for given rank L and orbital functions
            a, b, c and d at the given grid but without optimization. A value::Float64 is returned.
    """
    function XL_Coulomb_WO(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid)
        # Test for the triangular-delta conditions and calculate the reduced matrix elements of the C^L tensors
        la = Basics.subshell_l(a.subshell);    ja2 = Basics.subshell_2j(a.subshell)
        lb = Basics.subshell_l(b.subshell);    jb2 = Basics.subshell_2j(b.subshell)
        lc = Basics.subshell_l(c.subshell);    jc2 = Basics.subshell_2j(c.subshell)
        ld = Basics.subshell_l(d.subshell);    jd2 = Basics.subshell_2j(d.subshell)

        if  AngularMomentum.triangularDelta(ja2+1,jc2+1,L+L+1) * AngularMomentum.triangularDelta(jb2+1,jd2+1,L+L+1) == 0   ||
            rem(la+lc+L,2) == 1   ||   rem(lb+ld+L,2) == 1
            return( 0. )
        end
        xc = AngularMomentum.CL_reduced_me(a.subshell, L, c.subshell) * AngularMomentum.CL_reduced_me(b.subshell, L, d.subshell)
        if   rem(L,2) == 1    xc = - xc    end

        XL_Coulomb = xc * RadialIntegrals.SlaterRk_2dim_WO(L, a, b, c, d, grid)
        ##x XL_Coulomb = xc * RadialIntegrals.SlaterRk_new(L, a, b, c, d, grid)
        return( XL_Coulomb )
    end

    """
    `InteractionStrength.XL_Coulomb_reset_storage(keep::Bool; printout::Bool=false)`
        ... resets the global storage of XL_Coulomb interaction strength; nothing is returned.
    """
    function XL_Coulomb_reset_storage(keep::Bool; printout::Bool=false)
        if  keep
            if printout     println(">> Reset GBL_Storage_XL_Coulomb storage.")     end
            global GBL_Storage_XL_Coulomb = Dict{String, Float64}()
        else
        end
        return( nothing )
    end

    """
    `InteractionStrength.XL_Coulomb(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid; keep::Bool=false)`
        ... computes the the effective Coulomb interaction strengths X^L_Coulomb (abcd) for given rank L and orbital functions
            a, b, c and d at the given grid. For keep=true, the procedure looks up the (global) directory GBL_Storage_XL_Coulomb
            and returns the corresponding value without re-calculation of the interaction strength; it also 'stores' the calculated
            value if not yet included. For keep=false, the interaction strength is always computed on-fly. A value::Float64 is
            returned.
    """
    function XL_Coulomb(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid; keep::Bool=false)
        global GBL_Storage_XL_Coulomb
        # Test for the triangular-delta conditions and calculate the reduced matrix elements of the C^L tensors
        la = Basics.subshell_l(a.subshell);    ja2 = Basics.subshell_2j(a.subshell)
        lb = Basics.subshell_l(b.subshell);    jb2 = Basics.subshell_2j(b.subshell)
        lc = Basics.subshell_l(c.subshell);    jc2 = Basics.subshell_2j(c.subshell)
        ld = Basics.subshell_l(d.subshell);    jd2 = Basics.subshell_2j(d.subshell)

        if  AngularMomentum.triangularDelta(ja2+1,jc2+1,L+L+1) * AngularMomentum.triangularDelta(jb2+1,jd2+1,L+L+1) == 0   ||
            rem(la+lc+L,2) == 1   ||   rem(lb+ld+L,2) == 1
            return( 0. )
        end

        # Now distiguish due to the optional argument keep
        if  keep
            sa = "XL" * string(L) * " " * string(a.subshell) * string(b.subshell) * string(c.subshell) * string(d.subshell)
            if haskey(GBL_Storage_XL_Coulomb, sa )
                XL_Coulomb = GBL_Storage_XL_Coulomb[sa]
            else
                XL_Coulomb = InteractionStrength.XL_Coulomb(L::Int64, a, b, c, d, grid)
                ## global GBL_Storage_XL_Coulomb = Base.merge(GBL_Storage_XL_Coulomb, Dict( sa => XL_Coulomb))
                global GBL_Storage_XL_Coulomb[sa] = XL_Coulomb
            end
        else
            xc = AngularMomentum.CL_reduced_me(a.subshell, L, c.subshell) * AngularMomentum.CL_reduced_me(b.subshell, L, d.subshell)
            if   rem(L,2) == 1    xc = - xc    end
            XL_Coulomb = xc * RadialIntegrals.SlaterRk_2dim(L, a, b, c, d, grid)

            ##x XL_Coulomb = XL_Coulomb* (-1)^( (ja2+jb2+jc2+jd2)/2 )
        end

        return( XL_Coulomb )
    end

    """
    `InteractionStrength.XL_CoulombDamped(tau::Float64, L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid)`
        ... computes the the effective Coulomb interaction strengths X^L_Coulomb (abcd) for given rank L and orbital functions
            a, b, c and d at the given grid. A value::Float64 is returned.
    """
    function XL_CoulombDamped(tau::Float64, L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid)
        # Test for the triangular-delta conditions and calculate the reduced matrix elements of the C^L tensors
        la = Basics.subshell_l(a.subshell);    ja2 = Basics.subshell_2j(a.subshell)
        lb = Basics.subshell_l(b.subshell);    jb2 = Basics.subshell_2j(b.subshell)
        lc = Basics.subshell_l(c.subshell);    jc2 = Basics.subshell_2j(c.subshell)
        ld = Basics.subshell_l(d.subshell);    jd2 = Basics.subshell_2j(d.subshell)

        if  AngularMomentum.triangularDelta(ja2+1,jc2+1,L+L+1) * AngularMomentum.triangularDelta(jb2+1,jd2+1,L+L+1) == 0   ||
            rem(la+lc+L,2) == 1   ||   rem(lb+ld+L,2) == 1
            return( 0. )
        end
        xc = AngularMomentum.CL_reduced_me(a.subshell, L, c.subshell) * AngularMomentum.CL_reduced_me(b.subshell, L, d.subshell)
        if   rem(L,2) == 1    xc = - xc    end

        XL_Coulomb = xc * RadialIntegrals.SlaterRk_2dim_Damped(tau::Float64, L, a, b, c, d, grid)
        ##x XL_Coulomb = xc * RadialIntegrals.SlaterRk_new(L, a, b, c, d, grid)
        return( XL_Coulomb )
    end

    """
    `InteractionStrength.XL_Coulomb_DH(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid, lambda::Float64)`
        ... computes the the effective Coulomb-Debye-Hückel interaction strengths X^L_Coulomb_DH (abcd) for given rank L and
            orbital functions a, b, c and d at the given grid and for the given screening parameter lambda. A value::Float64 is
            returned.
    """
    function XL_Coulomb_DH(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, grid::Radial.Grid, lambda::Float64)
        # Test for the triangular-delta conditions and calculate the reduced matrix elements of the C^L tensors
        la = Basics.subshell_l(a.subshell);    ja2 = Basics.subshell_2j(a.subshell)
        lb = Basics.subshell_l(b.subshell);    jb2 = Basics.subshell_2j(b.subshell)
        lc = Basics.subshell_l(c.subshell);    jc2 = Basics.subshell_2j(c.subshell)
        ld = Basics.subshell_l(d.subshell);    jd2 = Basics.subshell_2j(d.subshell)

        if  AngularMomentum.triangularDelta(ja2+1,jc2+1,L+L+1) * AngularMomentum.triangularDelta(jb2+1,jd2+1,L+L+1) == 0   ||
            rem(la+lc+L,2) == 1   ||   rem(lb+ld+L,2) == 1
            return( 0. )
        end
        xc = AngularMomentum.CL_reduced_me(a.subshell, L, c.subshell) * AngularMomentum.CL_reduced_me(b.subshell, L, d.subshell)
        if   rem(L,2) == 1    xc = - xc    end

        XL_Coulomb_DH = xc * RadialIntegrals.SlaterRk_DebyeHueckel_2dim(L, a, b, c, d, grid, lambda)
        return( XL_Coulomb_DH )
    end

    """
    `InteractionStrength.XL_plasma_ionSphere(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, lambda::Float64)`
        ... computes the effective interaction strengths X^L_ion-sphere (abcd) for given rank L and orbital functions
            a, b, c and d and for the plasma parameter lambda. A value::Float64 is returned.  **Not yet implemented !**
    """
    function XL_plasma_ionSphere(L::Int64, a::Orbital, b::Orbital, c::Orbital, d::Orbital, lambda::Float64)
        error("Not yet implemented")
    end

    """
    `InteractionStrength.X_smsA(a::Orbital, b::Orbital, c::Orbital, d::Orbital, nm::Nuclear.Model, grid::Radial.Grid)`
        ... computes the the effective interaction strengths X^1_sms,A (abcd) for fixed rank 1 and orbital functions
            a, b, c and d at the given grid. A value::Float64 is returned.
    """
    function X_smsA(a::Orbital, b::Orbital, c::Orbital, d::Orbital, nm::Nuclear.Model, grid::Radial.Grid)
        wa = AngularMomentum.CL_reduced_me_sms(a.subshell, 1, c.subshell) *
             AngularMomentum.CL_reduced_me_sms(b.subshell, 1, d.subshell) *
             RadialIntegrals.Vinti(a, c, grid) * RadialIntegrals.Vinti(b, d, grid) / 2
        ## println("**  <$(a.subshell) || Vinti || $(c.subshell)>  = $(RadialIntegrals.Vinti(a, c, grid)) " )
        ## println("**  <$(b.subshell) || Vinti || $(d.subshell)>  = $(RadialIntegrals.Vinti(b, d, grid)) " )
        return( wa )
    end

    """
    `InteractionStrength.X_smsB(a::Orbital, b::Orbital, c::Orbital, d::Orbital, nm::Nuclear.Model, grid::Radial.Grid)`
        ... computes the the effective interaction strengths X^1_sms,B (abcd) for fixed rank 1 and orbital functions
            a, b, c and d at the given grid. A value::Float64 is returned.
    """
    function X_smsB(a::Orbital, b::Orbital, c::Orbital, d::Orbital, nm::Nuclear.Model, grid::Radial.Grid)
        ##x println("")
        ##x @show AngularMomentum.CL_reduced_me_sms(b.subshell, 1, d.subshell)
        ##x @show RadialIntegrals.Vinti(b, d, grid)
        ##x @show RadialIntegrals.isotope_smsB(a, c, nm.Z, grid)
        wa = - AngularMomentum.CL_reduced_me_sms(b.subshell, 1, d.subshell) * RadialIntegrals.Vinti(b, d, grid) *
               RadialIntegrals.isotope_smsB(a, c, nm.Z, grid) / 2
        return( wa )
    end

    """
    `InteractionStrength.X_smsC(a::Orbital, b::Orbital, c::Orbital, d::Orbital, nm::Nuclear.Model, grid::Radial.Grid)`
        ... computes the the effective interaction strengths X^1_sms,C (abcd) for fixed rank 1 and orbital functions
            a, b, c and d at the given grid. A value::Float64 is returned.
    """
    function X_smsC(a::Orbital, b::Orbital, c::Orbital, d::Orbital, nm::Nuclear.Model, grid::Radial.Grid)
        wa = - AngularMomentum.CL_reduced_me_sms(b.subshell, 1, d.subshell) *
               AngularMomentum.CL_reduced_me_sms(a.subshell, 1, c.subshell) *
               RadialIntegrals.Vinti(b, d, grid) * RadialIntegrals.isotope_smsC(a, c, nm.Z, grid) / 2
        return( wa )
    end

    """
    `InteractionStrength.zeeman_Delta_n1(a::Orbital, b::Orbital, grid::Radial.Grid)`
        ... computes the <a|| Delta n^(1) ||b> reduced matrix element for the Zeeman-Schwinger contribution to the coupling
            to an external magnetic field for orbital functions a, b. A value::Float64 is returned.
    """
    function zeeman_Delta_n1(a::Orbital, b::Orbital, grid::Radial.Grid)
        # Use Andersson, Jönson (2008), CPC ... test for the proper definition of the C^L tensors.
        minusa = Subshell(1, -a.subshell.kappa)
        wa = (Defaults.getDefaults("electron g-factor") - 2) / 2. * (a.subshell.kappa + b.subshell.kappa + 1) *
             AngularMomentum.CL_reduced_me(minusa, 1, b.subshell) * RadialIntegrals.rkDiagonal(0, a, b, grid)
        return( wa )
    end

    """
    `InteractionStrength.zeeman_n1(a::Orbital, b::Orbital, grid::Radial.Grid)`
        ... computes the <a|| n^(1) ||b> reduced matrix element for the Zeeman coupling to an external magnetic field for
            orbital functions a, b. A value::Float64 is returned.
    """
    function zeeman_n1(a::Orbital, b::Orbital, grid::Radial.Grid)
        # Use Andersson, Jönson (2008), CPC ... test for the proper definition of the C^L tensors.
        minusa = Subshell(1, -a.subshell.kappa)
        wa = - AngularMomentum.CL_reduced_me(minusa, 1, b.subshell) / (2 * Defaults.getDefaults("alpha")) *
               RadialIntegrals.rkNonDiagonal(1, a, b, grid)
        return( wa )
    end

end # module


text:
using SDD
const SDD = SententialDecisionDiagrams

# set up vtree and manager
vtree = SDD.read_vtree("$(@__DIR__)/input/rotate-left.vtree")
manager = SDD.sdd_manager(vtree)

# construct the term X_1 ^ X_2 ^ X_3 ^ X_4
x = [SDD.literal(i,manager) for i in 1:5]
α = x[1] & x[2] & x[3] & x[4]

# to perform a rotate, we need the manager's vtree
manager_vtree = SDD.vtree(manager)
manager_vtree_right = SDD.right(manager_vtree)

println("saving vtree & sdd ...")
SDD.dot("$(@__DIR__)/output/before-rotate-vtree.dot", α)

# ref alpha (so it is not gc'd)
SDD.ref(α)

# garbage collect (no dead nodes when performing vtree operations)
println("dead sdd nodes = $(SDD.dead_count(manager))")
println("garbage collection ...")
SDD.garbage_collect(manager)
println("dead sdd nodes = $(SDD.dead_count(manager))")

println("left rotating ... ")
succeeded = SDD.rotate_left(manager_vtree_right,manager,0)
if succeeded == 1; println("succeeded!") else println("did not succeed!") end

# deref alpha, since ref's are no longer needed
SDD.deref(α)

# the root changed after rotation, so get the manager's vtree again
# this time using root_location
manager_vtree = SDD.vtree(manager)

println("saving vtree & sdd ...")
SDD.dot("$(@__DIR__)/output/after-rotate-vtree.dot", manager_vtree)
SDD.dot("$(@__DIR__)/output/after-rotate-sdd.dot", α)


text:
"""
Type stable methods for small tuples
"""
module TupleTools

using Base: tuple_type_head, tuple_type_tail, tuple_type_cons, tail, front, setindex
import Base: getindex, +, -, *, /
using StaticArrays: SVector, MVector

"""
    struct StaticLength{N} end

Like `Val{N}`, `StaticLength` can be used to construct a tuple of inferrable length
using `ntuple(f, StaticLength(N))`. Here, `StaticLength(N)` creates `StaticLength{N}()`
using a `Base.@pure` constructor. Furthermore, one can add and subtract `StaticLength`
objects, such that
```
StaticLength(N₁) + StaticLength(N₂) == StaticLength(N₁+N₂)
```
and
```
StaticLength(N₁) - StaticLength(N₂) == StaticLength(max(0, N₁-N₂))
```
"""
struct StaticLength{N}
end
Base.@pure StaticLength(N::Int) = StaticLength{N}()
Base.@pure Base.:+(::StaticLength{N₁}, ::StaticLength{N₂}) where {N₁,N₂} = StaticLength(N₁+N₂)
Base.@pure Base.:-(::StaticLength{N₁}, ::StaticLength{N₂}) where {N₁,N₂} = StaticLength(max(0,N₁-N₂))

@inline Base.ntuple(f, ::StaticLength{N}) where {N} = ntuple(f, Val{N}())

"""
	ntuple(f, n, T)

Create an `NTuple{n,T}` with values `(f(1), ..., f(n))`. For n>10, this is much faster
than `ntuple(f, n)`.  An `InexactError` will occur if any of `f(1),...,f(n)` cannot be
converted to type `T`.
"""
function Base.ntuple(f, n::Integer, ::Type{T}) where {T}
	t = n == 0  ? () :
		n == 1  ? NTuple{n,T}((f(1),)) :
		n == 2  ? NTuple{n,T}((f(1), f(2))) :
		n == 3  ? NTuple{n,T}((f(1), f(2), f(3))) :
		n == 4  ? NTuple{n,T}((f(1), f(2), f(3), f(4))) :
		n == 5  ? NTuple{n,T}((f(1), f(2), f(3), f(4), f(5))) :
		n == 6  ? NTuple{n,T}((f(1), f(2), f(3), f(4), f(5), f(6))) :
		n == 7  ? NTuple{n,T}((f(1), f(2), f(3), f(4), f(5), f(6), f(7))) :
		n == 8  ? NTuple{n,T}((f(1), f(2), f(3), f(4), f(5), f(6), f(7), f(8))) :
		n == 9  ? NTuple{n,T}((f(1), f(2), f(3), f(4), f(5), f(6), f(7), f(8), f(9))) :
		n == 10 ? NTuple{n,T}((f(1), f(2), f(3), f(4), f(5), f(6), f(7), f(8), f(9), f(10))) :
		begin
			v = MVector{n, T}(undef)
			@inbounds for i in 1:n
				v[i] = f(i)
			end
			Tuple(v)
		end
end

"""
	oneto(n::Integer)

Construct the tuple `(1,2,...,n)`.
"""
oneto(n::T) where {T<:Integer} = ntuple(identity, n, T)
# When N can be inferred, the following are compile-time generated.
oneto(::Val{0}) = ()
oneto(::Val{N}) where {N} = (oneto(Val(N-1))..., N)

@inline argtail2(a, b, c...) = c

"""
    tail2(t::Tuple) -> ::Tuple

Returns a tuple with the first two elements stripped, equivalent to `tail(tail(t))`
"""
@inline tail2(t::Tuple{Any,Any,Vararg{Any}}) = argtail2(t...)

"""
    unsafe_tail(t::Tuple) -> ::Tuple

Returns a tuple with the first element stripped, similar to `tail(t)`, but does
not error on an empty tuple (instead returning an empty tuple again). An empty tuple
is thus the fixed point of this function.
"""
unsafe_tail(t::Tuple{}) = t
unsafe_tail(t::Tuple) = tail(t)

"""
    unsafe_front(t::Tuple) -> ::Tuple

Returns a tuple with the last element stripped, similar to `front(t)`, but does
not error on an empty tuple (instead returning an empty tuple again). An empty tuple
is thus the fixed point of this function.
"""
unsafe_front(t::Tuple{}) = t
unsafe_front(t::Tuple) = front(t)

"""
    vcat(args...) -> ::Tuple

Like `vcat` for tuples, concatenates a combination of tuple arguments and non-tuple
arguments into a single tuple. Only works one level deep, i.e. tuples in tuples are
not expanded.
"""
vcat(t::Tuple) = t
vcat() = ()
vcat(t) = (t,)
vcat(a, args...) = (vcat(a)..., vcat(args...)...)

"""
    flatten(args...) -> ::Tuple

Flatten one or more tuples into a single tuple, such that every element of that tuple is itself not a tuple, otherwise it would also be expanded (i.e. flattened).
"""
flatten(x::Any) = (x,)
flatten(t::Tuple{}) = ()
flatten(t::Tuple) = (flatten(t[1])..., flatten(tail(t))...)
flatten(x, r...) = (flatten(x)..., flatten(r)...)

"""
    deleteat(t::Tuple, i::Int) -> ::Tuple
    deleteat(t::Tuple, I::Tuple{Vararg{Int}}) -> ::Tuple

Delete the element at location `i` in `t`; if a list `I` of indices is specified
(again as a tuple), the elements of these different positions are deleted.
"""
deleteat(t::Tuple, I::Tuple{Int}) = deleteat(t, I[1])
function deleteat(t::Tuple, I::Tuple{Int, Int, Vararg{Int}})
    any(i->!(1 <= i <= length(t)), I) && throw(BoundsError(t, I))
    _deleteat(t, sort(I, rev = true))
end
deleteat(t::Tuple, i::Int) =
    1 <= i <= length(t) ? _deleteat(t, i) : throw(BoundsError(t, i))
@inline _deleteat(t::Tuple, i::Int) = i == 1 ? tail(t) : (t[1], _deleteat(tail(t), i-1)...)

@inline _deleteat(t::Tuple, I::Tuple{Int}) = _deleteat(t, I[1])
@inline _deleteat(t::Tuple, I::Tuple{Int,Int,Vararg{Int}}) =
    _deleteat(_deleteat(t, I[1]), tail(I)) # assumes sorted from big to small

"""
    insertat(t::Tuple, i::Int, t2::Tuple) -> ::Tuple

Insert the elements of tuple `t2` at location `i` in `t`, i.e. the output tuple will
look as (t[1:i-1]..., t2..., t[i+1:end]). Note that element `t[i]` is deleted. Use
`setindex` for setting a single value at position `i`, or `insertafter(t, i, t2)` to
insert the contents of `t2` in between element `i` and `i+1` in `t`.
"""
insertat(t::Tuple, i::Int, t2::Tuple) =
    1 <= i <= length(t) ? _insertat(t, i, t2) : throw(BoundsError(t, i))
@inline _insertat(t::Tuple, i::Int, t2::Tuple) =
    i == 1 ? (t2..., tail(t)...) : (t[1], _insertat(tail(t), i-1, t2)...)

"""
    insertafter(t::Tuple, i::Int, t2::Tuple) -> ::Tuple

Insert the elements of tuple `t2` after location `i` in `t`, i.e. the output tuple will
look as (t[1:i]..., t2..., t[i+1:end]). Use index `i=0` or just `(t2..., t...)` to insert
`t2` in front of `t`; also see `insertat` to overwrite the element at position `i`.
"""
insertafter(t::Tuple, i::Int, t2::Tuple) =
    0 <= i <= length(t) ? _insertafter(t, i, t2) : throw(BoundsError(t, i))
@inline _insertafter(t::Tuple, i::Int, t2::Tuple) =
    i == 0 ? (t2..., t...) : (t[1], _insertafter(tail(t), i-1, t2)...)

"""
    sum(t::Tuple)

Returns the sum of the element of a tuple, or `0` for an empty tuple.
"""
sum(t::Tuple{}) = 0
sum(t::Tuple{Any}) = t[1]
sum(t::Tuple) = t[1]+sum(tail(t))

"""
    cumsum(t::Tuple)

Returns the cumulative sum of the elements of a tuple, or `()` for an empty tuple.
"""
function cumsum(t::Tuple)
    t_1, t_tail = first(t), tail(t)
    return (t_1, cumsum((t_1 + first(t_tail), tail(t_tail)...))...)
end
cumsum(t::Tuple{Any}) = t
cumsum(t::Tuple{}) = t

"""
    prod(t::Tuple)

Returns the product of the elements of a tuple, or `1` for an empty tuple.
"""
prod(t::Tuple{}) = 1
prod(t::Tuple{Any}) = t[1]
prod(t::Tuple) = t[1]*prod(tail(t))

"""
    cumprod(t::Tuple)

Returns the cumulative product of the elements of a tuple, or `()` for an empty tuple.
"""
function cumprod(t::Tuple)
    t_1, t_tail = first(t), tail(t)
    return (t_1, cumprod((t_1*first(t_tail), tail(t_tail)...))...)
end
cumprod(t::Tuple{Any}) = t
cumprod(t::Tuple{}) = t

"""
    minimum(t::Tuple)

Returns the smallest element of a tuple
"""
minimum(t::Tuple{Any}) = t[1]
minimum(t::Tuple) = min(t[1], minimum(tail(t)))

"""
    maximum(t::Tuple)

Returns the largest element of a tuple
"""
maximum(t::Tuple{Any}) = t[1]
maximum(t::Tuple) = max(t[1], maximum(tail(t)))

"""
    argmin(t::Tuple)

Returns the index of the minimum element in a tuple. If there are multiple
minimal elements, then the first one will be returned.
"""
argmin(t::Tuple) = findmin(t)[2]

"""
    argmax(t::Tuple)

Returns the index of the maximum element in a tuple. If there are multiple
minimal elements, then the first one will be returned.
"""
argmax(t::Tuple) = findmax(t)[2]

"""
    findmin(t::Tuple)

Returns the value and index of the minimum element in a tuple. If there are multiple
minimal elements, then the first one will be returned.
"""
findmin(t::Tuple) = Base.findmin(t)

"""
    findmax(t::Tuple)

Returns the value and index of the maximum element in a tuple. If there are multiple
maximal elements, then the first one will be returned.
"""
findmax(t::Tuple) = Base.findmax(t)

"""
    sort(t::Tuple; lt=isless, by=identity, rev::Bool=false) -> ::Tuple

Sorts the tuple `t`.
"""
#sort(t::Tuple; lt=isless, by=identity, rev::Bool=false) = _sort(t, lt, by, rev)
sort(t::Tuple; lt=isless, by=identity, rev::Bool=false) = length(t)<=10 ? _sort(t, lt, by, rev) : _sort_long(t, lt, by, rev)

@inline function _sort(t::Tuple, lt=isless, by=identity, rev::Bool=false)
    t1, t2 = _split(t)
    t1s = _sort(t1, lt, by, rev)
    t2s = _sort(t2, lt, by, rev)
    return _merge(t1s, t2s, lt, by, rev)
end
_sort(t::Tuple{Any}, lt=isless, by=identity, rev::Bool=false) = t
_sort(t::Tuple{}, lt=isless, by=identity, rev::Bool=false) = t

function _split(t::NTuple{N}) where N
    M = N>>1
    ntuple(i->t[i], M), ntuple(i->t[i+M], N-M)
end

function _merge(t1::Tuple, t2::Tuple, lt, by, rev)
    if lt(by(first(t1)), by(first(t2))) != rev
        return (first(t1), _merge(tail(t1), t2, lt, by, rev)...)
    else
        return (first(t2), _merge(t1, tail(t2), lt, by, rev)...)
    end
end
_merge(::Tuple{}, t2::Tuple, lt, by, rev) = t2
_merge(t1::Tuple, ::Tuple{}, lt, by, rev) = t1
_merge(::Tuple{}, ::Tuple{}, lt, by, rev) = ()

function _sort_long(t, lt=isless, by=identity, rev::Bool=false)
	s = MVector(t)
	sort!(s; lt =lt, by = by, rev = rev, alg = QuickSort)
	Tuple(s)
end

"""
    sortperm(t::Tuple; lt=isless, by=identity, rev::Bool=false) -> ::Tuple

Computes a tuple that contains the permutation required to sort `t`.
"""
#sortperm(t::Tuple; lt=isless, by=identity, rev::Bool=false) = _sortperm(t, lt, by, rev)
sortperm(t::Tuple; lt=isless, by=identity, rev::Bool=false) = length(t)<=10 ? _sortperm(t, lt, by, rev) : _sortperm_long(t, lt, by, rev)

function _sortperm(t::Tuple, lt=isless, by=identity, rev::Bool=false)
    map(first, _sort(ntuple(n->(n,by(t[n])), length(t)), lt, last, rev))
end

function _sortperm_long(t, lt=isless, by=identity, rev::Bool=false)
	s = MVector(t)
	sortperm!(s, SVector(t); lt =lt, by = by, rev = rev, alg = QuickSort)
	Tuple(s)
end

"""
    getindices(t::Tuple, I::Tuple{Vararg{Int}}) -> ::Tuple

Get the indices `t[i] for i in I`, again as tuple.
"""
getindices(t::Tuple, ind::Tuple{Vararg{Int}}) =
    (t[ind[1]], getindices(t, tail(ind))...)
getindices(t::Tuple, ind::Tuple{}) = ()

getindex(t::Tuple, b::Tuple{Vararg{Bool}}) = length(b) == length(t) ? getindex(t, findall(b)) : throw(BoundsError(t, b))

"""
    permute(t::Tuple, p) -> ::Tuple

Permute the elements of tuple `t` according to the permutation in `p`.
"""
permute(t::NTuple{N,Any}, p::NTuple{N,Int}) where {N} =
    isperm(p) ? _permute(t, p) : throw(ArgumentError("not a valid permutation: $p"))
permute(t::NTuple{N,Any}, p) where {N} =
    isperm(p) && length(p) == N ? _permute(t, p) :
        throw(ArgumentError("not a valid permutation: $p"))

_permute(t::NTuple{N,Any}, p::NTuple{N,Int}) where {N} = getindices(t, p)
_permute(t::NTuple{N,Any}, p) where {N} = ntuple(n->t[p[n]], StaticLength(N))

"""
    isperm(p) -> ::Bool

A non-allocating alternative to Base.isperm(p) that is much faster for small permutations.
"""
function isperm(p)
    N = length(p)
    @inbounds for i = 1:N
        1 <= p[i] <= N || return false
        for j = i+1:N
            p[i] == p[j] && return false
        end
    end
    return true
end

"""
    invperm(p::NTuple{N,Int}) -> ::NTuple{N,Int}

Inverse permutation of a permutation `p`.
"""
invperm(p::Tuple{Vararg{Int}}) = length(p)<=6 ? _sortperm(p) : _invperm_long(p)

# Note:  Whether t is a permutation is not checked.
function _invperm_long(t::Tuple{Vararg{<:Number}})
	p = MVector{length(t), Int}(undef)
	for i = 1:length(t)
		p[t[i]] = i
	end
	Tuple(p)
end

"""
	invpermute(t::Tuple, p::Tuple)

Returns `s` such that `s[p] = t`.  `p` is not checked.
"""
function invpermute(t::NTuple{N}, p::NTuple{N,<:Number}) where {N}
	s = MVector{length(t), Int}(undef)
	for i = 1:length(t)
		s[p[i]] = t[i]
	end
	Tuple(s)
end

"""
    diff(v::Tuple) -> Tuple

Finite difference operator of tuple `v`.
"""
diff(v::Tuple{}) = () # similar to diff([])
diff(v::Tuple{Any}) = ()
diff(v::Tuple) = (v[2]-v[1], diff(Base.tail(v))...)

end # module


text:
using Cpp

# Find the test library (this is not normally needed)
fnames = ["libdemo.so", "libdemo.dylib", "libdemo.dll"]
paths = [pwd(), joinpath(Pkg.dir(), "Cpp", "deps")]
global libname
found = false
for path in paths
    if !found
        for fname in fnames
            libname = Base.find_in_path(joinpath(path, fname))
            if isfile(libname)
                found = true
                break
            end
        end
    end
end
if !isfile(libname)
    error("Library cannot be found; it may not have been built correctly.\n Try include(\"build.jl\") from within the deps directory.")
end
const libdemo = libname

x = 3.5
x2 = @cpp ccall((:timestwo, libdemo), Float64, (Float64,), x)
@assert x2 == 2x
y = 3
y2 = @cpp ccall((:timestwo, libdemo), Cint, (Cint,), y)
@assert 2y == y2

# Test substitution mangling
val = @cpp ccall((:testsub, libdemo), Cint, (Ptr{Cint}, Ptr{Float64}, Ptr{Float64}), [int32(5)], [7.8], [-1.4])
@assert val == 1