diff --git a/astra.yaml b/astra.yaml index 87c844909..b823f803b 100644 --- a/astra.yaml +++ b/astra.yaml @@ -1955,15 +1955,18 @@ analyses: rationale: >- [HARDCODED] detections with no ngmix row (no surviving epoch, or an exception during the fit) stay in the catalogue with sentinels: - sizes, fluxes, magnitudes and flags 0, flux and magnitude errors -1, + sizes, fluxes and magnitudes 0, flux and magnitude errors -1, ellipticities and their errors -10, size errors 1e30, - NGMIX_N_EPOCH 0. A failed object's NGMIX_MCAL_FLAGS reads 0, the - success value, so a cut on flags alone keeps it; NGMIX_N_EPOCH > 0 - removes it. + NGMIX_N_EPOCH 0. Their flag columns carry ngmix's own + LM_FUNC_NOTFINITE bit (2**12; FLAGS_ and NGMIX_MCAL_FLAGS) + and NGMIX_MCAL_TYPES_FAIL is 5, as for a metacal type that reports + success without a finite shear, so the MCAL_FLAGS == 0 and + MCAL_TYPES_FAIL == 0 cut rejects them; NGMIX_N_EPOCH == 0 tells + them apart from fitted objects. default: sentinel_values options: sentinel_values: - label: Keep with sentinels (flags 0, e -10, T_ERR 1e30) + label: Keep with sentinels (flags LM_FUNC_NOTFINITE, e -10, T_ERR 1e30) drop_unmatched: label: Drop objects without shapes excluded: true diff --git a/src/shapepipe/modules/make_cat_package/make_cat.py b/src/shapepipe/modules/make_cat_package/make_cat.py index ee545626a..18de1474e 100644 --- a/src/shapepipe/modules/make_cat_package/make_cat.py +++ b/src/shapepipe/modules/make_cat_package/make_cat.py @@ -15,6 +15,11 @@ from astropy.wcs import WCS from sqlitedict import SqliteDict +from shapepipe.modules.ngmix_package.ngmix import ( + get_mcal_flags, + get_mcal_types_fail, + get_type_flags, +) from shapepipe.pipeline import file_io from shapepipe.utilities import mask_query @@ -423,11 +428,16 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False): If True, write the parallel ``NGMIXm_*`` (moments-branch) columns. @sc [decision:catalogue_assembly.failure_sentinels,label:coupling] failure-sentinel-cut-semantics - Missing-row ``T``, ``SNR``, flux, magnitude, PSF-size and flag values - initialize to 0, which is in range and cannot identify failure. Use - ``NGMIX_N_EPOCH == 0`` for that cut; the -10 ellipticity, -1 - flux/magnitude-error and 1e30 size-error sentinels are out of range and - can also identify missing fits. + An object absent from the ngmix catalogue was never fit (no usable + epoch, or its fit raised). Its flag columns take what + :func:`ngmix.get_type_flags` derives for an empty metacal result: + ``LM_FUNC_NOTFINITE`` in ``FLAGS_`` and ``MCAL_FLAGS``, and + ``MCAL_TYPES_FAIL`` 5, never 0, so sp_validation's ``MCAL_FLAGS == 0`` + and ``MCAL_TYPES_FAIL == 0`` cut rejects it. Its ``T``, ``SNR``, flux, + magnitude and PSF-size values initialize to 0, which is in range and + cannot identify failure; ``NGMIX_N_EPOCH == 0`` identifies a + never-fit row, and the -10 ellipticity, -1 flux/magnitude-error and + 1e30 size-error sentinels are out of range. """ self._key_ends = ["1M", "1P", "2M", "2P", "NOSHEAR"] @@ -462,13 +472,19 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False): n_obj = len(self._obj_id) self._w_log.info(f"writing ngmix info for {n_obj} objects") + # An object ngmix never fit has no metacal result at all. + never_fit = {} + if moments: m = "m" else: m = "" self._add2dict("NGMIX_N_EPOCH", np.zeros(n_obj)) - self._add2dict("NGMIX_MCAL_TYPES_FAIL", np.zeros(n_obj)) + self._add2dict( + "NGMIX_MCAL_TYPES_FAIL", + np.full(n_obj, get_mcal_types_fail(never_fit), dtype=float), + ) self._add2dict("NGMIX_NEIGHBOUR_FLAG", np.zeros(n_obj)) prefix = f"NGMIX{m}" @@ -477,18 +493,22 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False): # reconvolution kernel (PSF_RECONV); see ngmix.average_original_psf / # average_multiepoch_psf for what each PSF family is. G1/G2 are scalar # reduced-shear components, not a 2-vector. Sentinels: - # sizes/fluxes/mags/flags 0, *_ERR fluxes/mags -1, ellipticities -10, - # *_ERR sizes 1e30. + # sizes/fluxes/mags 0, *_ERR fluxes/mags -1, ellipticities -10, + # *_ERR sizes 1e30; flags as ngmix derives for an empty metacal result + # (never_fit). for key_str in ( f"{prefix}_T_", f"{prefix}_SNR_", f"{prefix}_FLUX_", f"{prefix}_MAG_", - f"{prefix}_FLAGS_", f"{prefix}_T_PSF_ORIG_", f"{prefix}_T_PSF_RECONV_", ): self._update_dict(key_str, np.zeros(n_obj)) + self._update_dict( + f"{prefix}_FLAGS_", + np.full(n_obj, get_type_flags(never_fit), dtype=float), + ) for key_str in ( f"{prefix}_FLUX_ERR_", f"{prefix}_MAG_ERR_", @@ -515,7 +535,10 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False): f"{prefix}_T_ERR_PSF_RECONV_", ): self._update_dict(key_str, np.ones(n_obj) * 1e30) - self._add2dict(f"{prefix}_MCAL_FLAGS", np.zeros(n_obj)) + self._add2dict( + f"{prefix}_MCAL_FLAGS", + np.full(n_obj, get_mcal_flags(never_fit), dtype=float), + ) for idx, id_tmp in enumerate(self._obj_id): ind = np.where(id_tmp == ngmix_id)[0] diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 62502c4cd..707985b80 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -65,12 +65,57 @@ class MetacalResult(NamedTuple): orig: dict +def get_type_flags(fit): + """Get Type Flags. + + Fit flags of one metacal type, reading absence of evidence of success + as failure. + + @sc [label:convention] mcal-flags-zero-means-measured + A flag of 0 means the fit ran, reported success and returned a finite + shear; no default or fallback may produce 0. FLAGS_, MCAL_FLAGS + (OR) and MCAL_TYPES_FAIL (count) all derive from this function, and + sp_validation selects galaxies on MCAL_FLAGS == 0 and + MCAL_TYPES_FAIL == 0 as "measured". Every failure carries one of + ngmix's own bits (``ngmix.flags``); ShapePipe adds none: + + - the fitter reported failure: its own ``flags``, unchanged; + - the fitter reported success (``flags == 0``) without a finite shear + ``g``, or the result has no ``flags``, or the type is absent: + ``LM_FUNC_NOTFINITE`` (2**12). The shear the catalogue holds for + such a type is NaN or a sentinel, never a finite measurement. + + An object ngmix never fit (no usable epoch, or its fit raised) has no + metacal result at all; make_cat evaluates it as ``{}`` for every type, + so it carries ``LM_FUNC_NOTFINITE`` in every flag column and fails all + five types. ``NGMIX_N_EPOCH == 0`` tells it apart from a fitted object + whose LM fit set the same bit. + + Parameters + ---------- + fit : dict + One metacal type's fit result; ``{}`` when the type is absent. + + Returns + ------- + int + The fit's own ``flags``, or ``ngmix.flags.LM_FUNC_NOTFINITE`` when + the result is absent, lacks ``flags``, or claims success + (``flags == 0``) without a finite shear ``g``. + """ + notfinite = ngmix.flags.LM_FUNC_NOTFINITE + flags = int(fit.get('flags', notfinite)) + g = np.asarray(fit.get('g', (np.nan, np.nan)), dtype=float) + if flags == 0 and not np.all(np.isfinite(g)): + return notfinite + return flags + + def get_mcal_flags(res): """Get Metacal Flags. - Bitwise OR of the per-type metacal fit flags, the v1 contract for the - downstream NGMIX_MCAL_FLAGS column: nonzero whenever any metacal - type's galaxy fit failed. + Object-level metacal flags: the bitwise OR of :func:`get_type_flags` + over :data:`METACAL_TYPES` (the NGMIX_MCAL_FLAGS column). Parameters ---------- @@ -80,13 +125,259 @@ def get_mcal_flags(res): Returns ------- int - OR of all per-type ``flags``. + OR of all per-type flags; 0 only if every type was measured. """ return int(np.bitwise_or.reduce( - [res.get(name, {}).get('flags', 0) for name in METACAL_TYPES] + [get_type_flags(res.get(name, {})) for name in METACAL_TYPES] )) +def get_mcal_types_fail(res): + """Get Metacal Types Fail. + + Number of metacal types (0-5) with nonzero :func:`get_type_flags` (the + NGMIX_MCAL_TYPES_FAIL column). + + Parameters + ---------- + res : dict + MetacalBootstrapper result dict with one entry per metacal type. + + Returns + ------- + int + Count of failed metacal types. + """ + return sum( + get_type_flags(res.get(name, {})) != 0 for name in METACAL_TYPES + ) + + +def log_run_health(w_log, count, n_fitted, n_flagged): + """Log Run Health. + + Log an error when a run's metacal fits failed wholesale: either no + object fitted at all, or every fitted object carries nonzero + ``mcal_flags``. + + @sc [label:operations] run-health-logs-not-raises + Wholesale metacal failure is logged at error level, never raised: one + empty edge tile or broken input must not abort a multi-tile campaign + job, and the error line is the signal to catch in review. + + Parameters + ---------- + w_log : logging.Logger + Logging instance + count : int + Number of objects considered for fitting + n_fitted : int + Number of objects that were fitted (present in the results list) + n_flagged : int + Number of fitted objects whose ``mcal_flags`` ended up nonzero + + """ + if count > 0 and n_fitted == 0: + w_log.error( + f'ngmix: all {count} objects failed the metacal fit' + ' (0 fitted); writing an empty catalogue. Expected only for a' + ' tile with no usable epochs; otherwise check the vignettes,' + ' PSFs and ngmix installation.' + ) + if n_fitted > 0 and n_flagged == n_fitted: + w_log.error( + f'ngmix: 100% of {n_fitted} fitted objects carry nonzero' + ' mcal_flags -- the metacal fit failed wholesale; outputs' + ' are unusable.' + ) + + +def empty_metacal_output(): + """Empty Metacal Output. + + The five-HDU-shaped output dict :meth:`Ngmix.compile_results` returns + for zero fitted objects: one empty list per column, per metacal type. + Factored out so a tile with no measurable objects gets the identical + catalogue shape whether it is discovered by :meth:`Ngmix.process` + (partly-empty tile, one skipped object at a time) or by one of + ``ngmix_runner``'s all-empty-store guards (wholesale-empty tile, before + any object is read). + + Returns + ------- + dict + ``{metacal_type: {column: []}}``, matching an empty + :meth:`Ngmix.compile_results` call. + + Raises + ------ + ValueError + If the hardcoded HDU name list drifts out of sync with + :data:`METACAL_TYPES`. + """ + # Output HDU order. Same set as METACAL_TYPES, but kept in this + # fixed order so output catalogues stay byte-reproducible; the check + # below guards against the two lists silently diverging. + names = ["1m", "1p", "2m", "2p", "noshear"] + if set(names) != set(METACAL_TYPES): + raise ValueError( + "compile_results metacal type list is out of sync with" + + " METACAL_TYPES" + ) + names2 = [ + 'id', + 'n_epoch_model', + 'mcal_types_fail', + 'neighbour_flag', + 'nfev_fit', + # galaxy + 'g1', + 'g1_err', + 'g2', + 'g2_err', + 'T', + 'T_err', + 'flux', + 'flux_err', + 's2n', + 'mag', + 'mag_err', + 'flags', + 'mcal_flags', + # original image PSF (psfex/mccd), fit by average_original_psf + 'g1_psf_orig', + 'g2_psf_orig', + 'g1_err_psf_orig', + 'g2_err_psf_orig', + 'T_psf_orig', + 'T_err_psf_orig', + # metacal reconvolution kernel, fit by average_multiepoch_psf + 'g1_psf_reconv', + 'g2_psf_reconv', + 'g1_err_psf_reconv', + 'g2_err_psf_reconv', + 'T_psf_reconv', + 'T_err_psf_reconv', + ] + return {k: {kk: [] for kk in names2} for k in names} + + +def write_ngmix_fits(output_path, output_dict): + """Write Ngmix Fits. + + Write a compiled ngmix results dict to a fresh output FITS file, one + HDU per metacal type. The file must not already exist; an existing + output is appended to by :meth:`Ngmix.save_results`, not this function. + + Parameters + ---------- + output_path : str + Path of the FITS file to create + output_dict : dict + Compiled results, as returned by :meth:`Ngmix.compile_results` or + :func:`empty_metacal_output` + + Raises + ------ + IndexError + If ``output_dict`` does not have exactly five HDUs + """ + n_hdu = len(output_dict.keys()) + if n_hdu != 5: + raise IndexError( + f"FITS output file data has {n_hdu} HDUs," + + " expected are 5" + ) + f_out = file_io.FITSCatalogue( + output_path, open_mode=file_io.BaseCatalogue.OpenMode.ReadWrite + ) + for key in output_dict.keys(): + f_out.save_as_fits(output_dict[key], ext_name=key.upper()) + + +def write_empty_tile_output(output_dir, file_number_string, w_log, count): + """Write Empty Tile Output. + + Write ngmix's empty-tile product for a guard that fires before any + stamp is read: the same five-HDU empty catalogue and run-health error + line that :meth:`Ngmix.process` writes once every object in a tile has + been skipped, without constructing an ``Ngmix`` instance (which would + open the galaxy vignette store) or reading any vignette. + + Used by ``ngmix_runner``'s all-empty-store guards, which must return + before the galaxy vignette store is opened at all -- reading its + many-epoch, many-object arrays is what triggers a C-level malloc crash + when the store is large (see the PSF-empty guard's own comment). + + Parameters + ---------- + output_dir : str + Output directory + file_number_string : str + File numbering scheme + w_log : logging.Logger + Logging instance + count : int + Number of objects considered for fitting (see :func:`log_run_health`); + for these guards, the number of entries in the empty store. + """ + log_run_health(w_log, count, n_fitted=0, n_flagged=0) + output_path = f"{output_dir}/ngmix{file_number_string}.fits" + if not os.path.exists(output_path): + write_ngmix_fits(output_path, empty_metacal_output()) + + +def check_wcs_centroid_offset(centroid_source, tile_cat, gal_vign_cat): + """Check WCS Centroid Offset. + + Fail once, up front, when ``centroid_source="wcs"`` would have no + coadd-centroid offset to place the galaxy Jacobian at. + + @sc [label:coupling] wcs-centroid-needs-offset + ``centroid_source="wcs"`` reads the ``OFFSET`` the stamp extractor + (:func:`shapepipe.modules.vignetmaker_package.vignetmaker.get_stamps`) + writes into every vignette epoch entry; vignettes cut before that + extractor carry none. Left unchecked, :func:`make_ngmix_observation` + raises for every object in turn and :meth:`Ngmix.process`'s per-object + exception handling turns the whole tile into a silently empty + catalogue. OFFSET is a property of the extraction run, not of any one + object, so the first object with epochs speaks for the whole vignette + file: checking it is enough, and scanning every object would only cost + more sqlitedict unpickling for the same answer. + + Parameters + ---------- + centroid_source : {"wcs", "hsm"} + The configured centroid source; a no-op unless it is ``"wcs"``. + tile_cat : Tile_cat + Tile catalogue, read for its object ID order. + gal_vign_cat : Mapping + Galaxy vignette store, keyed by ``str(obj_id)``. + + Raises + ------ + ValueError + If ``centroid_source == "wcs"`` and the first object with epochs + has an epoch entry with no ``OFFSET``. + """ + if centroid_source != "wcs": + return + for obj_id in tile_cat.obj_id: + gal_obj = gal_vign_cat[str(obj_id)] + if gal_obj == 'empty' or not gal_obj: + continue + first_epoch = next(iter(gal_obj.values())) + if 'OFFSET' not in first_epoch: + raise ValueError( + "centroid_source='wcs' requires the coadd-centroid OFFSET" + " the stamp extractor writes into every vignette epoch," + " but this tile's vignettes carry none: re-extract the" + " stamps with the current vignetmaker, or set" + " centroid_source='hsm'." + ) + return + + def get_prior(pixel_scale, rng, T_range=None, F_range=None): """Build ngmix joint prior for a 6-parameter galaxy model. @@ -631,60 +922,26 @@ def compile_results(self, results): If SNR key not found """ - # Output HDU order. Same set as METACAL_TYPES, but kept in this - # fixed order so output catalogues stay byte-reproducible; the check - # below guards against the two lists silently diverging. - names = ["1m", "1p", "2m", "2p", "noshear"] - if set(names) != set(METACAL_TYPES): - raise ValueError( - "compile_results metacal type list is out of sync with" - + " METACAL_TYPES" - ) - names2 = [ - 'id', - 'n_epoch_model', - 'mcal_types_fail', - 'neighbour_flag', - 'nfev_fit', - # galaxy - 'g1', - 'g1_err', - 'g2', - 'g2_err', - 'T', - 'T_err', - 'flux', - 'flux_err', - 's2n', - 'mag', - 'mag_err', - 'flags', - 'mcal_flags', - # original image PSF (psfex/mccd), fit by average_original_psf - 'g1_psf_orig', - 'g2_psf_orig', - 'g1_err_psf_orig', - 'g2_err_psf_orig', - 'T_psf_orig', - 'T_err_psf_orig', - # metacal reconvolution kernel, fit by average_multiepoch_psf - 'g1_psf_reconv', - 'g2_psf_reconv', - 'g1_err_psf_reconv', - 'g2_err_psf_reconv', - 'T_psf_reconv', - 'T_err_psf_reconv', - ] - output_dict = {k: {kk: [] for kk in names2} for k in names} + # Column layout (HDU names and per-type columns) lives in + # empty_metacal_output, shared with the runner's all-empty-store + # guards so every zero-object catalogue has the identical shape. + output_dict = empty_metacal_output() + names = list(output_dict.keys()) for idx in range(len(results)): + # Object-level quality columns, derived from the same per-type + # flags as the ``flags`` column below (see get_type_flags). + mcal_flags = get_mcal_flags(results[idx]) + mcal_types_fail = get_mcal_types_fail(results[idx]) for name in names: - fit = results[idx][name] + fit = results[idx].get(name, {}) + flags = get_type_flags(fit) # ngmix 2.x does not raise on fit failure: after ntry the # result keeps flags != 0 and carries none of the # measurement keys (g, g_cov, T, T_err, flux, flux_err, - # s2n). NaN-fill those so failed types are recorded with - # their flags instead of crashing the tile on a KeyError. + # s2n). NaN-fill those (and an absent type) so failed types + # are recorded with their flags instead of crashing the tile + # on a KeyError. flux = fit.get("flux", np.nan) flux_err = fit.get("flux_err", np.nan) g = np.asarray(fit.get("g", (np.nan, np.nan))) @@ -701,9 +958,7 @@ def compile_results(self, results): output_dict[name]["n_epoch_model"].append( results[idx]["n_epoch_model"] ) - output_dict[name]["mcal_types_fail"].append( - results[idx]["mcal_types_fail"] - ) + output_dict[name]["mcal_types_fail"].append(mcal_types_fail) # Per-object blend flag (see process()); replicated across all # shear types like id / n_epoch_model / mcal_types_fail. output_dict[name]["neighbour_flag"].append( @@ -750,15 +1005,13 @@ def compile_results(self, results): output_dict[name]["s2n"].append(fit["s2n"]) elif "s2n_r" in fit: output_dict[name]["s2n"].append(fit["s2n_r"]) - elif fit["flags"] != 0: + elif flags != 0: output_dict[name]["s2n"].append(np.nan) else: raise KeyError("No SNR key (s2n, s2n_r) found in results") - output_dict[name]["flags"].append(fit["flags"]) - output_dict[name]["mcal_flags"].append( - results[idx].get("mcal_flags", 0) - ) + output_dict[name]["flags"].append(flags) + output_dict[name]["mcal_flags"].append(mcal_flags) return output_dict @@ -801,11 +1054,7 @@ def save_results(self, output_dict): output_name = self.get_output_path(self._output_dir) if not os.path.exists(output_name): - f_out = file_io.FITSCatalogue( - output_name, open_mode=file_io.BaseCatalogue.OpenMode.ReadWrite - ) - for key in output_dict.keys(): - f_out.save_as_fits(output_dict[key], ext_name=key.upper()) + write_ngmix_fits(output_name, output_dict) return with fits.open(output_name, mode='update') as hdul: @@ -978,11 +1227,22 @@ def process(self): dict Dictionary containing the NGMIX metacal results + Raises + ------ + ValueError + If ``centroid_source == "wcs"`` and the vignette catalogue + carries no coadd-centroid OFFSET (see + :func:`check_wcs_centroid_offset`). + @sc [decision:shape_measurement.fit_initialisation,decision:shape_measurement.ngmix_seed_mode] """ tile_cat = Tile_cat(self._tile_cat_path, self._seg_cat_path) vignet_cat = self._vignet_cat + check_wcs_centroid_offset( + self._centroid_source, tile_cat, vignet_cat.gal_vign_cat + ) + final_res = [] count = 0 @@ -990,6 +1250,7 @@ def process(self): n_no_epoch = 0 n_ngmix_fail = 0 n_fitted = 0 + n_flagged = 0 id_first = -1 id_last = -1 count_batch = 0 @@ -1009,8 +1270,12 @@ def process(self): # Read each store once here and pass the dicts down: every # sqlitedict access unpickles the object's whole all-epoch dict. psf_obj = vignet_cat.psf_vign_cat[str(obj_id)] + # Avoid allocating galaxy stamp arrays when there is no PSF coverage. + if psf_obj == 'empty' or not psf_obj: + n_empty_cat += 1 + continue gal_obj = vignet_cat.gal_vign_cat[str(obj_id)] - if psf_obj == 'empty' or gal_obj == 'empty': + if gal_obj == 'empty' or not gal_obj: n_empty_cat += 1 continue @@ -1086,15 +1351,10 @@ def process(self): # epochs that survived the PSF fit and entered the model, # not the number of epochs submitted (v1 contract) res['n_epoch_model'] = psf_res['n_epoch'] - # Count of metacal fit types (0-5) with nonzero fit flags. - # (In ngmix v1 the same-named column counted moments-initial-guess - # failures from get_guess, which no longer exists — hence the - # rename to mcal_types_fail / NGMIX_MCAL_TYPES_FAIL.) - res['mcal_types_fail'] = sum( - 1 for k in METACAL_TYPES - if res.get(k, {}).get('flags', 0) != 0 - ) - res['mcal_flags'] = get_mcal_flags(res) + # The mcal flag columns are derived from the per-type results in + # compile_results; here they only feed the run-health count. + if get_mcal_flags(res) != 0: + n_flagged += 1 # Two distinct PSF families (shapepipe#749), each carrying its own # ellipticity AND size: the metacal reconvolution kernel (psf_res) # and the original image PSF (psf_orig_res). Tag both into res from @@ -1139,6 +1399,8 @@ def process(self): + f" {n_fitted} fitted" ) + log_run_health(self._w_log, count, n_fitted, n_flagged) + vignet_cat.close() # Put all results together diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index 63f03872b..1509803ef 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -11,7 +11,7 @@ from sqlitedict import SqliteDict from shapepipe.modules.module_decorator import module_runner -from shapepipe.modules.ngmix_package.ngmix import Ngmix +from shapepipe.modules.ngmix_package.ngmix import Ngmix, write_empty_tile_output @module_runner( @@ -43,7 +43,18 @@ def ngmix_runner( w_log, ): """Define The Ngmix Runner. + @sc [decision:shape_measurement.blend_handling,decision:shape_measurement.centroid_source,decision:shape_measurement.defect_fill,decision:shape_measurement.metacal_scheme] + + @sc [label:operations] empty-tile-product + A tile whose PSF or galaxy vignette store is entirely empty never + reaches ``Ngmix``: an early guard in this runner writes the empty + catalogue and logs the zero-fitted run itself, before either store is + opened for its stamps. A tile that is only partly empty runs + ``Ngmix.process`` as usual, which skips each empty object individually + and writes the same empty-catalogue product if every object ends up + skipped. Either way, campaign completeness checks see one catalogue per + tile. """ # Read config file entries @@ -158,6 +169,9 @@ def ngmix_runner( f"All {len(psf_keys)} PSF vignet entries are empty in " f"{psf_vignet_path} — no PSF coverage for this tile. Skipping ngmix." ) + write_empty_tile_output( + run_dirs["output"], file_number_string, w_log, len(psf_keys) + ) return None, None # Check that image vignets are not all empty before initialising ngmix @@ -170,6 +184,9 @@ def ngmix_runner( f"All {len(keys)} image vignets are 'empty' in {image_vignet_path} " "— no valid CCD coverage for this tile. Skipping ngmix." ) + write_empty_tile_output( + run_dirs["output"], file_number_string, w_log, len(keys) + ) return None, None # Metacal reconvolution-kernel scheme (metacal_pars['psf']): "fitgauss" diff --git a/tests/module/test_make_cat.py b/tests/module/test_make_cat.py index 576ad2c95..520a82fa5 100644 --- a/tests/module/test_make_cat.py +++ b/tests/module/test_make_cat.py @@ -14,10 +14,14 @@ the pre-#749 code. """ +import functools +import operator + import numpy as np import numpy.testing as npt import pytest from astropy.io import fits +from ngmix.flags import LM_FUNC_NOTFINITE, NAME_MAP from sqlitedict import SqliteDict from shapepipe.modules.make_cat_package.make_cat import SaveCatalogue @@ -219,8 +223,9 @@ def test_save_ngmix_data_fills_sentinels_for_absent_objects(tmp_path): make_cat pre-fills every column with a type-specific sentinel and only overwrites the rows whose ``NUMBER`` matches an ngmix ``id``. An object SExtractor saw but ngmix never fit (no matching id) must therefore keep - the sentinels: 0 for sizes/flags, -10 for ellipticities, 1e30 for - ``T_ERR``, -1 for flux/mag errors. + the sentinels: 0 for sizes, -10 for ellipticities, 1e30 for ``T_ERR``, + -1 for flux/mag errors. (Its flag columns are pinned by + ``test_galaxy_cut_admits_only_measured_objects``.) """ ngmix_path = tmp_path / "ngmix-2.fits" # ngmix fit only object 22; the final cat also carries 11 and 99. @@ -255,6 +260,47 @@ def test_save_ngmix_data_fills_sentinels_for_absent_objects(tmp_path): npt.assert_allclose(n_epoch[absent], [0.0, 0.0]) +def _metacal_result(obj_id): + """One clean object's result as ``Ngmix.process`` hands it on. + + Every metacal type carries a successful fit; the PSF families are + object-level, copied from ``_ngmix_row``. + """ + row = _ngmix_row(obj_id) + fit = { + "nfev": row["nfev_fit"], + "g": [row["g1"], row["g2"]], + "g_cov": np.diag([row["g1_err"] ** 2, row["g2_err"] ** 2]), + "T": row["T"], "T_err": row["T_err"], + "flux": row["flux"], "flux_err": row["flux_err"], + "s2n": row["s2n"], "flags": 0, + } + res = { + "obj_id": obj_id, + "n_epoch_model": row["n_epoch_model"], + "neighbour_flag": row["neighbour_flag"], + } + for key in NGMIX_KEYS: + if key.endswith("_psf_orig") or key.endswith("_psf_reconv"): + res[key] = row[key] + res.update({name: dict(fit) for name in SHEAR_EXTS_LOWER}) + return res + + +def _serialise_then_merge(tmp_path, results, cat_ids): + """Write ``results`` with ngmix's own write path, read via make_cat. + + ``compile_results`` + ``save_results`` produce the ngmix catalogue; + ``_save_ngmix_data`` merges it into a final catalogue of ``cat_ids``. + """ + ngmix_inst = object.__new__(Ngmix) + ngmix_inst._zero_point = 30.0 + ngmix_inst._output_dir = str(tmp_path) + ngmix_inst._file_number_string = "-0" + ngmix_inst.save_results(ngmix_inst.compile_results(results)) + return _run_save_ngmix(ngmix_inst.get_output_path(str(tmp_path)), cat_ids) + + def test_save_ngmix_data_matches_module_serialised_catalogue(tmp_path): """End-to-end: make_cat reads a catalogue ngmix itself serialised. @@ -264,39 +310,9 @@ def test_save_ngmix_data_matches_module_serialised_catalogue(tmp_path): drift in the ngmix key set: the keys must line up end to end. """ obj_ids = [5, 7] - results = [] - for oid in obj_ids: - row = _ngmix_row(oid) - per_type = { - "nfev": row["nfev_fit"], - "g": [row["g1"], row["g2"]], - "g_cov": np.diag([row["g1_err"] ** 2, row["g2_err"] ** 2]), - "T": row["T"], "T_err": row["T_err"], - "flux": row["flux"], "flux_err": row["flux_err"], - "s2n": row["s2n"], "flags": 0, - } - res = { - "obj_id": oid, - "n_epoch_model": row["n_epoch_model"], - "mcal_types_fail": row["mcal_types_fail"], - "neighbour_flag": row["neighbour_flag"], - "mcal_flags": row["mcal_flags"], - } - for key in NGMIX_KEYS: - if key.endswith("_psf_orig") or key.endswith("_psf_reconv"): - res[key] = row[key] - res.update({name: dict(per_type) for name in SHEAR_EXTS_LOWER}) - results.append(res) - - ngmix_inst = object.__new__(Ngmix) - ngmix_inst._zero_point = 30.0 - ngmix_inst._output_dir = str(tmp_path) - ngmix_inst._file_number_string = "-0" - out_dict = ngmix_inst.compile_results(results) - ngmix_inst.save_results(out_dict) - - ngmix_path = ngmix_inst.get_output_path(str(tmp_path)) - out = _run_save_ngmix(ngmix_path, obj_ids) + out = _serialise_then_merge( + tmp_path, [_metacal_result(oid) for oid in obj_ids], obj_ids + ) # Both PSF families survive the round trip, distinct, on every object. npt.assert_allclose( @@ -456,6 +472,81 @@ def test_save_psf_data_fills_sentinel_for_absent_epochs(tmp_path): assert out[col][1] == -1, col +@pytest.mark.parametrize("shear", SHEAR_EXTS) +@pytest.mark.parametrize("component", [0, 1], ids=["g1", "g2"]) +@pytest.mark.parametrize("nonfinite", [np.nan, np.inf, -np.inf]) +def test_galaxy_cut_rejects_each_nonfinite_shear_component( + tmp_path, shear, component, nonfinite, +): + """A non-finite component in any metacal type cannot pass the galaxy cut.""" + result = _metacal_result(1) + result[shear.lower()]["g"] = [0.1, 0.2] + result[shear.lower()]["g"][component] = nonfinite + out = _serialise_then_merge(tmp_path, [result], np.array([1])) + + assert out["NGMIX_MCAL_FLAGS"][0] == LM_FUNC_NOTFINITE + assert out["NGMIX_MCAL_TYPES_FAIL"][0] == 1 + assert out[f"NGMIX_FLAGS_{shear}"][0] == LM_FUNC_NOTFINITE + + +def test_galaxy_cut_admits_only_measured_objects(tmp_path): + """Contracts mcal-flags-zero-means-measured and failure-sentinel-cut-semantics. + + Consumer-side invariant through the real write path: synthetic metacal + results -> ``compile_results`` / ``save_results`` -> ``_save_ngmix_data`` + -> sp_validation's galaxy cut ``MCAL_FLAGS == 0 & MCAL_TYPES_FAIL == 0``. + One object per way a fit can fail to be a measurement, plus objects + ngmix never fit. Only the two clean objects may pass, and every passing + row must carry a fitted shape in all five metacal types. + """ + nan, inf = float("nan"), float("inf") + results = {oid: _metacal_result(oid) for oid in range(1, 9)} + # 1, 8: clean. + results[2]["1p"] = {"flags": 0x8, "nfev": 5} # fitter reported failure + del results[3]["noshear"] # type absent from the result + del results[4]["2m"]["flags"] # no flags key, shape present + results[5]["noshear"]["g"] = [nan, nan] # flags 0, non-finite shear + results[6]["1m"]["g"] = [inf, 0.1] # flags 0, infinite shear + del results[7]["2p"]["g"] # flags 0, no shear at all + # 101, 102: in the final catalogue but never fit by ngmix. + cat_ids = np.array([101, 1, 2, 3, 102, 4, 5, 6, 7, 8]) + + out = _serialise_then_merge(tmp_path, list(results.values()), cat_ids) + + mcal_flags = np.asarray(out["NGMIX_MCAL_FLAGS"]).astype(np.int64) + types_fail = np.asarray(out["NGMIX_MCAL_TYPES_FAIL"]).astype(np.int64) + passed = (mcal_flags == 0) & (types_fail == 0) + + assert set(cat_ids[passed]) == {1, 8}, ( + "mcal-flags-zero-means-measured / failure-sentinel-cut-semantics: the cut" + f" admitted {sorted(set(cat_ids[passed]) - {1, 8})}" + ) + assert np.all(np.asarray(out["NGMIX_N_EPOCH"])[passed] > 0) + for shear in SHEAR_EXTS: + for comp in ("G1", "G2"): + shape = np.asarray(out[f"NGMIX_{comp}_{shear}"])[passed] + assert np.all(np.isfinite(shape) & (shape != -10.0)), ( + f"NGMIX_{comp}_{shear}: a row passing the cut has no fitted shape" + ) + + # The three flag columns agree on every row, fitted or never fit: + # MCAL_FLAGS is the OR and MCAL_TYPES_FAIL the count of FLAGS_. + type_flags = np.array( + [np.asarray(out[f"NGMIX_FLAGS_{shear}"]) for shear in SHEAR_EXTS] + ).astype(np.int64) + npt.assert_array_equal(mcal_flags, np.bitwise_or.reduce(type_flags)) + npt.assert_array_equal(types_fail, np.count_nonzero(type_flags, axis=0)) + + # Failures carry ngmix's own bits: the fitter's flags pass through, and + # every no-finite-shear case, never-fit objects included, reads + # LM_FUNC_NOTFINITE. No bit outside ngmix.flags is ever set. + expected = {oid: LM_FUNC_NOTFINITE for oid in (3, 4, 5, 6, 7, 101, 102)} + expected.update({1: 0, 8: 0, 2: 0x8}) + npt.assert_array_equal(mcal_flags, [expected[oid] for oid in cat_ids]) + ngmix_bits = functools.reduce(operator.or_, NAME_MAP) + assert not np.any(type_flags & ~ngmix_bits) + + # --- _save_psf_data: fixed per-epoch slot count (N_EPOCH_SLOTS) --- # Per-family empty-slot sentinel: what a slot holds when no epoch fills it. diff --git a/tests/module/test_ngmix.py b/tests/module/test_ngmix.py index 17ff692ba..ed87b4be2 100644 --- a/tests/module/test_ngmix.py +++ b/tests/module/test_ngmix.py @@ -512,7 +512,10 @@ def test_get_mcal_flags_ors_per_type_fit_flags(): """ from shapepipe.modules.ngmix_package.ngmix import get_mcal_flags - res = {name: {"flags": 0} for name in ("noshear", "1p", "1m", "2p", "2m")} + res = { + name: {"flags": 0, "g": [0.01, -0.02]} + for name in ("noshear", "1p", "1m", "2p", "2m") + } assert get_mcal_flags(res) == 0 res["1p"]["flags"] = 0x8 @@ -520,6 +523,241 @@ def test_get_mcal_flags_ors_per_type_fit_flags(): assert get_mcal_flags(res) == 0xA +class _RecordingLogger: + """Records ``error`` calls; drops ``info`` and ``warning``.""" + + def __init__(self): + self.errors = [] + + def error(self, msg): + self.errors.append(msg) + + def info(self, *_args, **_kwargs): + pass + + def warning(self, *_args, **_kwargs): + pass + + +@pytest.mark.parametrize( + "count, n_fitted, n_flagged, n_errors", + [ + (10, 0, 0, 1), # nothing fitted: an error line, not an exception + (10, 10, 10, 1), # every fit flagged: an error line + (10, 9, 1, 0), # healthy run: no false alarm + ], +) +def test_log_run_health_logs_wholesale_failure_without_raising( + count, n_fitted, n_flagged, n_errors +): + """Contract run-health-logs-not-raises. + + Failure modes: raising (one empty tile aborts a campaign job), staying + silent on a wholesale failure, and alarming on a healthy run. + """ + from shapepipe.modules.ngmix_package.ngmix import log_run_health + + w_log = _RecordingLogger() + log_run_health(w_log, count=count, n_fitted=n_fitted, n_flagged=n_flagged) + + assert len(w_log.errors) == n_errors + + +def test_process_survives_a_tile_with_nothing_to_fit(tmp_path): + """Contract run-health-logs-not-raises, through ``Ngmix.process``. + + A tile whose every object has no stamps (an empty edge tile) fits + nothing. ``process`` must log the 0-fitted error and still write its + (empty) catalogue, not abort the campaign job. + """ + n_obj = 3 + tile_cat = tmp_path / "tile_cat.fits" + objects = fits.BinTableHDU.from_columns( + [ + fits.Column(name="NUMBER", format="J", array=np.arange(1, n_obj + 1)), + fits.Column(name="XWIN_WORLD", format="D", array=np.zeros(n_obj)), + fits.Column(name="YWIN_WORLD", format="D", array=np.zeros(n_obj)), + ], + name="LDAC_OBJECTS", + ) + imhead = fits.BinTableHDU.from_columns( + [fits.Column(name="Field Header Card", format="1A", array=["x"])], + name="LDAC_IMHEAD", + ) + fits.HDUList([fits.PrimaryHDU(), imhead, objects]).writeto(tile_cat) + + sqlite_paths = [] + for name in ("gal", "bkg", "psf", "weight", "flag", "headers"): + path = str(tmp_path / f"{name}.sqlite") + db = SqliteDict(path) + if name in ("gal", "psf"): + for obj_id in range(1, n_obj + 1): + db[str(obj_id)] = "empty" + db.commit() + db.close() + sqlite_paths.append(path) + + w_log = _RecordingLogger() + ngmix = Ngmix( + [str(tile_cat)] + sqlite_paths[:5], + str(tmp_path), + "-001-001", + 30.0, + 0.186, + sqlite_paths[5], + w_log, + ) + ngmix.process() + + assert any("0 fitted" in msg for msg in w_log.errors) + with fits.open(ngmix.get_output_path(str(tmp_path))) as hdul: + assert len(hdul["NOSHEAR"].data) == 0 + + +@pytest.mark.parametrize("flags", [(8, 8), (8, 0), (0, 8)]) +def test_process_counts_flagged_fits_across_batches(tmp_path, monkeypatch, flags): + """Contract run-health-logs-not-raises includes every fitted batch. + + Known fitter outcomes feed the real process loop and FITS writer. Only + an all-flagged run should log an error, even when each fit is saved in + its own batch and no results remain in memory at the end. + """ + from types import SimpleNamespace + from shapepipe.modules.ngmix_package import ngmix as module + + tile = SimpleNamespace(obj_id=[1, 2], flux=None, seg=None) + galaxies = {str(i): {"exp-1": {"OFFSET": [0., 0.]}} for i in tile.obj_id} + stamp = SimpleNamespace( + gals=[np.ones((5, 5))], ra=[42.], dec=[30.], ccd=20, + ) + psf = dict( + n_epoch=1, g_psf=[.01, -.01], g_psf_err=[.001, .001], + T_psf=.1, T_psf_err=.01, + ) + results = [] + for flag in flags: + result = _fake_metacal_result(.18, .02, .09, .001) + result["1p"]["flags"] = flag + results.append((result, psf, psf)) + fits_to_return = iter(results) + monkeypatch.setattr(module, "Tile_cat", lambda *args: tile) + monkeypatch.setattr(module, "prepare_postage_stamps", lambda *args: stamp) + monkeypatch.setattr( + module, "do_ngmix_metacal", lambda *args, **kwargs: next(fits_to_return), + ) + inst = object.__new__(Ngmix) + inst._tile_cat_path = "in-memory-tile" + inst._seg_cat_path = None + inst._vignet_cat = SimpleNamespace( + gal_vign_cat=galaxies, psf_vign_cat=galaxies, close=lambda: None, + ) + inst._centroid_source = "wcs" + inst._id_obj_min = inst._id_obj_max = -1 + inst._bkg_sub = True + inst._pixel_scale = .186 + inst._blend_handling = "noisefill" + inst._dilate_neighbour = 1 + inst._metacal_psf = "fitgauss" + inst._save_batch = 1 + inst._zero_point = 30. + inst._output_dir = str(tmp_path) + inst._file_number_string = "-001-001" + inst._w_log = _RecordingLogger() + + inst.process() + + with fits.open(inst.get_output_path(str(tmp_path))) as hdul: + npt.assert_array_equal(hdul["NOSHEAR"].data["mcal_flags"], flags) + if all(flags): + assert len(inst._w_log.errors) == 1 + assert "100% of 2 fitted objects carry nonzero mcal_flags" in inst._w_log.errors[0] + else: + assert inst._w_log.errors == [] + + +def _write_tile_cat_with_one_object(tmp_path): + """A tile catalogue with a single object, for the OFFSET-check tests.""" + tile_cat = tmp_path / "tile_cat.fits" + objects = fits.BinTableHDU.from_columns( + [ + fits.Column(name="NUMBER", format="J", array=np.array([1])), + fits.Column(name="XWIN_WORLD", format="D", array=np.zeros(1)), + fits.Column(name="YWIN_WORLD", format="D", array=np.zeros(1)), + ], + name="LDAC_OBJECTS", + ) + imhead = fits.BinTableHDU.from_columns( + [fits.Column(name="Field Header Card", format="1A", array=["x"])], + name="LDAC_IMHEAD", + ) + fits.HDUList([fits.PrimaryHDU(), imhead, objects]).writeto(tile_cat) + return tile_cat + + +def _ngmix_with_offsetless_vignette(tmp_path, centroid_source): + """One object whose galaxy vignette has an epoch entry with no OFFSET. + Its PSF is marked 'empty', so a run that reaches the per-object loop + skips the object cleanly instead of failing for some other reason. + """ + tile_cat = _write_tile_cat_with_one_object(tmp_path) + + sqlite_paths = [] + for name in ("gal", "bkg", "psf", "weight", "flag", "headers"): + path = str(tmp_path / f"{name}.sqlite") + db = SqliteDict(path) + if name == "gal": + db["1"] = {"expA-1": {"VIGNET": np.ones((5, 5))}} + if name == "psf": + db["1"] = "empty" + db.commit() + db.close() + sqlite_paths.append(path) + + return Ngmix( + [str(tile_cat)] + sqlite_paths[:5], + str(tmp_path), + "-001-001", + 30.0, + 0.186, + sqlite_paths[5], + _RecordingLogger(), + centroid_source=centroid_source, + ) + + +def test_process_raises_before_any_fit_when_wcs_offset_is_missing(tmp_path): + """Contract wcs-centroid-needs-offset, through ``Ngmix.process``. + + Vignettes cut before the OFFSET-writing stamp extractor carry no + OFFSET. Under the default ``centroid_source="wcs"`` this must raise + once, up front -- not disappear into the per-object try/except that + would otherwise turn a wholesale failure into a silently empty + catalogue. The fixture's PSF store marks the object 'empty', so if the + check did not run first, ``process`` would simply skip the object and + return cleanly rather than raising at all. + """ + ngmix = _ngmix_with_offsetless_vignette(tmp_path, centroid_source="wcs") + + with pytest.raises(ValueError, match="OFFSET"): + ngmix.process() + + +def test_process_ignores_missing_offset_under_hsm(tmp_path): + """Contract wcs-centroid-needs-offset: a no-op under ``centroid_source= + "hsm"``, which never reads OFFSET. + + The same offset-less vignette that raises under "wcs" must not trip the + check under "hsm": the object is simply skipped (its PSF is marked + 'empty') and the run completes. + """ + ngmix = _ngmix_with_offsetless_vignette(tmp_path, centroid_source="hsm") + + ngmix.process() + + with fits.open(ngmix.get_output_path(str(tmp_path))) as hdul: + assert len(hdul["NOSHEAR"].data) == 0 + + def test_average_multiepoch_psf_skips_failed_psf_epochs(): """A failed-PSF epoch must be skipped, not KeyError the whole object. diff --git a/tests/module/test_ngmix_empty_tile.py b/tests/module/test_ngmix_empty_tile.py new file mode 100644 index 000000000..c581b8970 --- /dev/null +++ b/tests/module/test_ngmix_empty_tile.py @@ -0,0 +1,177 @@ +"""Empty edge tiles publish the same catalogue product as populated tiles.""" + +import subprocess +import sys + +from astropy.io import fits +import numpy as np +import pytest +from sqlitedict import SqliteDict + +from shapepipe.modules.ngmix_runner import ngmix_runner +from shapepipe.pipeline.config import CustomParser + + +@pytest.mark.parametrize("empty_store", ["galaxy", "psf"]) +def test_empty_tile_cli_writes_catalogue_and_health_log(tmp_path, empty_store): + """Contract empty-tile-product through the shipped launcher. + + A successful exit without a catalogue fails the campaign's per-tile + completeness check. Both empty-galaxy and empty-PSF stores must publish + one catalogue with five empty metacal HDUs and log the zero-fitted run. + """ + inputs = tmp_path / "input" + outputs = tmp_path / "output" + inputs.mkdir() + outputs.mkdir() + ids = np.arange(1, 4) + objects = fits.BinTableHDU.from_columns( + [ + fits.Column(name="NUMBER", format="J", array=ids), + fits.Column(name="XWIN_WORLD", format="D", array=np.zeros(3)), + fits.Column(name="YWIN_WORLD", format="D", array=np.zeros(3)), + ], + name="LDAC_OBJECTS", + ) + imhead = fits.BinTableHDU.from_columns( + [fits.Column(name="Field Header Card", format="1A", array=["x"])], + name="LDAC_IMHEAD", + ) + fits.HDUList([fits.PrimaryHDU(), imhead, objects]).writeto( + inputs / "tile_sexcat-001-001.fits" + ) + epoch = {"exp-1": {"VIGNET": np.ones((5, 5)), "OFFSET": [0.0, 0.0]}} + stores = { + "image": "empty" if empty_store == "galaxy" else epoch, + "galaxy_psf": {} if empty_store == "psf" else epoch, + "exp_background": {}, + "weight": {}, + "flag": {}, + "log_exp_headers": {}, + } + for name, value in stores.items(): + with SqliteDict(str(inputs / f"{name}-001-001.sqlite")) as db: + for obj_id in ids: + db[str(obj_id)] = value + db.commit() + + config = tmp_path / "config.ini" + config.write_text(f"""[DEFAULT] +RUN_NAME = run_empty +RUN_DATETIME = False +VERBOSE = False +[EXECUTION] +MODULE = ngmix_runner +MODE = smp +[FILE] +INPUT_DIR = {inputs} +OUTPUT_DIR = {outputs} +NUMBERING_SCHEME = -000-000 +[JOB] +SMP_BATCH_SIZE = 1 +TIMEOUT = 00:02:00 +[NGMIX_RUNNER] +MAG_ZP = 30 +PIXEL_SCALE = .186 +ID_OBJ_MIN = -1 +ID_OBJ_MAX = -1 +""") + run = subprocess.run( + [sys.executable, "-m", "shapepipe.shapepipe_run", "-c", str(config)], + cwd=tmp_path, capture_output=True, text=True, timeout=120, + ) + assert run.returncode == 0, run.stdout + run.stderr + product_dir = outputs / "run_empty" / "ngmix_runner" / "output" + products = list(product_dir.iterdir()) + assert [path.name for path in products] == ["ngmix-001-001.fits"] + with fits.open(products[0]) as hdul: + assert {hdu.name for hdu in hdul[1:]} == { + "NOSHEAR", "1P", "1M", "2P", "2M", + } + assert all(len(hdu.data) == 0 for hdu in hdul[1:]) + log = "\n".join(path.read_text() for path in outputs.rglob("process*.log")) + assert log.count("all 3 objects failed the metacal fit (0 fitted)") == 1 + + +class _RecordingLogger: + """Minimal logger that records messages passed to warning/error/info.""" + + def __init__(self): + self.messages = [] + + def info(self, msg, *_args, **_kwargs): + self.messages.append(msg) + + def warning(self, msg, *_args, **_kwargs): + self.messages.append(msg) + + def error(self, msg, *_args, **_kwargs): + self.messages.append(msg) + + +def _runner_config(): + config = CustomParser() + config.read_string( + "[NGMIX_RUNNER]\n" + "MAG_ZP = 30\n" + "PIXEL_SCALE = .186\n" + "ID_OBJ_MIN = -1\n" + "ID_OBJ_MAX = -1\n" + ) + return config + + +def test_psf_all_empty_guard_never_opens_galaxy_store(tmp_path): + """The PSF-all-empty guard returns before the galaxy vignette store is + opened at all -- the crash it guards against comes from unpickling that + store's many-epoch arrays, so the check must fire without ever touching + it, not just without acting on its contents. + + The galaxy ("image") store here is a corrupted, non-sqlite file: + ``ngmix_runner`` must still succeed. Moving the PSF-empty guard after the + image-vignet read -- or dropping it -- makes this test fail with a + ``sqlite3`` error instead of the assertions below (checked by hand: with + the two guards' order swapped, this test goes red). + """ + ids = [1, 2, 3] + psf_path = tmp_path / "galaxy_psf-001-001.sqlite" + with SqliteDict(str(psf_path)) as db: + for obj_id in ids: + db[str(obj_id)] = {} + db.commit() + + image_path = tmp_path / "image-001-001.sqlite" + image_path.write_bytes(b"not a sqlite database") + + outputs = tmp_path / "output" + outputs.mkdir() + log = _RecordingLogger() + + input_file_list = [ + "tile_sexcat-001-001.fits", # never opened: guard fires before Tile_cat + str(image_path), # galaxy store: corrupted, must stay unopened + "exp_background-001-001.sqlite", + str(psf_path), + "weight-001-001.sqlite", + "flag-001-001.sqlite", + "log_exp_headers-001-001.sqlite", + ] + result = ngmix_runner( + input_file_list, + {"output": str(outputs)}, + "-001-001", + _runner_config(), + "NGMIX_RUNNER", + log, + ) + + assert result == (None, None) + with fits.open(outputs / "ngmix-001-001.fits") as hdul: + assert {hdu.name for hdu in hdul[1:]} == { + "NOSHEAR", "1P", "1M", "2P", "2M", + } + assert all(len(hdu.data) == 0 for hdu in hdul[1:]) + assert any( + "all 3 objects failed the metacal fit (0 fitted)" in msg + for msg in log.messages + ) diff --git a/tests/module/test_psf_averaging_properties.py b/tests/module/test_psf_averaging_properties.py index a6ce2b782..d22435719 100644 --- a/tests/module/test_psf_averaging_properties.py +++ b/tests/module/test_psf_averaging_properties.py @@ -24,9 +24,13 @@ from astropy.io import fits from hypothesis import given, settings from hypothesis import strategies as st +from ngmix.flags import LM_FUNC_NOTFINITE from shapepipe.modules.make_cat_package.make_cat import SaveCatalogue -from shapepipe.modules.ngmix_package.ngmix import _average_psf_fits +from shapepipe.modules.ngmix_package.ngmix import ( + METACAL_TYPES, + _average_psf_fits, +) # --------------------------------------------------------------------------- # @@ -235,7 +239,7 @@ def test_all_epochs_failed_raises_zero_division(flagged_specs): _SENTINELS = { "NGMIX_T_NOSHEAR": 0.0, "NGMIX_SNR_NOSHEAR": 0.0, - "NGMIX_FLAGS_NOSHEAR": 0.0, + "NGMIX_FLAGS_NOSHEAR": LM_FUNC_NOTFINITE, "NGMIX_T_PSF_ORIG_NOSHEAR": 0.0, "NGMIX_T_PSF_RECONV_NOSHEAR": 0.0, "NGMIX_FLUX_ERR_NOSHEAR": -1.0, @@ -248,8 +252,11 @@ def test_all_epochs_failed_raises_zero_division(flagged_specs): "NGMIX_T_ERR_PSF_ORIG_NOSHEAR": 1e30, "NGMIX_T_ERR_PSF_RECONV_NOSHEAR": 1e30, "NGMIX_N_EPOCH": 0.0, - "NGMIX_MCAL_FLAGS": 0.0, - "NGMIX_MCAL_TYPES_FAIL": 0.0, + # Never fit reads as an empty metacal result, not a clean fit + # (contract failure-sentinel-cut-semantics): LM_FUNC_NOTFINITE, all types + # failed. + "NGMIX_MCAL_FLAGS": LM_FUNC_NOTFINITE, + "NGMIX_MCAL_TYPES_FAIL": len(METACAL_TYPES), "NGMIX_NEIGHBOUR_FLAG": 0.0, } @@ -277,8 +284,14 @@ def info(self, *_args, **_kwargs): def _measured_row(obj_id): - """One fit object whose every value is far from any sentinel (5 / 0.5).""" + """One fit object whose every value is far from any sentinel (5 / 0.5). + + ``mcal_types_fail`` is the one int column whose sentinel (absent -> + ``len(METACAL_TYPES)`` == 5, shapepipe#889) coincides with the generic + placeholder, so it gets its own distinct in-range value. + """ row = {key: (5 if key in _INT_KEYS else 0.5) for key in _NGMIX_KEYS} + row["mcal_types_fail"] = 2 row["id"] = obj_id return row