text:
@testset "singular case" begin
# From Nocedal & Wright, p. 288-289

# Jacobian is singular at the starting point.
# Used to test the behavior of algorithms in that context.

function f_sinj!(F, x)
    F[1] = x[1]
    F[2] = 10*x[1]/(x[1]+convert(eltype(x), 0.1))+2*x[2]^2
end

function g_sinj!(J, x)
    J[1, 1] = 1
    J[1, 2] = 0
    J[2, 1] = 1/(x[1]+convert(eltype(x), 0.1))^2
    J[2, 2] = 4*x[2]
end

df = OnceDifferentiable(f_sinj!, g_sinj!, [3.0, 0.0], [3.0, 0.0])
df32 = OnceDifferentiable(f_sinj!, g_sinj!, [3.0f0, 0.0f0], [3.0f0, 0.0f0])

# Test disabled, not stable across runs
#r = nlsolve(df, [ 3.0; 0], method = :newton, ftol = 1e-5)
#@assert converged(r)
#@assert norm(r.zero) < 1e-5

r = nlsolve(df, [ 3.0; 0.0], method = :trust_region)
@assert converged(r)
@assert norm(r.zero) < 1e-6

r = nlsolve(df32, [3.0f0; 0.0f0], method = :trust_region)
@assert converged(r)
@assert norm(r.zero) < 1e-6

let a = rand(10)
    A = a*a'
    global f_let!, g_let!
    function f_let!(fvec, x)
        copyto!(fvec, A*x)
    end

    function g_let!(fjac, x)
        copyto!(fjac, A)
    end
end

df = OnceDifferentiable(f_let!, g_let!, rand(10), rand(10))
r = nlsolve(df, rand(10), method = :trust_region)
end


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

name = "Libassuan"
version = v"2.5.5"

# Collection of sources required to build libgcrypt
sources = [
    ArchiveSource("https://gnupg.org/ftp/gcrypt/libassuan/libassuan-$(version).tar.bz2",
                  "8e8c2fcc982f9ca67dcbb1d95e2dc746b1739a4668bc20b3a3c5be632edb34e4"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms

# Tried -no-undefined but still couldn't build for windows
script = raw"""
cd $WORKSPACE/srcdir/libassuan-*/

if [[ "${target}" == x86_64-*-mingw* ]]; then
    # `gpgrt-config` for this platform returns garbage results.  We replace it with
    # a simple wrapper around `pkg-config`, so that we can easily build the shared library.
    FLAGS=(GPG_ERROR_CONFIG="../gpgrt-config.sh" ac_cv_path_GPGRT_CONFIG="../gpgrt-config.sh")
fi
export CPPFLAGS="-I${includedir}"
./configure --prefix=${prefix} --host=${target} --build=${MACHTYPE} "${FLAGS[@]}"
make -j${nproc}
make install

install_license ${WORKSPACE}/srcdir/libassuan-*/COPYING.LIB
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.  We are manually disabling
# many platforms that do not seem to work.
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libassuan", "libassuan6"], :libassuan),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # Future versions of `Libgpg_error_jll` maybe can have a more lax compat,
    # but the move 1.36 -> 1.42 changed the soname for FreeBSD and Windows.
    Dependency("Libgpg_error_jll", v"1.36.0"; compat="=1.36.0"),
]

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

ext: jl

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
__precompile__(false)

module PCLLite

const VERBOSE = Bool(parse(Int, get(ENV, "PCLLITEJL_VERBOSE", "0")))

const BOOST_INCLUDE_DIR = get(ENV, "BOOST_INCLUDE_DIR", "/usr/local/include")
const FLANN_INCLUDE_DIR = get(ENV, "FLANN_INCLUDE_DIR", "/usr/local/include")
const EIGEN_INCLUDE_DIR = get(ENV, "EIGEN_INCLUDE_DIR",
    "/usr/local/include/eigen3")

searchdir(path, key) = filter(x->contains(x, key), readdir(path))

const VTK_INCLUDE_DIR_PARENT = get(ENV, "VTK_INCLUDE_DIR_PARENT",
    "/usr/local/include")

# Search GUI backend
vtk_dirs = searchdir(VTK_INCLUDE_DIR_PARENT, "vtk-")
const VTK_INCLUDE_DIR = get(ENV, "VTK_INCLUDE_DIR",
    isempty(vtk_dirs) ? "" : joinpath(VTK_INCLUDE_DIR_PARENT, vtk_dirs[1]))

const has_vtk_backend = isfile(joinpath(VTK_INCLUDE_DIR, "vtkVersion.h"))
if has_vtk_backend
    VERBOSE && info("vtk include directory found: $VTK_INCLUDE_DIR")
end

using BinDeps

# Load dependency
deps = joinpath(Pkg.dir("PCLLite"), "deps", "deps.jl")
if isfile(deps)
    include(deps)
else
    error("PCLLite not properly installed. Please run Pkg.build(\"PCLLite\")")
end

VERBOSE && info("Loading Cxx.jl...")
using Cxx
using CxxStd

macro timevb(expr)
    if VERBOSE
        return Expr(:macrocall, Symbol("@time"), expr)
    else
        return expr
    end
end

VERBOSE && info("dlopen...")
for lib in [
        libpcl_common,
        libpcl_io,
        libpcl_features,
        libpcl_filters,
        libpcl_kdtree,
        libpcl_keypoints,
        libpcl_segmentation,
        libpcl_visualization,
        libpcl_recognition,
        libpcl_registration,
        libpcl_octree,
        libpcl_surface,
        libpcl_people,
        libpcl_search,
        libpcl_tracking,
        ]
    p = Libdl.dlopen_e(lib, Libdl.RTLD_GLOBAL)
    p == C_NULL && warn("Failed to load: $lib")
end

# Make sure vtk libraries are loaded before calling @cxx vtkVersion::xxx()
if has_vtk_backend
    cxxinclude(joinpath(VTK_INCLUDE_DIR, "vtkVersion.h"))
    global const vtk_version = bytestring(@cxx vtkVersion::GetVTKVersion())
    VERBOSE && info("vtk version: $vtk_version")
end

const system_include_top = "/usr/local/include"
const local_include_top = joinpath(Pkg.dir("PCL", "deps", "usr", "include"))

function get_pcl_version(top)
    dirs = searchdir(top, "pcl-")
    isempty(dirs) && error("could not find pcl directory")
    pcl_dir = dirs[1]
    return pcl_dir[5:end]
end

topdir_to_be_included = local_include_top

if isdir(local_include_top)
    VERBOSE && info("Including headers from local path: $local_include_top")
    pcl_version = get_pcl_version(local_include_top)
elseif !isempty(searchdir(joinpath(system_include_top), "pcl-"))
    VERBOSE && info("Including headers from system path: $system_include_top")
    pcl_version = get_pcl_version(system_include_top)
    topdir_to_be_included = system_include_top
else
    error("Cannot find PCL headers")
end

function add_header_dirs(top)
    # Boost (required)
    addHeaderDir(BOOST_INCLUDE_DIR, kind=C_System)

    # Eigen (required)
    addHeaderDir(EIGEN_INCLUDE_DIR, kind=C_System)

    # FLANN (required)
    addHeaderDir(FLANN_INCLUDE_DIR, kind=C_System)

    # PCL top directory
    addHeaderDir(top, kind=C_System)
    addHeaderDir(joinpath(top, "pcl"), kind=C_System)
end

function include_headers(top)
    # This is necesarry, but not sure why...
    cxxinclude("iostream")

    # top level
    VERBOSE && info("Include pcl top-level headers")
    @timevb for name in ["pcl_base.h", "point_cloud.h", "point_types.h"]
        cxxinclude(joinpath(top, "pcl", name))
    end

    # boost
    cxxinclude("boost/version.hpp")
end

VERBOSE && info("pcl_version: $pcl_version")
add_header_dirs(joinpath(topdir_to_be_included, "pcl-$pcl_version"))
include_headers(joinpath(topdir_to_be_included, "pcl-$pcl_version"))

# Check boost version
const _BOOST_VERSION = icxx"BOOST_VERSION;"
const BOOST_VERSION_MAJOR = trunc(Int, _BOOST_VERSION / 100000)
const BOOST_VERSION_MINOR = trunc(Int, _BOOST_VERSION / 100 % 1000)
const BOOST_VERSION_PATCH = trunc(Int, _BOOST_VERSION % 100)
const BOOST_VERSION = join([BOOST_VERSION_MAJOR, BOOST_VERSION_MINOR, BOOST_VERSION_PATCH], ".")
VERBOSE && info("boost version: $BOOST_VERSION")

# Check FLANN vesion
# make sure FLANN_INCLUDE_DIR is addded as kind=C_System
cxxinclude(joinpath(FLANN_INCLUDE_DIR, "flann/flann.h"))
cxx"""
namespace pcllite {
std::string getFLANNVersion() { return FLANN_VERSION_; }
}
"""
getFLANNVersion() = bytestring(@cxx pcllite::getFLANNVersion())
VERBOSE && info("FLANN version: $(getFLANNVersion())")

end # module PCLLite


text:
using Printf
using BenchmarkTools
using InteractiveUtils

"""
    compute_annihilation(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix, spin_case)

Compute an creation operator between states `bra_v` and `ket_v`
# Arguments
- `no`: number of orbitals
- `bra_na`: number of alpha electrons in bra
- `bra_nb`: number of beta electrons in bra
- `ket_na`: number of alpha electrons in ket
- `ket_nb`: number of beta electrons in ket
- `bra_v`: basis vectors in bra
- `ket_v`: basis vectors ket
- `spin_case`: either (alpha|beta)
"""
function compute_annihilation(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix, spin_case)
#={{{=#
    bra_a = DeterminantString(no,bra_na)
    bra_b = DeterminantString(no,bra_nb)

    ket_a = DeterminantString(no,ket_na)
    ket_b = DeterminantString(no,ket_nb)

    if spin_case == "alpha"
        bra_na == ket_na+1 || throw(DimensionMismatch)
        bra_nb == ket_nb || throw(DimensionMismatch)
    elseif spin_case == "beta"
        bra_nb == ket_nb+1 || throw(DimensionMismatch)
        bra_na == ket_na || throw(DimensionMismatch)
    end

    ket = ket_a
    bra = bra_a
    if spin_case == "beta"
        ket = ket_b
        bra = bra_b
    end
    virt = zeros(typeof(ket.no), ket.no-ket.ne)

    #
    # precompute lookup table for binomial coefficients
    _binomial = Array{Int,2}(undef,no+1,no+1)
    for i in 0:no
        for j in i:no
            _binomial[j+1,i+1] = binomial(j,i)
        end
    end

    bra_M = size(bra_v,2)
    ket_M = size(ket_v,2)

    #println((no, bra_na, bra_nb, bra_a.max, bra_b.max))
    #println(size(bra_v), bra_a.max, " ", bra_b.max, " ", bra_M)

    v1 = reshape(bra_v, bra_a.max, bra_b.max, bra_M)
    v2 = reshape(ket_v, ket_a.max, ket_b.max, ket_M)
    if spin_case == "alpha"
        v1 = permutedims(v1, [2,3,1])
        v2 = permutedims(v2, [2,3,1])
    elseif spin_case == "beta"
        v1 = permutedims(v1, [1,3,2])
        v2 = permutedims(v2, [1,3,2])
    else
        throw(Exception)
    end

    #
    #   TDM[p,s,t] =
    tdm = zeros(Float64, bra_M, ket_M, no)
    reset!(ket)
    reset!(bra)
    for J in 1:ket.max
        get_unoccupied!(virt, ket)
        for a in virt
            copy!(bra,ket)
            apply_creation!(bra,a)
            issorted(bra.config) || throw(Exception)
            calc_linear_index!(bra,_binomial)
            I = bra.lin_index
            @views bra_vI = v1[:,:,I]
            @views ket_vJ = v2[:,:,J]
            @views tdma = tdm[:,:,a]
            sgn = bra.sign
            if spin_case == "beta"
                if ket_a.ne%2 != 0
                    sgn = -sgn
                end
            end

            if sgn == 1
                @tensor begin
                    tdma[s,t] += bra_vI[I,s] * ket_vJ[I,t]
                end
                #tdma += bra_vI' * ket_vJ
            elseif sgn == -1
                #tdma -= bra_vI' * ket_vJ
                @tensor begin
                    tdma[s,t] -= bra_vI[I,s] * ket_vJ[I,t]
                end
            else
                throw(Exception)
            end
        end
        incr!(ket)
    end
    tdm = permutedims(tdm, [3,1,2])
    #println(size(tdm))
    return tdm
#=}}}=#
end

"""
    compute_AA(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix, spin_case)

Compute represntation of 2creation operators between states `bra_v` and `ket_v`
# Arguments
- `no`: number of orbitals
- `bra_na`: number of alpha electrons in bra
- `bra_nb`: number of beta electrons in bra
- `ket_na`: number of alpha electrons in ket
- `ket_nb`: number of beta electrons in ket
- `bra_v`: basis vectors in bra
- `ket_v`: basis vectors ket
- `spin_case`: either (alpha|beta)
"""
function compute_AA(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix, spin_case)
#={{{=#
    bra_a = DeterminantString(no,bra_na)
    bra_b = DeterminantString(no,bra_nb)

    ket_a = DeterminantString(no,ket_na)
    ket_b = DeterminantString(no,ket_nb)

    if spin_case == "alpha"
        bra_na == ket_na+2 || throw(DimensionMismatch)
        bra_nb == ket_nb || throw(DimensionMismatch)
    elseif spin_case == "beta"
        bra_nb == ket_nb+2 || throw(DimensionMismatch)
        bra_na == ket_na || throw(DimensionMismatch)
    end

    ket = ket_a
    bra = bra_a
    if spin_case == "beta"
        ket = ket_b
        bra = bra_b
    end
    virt = zeros(typeof(ket.no), ket.no-ket.ne)

    #
    # precompute lookup table for binomial coefficients
    _binomial = Array{Int,2}(undef,no+1,no+1)
    for i in 0:no
        for j in i:no
            _binomial[j+1,i+1] = binomial(j,i)
        end
    end

    bra_M = size(bra_v,2)
    ket_M = size(ket_v,2)

    #println((no, bra_na, bra_nb, bra_a.max, bra_b.max))
    #println(size(bra_v), bra_a.max, " ", bra_b.max, " ", bra_M)

    v1 = reshape(bra_v, bra_a.max, bra_b.max, bra_M)
    v2 = reshape(ket_v, ket_a.max, ket_b.max, ket_M)
    if spin_case == "alpha"
        v1 = permutedims(v1, [2,3,1])
        v2 = permutedims(v2, [2,3,1])
    elseif spin_case == "beta"
        v1 = permutedims(v1, [1,3,2])
        v2 = permutedims(v2, [1,3,2])
    else
        throw(Exception)
    end

    #
    #   TDM[p,s,t] =
    tdm = zeros(Float64, bra_M, ket_M, no, no)
    reset!(ket)
    reset!(bra)
    for J in 1:ket.max
        get_unoccupied!(virt, ket)
        for a in virt
            for b in virt
                a != b || continue
                copy!(bra,ket)
                apply_creation!(bra,b)
                apply_creation!(bra,a)
                calc_linear_index!(bra,_binomial)
                I = bra.lin_index
                @views bra_vI = v1[:,:,I]
                @views ket_vJ = v2[:,:,J]
                @views tdma = tdm[:,:,a,b]
                sgn = bra.sign

                if sgn == 1
                    @tensor begin
                        tdma[s,t] += bra_vI[I,s] * ket_vJ[I,t]
                    end
                elseif sgn == -1
                    @tensor begin
                        tdma[s,t] -= bra_vI[I,s] * ket_vJ[I,t]
                    end
                else
                    throw(Exception)
                end
            end
        end
        incr!(ket)
    end
    tdm = permutedims(tdm, [3,4,1,2])
    return tdm
#=}}}=#
end

"""
    compute_Aa(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix, spin_case)

Compute representation of a'a operators between states `bra_v` and `ket_v`
# Arguments
- `no`: number of orbitals
- `bra_na`: number of alpha electrons in bra
- `bra_nb`: number of beta electrons in bra
- `ket_na`: number of alpha electrons in ket
- `ket_nb`: number of beta electrons in ket
- `bra_v`: basis vectors in bra
- `ket_v`: basis vectors ket
- `spin_case`: either (alpha|beta)
"""
function compute_Aa(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix, spin_case)
#={{{=#
    bra_a = DeterminantString(no,bra_na)
    bra_b = DeterminantString(no,bra_nb)

    ket_a = DeterminantString(no,ket_na)
    ket_b = DeterminantString(no,ket_nb)

    if spin_case == "alpha"
        bra_na == ket_na || throw(DimensionMismatch)
        bra_nb == ket_nb || throw(DimensionMismatch)
    elseif spin_case == "beta"
        bra_nb == ket_nb || throw(DimensionMismatch)
        bra_na == ket_na || throw(DimensionMismatch)
    end

    ket = ket_a
    bra = bra_a
    if spin_case == "beta"
        ket = ket_b
        bra = bra_b
    end
    virt = zeros(typeof(ket.no), ket.no-ket.ne)

    #
    # precompute lookup table for binomial coefficients
    _binomial = Array{Int,2}(undef,no+1,no+1)
    for i in 0:no
        for j in i:no
            _binomial[j+1,i+1] = binomial(j,i)
        end
    end

    bra_M = size(bra_v,2)
    ket_M = size(ket_v,2)

    v1 = reshape(bra_v, bra_a.max, bra_b.max, bra_M)
    v2 = reshape(ket_v, ket_a.max, ket_b.max, ket_M)
    if spin_case == "alpha"
        v1 = permutedims(v1, [2,3,1])
        v2 = permutedims(v2, [2,3,1])
    elseif spin_case == "beta"
        v1 = permutedims(v1, [1,3,2])
        v2 = permutedims(v2, [1,3,2])
    else
        throw(Exception)
    end

    #
    #   TDM[p,s,t] =
    tdm = zeros(Float64, bra_M, ket_M, no, no)
    reset!(ket)
    reset!(bra)
    for J in 1:ket.max
        get_unoccupied!(virt, ket)
        for q in 1:no
            for p in 1:no
                copy!(bra,ket)
                apply_annihilation!(bra,q)
                apply_creation!(bra,p)
                bra.sign != 0 || continue
                calc_linear_index!(bra,_binomial)
                I = bra.lin_index
                @views bra_vI = v1[:,:,I]
                @views ket_vJ = v2[:,:,J]
                @views tdma = tdm[:,:,p,q]
                sgn = bra.sign

                if sgn == 1
                    @tensor begin
                        tdma[s,t] += bra_vI[I,s] * ket_vJ[I,t]
                    end
                elseif sgn == -1
                    @tensor begin
                        tdma[s,t] -= bra_vI[I,s] * ket_vJ[I,t]
                    end
                else
                    throw(DomainError(sgn))
                end
            end
        end
        incr!(ket)
    end
    tdm = permutedims(tdm, [3,4,1,2])
    return tdm
#=}}}=#
end

"""
    compute_Ab(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix)

Compute representation of p'q operators between states `bra_v` and `ket_v`, where
`p` is alpha and `q` is beta.
# Arguments
- `no`: number of orbitals
- `bra_na`: number of alpha electrons in bra
- `bra_nb`: number of beta electrons in bra
- `ket_na`: number of alpha electrons in ket
- `ket_nb`: number of beta electrons in ket
- `bra_v`: basis vectors in bra
- `ket_v`: basis vectors ket

    G(pqst) = v(IJs) <IJ|p'q|KL> v(KLt)
            = v(IJs) <J|<I|p'q|K>|L> v(KLt)
            = v(IJs) <J|<I|p'|K>q|L> v(KLt) (-1)^N(K)
            = v(IJs) <I|p'|K> <J|q|L> v(KLt) (-1)^N(K)

            this can be vectorized if needed. for now:
            loop over p,q
                loop over K,L
                    sign = -1^N(K)
                    find indices of I,J
                    G(pq,st) += v(IJ,s) v(KL,t) * sign
                end
            end
"""
function compute_Ab(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix)
#={{{=#
    bra_a = DeterminantString(no,bra_na)
    bra_b = DeterminantString(no,bra_nb)

    ket_a = DeterminantString(no,ket_na)
    ket_b = DeterminantString(no,ket_nb)

    bra_na == ket_na+1 || throw(DimensionMismatch)
    bra_nb == ket_nb-1 || throw(DimensionMismatch)

    #
    # precompute lookup table for binomial coefficients
    _binomial = Array{Int,2}(undef,no+1,no+1)
    for i in 0:no
        for j in i:no
            _binomial[j+1,i+1] = binomial(j,i)
        end
    end

    bra_M = size(bra_v,2)
    ket_M = size(ket_v,2)

    v1 = reshape(bra_v, bra_a.max, bra_b.max, bra_M)
    v2 = reshape(ket_v, ket_a.max, ket_b.max, ket_M)
    v1 = permutedims(v1, [3,1,2])
    v2 = permutedims(v2, [3,1,2])

    #
    #   TDM[p,s,t] =
    tdm = zeros(Float64, bra_M, ket_M, no, no)
    bra_a = deepcopy(ket_a)
    bra_b = deepcopy(ket_b)

    sgnK = 1
    if ket_a.ne % 2 != 0
        sgnK = -sgnK
    end

    reset!(ket_a)
    for K in 1:ket_a.max
        for p in 1:no
            copy!(bra_a,ket_a)
            apply_creation!(bra_a,p)
            bra_a.sign != 0 || continue
            calc_linear_index!(bra_a,_binomial)
            I = bra_a.lin_index

            reset!(ket_b)
            for L in 1:ket_b.max
                for q in 1:no
                    copy!(bra_b,ket_b)
                    apply_annihilation!(bra_b,q)
                    bra_b.sign != 0 || continue
                    calc_linear_index!(bra_b,_binomial)
                    J = bra_b.lin_index
                    @views tdm_pq = tdm[:,:,p,q]
                    @views v1_IJ = v1[:,I,J]
                    @views v2_KL = v2[:,K,L]
                    sgn = sgnK * bra_a.sign * bra_b.sign

                    if sgn == 1
                        @tensor begin
                            tdm_pq[s,t] += v1_IJ[s] * v2_KL[t]
                        end
                    elseif sgn == -1
                        @tensor begin
                            tdm_pq[s,t] -= v1_IJ[s] * v2_KL[t]
                        end
                    else
                        throw(Exception)
                    end
                end #q
                incr!(ket_b)
            end #L
        end #p
        incr!(ket_a)
    end #K
    tdm = permutedims(tdm, [3,4,1,2])
    return tdm
#=}}}=#
end

"""
    compute_AB(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix)

Compute representation of p'q' operators between states `bra_v` and `ket_v`, where
`p` is alpha and `q` is beta.
# Arguments
- `no`: number of orbitals
- `bra_na`: number of alpha electrons in bra
- `bra_nb`: number of beta electrons in bra
- `ket_na`: number of alpha electrons in ket
- `ket_nb`: number of beta electrons in ket
- `bra_v`: basis vectors in bra
- `ket_v`: basis vectors ket

    G(pqst) = v(IJs) <IJ|p'q'|KL> v(KLt)
            = v(IJs) <J|<I|p'q'|K>|L> v(KLt)
            = v(IJs) <J|<I|p'|K>q'|L> v(KLt) (-1)^N(K)
            = v(IJs) <I|p'|K> <J|q'|L> v(KLt) (-1)^N(K)

            this can be vectorized if needed. for now:
            loop over p,q
                loop over K,L
                    sign = -1^N(K)
                    find indices of I,J
                    G(pq,st) += v(IJ,s) v(KL,t) * sign
                end
            end
"""
function compute_AB(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix)
#={{{=#
    bra_a = DeterminantString(no,bra_na)
    bra_b = DeterminantString(no,bra_nb)

    ket_a = DeterminantString(no,ket_na)
    ket_b = DeterminantString(no,ket_nb)

    bra_na == ket_na+1 || throw(DimensionMismatch)
    bra_nb == ket_nb+1 || throw(DimensionMismatch)

    #
    # precompute lookup table for binomial coefficients
    _binomial = Array{Int,2}(undef,no+1,no+1)
    for i in 0:no
        for j in i:no
            _binomial[j+1,i+1] = binomial(j,i)
        end
    end

    bra_M = size(bra_v,2)
    ket_M = size(ket_v,2)

    v1 = reshape(bra_v, bra_a.max, bra_b.max, bra_M)
    v2 = reshape(ket_v, ket_a.max, ket_b.max, ket_M)
    v1 = permutedims(v1, [3,1,2])
    v2 = permutedims(v2, [3,1,2])

    #
    #   TDM[p,s,t] =
    tdm = zeros(Float64, bra_M, ket_M, no, no)
    bra_a = deepcopy(ket_a)
    bra_b = deepcopy(ket_b)

    sgnK = 1
    if ket_a.ne % 2 != 0
        sgnK = -sgnK
    end

    reset!(ket_a)
    for K in 1:ket_a.max
        for p in 1:no
            copy!(bra_a,ket_a)
            apply_creation!(bra_a,p)
            bra_a.sign != 0 || continue
            calc_linear_index!(bra_a,_binomial)
            I = bra_a.lin_index

            reset!(ket_b)
            for L in 1:ket_b.max
                for q in 1:no
                    copy!(bra_b,ket_b)
                    apply_creation!(bra_b,q)
                    bra_b.sign != 0 || continue
                    calc_linear_index!(bra_b,_binomial)
                    J = bra_b.lin_index
                    @views tdm_pq = tdm[:,:,p,q]
                    @views v1_IJ = v1[:,I,J]
                    @views v2_KL = v2[:,K,L]
                    sgn = sgnK * bra_a.sign * bra_b.sign

                    if sgn == 1
                        @tensor begin
                            tdm_pq[s,t] += v1_IJ[s] * v2_KL[t]
                        end
                    elseif sgn == -1
                        @tensor begin
                            tdm_pq[s,t] -= v1_IJ[s] * v2_KL[t]
                        end
                    else
                        throw(Exception)
                    end
                end #q
                incr!(ket_b)
            end #L
        end #p
        incr!(ket_a)
    end #K
    tdm = permutedims(tdm, [3,4,1,2])
    return tdm
#=}}}=#
end

"""
    compute_AAa(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix, spin_case)

Compute represntation of 2creation operators between states `bra_v` and `ket_v`
# Arguments
- `no`: number of orbitals
- `bra_na`: number of alpha electrons in bra
- `bra_nb`: number of beta electrons in bra
- `ket_na`: number of alpha electrons in ket
- `ket_nb`: number of beta electrons in ket
- `bra_v`: basis vectors in bra
- `ket_v`: basis vectors ket
- `spin_case`: either (alpha|beta)
"""
function compute_AAa(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix, spin_case)
#={{{=#
    bra_a = DeterminantString(no,bra_na)
    bra_b = DeterminantString(no,bra_nb)

    ket_a = DeterminantString(no,ket_na)
    ket_b = DeterminantString(no,ket_nb)

    if spin_case == "alpha"
        bra_na == ket_na+1 || throw(DimensionMismatch)
        bra_nb == ket_nb || throw(DimensionMismatch)
    elseif spin_case == "beta"
        bra_nb == ket_nb+1 || throw(DimensionMismatch)
        bra_na == ket_na || throw(DimensionMismatch)
    end

    ket = ket_a
    bra = bra_a
    if spin_case == "beta"
        ket = ket_b
        bra = bra_b
    end
    virt = zeros(typeof(ket.no), ket.no-ket.ne)

    #
    # precompute lookup table for binomial coefficients
    _binomial = Array{Int,2}(undef,no+1,no+1)
    for i in 0:no
        for j in i:no
            _binomial[j+1,i+1] = binomial(j,i)
        end
    end

    bra_M = size(bra_v,2)
    ket_M = size(ket_v,2)

    #println((no, bra_na, bra_nb, bra_a.max, bra_b.max))
    #println(size(bra_v), bra_a.max, " ", bra_b.max, " ", bra_M)

    v1 = reshape(bra_v, bra_a.max, bra_b.max, bra_M)
    v2 = reshape(ket_v, ket_a.max, ket_b.max, ket_M)
    if spin_case == "alpha"
        v1 = permutedims(v1, [2,3,1])
        v2 = permutedims(v2, [2,3,1])
    elseif spin_case == "beta"
        v1 = permutedims(v1, [1,3,2])
        v2 = permutedims(v2, [1,3,2])
    else
        throw(Exception)
    end

    #
    #   TDM[s,t,p,q,r] =
    tdm = zeros(Float64, bra_M, ket_M, no, no, no)
    reset!(ket)
    reset!(bra)
    for J in 1:ket.max
        get_unoccupied!(virt, ket)
        for p in 1:no
            for q in 1:no
                for r in 1:no
                copy!(bra,ket)
                apply_annihilation!(bra,r)
                apply_creation!(bra,q)
                apply_creation!(bra,p)
                bra.sign != 0 || continue
                calc_linear_index!(bra,_binomial)
                I = bra.lin_index
                @views bra_vI = v1[:,:,I]
                @views ket_vJ = v2[:,:,J]
                @views tdm_pqr = tdm[:,:,p,q,r]
                sgn = bra.sign
                if spin_case == "beta"
                    if ket_a.ne%2 != 0
                        sgn = -sgn
                    end
                end

                if sgn == 1
                    @tensor begin
                        tdm_pqr[s,t] += bra_vI[I,s] * ket_vJ[I,t]
                    end
                elseif sgn == -1
                    @tensor begin
                        tdm_pqr[s,t] -= bra_vI[I,s] * ket_vJ[I,t]
                    end
                else
                    throw(Exception)
                end
            end
            end
        end
        incr!(ket)
    end
    tdm = permutedims(tdm, [3,4,5,1,2])
    return tdm
#=}}}=#
end

"""
    compute_ABa(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix, spin_case)

Compute represntation of 2creation operators between states `bra_v` and `ket_v`
# Arguments
- `no`: number of orbitals
- `bra_na`: number of alpha electrons in bra
- `bra_nb`: number of beta electrons in bra
- `ket_na`: number of alpha electrons in ket
- `ket_nb`: number of beta electrons in ket
- `bra_v`: basis vectors in bra
- `ket_v`: basis vectors ket
- `spin_case`: either (alpha|beta) giving ABa or BAb
"""
function compute_ABa(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix)
#={{{=#
    bra_a = DeterminantString(no,bra_na)
    bra_b = DeterminantString(no,bra_nb)

    ket_a = DeterminantString(no,ket_na)
    ket_b = DeterminantString(no,ket_nb)

    bra_na == ket_na+0 || throw(DimensionMismatch)
    bra_nb == ket_nb+1 || throw(DimensionMismatch)

    #
    # precompute lookup table for binomial coefficients
    _binomial = Array{Int,2}(undef,no+1,no+1)
    for i in 0:no
        for j in i:no
            _binomial[j+1,i+1] = binomial(j,i)
        end
    end

    bra_M = size(bra_v,2)
    ket_M = size(ket_v,2)

    v1 = reshape(bra_v, bra_a.max, bra_b.max, bra_M)
    v2 = reshape(ket_v, ket_a.max, ket_b.max, ket_M)
    v1 = permutedims(v1, [3,1,2])
    v2 = permutedims(v2, [3,1,2])

    #
    #   TDM[pqr,s,t] =
    tdm = zeros(Float64, bra_M, ket_M, no, no, no)
    bra_a = deepcopy(ket_a)
    bra_b = deepcopy(ket_b)

    sgnK = -1 # start from -1 to account for ABa -> AaB
    if (ket_a.ne) % 2 != 0
        sgnK = -sgnK
    end

    # here, p'r are Alpha
    #  v(K,L,s) <KL|p'q'r|IJ> V(I,J,t)
    # -v(K,L,s) <L|<K|p'rq'|I>|J> V(I,J,t)
    # -v(K,L,s) <L|<K|p'r|I>q'|J> V(I,J,t) (-1)^N(I)
    # -v(K,L,s) <K|p'r|I> <L|q'|J> V(I,J,t) (-1)^N(I)
    #
    reset!(ket_a)
    for K in 1:ket_a.max

        reset!(ket_b)
        for L in 1:ket_b.max
            for p in 1:no
                for r in 1:no
                    copy!(bra_a,ket_a)
                    apply_annihilation!(bra_a,r)
                    bra_a.sign != 0 || continue
                    apply_creation!(bra_a,p)
                    bra_a.sign != 0 || continue
                    calc_linear_index!(bra_a,_binomial)
                    I = bra_a.lin_index

                    for q in 1:no

                        copy!(bra_b,ket_b)
                        apply_creation!(bra_b,q)
                        bra_b.sign != 0 || continue
                        calc_linear_index!(bra_b,_binomial)
                        J = bra_b.lin_index

                        @views tdm_pqr = tdm[:,:,p,q,r]
                        @views v1_IJ = v1[:,I,J]
                        @views v2_KL = v2[:,K,L]
                        #sgn = bra_a.sign * bra_b.sign
                        sgn = sgnK * bra_a.sign * bra_b.sign

                        if sgn == 1
                            @tensor begin
                                tdm_pqr[s,t] += v1_IJ[s] * v2_KL[t]
                            end
                        elseif sgn == -1
                            @tensor begin
                                tdm_pqr[s,t] -= v1_IJ[s] * v2_KL[t]
                            end
                        else
                            throw(Exception)
                        end
                    end #q
                end #r
            end #p
            incr!(ket_b)
        end #L
        incr!(ket_a)
    end #K
    tdm = permutedims(tdm, [3,4,5,1,2])
    return tdm
#=}}}=#
end

"""
    compute_ABb(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix, spin_case)

Compute represntation of 2creation operators between states `bra_v` and `ket_v`
# Arguments
- `no`: number of orbitals
- `bra_na`: number of alpha electrons in bra
- `bra_nb`: number of beta electrons in bra
- `ket_na`: number of alpha electrons in ket
- `ket_nb`: number of beta electrons in ket
- `bra_v`: basis vectors in bra
- `ket_v`: basis vectors ket
- `spin_case`: either (alpha|beta) giving ABa or BAb
"""
function compute_ABb(no::Integer, bra_na, bra_nb, ket_na, ket_nb, bra_v::Matrix, ket_v::Matrix)
#={{{=#
    bra_a = DeterminantString(no,bra_na)
    bra_b = DeterminantString(no,bra_nb)

    ket_a = DeterminantString(no,ket_na)
    ket_b = DeterminantString(no,ket_nb)

    bra_na == ket_na+1 || throw(DimensionMismatch)
    bra_nb == ket_nb+0 || throw(DimensionMismatch)

    #
    # precompute lookup table for binomial coefficients
    _binomial = Array{Int,2}(undef,no+1,no+1)
    for i in 0:no
        for j in i:no
            _binomial[j+1,i+1] = binomial(j,i)
        end
    end

    bra_M = size(bra_v,2)
    ket_M = size(ket_v,2)

    v1 = reshape(bra_v, bra_a.max, bra_b.max, bra_M)
    v2 = reshape(ket_v, ket_a.max, ket_b.max, ket_M)
    v1 = permutedims(v1, [3,1,2])
    v2 = permutedims(v2, [3,1,2])

    #
    #   TDM[pqr,s,t] =
    tdm = zeros(Float64, bra_M, ket_M, no, no, no)
    bra_a = deepcopy(ket_a)
    bra_b = deepcopy(ket_b)

    # here, q'r are Beta
    #  v(K,L,s) <KL|p'q'r|IJ> V(I,J,t)
    # -v(K,L,s) <L|<K|p'|I>q'r|J> V(I,J,t)
    # -v(K,L,s) <K|p'|I> <L|q'r|J> V(I,J,t)
    #
    reset!(ket_a)
    for K in 1:ket_a.max

        reset!(ket_b)
        for L in 1:ket_b.max
            for p in 1:no
                copy!(bra_a,ket_a)
                apply_creation!(bra_a,p)
                bra_a.sign != 0 || continue
                calc_linear_index!(bra_a,_binomial)
                I = bra_a.lin_index

                for q in 1:no
                    for r in 1:no

                        copy!(bra_b,ket_b)
                        apply_annihilation!(bra_b,r)
                        apply_creation!(bra_b,q)

                        bra_b.sign != 0 || continue
                        calc_linear_index!(bra_b,_binomial)
                        J = bra_b.lin_index

                        @views tdm_pqr = tdm[:,:,p,q,r]
                        @views v1_IJ = v1[:,I,J]
                        @views v2_KL = v2[:,K,L]
                        sgn = bra_a.sign * bra_b.sign

                        if sgn == 1
                            @tensor begin
                                tdm_pqr[s,t] += v1_IJ[s] * v2_KL[t]
                            end
                        elseif sgn == -1
                            @tensor begin
                                tdm_pqr[s,t] -= v1_IJ[s] * v2_KL[t]
                            end
                        else
                            error(" sign?: ", sgn)
                        end
                    end #r
                end #q
            end #p
            incr!(ket_b)
        end #L
        incr!(ket_a)
    end #K
    tdm = permutedims(tdm, [3,4,5,1,2])
    return tdm
#=}}}=#
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
@testset "numeric fluents" begin

# Test numeric fluent functionality
path = joinpath(dirname(pathof(PDDL)), "..", "test", "numeric")

domain = load_domain(joinpath(path, "domain.pddl"))
@test domain.name == Symbol("zeno-travel")
@test convert(Term, domain.functions[:distance]) == pddl"(distance ?c1 ?c2)"
@test domain.functions[:fuel].argtypes == (:aircraft,)

problem = load_problem(joinpath(path, "problem.pddl"))
@test problem.metric == pddl"(minimize (+ (* 4 (total-time)) (* 5 (total-fuel-used))))"

problem
# Test for static functions
static_fluents = infer_static_fluents(domain)
@test :capacity in static_fluents
@test length(static_fluents) == 5

state = initstate(domain, problem)
implementations = [
    "concrete interpreter" => (domain, state),
    "abstract interpreter" => abstracted(domain, state),
    "ground interpreter" => (ground(domain, state), state),
    "concrete compiler" => compiled(domain, state),
    "abstract compiler" => compiled(abstracted(domain), state)
]

@testset "numeric fluents ($name)" for (name, (domain, _)) in implementations

    # Execute plan to goal
    state = initstate(domain, problem)

    # Person 1 boards plane 1
    state = execute(domain, state, pddl"(board person1 plane1 city0)")
    @test domain[state => pddl"(onboard plane1)"] ≃ 1

    # Plane 1 flies from city 0 to city 2
    state = execute(domain, state, pddl"(fly plane1 city0 city2)")
    @test domain[state => pddl"(total-fuel-used)"] ≃ 3100

    # Person 1 debarks at city 2
    state = execute(domain, state, pddl"(debark person1 plane1 city2)")
    @test domain[state => pddl"(at person1 city2)"] ≃ true

    # Plane 1 refuels at city 2
    state = execute(domain, state, pddl"(refuel plane1 city2)")
    @test domain[state => pddl"(fuel plane1)"] ≃ 10232

    # Person 2 boards plane 1 at city2
    state = execute(domain, state, pddl"(board person2 plane1 city2)")
    @test domain[state => pddl"(in person2 plane1)"] ≃ true

    # Plane 1 flies from city 2 to city 0
    state = execute(domain, state, pddl"(fly plane1 city2 city0)")
    @test domain[state => pddl"(total-fuel-used)"] ≃ 6200

    # Person 2 debarks at city 0
    state = execute(domain, state, pddl"(debark person2 plane1 city0)")
    @test domain[state => pddl"(at person2 city0)"] ≃ true

    # Plane 1 refuels at city 0
    state = execute(domain, state, pddl"(refuel plane1 city0)")
    @test domain[state => pddl"(fuel plane1)"] ≃ 10232

    # Plane 1 zooms from city 0 to city 1
    state = execute(domain, state, pddl"(zoom plane1 city0 city1)")
    @test domain[state => pddl"(total-fuel-used)"] ≃ 16370

    # The whole plan took 9 steps
    @test domain[state => pddl"(total-time)"] ≃ 9

    # Check if goal is satisfied
    @test satisfy(domain, state, problem.goal) == true

    # Test execution of entire plan
    state = initstate(domain, problem)
    plan = @pddl(
        "(board person1 plane1 city0)",
        "(fly plane1 city0 city2)",
        "(debark person1 plane1 city2)",
        "(refuel plane1 city2)",
        "(board person2 plane1 city2)",
        "(fly plane1 city2 city0)",
        "(debark person2 plane1 city0)",
        "(refuel plane1 city0)",
        "(zoom plane1 city0 city1)"
    )
    sim = EndStateSimulator()
    state = sim(domain, state, plan)
    @test satisfy(domain, state, problem.goal) == true

end

end # numeric fluents


text:
#=##############################################################################################
Filename: solve_vlm_bem.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: `<: Action` function calculates aerodynamic forces on a blade-element momentum-vortex
        lattice system; only accounts for rotor-on-wing effects
=###############################################################################################

"""
    solve_vlm_bem(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol) <: Action

Solves for the aerodynamic forces at each step.

# Arguments:

* `aircraft::Aircraft` : `Aircraft` system object
* `parameters<:Parameters` `Parameters` struct
* `freestream::Freestream` : `Freestream` object
* `environment::Environment` `Environment` object
* `step_range::AbstractArray` : array of steps for which the simulation is run
* `stepi::Int` : index of the current step
* `step_symbol::String` : defines the step, e.g. `alpha` or `time`

`parameters <: Parameters` requires the following elements:

* `omegas::Vector{Float64}` : a vector of length `length(step_range)` containing a vector of rotational velocities for each rotor
* `Js::Array{Float64,2}` : a vector of length `length(step_range)` containing a vector of advance ratios for each rotor
* `Ts::Array{Float64,2}` : a vector of length `length(step_range)` containing a vector of dimensional thrust values for each rotor
* `Qs::Array{Float64,2}` : a vector of length `length(step_range)` containing a vector of dimensional torque values for each rotor
* `us::Vector{Vector{Vector{Float64}}}` : each [i][j][k]th element is the axial induced velocity at ith step of the jth rotor at the kth radial section
* `vs::Vector{Vector{Vector{Float64}}}` : each [i][j][k]th element is the swirl induced velocity at ith step of the jth rotor at the kth radial section
* `wake_function::Vector{Function}`: [i]th element is a function vwake(X::Vector{Float64}) describing the wake induced velocity at `X` at the ith step
* `wake_shape_functions::Vector{Function}` : [i]th element is a function f(Rtip, x) describing the radial distance from the rotor axis to the boundary of the wake of the ith rotor
* `axial_interpolations::Vector{Function}` : [i]th element is a function f(rs, us, r, Rtip) that returns the axial component of rotor-induced velocity at distance r from the rotor axis based on the calculated axial induced velocities output from CCBlade of the ith rotor
* `swirl_interpolations::Vector{Function}` : [i]th element is a function f(rs, us, r, Rtip) that returns the swirl component of rotor-induced velocity at distance r from the rotor axis based on the calculated axial induced velocities output from CCBlade of the ith rotor
* `axial_multipliers::Vector{Function}` : [i]th element is a function f(distance2plane, Rtip) that is multiplied by the axial induced velocity function of the ith rotor
* `swirl_multipliers::Vector{Function}` : [i]th element is a function f(distance2plane, Rtip) that is multiplied by the swirl induced velocity function of the ith rotor
* `CFs::Array{Float64,2}` : [:,j]th element is the [CD,CY,CL] force coefficients of the aircraft at the jth step
* `CMs::Array{Float64,2}` : [:,j]th element is the [CMx,CMy,CMz] moment coefficients of the aircraft at the jth step
* `cfs::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of force coefficients cd, cy, cl corresponding to the ith step, jth lifting surface
* `cms::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of moment coefficients cmx, cmy, cmz corresponding to the ith step, jth lifting surface

"""
function solve_vlm_bem(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)

    flags = Vector{Bool}(undef,4)
    # solves rotors
    flags[1] = solve_rotor_nondimensional(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)
    # get wake function
    flags[2] = solve_rotor_wake(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)
    # solve VLM
    flags[3] = solve_wing_CF_CM(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)
    # extract lift and moment distribution
    flags[4] = lift_moment_distribution_blownwing(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)
    # flags[4] = lift_moment_distribution_blownwing(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)

    return prod(flags)
end

"""
    solve_vlm_bem(aircraft, step_range)

Method returns initialized elements required for the `parameters <: Parameters` struct during simulation.

# Arguments:

* `aircraft::Aircraft` : system to be simulated
* `step_range::AbstractArray` : defines each step of the simulation

# Returns:

* `omegas::Vector{Float64}` : a vector of rotational speeds in rad/s at the current step
* `Js::Array{Float64,2}` : each [i,j]th element is the advance ratio of the ith rotor at the jth step
* `Ts::Array{Float64,2}` : each [i,j]th element is the thrust of the ith rotor at the jth step
* `Qs::Array{Float64,2}` : each [i,j]th element is the torque of the ith rotor at the jth step
* `CTs::Array{Float64,2}` : each [i,j]th element is the thrust coefficient of the ith rotor at the jth step
* `CQs::Array{Float64,2}` : each [i,j]th element is the torque coefficient of the ith rotor at the jth step
* `ηs::Array{Float64,2}` : each [i,j]th element is the propulsive efficiency of the ith rotor at the jth step
* `us::Vector{Vector{Vector{Float64}}}` : each [i][j][k]th element is the axial induced velocity at ith step of the jth rotor at the kth radial section
* `vs::Vector{Vector{Vector{Float64}}}` : each [i][j][k]th element is the swirl induced velocity at ith step of the jth rotor at the kth radial section
* `wake_function::Vector{Function}`: [i]th element is a function vwake(X::Vector{Float64}) describing the wake induced velocity at `X` at the ith step
* `wake_shape_functions::Vector{Function}` : [i]th element is a function f(Rtip, x) describing the radial distance from the rotor axis to the boundary of the wake of the ith rotor
* `axial_interpolations::Vector{Function}` : [i]th element is a function f(rs, us, r, Rtip) that returns the axial component of rotor-induced velocity at distance r from the rotor axis based on the calculated axial induced velocities output from CCBlade of the ith rotor
* `swirl_interpolations::Vector{Function}` : [i]th element is a function f(rs, us, r, Rtip) that returns the swirl component of rotor-induced velocity at distance r from the rotor axis based on the calculated axial induced velocities output from CCBlade of the ith rotor
* `axial_multipliers::Vector{Function}` : [i]th element is a function f(distance2plane, Rtip) that is multiplied by the axial induced velocity function of the ith rotor
* `swirl_multipliers::Vector{Function}` : [i]th element is a function f(distance2plane, Rtip) that is multiplied by the swirl induced velocity function of the ith rotor
* `CFs::Array{Float64,2}` : [:,j]th element is the [CD,CY,CL] force coefficients of the aircraft at the jth step
* `CMs::Array{Float64,2}` : [:,j]th element is the [CMx,CMy,CMz] moment coefficients of the aircraft at the jth step
* `cfs::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of force coefficients cd, cy, cl corresponding to the ith step, jth lifting surface
* `cms::Vector{Vector{Array{Float64,2}}}`: [i][j]th element is an array of size (3,nspanwisepanels) of moment coefficients cmx, cmy, cmz corresponding to the ith step, jth lifting surface

"""
function solve_vlm_bem(aircraft, step_range)

    params_solve_rotor_nondimensional = solve_rotor_nondimensional(aircraft, step_range) # omegas, Js, Ts, Qs, CTs, CQs, ηs, us, vs
    params_solve_rotor_wake = solve_rotor_wake(aircraft, step_range) # wake_function, us, vs
    params_solve_wing_CF_CM = solve_wing_CF_CM(aircraft, step_range) # wake_function, CFs, CMs
    params_lift_moment_distribution = lift_moment_distribution(aircraft, step_range) # cfs, cms

    return params_solve_rotor_nondimensional..., params_solve_rotor_wake[[1,4,5,6,7,8]]..., params_solve_wing_CF_CM[2:end]..., params_lift_moment_distribution...
end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
########
# Set up parameters
########
parameters = (
    xmax = 10000,   # horizontal extent of the SingleStack in the x direction [m]
    ymax = 500,     # horizontal extent of the SingleStack in the y direction [m]
    zmax = 10000,   # vertical extent of the SingleStack [m]
    g = 9.81,    # gravitation acceleration [m /s^2]
    pₒ = 1.01325e5,                                               # initial surface pressure [Pa]
    cp_d = (8.3144598 / 28.97e-3) / (2 / 7),                          # specific heat for ideal gas at constant pressure [J / K / kg]
    cv_d = (8.3144598 / 28.97e-3) / (2 / 7) - 8.3144598 / 28.97e-3,   # specific heat for ideal gas at constant volume [J / K / kg]
    R_d = 8.3144598 / 28.97e-3,                                    # gas constant for dry air [J / K / kg]
    κ = 2.0 / 7.0, # R_d/cp_d
    cₛ = 340,     # speed of sound [m / s]
    c_s = 800,     # specific heat for land (soil)  [J / K / kg]
    κ_s = 0.0,     # soil conductivity [W / m / K] (set to 0 for energy conservation checks)
    h_s = 100,     # depth of the modelled soil model [m]
    T_h = 280,     # temperature of soil at depth h [K]
    ρ_s = 1500,    # density for land (soil) [kg / m^3]
    τ = 0.9,     # atmospheric transmissivity
    α = 0.5,     # surface albedo
    σ = 5.67e-8, # Steffan Boltzmann constant [kg / s^3 / K^4]
    g_a = 0.06,    # aerodynamic conductance for heat transfer [kg / m^2 / s]
    ϵ = 0.98,    # broadband emissivity / absorptivity of the surface
    F_a = 0.0,     # downward LW flux from the atmosphere [W / m^2]
    F_sol = 1361,   # incoming solar TOA radiation [W / m^2]
    τ_d = 10,     # idealized daily cycle period [s]
)

########
# Set up inital conditions
########

# Land (slab) initial condition
T_sfc₀(p, x, y, z) = p.T_h

# Atmos (single stack) initial conditions
θ₀Atmos(p, x, y, z) = 300.0 + 10.0 * (z / p.zmax)
π_exnerAtmos(p, x, y, z) = 1.0 - p.g / (p.cp_d * θ₀Atmos(p, x, y, z)) * z
ρ₀Atmos(p, x, y, z) = p.pₒ / (p.R_d * θ₀Atmos(p, x, y, z)) * (π_exnerAtmos(p, x, y, z))^(p.cv_d / p.R_d)
ρθ₀Atmos(p, x, y, z) = ρ₀Atmos(p, x, y, z) * θ₀Atmos(p, x, y, z)
ρu₀Atmos(p, x, y, z) = @SVector [1, 0, 0]

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using Test
using HMCResearchRandomWalks.Enz

"testing functions for locations"
function rand_x() return Enz.ESCAPE_X * rand() * 2 - Enz.ESCAPE_X end
function rand_y() return Enz.ESCAPE_Y * rand() end # greater than 0.0

@testset "constant locations" begin
    @test Enz.inwater([rand_x(), 5.0])
    @test Enz.inwater([rand_x(), 149.0])

    @test Enz.insensor([0.0, -0.001])
    @test Enz.insensor([140.0, -0.001])
    @test Enz.insensor([-140.0, -0.001])
    @test !Enz.insensor([140.1, -0.001])
    @test !Enz.insensor([-140.1, -0.001])
    @test !Enz.insensor([rand_x(), 0.001])

    @test Enz.inwalls(140.1, -0.001)
    @test Enz.inwalls(-140.1, -0.001)
    @test Enz.inwalls(280.0, -0.001)
    @test Enz.inwalls(-280.0, -0.001)
    @test !Enz.inwalls(rand_x(), rand_y())

    @test !Enz.inppd([rand_x(), rand_y()])
end

@testset "thick/thin locations" begin
    Enz.set_THICK_ENZ(true)
    @test Enz.inenz([0.0, 1.5])
    @test !Enz.inwater([0.0, 1.5])

    Enz.set_THICK_ENZ(false)
    @test !Enz.inenz([0.0, 1.5])
    @test Enz.inwater([0.0, 1.5])
    @test Enz.inenz([0.0, 0.15])
    @test !Enz.inwater([0.0, 0.15])
end

@testset "boundarycheck cases" begin
    @test Enz.sensorcases(rand_x()) == "sensor"

    initxy = [150.0, 0.003]
    dx = 0.0
    dy = -1.0
    ss = 0.005
    output = Enz.wallcases(initxy, dx, dy, ss)
    @test output[1] == [150.0, 0.0005]
    @test output[2] == "no collision"
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
@testset "BLSPeriodogram printing" begin
    t, y, dy, params = make_data()
    result = BLS(t, y, dy; params.duration)

    N_in = length(t)
    N_out = length(result.periods)
    T = eltype(t)
    V = eltype(y)
    # shortform
    @test sprint(show, result) == "$N_out-element BLSPeriodogram{$T,$V}"

    # longform
    minP, maxP = extrema(result.periods)
    pars = BoxLeastSquares.params(result)
    @test repr("text/plain", result) == """
    BLSPeriodogram
    ==============
    input dim: $N_in
    output dim: $N_out
    period range: $minP - $maxP
    duration range: $(params.duration) - $(params.duration)
    objective: $(string(result.objective))

    parameters
    ----------
    index: $(pars.index)
    period: $(pars.period)
    duration: $(pars.duration)
    t0: $(pars.t0)
    depth: $(pars.depth) ± $(pars.depth_err)
    snr: $(pars.snr)
    log-likelihood: $(pars.loglike)"""
end


text:
module FourPlayer

import Cribbage: CribbageGame, GameState, Player, nextplayer, UnexpectedPlayerException
import Cribbage.TwoPlayer: incrementscore, numplayers, cardsdealt, cardstocrib, players, reportscore, isgameover, rotatedealer

export incrementscore, isgameover, players, reportscore, rotatedealer, FourPlayerGame, FourPlayerGameState

struct FourPlayerGame <: CribbageGame
    p₁ :: Player
    p₂ :: Player
    p₃ :: Player
    p₄ :: Player
    winningscore :: Int8
    function FourPlayerGame(
        p₁::Player, p₂::Player, p₃::Player, p₄::Player, winningscore::Integer=121)
        @assert(winningscore == 121 || winningscore == 61, "winningscore is either 61 or 121")
        @assert(p₁ ≢ p₂)
        @assert(p₁ ≢ p₃)
        @assert(p₁ ≢ p₄)
        @assert(p₂ ≢ p₃)
        @assert(p₂ ≢ p₄)
        @assert(p₃ ≢ p₄)
        new(p₁, p₂, p₃, p₄, winningscore)
    end
end

struct FourPlayerGameState <: GameState
    game :: FourPlayerGame
    dealer :: Player
    s₁ :: Int8
    s₂ :: Int8
    function FourPlayerGameState(
        game::FourPlayerGame, dealer::Player, s₁::Integer=0, s₂::Integer=0)
        @assert(dealer ∈ players(game), "dealer is not one of the players of the game")
        @assert(0 ≤ s₁ ≤ game.winningscore, "s₁ is an invalid score")
        @assert(0 ≤ s₂ ≤ game.winningscore, "s₂ is an invalid score")
        new(game, dealer, s₁, s₂)
    end
end

function numplayers(game::FourPlayerGame)::Int8
    return 4
end

function cardsdealt(game::FourPlayerGame)::Int8
    return 5
end

function cardstocrib(game::FourPlayerGame)::Int8
    return 0
end

function players(game::FourPlayerGame)::Array{Player,1}
    return [ game.p₁, game.p₂, game.p₃, game.p₄ ]
end

function incrementscore(
        gamestate::FourPlayerGameState,
        player::Player,
        amount::Integer
    )
    @assert(amount ≥ 0)
    if player === gamestate.game.p₁ || player === gamestate.game.p₃
        return FourPlayerGameState(
            gamestate.game,
            gamestate.dealer,
            min(gamestate.s₁ + amount, gamestate.game.winningscore),
            gamestate.s₂
        )
    elseif player === gamestate.game.p₂ || player === gamestate.game.p₄
        return FourPlayerGameState(
            gamestate.game,
            gamestate.dealer,
            gamestate.s₁,
            min(gamestate.s₂ + amount, gamestate.game.winningscore)
        )
    else
        throw(UnexpectedPlayerException("$(player.name) is not one of the players in this game"))
    end
end

function isgameover(gamestate::FourPlayerGameState)::Bool
    return gamestate.s₁ == gamestate.game.winningscore || gamestate.s₂ == gamestate.game.winningscore
end

function reportscore(gamestate::FourPlayerGameState)::String
    return "$(gamestate.s₁) - $(gamestate.s₂)"
end

function rotatedealer(gamestate::FourPlayerGameState)::FourPlayerGameState
    return FourPlayerGameState(
        gamestate.game,
        nextplayer(players(gamestate.game), gamestate.dealer),
        gamestate.s₁,
        gamestate.s₂
    )
end

end  # module FourPlayer


text:
module SpringCollab2020BlackholeCustomColors

using ..Ahorn, Maple

@mapdef Effect "SpringCollab2020/BlackholeCustomColors" BlackholeCustomColors(only::String="*", exclude::String="", colorsMild::String="6e3199,851f91,3026b0", colorsWild::String="ca4ca7,b14cca,ca4ca7",
    bgColorInner::String="000000", bgColorOuterMild::String="512a8b", bgColorOuterWild::String="bd2192", alpha::Number=1.0)

placements = BlackholeCustomColors

function Ahorn.canFgBg(effect::BlackholeCustomColors)
    return true, true
end

end


text:
# Correlation analysis of signals
#
#  autocorrelation
#  cross-correlation
#  partial autocorrelation
#

#######################################
#
#   Helper functions
#
#######################################

default_laglen(lx::Int) = min(lx-1, round(Int,10*log10(lx)))
check_lags(lx::Int, lags::AbstractVector) = (maximum(lags) < lx || error("lags must be less than the sample length."))

function demean_col!(z::AbstractVector{T}, x::AbstractMatrix{T}, j::Int, demean::Bool) where T<:RealFP
    m = size(x, 1)
    @assert m == length(z)
    b = m * (j-1)
    if demean
        s = zero(T)
        for i = 1 : m
            s += x[b + i]
        end
        mv = s / m
        for i = 1 : m
            z[i] = x[b + i] - mv
        end
    else
        copyto!(z, 1, x, b+1, m)
    end
    z
end

#######################################
#
#   Auto-correlations
#
#######################################

default_autolags(lx::Int) = 0 : default_laglen(lx)

_autodot(x::AbstractVector{<:RealFP}, lx::Int, l::Int) = dot(x, 1:lx-l, x, 1+l:lx)

## autocov
"""
    autocov!(r, x, lags; demean=true)

Compute the autocovariance of a vector or matrix `x` at `lags` and store the result
in `r`. `demean` denotes whether the mean of `x` should be subtracted from `x`
before computing the autocovariance.

If `x` is a vector, `r` must be a vector of the same length as `x`.
If `x` is a matrix, `r` must be a matrix of size `(length(lags), size(x,2))`, and
where each column in the result will correspond to a column in `x`.

The output is not normalized. See [`autocor!`](@ref) for a method with normalization.
"""
function autocov!(r::RealVector, x::AbstractVector{T}, lags::IntegerVector; demean::Bool=true) where T<:RealFP
    lx = length(x)
    m = length(lags)
    length(r) == m || throw(DimensionMismatch())
    check_lags(lx, lags)

    z::Vector{T} = demean ? x .- mean(x) : x
    for k = 1 : m  # foreach lag value
        r[k] = _autodot(z, lx, lags[k]) / lx
    end
    return r
end

function autocov!(r::RealMatrix, x::AbstractMatrix{T}, lags::IntegerVector; demean::Bool=true) where T<:RealFP
    lx = size(x, 1)
    ns = size(x, 2)
    m = length(lags)
    size(r) == (m, ns) || throw(DimensionMismatch())
    check_lags(lx, lags)

    z = Vector{T}(undef, lx)
    for j = 1 : ns
        demean_col!(z, x, j, demean)
        for k = 1 : m
            r[k,j] = _autodot(z, lx, lags[k]) / lx
        end
    end
    return r
end

"""
    autocov(x, [lags]; demean=true)

Compute the autocovariance of a vector or matrix `x`, optionally specifying
the `lags` at which to compute the autocovariance. `demean` denotes whether
the mean of `x` should be subtracted from `x` before computing the autocovariance.

If `x` is a vector, return a vector of the same length as `x`.
If `x` is a matrix, return a matrix of size `(length(lags), size(x,2))`,
where each column in the result corresponds to a column in `x`.

When left unspecified, the lags used are the integers from 0 to
`min(size(x,1)-1, 10*log10(size(x,1)))`.

The output is not normalized. See [`autocor`](@ref) for a function with normalization.
"""
function autocov(x::AbstractVector{T}, lags::IntegerVector; demean::Bool=true) where T<:Real
    autocov!(Vector{fptype(T)}(undef, length(lags)), float(x), lags; demean=demean)
end

function autocov(x::AbstractMatrix{T}, lags::IntegerVector; demean::Bool=true) where T<:Real
    autocov!(Matrix{fptype(T)}(undef, length(lags), size(x,2)), float(x), lags; demean=demean)
end

autocov(x::AbstractVecOrMat{<:Real}; demean::Bool=true) = autocov(x, default_autolags(size(x,1)); demean=demean)

## autocor

"""
    autocor!(r, x, lags; demean=true)

Compute the autocorrelation function (ACF) of a vector or matrix `x` at `lags`
and store the result in `r`. `demean` denotes whether the mean of `x` should
be subtracted from `x` before computing the ACF.

If `x` is a vector, `r` must be a vector of the same length as `x`.
If `x` is a matrix, `r` must be a matrix of size `(length(lags), size(x,2))`, and
where each column in the result will correspond to a column in `x`.

The output is normalized by the variance of `x`, i.e. so that the lag 0
autocorrelation is 1. See [`autocov!`](@ref) for the unnormalized form.
"""
function autocor!(r::RealVector, x::AbstractVector{T}, lags::IntegerVector; demean::Bool=true) where T<:RealFP
    lx = length(x)
    m = length(lags)
    length(r) == m || throw(DimensionMismatch())
    check_lags(lx, lags)

    z::Vector{T} = demean ? x .- mean(x) : x
    zz = dot(z, z)
    for k = 1 : m  # foreach lag value
        r[k] = _autodot(z, lx, lags[k]) / zz
    end
    return r
end

function autocor!(r::RealMatrix, x::AbstractMatrix{T}, lags::IntegerVector; demean::Bool=true) where T<:RealFP
    lx = size(x, 1)
    ns = size(x, 2)
    m = length(lags)
    size(r) == (m, ns) || throw(DimensionMismatch())
    check_lags(lx, lags)

    z = Vector{T}(undef, lx)
    for j = 1 : ns
        demean_col!(z, x, j, demean)
        zz = dot(z, z)
        for k = 1 : m
            r[k,j] = _autodot(z, lx, lags[k]) / zz
        end
    end
    return r
end

"""
    autocor(x, [lags]; demean=true)

Compute the autocorrelation function (ACF) of a vector or matrix `x`,
optionally specifying the `lags`. `demean` denotes whether the mean
of `x` should be subtracted from `x` before computing the ACF.

If `x` is a vector, return a vector of the same length as `x`.
If `x` is a matrix, return a matrix of size `(length(lags), size(x,2))`,
where each column in the result corresponds to a column in `x`.

When left unspecified, the lags used are the integers from 0 to
`min(size(x,1)-1, 10*log10(size(x,1)))`.

The output is normalized by the variance of `x`, i.e. so that the lag 0
autocorrelation is 1. See [`autocov`](@ref) for the unnormalized form.
"""
function autocor(x::AbstractVector{T}, lags::IntegerVector; demean::Bool=true) where T<:Real
    autocor!(Vector{fptype(T)}(undef, length(lags)), float(x), lags; demean=demean)
end

function autocor(x::AbstractMatrix{T}, lags::IntegerVector; demean::Bool=true) where T<:Real
    autocor!(Matrix{fptype(T)}(undef, length(lags), size(x,2)), float(x), lags; demean=demean)
end

autocor(x::AbstractVecOrMat{<:Real}; demean::Bool=true) = autocor(x, default_autolags(size(x,1)); demean=demean)

#######################################
#
#   Cross-correlations
#
#######################################

default_crosslags(lx::Int) = (l=default_laglen(lx); -l:l)

_crossdot(x::AbstractVector{T}, y::AbstractVector{T}, lx::Int, l::Int) where {T<:RealFP} =
    (l >= 0 ? dot(x, 1:lx-l, y, 1+l:lx) : dot(x, 1-l:lx, y, 1:lx+l))

## crosscov

"""
    crosscov!(r, x, y, lags; demean=true)

Compute the cross covariance function (CCF) between real-valued vectors or matrices
`x` and `y` at `lags` and store the result in `r`. `demean` specifies whether the
respective means of `x` and `y` should be subtracted from them before computing their
CCF.

If both `x` and `y` are vectors, `r` must be a vector of the same length as
`lags`. If either `x` is a matrix and `y` is a vector, `r` must be a matrix of size
`(length(lags), size(x, 2))`; if `x` is a vector and `y` is a matrix, `r` must be a matrix
of size `(length(lags), size(y, 2))`. If both `x` and `y` are matrices, `r` must be a
three-dimensional array of size `(length(lags), size(x, 2), size(y, 2))`.

The output is not normalized. See [`crosscor!`](@ref) for a function with normalization.
"""
function crosscov!(r::RealVector, x::AbstractVector{T}, y::AbstractVector{T}, lags::IntegerVector; demean::Bool=true) where T<:RealFP
    lx = length(x)
    m = length(lags)
    (length(y) == lx && length(r) == m) || throw(DimensionMismatch())
    check_lags(lx, lags)

    zx::Vector{T} = demean ? x .- mean(x) : x
    zy::Vector{T} = demean ? y .- mean(y) : y
    for k = 1 : m  # foreach lag value
        r[k] = _crossdot(zx, zy, lx, lags[k]) / lx
    end
    return r
end

function crosscov!(r::RealMatrix, x::AbstractMatrix{T}, y::AbstractVector{T}, lags::IntegerVector; demean::Bool=true) where T<:RealFP
    lx = size(x, 1)
    ns = size(x, 2)
    m = length(lags)
    (length(y) == lx && size(r) == (m, ns)) || throw(DimensionMismatch())
    check_lags(lx, lags)

    zx = Vector{T}(undef, lx)
    zy::Vector{T} = demean ? y .- mean(y) : y
    for j = 1 : ns
        demean_col!(zx, x, j, demean)
        for k = 1 : m
            r[k,j] = _crossdot(zx, zy, lx, lags[k]) / lx
        end
    end
    return r
end

function crosscov!(r::RealMatrix, x::AbstractVector{T}, y::AbstractMatrix{T}, lags::IntegerVector; demean::Bool=true) where T<:RealFP
    lx = length(x)
    ns = size(y, 2)
    m = length(lags)
    (size(y, 1) == lx && size(r) == (m, ns)) || throw(DimensionMismatch())
    check_lags(lx, lags)

    zx::Vector{T} = demean ? x .- mean(x) : x
    zy = Vector{T}(undef, lx)
    for j = 1 : ns
        demean_col!(zy, y, j, demean)
        for k = 1 : m
            r[k,j] = _crossdot(zx, zy, lx, lags[k]) / lx
        end
    end
    return r
end

function crosscov!(r::AbstractArray{T,3}, x::AbstractMatrix{T}, y::AbstractMatrix{T}, lags::IntegerVector; demean::Bool=true) where T<:RealFP
    lx = size(x, 1)
    nx = size(x, 2)
    ny = size(y, 2)
    m = length(lags)
    (size(y, 1) == lx && size(r) == (m, nx, ny)) || throw(DimensionMismatch())
    check_lags(lx, lags)

    # cached (centered) columns of x
    zxs = Vector{T}[]
    sizehint!(zxs, nx)
    for j = 1 : nx
        xj = x[:,j]
        if demean
            mv = mean(xj)
            for i = 1 : lx
                xj[i] -= mv
            end
        end
        push!(zxs, xj)
    end

    zx = Vector{T}(undef, lx)
    zy = Vector{T}(undef, lx)
    for j = 1 : ny
        demean_col!(zy, y, j, demean)
        for i = 1 : nx
            zx = zxs[i]
            for k = 1 : m
                r[k,i,j] = _crossdot(zx, zy, lx, lags[k]) / lx
            end
        end
    end
    return r
end

"""
    crosscov(x, y, [lags]; demean=true)

Compute the cross covariance function (CCF) between real-valued vectors or
matrices `x` and `y`, optionally specifying the `lags`. `demean` specifies
whether the respective means of `x` and `y` should be subtracted from them
before computing their CCF.

If both `x` and `y` are vectors, return a vector of the same length as
`lags`. Otherwise, compute cross covariances between each pairs of columns in `x` and `y`.

When left unspecified, the lags used are the integers from
`-min(size(x,1)-1, 10*log10(size(x,1)))` to `min(size(x,1), 10*log10(size(x,1)))`.

The output is not normalized. See [`crosscor`](@ref) for a function with normalization.
"""
function crosscov(x::AbstractVector{T}, y::AbstractVector{T}, lags::IntegerVector; demean::Bool=true) where T<:Real
    crosscov!(Vector{fptype(T)}(undef, length(lags)), float(x), float(y), lags; demean=demean)
end

function crosscov(x::AbstractMatrix{T}, y::AbstractVector{T}, lags::IntegerVector; demean::Bool=true) where T<:Real
    crosscov!(Matrix{fptype(T)}(undef, length(lags), size(x,2)), float(x), float(y), lags; demean=demean)
end

function crosscov(x::AbstractVector{T}, y::AbstractMatrix{T}, lags::IntegerVector; demean::Bool=true) where T<:Real
    crosscov!(Matrix{fptype(T)}(undef, length(lags), size(y,2)), float(x), float(y), lags; demean=demean)
end

function crosscov(x::AbstractMatrix{T}, y::AbstractMatrix{T}, lags::IntegerVector; demean::Bool=true) where T<:Real
    crosscov!(Array{fptype(T),3}(undef, length(lags), size(x,2), size(y,2)), float(x), float(y), lags; demean=demean)
end

crosscov(x::AbstractVecOrMat{T}, y::AbstractVecOrMat{T}; demean::Bool=true) where {T<:Real} = crosscov(x, y, default_crosslags(size(x,1)); demean=demean)

## crosscor
"""
    crosscor!(r, x, y, lags; demean=true)

Compute the cross correlation between real-valued vectors or matrices `x` and `y` at
`lags` and store the result in `r`. `demean` specifies whether the respective means of
`x` and `y` should be subtracted from them before computing their cross correlation.

If both `x` and `y` are vectors, `r` must be a vector of the same length as
`lags`. If either `x` is a matrix and `y` is a vector, `r` must be a matrix of size
`(length(lags), size(x, 2))`; if `x` is a vector and `y` is a matrix, `r` must be a matrix
of size `(length(lags), size(y, 2))`. If both `x` and `y` are matrices, `r` must be a
three-dimensional array of size `(length(lags), size(x, 2), size(y, 2))`.

The output is normalized by `sqrt(var(x)*var(y))`. See [`crosscov!`](@ref) for the
unnormalized form.
"""
function crosscor!(r::RealVector, x::AbstractVector{T}, y::AbstractVector{T}, lags::IntegerVector; demean::Bool=true) where T<:RealFP
    lx = length(x)
    m = length(lags)
    (length(y) == lx && length(r) == m) || throw(DimensionMismatch())
    check_lags(lx, lags)

    zx::Vector{T} = demean ? x .- mean(x) : x
    zy::Vector{T} = demean ? y .- mean(y) : y
    sc = sqrt(dot(zx, zx) * dot(zy, zy))
    for k = 1 : m  # foreach lag value
        r[k] = _crossdot(zx, zy, lx, lags[k]) / sc
    end
    return r
end

function crosscor!(r::RealMatrix, x::AbstractMatrix{T}, y::AbstractVector{T}, lags::IntegerVector; demean::Bool=true) where T<:RealFP
    lx = size(x, 1)
    ns = size(x, 2)
    m = length(lags)
    (length(y) == lx && size(r) == (m, ns)) || throw(DimensionMismatch())
    check_lags(lx, lags)

    zx = Vector{T}(undef, lx)
    zy::Vector{T} = demean ? y .- mean(y) : y
    yy = dot(zy, zy)
    for j = 1 : ns
        demean_col!(zx, x, j, demean)
        sc = sqrt(dot(zx, zx) * yy)
        for k = 1 : m
            r[k,j] = _crossdot(zx, zy, lx, lags[k]) / sc
        end
    end
    return r
end

function crosscor!(r::RealMatrix, x::AbstractVector{T}, y::AbstractMatrix{T}, lags::IntegerVector; demean::Bool=true) where T<:RealFP
    lx = length(x)
    ns = size(y, 2)
    m = length(lags)
    (size(y, 1) == lx && size(r) == (m, ns)) || throw(DimensionMismatch())
    check_lags(lx, lags)

    zx::Vector{T} = demean ? x .- mean(x) : x
    zy = Vector{T}(undef, lx)
    xx = dot(zx, zx)
    for j = 1 : ns
        demean_col!(zy, y, j, demean)
        sc = sqrt(xx * dot(zy, zy))
        for k = 1 : m
            r[k,j] = _crossdot(zx, zy, lx, lags[k]) / sc
        end
    end
    return r
end

function crosscor!(r::AbstractArray{T,3}, x::AbstractMatrix{T}, y::AbstractMatrix{T}, lags::IntegerVector; demean::Bool=true) where T<:RealFP
    lx = size(x, 1)
    nx = size(x, 2)
    ny = size(y, 2)
    m = length(lags)
    (size(y, 1) == lx && size(r) == (m, nx, ny)) || throw(DimensionMismatch())
    check_lags(lx, lags)

    # cached (centered) columns of x
    zxs = Vector{T}[]
    sizehint!(zxs, nx)
    xxs = Vector{T}(undef, nx)

    for j = 1 : nx
        xj = x[:,j]
        if demean
            mv = mean(xj)
            for i = 1 : lx
                xj[i] -= mv
            end
        end
        push!(zxs, xj)
        xxs[j] = dot(xj, xj)
    end

    zx = Vector{T}(undef, lx)
    zy = Vector{T}(undef, lx)
    for j = 1 : ny
        demean_col!(zy, y, j, demean)
        yy = dot(zy, zy)
        for i = 1 : nx
            zx = zxs[i]
            sc = sqrt(xxs[i] * yy)
            for k = 1 : m
                r[k,i,j] = _crossdot(zx, zy, lx, lags[k]) / sc
            end
        end
    end
    return r
end

"""
    crosscor(x, y, [lags]; demean=true)

Compute the cross correlation between real-valued vectors or matrices `x` and `y`,
optionally specifying the `lags`. `demean` specifies whether the respective means of
`x` and `y` should be subtracted from them before computing their cross correlation.

If both `x` and `y` are vectors, return a vector of the same length as
`lags`. Otherwise, compute cross covariances between each pairs of columns in `x` and `y`.

When left unspecified, the lags used are the integers from
`-min(size(x,1)-1, 10*log10(size(x,1)))` to `min(size(x,1), 10*log10(size(x,1)))`.

The output is normalized by `sqrt(var(x)*var(y))`. See [`crosscov`](@ref) for the
unnormalized form.
"""
function crosscor(x::AbstractVector{T}, y::AbstractVector{T}, lags::IntegerVector; demean::Bool=true) where T<:Real
    crosscor!(Vector{fptype(T)}(undef, length(lags)), float(x), float(y), lags; demean=demean)
end

function crosscor(x::AbstractMatrix{T}, y::AbstractVector{T}, lags::IntegerVector; demean::Bool=true) where T<:Real
    crosscor!(Matrix{fptype(T)}(undef, length(lags), size(x,2)), float(x), float(y), lags; demean=demean)
end

function crosscor(x::AbstractVector{T}, y::AbstractMatrix{T}, lags::IntegerVector; demean::Bool=true) where T<:Real
    crosscor!(Matrix{fptype(T)}(undef, length(lags), size(y,2)), float(x), float(y), lags; demean=demean)
end

function crosscor(x::AbstractMatrix{T}, y::AbstractMatrix{T}, lags::IntegerVector; demean::Bool=true) where T<:Real
    crosscor!(Array{fptype(T),3}(undef, length(lags), size(x,2), size(y,2)), float(x), float(y), lags; demean=demean)
end

crosscor(x::AbstractVecOrMat{T}, y::AbstractVecOrMat{T}; demean::Bool=true) where {T<:Real} = crosscor(x, y, default_crosslags(size(x,1)); demean=demean)

#######################################
#
#   Partial auto-correlations
#
#   TODO: the codes below need cleanup.
#
#######################################

function pacf_regress!(r::RealMatrix, X::AbstractMatrix{T}, lags::IntegerVector, mk::Integer) where T<:RealFP
    lx = size(X, 1)
    tmpX = ones(T, lx, mk + 1)
    for j = 1 : size(X,2)
        for l = 1 : mk
            for i = 1+l:lx
                tmpX[i,l+1] = X[i-l,j]
            end
        end
        for i = 1 : length(lags)
            l = lags[i]
            sX = view(tmpX, 1+l:lx, 1:l+1)
            r[i,j] = l == 0 ? 1 : (cholesky!(sX'sX, Val(false)) \ (sX'view(X, 1+l:lx, j)))[end]
        end
    end
    r
end

function pacf_yulewalker!(r::RealMatrix, X::AbstractMatrix{T}, lags::IntegerVector, mk::Integer) where T<:RealFP
    tmp = Vector{T}(undef, mk)
    for j = 1 : size(X,2)
        acfs = autocor(X[:,j], 1:mk)
        for i = 1 : length(lags)
            l = lags[i]
            r[i,j] = l == 0 ? 1 : l == 1 ? acfs[i] : -durbin!(view(acfs, 1:l), tmp)[l]
        end
    end
end

"""
    pacf!(r, X, lags; method=:regression)

Compute the partial autocorrelation function (PACF) of a matrix `X` at `lags` and
store the result in `r`. `method` designates the estimation method. Recognized values
are `:regression`, which computes the partial autocorrelations via successive
regression models, and `:yulewalker`, which computes the partial autocorrelations
using the Yule-Walker equations.

`r` must be a matrix of size `(length(lags), size(x, 2))`.
"""
function pacf!(r::RealMatrix, X::AbstractMatrix{T}, lags::IntegerVector; method::Symbol=:regression) where T<:RealFP
    lx = size(X, 1)
    m = length(lags)
    minlag, maxlag = extrema(lags)
    (0 <= minlag && 2maxlag < lx) || error("Invalid lag value.")
    size(r) == (m, size(X,2)) || throw(DimensionMismatch())

    if method == :regression
        pacf_regress!(r, X, lags, maxlag)
    elseif method == :yulewalker
        pacf_yulewalker!(r, X, lags, maxlag)
    else
        error("Invalid method: $method")
    end
    return r
end

"""
    pacf(X, lags; method=:regression)

Compute the partial autocorrelation function (PACF) of a real-valued vector
or matrix `X` at `lags`. `method` designates the estimation method. Recognized
values are `:regression`, which computes the partial autocorrelations via successive
regression models, and `:yulewalker`, which computes the partial autocorrelations
using the Yule-Walker equations.

If `x` is a vector, return a vector of the same length as `lags`.
If `x` is a matrix, return a matrix of size `(length(lags), size(x, 2))`,
where each column in the result corresponds to a column in `x`.
"""
function pacf(X::AbstractMatrix{T}, lags::IntegerVector; method::Symbol=:regression) where T<:Real
    pacf!(Matrix{fptype(T)}(undef, length(lags), size(X,2)), float(X), lags; method=method)
end

function pacf(x::AbstractVector{T}, lags::IntegerVector; method::Symbol=:regression) where T<:Real
    vec(pacf(reshape(x, length(x), 1), lags, method=method))
end


text:
module PropertyDicts

export PropertyDict

struct PropertyDict{K, V, T <: AbstractDict{K, V}} <: AbstractDict{K, V}
    d::T

    PropertyDict(d::T) where {T <: AbstractDict} = new{keytype(d), valtype(d), T}(d)
    PropertyDict(args...) = PropertyDict(Dict(args...))
end

unwrap(d::PropertyDict) = getfield(d, :d)

function Base.sizehint!(d::PropertyDict, n::Integer)
    sizehint!(unwrap(d), n)
    return d
end

Base.push!(d::PropertyDict, p::Pair) = push!(unwrap(d), p)
Base.pop!(d::PropertyDict, args...) = pop!(unwrap(d), args...)
function Base.empty!(d::PropertyDict)
    empty!(unwrap(d))
    return d
end
function Base.delete!(d::PropertyDict, key)
    delete!(unwrap(d), key)
    return d
end

Base.getproperty(d::PropertyDict, n::Symbol) = getindex(d, n)
Base.getproperty(d::PropertyDict, n::String) = getindex(d, n)

Base.setproperty!(d::PropertyDict, n::Symbol, v) = setindex!(d, v, n)
Base.setproperty!(d::PropertyDict, n::String, v) = setindex!(d, v, n)

Base.convert(::Type{Any}, d::PropertyDict) = d
Base.convert(::Type{PropertyDict{K,V,T}}, d::PropertyDict{K,V,T}) where {K,V,T<:AbstractDict{K,V}} = d
Base.convert(::Type{T}, d::PropertyDict) where T <: AbstractDict = T === AbstractDict ? d : convert(T, PropertyDicts.unwrap(d))
Base.convert(::Type{T}, d::PropertyDict) where T = convert(T, PropertyDicts.unwrap(d))

Base.get(d::PropertyDict, k, default) = get(unwrap(d), k, default)
Base.get(d::PropertyDict{Symbol}, k::AbstractString, default) = get(d, Symbol(k), default)
Base.get(d::PropertyDict{<:AbstractString}, k::Symbol, default) = get(d, String(k), default)
function Base.get(d::PropertyDict{Any}, k::Symbol, default)
    out = get(unwrap(d), k, Base.secret_table_token)
    if out === Base.secret_table_token
        return get(unwrap(d), String(k), default)
    else
        return out
    end
end
function Base.get(d::PropertyDict{Any}, k::String, default)
    out = get(unwrap(d), k, Base.secret_table_token)
    if out === Base.secret_table_token
        return get(unwrap(d), Symbol(k), default)
    else
        return out
    end
end

function Base.get(f::Union{Function,Type}, d::PropertyDict, k)
    out = get(d, k, Base.secret_table_token)
    if out === Base.secret_table_token
        return f()
    else
        return out
    end
end

Base.get!(d::PropertyDict, k, default) = get!(unwrap(d), k, default)
Base.get!(d::PropertyDict{Symbol}, k::AbstractString, default) = get!(d, Symbol(k), default)
Base.get!(d::PropertyDict{<:AbstractString}, k::Symbol, default) = get!(d, String(k), default)
function Base.get!(f::Union{Function,Type}, d::PropertyDict, k)
    out = get(d, k, Base.secret_table_token)
    if out === Base.secret_table_token
        default = f()
        setindex!(d, default, k)
        return default
    else
        return out
    end
end

function Base.getindex(d::PropertyDict, k)
    out = get(d, k, Base.secret_table_token)
    out === Base.secret_table_token && throw(KeyError(k))
    return out
end

Base.setindex!(d::PropertyDict, v, i) = setindex!(unwrap(d), v, i)
Base.setindex!(d::PropertyDict{<:AbstractString}, v, i::Symbol) = setindex!(d, v, String(i))
Base.setindex!(d::PropertyDict{Symbol}, v, i::AbstractString) = setindex!(d, v, Symbol(i))

Base.iterate(d::PropertyDict) = iterate(unwrap(d))
Base.iterate(d::PropertyDict, i) = iterate(unwrap(d), i)

Base.IteratorSize(::Type{PropertyDict{K,V,T}}) where {K,V,T} = Base.IteratorSize(T)
Base.IteratorEltype(::Type{PropertyDict{K,V,T}}) where {K,V,T} = Base.IteratorEltype(T)

Base.length(d::PropertyDict) = length(unwrap(d))

Base.string(d::PropertyDict) = string(unwrap(d))

@static if isdefined(Base, :hasproperty)
    Base.hasproperty(d::PropertyDict, key::Symbol) = haskey(d, key)
    Base.hasproperty(d::PropertyDict, key) = haskey(d, key)
end
Base.haskey(d::PropertyDict, key) = haskey(unwrap(d), key)
Base.haskey(d::PropertyDict{<:AbstractString}, key::Symbol) = haskey(d, String(key))
Base.haskey(d::PropertyDict{Symbol}, key::AbstractString) = haskey(d, Symbol(key))
function Base.haskey(d::PropertyDict{Any}, key::AbstractString)
    haskey(unwrap(d), key) || haskey(unwrap(d), Symbol(key))
end
function Base.haskey(d::PropertyDict{Any}, key::Symbol)
    haskey(unwrap(d), key) || haskey(unwrap(d), String(key))
end

# a handful of dictionaries aren't just wrapped in `KeySet` and `ValueIterator`
Base.keys(d::PropertyDict) = keys(unwrap(d))
Base.values(d::PropertyDict) = values(unwrap(d))

Base.propertynames(d::PropertyDict) = keys(d)

end # module PropertyDicts


text:
"""
Neighborhoods define how surrounding cells are related to the current cell.
The `neighbors` function returns the sum of surrounding cells, as defined
by the neighborhood.
"""
abstract type Neighborhood{R} end

radius(hood::Neighborhood{R}) where R = R

"""
A Moore-style neighborhood where a square are with a center radius `(D - 1) / 2`
where D is the diameter.
"""
abstract type AbstractRadialNeighborhood{R} <: Neighborhood{R} end

"""
Radial neighborhoods calculate the surrounding neighborood
from the radius around the central cell. The central cell
is ommitted.
"""
struct RadialNeighborhood{R} <: AbstractRadialNeighborhood{R} end

"""
    neighbors(hood::AbstractRadialNeighborhood, buf, state)

Sums radial Moore nieghborhoods of any dimension.
"""
neighbors(hood::AbstractRadialNeighborhood, model, buf, state) = sum(buf) - state

@inline mapreduceneighbors(f, data, hood::AbstractRadialNeighborhood, rule, state, index) = begin
    r = radius(hood)
    sum = zero(state)
    # Loop over dispersal kernel grid dimensions
    for x = one(r):2r + one(r)
        xdest = x + index[2] - r - one(r)
        for y = one(r):2r + one(r)
            x == (r + one(r)) && y == (r + one(r)) && continue
            ydest = y + index[1] - r - one(r)
            hood_index = (y, x)
            dest_index = (ydest, xdest)
            sum += f(data, hood, rule, state, hood_index, dest_index)
        end
    end
    sum
end

"""
Custom neighborhoods are tuples of custom coordinates (also tuples) specified in relation
to the central point of the current cell. They can be any arbitrary shape or size, but
should be listed in column-major order for performance.
"""
abstract type AbstractCustomNeighborhood{R} <: Neighborhood{R} end

"""
Allows completely arbitrary neighborhood shapes by specifying each coordinate
specifically.
"""
@description struct CustomNeighborhood{R,C} <: AbstractCustomNeighborhood{R}
    coords::C | "A tuple of tuples of Int, containing 2-D coordinates relative to the central point"
end
CustomNeighborhood(args...) = CustomNeighborhood(args)
CustomNeighborhood(coords) = CustomNeighborhood{absmaxcoord(coords), typeof(coords)}(coords)

coords(hood::CustomNeighborhood) = hood.coords
# Calculate the maximum absolute value in the coords to use as the radius
absmaxcoord(coords) = maximum((x -> maximum(abs.(x))).(coords))

neighbors(hood::CustomNeighborhood, rule, buf, state) = begin
    r = radius(hood); sum = zero(state)
    for coord in coords(hood)
        sum += buf[(coord .+ r .+ 1)...]
    end
    sum
end

@inline mapreduceneighbors(f, data, hood::CustomNeighborhood, rule, state, index) = begin
    r = radius(hood); sum = zero(state)
    # Loop over dispersal kernel grid dimensions
    for coord in coords(hood)
        hood_index = coord .+ r
        dest_index = index .+ coord
        sum += f(data, hood, rule, state, hood_index, dest_index)
    end
    sum
end

"""
Sets of custom neighborhoods that can have separate rules for each set.
"""
@description struct LayeredCustomNeighborhood{R,C} <: AbstractCustomNeighborhood{R}
    layers::C | "A tuple of tuple of custom neighborhoods"
end
LayeredCustomNeighborhood(l::Tuple) =
    LayeredCustomNeighborhood{maximum(absmaxcoord.(coords.(l))), typeof(l)}(l)

@inline neighbors(hood::LayeredCustomNeighborhood, rule, buf, state) =
    map(layer -> neighbors(layer, rule, buf, state), hood.layers)

@inline mapreduceneighbors(f, data, hood::LayeredCustomNeighborhood, rule, state, index) =
    map(layer -> mapreduceneighbors(f, data, layer, rule, state, index), hood.layers)

"""
A convenience wrapper to build a VonNeumann neighborhoods as a `CustomNeighborhood`.

# TODO: variable radius
"""
VonNeumannNeighborhood() = CustomNeighborhood(((0,-1), (-1,0), (1,0), (0,1)))

"""
Find the largest radius present in the passed in rules.
"""
radius(set::Ruleset) = max(radius(rules(ruleset)))
radius(set::MultiRuleset) = begin
    ruleradius = map(radius, map(rules, ruleset(set)))
    intradius = map(key -> radius(interactions(set), key), map(Val, keys(set)))
    map(max, Tuple(ruleradius), Tuple(intradius))
end
radius(ruleset::Ruleset) = radius(rules(ruleset))
radius(rules::Tuple) = mapreduce(radius, max, rules)
radius(rules::Tuple, key) = mapreduce(rule -> radius(rule, key), max, rules)
radius(rules::Tuple{}, args...) = 0

radius(rule::NeighborhoodRule) = radius(neighborhood(rule))
radius(rule::PartialNeighborhoodRule) = radius(neighborhood(rule))
radius(rule::Rule, args...) = 0
# Only the first rule in a chain can have a radius larger than zero.
radius(chain::Chain) = radius(chain[1])


text:
using Distributions

"""Equivalent gaussian sigma for small log-probability.

Return the equivalent gaussian sigma corresponding to the natural log of
the cumulative gaussian probability logp. In other words, return x, such
that Q(x) = p, where Q(x) is the cumulative normal distribution. This
version uses the rational approximation from Abramowitz and Stegun,
eqn 26.2.23, that claims to be precise to ~1e-4. Using the log(P) as input
gives a much extended range.

The parameters here are the result of a best-fit, with no physical meaning.

Translated from Scott Ransom's PRESTO
"""
function extended_equiv_gaussian_Nsigma(logp::Number)
    t = sqrt(-2.0 * logp)
    num = 2.515517 + t * (0.802853 + t * 0.010328)
    denom = 1.0 + t * (1.432788 + t * (0.189269 + t * 0.001308))
    return t - num / denom
end

"""Number of Gaussian sigmas corresponding to tail log-probability.

This function computes the value of the characteristic function of a
standard Gaussian distribution for the tail probability equivalent to the
provided p-value, and turns this value into units of standard deviations
away from the Gaussian mean. This allows the user to make a statement
about the signal such as “I detected this pulsation at 4.1 sigma

The example values below are obtained by brute-force integrating the
Gaussian probability density function using the mpmath library
between Nsigma and +inf.

# Examples
```jldoctest
julia> using PulsarSearch

julia> pvalues = Array([0.15865525393145707, 0.0013498980316301035, 9.865877e-10, 6.22096e-16, 3.0567e-138]);

julia> log_pvalues = log.(pvalues);

julia> sigmas = Array([1, 3, 6, 8, 25]);

julia> all(isapprox(equivalent_gaussian_Nsigma_from_logp.(log_pvalues), sigmas; atol=0.1))
true
```
"""
function equivalent_gaussian_Nsigma_from_logp(logp)

    if logp < -300
        # print("Extended")
        return extended_equiv_gaussian_Nsigma(logp)
    end
    d = Normal()
    return cquantile(d, exp(logp))
end

export equivalent_gaussian_Nsigma_from_logp

"""Number of Gaussian sigmas corresponding to tail probability.

This function computes the value of the characteristic function of a
standard Gaussian distribution for the tail probability equivalent to the
provided p-value, and turns this value into units of standard deviations
away from the Gaussian mean. This allows the user to make a statement
about the signal such as “I detected this pulsation at 4.1 sigma
"""
function equivalent_gaussian_Nsigma(p)
    return equivalent_gaussian_Nsigma_from_logp(log(p))
end

export equivalent_gaussian_Nsigma

"""Asymptotic natural log of incomplete gamma function.

Return the natural log of the incomplete gamma function in
its asymptotic limit as z->infty.  This is from Abramowitz
and Stegun eqn 6.5.32.

Translated from Scott Ransom's PRESTO

# Examples

```jldoctest
julia> using PulsarSearch

julia> pvalues = Array([0.15865525393145707, 0.0013498980316301035]);

julia> sigmas = Array([1, 3]);

julia> all(isapprox(equivalent_gaussian_Nsigma.(pvalues), sigmas; atol=0.1))
true

```

"""
function log_asymptotic_incomplete_gamma(a, z)
    x = 1.0
    newxpart = 1.0
    term = 1.0
    ii = 1

    while (abs(newxpart) > 1e-15)
        term *= (a - ii)
        newxpart = term / z^ii
        x += newxpart
        ii += 1
    end

    return (a - 1.0) * log(z) - z + log(x)
end

export log_asymptotic_incomplete_gamma

"""Natural log of the Gamma function in its asymptotic limit.

Return the natural log of the gamma function in its asymptotic limit
as z->infty.  This is from Abramowitz and Stegun eqn 6.1.41.

Translated from Scott Ransom's PRESTO
"""
function log_asymptotic_gamma(z)
    half_log_twopi = 0.91893853320467267  # (1/2)*log(2*pi)
    one_twelfth = 8.3333333333333333333333e-2
    one_degree = 2.7777777777777777777778e-3  # 1 / 360
    one_over_1680 = 5.9523809523809529e-4
    one_over_1260 = 7.9365079365079365079365e-4
    x = (z - 0.5) * log(z) - z + half_log_twopi
    y = 1.0 / (z * z)
    x +=
        (
            ((-one_over_1680 * y + one_over_1260) * y - one_degree) * y +
            one_twelfth
        ) / z
    return x
end

"""Log survival function of the chi-squared distribution.

# Examples

```jldoctest chi2_logp
julia> using PulsarSearch;

julia> using Distributions;

julia> chi2 = 31;

julia> d = Chisq(2);

julia> isapprox(chi2_logp(chi2, 2), logccdf(d, chi2), atol=0.1)
true

julia> chi2 = Array([5, 32]);

julia> all(isapprox.(chi2_logp.(chi2, 2), logccdf.(d, chi2), atol=0.1))
true
```
"""
function chi2_logp(chi2, dof)

    # If very large reduced chi squared, use approximation. This is an
    # eyeballed limit parameter space where the difference between the
    # approximation and the scipy version is tiny, but above which the scipy
    # version starts failing.
    if (chi2 / dof > 15.0) || ((dof > 150) && (chi2 / dof > 6.0))
        return log_asymptotic_incomplete_gamma(0.5 * dof, 0.5 * chi2) -
               log_asymptotic_gamma(0.5 * dof)
    end
    d = Chisq(dof)
    return logccdf(d, chi2)
end

export chi2_logp

"""Calculate a multi-trial p-value from the log of a single-trial one.

This allows to work around Numba's limitation on longdoubles, a way to
vectorize the computation when we need longdouble precision.

Parameters
----------
logp1 : float
    The natural logarithm of the significance at which we reject the null
    hypothesis on each single trial.
n : int
    The number of trials

Returns
-------
logpn : float
    The log of the significance at which we reject the null hypothesis
    after multiple trials
"""
function logp_multitrial_from_single_logp(logp1, n)
    # If the the probability is very small (p1 * n) < 1e-6, use Bonferroni
    # approximation.
    logn = log(n)
    if logp1 + logn < -7
        return logp1 + logn
    end

    return log(1 - (1 - exp(logp1))^n)
end

export logp_multitrial_from_single_logp

raw"""Calculate a multi-trial p-value from a single-trial one.

Calling *p* the probability of a single success, the Binomial
distributions says that the probability *at least* one outcome
in n trials is

``P(k\geq 1) = \sum_{k\geq 1} \binom{n}{k} p^k (1-p)^{(n-k)}``

or more simply, using P(k ≥ 0) = 1

P(k\geq 1) = 1 - \binom{n}{0} (1-p)^n = 1 - (1-p)^n``

Parameters
----------
p1 : float
    The significance at which we reject the null hypothesis on
    each single trial.
n : int
    The number of trials

Returns
-------
pn : float
    The significance at which we reject the null hypothesis
    after multiple trials
"""
function p_multitrial_from_single_trial(p1, n)
    logpn = logp_multitrial_from_single_logp(log(p1), n)
    return exp(logpn)
end

export p_multitrial_from_single_trial

"""Calculate a multi-trial p-value from the log of a single-trial one.

This allows to work around Numba's limitation on longdoubles, a way to
vectorize the computation when we need longdouble precision.

Parameters
----------
logpn : float
    The natural logarithm of the significance at which we want to reject
    the null hypothesis after multiple trials
n : int
    The number of trials

Returns
-------
logp1 : float
    The log of the significance at which we reject the null hypothesis on
    each single trial.
"""
function logp_single_trial_from_logp_multitrial(logpn, n)
    logn = log(n)
    # If the the probability is very small, use Bonferroni approximation.
    if logpn < -7
        return logpn - logn
    end
    # Numerical errors arise when pn is very close to 1.
    p1 = 1 - (1 - exp(logpn))^(1 / n)
    return log(p1)
end

export logp_single_trial_from_logp_multitrial

raw"""Calculate the single-trial p-value from a total p-value

Let us say that we want to reject a null hypothesis at the
``pn`` level, after executing ``n`` different measurements.
This might be the case because, e.g., we
want to have a 1% probability of detecting a signal in an
entire power spectrum, and we need to correct the detection
level accordingly.

The typical procedure is dividing the initial probability
(often called _epsilon_) by the number of trials. This is
called the Bonferroni correction and it is often a good
approximation, when ``pn`` is low: ``p1 = pn / n``.

However, if ``pn`` is close to 1, this approximation gives
incorrect results.

Here we calculate this probability by inverting the Binomial
problem. Given that (see ``p_multitrial_from_single_trial``)
the probability of getting more than one hit in n trials,
given the single-trial probability *p*, is

``P (k \geq 1) =  1 - (1 - p)^n``

we get the single trial probability from the multi-trial one
from

``p = 1 - (1 - P)^{(1/n)}``

This is also known as Šidák correction.

Parameters
----------
pn : float
    The significance at which we want to reject the null
    hypothesis after multiple trials
n : int
    The number of trials

Returns
-------
p1 : float
    The significance at which we reject the null hypothesis on
    each single trial.
"""
function p_single_trial_from_p_multitrial(pn, n)
    logp = logp_single_trial_from_logp_multitrial(log(pn), n)

    return exp(logp)
end

export p_single_trial_from_p_multitrial

"""Calculate the probability of a certain folded profile, due to noise.

Parameters
----------
z2 : float
    A ``Z^2_n`` statistics value
n : int, default 2
    The ``n`` in ``Z^2_n`` (number of harmonics, including the fundamental)

Other Parameters
----------------
ntrial : int
    The number of trials executed to find this profile
n_summed_spectra : int
    Number of ``Z_2^n`` periodograms that were averaged to obtain z2

Returns
-------
p : float
    The probability that the ``Z^2_n`` value has been produced by noise
"""
function z2_n_probability(z2, n; ntrial = 1, n_summed_spectra = 1)
    d = Chisq(2 * n * n_summed_spectra)
    epsilon_1 = ccdf(d, z2 * n_summed_spectra)
    epsilon = p_multitrial_from_single_trial(epsilon_1, ntrial)
    return epsilon
end

export z2_n_probability

"""Calculate the probability of a certain folded profile, due to noise.

Parameters
----------
z2 : float
    A ``Z^2_n`` statistics value
n : int, default 2
    The ``n`` in ``Z^2_n`` (number of harmonics, including the fundamental)

Other Parameters
----------------
ntrial : int
    The number of trials executed to find this profile
n_summed_spectra : int
    Number of ``Z_2^n`` periodograms that were averaged to obtain z2

Returns
-------
p : float
    The probability that the ``Z^2_n`` value has been produced by noise
"""
function z2_n_logprobability(z2, n; ntrial = 1, n_summed_spectra = 1)

    epsilon_1 = chi2_logp(z2 * n_summed_spectra, 2 * n * n_summed_spectra)
    epsilon = logp_multitrial_from_single_logp(epsilon_1, ntrial)
    return epsilon
end

export z2_n_logprobability

"""Return the detection level for the ``Z^2_n`` statistics.

See Buccheri et al. (1983), Bendat and Piersol (1971).

Parameters
----------
n : int, default 2
    The ``n`` in ``Z^2_n`` (number of harmonics, including the fundamental)
epsilon : float, default 0.01
    The fractional probability that the signal has been produced by noise

Other Parameters
----------------
ntrial : int
    The number of trials executed to find this profile
n_summed_spectra : int
    Number of Z_2^n periodograms that are being averaged

Returns
-------
detlev : float
    The epoch folding statistics corresponding to a probability
    epsilon * 100 % that the signal has been produced by noise
"""
function z2_n_detection_level(
    n = 2,
    epsilon = 0.01;
    ntrial = 1,
    n_summed_spectra = 1,
)
    epsilon = p_single_trial_from_p_multitrial(epsilon, ntrial)
    d = Chisq(2 * n_summed_spectra * n)
    retlev = cquantile(d, epsilon) / (n_summed_spectra)
    return retlev
end

export z2_n_detection_level


text:
"""
```
Reaction(name, rate, formula)
```

Define a `Reaction` with a `name`, a stochastic `rate` constant, and a `formula`.

The `formula` should use `Species` from a `Network`. Multiple reactants are separated by a `+` symbol; similarly for products. Reactants and products are separated by a `-->` symbol. Refer to the examples below:

# Examples

```
Reaction("birth", 1.0, "X --> X + X")   # equiv. X --> 2*X
Reaction("death", 0.5, "X --> 0")       # no products
Reaction("immigration", 0.5, "0 --> X") # no reactants
```

"""
struct Reaction
    identifier :: Symbol
    rate :: Float64

    reactants :: OrderedDict{Symbol,Int}
    products  :: OrderedDict{Symbol,Int}

    affectedby :: Vector{Symbol}
    affects    :: Vector{Symbol}

    origex :: Expr

    function Reaction(name, rate, formula)
        reactants, products, origex = parse_reaction(formula)

        if isempty(reactants) && isempty(products)
            error("Reaction needs at least one reactant or product.")
        end

        if any(x -> x < 0, values(reactants)) || any(x -> x < 0, values(products))
            error("Coefficients must be positive.")
        end

        affectedby = Symbol[]
        affects = collect(keys(reactants))

        for s in union(keys(reactants), keys(products))
          vm = get(reactants, s, 0)
          vp = get(products, s, 0)

          v = vp - vm

          v != 0 && push!(affectedby, s)
        end

        return new(Symbol(replace(name, " " => "_")), rate, reactants, products, affectedby, affects, origex)
    end
end

function Base.show(io::IO, x::Reaction)
    print_participants(io, x.reactants)
    print(io, " -> ")
    print_participants(io, x.products)
    println(io)
    println(io, "  affected by: $(join(x.affectedby, ", "))")
    println(io, "  affects: $(join(x.affects, ", "))")
end

function print_participants(io, participants)
  n = length(participants)
  if n == 0
    print(io, "∅")
  else
    i = 1
    for (id, coeff) in participants
      if coeff == 1
        print(io, id)
      else
        print(io, coeff, " ", id)
      end

      if i < n
        print(io, " + ")
      end
      i = i + 1
    end
  end
end

parse_reaction(formula::String) = parse_reaction(Meta.parse(formula))

function parse_reaction(ex::Expr)
  reactants = OrderedDict{Symbol,Int}()
  products  = OrderedDict{Symbol,Int}()

  if ex.head == :-->
    exr = ex.args[1]
    exp = ex.args[2]

    add_participants!(reactants, exr)
    add_participants!(products,  exp)
  else
    throw("malformed reaction")
  end
  return reactants, products, ex
end

function add_participants!(dict, ex)
  if isa(ex, Symbol)
    #push!(dict, ex)
    val = get(dict, ex, 0)
    dict[ex] = val + 1
  elseif isa(ex, Expr)
    if ex.args[1] == :*
      #push!(dict, ex.args[3], ex.args[2])
      val = get(dict, ex.args[3], 0)
      dict[ex.args[3]] = val + ex.args[2]
    else
      for i in 2:length(ex.args)
        add_participants!(dict, ex.args[i])
      end
    end
  end
end


text:
s = CuStream()
synchronize(s)
@test CUDA.query(s) == true

let s2 = CuStream()
    @test s != s2
    @test !(s == s2)
end

let s3 = CuStream(; flags=CUDA.STREAM_NON_BLOCKING)
    @test s != s3
    @test !(s == s3)
end

prio = priority_range()
let s = CuStream(; priority=first(prio))
    @test priority(s) == first(prio)
end
let s = CuStream(; priority=last(prio))
    @test priority(s) == last(prio)
end

synchronize(CuDefaultStream())