Basic Usage

Points and PointContainers

The main objects in CatAna are points, coordinates in 3d space stored in spherical components. They can either be directly obtained from python arrays or read from files using the catana.io module.

Points are collected in point containers, which can then be given to a spherical Fourier-Bessel transform routine to analyze.

import numpy as np
import catana

# Some random particle positions in a sphere with radius 100
python_particles = np.random.uniform(-100,100,(400000,3))
python_particles = python_particles[np.sum(python_particles**2, axis=1) <= 100**2]

# Load into PointContainer
point_container = catana.PointContainer(python_particles, 'cartesian')

# If we want to retreive a numpy array from the point_container again
# (note that the coordinates will be spherical, with columns [r, theta, phi])
python_particles_spherical = np.array(point_container)

Computing the SFB transform will be faster if the points are binned in angular pixels. CatAna uses the HEALPix pixelization scheme. The resolution of the angular grid is given by the nside parameter, which determines how many times each base pixel is subdivided and has to be a power of 2.

CatAna contains a specialized container, the PixelizedPointContainer. It can be used as follows:

# Create a PixelizedPointContainer with nside 64
pixelized_point_container = catana.PixelizedPointContainer(64, point_container)

# We can create a HEALPix map with each pixel containing the number of points it contains
import healpy as hp
countmap = pixelized_point_container.get_countmap()
hp.mollview(countmap)
_images/countmap.png

Computing the SFB transform

The SFB decomposition can either be directly done on PointContainers (brute-force) or on PixelizedObjectContainers which is much faster.

# Parameters
params = {
    'lmax': 100,  # maximum multipole to compute
    'nmax': 50,  # maximum k-index to compute
    'rmax': 100,  # largest distance to point (field is supposed to be 0 outside this radius)
    'store_flmn': False,  # we're only interested in the C_ln components here
    'verbose': False
}

# brute-force from PointContainer (not recommended, will take a while)
%time kclkk_1 = catana.sfb_decomposition(point_container, **params)

# using the pixelized container, much faster
%time kclkk_2 = catana.sfb_decomposition(pixelized_point_container, **params)

On my machine, the recorded times are something like

CPU times: user 40min 45s, sys: 7.53 s, total: 40min 53s
Wall time: 5min 29s

CPU times: user 45.4 s, sys: 1.55 s, total: 46.9 s
Wall time: 8.93 s

As you can see, the speed-up from the fast method is big. We used a very low pixelization resolution (nside), hence the errors in the \(f_{lmn}\) coefficients and even in the \(C_{ln}\) coefficients are quite large. Here is from the previous example:

fig, axes = plt.subplots(2, 1, figsize=(6,6))
axes[1].axhline(0, color='black')
for l in [1, 10, 20]:
    g, = axes[0].plot(kclkk_1.k_ln[l], kclkk_1.c_ln[l], label='l={} (brute-force)'.format(l))
    axes[0].plot(kclkk_2.k_ln[l], kclkk_2.c_ln[l], label='l={} (pixel)'.format(l),
        color=g.get_color(), linestyle='dashed')
    axes[1].plot(kclkk_1.k_ln[l], kclkk_2.c_ln[l]/kclkk_1.c_ln[l]-1,
        label='l={}', color=g.get_color())
axes[0].set(ylabel=r'$C_l(k,k)$', xticklabels=[])
axes[1].set(xlabel=r'$k$', ylabel=r'error', yticks=[-0.1, -0.05, 0, 0.05, 0.1])
axes[0].legend()
fig.tight_layout()
_images/pixel_accuracy.png

By choosing a larger NSide, the errors will drop significantly while we can still profit from a significant speed increase over the brute-force method.

Normalization and shot noise

The \(f_{lmn}\) components are computed as

\[f_{lmn} = \sqrt{\frac{2}{\pi}} \frac{V}{N} \sum_p k j_l(k r_p) Y_{lm}^{*}(\hat{r}_p),\]

where \(V=\frac{4}{3} \pi r_{max}^3\) is the volume of the supporting tophat with radius \(r_{max}\). In case the survey window function is not a tophat (e.g. a redshift depending window function or an angular mask) you need to normalize the SFB decomposition results by multiplying with the fraction that the true volume occupies in the tophat for \(f_{lmn}\) and by the square of it for \(C_{ln}\).

Due to the discrete nature of our point objects, we add an additional term to the spherical power spectrum of the underlying smooth distribution (\(C_l(k,k)\)), the so-called shot-noise (\(C_l^{sn}(k,k)\)).

\[ \begin{align}\begin{aligned}C_l^d(k, k) = C_l(k, k) + C_l^{sn}(k,k)\\C_l^{sn}(k,k) = \frac{V}{N} \frac{2}{\pi} k^2 \int_0^\infty dr \; r^2 \; w(r) \; j_l(k r)^2\end{aligned}\end{align} \]

where w(r) is the radial window function. This contribution has to be either computed numerically from a random catalog (basically our example above) or analytically. CatAna provides a function to compute these spherical Bessel function integrals, see double_sbessel_integrator.

Reading Data

The most convenient way to use CatAna in Python is to directly use data from numpy arrays as shown above. CatAna can however also deal with files in (space delimited) text files and Gadget files. For detailed documentation take a look at catana.io.

# Save data so that we can read it again (uses [r theta phi] columns)
sink = catana.io.SphericalTextSink("points.txt")
sink.write(point_container)

# Read data from a text file
source = catana.io.SphericalTextSource("points.txt")
point_container = source.get_point_container()

# Other supported text sources and sinks:
#     - CartesianTextSource / CartesianTextSink
#     - SphericalTextSource / SphericalTextSink
#     - SphericalTextSource_3dex / SphericalTextSink_3dex (uses [theta phi r] columns)
#     - GadgetSource

Filtering Data

Masking and removing points prior to decomposition can be done directly on the data array. CatAna provides some additional functionality to filter PointContainers which are illustrated below.

# Gaussian Filter with scale factor 20
gaussian_filter = catana.io.GaussianRadialWindowFunctionFilter(50)

# Apply to our point_container
gaussian_filter(point_container_2)

# Plot point radii histogram normalized by the surface at the given radii
radii = np.array(point_container_2)[:,0]
fig, ax = plt.subplots(1, 1, figsize=(6,4))
ax.hist(radii, weights=1/radii**2, bins=20, normed=True)
ax.set(yticks=[], xlabel='r')
_images/gaussian_filter_hist.png 
# Angular mask from HEALPix map
mask = np.zeros(12)
mask[3] = 1
filt_mask = catana.io.AngularMaskFilter(mask)
filt_mask(point_container_2)

pixelized_point_container_2 = catana.PixelizedPointContainer(64, point_container_2)
hp.mollview(pixelized_point_container_2.get_countmap())
_images/countmap_angmask.png

Note

For more filters, take a look at Filters.