NumCosmoGalaxyShapeFactorQuad

Exact quadrature evaluation of the intrinsic-ellipticity marginal.

Computes $$P(\epsilon_\mathrm{obs} \mid g) = \int_{|\chi_I|<1} \mathrm{d}^2\chi_I\, P_\mathrm{pop}(\chi_I)\, N_2\big(\epsilon_\mathrm{obs} - f_g(\chi_I); \sigma_\mathrm{noise}^2\big)$$ exactly, with no linearization of the shear map and no truncation of the intrinsic population to an untruncated Gaussian (contrast NcGalaxyShapeFactorVarAdd, which approximates both).

The substitution variable is the LENSED (noiseless, pre-noise) ellipticity $\chi_L=f_g(\chi_I)$, not the intrinsic $\chi_I$: change variables in the integral above to $\chi_L$, $$P(\epsilon_\mathrm{obs} \mid g) = \int_{|\chi_L|<1} \mathrm{d}^2\chi_L\, P_\mathrm{pop}\big(f_g^{-1}(\chi_L)\big)\, \left|\det J_{f_g^{-1}}(\chi_L)\right|\, N_2\big(\epsilon_\mathrm{obs} - \chi_L; \sigma_\mathrm{noise}^2\big),$$ where $\det J_{f_g^{-1}}$ is supplied by the existing nc_galaxy_shape_factor_lndet_jac() machinery (already used by VarAdd), so no new map algebra is needed. $\chi_L$ is turned into a plane integral by the usual bijection, $$\chi_L(u,v) = (u,v)/\sqrt{1+u^2+v^2}, \qquad \left|\det\frac{\partial\chi_L}{\partial(u,v)}\right| = \frac{1}{(1+u^2+v^2)^2}.$$ The point of using $\chi_L$ rather than $\chi_I$ is that the noise kernel $N_2(\epsilon_\mathrm{obs}-\chi_L;\sigma_\mathrm{noise}^2)$ becomes an ordinary Gaussian of EXACTLY known width $\sigma_\mathrm{noise}$ centered at EXACTLY $\epsilon_\mathrm{obs}$ — no shear map, no approximation. Because the population is evaluated at $\chi_I=f_g^{-1}(\chi_L)$ instead of directly at the substitution point, nc_galaxy_shape_pop_eval_p() is used (not nc_galaxy_shape_pop_eval_p_rho2(), whose $\rho^2$ contract is specifically $\rho^2=u^2+v^2$ for the point being substituted for, which is $\chi_L$ here, not $\chi_I$); $x_I=\lvert\chi_I\rvert^2$ is computed directly from the complex division (a sum of two squares, not itself subtraction-sensitive), only losing the extra conditioning nc_galaxy_shape_pop_eval_p_rho2() offers for the $(1-x)^a$ factor right at $x\to1^-$ (relevant mainly for NcGalaxyShapePopBeta at the disc boundary).

The plane integral is evaluated by the Divonne cubature algorithm from the Cuba library (ncm_integrate_2dim_divonne()), over a FIXED box of half-width bound centered on the exact $(u,v)$ preimage of $\epsilon_\mathrm{obs}$ (no shear map needed: $\chi_L$’s noise kernel peaks exactly there), seeded with two explicit peak hints: that same point (the noise peak, exact) and the preimage of $f_g(0)$ (where $\chi_L$ sits when $\chi_I=0$, i.e. the population’s peak location under the common assumption that $P_\mathrm{pop}$ is radially symmetric about $\chi_I=0$, e.g. NcGalaxyShapePopGauss / NcGalaxyShapePopGaussLocal — for a population peaked elsewhere, e.g. NcGalaxyShapePopBeta with $\mu$ away from 0, this hint is only approximate, but Divonne’s own stratified search still reliably finds the true peak from there in every case tested, see below). $f_g(0)$ is convention-dependent — exactly $g$ in the TRACE_DET (ellipticity) convention, but $2g/(1+\lvert g\rvert^2)$ in the TRACE (distortion) convention — so it is computed via the actual forward map rather than hardcoded as $g$, since hardcoding it would miss the peak for narrow, off-center TRACE-convention populations.

Divonne’s explicit peak hints matter: a fixed-degree base rule with no explicit hints (e.g. Cuba’s Cuhre) has no mechanism to notice an isolated feature it never samples near, and can return a confidently wrong result (off by orders of magnitude, at tight reltol) whenever the population is narrow relative to the integration box and off-center. Because Divonne is told where to look, bound can be a single generous fixed constant: shrinking it does not help an already-hint-found narrow feature, but does progressively truncate a broad population’s disc-spanning support (a few units already covers the disc for any population). Verified against an independent scipy reference across populations from $\sigma=0.30$ down to $\sigma=0.001$, five different $g$/$\epsilon_\mathrm{obs}$/ $\sigma_\mathrm{noise}$ geometries (including near the disc boundary), and NcGalaxyShapePopBeta concentrations up to $\nu=10^5$: every case matches to $\sim!10^{-6}$ relative accuracy or better. See docs/theory/wl_shape_factor_history.md for earlier implementations of this class that were tried and rejected.

The integrand clamps any non-finite evaluation to zero before returning: Cuba can segfault outright on a NaN/Inf sample (reproduced directly), and some populations have a genuine, mathematically correct divergence at a disc point (e.g. NcGalaxyShapePopBeta with $\alpha<1$ diverges at $x=0$, which can coincide with a peak hint).

This scheme is exact but substantially more expensive per evaluation than VarAdd (a full 2D cubature vs. one closed-form expression): Divonne with two hints typically costs $\sim!100\,\mathrm{ms}$ per evaluation, rising to several seconds for extremely concentrated populations. It is meant as an accuracy reference / fallback for regimes where the variance-add approximation is not trusted, not as a routine replacement in large-catalog likelihood evaluations. Cuba’s own fork()-based internal parallelism is disabled globally by ncm_cfg_init() (cubacores(0, 0)), so this is safe to call concurrently from multiple OpenMP threads.