<jupyter_start><jupyter_text>API demonstration for paper of v1.0_the LSST-DESC CLMM team_Here we demonstrate how to use `clmm` to estimate a WL halo mass from observations of a galaxy cluster when source galaxies follow a given distribution (The LSST DESC Science Requirements Document - arXiv:1809.01669,  implemented in `clmm`). It uses several functionalities of the support `mock_data` module to produce mock datasets.- Setting things up, with the proper imports.- Computing the binned reduced tangential shear profile, for the 2 datasets, using logarithmic binning.- Setting up a model accounting for the redshift distribution.- Perform a simple fit using `scipy.optimize.curve_fit` included in `clmm` and visualize the results.  Setup First, we import some standard packages.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
plt.rcParams['font.family'] = ['gothambook','gotham','gotham-book','serif']
from clmm.dataops import _compute_lensing_angles_flatsky as clmm_flat
from clmm.dataops import _compute_lensing_angles_astropy as clmm_full
ra1, dec1 = 0, 0
t = np.linspace(0, np.pi*2, 100)[:-1]
ra2 = ra1+np.cos(t)
dec2 = dec1+np.sin(t)

flat_vals = clmm_flat(ra1, dec1, ra2, dec2)
full_vals = clmm_full(ra1, dec1, ra2, dec2)
f, axes = plt.subplots(2, sharex=True)

axes[0].plot(full_vals[1], flat_vals[0]-full_vals[0])
axes[0].set_ylabel('$\\theta$ [Flat-Curve]')

axes[1].plot(full_vals[1], flat_vals[1]-full_vals[1])
axes[1].set_xlabel('Azimuthal angle')
axes[1].set_ylabel('$\phi$ [Flat-Curve]')

for ax in axes:
    ax.set_xticks(np.arange(-1, 1.01, 0.5)*np.pi)
    ax.xaxis.grid(True, which='major', lw=.5)
    ax.yaxis.grid(True, which='major', lw=.5)
    ax.xaxis.grid(True, which='minor', lw=.1)
    ax.yaxis.grid(True, which='minor', lw=.1)
    
axes[1].set_xticklabels(['$-\pi$', '$-\\frac{\pi}{2}$', '$0$', '$\\frac{\pi}{2}$', '$\pi$'])<jupyter_output><empty_output><jupyter_text>Generating mock data `clmm` has a support code to generate a mock catalog given a input cosmology and cluster parameters. We will use this to generate a data sample to be used in this example:<jupyter_code>from clmm import Cosmology
import clmm.support.mock_data as mock
np.random.seed(14) # For reproducibility

# Set cosmology of mock data
cosmo = Cosmology(H0=70.0, Omega_dm0=0.27-0.045, Omega_b0=0.045, Omega_k0=0.0)

# Cluster info
cluster_m = 1.e15 # Cluster mass - ($M200_m$) [Msun]
concentration = 4  # Cluster concentration
cluster_z = 0.3 # Cluster redshift
cluster_ra = 20. # Cluster Ra in deg
cluster_dec = 85. # Cluster Dec in deg

# Catalog info
field_size = 10 # i.e. 10 x 10 Mpc field at the cluster redshift, cluster in the center

# Make mock galaxies
mock_galaxies = mock.generate_galaxy_catalog(
    cluster_m=cluster_m, cluster_z=cluster_z, cluster_c=concentration, cluster_ra=cluster_ra, cluster_dec=cluster_dec,# Cluster data
    cosmo=cosmo, # Cosmology object
    zsrc='desc_srd', # Galaxy redshift distribution, 
    zsrc_min=0.4, # Minimum redshift of the galaxies
    shapenoise=0.05, # Gaussian shape noise to the galaxy shapes
    photoz_sigma_unscaled=0.05, # Photo-z errors to source redshifts
    field_size=field_size,
    ngal_density=20 # number of gal/arcmin2 for z in [0, infty]
)['ra', 'dec', 'e1', 'e2', 'z', 'ztrue', 'pzbins', 'pzpdf', 'id']
print(f'Catalog table with the columns: {", ".join(mock_galaxies.colnames)}')

ngals_init = len(mock_galaxies)
print(f'Initial number of galaxies: {ngals_init:,}')

# Keeping only galaxies with "measured" redshift greater than cluster redshift
mock_galaxies = mock_galaxies[(mock_galaxies['z']>cluster_z)]
ngals_good = len(mock_galaxies)

if ngals_good < ngals_init:
    print(f'Number of excluded galaxies (with photoz < cluster_z): {ngals_init-ngals_good:,}')
    # reset galaxy id for later use
    mock_galaxies['id'] = np.arange(ngals_good)

# Check final density
from clmm.utils import convert_units
field_size_arcmin = convert_units(field_size, 'Mpc', 'arcmin', redshift=cluster_z, cosmo=cosmo)
print(f'Background galaxy density = {ngals_good/field_size_arcmin**2:.2f} gal/arcmin2\n')<jupyter_output><empty_output><jupyter_text>Shear values<jupyter_code>import clmm
import copy
# Create a GalaxyCluster
cluster_flat = clmm.GalaxyCluster("Name of cluster", cluster_ra, cluster_dec,
                                   cluster_z, mock_galaxies)
cluster = clmm.GalaxyCluster("Name of cluster", cluster_ra, cluster_dec,
                                   cluster_z, copy.deepcopy(mock_galaxies))

# Convert elipticities into shears for the members
cluster_flat.compute_tangential_and_cross_components(geometry="flat")
cluster.compute_tangential_and_cross_components(geometry="curve")
pass
def colorp(x, y, c, **kwargs):
    isort = np.argsort(c)
    return plt.scatter(x[isort], y[isort], c=c[isort], **kwargs)
sc = colorp(mock_galaxies['ra'], mock_galaxies['dec'], 
            c=abs(cluster_flat.galcat['et']-cluster.galcat['et']), s=1)
plt.colorbar(sc, label='$\Delta g_t$')
plt.show()
sc = colorp(mock_galaxies['ra'], mock_galaxies['dec'], 
            c=abs(cluster_flat.galcat['ex']-cluster.galcat['ex']), s=1)
plt.colorbar(sc, label='$\Delta g_x$')
plt.show()<jupyter_output><empty_output><jupyter_text>Profile values<jupyter_code>import clmm.dataops as da
# Measure profile and add profile table to the cluster
cluster_flat.make_radial_profile(bins=da.make_bins(0.1, field_size/2., 25, method='evenlog10width'),
                            bin_units="Mpc",
                            cosmo=cosmo,
                            include_empty_bins=False,
                            gal_ids_in_bins=True,
                           )
cluster.make_radial_profile(bins=da.make_bins(0.1, field_size/2., 25, method='evenlog10width'),
                            bin_units="Mpc",
                            cosmo=cosmo,
                            include_empty_bins=False,
                            gal_ids_in_bins=True,
                           )
pass
#from paper_formating import prep_plot
#prep_plot(figsize=(9, 9))
f, axes = plt.subplots(2, sharex=True)

errorbar_kwargs = dict(linestyle='', marker='o',
    markersize=1, elinewidth=.5, capthick=.5)

axes[0].errorbar(cluster.profile['radius'], cluster.profile['gt'],
             cluster.profile['gt_err'], c='k', label='curve', **errorbar_kwargs)
axes[0].errorbar(cluster_flat.profile['radius'], cluster_flat.profile['gt'],
             cluster_flat.profile['gt_err'], c='r', label='flat', **errorbar_kwargs)
axes[1].errorbar(cluster.profile['radius'],
                (cluster_flat.profile['gt']-cluster.profile['gt']),
                 cluster.profile['gt_err'], **errorbar_kwargs)
axes[1].set_xlabel('r [Mpc]', fontsize = 10)
axes[0].set_ylabel(r'$g_t$', fontsize = 10)
axes[1].set_ylabel(r'$\Delta g_t$', fontsize = 10)
axes[1].set_xscale('log')
axes[0].set_yscale('log')
axes[0].legend()

for ax in axes:
    ax.xaxis.grid(True, which='major', lw=.5)
    ax.yaxis.grid(True, which='major', lw=.5)
    ax.xaxis.grid(True, which='minor', lw=.1)
    ax.yaxis.grid(True, which='minor', lw=.1)<jupyter_output><empty_output><jupyter_text>Theoretical predictionsWe consider 3 models:1. One model where all sources are considered at the same redshift2. One model using the overall source redshift distribution to predict the reduced tangential shear3. A more accurate model, relying on the fact that we have access to the individual redshifts of the sources, where the average reduced tangential shear is averaged independently in each bin, accounting for the acutal population of sources in each bin.All models rely on `clmm.predict_reduced_tangential_shear` to make a prediction that accounts for the redshift distribution of the galaxies in each radial bin:  Model considering all sources located at the average redshift\begin{equation}     g_{t,i}^{\rm{avg(z)}} = g_t(R_i, \langle z \rangle)\;, \label{eq:wrong_gt_model} \end{equation}<jupyter_code>def predict_reduced_tangential_shear_mean_z(profile, logm):
    return clmm.compute_reduced_tangential_shear(
            r_proj=profile['radius'], # Radial component of the profile
            mdelta=10**logm, # Mass of the cluster [M_sun]
            cdelta=4, # Concentration of the cluster
            z_cluster=cluster_z, # Redshift of the cluster
            z_source=np.mean(cluster.galcat['z']), # Mean value of source galaxies redshift
            cosmo=cosmo,
            delta_mdef=200,
            halo_profile_model='nfw'
        )<jupyter_output><empty_output><jupyter_text>Model relying on the overall redshift distribution of the sources N(z), not using individual redshift information (eq. (6) from Applegate et al. 2014, MNRAS, 439, 48) \begin{equation}     g_{t,i}^{N(z)} = \frac{\langle\beta_s\rangle \gamma_t(R_i, z\rightarrow\infty)}{1-\frac{\langle\beta_s^2\rangle}{\langle\beta_s\rangle}\kappa(R_i, z\rightarrow\infty)}     \label{eq:approx_model} \end{equation}<jupyter_code>z_inf = 1000
dl_inf = cosmo.eval_da_z1z2(cluster_z, z_inf)
d_inf = cosmo.eval_da(z_inf)

def betas(z):
    dls = cosmo.eval_da_z1z2(cluster_z, z)
    ds = cosmo.eval_da(z)
    return dls * d_inf / (ds * dl_inf)

def predict_reduced_tangential_shear_approx(profile, logm):

    bs_mean = np.mean(betas(cluster.galcat['z'])) 
    bs2_mean = np.mean(betas(cluster.galcat['z'])**2)

    gamma_t_inf = clmm.compute_tangential_shear(
            r_proj=profile['radius'], # Radial component of the profile
            mdelta=10**logm, # Mass of the cluster [M_sun]
            cdelta=4, # Concentration of the cluster
            z_cluster=cluster_z, # Redshift of the cluster
            z_source=z_inf, # Redshift value at infinity
            cosmo=cosmo,
            delta_mdef=200,
            halo_profile_model='nfw')
    convergence_inf = clmm.compute_convergence(
            r_proj=profile['radius'], # Radial component of the profile
            mdelta=10**logm, # Mass of the cluster [M_sun]
            cdelta=4, # Concentration of the cluster
            z_cluster=cluster_z, # Redshift of the cluster
            z_source=z_inf, # Redshift value at infinity
            cosmo=cosmo,
            delta_mdef=200,
            halo_profile_model='nfw')
        
    return bs_mean*gamma_t_inf/(1-(bs2_mean/bs_mean)*convergence_inf)<jupyter_output><empty_output><jupyter_text>Model using individual redshift and radial information, to compute the averaged shear in each radial bin, based on the galaxies actually present in that bin.\begin{equation}    g_{t,i}^{z, R} = \frac{1}{N_i}\sum_{{\rm gal\,}j\in {\rm bin\,}i} g_t(R_j, z_j)    \label{eq:exact_model} \end{equation}<jupyter_code>cluster.galcat['theta_mpc'] = convert_units(cluster.galcat['theta'], 'radians', 'mpc',cluster.z, cosmo)

def predict_reduced_tangential_shear_exact(profile, logm):
    return np.array([np.mean(
        clmm.compute_reduced_tangential_shear(
            # Radial component of each source galaxy inside the radial bin
            r_proj=cluster.galcat[radial_bin['gal_id']]['theta_mpc'],
            mdelta=10**logm, # Mass of the cluster [M_sun]
            cdelta=4, # Concentration of the cluster
            z_cluster=cluster_z, # Redshift of the cluster
            # Redshift value of each source galaxy inside the radial bin
            z_source=cluster.galcat[radial_bin['gal_id']]['z'],
            cosmo=cosmo,
            delta_mdef=200,
            halo_profile_model='nfw'
        )) for radial_bin in profile])<jupyter_output><empty_output><jupyter_text>Mass fitting We estimate the best-fit mass using `scipy.optimize.curve_fit`. The choice of fitting $\log M$ instead of $M$ lowers the range of pre-defined fitting bounds from several order of magnitude for the mass to unity. From the associated error $\sigma_{\log M}$ we calculate the error to mass as $\sigma_M = M_{fit}\ln(10)\sigma_{\log M}$.  First, identify bins with sufficient galaxy statistics to be kept for the fitFor small samples, error bars should not be computed using the simple error on the mean approach available so far in CLMM)<jupyter_code>data_for_fit = cluster.profile[cluster.profile['n_src'] > 5]
data_for_fit_flat = cluster_flat.profile[cluster_flat.profile['n_src'] > 5]<jupyter_output><empty_output><jupyter_text>Perform the fits<jupyter_code>from clmm.support.sampler import fitters
def fit_mass(predict_function, data_for_fit=data_for_fit):
    popt, pcov = fitters['curve_fit'](predict_function,
        data_for_fit, 
        data_for_fit['gt'], 
        data_for_fit['gt_err'], bounds=[10.,17.])
    logm, logm_err = popt[0], np.sqrt(pcov[0][0])
    return {'logm':logm, 'logm_err':logm_err,
            'm': 10**logm, 'm_err': (10**logm)*logm_err*np.log(10)}
fit_mean_z = fit_mass(predict_reduced_tangential_shear_mean_z)
fit_approx = fit_mass(predict_reduced_tangential_shear_approx)
fit_exact = fit_mass(predict_reduced_tangential_shear_exact)
fit_mean_z_flat = fit_mass(predict_reduced_tangential_shear_mean_z, data_for_fit_flat)
fit_approx_flat = fit_mass(predict_reduced_tangential_shear_approx, data_for_fit_flat)
fit_exact_flat = fit_mass(predict_reduced_tangential_shear_exact, data_for_fit_flat)
print(f'Input mass = {cluster_m:.2e} Msun\n')
print('Curve sky')
print(f'Best fit mass for average redshift               = {fit_mean_z["m"]:.3e} +/- {fit_mean_z["m_err"]:.3e} Msun')
print(f'Best fit mass for N(z) model                     = {fit_approx["m"]:.3e} +/- {fit_approx["m_err"]:.3e} Msun')
print(f'Best fit mass for individual redshift and radius = {fit_exact["m"]:.3e} +/- {fit_exact["m_err"]:.3e} Msun')
print()
print('Flat sky')
print(f'Best fit mass for average redshift               = {fit_mean_z_flat["m"]:.3e} +/- {fit_mean_z_flat["m_err"]:.3e} Msun')
print(f'Best fit mass for N(z) model                     = {fit_approx_flat["m"]:.3e} +/- {fit_approx_flat["m_err"]:.3e} Msun')
print(f'Best fit mass for individual redshift and radius = {fit_exact_flat["m"]:.3e} +/- {fit_exact_flat["m_err"]:.3e} Msun')
print()
print('Difference (Curve-Flat)')
print(f'Best fit mass for average redshift               = {fit_mean_z_flat["m"]-fit_mean_z["m"]:.3e} Msun')
print(f'Best fit mass for N(z) model                     = {fit_approx_flat["m"]-fit_approx["m"]:.3e} Msun')
print(f'Best fit mass for individual redshift and radius = {fit_exact_flat["m"]-fit_exact["m"]:.3e} Msun')<jupyter_output><empty_output><jupyter_text>Reduced $\chi$2<jupyter_code>def get_predicted_shear(predict_function, fit_values):
    gt_est = predict_function(data_for_fit, fit_values['logm'])
    gt_est_err = [predict_function(data_for_fit, fit_values['logm']+i*fit_values['logm_err'])
                          for i in (-3, 3)]
    return gt_est, gt_est_err
def chi2(predict_function, fit_values):
    gt_mean_z, gt_err_mean_z =  get_predicted_shear(predict_function, fit_values)
    return np.sum((gt_mean_z-data_for_fit['gt'])**2/(data_for_fit['gt_err'])**2)/(len(data_for_fit)-1)

print('Curve sky')
print(f'Reduced chi2 (mean z model) = {chi2(predict_reduced_tangential_shear_mean_z, fit_mean_z)}')
print(f'Reduced chi2 (N(z) model) = {chi2(predict_reduced_tangential_shear_approx, fit_approx)}')
print(f'Reduced chi2 (individual (R,z) model) = {chi2(predict_reduced_tangential_shear_exact, fit_exact)}')
print()
print('Flat sky')
print(f'Reduced chi2 (mean z model) = {chi2(predict_reduced_tangential_shear_mean_z, fit_mean_z_flat)}')
print(f'Reduced chi2 (N(z) model) = {chi2(predict_reduced_tangential_shear_approx, fit_approx_flat)}')
print(f'Reduced chi2 (individual (R,z) model) = {chi2(predict_reduced_tangential_shear_exact, fit_exact_flat)}')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12634753593429157, 'fraction_numerical': 0.016940451745379878, 'mean_word_length': 4.342817963661296, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 44, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '19447', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5651, 'n_tokens_neox': 5330, 'n_words': 1219, 'path': 'examples/Flat_sky_limitations.ipynb', 'repo_name': 'nicolahunfeld/CLMM', 'size': '15584'}