%matplotlib inline
%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import scipy.stats as stats
from matplotlib import cm
from astropy.table import Table
import emcee

import clmm
from clmm import GCData, GalaxyCluster
from clmm import Cosmology, utils as clmm_utils
from clmm.support.sampler import fitters
cosmo = clmm.Cosmology(H0=70.0, Omega_dm0=0.3 - 0.045, Omega_b0=0.045, Omega_k0=0.0)

from utils import (
    anacal_get_tang_cross, calibrate_shapes,
    build_flat_wcs_grid, compute_mass_map, schirmer_filter,
)

z_cl = 0.22

# Assume sources follow the DESC SRD redshift distribution.
z_distrib_func = clmm_utils.redshift_distributions.desc_srd
beta_kwargs = {
    "z_cl": z_cl,
    "z_inf": 10.0,
    "cosmo": cosmo,
    "z_distrib_func": z_distrib_func,
}
beta_s_mean = clmm_utils.compute_beta_s_mean_from_distribution(**beta_kwargs)
beta_s_square_mean = clmm_utils.compute_beta_s_square_mean_from_distribution(**beta_kwargs)
rproj = np.logspace(np.log10(0.3), np.log10(7.), 100)

def Da(z):
    return cosmo.eval_da(z)

import os
import tables_io
from astropy.io import fits

cluster_coords = (37.865017, 6.982205)  # BCG of Abell 360
ra_bcg, dec_bcg = cluster_coords

Dir = "/gpfs/mnt/gpfs02/astro/workarea/xli6/work/DP1"
anacal_table = Table.read(os.path.join(Dir, "catalogs/anacal_catalog_a360.fits"))
anacal_table = anacal_table[anacal_table["is_primary"]]

# Cached photo-z point estimates + per-source PDFs (row-aligned with the
# is_primary anacal table).
a360_redshifts = tables_io.read(
    os.path.join(Dir, "rail_data/data_a360_redshift.hdf5"),
    tType="astropyTable",
)
with fits.open(os.path.join(Dir, "rail_data/data_a360_redshift_pdfs.fits")) as hdul:
    a360_pdfs = hdul[0].data           # shape (N, 501); z-grid = arange(0, 5.01, 0.01)
anacal_table["zbest"] = a360_redshifts["zmode"]
anacal_table["index"] = np.arange(len(anacal_table))

# Derive AB-zero=31.4 mags from the gauss2 fluxes.
mag_zero = 31.4
for b in ("g", "r", "i", "z"):
    anacal_table[f"{b}_mag"] = mag_zero - 2.5 * np.log10(
        anacal_table[f"{b}_flux_gauss2"]
    )

anacal_table0 = anacal_table

/tmp/ipykernel_147858/699198058.py:26: RuntimeWarning: invalid value encountered in log10
  anacal_table[f"{b}_mag"] = mag_zero - 2.5 * np.log10(

plt.scatter(anacal_table["ra"], anacal_table["dec"], s=0.1)

<matplotlib.collections.PathCollection at 0x14eb4cce54f0>

print(len(anacal_table0))

113945

plt.hist(anacal_table0["z_mag"], bins=50)

(array([2.0000e+00, 2.8000e+01, 5.4000e+01, 1.1700e+02, 1.8800e+02,
        2.1200e+02, 2.5900e+02, 3.5400e+02, 4.4800e+02, 6.2700e+02,
        7.8700e+02, 1.1260e+03, 1.5450e+03, 2.0920e+03, 2.9340e+03,
        4.1020e+03, 5.5210e+03, 7.6640e+03, 1.0901e+04, 1.4654e+04,
        1.7363e+04, 1.5074e+04, 1.0107e+04, 5.3420e+03, 3.1340e+03,
        1.7630e+03, 1.0820e+03, 6.7700e+02, 4.4800e+02, 2.9300e+02,
        1.9100e+02, 9.7000e+01, 9.4000e+01, 6.3000e+01, 4.2000e+01,
        3.9000e+01, 1.9000e+01, 1.1000e+01, 7.0000e+00, 8.0000e+00,
        2.0000e+00, 2.0000e+00, 0.0000e+00, 2.0000e+00, 0.0000e+00,
        0.0000e+00, 0.0000e+00, 1.0000e+00, 0.0000e+00, 1.0000e+00]),
 array([15.12294085, 15.53348562, 15.94403039, 16.35457516, 16.76511993,
        17.1756647 , 17.58620946, 17.99675423, 18.407299  , 18.81784377,
        19.22838854, 19.63893331, 20.04947808, 20.46002285, 20.87056762,
        21.28111239, 21.69165716, 22.10220193, 22.5127467 , 22.92329147,
        23.33383623, 23.744381  , 24.15492577, 24.56547054, 24.97601531,
        25.38656008, 25.79710485, 26.20764962, 26.61819439, 27.02873916,
        27.43928393, 27.8498287 , 28.26037347, 28.67091823, 29.081463  ,
        29.49200777, 29.90255254, 30.31309731, 30.72364208, 31.13418685,
        31.54473162, 31.95527639, 32.36582116, 32.77636593, 33.1869107 ,
        33.59745547, 34.00800024, 34.418545  , 34.82908977, 35.23963454,
        35.65017931]),
 <BarContainer object of 50 artists>)

msk = np.logical_and.reduce((
    (anacal_table0["zbest"] > 0.40),
    (anacal_table0["zbest"] < 2.0),
    (anacal_table0["wsel"] > 1e-5),
    (anacal_table0["g_mag"] < 25.5),
    (anacal_table0["r_mag"] < 25.0),
    (anacal_table0["i_mag"] < 23.5),
    (anacal_table0["z_mag"] < 24.5),
    (((anacal_table0["i_fpfs1_m20"] + anacal_table0["i_fpfs1_m00"]) / anacal_table0["i_fpfs1_m00"]) > 0.1),
))
print(np.sum(msk))
anacal_table = anacal_table0[msk]

26591

e1, e2, res = calibrate_shapes(anacal_table)
anacal_table["response"] = res

source_phi = np.arctan2(
    anacal_table['dec'] - cluster_coords[1],
    (cluster_coords[0] - anacal_table['ra']) * np.cos(np.deg2rad(cluster_coords[1])),
)
ang_dist = np.sqrt(
    ((anacal_table['ra'] - cluster_coords[0]) * np.cos(np.deg2rad(cluster_coords[1]))) ** 2
    + (anacal_table['dec'] - cluster_coords[1]) ** 2
)
sky_distance = Da(z_cl) * ang_dist * (np.pi / 180)

trial_shear = e1 + 1j * e2
cl_shear = trial_shear * -np.exp(-2j * source_phi)

bins_mpc = np.array([0.3, 0.7, 1.06495522, 1.62018517, 2.46489237, 3.75, 5, 6])
bin_mids = 0.5 * (bins_mpc[1:] + bins_mpc[:-1])
shear_cl = anacal_get_tang_cross(
    cl_shear, sky_distance, bins_mpc, res, ci_level=.68,
)

moo = clmm.Modeling(massdef="critical", delta_mdef=500, halo_profile_model="nfw")
moo.set_cosmo(cosmo)
moo.set_concentration(4)
moo.set_mass(6e14)
gt_z = moo.eval_reduced_tangential_shear(
    rproj, z_cl, [beta_s_mean, beta_s_square_mean],
    z_src_info="beta", approx="order2",
)

plt.close()
cmap = cm.coolwarm

plt.plot(rproj, gt_z, label=r'NFW, $M_{500c}$=$6 \times 10^{14} M_{\odot}$, c=4, n(z)=SRD', ls=':')
plt.plot(bin_mids, shear_cl[0], '.', label='tangential', color=cmap(.05))
plt.vlines(bin_mids, shear_cl[2][:,0], shear_cl[2][:,1], color=cmap(0.05))

plt.plot(1.05*bin_mids, shear_cl[1], '.', label='cross', color=cmap(.95))
plt.vlines(1.05*bin_mids, shear_cl[3][:,0], shear_cl[3][:,1], color=cmap(0.95))
chi2_t = np.sum(shear_cl[0] ** 2.0 / ((shear_cl[2][:,0] - shear_cl[2][:,1]) / 2.0) ** 2.0)
chi2_x = np.sum(shear_cl[1] ** 2.0 / ((shear_cl[3][:,0] - shear_cl[3][:,1]) / 2.0) ** 2.0)

plt.semilogx()
plt.axhline(0, ls='--', color='k', alpha=.5)

plt.ylim([-0.05,0.10])
plt.xlim([0.45,7.5])
plt.ylabel("Reduced shear")
plt.xlabel("R (Mpc)")
plt.legend(frameon=False, loc='upper right', fontsize=15)

<matplotlib.legend.Legend at 0x14eb2a74f710>

print(np.sqrt(chi2_t))
print(np.sqrt(chi2_x))

4.581021216262394
1.76505803184541

test_zs = np.arange(0, 5.01, 0.01)

nz_ndxs = anacal_table["index"]
response_stacked_nz = np.dot(anacal_table["response"], a360_pdfs[nz_ndxs])
nz_normalization = integrate.trapezoid(response_stacked_nz, test_zs)
normalized_nz = response_stacked_nz / nz_normalization

full_nz = np.sum(a360_pdfs, axis=0)
full_norm_nz = full_nz / integrate.trapezoid(full_nz, test_zs)

cmap = cm.coolwarm
plt.plot(test_zs, normalized_nz, '-', label='Source Sample', color=cmap(0.95))
plt.plot(test_zs, full_norm_nz, '-', label='Full Sample', color=cmap(0.05))
plt.xlabel("Redshift")
plt.ylabel("Normalized $N(z)$")
plt.title("Abell 360 Redshift Distribution")
plt.axvline(z_cl, ls='--', color='k', alpha=0.3)
plt.legend()

<matplotlib.legend.Legend at 0x14eb2a3c4920>

def Daz1z2(z1, z2):
    return cosmo.eval_da_z1z2(z1, z2)

def beta_s_of_zs(zl, zs):
    if zs <= zl:
        return 0.0
    return Daz1z2(zl, zs) / cosmo.eval_da(zs) * cosmo.eval_da(1e3) / Daz1z2(zl, 1e3)


def get_beta_betasqr_pdf(nz, beta_s, beta_sqr, zs):
    norm = integrate.trapezoid(nz, zs)
    return (
        integrate.trapezoid(beta_s * nz, zs) / norm,
        integrate.trapezoid(beta_sqr * nz, zs) / norm,
    )


beta_s_vals = np.array([beta_s_of_zs(z_cl, zz) for zz in test_zs])
beta_sqr_vals = beta_s_vals ** 2
cluster_beta_s, cluster_beta_s_sqr = get_beta_betasqr_pdf(
    normalized_nz, beta_s_vals, beta_sqr_vals, test_zs,
)
print(f"cluster_beta_s     = {cluster_beta_s:.4f}")
print(f"cluster_beta_s_sqr = {cluster_beta_s_sqr:.4f}")

cluster_beta_s     = 0.6577
cluster_beta_s_sqr = 0.4609

def fit_mass(predict_function, fit_data, **kwargs):
    popt, pcov = fitters["curve_fit"](
        predict_function,
        fit_data["radius"],
        fit_data["gt"],
        fit_data["gt_err"],
        **kwargs,
    )
    return {"vals": popt, "err": np.sqrt(np.diag(pcov)), "cov": pcov}


def fixed_concentration_fit(profile, beta, betasqr, cdelta, halo_profile="nfw"):
    fn = lambda radius, logm: clmm.compute_reduced_tangential_shear(
        r_proj=radius,
        mdelta=10 ** logm,
        cdelta=cdelta,
        z_cluster=z_cl,
        z_src=(beta, betasqr),
        z_src_info="beta",
        approx="order2",
        cosmo=cosmo,
        delta_mdef=500,
        massdef="critical",
        halo_profile_model=halo_profile,
    )
    mass_fit = fit_mass(fn, profile.profile, bounds=[10.0, 17.0])
    return fn, mass_fit


def mcmc_fit_1d(profile, fn, init, nwalkers=32, nsteps=100):
    fit_data = profile.profile

    def log_prob(pars):
        if pars[0] < 10 or pars[0] > 16:
            return -np.inf
        model = fn(fit_data["radius"], *pars)
        diff = model - fit_data["gt"]
        sig = fit_data["gt_err"]
        return float(-0.5 * np.sum(diff ** 2 / sig ** 2))

    sampler = emcee.EnsembleSampler(nwalkers, 1, log_prob)
    p0 = np.random.normal(init["vals"], init["err"] * 0.3, (nwalkers, 1))
    p0[p0 < 0] *= -1
    sampler.run_mcmc(p0, nsteps, progress=True)
    return {
        "mcmc_chain": sampler.get_chain().flatten(),
        "mcmc_loglike": sampler.get_log_prob().flatten(),
    }

global_R_fit = float(np.mean(res))
print(f"global R (mass fit) = {global_R_fit:.4f}")

galcat = GCData()
galcat['ra'] = anacal_table['ra']
galcat['dec'] = anacal_table['dec']
galcat['e1'] = e1 / global_R_fit   # AnaCal response correction
galcat['e2'] = e2 / global_R_fit
galcat['id'] = np.arange(len(anacal_table))
galcat['z'] = anacal_table["zbest"]

gc_object = GalaxyCluster("A360", ra_bcg, dec_bcg, z_cl, galcat)
gc_object.compute_tangential_and_cross_components(add=True);
gc_object.make_radial_profile(
    bins=bins_mpc, bin_units='Mpc', add=True,
    cosmo=cosmo, overwrite=True,
    use_weights=False, gal_ids_in_bins=False,
);
profile = gc_object.profile

global R (mass fit) = 0.3460

/gpfs/mnt/gpfs02/astro/workarea/xli6/superonionSDCC/miniconda3/envs/analysis/lib/python3.12/site-packages/clmm/gcdata.py:30: UserWarning: coordinate_system not set, defaulting to 'euclidean'.
  warnings.warn("coordinate_system not set, defaulting to 'euclidean'.")

fn, init = fixed_concentration_fit(
    gc_object, cluster_beta_s, cluster_beta_s_sqr, 4, "nfw",
)
print(f"curve_fit: logM = {init['vals'][0]:.3f} ± {init['err'][0]:.3f}")

mcmc_results = mcmc_fit_1d(gc_object, fn, init, nwalkers=32, nsteps=100)
chain = mcmc_results['mcmc_chain']
ll = mcmc_results['mcmc_loglike']
fin = np.isfinite(ll)
chain = chain[fin]
ll = ll[fin]

curve_fit: logM = 14.630 ± 0.145

100%|██████████| 100/100 [00:08<00:00, 12.19it/s]

  2%|▏         | 2/100 [00:00<00:08, 12.13it/s]

  4%|▍         | 4/100 [00:00<00:09, 10.66it/s]

  6%|▌         | 6/100 [00:00<00:08, 11.52it/s]

  8%|▊         | 8/100 [00:00<00:07, 11.70it/s]

 10%|█         | 10/100 [00:00<00:07, 11.99it/s]

 12%|█▏        | 12/100 [00:01<00:07, 12.14it/s]

 14%|█▍        | 14/100 [00:01<00:06, 12.39it/s]

 16%|█▌        | 16/100 [00:01<00:06, 12.54it/s]

fig = plt.figure(figsize=(7, 3))
cmap = cm.coolwarm
mass_bins = np.arange(13.1, 15.5, 0.05)
plt.hist(chain, bins=mass_bins, weights=np.exp(ll), density=True,
         linewidth=2, histtype='step', color=cmap(0.05))
plt.hist(chain, bins=mass_bins, weights=np.exp(ll), density=True,
         linewidth=2, histtype='stepfilled', color=cmap(0.05), alpha=0.2)
plt.xlim(13, 16)
plt.xlabel(r"$\log{M_{500c}/M_\odot}$")
plt.ylabel("Normalized counts")
plt.title("Abell 360 — NFW mass posterior (fixed c=4)")

Text(0.5, 1.0, 'Abell 360 — NFW mass posterior (fixed c=4)')

N, rr = 121, 0.5
x, y, x_grid, y_grid = build_flat_wcs_grid(
    ra_bcg, dec_bcg, anacal_table['ra'], anacal_table['dec'],
    n=N, half_extent_deg=rr,
)
global_R = float(np.mean(res))
print(f"global R = {global_R:.4f}")
g1 = e1 / global_R
g2 = e2 / global_R
weights = np.ones(len(x))

e_ap, b_ap, v_ap = compute_mass_map(
    x_grid, y_grid, x, y, g1, g2, weights,
    schirmer_filter,
    filter_kwargs={'aperture_size': 0.4, 'a': 15, 'b': 150},
)

global R = 0.3460

test_xs = np.linspace(-5, 5, num=1001)
normal_dist = stats.norm(0, 1).pdf(test_xs)
plt.hist((b_ap / np.sqrt(v_ap)).flatten(), bins=51, label='B Mode',
         histtype='step', density=True)
plt.hist((e_ap / np.sqrt(v_ap)).flatten(), bins=51, label='E Mode',
         histtype='step', density=True)
plt.plot(test_xs, normal_dist, '--')
plt.yscale('log')
plt.ylim(1e-4, 5e-1)
plt.legend()

<matplotlib.legend.Legend at 0x14eb2a2f9790>

fig, axs = plt.subplots(ncols=2, figsize=(10, 6))
extent = (rr, -rr, -rr, rr)
sn_v = np.sqrt(np.mean(v_ap))

axs[0].set_title('E-Mode SN')
MapE = axs[0].imshow(e_ap / sn_v, vmin=-2, vmax=5, origin='lower', extent=extent)
y_p, x_p = np.where(e_ap / sn_v == np.max(e_ap / sn_v))
print((e_ap / sn_v)[y_p, x_p])
x_vals = np.linspace(extent[0], extent[1], N)
y_vals = np.linspace(extent[2], extent[3], N)
axs[0].scatter(x_vals[x_p], y_vals[y_p], s=10.0, color='red')

axs[1].set_title('B-Mode SN')
axs[1].imshow(b_ap / sn_v, vmin=-2, vmax=5, origin='lower', extent=extent)

fig.colorbar(MapE, ax=axs, fraction=0.025)
fig.supxlabel('RA',  y=0.14)
fig.supylabel('DEC', x=0.05)

[5.35517838]

Text(0.05, 0.5, 'DEC')

3.2 — Abell 360 cluster¶

Purpose¶

Steps¶

Inputs¶

Outputs¶

Source-sample n(z) and β_s¶

CLMM mass fit¶

Schirmer aperture-mass S/N map¶