From 300b27b6324fb73d825d563dcbef92e11ce7a9cd Mon Sep 17 00:00:00 2001 From: Jim Bosch Date: Fri, 25 Aug 2023 12:13:29 -0400 Subject: [PATCH 01/18] Add first-draft shell of metadetection PipelineTask. --- python/lsst/drp/tasks/metadetection_shear.py | 286 +++++++++++++++++++ 1 file changed, 286 insertions(+) create mode 100644 python/lsst/drp/tasks/metadetection_shear.py diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py new file mode 100644 index 00000000..16219665 --- /dev/null +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -0,0 +1,286 @@ +# This file is part of drp_tasks. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +from __future__ import annotations + +from lsst.daf.butler import DataCoordinate, DatasetRef +from lsst.pipe.base.connectionTypes import BaseInput, Output + +__all__ = () + +from collections.abc import Collection, Mapping, Sequence +from typing import Any, ClassVar + +import pyarrow as pa + +from lsst.pex.config import ListField +from lsst.pipe.base import ( + InputQuantizedConnection, + NoWorkFound, + OutputQuantizedConnection, + QuantumContext, + PipelineTask, + PipelineTaskConfig, + PipelineTaskConnections, + Struct, +) +import lsst.pipe.base.connectionTypes as cT +from lsst.cell_coadds import MultipleCellCoadd, SingleCellCoadd + + +class MetadetectionShearConnections(PipelineTaskConnections, dimensions={"patch"}): + """Definitions of inputs and outputs for MetadetectionShearTask.""" + + input_coadds = cT.Input( + "DeepCoadd", + storageClass="MultipleCellCoadd", + doc="Per-band deep coadds.", + multiple=True, + dimensions={"patch", "band"}, + ) + + # TODO: If there are image-like or other non-catalog output products (e.g. + # detection masks), add them here. + + object_catalog = cT.Output( + "ShearObject", + storageClass="ArrowTable", + doc="Output catalog with all quantities measured inside the metacalibration loop.", + multiple=False, + dimensions={"patch"}, + ) + object_schema = cT.InitOutput( + "ShearObject_schema", + # TODO: It's not currently possible to save ArrowSchema objects on + # their own, but some combination of Eli and Jim can figure out how to + # fix that. + storageClass="ArrowSchema", + doc="Schema of the output catalog.", + ) + + # TODO: if we want a per-cell output catalog instead of just denormalizing + # everything into per-object catalogs, add it and its schema here. + + config: MetadetectionShearConfig + + def adjustQuantum( + self, + inputs: dict[str, tuple[BaseInput, Collection[DatasetRef]]], + outputs: dict[str, tuple[Output, Collection[DatasetRef]]], + label: str, + data_id: DataCoordinate, + ) -> tuple[ + Mapping[str, tuple[BaseInput, Collection[DatasetRef]]], + Mapping[str, tuple[Output, Collection[DatasetRef]]], + ]: + # Docstring inherited. + # This is a hook for customizing what is input and output to each + # invocation of the task as early as possible, which we override here + # to make sure we have exactly the required bands, no more, no less. + connection, original_input_coadds = inputs["input_coadds"] + bands_missing = set(self.config.required_bands) + adjusted_input_coadds = [] + for ref in original_input_coadds: + if ref.dataId["band"] in bands_missing: + adjusted_input_coadds.append(ref) + bands_missing.remove(ref.dataId["band"]) + if bands_missing: + raise NoWorkFound( + f"Required bands {bands_missing} not present for {label}@{data_id})." + ) + adjusted_inputs = {"input_coadds": (connection, adjusted_input_coadds)} + inputs.update(adjusted_inputs) + super().adjustQuantum(inputs, outputs, label, data_id) + return adjusted_inputs, {} + + +class MetadetectionShearConfig( + PipelineTaskConfig, pipelineConnections=MetadetectionShearConnections +): + """Configuration definition for MetadetectionShearTask.""" + + required_bands = ListField[str]( + "Bands expected to be present. Cells with one or more of these bands " + "missing will be skipped. Bands other than those listed here will " + "not be processed.", + default=["g", "r", "i", "z"], + optional=False, + ) + + # TODO: expose more configuration options here. + + +class MetadetectionShearTask(PipelineTask): + """A PipelineTask that measures shear using metadetection.""" + + _DefaultName: ClassVar[str] = "metadetectionShear" + ConfigClass: ClassVar[type[MetadetectionShearConfig]] = MetadetectionShearConfig + + config: MetadetectionShearConfig + + def __init__(self, *, initInputs: dict[str, Any] | None = None, **kwargs: Any): + super().__init__(initInputs=initInputs, **kwargs) + self.object_schema = self.make_object_schema(self.config) + + @classmethod + def make_object_schema(cls, config: MetadetectionShearConfig) -> pa.Schema: + """Construct a PyArrow Schema for this task's main output catalog. + + Parameters + ---------- + config : `MetadetectionShearConfig` + Configuration that may be used to control details of the schema. + + Returns + ------- + object_schema : `pyarrow.Schema` + Schema for the object catalog produced by this task. Each field's + metadata should include both a 'doc' entry and a 'unit' entry. + """ + return pa.schema( + [ + pa.field( + "id", + pa.uint64(), + nullable=False, + metadata={ + "doc": ( + "Unique identifier for a ShearObject, specific " + "to a single metacalibration counterfactual image." + ), + "unit": "", + }, + ), + pa.field( + "tract", + pa.uint64(), + nullable=False, + metadata={ + "doc": "ID of the tract on which this measurement was made.", + "unit": "", + }, + ), + pa.field( + "patch_x", + pa.uint64(), + nullable=False, + metadata={ + "doc": "Column within the tract of the patch on which this measurement was made.", + "unit": "", + }, + ), + pa.field( + "patch_y", + pa.uint64(), + nullable=False, + metadata={ + "doc": "Row within the tract of the patch on which this measurement was made.", + "unit": "", + }, + ), + pa.field( + "cell_x", + pa.uint64(), + nullable=False, + metadata={ + "doc": "Column within the patch of the cell on which this measurement was made.", + "unit": "", + }, + ), + pa.field( + "cell_y", + pa.uint64(), + nullable=False, + metadata={ + "doc": "Row within the patch of the cell on which this measurement was made.", + "unit": "", + }, + ), + # TODO: add more field definitions here + ] + ) + + def runQuantum( + self, + qc: QuantumContext, + inputRefs: InputQuantizedConnection, + outputRefs: OutputQuantizedConnection, + ) -> None: + # Docstring inherited. + # Read the coadds and put them in the order defined by + # config.required_bands (note that each MultipleCellCoadd object also + # knows its own band, if that's needed). + coadds_by_band = { + ref.dataId["band"]: qc.get(ref) for ref in inputRefs.input_coadds + } + outputs = self.run([coadds_by_band[b] for b in self.config.required_bands]) + qc.put(outputs, outputRefs) + + def run(self, patch_coadds: Sequence[MultipleCellCoadd]) -> Struct: + """Run metadetection on a patch. + + Parameters + ---------- + patch_coadds : `~collections.abc.Sequence` [ \ + `~lsst.cell_coadds.MultipleCellCoadd` ] + Per-band, per-patch coadds, in the order specified by + `MetadetectionShearConfig.required_bands`. + + Returns + ------- + results : `lsst.pipe.base.Struct` + Structure with the following attributes: + + - ``object_catalog`` [ `pyarrow.Table` ]: the output object + catalog for the patch, with schema equal to `object_schema`. + """ + single_cell_tables: list[pa.Table] = [] + for single_cell_coadds in zip(*patch_coadds, strict=True): + single_cell_tables.append(self.process_cell(single_cell_coadds)) + # TODO: if we need to do any cell-overlap-region deduplication here + # (instead of purely in science analysis code), this is where'd it'd + # happen. + return Struct( + object_catalog=pa.concat_tables(single_cell_tables), + ) + + def process_cell(self, single_cell_coadds: Sequence[SingleCellCoadd]) -> pa.Table: + """Run metadetection on a single cell. + + Parameters + ---------- + single_cell_coadds : `~collections.abc.Sequence` [ \ + `~lsst.cell_coadds.SingleCellCoadd` ] + Per-band, per-cell coadds, in the order specified by + `MetadetectionShearConfig.required_bands`. + + Returns + ------- + object_catalog : `pyarrow.Table` + Output object catalog for the cell, with schema equal to + `object_schema`. + """ + rows: list[dict[str, Any]] = [] + # TODO: run metadetection on the cell, filling in 'rows' with + # measurements. Or replace 'rows' with a 'columns' dict of numpy array + # columns and call 'from_pydict' instead of 'from_pylist' below, if + # that's more convenient. + return pa.Table.from_pylist(rows, self.object_schema) From a25beb144b2adeee13e85bb547e65e72a7c51d35 Mon Sep 17 00:00:00 2001 From: Jim Bosch Date: Fri, 25 Aug 2023 14:03:27 -0400 Subject: [PATCH 02/18] Fix zip-iteration over cells. MultipleCellCoadd isn't iterable, but its .cells attribute is a mapping. --- python/lsst/drp/tasks/metadetection_shear.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index 16219665..a3dccbec 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -253,7 +253,9 @@ def run(self, patch_coadds: Sequence[MultipleCellCoadd]) -> Struct: catalog for the patch, with schema equal to `object_schema`. """ single_cell_tables: list[pa.Table] = [] - for single_cell_coadds in zip(*patch_coadds, strict=True): + for single_cell_coadds in zip( + *[patch_coadd.cells.values() for patch_coadd in patch_coadds], strict=True + ): single_cell_tables.append(self.process_cell(single_cell_coadds)) # TODO: if we need to do any cell-overlap-region deduplication here # (instead of purely in science analysis code), this is where'd it'd From a51cb94d0b1d87ace02601b80cec8e2c17e48bb5 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Mon, 28 Aug 2023 10:32:44 -0700 Subject: [PATCH 03/18] correct name "cellCoadd" and temp fixes The Input name is set correctly to cellCoadd temporarily commented object schema "cT.InitOutput" because this doesn't work yet temporarilyl set required_bands to r, need to learn how to set with config --- python/lsst/drp/tasks/metadetection_shear.py | 23 +++++++++++--------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index a3dccbec..ed03a477 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -50,7 +50,7 @@ class MetadetectionShearConnections(PipelineTaskConnections, dimensions={"patch" """Definitions of inputs and outputs for MetadetectionShearTask.""" input_coadds = cT.Input( - "DeepCoadd", + "cellCoadd", storageClass="MultipleCellCoadd", doc="Per-band deep coadds.", multiple=True, @@ -67,14 +67,15 @@ class MetadetectionShearConnections(PipelineTaskConnections, dimensions={"patch" multiple=False, dimensions={"patch"}, ) - object_schema = cT.InitOutput( - "ShearObject_schema", - # TODO: It's not currently possible to save ArrowSchema objects on - # their own, but some combination of Eli and Jim can figure out how to - # fix that. - storageClass="ArrowSchema", - doc="Schema of the output catalog.", - ) + # TODO make this exist + # object_schema = cT.InitOutput( + # "ShearObject_schema", + # # TODO: It's not currently possible to save ArrowSchema objects on + # # their own, but some combination of Eli and Jim can figure out how to + # # fix that. + # storageClass="ArrowSchema", + # doc="Schema of the output catalog.", + # ) # TODO: if we want a per-cell output catalog instead of just denormalizing # everything into per-object catalogs, add it and its schema here. @@ -121,7 +122,9 @@ class MetadetectionShearConfig( "Bands expected to be present. Cells with one or more of these bands " "missing will be skipped. Bands other than those listed here will " "not be processed.", - default=["g", "r", "i", "z"], + # TODO learn how to set in a config file + # default=["g", "r", "i", "z"], + default=["r"], optional=False, ) From 29d1169e9f9f306bf6552173e6c3f214d3ba0c93 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Mon, 28 Aug 2023 10:42:24 -0700 Subject: [PATCH 04/18] satisfy the flake8 --- python/lsst/drp/tasks/metadetection_shear.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index ed03a477..34401eb4 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -71,8 +71,8 @@ class MetadetectionShearConnections(PipelineTaskConnections, dimensions={"patch" # object_schema = cT.InitOutput( # "ShearObject_schema", # # TODO: It's not currently possible to save ArrowSchema objects on - # # their own, but some combination of Eli and Jim can figure out how to - # # fix that. + # # their own, but some combination of Eli and Jim can figure out + # # how to # fix that. # storageClass="ArrowSchema", # doc="Schema of the output catalog.", # ) From 790bf5f2b9d336f62c7e8123fef18b04a71f991e Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Wed, 30 Aug 2023 10:33:24 -0700 Subject: [PATCH 05/18] running a sim to get something we can process Currently the real cell coadds don't have everything we need --- python/lsst/drp/tasks/metadetection_shear.py | 299 ++++++++++++++++++- 1 file changed, 293 insertions(+), 6 deletions(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index 34401eb4..9653bf44 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -20,6 +20,7 @@ # along with this program. If not, see . from __future__ import annotations +import numpy as np from lsst.daf.butler import DataCoordinate, DatasetRef from lsst.pipe.base.connectionTypes import BaseInput, Output @@ -118,6 +119,8 @@ class MetadetectionShearConfig( ): """Configuration definition for MetadetectionShearTask.""" + from lsst.meas.base import SkyMapIdGeneratorConfig + required_bands = ListField[str]( "Bands expected to be present. Cells with one or more of these bands " "missing will be skipped. Bands other than those listed here will " @@ -128,6 +131,8 @@ class MetadetectionShearConfig( optional=False, ) + idGenerator = SkyMapIdGeneratorConfig.make_field() + # TODO: expose more configuration options here. @@ -231,13 +236,20 @@ def runQuantum( # Read the coadds and put them in the order defined by # config.required_bands (note that each MultipleCellCoadd object also # knows its own band, if that's needed). + + # TODO make this work + idGenerator = self.config.idGenerator.apply(qc.quantum.dataId) + seed = idGenerator.catalog_id + coadds_by_band = { ref.dataId["band"]: qc.get(ref) for ref in inputRefs.input_coadds } - outputs = self.run([coadds_by_band[b] for b in self.config.required_bands]) + outputs = self.run( + [coadds_by_band[b] for b in self.config.required_bands], seed, + ) qc.put(outputs, outputRefs) - def run(self, patch_coadds: Sequence[MultipleCellCoadd]) -> Struct: + def run(self, patch_coadds: Sequence[MultipleCellCoadd], seed: int) -> Struct: """Run metadetection on a patch. Parameters @@ -246,6 +258,9 @@ def run(self, patch_coadds: Sequence[MultipleCellCoadd]) -> Struct: `~lsst.cell_coadds.MultipleCellCoadd` ] Per-band, per-patch coadds, in the order specified by `MetadetectionShearConfig.required_bands`. + seed: int + A seed for random number generator, used for simulation, and/or + fitting. Returns ------- @@ -255,11 +270,28 @@ def run(self, patch_coadds: Sequence[MultipleCellCoadd]) -> Struct: - ``object_catalog`` [ `pyarrow.Table` ]: the output object catalog for the patch, with schema equal to `object_schema`. """ + + # TODO figure out how to use a config to switch to simulation mode + # in that case also need to fake the id info + rng = np.random.RandomState(seed) + idstart = 0 + single_cell_tables: list[pa.Table] = [] for single_cell_coadds in zip( *[patch_coadd.cells.values() for patch_coadd in patch_coadds], strict=True ): - single_cell_tables.append(self.process_cell(single_cell_coadds)) + self._log_ids(single_cell_coadds[0]) + + res = self.process_cell(single_cell_coadds, rng, simulate=True) + # TODO figure out how to make the object ids + if len(res) > 0: + res['id'] = np.arange(idstart, idstart + len(res)) + da = _dictify(res) + table = pa.Table.from_pydict(da, self.object_schema) + + single_cell_tables.append(table) + idstart += len(res) + # TODO: if we need to do any cell-overlap-region deduplication here # (instead of purely in science analysis code), this is where'd it'd # happen. @@ -267,7 +299,10 @@ def run(self, patch_coadds: Sequence[MultipleCellCoadd]) -> Struct: object_catalog=pa.concat_tables(single_cell_tables), ) - def process_cell(self, single_cell_coadds: Sequence[SingleCellCoadd]) -> pa.Table: + def process_cell( + self, single_cell_coadds: Sequence[SingleCellCoadd], rng, + simulate=False, + ) -> pa.Table: """Run metadetection on a single cell. Parameters @@ -277,15 +312,267 @@ def process_cell(self, single_cell_coadds: Sequence[SingleCellCoadd]) -> pa.Tabl Per-band, per-cell coadds, in the order specified by `MetadetectionShearConfig.required_bands`. + rng: `np.random.RandomState` + Random number generator. + simulate: bool, optional + If set to True, simulate the data + Returns ------- object_catalog : `pyarrow.Table` Output object catalog for the cell, with schema equal to `object_schema`. """ - rows: list[dict[str, Any]] = [] + from metadetect.lsst.metadetect import run_metadetect + from metadetect.lsst.metadetect import get_config as get_mdet_config + + # rows: list[dict[str, Any]] = [] # TODO: run metadetection on the cell, filling in 'rows' with # measurements. Or replace 'rows' with a 'columns' dict of numpy array # columns and call 'from_pydict' instead of 'from_pylist' below, if # that's more convenient. - return pa.Table.from_pylist(rows, self.object_schema) + + coadd_data = _simulate_coadd(rng) + + mask_frac = _get_mask_frac( + coadd_data['mfrac_mbexp'], + trim_pixels=0, + ) + + mdet_config = get_mdet_config() + # TODO how to get the id/tract/patch_x etc. currently in schema + res = run_metadetect( + rng=rng, + config=mdet_config, + **coadd_data, + ) + comb_res = _make_comb_data( + single_cell_coadd=single_cell_coadds[0], + res=res, + meas_type=mdet_config['meas_type'], + mask_frac=mask_frac, + full_output=True, + ) + + return comb_res + if len(comb_res) == 0: + return pa.Table.from_pylist(comb_res, self.object_schema) + else: + da = _dictify(comb_res) + return pa.Table.from_pydict(da, self.object_schema) + + def _log_ids(self, single_cell_coadd): + idinfo = single_cell_coadd.identifiers + tract = idinfo.tract + patch_x = idinfo.patch.x + patch_y = idinfo.patch.y + cell_x = idinfo.cell.x + cell_y = idinfo.cell.y + mess = ( + f'tract: {tract} patch: {patch_x},{patch_y} cell: {cell_x},{cell_y}' + ) + self.log.info(mess) + + +def _dictify(data): + output = {} + for name in data.dtype.names: + output[name] = data[name] + return output + + +def _make_comb_data( + single_cell_coadd, + res, + meas_type, + mask_frac, + full_output=False, +): + import esutil as eu + + idinfo = single_cell_coadd.identifiers + + add_dt = [ + ('id', 'u8'), + ('tract', 'u4'), + ('patch_x', 'u1'), + ('patch_y', 'u1'), + ('cell_x', 'u1'), + ('cell_y', 'u1'), + ('shear_type', 'U2'), + ('mask_frac', 'f4'), + ('primary', bool), + ] + + if not hasattr(res, 'keys'): + res = {'noshear': res} + + dlist = [] + for stype in res.keys(): + data = res[stype] + if data is not None: + + if not full_output: + data = _trim_output_columns(data, meas_type) + + newdata = eu.numpy_util.add_fields(data, add_dt) + + newdata['tract'] = idinfo.tract + newdata['patch_x'] = idinfo.patch.x + newdata['patch_y'] = idinfo.patch.y + newdata['cell_x'] = idinfo.cell.x + newdata['cell_y'] = idinfo.cell.y + + if stype == 'noshear': + newdata['shear_type'] = 'ns' + else: + newdata['shear_type'] = stype + + dlist.append(newdata) + + if len(dlist) > 0: + output = eu.numpy_util.combine_arrlist(dlist) + # output['id'] = np.arange(output.size) + else: + output = [] + + return output + + +def _trim_output_columns(data, meas_type): + # TODO decide what to keep + raise NotImplementedError('implement trim output') + + # if meas_type == 'admom': + # meas_type = 'am' + # + # # note the bmask/ormask compress to nothing + # cols2keep_orig = [ + # get_flags_name(data=data, meas_type=meas_type), + # 'bmask', + # 'ormask', + # 'row', 'row0', + # 'col', 'col0', + # 'mfrac', + # '%s_s2n' % meas_type, + # '%s_T_ratio' % meas_type, + # '%s_g' % meas_type, + # '%s_g_cov' % meas_type, + # ] + # + # cols2keep = [] + # for col in cols2keep_orig: + # if col in data.dtype.names: + # cols2keep.append(col) + # + # return eu.numpy_util.extract_fields(data, cols2keep) + + +def _get_mask_frac(mfrac_mbexp, trim_pixels=0): + """ + get the average mask frac for each band and then return the max of those + """ + + mask_fracs = [] + for mfrac_exp in mfrac_mbexp: + mfrac = mfrac_exp.image.array + dim = mfrac.shape[0] + mfrac = mfrac[ + trim_pixels:dim - trim_pixels - 1, + trim_pixels:dim - trim_pixels - 1, + ] + mask_fracs.append(mfrac.mean()) + + return max(mask_fracs) + + +def _simulate_coadd(rng): + from descwl_shear_sims.sim import ( + make_sim, + get_sim_config, + get_se_dim, + ) + from descwl_shear_sims.galaxies import make_galaxy_catalog + from descwl_shear_sims.psfs import make_fixed_psf, make_ps_psf, make_rand_psf + # from descwl_shear_sims.stars import make_star_catalog + from descwl_coadd.coadd_nowarp import make_coadd_nowarp + from metadetect.lsst.util import extract_multiband_coadd_data + + g1, g2 = 0.02, 0.00 + + sim_config = get_sim_config() + if sim_config['se_dim'] is None: + sim_config['se_dim'] = get_se_dim( + coadd_dim=sim_config['coadd_dim'], + dither=sim_config['dither'], + rotate=sim_config['rotate'], + ) + + if sim_config['gal_type'] != 'wldeblend': + gal_config = sim_config.get('gal_config', None) + else: + sim_config["layout"] = None + gal_config = None + + galaxy_catalog = make_galaxy_catalog( + rng=rng, + gal_type=sim_config['gal_type'], + coadd_dim=sim_config['coadd_dim'], + buff=sim_config['buff'], + layout=sim_config['layout'], + sep=sim_config['sep'], # for layout='pair' + gal_config=gal_config, + ) + if sim_config['psf_type'] == 'ps': + psf = make_ps_psf( + rng=rng, + dim=sim_config['se_dim'], + variation_factor=sim_config['psf_variation_factor'], + ) + elif sim_config['randomize_psf']: + psf = make_rand_psf( + psf_type=sim_config["psf_type"], rng=rng, + ) + else: + psf = make_fixed_psf(psf_type=sim_config["psf_type"]) + + sim_data = make_sim( + rng=rng, + galaxy_catalog=galaxy_catalog, + coadd_dim=sim_config['coadd_dim'], + se_dim=sim_config['se_dim'], + g1=g1, + g2=g2, + psf=psf, + draw_stars=sim_config['draw_stars'], + psf_dim=sim_config['psf_dim'], + dither=sim_config['dither'], + rotate=sim_config['rotate'], + bands=sim_config['bands'], + epochs_per_band=sim_config['epochs_per_band'], + noise_factor=sim_config['noise_factor'], + cosmic_rays=sim_config['cosmic_rays'], + bad_columns=sim_config['bad_columns'], + star_bleeds=sim_config['star_bleeds'], + sky_n_sigma=sim_config['sky_n_sigma'], + ) + + bands = list(sim_data['band_data'].keys()) + exps = sim_data['band_data'][bands[0]] + assert ( + not sim_config['dither'] + and not sim_config['rotate'] + and len(exps) == 1 + ) + + coadd_data_list = [ + make_coadd_nowarp( + exp=exps[0], + psf_dims=sim_data['psf_dims'], + rng=rng, + remove_poisson=False, + ) + for band in bands + ] + + return extract_multiband_coadd_data(coadd_data_list) From 8f18935011b425750c44ecdc7174ec017ebd9710 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Wed, 30 Aug 2023 11:54:20 -0700 Subject: [PATCH 06/18] add more schema --- python/lsst/drp/tasks/metadetection_shear.py | 326 ++++++++++++++++++- 1 file changed, 319 insertions(+), 7 deletions(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index 9653bf44..f89615c5 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -222,7 +222,300 @@ def make_object_schema(cls, config: MetadetectionShearConfig) -> pa.Schema: "unit": "", }, ), - # TODO: add more field definitions here + + pa.field( + "wmom_flags", + pa.uint32(), + nullable=False, + metadata={ + "doc": "Overall flags for wmom measurement.", + "unit": "", + }, + ), + + # Original PSF measurements + pa.field( + "psfrec_flags", + pa.uint32(), + nullable=False, + metadata={ + "doc": "Flags for admom PSF measurement.", + "unit": "", + }, + ), + pa.field( + "psfrec_g_1", + pa.float32(), + nullable=False, + metadata={ + "doc": "admom g1 measurement for PSF.", + "unit": "", + }, + ), + pa.field( + "psfrec_g_2", + pa.float32(), + nullable=False, + metadata={ + "doc": "admom g2 measurement for PSF.", + "unit": "", + }, + ), + pa.field( + "psfrec_T", + pa.float32(), + nullable=False, + metadata={ + "doc": "admom wmom T ( + ) measurement for PSF.", + "unit": "", + }, + ), + + # reconvolved PSF measurements + pa.field( + "wmom_psf_flags", + pa.uint32(), + nullable=False, + metadata={ + "doc": "Flags for wmom reconvolved PSF measurement.", + "unit": "", + }, + ), + pa.field( + "wmom_psf_g_1", + pa.float32(), + nullable=False, + metadata={ + "doc": "wmom g1 measurement for reconvolved PSF.", + "unit": "", + }, + ), + pa.field( + "wmom_psf_g_2", + pa.float32(), + nullable=False, + metadata={ + "doc": "wmom g2 measurement for reconvolved PSF.", + "unit": "", + }, + ), + pa.field( + "wmom_psf_T", + pa.float32(), + nullable=False, + metadata={ + "doc": "wmom T ( + ) measurement for reconvolved PSF.", + "unit": "", + }, + ), + + # Object measurements + pa.field( + "wmom_obj_flags", + pa.uint32(), + nullable=False, + metadata={ + "doc": "Flags for wmom object measurement.", + "unit": "", + }, + ), + pa.field( + "wmom_s2n", + pa.float32(), + nullable=False, + metadata={ + "doc": "wmom object s2n measurement.", + "unit": "", + }, + ), + + pa.field( + "wmom_g_1", + pa.float32(), + nullable=False, + metadata={ + "doc": "wmom object g1 measurement.", + "unit": "", + }, + ), + pa.field( + "wmom_g_2", + pa.float32(), + nullable=False, + metadata={ + "doc": "wmom object g2 measurement.", + "unit": "", + }, + ), + + pa.field( + "wmom_T_flags", + pa.uint32(), + nullable=False, + metadata={ + "doc": "Flags for wmom T measurement for object.", + "unit": "", + }, + ), + + pa.field( + "wmom_T", + pa.float32(), + nullable=False, + metadata={ + "doc": "wmom T ( + ) measurement for object.", + "unit": "", + }, + ), + pa.field( + "wmom_T_err", + pa.float32(), + nullable=False, + metadata={ + "doc": "wmom T uncertainty for object.", + "unit": "", + }, + ), + pa.field( + "wmom_T_ratio", + pa.float32(), + nullable=False, + metadata={ + "doc": "wmom T/Tpsf for object.", + "unit": "", + }, + ), + + pa.field( + "wmom_band_flux_flags_1", + pa.uint32(), + nullable=False, + metadata={ + "doc": "wmom flux measurement flags for object.", + "unit": "", + }, + ), + pa.field( + "wmom_band_flux_1", + pa.float32(), + nullable=False, + metadata={ + "doc": "wmom flux for object.", + "unit": "", + }, + ), + pa.field( + "wmom_band_flux_err_1", + pa.float32(), + nullable=False, + metadata={ + "doc": "wmom flux uncertainty for object.", + "unit": "", + }, + ), + # pa.field( + # "row0", + # pa.float32(), + # nullable=False, + # metadata={ + # "doc": "row start for stamp.", + # "unit": "", + # }, + # ), + # pa.field( + # "col0", + # pa.float32(), + # nullable=False, + # metadata={ + # "doc": "column start for stamp.", + # "unit": "", + # }, + # ), + pa.field( + "row", + pa.float32(), + nullable=False, + metadata={ + "doc": "detected row for object, within image.", + "unit": "", + }, + ), + pa.field( + "col", + pa.float32(), + nullable=False, + metadata={ + "doc": "detected column for object, within image.", + "unit": "", + }, + ), + pa.field( + "row_diff", + pa.float32(), + nullable=False, + metadata={ + "doc": "difference of measured row from detected row.", + "unit": "", + }, + ), + pa.field( + "col_diff", + pa.float32(), + nullable=False, + metadata={ + "doc": "difference of measured column from detected row.", + "unit": "", + }, + ), + + pa.field( + "ra", + pa.float64(), + nullable=False, + metadata={ + "doc": "detected ra for object.", + "unit": "degrees", + }, + ), + pa.field( + "dec", + pa.float64(), + nullable=False, + metadata={ + "doc": "detected dec for object.", + "unit": "degrees", + }, + ), + + pa.field( + "bmask", + pa.uint32(), + nullable=False, + metadata={ + "doc": "bmask flags for object", + "unit": "", + }, + ), + pa.field( + "ormask", + pa.uint32(), + nullable=False, + metadata={ + "doc": "ored mask flags for object", + "unit": "", + }, + ), + + pa.field( + "mfrac", + pa.float32(), + nullable=False, + metadata={ + "doc": "gaussian weighted masked fraction for object.", + "unit": "", + }, + ), + ] ) @@ -283,6 +576,7 @@ def run(self, patch_coadds: Sequence[MultipleCellCoadd], seed: int) -> Struct: self._log_ids(single_cell_coadds[0]) res = self.process_cell(single_cell_coadds, rng, simulate=True) + # TODO figure out how to make the object ids if len(res) > 0: res['id'] = np.arange(idstart, idstart + len(res)) @@ -291,6 +585,7 @@ def run(self, patch_coadds: Sequence[MultipleCellCoadd], seed: int) -> Struct: single_cell_tables.append(table) idstart += len(res) + break # TODO: if we need to do any cell-overlap-region deduplication here # (instead of purely in science analysis code), this is where'd it'd @@ -355,11 +650,6 @@ def process_cell( ) return comb_res - if len(comb_res) == 0: - return pa.Table.from_pylist(comb_res, self.object_schema) - else: - da = _dictify(comb_res) - return pa.Table.from_pydict(da, self.object_schema) def _log_ids(self, single_cell_coadd): idinfo = single_cell_coadd.identifiers @@ -392,6 +682,19 @@ def _make_comb_data( idinfo = single_cell_coadd.identifiers + copy_dt = [ + # we will copy out of arrays to these + ('psfrec_g_1', 'f4'), + ('psfrec_g_2', 'f4'), + ('wmom_psf_g_1', 'f4'), + ('wmom_psf_g_2', 'f4'), + ('wmom_g_1', 'f4'), + ('wmom_g_2', 'f4'), + ('wmom_band_flux_flags_1', 'i4'), + ('wmom_band_flux_1', 'f4'), + ('wmom_band_flux_err_1', 'f4'), + ] + add_dt = [ ('id', 'u8'), ('tract', 'u4'), @@ -402,7 +705,7 @@ def _make_comb_data( ('shear_type', 'U2'), ('mask_frac', 'f4'), ('primary', bool), - ] + ] + copy_dt if not hasattr(res, 'keys'): res = {'noshear': res} @@ -416,6 +719,15 @@ def _make_comb_data( data = _trim_output_columns(data, meas_type) newdata = eu.numpy_util.add_fields(data, add_dt) + newdata['psfrec_g_1'] = newdata['psfrec_g'][:, 0] + newdata['psfrec_g_2'] = newdata['psfrec_g'][:, 1] + newdata['wmom_psf_g_1'] = newdata['wmom_psf_g'][:, 0] + newdata['wmom_psf_g_2'] = newdata['wmom_psf_g'][:, 1] + newdata['wmom_g_1'] = newdata['wmom_g'][:, 0] + newdata['wmom_g_2'] = newdata['wmom_g'][:, 1] + newdata['wmom_band_flux_flags_1'] = newdata['wmom_band_flux_flags'] + newdata['wmom_band_flux_1'] = newdata['wmom_band_flux'] + newdata['wmom_band_flux_err_1'] = newdata['wmom_band_flux_err'] newdata['tract'] = idinfo.tract newdata['patch_x'] = idinfo.patch.x From 6364e095dca52f70be093953d71fb14cd0992b23 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Wed, 30 Aug 2023 12:13:26 -0700 Subject: [PATCH 07/18] add shear_bands --- python/lsst/drp/tasks/metadetection_shear.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index f89615c5..1be96524 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -413,6 +413,16 @@ def make_object_schema(cls, config: MetadetectionShearConfig) -> pa.Schema: "unit": "", }, ), + pa.field( + "shear_bands", + pa.string(), + nullable=False, + metadata={ + "doc": "bands used for shear measurement.", + "unit": "", + }, + ), + # pa.field( # "row0", # pa.float32(), From 9f184a2e58b46afd0c0a20ff68a7b754364f90fb Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Thu, 31 Aug 2023 09:14:49 -0700 Subject: [PATCH 08/18] add code to convert cell coadd data type The metadetect code expects MultiBandExposure --- python/lsst/drp/tasks/metadetection_shear.py | 120 ++++++++++++++++--- 1 file changed, 104 insertions(+), 16 deletions(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index 1be96524..dfdb9d05 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -46,6 +46,11 @@ import lsst.pipe.base.connectionTypes as cT from lsst.cell_coadds import MultipleCellCoadd, SingleCellCoadd +from lsst.afw.image import ImageF, MaskedImageF, ExposureF, FilterLabel +from lsst.afw.math import FixedKernel +from lsst.meas.algorithms import KernelPsf +from metadetect.lsst.util import extract_multiband_coadd_data + class MetadetectionShearConnections(PipelineTaskConnections, dimensions={"patch"}): """Definitions of inputs and outputs for MetadetectionShearTask.""" @@ -576,16 +581,16 @@ def run(self, patch_coadds: Sequence[MultipleCellCoadd], seed: int) -> Struct: # TODO figure out how to use a config to switch to simulation mode # in that case also need to fake the id info - rng = np.random.RandomState(seed) + self.rng = np.random.RandomState(seed) idstart = 0 single_cell_tables: list[pa.Table] = [] - for single_cell_coadds in zip( + for cell_coadds in zip( *[patch_coadd.cells.values() for patch_coadd in patch_coadds], strict=True ): - self._log_ids(single_cell_coadds[0]) + self._log_ids(cell_coadds[0]) - res = self.process_cell(single_cell_coadds, rng, simulate=True) + res = self.process_cell(cell_coadds, simulate=False) # TODO figure out how to make the object ids if len(res) > 0: @@ -605,20 +610,18 @@ def run(self, patch_coadds: Sequence[MultipleCellCoadd], seed: int) -> Struct: ) def process_cell( - self, single_cell_coadds: Sequence[SingleCellCoadd], rng, + self, cell_coadds: Sequence[SingleCellCoadd], simulate=False, ) -> pa.Table: """Run metadetection on a single cell. Parameters ---------- - single_cell_coadds : `~collections.abc.Sequence` [ \ + cell_coadds : `~collections.abc.Sequence` [ \ `~lsst.cell_coadds.SingleCellCoadd` ] Per-band, per-cell coadds, in the order specified by `MetadetectionShearConfig.required_bands`. - rng: `np.random.RandomState` - Random number generator. simulate: bool, optional If set to True, simulate the data @@ -637,7 +640,10 @@ def process_cell( # columns and call 'from_pydict' instead of 'from_pylist' below, if # that's more convenient. - coadd_data = _simulate_coadd(rng) + if simulate: + coadd_data = _simulate_coadd(self.rng) + else: + coadd_data = self._cell_to_coadd_data(cell_coadds) mask_frac = _get_mask_frac( coadd_data['mfrac_mbexp'], @@ -647,12 +653,12 @@ def process_cell( mdet_config = get_mdet_config() # TODO how to get the id/tract/patch_x etc. currently in schema res = run_metadetect( - rng=rng, + rng=self.rng, config=mdet_config, **coadd_data, ) comb_res = _make_comb_data( - single_cell_coadd=single_cell_coadds[0], + cell_coadd=cell_coadds[0], res=res, meas_type=mdet_config['meas_type'], mask_frac=mask_frac, @@ -661,8 +667,8 @@ def process_cell( return comb_res - def _log_ids(self, single_cell_coadd): - idinfo = single_cell_coadd.identifiers + def _log_ids(self, cell_coadd): + idinfo = cell_coadd.identifiers tract = idinfo.tract patch_x = idinfo.patch.x patch_y = idinfo.patch.y @@ -673,6 +679,89 @@ def _log_ids(self, single_cell_coadd): ) self.log.info(mess) + # these are implemented as methods to the class to access task logging + # in _get_cell_psf + + def _cell_to_coadd_data(self, cell_coadds: Sequence[SingleCellCoadd]): + + coadd_data_list = [] + for cell_coadd in cell_coadds: + coadd_data = {} + coadd_data['coadd_exp'] = self._make_main_exposure(cell_coadd) + coadd_data['coadd_noise_exp'] = self._make_noise_exposure(cell_coadd, index=0) + coadd_data['coadd_mfrac_exp'] = self._make_mfrac_exposure(cell_coadd) + coadd_data_list.append(coadd_data) + + return extract_multiband_coadd_data(coadd_data_list) + + def _make_main_exposure(self, cell_coadd: SingleCellCoadd) -> ExposureF: + return self._make_cell_exposure( + cell_coadd.outer.image, cell_coadd, + ) + + def _make_noise_exposure(self, cell_coadd: SingleCellCoadd, index: int) -> ExposureF: + # TODO: cell coadds will have real noise realization + fake_noise_image = ImageF(cell_coadd.outer.image, True) + noise = np.median(cell_coadd.outer.variance.array[:, :]) + fake_noise_image.array[:, :] = self.rng.normal( + scale=noise, + size=fake_noise_image.array.shape, + ) + return self._make_cell_exposure( + fake_noise_image, cell_coadd, + ) + # return self._make_cell_exposure( + # cell_coadd.outer.noise_realizations[index], cell_coadd, + # ) + + def _make_mfrac_exposure(self, cell_coadd: SingleCellCoadd) -> ExposureF: + # TODO: cell coadds will have a real mfrac image + # True means deep copy + fake_mfrac_image = ImageF(cell_coadd.outer.image, True) + fake_mfrac_image.array[:, :] = 0 + return self._make_cell_exposure( + fake_mfrac_image, cell_coadd, + ) + + def _make_cell_exposure( + self, image: ImageF, cell_coadd: SingleCellCoadd, + ) -> ExposureF: + + masked_image = MaskedImageF( + image, cell_coadd.outer.mask, cell_coadd.outer.variance, + ) + + exposure = ExposureF(masked_image) + exposure.setWcs(cell_coadd.common.wcs) + + psf = self._get_cell_psf(cell_coadd) + exposure.setPsf(psf) + + band = cell_coadd.common.band + filter_label = FilterLabel(band=band, physical=band) + exposure.setFilter(filter_label) + + return exposure + + def _get_cell_psf(self, cell_coadd: SingleCellCoadd) -> KernelPsf: + + # this makes a deep copy + psf_image = cell_coadd.psf_image.convertD() + psf_array = psf_image.array + + wbad = np.where(~np.isfinite(psf_array)) + if wbad[0].size == psf_array.size: + raise ValueError('no good pixels in the psf') + + if wbad[0].size > 0: + self.log.debug('zeroing %d bad psf pixels' % wbad[0].size) + psf_array[wbad] = 0.0 + + psf_array *= 1/psf_array.sum() + # assert np.all(psf_image.array[:, :] == psf_array) + kernel = FixedKernel(psf_image) + return KernelPsf(kernel) + def _dictify(data): output = {} @@ -682,7 +771,7 @@ def _dictify(data): def _make_comb_data( - single_cell_coadd, + cell_coadd, res, meas_type, mask_frac, @@ -690,7 +779,7 @@ def _make_comb_data( ): import esutil as eu - idinfo = single_cell_coadd.identifiers + idinfo = cell_coadd.identifiers copy_dt = [ # we will copy out of arrays to these @@ -818,7 +907,6 @@ def _simulate_coadd(rng): from descwl_shear_sims.psfs import make_fixed_psf, make_ps_psf, make_rand_psf # from descwl_shear_sims.stars import make_star_catalog from descwl_coadd.coadd_nowarp import make_coadd_nowarp - from metadetect.lsst.util import extract_multiband_coadd_data g1, g2 = 0.02, 0.00 From df3e1a2c63f5c554add6b3cb2ccf2ca13f72665d Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Thu, 31 Aug 2023 09:21:48 -0700 Subject: [PATCH 09/18] a little cleanup --- python/lsst/drp/tasks/metadetection_shear.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index dfdb9d05..f7ce3c0d 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -634,12 +634,6 @@ def process_cell( from metadetect.lsst.metadetect import run_metadetect from metadetect.lsst.metadetect import get_config as get_mdet_config - # rows: list[dict[str, Any]] = [] - # TODO: run metadetection on the cell, filling in 'rows' with - # measurements. Or replace 'rows' with a 'columns' dict of numpy array - # columns and call 'from_pydict' instead of 'from_pylist' below, if - # that's more convenient. - if simulate: coadd_data = _simulate_coadd(self.rng) else: @@ -651,7 +645,7 @@ def process_cell( ) mdet_config = get_mdet_config() - # TODO how to get the id/tract/patch_x etc. currently in schema + res = run_metadetect( rng=self.rng, config=mdet_config, From a94f322e6dc0d448c00aac56c17e0fac2087934f Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Wed, 6 Sep 2023 07:43:11 -0700 Subject: [PATCH 10/18] units on T --- python/lsst/drp/tasks/metadetection_shear.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index f7ce3c0d..558b9630 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -272,7 +272,7 @@ def make_object_schema(cls, config: MetadetectionShearConfig) -> pa.Schema: nullable=False, metadata={ "doc": "admom wmom T ( + ) measurement for PSF.", - "unit": "", + "unit": "arcseconds squared", }, ), @@ -310,7 +310,7 @@ def make_object_schema(cls, config: MetadetectionShearConfig) -> pa.Schema: nullable=False, metadata={ "doc": "wmom T ( + ) measurement for reconvolved PSF.", - "unit": "", + "unit": "arcseconds squared", }, ), @@ -369,7 +369,7 @@ def make_object_schema(cls, config: MetadetectionShearConfig) -> pa.Schema: nullable=False, metadata={ "doc": "wmom T ( + ) measurement for object.", - "unit": "", + "unit": "arcseconds squared", }, ), pa.field( @@ -378,7 +378,7 @@ def make_object_schema(cls, config: MetadetectionShearConfig) -> pa.Schema: nullable=False, metadata={ "doc": "wmom T uncertainty for object.", - "unit": "", + "unit": "arcseconds squared", }, ), pa.field( From 545f21d25c3104f69cee4e01741f7d2ba52362a1 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Wed, 6 Sep 2023 08:47:55 -0700 Subject: [PATCH 11/18] field doc --- python/lsst/drp/tasks/metadetection_shear.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index 558b9630..29de5210 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -396,7 +396,7 @@ def make_object_schema(cls, config: MetadetectionShearConfig) -> pa.Schema: pa.uint32(), nullable=False, metadata={ - "doc": "wmom flux measurement flags for object.", + "doc": "wmom flux measurement flags for object in filter 1.", "unit": "", }, ), @@ -405,7 +405,7 @@ def make_object_schema(cls, config: MetadetectionShearConfig) -> pa.Schema: pa.float32(), nullable=False, metadata={ - "doc": "wmom flux for object.", + "doc": "wmom flux for object in filter 1.", "unit": "", }, ), @@ -414,7 +414,7 @@ def make_object_schema(cls, config: MetadetectionShearConfig) -> pa.Schema: pa.float32(), nullable=False, metadata={ - "doc": "wmom flux uncertainty for object.", + "doc": "wmom flux uncertainty for object in filter 1.", "unit": "", }, ), From 93b77862d191cece712c6f6dbff222895729d30c Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Wed, 13 Sep 2023 08:33:36 -0700 Subject: [PATCH 12/18] restore griz as default bands We currently run with a config restricting to r for real data, as other bands don't yet exist for the testbed Note we will want to process g separately, this is still TODO --- python/lsst/drp/tasks/metadetection_shear.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index 29de5210..85eb19fe 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -131,8 +131,8 @@ class MetadetectionShearConfig( "missing will be skipped. Bands other than those listed here will " "not be processed.", # TODO learn how to set in a config file - # default=["g", "r", "i", "z"], - default=["r"], + default=["g", "r", "i", "z"], + # default=["r"], optional=False, ) From 1dc54f984de1bce5fe1d94e295f5ff2993244d19 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Mon, 16 Oct 2023 07:23:17 -0700 Subject: [PATCH 13/18] add bright info for masking currently the sims cannot generate the bright objects because the catalog is not available Also when using real data we need to pass in the info for bright objects to be masked --- python/lsst/drp/tasks/metadetection_shear.py | 35 ++++++++++++++++++-- 1 file changed, 33 insertions(+), 2 deletions(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index 85eb19fe..05a9c0c8 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -633,11 +633,25 @@ def process_cell( """ from metadetect.lsst.metadetect import run_metadetect from metadetect.lsst.metadetect import get_config as get_mdet_config + from metadetect.lsst.masking import ( + apply_apodized_edge_masks_mbexp, + apply_apodized_bright_masks_mbexp, + ) if simulate: - coadd_data = _simulate_coadd(self.rng) + coadd_data, bright_info = _simulate_coadd(self.rng) else: coadd_data = self._cell_to_coadd_data(cell_coadds) + # TODO get bright star etc. info as input + bright_info = [] + + apply_apodized_edge_masks_mbexp(**coadd_data) + + if len(bright_info) > 0: + apply_apodized_bright_masks_mbexp( + bright_info=bright_info, + **coadd_data + ) mask_frac = _get_mask_frac( coadd_data['mfrac_mbexp'], @@ -892,6 +906,12 @@ def _get_mask_frac(mfrac_mbexp, trim_pixels=0): def _simulate_coadd(rng): + """ + TODO fake bright info + + We can't use the usual star catalog because it is not available + to the science pipelines + """ from descwl_shear_sims.sim import ( make_sim, get_sim_config, @@ -901,6 +921,8 @@ def _simulate_coadd(rng): from descwl_shear_sims.psfs import make_fixed_psf, make_ps_psf, make_rand_psf # from descwl_shear_sims.stars import make_star_catalog from descwl_coadd.coadd_nowarp import make_coadd_nowarp + from metadetect.masking import get_ap_range + from metadetect.lsst.masking import AP_RAD g1, g2 = 0.02, 0.00 @@ -979,4 +1001,13 @@ def _simulate_coadd(rng): for band in bands ] - return extract_multiband_coadd_data(coadd_data_list) + coadd_data = extract_multiband_coadd_data(coadd_data_list) + bright_info = sim_data['bright_info'] + if len(bright_info) > 0: + # Note padding due to apodization, otherwise we get donuts the + # radii coming out of the sim code are not super conservative, + # just going to the noise level + ap_padding = get_ap_range(AP_RAD) + sim_data['bright_info']['radius_pixels'] += ap_padding + + return coadd_data, bright_info From 73fb34bd9d37935b7e12be4f2eeda8d0a0a80cc0 Mon Sep 17 00:00:00 2001 From: Arun Kannawadi Date: Mon, 30 Oct 2023 11:00:09 -0700 Subject: [PATCH 14/18] Pass in a value for filterName --- python/lsst/drp/tasks/metadetection_shear.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index aeec4b3b..9cd94528 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -559,7 +559,8 @@ def runQuantum( config=self.config.ref_loader, log=self.log, ) - ref_cat = ref_loader.loadRegion(qc.quantum.dataId.region) + # What should decide the filterName? + ref_cat = ref_loader.loadRegion(qc.quantum.dataId.region, filterName="r") # THIS IS A HACK. coadds_by_band = { ref.dataId["band"]: qc.get(ref) for ref in inputRefs.input_coadds From 4b43b94b2c87a89b19c31b8d49e8673aedad6c3d Mon Sep 17 00:00:00 2001 From: Arun Kannawadi Date: Mon, 30 Oct 2023 11:00:36 -0700 Subject: [PATCH 15/18] DROP: Set a default filterMap --- python/lsst/drp/tasks/metadetection_shear.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index 9cd94528..9f9e032d 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -156,6 +156,11 @@ class MetadetectionShearConfig( # TODO: expose more configuration options here. + def setDefaults(self): + super().setDefaults() + # This is a DC2/cal_ref_cat_2_2 specific hack. This should be ideally specified in a config file + # To be removed in the cleanup before merging to main + self.ref_loader.filterMap = {band: f"lsst_{band}_smeared" for band in self.required_bands} class MetadetectionShearTask(PipelineTask): """A PipelineTask that measures shear using metadetection.""" From 70d4938a667d50fc0d0334356ae2a0c7f6c8c328 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Wed, 1 Nov 2023 06:37:37 -0700 Subject: [PATCH 16/18] set ref cat name to cal_ref_cat_2_2 for DC2 We need to make this configurable for different repos. E.g. for /repo/main we would need gaia_dr3_20230707 Removed temporary break statement added for testing Fixed flake8 complaints --- python/lsst/drp/tasks/metadetection_shear.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index 9f9e032d..8263da9f 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -65,9 +65,14 @@ class MetadetectionShearConnections(PipelineTaskConnections, dimensions={"patch" dimensions={"patch", "band"}, ) + # TODO: make "name" configurable, as it will depend on + # the repo we are using ref_cat = cT.PrerequisiteInput( doc="Reference catalog used to mask bright objects.", - name="ref_cat", + # when using /repo/main + # name="gaia_dr3_20230707", + # when using /repo/dc2 + name="cal_ref_cat_2_2", storageClass="SimpleCatalog", dimensions=("skypix",), deferLoad=True, @@ -141,9 +146,7 @@ class MetadetectionShearConfig( "Bands expected to be present. Cells with one or more of these bands " "missing will be skipped. Bands other than those listed here will " "not be processed.", - # TODO learn how to set in a config file default=["g", "r", "i", "z"], - # default=["r"], optional=False, ) @@ -158,10 +161,12 @@ class MetadetectionShearConfig( def setDefaults(self): super().setDefaults() - # This is a DC2/cal_ref_cat_2_2 specific hack. This should be ideally specified in a config file + # This is a DC2/cal_ref_cat_2_2 specific hack. This should be ideally + # specified in a config file # To be removed in the cleanup before merging to main self.ref_loader.filterMap = {band: f"lsst_{band}_smeared" for band in self.required_bands} + class MetadetectionShearTask(PipelineTask): """A PipelineTask that measures shear using metadetection.""" @@ -626,7 +631,6 @@ def run( single_cell_tables.append(table) idstart += len(res) - break # TODO: if we need to do any cell-overlap-region deduplication here # (instead of purely in science analysis code), this is where'd it'd From 13373cff789b3027697a8511d2acdcca347fafe6 Mon Sep 17 00:00:00 2001 From: Erin Sheldon Date: Wed, 1 Nov 2023 08:12:45 -0700 Subject: [PATCH 17/18] add explanatory comment from Arun K. about ref cat name --- python/lsst/drp/tasks/metadetection_shear.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/python/lsst/drp/tasks/metadetection_shear.py b/python/lsst/drp/tasks/metadetection_shear.py index 8263da9f..04c95229 100644 --- a/python/lsst/drp/tasks/metadetection_shear.py +++ b/python/lsst/drp/tasks/metadetection_shear.py @@ -65,8 +65,10 @@ class MetadetectionShearConnections(PipelineTaskConnections, dimensions={"patch" dimensions={"patch", "band"}, ) - # TODO: make "name" configurable, as it will depend on - # the repo we are using + # from Arun K.: The "name" field is configurable from the pipeline yaml + # file which are specific to what repo we are running this against. + # cal_ref_cat_2_2 is just the default value in the absence of a pipeline + # file ref_cat = cT.PrerequisiteInput( doc="Reference catalog used to mask bright objects.", # when using /repo/main From f117d77c710d6ec0abe765243dea5d2a22073247 Mon Sep 17 00:00:00 2001 From: Arun Kannawadi Date: Tue, 9 Jan 2024 13:14:00 -0500 Subject: [PATCH 18/18] Add metadetect to ups table --- ups/drp_tasks.table | 1 + 1 file changed, 1 insertion(+) diff --git a/ups/drp_tasks.table b/ups/drp_tasks.table index 2dcc831a..67b82f60 100644 --- a/ups/drp_tasks.table +++ b/ups/drp_tasks.table @@ -15,6 +15,7 @@ setupRequired(pipe_base) setupRequired(utils) setupRequired(meas_algorithms) setupRequired(pipe_tasks) +setupRequired(metadetect) # The following is boilerplate for all packages. # See https://dmtn-001.lsst.io for details on LSST_LIBRARY_PATH.