Density-Split statistics

Note : Due to a conflict of versions in the cosmodesienvironment, it is not possible to run the densitysplit code in a notebook at the moment. However, the code is correct and can be run in a python script (given a swap of the pyrecon module to a previous version).

import numpy as np
from pathlib import Path
from astropy.io import fits
import matplotlib.pyplot as plt

from acm.estimators.galaxy_clustering.densitysplit import DensitySplit

# Load the data
data_fn = Path('data/mock_data.fits')
with fits.open(data_fn) as hdul:
    data = hdul[1].data
data_positions = np.c_[data['X'], data['Y'], data['Z']]

# Densitysplit properties 
boxsize = 1500
smoothing_radius = 20
cellsize = 10.0
nquantiles=3

# Correlation function properties
sedges = np.linspace(0, 150, 151) # 1 Mpc/h bins
muedges = np.linspace(-1, 1, 120) # 240 bins
edges = (sedges, muedges)

ds = DensitySplit(boxsize=boxsize, boxcenter=boxsize/2, cellsize=cellsize)
ds.assign_data(positions=data_positions, wrap=True, clear_previous=True)
ds.set_density_contrast(smoothing_radius=smoothing_radius)

quantiles, quantiles_idx, _ = ds.set_quantiles(nquantiles=nquantiles, query_method='randoms')

ccf_los = ds.quantile_data_correlation(
    data_positions, 
    edges=edges, 
    los='z'
)
acf_los = ds.quantile_correlation(
    edges=edges, los='z', 
)

# Plot 3 quantiles
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
colors=['blue', 'gray', 'crimson']
labels = [fr'$Q_{{{i}}}$' for i in range(3)]

for i in range(3):
    cut = quantiles[i][:, 2] < 50

    ax.scatter(quantiles[i][cut, 0], quantiles[i][cut, 1],
        s=1.0, color=colors[i], label=labels[i])

ax.set_xlabel('x [Mpc/h]', fontsize=12)
ax.set_ylabel('x [Mpc/h]', fontsize=12)
ax.set_xlim(0, boxsize)
ax.set_ylim(0, boxsize)

ax.legend(fontsize=12, markerscale=5, loc='upper left');

ds.plot_quantiles()
ds.plot_quantile_data_correlation()
ds.plot_quantile_correlation()