Reference

calc_fsky(win, wmodes)

This functions returns a measure of the sky fraction covered by the survey with window win. The exact implementation is considered an implementation detail and can change in the future.

pixwin(cmodes)

Return the angular pixel window for application to a Clnn object.

xyz2rtp(xyz)

Convert the Cartesian positions in xyz to spherical coordinates rθϕ. The first dimension is of length 3, the second is the number of galaxies. Assumes a flat geometry.

GNLs

Module to define radial basis functions. The specific basis functions are defined in an enum, BoundaryConditions.

BoundaryConditions

This enum defines the radial basis functions.

GNL(nmax, lmax, rmin, rmax)
GNL(kmax, rmin, rmax; nmax=typemax(Int64), lmax=typemax(Int64))

Generate gnl(n,l,r). Returns a struct that can be called for calculating gnl. Note that the last argument is r, not kr.

SphericalBesselGNL(nmax, lmax, rmin, rmax)
SphericalBesselGNL(kmax, rmin, rmax; nmax=typemax(Int64), lmax=typemax(Int64))

Generate gnl(n,l,r). Returns a struct that can be called for calculating gnl. Note that the last argument is r, not kr.

calc_first_n_zeros(func, nmax; δ=π/20, xmin=0.0)

Calculate the first nmax zeros of the function func. Assumes that zeros are spaced more than δ apart. First zero could be xmin.

calc_knl(nmax, lmax, rmin, rmax; boundary=potential)
calc_knl(kmax, rmin, rmax; nmax=typemax(Int64), lmax=typemax(Int64), boundary=potential)

Calculate the knl for boundary conditions.

calc_zeros(func, xmin, xmax; δ=π/20)

Calculate all zeros of the function func between xmin and xmax. Assumes that zeros are spaced more than δ apart.

CryoGNLs

To generate the radial basis functions, first define the radial binning, e.g.,

rmin = 1000.0
rmax = 4000.0
nbins = 1000
Δr = (rmax - rmin) / nbins
rmid = range(rmin + Δr/2, rmax - Δr/2, length=nbins)

Then, to generate the radial basis functions at ℓ=3,

ell = 3
k, g, ginv = CryoGNLs.get_radial_cryofunks(ell, rmid, Δr)

The functions gnl(r) are now saved in g, so that g[i,n] is gnl(rmid[i]) and the corresponding k-mode is k[n].

AnlmModes(kmax, rmin, rmax; cache=true, nside=nothing, boundary=SFB.GNL.potential)
AnlmModes(nmax, lmax, rmin, rmax; cache=true, nside=nothing, boundary=SFB.GNL.potential)

This is where we define which modes are included. As our criterium, we set a maximum qmax = kmax * rmax, and we include all modes below that.

The modes are arranged in the following order. The fastest loop is through 'm', then 'l', finally 'n', from small number to larger number. We restrict 'm' to m >= 0, and we assume a real field.

Example:

julia> kmax = 0.05
julia> rmin = 500.0
julia> rmax = 1000.0
julia> modes = AnlmModes(kmax, rmin, rmax)

ClnnModes(::AnlmModes)
ClnnModes(::AnlmModes, ::AnlmModes)

This is where we define which modes are included in the power spectrum, given a AnlmModes struct.

The modes are arranged in the following order. The fastest loop is through 'n̄', then 'Δn', finally 'ℓ', from small number to larger number. We make the convention that n̄ is the smaller of the k-modes, and Δn >= 0.

More useful is the labeling by ℓ, n₁, n₂. In that case we make the convention that n₁ = n̄ and n₂ = n̄ + Δn.

If two AnlmModes are given, then we assume a cross-correlation is desired. The S parameter specifies whether the resulting modes are symmetric (S=true) on interchange of k1 and k2, or not (S=false). This is useful for auto- and cross-correlation, respectively.

There is some ambiguity in what exactly the S parameter encodes. First, it could encode purely the symmetry of the Clnn modes, or it could encode whether we are doing a cross-correlation or an auto-correlation. Are these two equivalent? Can we ever have a symmetric cross-correlation? We define the cross-correlation as <delta^A delta^{B,*}>, so that for a fixedl,Clnn^{AB}is the transpose ofClnn^{BA}`.

For the purposes of ClnnModes, the S parameter shall encode whether the two amodes are equal or not. That is, S = (amodesA==amodesB). Whether it encodes any symmetry A==B beyond that will be left to the application using ClnnModes. The recommendation is to indeed use it to encode A==B.

getlkk(::ClnnModes, [i])
getlkk(::ClnnBinnedModes, [i])

Get the physical modes ℓ, k, and k' corresponding to the index i. If i is left out, an array lkk of all modes is returned so that lkk[1,:] are all the ℓ-values, lkk[2,:] all the k-values, and lkk[3,:] are all the k'-values.

exponentiate(s::SeparableArray, exp::Number)

This function is used for elementwise exponentiation of the array 's'. It could be made more elegant by extending the broadcast syntax. PRs welcome.

cat2amln(rθϕ, amodes, nbar, win_rhat_ln, weights)

Computes the spherical Fourier-Bessel decomposition coefficients from a catalogue of sources. The number density is measured from the survey as $\bar n = N_\mathrm{gals} / V_\mathrm{eff}$.

weights is an array containing a weight for each galaxy.

Example

julia> using SphericalFourierBesselDecompositions
julia> cat2amln(rθϕ, ...)

winweights2galweights(weights, wmodes, rθϕ)

Returns an array with the weight for each galaxy. weights is a 2D-array where the first index goes over r, the second over healpix pixel p.

ConfigurationSpaceModes(rmin, rmax, nr, nside)

A struct to describe and define the voxelization scheme.

apply_window(rθϕ, win, wmodes::ConfigurationSpaceModes; rng=Random.GLOBAL_RNG)
apply_window(rθϕ::AbstractArray{T}, win, rmin, rmax, win_r, win_Δr; rng=Random.GLOBAL_RNG) where {T<:Real}

The function apply_window() takes a sample of points in rθϕ and filters out points with probability specified by 1-win/maximum(win). Thus, all points are retained where win == maximum(win), and points are filtered out with proportional probability so that none are kept where win <= 0.

power_win_mix(win1, win2, wmodes, cmodes; div2Lp1=false, interchange_NN′=false)
power_win_mix(win1, win2, w̃mat, vmat, wmodes, bcmodes; div2Lp1=false, interchange_NN′=false)
power_win_mix(wmix, wmix_negm, cmodes)

This function is used to calculate the coupling matrix $\mathcal{M}_{\ell nn'}^{LNN'}$ – the first version without any binning, the second version also takes the binning matrices w̃mat and vmat to calculate the coupling matrix of the binned modes $\mathcal{N}_{LNN'}^{lnn'}$. These assume the symmetry between $N$ and $N'$.

The last version is probably not useful except for testing. It takes a fully calculated window mixing matrix to calculate the coupling matrix brute-force.

If div2Lp1=true then the whole matrix is divided by $2L+1$.

If interchange_NN′=true then calculate the same, but with $N$ and $N'$ interchanged, which might be useful for the covariance matrix.

Either version of power_win_mix() will specialize to a separable window function if win is a SeparableArray.

The basic usage is to multiply the power spectrum Clnn by this matrix, and the assuption is that there is symmetry in the exchange of k_n and k_n′. (Note that this assumed symmetry, however, destroyes the symmetry in the coupling matrix.)

window_r(wmodes::ConfigurationSpaceModes)

Get the $r$-values of the radial bins and corresponding widths $\Delta r$, e.g.,

r, Δr = SFB.window_r(wmodes)

WindowChains

This module exposes several ways to calculate a chain of window functions. All require a cache to be created that speeds up subsequent calls. The type of cache determines which method is used.

To create the cache, simply call one of

julia> cache = SFB.WindowChains.WindowChainsCacheWignerChain(win, wmodes, amodes)
julia> cache = SFB.WindowChains.WindowChainsCacheFullWmix(win, wmodes, amodes)
julia> cache = SFB.WindowChains.WindowChainsCacheFullWmixOntheflyWmix(win, wmodes, amodes)
julia> cache = SFB.WindowChains.WindowChainsCacheSeparableWmix(win, wmodes, amodes)
julia> cache = SFB.WindowChains.WindowChainsCacheSeparableWmixOntheflyWlmlm(win, wmodes, amodes)
julia> cache = SFB.WindowChainsCache(win, wmodes, amodes)

The last one will automatically select the typically fastest algorithm.

Then to calculate an element of a chain of windows, use window_chain(ell, n1, n2, cache). For example, for a chain of four windows,

julia> ell = [0, 1, 2, 3]
julia> n1 = [1, 2, 3, 4]
julia> n2 = [5, 6, 7, 8]
julia> Wk = SFB.window_chain(ell, n1, n2, cache)

window_chain(ell, n1, n2, cache, symmetries)

This version adds up several window chains taking into account the symmetries. symmetries is an array of pairs of numbers specifying the symmetries to consider. Each pair specifies the ℓnn′ index and the type of symmetry. For example, when k≥3, then symmetries = [1=>0, 2=>1, 3=>2] would specify that no symmetries are taken into account for the first lnn′ combination, the second symmetry will flip n and n′ and add the result only when n ≠ n′, and the third will add the result regardless whether the n equal or not.

calc_C4(C_th, cmix_W, cmix_wW, Veff, cmodes)

Calculate the observed power spectrum including the local average effect for a constant nbar.

The method is by default an approximate formula.

calc_NobsA(NwW_th, NW_th, cmix_wW, nbar, Veff, cmodes)

Calculate the observed shot noise including the local average effect for a constant nbar.

calc_NobsA_z(NwW_th, nbar, cmodes, amodes_red, wWmix, wWmix_negm, wWtildemix, wWtildemix_negm)

Calculate the observed shot noise including the local average effect for measured nbar(z).

MyBroadcast

This module defines the function mybroadcast. It behave similarly to a threaded broadcast, except that it tries to batch iterations such that each batch takes about 0.5 seconds to perform.

The idea is to automatically adjust the number of iterations per batch so that overhead per iteration is low and batch size is small so that the threads keep getting scheduled.

For example, imagine that the execution time per iteration increases. With a static scheduler, this would mean that the first threads finish long before the last thread. This avoids that by adjusting the number of iterations so that each batch should take approximately 0.5 seconds.

So why batch iterations? Imagine you need to allocate a buffer for each iteration, and this buffer can be shared for sequentially run iterations. Allocating a separate buffer would add a lot of overhead, so that traditional map() can take longer than the serial implementation. Batching avoids that pitfall.

mybroadcast(fn, x...; num_threads=Threads.nthreads())

MeshedArrays

A 'MeshedArray' struct is a kind of view into a smaller array 'x'. It is useful for treating a 1D array as a 2D array ignoring one of the dimensions. A 'MeshedArray' is an 'AbstractArray'.