From d6531b1bb70b6a3bc612f883f5f5ac019f21f838 Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Wed, 10 Jun 2026 03:56:14 +0200 Subject: [PATCH] fix(ngmix): emit true half-light radii in r50 columns; dedupe PSF size columns The v2.0 size columns were mislabeled: galaxy r50/r50_err were bit-for-bit copies of the area T = 2 sigma^2 (pars[4]), and the PSF r50 columns held bare sigma, missing the sqrt(2 ln 2) = 1.17741 factor. No column in the file was a true half-light radius, and the same name root meant an area on the galaxy side and a length on the PSF side. Now every r50 column is a genuine half-light radius (r50 = 1.17741 sigma = sqrt(ln 2 * T)) with propagated errors, matching the UNIONS-3500 WL I convention (r_h primary, T = 2 sigma^2 the derived DES area): - galaxy: r50 = sqrt(ln 2 * T), r50_err = r50 * T_err / (2 T); NaN when the fitted T is non-positive - PSF: r50psf = sqrt(2 ln 2) * sigma_psf, with the corrected error d sigma / d T = 1 / (4 sigma) (the previous T_err / (2 sigma) was a factor-2 over-estimate) - retire the duplicate columns T_psfo_ngmix (== Tpsf), T_err_psfo_ngmix, r50_psfo_ngmix, r50_err_psfo_ngmix (== r50psf and its error); add the missing Tpsf_err so areas carry their error on both sides - make_cat: read Tpsf instead of the retired T_psfo_ngmix (output column NGMIX_T_PSFo_ unchanged, so downstream consumers are unaffected) - fix the get_noise docstring: guess[4] is T, not r50 Tests pin the new semantics: both r50 columns on the same scale (their ratio is sqrt(T/Tpsf)), retired names absent, NaN for non-positive T. Co-Authored-By: Claude Fable 5 --- .../modules/make_cat_package/make_cat.py | 2 +- src/shapepipe/modules/ngmix_package/ngmix.py | 50 +++++---- src/shapepipe/tests/test_ngmix.py | 104 ++++++++++++++++++ 3 files changed, 135 insertions(+), 21 deletions(-) diff --git a/src/shapepipe/modules/make_cat_package/make_cat.py b/src/shapepipe/modules/make_cat_package/make_cat.py index 72024b3a4..890e24713 100644 --- a/src/shapepipe/modules/make_cat_package/make_cat.py +++ b/src/shapepipe/modules/make_cat_package/make_cat.py @@ -449,7 +449,7 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False): ) self._add2dict(f"NGMIX{m}_ELL_PSFo_{key}", g_psf, idx) - t_psfo = ncf_data["T_psfo_ngmix"][ind[0]] + t_psfo = ncf_data["Tpsf"][ind[0]] self._add2dict(f"NGMIX{m}_T_PSFo_{key}", t_psfo, idx) self._add2dict( diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 630e1ade8..e6f183875 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -19,6 +19,11 @@ from shapepipe.pipeline import file_io +# Gaussian half-light radius per unit sigma: r50 = sqrt(2 ln 2) * sigma +# = 1.17741 * sigma. With the ngmix area parameter T = 2 sigma^2 this +# gives r50 = SIGMA_TO_R50 * sqrt(T / 2) = sqrt(ln 2 * T). +SIGMA_TO_R50 = np.sqrt(2.0 * np.log(2.0)) + def get_prior(pixel_scale, rng, T_range=None, F_range=None): """Build ngmix joint prior for a 6-parameter galaxy model. @@ -331,12 +336,8 @@ def compile_results(self, results): 'ntry_fit', 'g1_psfo_ngmix', 'g2_psfo_ngmix', - 'T_psfo_ngmix', - 'T_err_psfo_ngmix', - 'r50_psfo_ngmix', 'g1_err_psfo_ngmix', 'g2_err_psfo_ngmix', - 'r50_err_psfo_ngmix', 'g1', 'g1_err', 'g2', @@ -344,9 +345,11 @@ def compile_results(self, results): 'T', 'T_err', 'Tpsf', + 'Tpsf_err', 'r50', 'r50_err', 'r50psf', + 'r50psf_err', 'g1_psf', 'g2_psf', 'flux', @@ -391,24 +394,28 @@ def compile_results(self, results): output_dict[name]["g2_err_psfo_ngmix"].append( results[idx]["g_err_PSFo"][1] ) - output_dict[name]["T_psfo_ngmix"].append(results[idx]["T_PSFo"]) - output_dict[name]["T_err_psfo_ngmix"].append( - results[idx]["T_err_PSFo"] - ) - output_dict[name]['r50_psfo_ngmix'].append( - results[idx]['r50_PSFo'] - ) - output_dict[name]['r50_err_psfo_ngmix'].append( - results[idx]['r50_err_PSFo'] - ) output_dict[name]["T"].append(results[idx][name]["T"]) output_dict[name]["T_err"].append(results[idx][name]["T_err"]) output_dict[name]["Tpsf"].append(results[idx]["T_PSFo"]) + output_dict[name]["Tpsf_err"].append(results[idx]["T_err_PSFo"]) output_dict[name]["g1_psf"].append(results[idx]["g_PSFo"][0]) output_dict[name]["g2_psf"].append(results[idx]["g_PSFo"][1]) - output_dict[name]['r50'].append(results[idx][name]['pars'][4]) - output_dict[name]['r50_err'].append(results[idx][name]['pars_err'][4]) + + # Galaxy half-light radius from the fitted area T = 2 sigma^2: + # r50 = sqrt(ln 2 * T), with d r50 / d T = r50 / (2 T) + T_gal = results[idx][name]["T"] + T_gal_err = results[idx][name]["T_err"] + if T_gal > 0: + r50_gal = SIGMA_TO_R50 * np.sqrt(T_gal / 2) + r50_gal_err = r50_gal * T_gal_err / (2 * T_gal) + else: + r50_gal = r50_gal_err = np.nan + output_dict[name]['r50'].append(r50_gal) + output_dict[name]['r50_err'].append(r50_gal_err) output_dict[name]['r50psf'].append(results[idx]["r50_PSFo"]) + output_dict[name]['r50psf_err'].append( + results[idx]["r50_err_PSFo"] + ) output_dict[name]["g1"].append(results[idx][name]["g"][0]) output_dict[name]["g2"].append(results[idx][name]["g"][1]) output_dict[name]["g1_err"].append( @@ -621,14 +628,17 @@ def process(self): 1 for k in ['noshear', '1p', '1m', '2p', '2m'] if res.get(k, {}).get('flags', 0) != 0 ) - r50_psfo = np.sqrt(max(psf_res['T_psf'], 0) / 2) + # PSF half-light radius r50 = sqrt(2 ln 2) * sigma with + # sigma = sqrt(T / 2); error from d sigma / d T = 1 / (4 sigma) + sigma_psfo = np.sqrt(max(psf_res['T_psf'], 0) / 2) res['g_PSFo'] = psf_res['g_psf'] res['g_err_PSFo'] = psf_res['g_psf_err'] res['T_PSFo'] = psf_res['T_psf'] res['T_err_PSFo'] = psf_res['T_psf_err'] - res['r50_PSFo'] = r50_psfo + res['r50_PSFo'] = SIGMA_TO_R50 * sigma_psfo res['r50_err_PSFo'] = ( - psf_res['T_psf_err'] / (2 * r50_psfo) if r50_psfo > 0 else np.nan + SIGMA_TO_R50 * psf_res['T_psf_err'] / (4 * sigma_psfo) + if sigma_psfo > 0 else np.nan ) res['mcal_flags'] = 0 final_res.append(res) @@ -838,7 +848,7 @@ def get_noise(gal, weight, guess, pixel_scale, thresh=1.2): Weight image guess : list Gaussian parameters fot the window function - ``[x0, y0, g1, g2, r50, flux]`` + ``[x0, y0, g1, g2, T, flux]`` pixel_scale : float Pixel scale of the galaxy image thresh : float, optional diff --git a/src/shapepipe/tests/test_ngmix.py b/src/shapepipe/tests/test_ngmix.py index aeb2274c7..c324dfa6e 100644 --- a/src/shapepipe/tests/test_ngmix.py +++ b/src/shapepipe/tests/test_ngmix.py @@ -94,3 +94,107 @@ def test_metacal_is_reproducible_with_fixed_seed(): irreproducible from one run to the next. """ npt.assert_array_equal(_metacal_noshear_g(42), _metacal_noshear_g(42)) + + +def _fake_metacal_result(T, T_err, T_psf, T_psf_err): + """Build one minimal metacal result as produced by the process loop. + + The PSF size entries mirror the conversion done at the source: + ``r50_PSFo = sqrt(2 ln 2) * sigma`` with ``sigma = sqrt(T_psf / 2)``. + """ + from shapepipe.modules.ngmix_package.ngmix import SIGMA_TO_R50 + + sigma_psf = np.sqrt(max(T_psf, 0) / 2) + per_type = { + "nfev": 1, + "g": [0.01, -0.02], + "g_cov": np.diag([1e-4, 1e-4]), + "T": T, + "T_err": T_err, + "flux": 100.0, + "flux_err": 1.0, + "s2n": 50.0, + "flags": 0, + } + res = { + "obj_id": 1, + "n_epoch_model": 1, + "moments_fail": 0, + "g_PSFo": [0.001, 0.002], + "g_err_PSFo": [1e-5, 1e-5], + "T_PSFo": T_psf, + "T_err_PSFo": T_psf_err, + "r50_PSFo": SIGMA_TO_R50 * sigma_psf, + "r50_err_PSFo": ( + SIGMA_TO_R50 * T_psf_err / (4 * sigma_psf) + if sigma_psf > 0 + else np.nan + ), + "mcal_flags": 0, + } + res.update( + {name: dict(per_type) for name in ("1m", "1p", "2m", "2p", "noshear")} + ) + return res + + +def test_compile_results_size_columns_are_half_light_radii(): + """Every r50 column is a true half-light radius, on both sides. + + Galaxy ``r50 = sqrt(ln 2 * T)`` (not the raw area ``pars[4]``), PSF + ``r50psf = sqrt(2 ln 2) * sigma_psf`` (not bare sigma), and the + redundant ``*_psfo_ngmix`` size duplicates are gone. + """ + from shapepipe.modules.ngmix_package.ngmix import Ngmix, SIGMA_TO_R50 + + T, T_err = 0.18, 0.02 + T_psf, T_psf_err = 0.09, 0.001 + + inst = object.__new__(Ngmix) + inst._zero_point = 30.0 + out = inst.compile_results([_fake_metacal_result(T, T_err, T_psf, T_psf_err)]) + + noshear = out["noshear"] + + # galaxy: r50 = sqrt(ln2 * T) = 1.17741 * sqrt(T / 2), error dT * r50/(2T) + r50_expected = np.sqrt(np.log(2) * T) + npt.assert_allclose(noshear["r50"], [r50_expected]) + npt.assert_allclose(noshear["r50_err"], [r50_expected * T_err / (2 * T)]) + + # PSF: r50psf = sqrt(2 ln2) * sigma, sigma = sqrt(T_psf / 2) + sigma_psf = np.sqrt(T_psf / 2) + npt.assert_allclose(noshear["r50psf"], [SIGMA_TO_R50 * sigma_psf]) + npt.assert_allclose( + noshear["r50psf_err"], [SIGMA_TO_R50 * T_psf_err / (4 * sigma_psf)] + ) + + # galaxy/PSF r50 are now commensurable: same convention on both sides + npt.assert_allclose( + np.array(noshear["r50"]) / np.array(noshear["r50psf"]), + [np.sqrt(T / T_psf)], + ) + + # areas pass through untouched, with the err asymmetry fixed + npt.assert_allclose(noshear["Tpsf"], [T_psf]) + npt.assert_allclose(noshear["Tpsf_err"], [T_psf_err]) + + # the *_psfo_ngmix size duplicates are retired + for retired in ( + "T_psfo_ngmix", + "T_err_psfo_ngmix", + "r50_psfo_ngmix", + "r50_err_psfo_ngmix", + ): + assert retired not in noshear + + +def test_compile_results_nonpositive_T_gives_nan_r50(): + """A non-positive fitted area cannot yield a real half-light radius.""" + from shapepipe.modules.ngmix_package.ngmix import Ngmix + + inst = object.__new__(Ngmix) + inst._zero_point = 30.0 + out = inst.compile_results([_fake_metacal_result(-0.05, 0.02, 0.09, 0.001)]) + + assert np.isnan(out["noshear"]["r50"]).all() + assert np.isnan(out["noshear"]["r50_err"]).all()