Advanced example: Higher-order statistics on the MICEv2 Catalog
!!Warning!! This notebook requires siginficant memory (~60GB) and CPU (~170h) resources!
[1]:
import numpy as np
from astropy.table import Table
import orpheus
from matplotlib import pyplot as plt
Create a mock catalog
Query and download data
To create the dataset head over to cosmohub and enter the following query selecting galaxies that are in regions which are ´DES-complete´.
SELECT ra_gal, dec_gal, kappa, gamma1, gamma2, z_cgal, z_desdm_mc, sdss_r_true, des_asahi_full_i_true FROM cosmohub.micecatv2_0_view WHERE des_asahi_full_i_true - 0.8 * (ATAN(1.5 * z_cgal) - 0.1489) < 24 AND (dec_gal>30.0 OR (dec_gal<30.0 AND ra_gal>30.0 AND ra_gal<60.0))
Next, download the catalog somewhere.
[2]:
path_to_cat = '/vol/euclidraid4/data/lporth/HigherOrderLensing/Mocks/MICEv2/micev2_descomplete_orpheustutorial.fits'
Create a tompgraphic setup
In this example we choose four equally populated bins. We further thin the catalog by a factor of ten.
[3]:
micecat = Table.read(path_to_cat)
[4]:
micecat.keys()
[4]:
['ra_gal',
'dec_gal',
'kappa',
'gamma1',
'gamma2',
'z_cgal',
'z_desdm_mc',
'sdss_r_true',
'des_asahi_full_i_true']
[5]:
fthin = 10
sel_ra = micecat['ra_gal']>0.
nbinsz = 4
zsort = np.argsort(micecat['z_desdm_mc'][sel_ra][::fthin])
ngal_tot = len(zsort)
[6]:
zbounds = np.zeros((nbinsz+1))
zsel = [None]*nbinsz
prevbound = 0
ras = np.zeros(ngal_tot)
decs = np.zeros(ngal_tot)
e1s = np.zeros(ngal_tot)
e2s = np.zeros(ngal_tot)
zinds = np.zeros(ngal_tot, dtype=np.int32)
for elz in range(nbinsz):
print(elz)
nextbound = int((elz+1)*0.25*ngal_tot)
zbounds[elz+1] = zsort[nextbound-1]
zsel[elz] = zsort[prevbound:nextbound]
ras[prevbound:nextbound] = micecat['ra_gal'][sel_ra][::fthin][zsel[elz]]
decs[prevbound:nextbound] = micecat['dec_gal'][sel_ra][::fthin][zsel[elz]]
e1s[prevbound:nextbound] = micecat['gamma1'][sel_ra][::fthin][zsel[elz]]
e2s[prevbound:nextbound] = micecat['gamma2'][sel_ra][::fthin][zsel[elz]]
zinds[prevbound:nextbound] = elz
prevbound = nextbound
plt.hist(micecat['z_cgal'][sel_ra][::fthin][zsel[elz]],bins=np.linspace(0,1.4,71),histtype='step',label='zbin %i'%elz)
plt.xlabel('Redshift $z$')
plt.ylabel("Number counts")
plt.show()
0
1
2
3
Allocate relevant orpheus instances
Keep in mind that when decomposing the catalogs into patches to use the other_cats parameter which enforces matching patches between the lens- and the source catalog.
[7]:
sel_lenses = zinds < 2 # Only select lenses from lowest two zbins
sel_sources = zinds > 0 # Only select sources from largest three zbins
thin_lens = 5 # Additionally thin the lens catalog
[8]:
sourcecat = orpheus.SpinTracerCatalog(
spin=2, pos1=ras[sel_sources], pos2=decs[sel_sources], tracer_1=-e1s[sel_sources], tracer_2=-e2s[sel_sources], zbins=zinds[sel_sources],
geometry='spherical', units_pos1='deg', units_pos2='deg')
lenscat = orpheus.ScalarTracerCatalog(
pos1=ras[sel_lenses][::thin_lens], pos2=decs[sel_lenses][::thin_lens], tracer=np.ones_like(ras[sel_lenses][::thin_lens]), zbins=zinds[sel_lenses][::thin_lens],
geometry='spherical', units_pos1='deg', units_pos2='deg')
print(sourcecat.ngal, lenscat.ngal)
19865593 2648746
[9]:
npatches = 100
sourcecat.topatches(npatches=npatches, patchextend_deg=3., other_cats=[lenscat], verbose=True)
Computing inner region of patches
Took 57.299 seconds
Mapping catalog to healpix grid
Took 1.593 seconds
Building index hash
Took 2.129 seconds
Building buffer around patches
100/100Took 49.146 seconds
Compute projected clustering correlation function
To debias the GNN estimator lateron we need to know \(w(\theta)\). As our mock catalog has a selection function with no angular dependence we can quickly create a corresponding random catalog and store it in a ScalarTracerCatalog instance.
[10]:
ngal_rand_target = 50*lenscat.ngal # Use 50 times more randoms than data
# Area of bounding box in MICE is subregion of octant covering pi/3 rad. We are lazy and first sample the full octant and
# then subselect galaxies that are in the interior of the footprint --> Need to oversample by A_oct/A_MICE = pi/2 / pi/3 = 1.5
_ngal_randsample = int(1.5*ngal_rand_target)
# Sample on octant
rand_ra = np.random.uniform(0, np.pi/2., _ngal_randsample)
rand_sindec = np.random.uniform(np.sin(0), np.sin(np.pi/2), _ngal_randsample)
rand_dec = np.arcsin(rand_sindec)
# Select data within MICE footprint
_lo = rand_dec<(np.pi/6.)
masked = np.logical_or(_lo*(rand_ra<(np.pi/6.)), _lo*(rand_ra>(np.pi/3.)))
# Assign zbins to random catalog based on relative nbar in zbins of lens catalog
_ngalz0 = int(np.mean(lenscat.zbins==0) * np.sum(~masked))
rand_zbins = np.zeros(np.sum(~masked), dtype=int)
rand_zbins[_ngalz0:] = 1
# Allocate orpheus instance
randcat = orpheus.ScalarTracerCatalog(
pos1=rand_ra[~masked]*180/np.pi, pos2=rand_dec[~masked]*180/np.pi, tracer=np.ones(np.sum(~masked)), zbins=rand_zbins,
geometry='spherical', units_pos1='deg', units_pos2='deg')
As a sanity check, let us make sure that the footprints of the two catalogs agree.
[11]:
plt.scatter(lenscat.pos1[::100],lenscat.pos2[::100],s=0.1, label='Data')
plt.scatter(randcat.pos1[::10000],randcat.pos2[::10000],s=0.1, label='Random')
plt.legend(loc="upper right", markerscale=10)
plt.show()
We now define the binning for the 2pcf. As this will be used for GNN we need to make sure to cover a wide range of scales. In order to not get extremely noisy results we choose a fairly coarse radial binning.
Note: There is a fairly subtle point hidden: As we have allowed the patches to only overlap by three degrees (i.e. the largest separation used for all 3PCFs) above, the same size is also chosen for the buffer region of the joint data+random catalog used to compute :math:`w(theta)`, which however is computed up until six degrees, meaning that we will only count some large-separation pairs once in the overlap. This should be what we want as this setup only counts those separations that are contained within the corresponding extended patches of the GNN correlator.
[12]:
min_sep = 0.03125
max_sep = 360.
binsize = 0.2
rmin_pixsize=10
tree_resos = [0, 0.5, 1., 2., 4]
nthreads=48
[13]:
nncorr = orpheus.NNCorrelation(min_sep=min_sep, max_sep=max_sep, binsize=binsize, verbosity=1,
tree_resos=tree_resos, rmin_pixsize=rmin_pixsize, nthreads=nthreads)
Processing the function we note that a significant amount of time is spent in the creation of the patch decomposition (specifically in allocating the buffer regions) and in the creation of the multihash. Part of this is due to the fact that these operations are not (yet) parallelised.
[14]:
%%time
# Warning, this takes ~4 CPUh!
nncorr.process(cat=lenscat, cat_random=randcat)
Doing patch 1/100
Doing patch 2/100
Doing patch 3/100
Doing patch 4/100
Doing patch 5/100
Doing patch 6/100
Doing patch 7/100
Doing patch 8/100
Doing patch 9/100
Doing patch 10/100
Doing patch 11/100
Doing patch 12/100
Doing patch 13/100
Doing patch 14/100
Doing patch 15/100
Doing patch 16/100
Doing patch 17/100
Doing patch 18/100
Doing patch 19/100
Doing patch 20/100
Doing patch 21/100
Doing patch 22/100
Doing patch 23/100
Doing patch 24/100
Doing patch 25/100
Doing patch 26/100
Doing patch 27/100
Doing patch 28/100
Doing patch 29/100
Doing patch 30/100
Doing patch 31/100
Doing patch 32/100
Doing patch 33/100
Doing patch 34/100
Doing patch 35/100
Doing patch 36/100
Doing patch 37/100
Doing patch 38/100
Doing patch 39/100
Doing patch 40/100
Doing patch 41/100
Doing patch 42/100
Doing patch 43/100
Doing patch 44/100
Doing patch 45/100
Doing patch 46/100
Doing patch 47/100
Doing patch 48/100
Doing patch 49/100
Doing patch 50/100
Doing patch 51/100
Doing patch 52/100
Doing patch 53/100
Doing patch 54/100
Doing patch 55/100
Doing patch 56/100
Doing patch 57/100
Doing patch 58/100
Doing patch 59/100
Doing patch 60/100
Doing patch 61/100
Doing patch 62/100
Doing patch 63/100
Doing patch 64/100
Doing patch 65/100
Doing patch 66/100
Doing patch 67/100
Doing patch 68/100
Doing patch 69/100
Doing patch 70/100
Doing patch 71/100
Doing patch 72/100
Doing patch 73/100
Doing patch 74/100
Doing patch 75/100
Doing patch 76/100
Doing patch 77/100
Doing patch 78/100
Doing patch 79/100
Doing patch 80/100
Doing patch 81/100
Doing patch 82/100
Doing patch 83/100
Doing patch 84/100
Doing patch 85/100
Doing patch 86/100
Doing patch 87/100
Doing patch 88/100
Doing patch 89/100
Doing patch 90/100
Doing patch 91/100
Doing patch 92/100
Doing patch 93/100
Doing patch 94/100
Doing patch 95/100
Doing patch 96/100
Doing patch 97/100
Doing patch 98/100
Doing patch 99/100
Doing patch 100/100
CPU times: user 3h 42min 25s, sys: 43min 29s, total: 4h 25min 54s
Wall time: 23min 56s
[15]:
plt.semilogx(nncorr.bin_centers_mean, nncorr.bin_centers_mean*nncorr.xi[0], 'b-', label="Z0 Z0")
plt.semilogx(nncorr.bin_centers_mean, nncorr.bin_centers_mean*nncorr.xi[1], 'g-', label="Z0 Z1")
plt.semilogx(nncorr.bin_centers_mean, nncorr.bin_centers_mean*nncorr.xi[2], 'c-', label="Z1 Z0")
plt.semilogx(nncorr.bin_centers_mean, nncorr.bin_centers_mean*nncorr.xi[3], 'r-', label="Z1 Z1")
plt.axhline(0,color='black',ls='--')
plt.xlabel(r"$\theta$ [arcmin]")
plt.ylabel(r"$\theta \ w(\theta)$")
plt.legend()
[15]:
<matplotlib.legend.Legend at 0x7c203c0f4710>
Compute GNN statistics
[16]:
min_sep = 0.25
max_sep = 180.
binsize=0.1
nmax=15
rmin_pixsize=15
tree_resos = [0,1.,2.,4.]
nthreads=48
apradii = np.geomspace(1,32,21)
[17]:
gnn = orpheus.GNNCorrelation(min_sep=min_sep, max_sep=max_sep, binsize=binsize, nmaxs=nmax, tree_resos=tree_resos, rmin_pixsize=rmin_pixsize,
verbosity=2, nthreads=nthreads)
Process the data
Processing the data works the same way as for the flat-sky case by using the .process method.
In the GNN case we can individually select whether we want to use the tomographic information in the lens- and/or source catalog.
We need to set the
rotsignflipparameter that relates to the convention on how the poles on the full-sky are defined. For simulated data I have always usedrotsignflip=False. A simple test to choose the sign forrotsignflipwould be to compute the shear2PCF and to make sure that xi- does not depend on the declination.Similar to the other statistics we can return the results from the individual patches using the
keep_patchresparameter.
[18]:
%%time
# Warning, this takes ~20 CPUh!
gnn_patches = gnn.process(cat_source=sourcecat, cat_lens=lenscat, dotomo_lens=True, dotomo_source=True, rotsignflip=False, keep_patchres=True)
Doing patch 1/100
Done 99.82 per cent
/users/lporth/anaconda3/envs/orpheus_devel/lib/python3.12/site-packages/orpheus/npcf_third.py:762: RuntimeWarning: invalid value encountered in sqrt
_patchnorm = np.nan_to_num(np.sqrt(_shelltriplets))
Doing patch 2/100
Done 99.86 per centDoing patch 3/100
Done 99.86 per centDoing patch 4/100
Done 99.86 per centDoing patch 5/100
Done 99.82 per centDoing patch 6/100
Done 99.82 per centDoing patch 7/100
Done 99.86 per centDoing patch 8/100
Done 99.86 per centDoing patch 9/100
Done 99.81 per centDoing patch 10/100
Done 99.86 per centDoing patch 11/100
Done 99.86 per centDoing patch 12/100
Done 99.86 per centDoing patch 13/100
Done 99.75 per centDoing patch 14/100
Done 99.77 per centDoing patch 15/100
Done 99.82 per centDoing patch 16/100
Done 99.86 per centDoing patch 17/100
Done 99.82 per centDoing patch 18/100
Done 99.86 per centDoing patch 19/100
Done 99.86 per centDoing patch 20/100
Done 99.84 per centDoing patch 21/100
Done 99.87 per centDoing patch 22/100
Done 99.86 per centDoing patch 23/100
Done 99.86 per centDoing patch 24/100
Done 99.81 per centDoing patch 25/100
Done 99.86 per centDoing patch 26/100
Done 99.86 per centDoing patch 27/100
Done 99.82 per centDoing patch 28/100
Done 99.86 per centDoing patch 29/100
Done 99.82 per centDoing patch 30/100
Done 99.77 per centDoing patch 31/100
Done 99.76 per centDoing patch 32/100
Done 99.87 per centDoing patch 33/100
Done 99.87 per centDoing patch 34/100
Done 99.85 per centDoing patch 35/100
Done 99.83 per centDoing patch 36/100
Done 99.86 per centDoing patch 37/100
Done 99.86 per centDoing patch 38/100
Done 99.86 per centDoing patch 39/100
Done 99.86 per centDoing patch 40/100
Done 99.86 per centDoing patch 41/100
Done 99.82 per centDoing patch 42/100
Done 99.86 per centDoing patch 43/100
Done 99.87 per centDoing patch 44/100
Done 99.81 per centDoing patch 45/100
Done 99.87 per centDoing patch 46/100
Done 99.81 per centDoing patch 47/100
Done 99.86 per centDoing patch 48/100
Done 99.87 per centDoing patch 49/100
Done 99.86 per centDoing patch 50/100
Done 99.81 per centDoing patch 51/100
Done 99.81 per centDoing patch 52/100
Done 99.86 per centDoing patch 53/100
Done 99.87 per centDoing patch 54/100
Done 99.86 per centDoing patch 55/100
Done 99.86 per centDoing patch 56/100
Done 99.86 per centDoing patch 57/100
Done 99.81 per centDoing patch 58/100
Done 99.82 per centDoing patch 59/100
Done 99.87 per centDoing patch 60/100
Done 99.82 per centDoing patch 61/100
Done 99.86 per centDoing patch 62/100
Done 99.86 per centDoing patch 63/100
Done 99.82 per centDoing patch 64/100
Done 99.86 per centDoing patch 65/100
Done 99.82 per centDoing patch 66/100
Done 99.82 per centDoing patch 67/100
Done 99.86 per centDoing patch 68/100
Done 99.81 per centDoing patch 69/100
Done 99.86 per centDoing patch 70/100
Done 99.82 per centDoing patch 71/100
Done 99.84 per centDoing patch 72/100
Done 99.86 per centDoing patch 73/100
Done 99.86 per centDoing patch 74/100
Done 99.82 per centDoing patch 75/100
Done 99.86 per centDoing patch 76/100
Done 99.80 per centDoing patch 77/100
Done 99.80 per centDoing patch 78/100
Done 99.86 per centDoing patch 79/100
Done 99.86 per centDoing patch 80/100
Done 99.75 per centDoing patch 81/100
Done 99.86 per centDoing patch 82/100
Done 99.82 per centDoing patch 83/100
Done 99.86 per centDoing patch 84/100
Done 99.86 per centDoing patch 85/100
Done 99.86 per centDoing patch 86/100
Done 99.82 per centDoing patch 87/100
Done 99.87 per centDoing patch 88/100
Done 99.82 per centDoing patch 89/100
Done 99.86 per centDoing patch 90/100
Done 99.87 per centDoing patch 91/100
Done 99.86 per centDoing patch 92/100
Done 99.85 per centDoing patch 93/100
Done 99.86 per centDoing patch 94/100
Done 99.86 per centDoing patch 95/100
Done 99.82 per centDoing patch 96/100
Done 99.86 per centDoing patch 97/100
Done 99.86 per centDoing patch 98/100
Done 99.86 per centDoing patch 99/100
Done 99.86 per centDoing patch 100/100
Done 99.86 per centCPU times: user 13h 11min 8s, sys: 47min 11s, total: 13h 58min 20s
Wall time: 23min 54s
[19]:
centers_patches_gnn, npcf_multipoles_patches_gnn, npcf_multipoles_norm_patches_gnn = gnn_patches
Now let us collect the measurements, compute MNN on the individual patches, and plot the results. (Note that in case you are solely interested in the MNN of the full footprint, the loop over the patches is not required). To assess the relevance of the clusering-induced correction we compute the MNN statistics for both, the corrected and for the uncorrected 3PCFs.
[20]:
%%time
import sys
# Create a copy of our GNN instance which will be used to compute MNN without the clustering correction applied
gnn_uncorr = orpheus.GNNCorrelation(
min_sep=gnn.min_sep,
max_sep=gnn.max_sep,
nbinsr=gnn.nbinsr,
nmaxs=gnn.nmaxs,
nbinsphi=gnn.nbinsphi,
nthreads=gnn.nthreads)
gnn_uncorr.nbinsz_lens = gnn.nbinsz_lens
gnn_uncorr.nbinsz_source = gnn.nbinsz_source
gnn_uncorr.projection = 'X'
gnn_uncorr.bin_centers_mean = 1.*gnn.bin_centers_mean
gnn_uncorr.npcf_multipoles = 1.*gnn.npcf_multipoles
gnn_uncorr.npcf_multipoles_norm = 1.*gnn.npcf_multipoles_norm
# Compute MNN of full 3pcf
fskyMNN = gnn.computeNNM(radii=apradii, xi=[nncorr.bin_centers_mean, nncorr.xi])
fskyMNN_uncorr = gnn_uncorr.computeNNM(radii=apradii)
# Compute MNN of patchwise 3pcfs
allMNN = np.zeros((sourcecat.npatches, *fskyMNN.shape), dtype=fskyMNN.dtype)
allMNN_uncorr = np.zeros((sourcecat.npatches, *fskyMNN_uncorr.shape), dtype=fskyMNN_uncorr.dtype)
for elp in range(sourcecat.npatches):
sys.stdout.write('%i '%elp)
_pinst = orpheus.GNNCorrelation(
min_sep=gnn.min_sep,
max_sep=gnn.max_sep,
nbinsr=gnn.nbinsr,
nmaxs=gnn.nmaxs,
nbinsphi=gnn.nbinsphi,
nthreads=gnn.nthreads)
_pinst.nbinsz_lens = gnn.nbinsz_lens
_pinst.nbinsz_source = gnn.nbinsz_source
_pinst.projection = 'X'
_pinst.bin_centers_mean = np.mean(centers_patches_gnn[elp],axis=(0,1))
_pinst.npcf_multipoles = npcf_multipoles_patches_gnn[elp]
_pinst.npcf_multipoles_norm = npcf_multipoles_norm_patches_gnn[elp]
_pmnn = _pinst.computeNNM(radii=apradii)
allMNN_uncorr[elp] += 1*_pmnn
_pinst.npcf = None # Wipe npcf from instance as this will trigger computeNNM to reallocate it in the next line
_pmnn = _pinst.computeNNM(radii=apradii, xi=[nncorr.bin_centers_mean, nncorr.xi])
allMNN[elp] += 1*_pmnn
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 CPU times: user 20min 49s, sys: 2min 2s, total: 22min 52s
Wall time: 10min 12s
[21]:
# Estimate of the sample variance for data with different weights. We choose as effective weights the sum over the n===0 multipole
MNN_effws = np.sum(npcf_multipoles_norm_patches_gnn[:,0], axis=(-1,-2)).real
MNN_effws /= np.mean(MNN_effws,axis=0)
_V1 = np.sum(MNN_effws,axis=0)
_V2 = np.sum(MNN_effws**2,axis=0)
MMN_samplestd_re = np.sqrt(np.sum(MNN_effws.real[...,np.newaxis]*(allMNN.real[:,0]-fskyMNN.real)**2,axis=0)/(_V1-_V2/_V1**2)[:,np.newaxis])
MMN_samplestd_im = np.sqrt(np.sum(MNN_effws.real[...,np.newaxis]*(allMNN.imag[:,0]-fskyMNN.imag)**2,axis=0)/(_V1-_V2/_V1**2)[:,np.newaxis])
MMN_uncorr_samplestd_re = np.sqrt(np.sum(MNN_effws.real[...,np.newaxis]*(allMNN_uncorr.real[:,0]-fskyMNN_uncorr.real)**2,axis=0)/(_V1-_V2/_V1**2)[:,np.newaxis])
MMN_uncorr_samplestd_im = np.sqrt(np.sum(MNN_effws.real[...,np.newaxis]*(allMNN_uncorr.imag[:,0]-fskyMNN_uncorr.imag)**2,axis=0)/(_V1-_V2/_V1**2)[:,np.newaxis])
for sbin in range(sourcecat.nbinsz):
for lbin1 in range(lenscat.nbinsz):
for lbin2 in range(lbin1,lenscat.nbinsz):
elzcombi = sbin*lenscat.nbinsz**2 + lbin1*lenscat.nbinsz + lbin2
plt.plot(MNN_effws[:,elzcombi],label="ZS%i ZL%i ZL%i"%(2+sbin,1+lbin1,1+lbin2))
plt.xlabel('Patch index')
plt.ylabel('Relative weight of patch')
plt.legend()
plt.show()
[22]:
# Plot the result. Solid lines are the corrected MNN while dashed lines are the uncorrected ones. We
# rescale the std by 0.1, approximating the error-on-the-mean statistics.
for sbin in range(sourcecat.nbinsz):
for lbin1 in range(lenscat.nbinsz):
for lbin2 in range(lbin1,lenscat.nbinsz):
elzcombi = sbin*lenscat.nbinsz**2 + lbin1*lenscat.nbinsz + lbin2
plt.axhline(0,color='black',ls='--',lw=2)
plt.errorbar(x=apradii,
y=apradii*fskyMNN[0,elzcombi].real,
yerr=.1*apradii*MMN_samplestd_re[elzcombi],color='blue',lw=3,
label=r'$\left\langle\mathcal{M}_{\rm ap} \, \mathcal{N}_{\rm ap} \, \mathcal{N}_{\rm ap}\right\rangle$')
plt.errorbar(x=apradii,
y=apradii*fskyMNN[0,elzcombi].imag,
yerr=.1*apradii*MMN_samplestd_im[elzcombi],color='red',lw=3,
label=r'$\left\langle\mathcal{M}_{\times} \, \mathcal{N}_{\rm ap} \, \mathcal{N}_{\rm ap}\right\rangle$')
plt.errorbar(x=apradii,
y=apradii*fskyMNN_uncorr[0,elzcombi].real,
yerr=.1*apradii*MMN_uncorr_samplestd_re[elzcombi],color='blue',ls='--',lw=3)
plt.errorbar(x=apradii,
y=apradii*fskyMNN_uncorr[0,elzcombi].imag,
yerr=.1*apradii*MMN_uncorr_samplestd_im[elzcombi],color='red',ls='--',lw=3)
for elp in range(sourcecat.npatches):
plt.semilogx(apradii, apradii*allMNN[elp,0,elzcombi],color='blue',alpha=0.1)
plt.xlabel('Aperture radius [arcmin]')
plt.ylabel(r'$\theta \ \times \ \mathcal{MNN}(z_%i,z_%i,z_%i)$'%(2+sbin,1+lbin1,1+lbin2))
plt.ylim(-0.2*np.max(apradii*fskyMNN[0,elzcombi].real),1.2*np.max(apradii*fskyMNN[0,elzcombi].real))
plt.xscale('log')
plt.legend(loc='upper right')
plt.show()
/users/lporth/anaconda3/envs/orpheus_devel/lib/python3.12/site-packages/matplotlib/cbook.py:1709: ComplexWarning: Casting complex values to real discards the imaginary part
return math.isfinite(val)
/users/lporth/anaconda3/envs/orpheus_devel/lib/python3.12/site-packages/matplotlib/cbook.py:1345: ComplexWarning: Casting complex values to real discards the imaginary part
return np.asarray(x, float)
Compute NGG statistics
[23]:
min_sep = 0.25
max_sep = 180.
binsize=0.1
nmax=15
rmin_pixsize=40 # We choose a large value here as a small one might produce an artificial B-Mode
tree_resos = [0.,1.,2.]
nthreads=48
[24]:
ngg = orpheus.NGGCorrelation(min_sep=min_sep, max_sep=max_sep, binsize=binsize, nmaxs=nmax,
tree_resos=tree_resos, rmin_pixsize=rmin_pixsize, verbosity=2, nthreads=nthreads)
[25]:
%%time
# Warning, this takes ~50 CPUh!
ngg_patches = ngg.process(cat_source=sourcecat, cat_lens=lenscat, dotomo_lens=True, dotomo_source=True, rotsignflip=False, keep_patchres=True)
Doing patch 1/100
.........|.........|.........|.........|.........|.........|.........|.........|.........|.........|Doing patch 2/100
..........................................................................................Doing patch 3/100
..........................................................................................Doing patch 4/100
..........................................................................................Doing patch 5/100
..........................................................................................Doing patch 6/100
..........................................................................................Doing patch 7/100
..........................................................................................Doing patch 8/100
..........................................................................................Doing patch 9/100
..........................................................................................Doing patch 10/100
..........................................................................................Doing patch 11/100
..........................................................................................Doing patch 12/100
..........................................................................................Doing patch 13/100
..........................................................................................Doing patch 14/100
..........................................................................................Doing patch 15/100
..........................................................................................Doing patch 16/100
..........................................................................................Doing patch 17/100
..........................................................................................Doing patch 18/100
..........................................................................................Doing patch 19/100
..........................................................................................Doing patch 20/100
..........................................................................................Doing patch 21/100
..........................................................................................Doing patch 22/100
..........................................................................................Doing patch 23/100
..........................................................................................Doing patch 24/100
..........................................................................................Doing patch 25/100
..........................................................................................Doing patch 26/100
..........................................................................................Doing patch 27/100
..........................................................................................Doing patch 28/100
..........................................................................................Doing patch 29/100
..........................................................................................Doing patch 30/100
..........................................................................................Doing patch 31/100
..........................................................................................Doing patch 32/100
..........................................................................................Doing patch 33/100
..........................................................................................Doing patch 34/100
..........................................................................................Doing patch 35/100
..........................................................................................Doing patch 36/100
..........................................................................................Doing patch 37/100
..........................................................................................Doing patch 38/100
..........................................................................................Doing patch 39/100
..........................................................................................Doing patch 40/100
..........................................................................................Doing patch 41/100
..........................................................................................Doing patch 42/100
..........................................................................................Doing patch 43/100
..........................................................................................Doing patch 44/100
..........................................................................................Doing patch 45/100
..........................................................................................Doing patch 46/100
..........................................................................................Doing patch 47/100
..........................................................................................Doing patch 48/100
..........................................................................................Doing patch 49/100
..........................................................................................Doing patch 50/100
..........................................................................................Doing patch 51/100
..........................................................................................Doing patch 52/100
..........................................................................................Doing patch 53/100
..........................................................................................Doing patch 54/100
..........................................................................................Doing patch 55/100
..........................................................................................Doing patch 56/100
..........................................................................................Doing patch 57/100
..........................................................................................Doing patch 58/100
..........................................................................................Doing patch 59/100
..........................................................................................Doing patch 60/100
..........................................................................................Doing patch 61/100
..........................................................................................Doing patch 62/100
..........................................................................................Doing patch 63/100
..........................................................................................Doing patch 64/100
..........................................................................................Doing patch 65/100
..........................................................................................Doing patch 66/100
..........................................................................................Doing patch 67/100
..........................................................................................Doing patch 68/100
..........................................................................................Doing patch 69/100
..........................................................................................Doing patch 70/100
..........................................................................................Doing patch 71/100
..........................................................................................Doing patch 72/100
..........................................................................................Doing patch 73/100
..........................................................................................Doing patch 74/100
..........................................................................................Doing patch 75/100
..........................................................................................Doing patch 76/100
..........................................................................................Doing patch 77/100
..........................................................................................Doing patch 78/100
..........................................................................................Doing patch 79/100
..........................................................................................Doing patch 80/100
..........................................................................................Doing patch 81/100
..........................................................................................Doing patch 82/100
..........................................................................................Doing patch 83/100
..........................................................................................Doing patch 84/100
..........................................................................................Doing patch 85/100
..........................................................................................Doing patch 86/100
..........................................................................................Doing patch 87/100
..........................................................................................Doing patch 88/100
..........................................................................................Doing patch 89/100
..........................................................................................Doing patch 90/100
..........................................................................................Doing patch 91/100
..........................................................................................Doing patch 92/100
..........................................................................................Doing patch 93/100
..........................................................................................Doing patch 94/100
..........................................................................................Doing patch 95/100
..........................................................................................Doing patch 96/100
..........................................................................................Doing patch 97/100
..........................................................................................Doing patch 98/100
..........................................................................................Doing patch 99/100
..........................................................................................Doing patch 100/100
..........................................................................................CPU times: user 1d 21h 22min 1s, sys: 1h 59min 7s, total: 1d 23h 21min 9s
Wall time: 1h 58s
[26]:
centers_patches_ngg, npcf_multipoles_patches_ngg, npcf_multipoles_norm_patches_ngg = ngg_patches
[27]:
import sys
apradii = np.geomspace(1,32,21)
# Compute Map3 of full 3pcf
fskyNMM = ngg.computeNMM(radii=apradii)
# Compute Map3 of patchwise 3pcfs
allNMM = np.zeros((sourcecat.npatches, *fskyNMM.shape), dtype=fskyNMM.dtype)
for elp in range(sourcecat.npatches):
sys.stdout.write('%i '%elp)
_pinst = orpheus.NGGCorrelation(
min_sep=ngg.min_sep,
max_sep=ngg.max_sep,
nbinsr=ngg.nbinsr,
nmaxs=ngg.nmaxs,
nbinsphi=ngg.nbinsphi,
nthreads=ngg.nthreads)
_pinst.nbinsz_lens = ngg.nbinsz_lens
_pinst.nbinsz_source = ngg.nbinsz_source
_pinst.projection = 'X'
_pinst.bin_centers_mean = np.mean(centers_patches_ngg[elp],axis=(0,1))
_pinst.npcf_multipoles = npcf_multipoles_patches_ngg[elp]
_pinst.npcf_multipoles_norm = npcf_multipoles_norm_patches_ngg[elp]
#if elp==0:
# filtercache = []
# for elr, r in enumerate(apradii):
# filtercache.append(_pinst._NMM_filtergrid(r, r, r))
_pnmm = _pinst.computeNMM(radii=apradii)
allNMM[elp] += 1*_pnmm
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
[28]:
# Choose as effective weights the sum over the n===0 multipole
NMM_effws = np.sum(npcf_multipoles_norm_patches_ngg[:,ngg.nmaxs[0]], axis=(-1,-2)).real
NMM_effws /= np.mean(NMM_effws,axis=0)
NMM_stds_re = np.std(allNMM.real[:,0],axis=0)
NMM_stds_im = np.std(allNMM.imag[:,0],axis=0)
_V1 = np.sum(NMM_effws,axis=0)
_V2 = np.sum(NMM_effws**2,axis=0)
NMM_samplestd_re = np.sqrt(np.sum(NMM_effws.real[:,np.newaxis,:,np.newaxis]*(allNMM.real-fskyNMM.real)**2,axis=0)/(_V1-_V2/_V1**2)[:,np.newaxis])
NMM_samplestd_im = np.sqrt(np.sum(NMM_effws.real[:,np.newaxis,:,np.newaxis]*(allNMM.imag-fskyNMM.imag)**2,axis=0)/(_V1-_V2/_V1**2)[:,np.newaxis])
[29]:
for lbin in range(lenscat.nbinsz):
for sbin1 in range(sourcecat.nbinsz):
for sbin2 in range(sbin1, sourcecat.nbinsz):
elzcombi = lbin*sourcecat.nbinsz**2 + sbin1*sourcecat.nbinsz + sbin2
plt.axhline(0,color='black',ls='--',lw=2)
plt.errorbar(x=apradii,
y=apradii*fskyNMM[0,elzcombi].real,
yerr=apradii*NMM_samplestd_re[0,elzcombi],color='blue',lw=3,
label=r'$\left\langle\mathcal{N}_{\rm ap} \mathcal{M}_{\rm ap} \, \mathcal{M}_{\rm ap}\right\rangle$')
plt.errorbar(x=apradii,
y=apradii*fskyNMM[1,elzcombi].real,
yerr=apradii*NMM_samplestd_re[1,elzcombi],color='red',lw=3,
label=r'$\left\langle\mathcal{N}_{\rm ap} \mathcal{M}_{\times} \, \mathcal{M}_{\times}\right\rangle$')
plt.errorbar(x=apradii,
y=apradii*fskyNMM[2,elzcombi].real,
yerr=apradii*NMM_samplestd_re[2,elzcombi],color='green',lw=3,
label=r'$\left\langle\mathcal{N}_{\rm ap} \mathcal{M}_{\times} \, \mathcal{M}_{\rm ap}\right\rangle$')
for elp in range(sourcecat.npatches):
plt.semilogx(apradii, apradii*allNMM[elp,0,elzcombi],color='blue',alpha=0.1)
plt.xlabel('Aperture radius [arcmin]')
plt.ylabel(r'$\theta \ \times \ \mathcal{NMM}(z_%i,z_%i,z_%i)$'%(1+lbin,2+sbin1,2+sbin2))
plt.ylim(-0.2*np.max(apradii*fskyNMM[0,elzcombi].real),1.2*np.max(apradii*fskyNMM[0,elzcombi].real))
plt.legend(loc='upper right')
plt.show()
/users/lporth/anaconda3/envs/orpheus_devel/lib/python3.12/site-packages/matplotlib/cbook.py:1709: ComplexWarning: Casting complex values to real discards the imaginary part
return math.isfinite(val)
/users/lporth/anaconda3/envs/orpheus_devel/lib/python3.12/site-packages/matplotlib/cbook.py:1345: ComplexWarning: Casting complex values to real discards the imaginary part
return np.asarray(x, float)
Compute GGG statistics
[30]:
min_sep = 0.25
max_sep = 180.
binsize=0.1
nmax=15
rmin_pixsize=15
tree_resos = [0,1.,2.,4.]
nthreads=48
[31]:
ggg = orpheus.GGGCorrelation(n_cfs=4,min_sep=min_sep, max_sep=max_sep, binsize=binsize, tree_resos=tree_resos, rmin_pixsize=rmin_pixsize,
verbosity=2, nthreads=nthreads)
[32]:
%%time
# Warning, this takes ~100 CPUh!
ggg_patches = ggg.process(cat=sourcecat, dotomo=True, rotsignflip=False, keep_patchres=True)
Doing patch 1/100
..........................................................................................
Doing patch 2/100
..........................................................................................
Doing patch 3/100
...........................................................................................
Doing patch 4/100
...........................................................................................
Doing patch 5/100
...........................................................................................
Doing patch 6/100
..........................................................................................
Doing patch 7/100
...........................................................................................
Doing patch 8/100
..........................................................................................
Doing patch 9/100
..........................................................................................
Doing patch 10/100
..........................................................................................
Doing patch 11/100
...........................................................................................
Doing patch 12/100
..........................................................................................
Doing patch 13/100
...........................................................................................
Doing patch 14/100
..........................................................................................
Doing patch 15/100
..........................................................................................
Doing patch 16/100
...........................................................................................
Doing patch 17/100
..........................................................................................
Doing patch 18/100
...........................................................................................
Doing patch 19/100
..........................................................................................
Doing patch 20/100
..........................................................................................
Doing patch 21/100
..........................................................................................
Doing patch 22/100
..........................................................................................
Doing patch 23/100
..........................................................................................
Doing patch 24/100
..........................................................................................
Doing patch 25/100
..........................................................................................
Doing patch 26/100
..........................................................................................
Doing patch 27/100
..........................................................................................
Doing patch 28/100
..........................................................................................
Doing patch 29/100
...........................................................................................
Doing patch 30/100
..........................................................................................
Doing patch 31/100
..........................................................................................
Doing patch 32/100
..........................................................................................
Doing patch 33/100
..........................................................................................
Doing patch 34/100
..........................................................................................
Doing patch 35/100
..........................................................................................
Doing patch 36/100
..........................................................................................
Doing patch 37/100
..........................................................................................
Doing patch 38/100
..........................................................................................
Doing patch 39/100
..........................................................................................
Doing patch 40/100
...........................................................................................
Doing patch 41/100
..........................................................................................
Doing patch 42/100
..........................................................................................
Doing patch 43/100
..........................................................................................
Doing patch 44/100
..........................................................................................
Doing patch 45/100
..........................................................................................
Doing patch 46/100
..........................................................................................
Doing patch 47/100
...........................................................................................
Doing patch 48/100
..........................................................................................
Doing patch 49/100
..........................................................................................
Doing patch 50/100
..........................................................................................
Doing patch 51/100
..........................................................................................
Doing patch 52/100
..........................................................................................
Doing patch 53/100
..........................................................................................
Doing patch 54/100
...........................................................................................
Doing patch 55/100
..........................................................................................
Doing patch 56/100
..........................................................................................
Doing patch 57/100
..........................................................................................
Doing patch 58/100
...........................................................................................
Doing patch 59/100
..........................................................................................
Doing patch 60/100
...........................................................................................
Doing patch 61/100
..........................................................................................
Doing patch 62/100
..........................................................................................
Doing patch 63/100
..........................................................................................
Doing patch 64/100
...........................................................................................
Doing patch 65/100
..........................................................................................
Doing patch 66/100
..........................................................................................
Doing patch 67/100
..........................................................................................
Doing patch 68/100
..........................................................................................
Doing patch 69/100
..........................................................................................
Doing patch 70/100
..........................................................................................
Doing patch 71/100
..........................................................................................
Doing patch 72/100
..........................................................................................
Doing patch 73/100
...........................................................................................
Doing patch 74/100
...........................................................................................
Doing patch 75/100
..........................................................................................
Doing patch 76/100
..........................................................................................
Doing patch 77/100
..........................................................................................
Doing patch 78/100
..........................................................................................
Doing patch 79/100
..........................................................................................
Doing patch 80/100
..........................................................................................
Doing patch 81/100
..........................................................................................
Doing patch 82/100
..........................................................................................
Doing patch 83/100
..........................................................................................
Doing patch 84/100
...........................................................................................
Doing patch 85/100
..........................................................................................
Doing patch 86/100
..........................................................................................
Doing patch 87/100
..........................................................................................
Doing patch 88/100
..........................................................................................
Doing patch 89/100
..........................................................................................
Doing patch 90/100
..........................................................................................
Doing patch 91/100
..........................................................................................
Doing patch 92/100
..........................................................................................
Doing patch 93/100
..........................................................................................
Doing patch 94/100
...........................................................................................
Doing patch 95/100
..........................................................................................
Doing patch 96/100
..........................................................................................
Doing patch 97/100
..........................................................................................
Doing patch 98/100
..........................................................................................
Doing patch 99/100
..........................................................................................
Doing patch 100/100
...........................................................................................
CPU times: user 4d 4h 34min 35s, sys: 4h 12min 17s, total: 4d 8h 46min 52s
Wall time: 2h 15min 24s
[33]:
centers_patches_ggg, npcf_multipoles_patches_ggg, npcf_multipoles_norm_patches_ggg = ggg_patches
[34]:
%%time
apradii = np.geomspace(1,32,21)
filtercache_ggg = [None]*len(apradii)
for elapr, apr in enumerate(apradii):
filtercache_ggg[elapr] = ggg._map3_filtergrid_singleR(apr,apr,apr)
# Compute Map3 of full 3pcf
fskyMMM = ggg.computeMap3(radii=apradii, filtercache=filtercache_ggg)
# Compute Map3 of patchwise 3pcfs
allMMM = np.zeros((sourcecat.npatches, *fskyMMM.shape), dtype=fskyMMM.dtype)
for elp in range(sourcecat.npatches):
sys.stdout.write('%i '%elp)
_pinst = orpheus.GGGCorrelation(n_cfs=4,
min_sep=ggg.min_sep,
max_sep=ggg.max_sep,
nbinsr=ggg.nbinsr,
nmaxs=ggg.nmaxs,
nbinsphi=ggg.nbinsphi,
nthreads=ggg.nthreads)
_pinst.nbinsz = ggg.nbinsz
_pinst.nzcombis = ggg.nzcombis
_pinst.projection = 'X'
_pinst.bin_centers_mean = np.mean(centers_patches_ggg[elp],axis=(0))
_pinst.npcf_multipoles = npcf_multipoles_patches_ggg[elp]
_pinst.npcf_multipoles_norm = npcf_multipoles_norm_patches_ggg[elp]
_pmmm = _pinst.computeMap3(radii=apradii, filtercache=filtercache_ggg)
del _pinst
allMMM[elp] += 1*_pmmm
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 CPU times: user 29min 52s, sys: 15min 20s, total: 45min 12s
Wall time: 23min 28s
[35]:
# Choose as effective weights the sum over the n===0 multipole
MMM_effws = np.sum(npcf_multipoles_norm_patches_ggg[:,0], axis=(-1,-2)).real
MMM_effws /= np.mean(MMM_effws,axis=0)
MMM_stds_re = np.std(allMMM.real[:,0],axis=0)
MMM_stds_im = np.std(allMMM.imag[:,0],axis=0)
_V1 = np.sum(MMM_effws,axis=0)
_V2 = np.sum(MMM_effws**2,axis=0)
MMM_samplestd_re = np.sqrt(np.sum(MMM_effws.real[:,np.newaxis,:,np.newaxis]*(allMMM.real-fskyMMM.real)**2,axis=0)/(_V1-_V2/_V1**2)[:,np.newaxis])
MMM_samplestd_im = np.sqrt(np.sum(MMM_effws.real[:,np.newaxis,:,np.newaxis]*(allMMM.imag-fskyMMM.imag)**2,axis=0)/(_V1-_V2/_V1**2)[:,np.newaxis])
[36]:
for sbin1 in range(sourcecat.nbinsz):
for sbin2 in range(sbin1,sourcecat.nbinsz):
for sbin3 in range(sbin2, sourcecat.nbinsz):
elzcombi = sbin1*sourcecat.nbinsz**2 + sbin2*sourcecat.nbinsz + sbin3
plt.axhline(0,color='black',ls='--',lw=2)
plt.errorbar(x=apradii,
y=apradii*fskyMMM[0,elzcombi].real,
yerr=apradii*MMM_samplestd_re[0,elzcombi],color='blue',lw=3,
label=r'$\left\langle\mathcal{M}_{\rm ap} \mathcal{M}_{\rm ap} \, \mathcal{M}_{\rm ap}\right\rangle$')
plt.errorbar(x=apradii,
y=apradii*fskyMMM[1,elzcombi].real,
yerr=apradii*MMM_samplestd_re[1,elzcombi],color='red',lw=3,
label=r'$\left\langle\mathcal{M}_{\rm ap} \mathcal{M}_{\rm ap} \, \mathcal{M}_{\times}\right\rangle$')
plt.errorbar(x=apradii,
y=apradii*fskyMMM[2,elzcombi].real,
yerr=apradii*MMM_samplestd_re[2,elzcombi],color='red',lw=3)
plt.errorbar(x=apradii,
y=apradii*fskyMMM[3,elzcombi].real,
yerr=apradii*MMM_samplestd_re[3,elzcombi],color='red',lw=3)
plt.errorbar(x=apradii,
y=apradii*fskyMMM[4,elzcombi].real,
yerr=apradii*MMM_samplestd_re[4,elzcombi],color='green',lw=3,
label=r'$\left\langle\mathcal{M}_{\rm ap} \mathcal{M}_{\times} \, \mathcal{M}_{\times}\right\rangle$')
plt.errorbar(x=apradii,
y=apradii*fskyMMM[5,elzcombi].real,
yerr=apradii*MMM_samplestd_re[5,elzcombi],color='green',lw=3)
plt.errorbar(x=apradii,
y=apradii*fskyMMM[6,elzcombi].real,
yerr=apradii*MMM_samplestd_re[6,elzcombi],color='green',lw=3)
plt.errorbar(x=apradii,
y=apradii*fskyMMM[7,elzcombi].real,
yerr=apradii*MMM_samplestd_re[7,elzcombi],color='cyan',lw=3,
label=r'$\left\langle\mathcal{M}_{\times} \mathcal{M}_{\times} \, \mathcal{M}_{\times}\right\rangle$')
for elp in range(sourcecat.npatches):
plt.semilogx(apradii, apradii*allMMM[elp,0,elzcombi],color='blue',alpha=0.1)
plt.xlabel('Aperture radius [arcmin]')
plt.ylabel(r'$\theta \ \times \ \mathcal{MMM}(z_%i,z_%i,z_%i)$'%(2+sbin1,2+sbin2,2+sbin3))
plt.ylim(-0.2*np.max(apradii*fskyMMM[0,elzcombi].real),1.2*np.max(apradii*fskyMMM[0,elzcombi].real))
plt.legend(loc='upper right')
plt.show()
/users/lporth/anaconda3/envs/orpheus_devel/lib/python3.12/site-packages/matplotlib/cbook.py:1709: ComplexWarning: Casting complex values to real discards the imaginary part
return math.isfinite(val)
/users/lporth/anaconda3/envs/orpheus_devel/lib/python3.12/site-packages/matplotlib/cbook.py:1345: ComplexWarning: Casting complex values to real discards the imaginary part
return np.asarray(x, float)
[ ]: