From 8ba4a2e45e25193a0182743a6200053778b584f0 Mon Sep 17 00:00:00 2001 From: Eric Charles Date: Wed, 14 Oct 2026 11:37:31 -0700 Subject: [PATCH 1/6] Add setReducedData, ShearMatch.splitAndReduce, and --joined-input CLI option MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Match.setReducedData() allows passing already-reduced DataFrames directly to a matcher without reading from disk or re-running reduceDataFrame(). ShearMatch.splitAndReduce() accepts a single joined catalog (all shear types concatenated with a type column), splits by type, reduces each subset, and stores via setReducedData — enabling an end-to-end in-memory workflow without intermediate per-type files. The shear match CLI command gains --joined-input and --shear-type-col options so the splitAndReduce path is accessible from the command line. Positional inputs and --joined-input are mutually exclusive. Generated with AI Co-Authored-By: SLAC AI --- src/hpmcm/cli/commands.py | 44 ++++++++++++++++++++++++++++++++----- src/hpmcm/cli/options.py | 21 ++++++++++++++++++ src/hpmcm/match.py | 24 ++++++++++++++++++++ src/hpmcm/shear_match.py | 38 ++++++++++++++++++++++++++++++++ tests/test_cli.py | 31 ++++++++++++++++++++++++++ tests/test_shear_match.py | 46 +++++++++++++++++++++++++++++++++++++++ tests/test_wcs_match.py | 41 ++++++++++++++++++++++++++++++++++ 7 files changed, 239 insertions(+), 6 deletions(-) diff --git a/src/hpmcm/cli/commands.py b/src/hpmcm/cli/commands.py index e8fa2f8..f86f01d 100644 --- a/src/hpmcm/cli/commands.py +++ b/src/hpmcm/cli/commands.py @@ -1,5 +1,6 @@ import click import numpy as np +import pandas import tables_io import hpmcm @@ -85,6 +86,8 @@ def shearGroup() -> None: @options.pixel_r2_cut() @options.pixel_match_scale() @options.deshear() +@options.joined_input() +@options.shear_type_col() def shearMatchCommand( inputs: list[str], output_file: click.Path, @@ -92,21 +95,50 @@ def shearMatchCommand( pixel_r2_cut: float | None, pixel_match_scale: int, deshear: bool, + joined_input: str | None, + shear_type_col: str, ) -> None: - """Match shear catalogs""" + """Match shear catalogs. + + Input catalogs can be supplied in two ways: + + \b + 1. As multiple positional INPUTS (one file per shear type, pre-split): + hpmcm shear match ns.pq 1p.pq 1m.pq 2p.pq 2m.pq --output-file out.pq + + \b + 2. As a single joined catalog via --joined-input (all types in one file): + hpmcm shear match --joined-input all_types.pq --output-file out.pq + """ + if inputs and joined_input is not None: + raise click.UsageError( + "Cannot use both positional INPUTS and --joined-input. Choose one." + ) + if not inputs and joined_input is None: + raise click.UsageError( + "Must provide either positional INPUTS or --joined-input." + ) + deshear_value: float | None = None if deshear: - assert shear is not None + if shear is None: + raise click.UsageError("--shear is required when --deshear is set.") deshear_value = -1 * shear matcher = hpmcm.ShearMatch.createShearMatch( pixelR2Cut=pixel_r2_cut, pixelMatchScale=pixel_match_scale, deshear=deshear_value, ) - input_files = list(inputs) - input_files.reverse() - visit_ids: list[int] = list(np.arange(len(input_files))) - matcher.reduceData(input_files, visit_ids) + + if joined_input is not None: + df = pandas.read_parquet(joined_input) + matcher.splitAndReduce(df, shear_type_col) + else: + input_files = list(inputs) + input_files.reverse() + visit_ids: list[int] = list(np.arange(len(input_files))) + matcher.reduceData(input_files, visit_ids) + print("Matching catalogs:") matcher.analysisLoop() print("Extracting stats.") diff --git a/src/hpmcm/cli/options.py b/src/hpmcm/cli/options.py index ada3935..b2c0c28 100644 --- a/src/hpmcm/cli/options.py +++ b/src/hpmcm/cli/options.py @@ -9,6 +9,7 @@ "catalog_type", "deshear", "inputs", + "joined_input", "output_file", "cell_size", "ra_ref", @@ -19,6 +20,7 @@ "pixel_match_scale", "pixel_r2_cut", "shear", + "shear_type_col", "snr_cut", "tract", ] @@ -169,3 +171,22 @@ def __call__(self, *args: Any, **kwargs: Any) -> Any: # pragma: no cover tract = PartialOption("--tract", help="Tract to select", type=int, required=True) inputs = PartialArgument("inputs", nargs=-1) + +joined_input = PartialOption( + "--joined-input", + type=click.Path(), + default=None, + help=( + "Path to a single joined shear catalog (all shear types concatenated). " + "Use instead of positional INPUTS to invoke splitAndReduce without " + "pre-splitting by type." + ), +) + +shear_type_col = PartialOption( + "--shear-type-col", + type=str, + default="meta_step", + show_default=True, + help="Column in --joined-input that identifies the shear catalog type.", +) diff --git a/src/hpmcm/match.py b/src/hpmcm/match.py index 80958bd..df002cc 100644 --- a/src/hpmcm/match.py +++ b/src/hpmcm/match.py @@ -164,6 +164,30 @@ def getIdOffset( cell_idx = self.getCellIdx(ix, iy) return int(self.cell_max_object * cell_idx) + def setReducedData( + self, + dataframes: list[pandas.DataFrame], + catalog_ids: list[int], + ) -> None: + """Populate red_data directly from in-memory DataFrames. + + Use instead of :meth:`reduceData` when the DataFrames are already + available in memory and need not be read from disk or passed through + :meth:`reduceDataFrame`. Each DataFrame must already have the + columns expected by the matcher (same schema as + :meth:`reduceDataFrame` output). + + Parameters + ---------- + dataframes: + One already-reduced DataFrame per catalog. + catalog_ids: + Catalog IDs corresponding to each entry in ``dataframes``. + """ + for idx, (df, cid) in enumerate(zip(dataframes, catalog_ids)): + self.catalog_id_map[cid] = idx + self.red_data[cid] = df + def reduceData( self, input_files: list[str], diff --git a/src/hpmcm/shear_match.py b/src/hpmcm/shear_match.py index 1136038..f6026fd 100644 --- a/src/hpmcm/shear_match.py +++ b/src/hpmcm/shear_match.py @@ -270,6 +270,44 @@ def _buildCellData( ) -> CellData: return ShearCellData(self, id_offset, corner, size, idx, self.cell_buffer) + def splitAndReduce( + self, + joined_df: pandas.DataFrame, + shear_type_col: str = "meta_step", + ) -> None: + """Split a joined shear catalog by type and populate red_data. + + Splits ``joined_df`` by shear catalog type, runs each subset through + :meth:`reduceDataFrame`, and stores the results via + :meth:`~hpmcm.match.Match.setReducedData`. No files are read or + written. + + The joined DataFrame must already contain all columns expected by + :meth:`reduceDataFrame`: ``x_cell_coadd``, ``y_cell_coadd``, + ``cell_idx_x``, ``cell_idx_y``, ``x_pix``, ``y_pix``, ``snr``, + ``g_1``, ``g_2``, ``id``. This is the schema produced by + :func:`~hpmcm.shear_utils.splitDESCMDByTypeAndClean` and + :func:`~hpmcm.shear_utils.splitRubinMDByTypeAndClean`. + + Parameters + ---------- + joined_df: + Single DataFrame with all shear types concatenated; must have a + column whose values are drawn from + :py:data:`~hpmcm.shear_utils.SHEAR_NAMES`. + shear_type_col: + Name of the column that identifies the shear catalog type. + Defaults to ``'meta_step'`` (the column written by the split + functions in :mod:`hpmcm.shear_utils`). + """ + dfs = [] + cids = [] + for cid, name in enumerate(self.shear_names): + sub = joined_df[joined_df[shear_type_col] == name] + dfs.append(self.reduceDataFrame(sub)) + cids.append(cid) + self.setReducedData(dfs, cids) + def extractShearStats(self, central_only: bool = True) -> dict[str, pandas.DataFrame]: """Extract shear stats diff --git a/tests/test_cli.py b/tests/test_cli.py index 3ae57a0..04aa9ac 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -32,3 +32,34 @@ def testCliShearMatch() -> None: result = runner.invoke(cli, "shear match --help") checkResult(result) + assert "--joined-input" in result.output + assert "--shear-type-col" in result.output + + +def testCliShearMatchMutualExclusion() -> None: + """Both positional inputs and --joined-input should be rejected.""" + runner = CliRunner() + + result = runner.invoke( + cli, + [ + "shear", "match", + "file1.pq", + "--joined-input", "joined.pq", + "--output-file", "out.pq", + ], + ) + assert result.exit_code != 0 + assert "Cannot use both" in result.output + + +def testCliShearMatchNoInput() -> None: + """No inputs at all should be rejected.""" + runner = CliRunner() + + result = runner.invoke( + cli, + ["shear", "match", "--output-file", "out.pq"], + ) + assert result.exit_code != 0 + assert "Must provide" in result.output diff --git a/tests/test_shear_match.py b/tests/test_shear_match.py index 069064e..7a48cac 100644 --- a/tests/test_shear_match.py +++ b/tests/test_shear_match.py @@ -2,6 +2,7 @@ import os import numpy as np +import pandas import hpmcm @@ -132,3 +133,48 @@ def testShearReport(setup_data: int) -> None: tract=TRACT, snr_cut=SNR_CUT, ) + + +def testShearMatchSplitAndReduce(setup_data: int) -> None: + """Verify splitAndReduce on a joined catalog produces valid results.""" + assert setup_data == 0 + + # Build a joined catalog by reading the per-type files and tagging each row. + per_type_files = sorted( + glob.glob( + os.path.join(DATADIR, f"shear_{CATALOG_TYPE}_{SHEAR_ST}_uncleaned_{TRACT}_*.pq") + ) + ) + dfs = [] + for f in per_type_files: + # Infer type name from file suffix (e.g., "_ns.pq" → "ns"). + type_ = os.path.splitext(f)[0].rsplit("_", 1)[-1] + df = pandas.read_parquet(f) + df = df.copy() + df["meta_step"] = type_ + dfs.append(df) + joined_df = pandas.concat(dfs, ignore_index=True) + + matcher = hpmcm.ShearMatch.createShearMatch( + geometry=OUR_GEOM, + pixel_R2_cut=PIXEL_R2CUT, + pixel_match_scale=PIXEL_MATCH_SCALE, + deshear=-1 * SHEAR, + ) + matcher.splitAndReduce(joined_df) + + x_range = range(50, 55) + y_range = range(170, 175) + matcher.analysisLoop(x_range, y_range) + + stats = matcher.extractStats() + assert stats is not None + + obj_lists = hpmcm.classify.classifyObjects(matcher, SNRCut=SNR_CUT) + n_good = len(obj_lists["ideal"]) + bad_list = [ + "edge_mixed", "edge_missing", "edge_extra", "orphan", + "missing", "two_missing", "many_missing", "extra", "caught", + ] + n_bad = np.sum([len(obj_lists[x]) for x in bad_list]) + assert n_good / (n_good + n_bad) > 0.90 diff --git a/tests/test_wcs_match.py b/tests/test_wcs_match.py index d62aed4..e0a6c20 100644 --- a/tests/test_wcs_match.py +++ b/tests/test_wcs_match.py @@ -144,3 +144,44 @@ def testWCSMatchPixelMatchScale(setup_data: int) -> None: n_bad = np.sum([len(obj_lists[x]) for x in bad_list]) effic = n_good / (n_good + n_bad) assert effic > 0.75 + + +def testWCSMatchSetReducedData(setup_data: int) -> None: + """Verify setReducedData produces the same result as reduceData.""" + assert setup_data == 0 + + source_tablesfiles = sorted( + glob.glob(os.path.join(DATADIR, f"shear_*_{SHEAR_ST}_cleaned_{TRACT}_ns.pq")) + ) + source_tablesfiles.append(os.path.join(DATADIR, f"object_{TRACT}.pq")) + source_tablesfiles.reverse() + source_tablesfiles = [source_tablesfiles[0], source_tablesfiles[1]] + catalog_ids: list[int] = list(np.arange(len(source_tablesfiles))) + + # Reduce via the normal file path to capture the reduced DataFrames. + ref_matcher = hpmcm.WcsMatch.create( + REF_DIR, REGION_SIZE, pixel_size=PIXEL_SIZE, pixel_R2_cut=PIXEL_R2CUT + ) + ref_matcher.reduceData(source_tablesfiles, catalog_ids) + + # Build a fresh matcher and supply the already-reduced frames directly. + matcher = hpmcm.WcsMatch.create( + REF_DIR, REGION_SIZE, pixel_size=PIXEL_SIZE, pixel_R2_cut=PIXEL_R2CUT + ) + matcher.setReducedData( + list(ref_matcher.red_data.values()), + list(ref_matcher.red_data.keys()), + ) + + matcher.analysisLoop(range(1, 2), range(1, 2)) + stats = matcher.extractStats() + assert stats is not None + + obj_lists = hpmcm.classify.classifyObjects(matcher, SNRCut=10.0) + n_good = len(obj_lists["ideal"]) + bad_list = [ + "edge_mixed", "edge_missing", "edge_extra", "orphan", + "missing", "two_missing", "many_missing", "extra", "caught", + ] + n_bad = np.sum([len(obj_lists[x]) for x in bad_list]) + assert n_good / (n_good + n_bad) > 0.85 From 3507023a1f5d3675bad4bf1d2e27688ef068f7f1 Mon Sep 17 00:00:00 2001 From: Eric Charles Date: Wed, 14 Oct 2026 16:16:47 -0700 Subject: [PATCH 2/6] Add catalog-format support to shear match CLI and prepare* functions Add prepareRubinMDForMatch and prepareDESCMDForMatch to shear_utils: these are in-memory counterparts to splitRubinMDByTypeAndClean and splitDESCMDByTypeAndClean that return a joined DataFrame (with meta_step column) instead of writing per-type files. Add --catalog-format option (pre-processed / rubin-md / desc-md) to the shear match CLI command. For rubin-md the tract must also be supplied via --tract. The format drives which preparation function is called before splitAndReduce; pre-processed (the default) passes the file straight to splitAndReduce, relying on --shear-type-col to identify the type column. Generated with AI Co-Authored-By: SLAC AI --- src/hpmcm/cli/commands.py | 36 +++++++- src/hpmcm/cli/options.py | 15 ++++ src/hpmcm/shear_utils.py | 184 ++++++++++++++++++++++++++++++++++++++ tests/test_cli.py | 19 ++++ 4 files changed, 251 insertions(+), 3 deletions(-) diff --git a/src/hpmcm/cli/commands.py b/src/hpmcm/cli/commands.py index f86f01d..7115cf1 100644 --- a/src/hpmcm/cli/commands.py +++ b/src/hpmcm/cli/commands.py @@ -87,7 +87,9 @@ def shearGroup() -> None: @options.pixel_match_scale() @options.deshear() @options.joined_input() +@options.catalog_format() @options.shear_type_col() +@options.tract(required=False, default=None) def shearMatchCommand( inputs: list[str], output_file: click.Path, @@ -96,7 +98,9 @@ def shearMatchCommand( pixel_match_scale: int, deshear: bool, joined_input: str | None, + catalog_format: str, shear_type_col: str, + tract: int | None, ) -> None: """Match shear catalogs. @@ -107,8 +111,22 @@ def shearMatchCommand( hpmcm shear match ns.pq 1p.pq 1m.pq 2p.pq 2m.pq --output-file out.pq \b - 2. As a single joined catalog via --joined-input (all types in one file): - hpmcm shear match --joined-input all_types.pq --output-file out.pq + 2. As a single joined catalog via --joined-input (all types in one file). + Use --catalog-format to specify the input format: + + \b + # already-processed catalog (x_cell_coadd etc. present): + hpmcm shear match --joined-input all.pq --output-file out.pq + + \b + # raw joined Rubin MD catalog (requires --tract): + hpmcm shear match --joined-input raw.pq --catalog-format rubin-md \\ + --tract 10463 --output-file out.pq + + \b + # raw joined DESC MD catalog: + hpmcm shear match --joined-input raw.pq --catalog-format desc-md \\ + --output-file out.pq """ if inputs and joined_input is not None: raise click.UsageError( @@ -118,6 +136,10 @@ def shearMatchCommand( raise click.UsageError( "Must provide either positional INPUTS or --joined-input." ) + if catalog_format == "rubin-md" and tract is None: + raise click.UsageError( + "--tract is required when --catalog-format is 'rubin-md'." + ) deshear_value: float | None = None if deshear: @@ -132,7 +154,15 @@ def shearMatchCommand( if joined_input is not None: df = pandas.read_parquet(joined_input) - matcher.splitAndReduce(df, shear_type_col) + if catalog_format == "rubin-md": + assert tract is not None # validated above + df = hpmcm.shear_utils.prepareRubinMDForMatch(df, tract) + matcher.splitAndReduce(df) + elif catalog_format == "desc-md": + df = hpmcm.shear_utils.prepareDESCMDForMatch(df) + matcher.splitAndReduce(df) + else: + matcher.splitAndReduce(df, shear_type_col) else: input_files = list(inputs) input_files.reverse() diff --git a/src/hpmcm/cli/options.py b/src/hpmcm/cli/options.py index b2c0c28..fd4e869 100644 --- a/src/hpmcm/cli/options.py +++ b/src/hpmcm/cli/options.py @@ -6,6 +6,7 @@ __all__: list[str] = [ "basefile", + "catalog_format", "catalog_type", "deshear", "inputs", @@ -172,6 +173,20 @@ def __call__(self, *args: Any, **kwargs: Any) -> Any: # pragma: no cover inputs = PartialArgument("inputs", nargs=-1) +catalog_format = PartialOption( + "--catalog-format", + type=click.Choice(["pre-processed", "rubin-md", "desc-md"]), + default="pre-processed", + show_default=True, + help=( + "Format of the --joined-input catalog. " + "'pre-processed': coordinates already computed; expects x_cell_coadd, " + "cell_idx_x/y, and --shear-type-col. " + "'rubin-md': raw joined Rubin MD catalog (requires --tract). " + "'desc-md': raw joined DESC MD catalog." + ), +) + joined_input = PartialOption( "--joined-input", type=click.Path(), diff --git a/src/hpmcm/shear_utils.py b/src/hpmcm/shear_utils.py index e260085..16a9045 100644 --- a/src/hpmcm/shear_utils.py +++ b/src/hpmcm/shear_utils.py @@ -281,6 +281,190 @@ def mergeShearReports( out_df.to_parquet(output_file) +def prepareRubinMDForMatch( + df: pandas.DataFrame, + tract: int, + *, + clean: bool = False, + geometry: ShearCellGeometry = DEFAULT_GEOMETRY, + shear_names: list[str] | None = None, +) -> pandas.DataFrame: # pragma: no cover + """Process a raw joined Rubin MD DataFrame for in-memory matching. + + Applies the same spatial filtering and coordinate computation as + :func:`splitRubinMDByTypeAndClean` but returns a single joined + DataFrame instead of writing per-type files. Pass the result + directly to :meth:`~hpmcm.shear_match.ShearMatch.splitAndReduce`. + + Parameters + ---------- + df: + Raw joined Rubin MD catalog with all shear types concatenated. + Required columns: ``metaStep``, ``cell_x``, ``cell_y``, + ``patch`` (or ``patch_x`` / ``patch_y``), ``x``, ``y``, + ``tract``, ``shearObjectId``, ``snr``, ``g_1``, ``g_2``, + ``ra``, ``dec``. + tract: + Tract to select. + clean: + When ``True``, apply strict cell bounds (no ``cell_inner_buffer`` + padding around the inner cell). + geometry: + Cell geometry configuration. + shear_names: + Subset of :py:data:`SHEAR_NAMES` to process. ``None`` → all five. + + Returns + ------- + Joined DataFrame with a ``meta_step`` column and the columns expected + by :meth:`~hpmcm.shear_match.ShearMatch.reduceDataFrame`: + ``id``, ``x_cell_coadd``, ``y_cell_coadd``, ``cell_idx_x``, + ``cell_idx_y``, ``x_pix``, ``y_pix``, ``snr``, ``g_1``, ``g_2``, + ``ra``, ``dec``. + """ + if shear_names is None: + shear_names = SHEAR_NAMES + buf = 0 if clean else geometry.cell_inner_buffer + + frames: list[pandas.DataFrame] = [] + for type_ in shear_names: + sub = df[df["metaStep"] == type_] + + right_tract = sub["tract"] == tract + central_to_patch = ( + np.fabs(sub["cell_x"].values - geometry.patch_offset) + < (geometry.n_cell_in_patch / 2) + ) & ( + np.fabs(sub["cell_y"].values - geometry.patch_offset) + < (geometry.n_cell_in_patch / 2) + ) + sub = sub[right_tract & central_to_patch].copy(deep=True) + + if "patch_x" not in sub.columns: + sub["patch_x"] = sub["patch"] % 10 + sub["patch_y"] = sub["patch"] // 10 + + cell_idx_x = ( + geometry.n_cell_in_patch * sub["patch_x"].values + sub["cell_x"].values + ).astype(int) + cell_idx_y = ( + geometry.n_cell_in_patch * sub["patch_y"].values + sub["cell_y"].values + ).astype(int) + + x_cell_coadd = ( + sub["x"].values + - (cell_idx_x - geometry.n_cell_patch_buffer) * geometry.cell_inner_size + ) + y_cell_coadd = ( + sub["y"].values + - (cell_idx_y - geometry.n_cell_patch_buffer) * geometry.cell_inner_size + ) + + central_to_cell = ( + (x_cell_coadd >= -buf) + & (x_cell_coadd < geometry.cell_inner_size + buf) + & (y_cell_coadd >= -buf) + & (y_cell_coadd < geometry.cell_inner_size + buf) + ) + out = sub[central_to_cell].copy(deep=True) + out["x_cell_coadd"] = x_cell_coadd[central_to_cell] + out["y_cell_coadd"] = y_cell_coadd[central_to_cell] + out["x_pix"] = out["x"] + out["y_pix"] = out["y"] + out["cell_idx_x"] = cell_idx_x[central_to_cell] + out["cell_idx_y"] = cell_idx_y[central_to_cell] + out["id"] = out["shearObjectId"] + out["meta_step"] = type_ + frames.append(out) + + return pandas.concat(frames, ignore_index=True) + + +def prepareDESCMDForMatch( + df: pandas.DataFrame, + *, + clean: bool = False, + geometry: ShearCellGeometry = DEFAULT_GEOMETRY, + shear_names: list[str] | None = None, +) -> pandas.DataFrame: # pragma: no cover + """Process a raw joined DESC MD DataFrame for in-memory matching. + + Applies the same spatial filtering and coordinate computation as + :func:`splitDESCMDByTypeAndClean` but returns a single joined + DataFrame instead of writing per-type files. Pass the result + directly to :meth:`~hpmcm.shear_match.ShearMatch.splitAndReduce`. + + Parameters + ---------- + df: + Raw joined DESC MD catalog with all shear types concatenated. + Required columns: ``mcal_step``, ``cell_i``, ``cell_j``, + ``patch_x``, ``patch_y``, ``xcell``, ``ycell``, ``x``, ``y``, + ``snr``, ``g_1``, ``g_2``, ``ra``, ``dec``. If ``id`` is + absent it is generated as ``arange(len(df))``. + clean: + When ``True``, apply strict cell bounds (no ``cell_inner_buffer`` + padding around the inner cell). + geometry: + Cell geometry configuration. + shear_names: + Subset of :py:data:`SHEAR_NAMES` to process. ``None`` → all five. + + Returns + ------- + Joined DataFrame with a ``meta_step`` column and the columns expected + by :meth:`~hpmcm.shear_match.ShearMatch.reduceDataFrame`. + """ + if shear_names is None: + shear_names = SHEAR_NAMES + buf = 0 if clean else geometry.cell_inner_buffer + + if "id" not in df.columns: + df = df.copy() + df["id"] = np.arange(len(df)) + + frames: list[pandas.DataFrame] = [] + for type_ in shear_names: + sub = df[df["mcal_step"] == type_] + + central_to_patch = ( + np.fabs(sub["cell_i"].values - geometry.patch_offset) + < (geometry.n_cell_in_patch / 2) + ) & ( + np.fabs(sub["cell_j"].values - geometry.patch_offset) + < (geometry.n_cell_in_patch / 2) + ) + sub = sub[central_to_patch].copy(deep=True) + + cell_idx_x = ( + geometry.n_cell_in_patch * sub["patch_x"].values + sub["cell_j"].values + ).astype(int) + cell_idx_y = ( + geometry.n_cell_in_patch * sub["patch_y"].values + sub["cell_i"].values + ).astype(int) + + x_cell_coadd = sub["xcell"].values - geometry.cell_buffer + y_cell_coadd = sub["ycell"].values - geometry.cell_buffer + + central_to_cell = ( + (x_cell_coadd >= -buf) + & (x_cell_coadd < geometry.cell_inner_size + buf) + & (y_cell_coadd >= -buf) + & (y_cell_coadd < geometry.cell_inner_size + buf) + ) + out = sub[central_to_cell].copy(deep=True) + out["x_cell_coadd"] = x_cell_coadd[central_to_cell] + out["y_cell_coadd"] = y_cell_coadd[central_to_cell] + out["x_pix"] = out["x"] + out["y_pix"] = out["y"] + out["cell_idx_x"] = cell_idx_x[central_to_cell] + out["cell_idx_y"] = cell_idx_y[central_to_cell] + out["meta_step"] = type_ + frames.append(out) + + return pandas.concat(frames, ignore_index=True) + + def splitRubinMDByTypeAndClean( basefile: str, tract: int, diff --git a/tests/test_cli.py b/tests/test_cli.py index 04aa9ac..f2b5cad 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -33,7 +33,9 @@ def testCliShearMatch() -> None: result = runner.invoke(cli, "shear match --help") checkResult(result) assert "--joined-input" in result.output + assert "--catalog-format" in result.output assert "--shear-type-col" in result.output + assert "--tract" in result.output def testCliShearMatchMutualExclusion() -> None: @@ -63,3 +65,20 @@ def testCliShearMatchNoInput() -> None: ) assert result.exit_code != 0 assert "Must provide" in result.output + + +def testCliShearMatchRubinMDRequiresTract() -> None: + """--catalog-format rubin-md without --tract should be rejected.""" + runner = CliRunner() + + result = runner.invoke( + cli, + [ + "shear", "match", + "--joined-input", "raw.pq", + "--catalog-format", "rubin-md", + "--output-file", "out.pq", + ], + ) + assert result.exit_code != 0 + assert "--tract" in result.output From e8369d7ad2d819a0be2519f67242317284a24f8b Mon Sep 17 00:00:00 2001 From: Eric Charles Date: Wed, 14 Oct 2026 16:43:18 -0700 Subject: [PATCH 3/6] Update docs for in-memory catalog loading and catalog-format CLI option - docs/source/cli.rst: add narrative for shear match --joined-input modes, including a table of --catalog-format values (pre-processed/rubin-md/desc-md), required extra options, and three example command lines; fix pre-existing RST title-overline length warnings on the split-rubin/split-desc sections - src/hpmcm/__init__.py: add two new bullet points to the Additional features section describing setReducedData/splitAndReduce and the prepare* functions Generated with AI Co-Authored-By: SLAC AI --- docs/source/cli.rst | 80 ++++++++++++++++++++++++++++++++++++++++--- src/hpmcm/__init__.py | 16 +++++++++ 2 files changed, 91 insertions(+), 5 deletions(-) diff --git a/docs/source/cli.rst b/docs/source/cli.rst index 293aefa..c0bd4fc 100644 --- a/docs/source/cli.rst +++ b/docs/source/cli.rst @@ -19,19 +19,89 @@ Matching using cell-based coadd frame :prog: hpmcm shear match :nested: full - -============================================= +Input catalog modes +------------------- + +The ``shear match`` command accepts catalogs in two ways: + +**Multiple positional inputs** — one pre-split parquet file per shear type, +passed as positional arguments. Files should be in reverse alphabetical order +(``ns`` last) to match the expected catalog-ID ordering:: + + hpmcm shear match ns.pq 1p.pq 1m.pq 2p.pq 2m.pq --output-file out.pq + +This is the classic workflow: split the raw catalog with ``split-rubin`` or +``split-desc`` first, then pass the resulting per-type files here. + +**Single joined catalog** — pass ``--joined-input`` with one file that contains +all shear types. Use ``--catalog-format`` to tell the command what format the +file is in: + +.. list-table:: + :header-rows: 1 + :widths: 20 45 35 + + * - ``--catalog-format`` + - Input expected + - Extra options + * - ``pre-processed`` *(default)* + - Coordinates already computed (``x_cell_coadd``, ``cell_idx_x/y`` + present); shear type identified by ``--shear-type-col`` (default + ``meta_step``). + - ``--shear-type-col`` if the column is named differently + * - ``rubin-md`` + - Raw joined Rubin MD catalog (``metaStep``, ``cell_x/y``, ``patch``, + ``x/y``, ``shearObjectId``, …). Tract/patch filtering and coordinate + computation are applied automatically. + - ``--tract`` *(required)* + * - ``desc-md`` + - Raw joined DESC MD catalog (``mcal_step``, ``cell_i/j``, + ``xcell/ycell``, ``x/y``, …). Patch filtering and coordinate + computation are applied automatically. + - — + +Example — raw Rubin MD catalog:: + + hpmcm shear match \ + --joined-input raw_rubin.parq \ + --catalog-format rubin-md \ + --tract 10463 \ + --deshear --shear 0.01 \ + --output-file out.pq + +Example — raw DESC MD catalog:: + + hpmcm shear match \ + --joined-input raw_desc.parq \ + --catalog-format desc-md \ + --deshear --shear 0.01 \ + --output-file out.pq + +Example — pre-processed joined catalog (e.g.\ files concatenated after +running ``split-rubin``, where the type column is ``meta_step``):: + + hpmcm shear match \ + --joined-input joined.parq \ + --output-file out.pq + +When ``--joined-input`` is used, the command internally calls +:func:`~hpmcm.shear_utils.prepareRubinMDForMatch` or +:func:`~hpmcm.shear_utils.prepareDESCMDForMatch` (for raw formats) and then +:meth:`~hpmcm.shear_match.ShearMatch.splitAndReduce` to populate the matcher +in memory without writing intermediate files. + +================================================ Splitting Rubin input catalogs for shear matching -============================================= +================================================ .. click:: hpmcm.cli.commands:shearSplitRubinCommand :prog: hpmcm shear split-rubin :nested: full -============================================= +================================================ Splitting DESC input catalogs for shear matching -============================================= +================================================ .. click:: hpmcm.cli.commands:shearSplitDESCCommand :prog: hpmcm shear split-desc diff --git a/src/hpmcm/__init__.py b/src/hpmcm/__init__.py index b56bd33..034bf13 100644 --- a/src/hpmcm/__init__.py +++ b/src/hpmcm/__init__.py @@ -56,6 +56,22 @@ :py:func:`buildJoinedClusterTable ` combine the association and statistics tables into a single flat table for downstream analysis. + +* **In-memory catalog loading** — :py:meth:`Match.setReducedData() + ` lets you pass already-reduced + DataFrames directly to any matcher, bypassing file I/O. + :py:meth:`ShearMatch.splitAndReduce() + ` accepts a single joined + DataFrame (all shear types concatenated with a type column), splits it by + type, and loads the result — no intermediate per-type files needed. + +* **Raw catalog preparation** — :py:func:`shear_utils.prepareRubinMDForMatch + ` and + :py:func:`shear_utils.prepareDESCMDForMatch + ` apply the same spatial + filtering and coordinate computation as the ``split-rubin`` / ``split-desc`` + CLI commands but return an in-memory joined DataFrame instead of writing + files. Pass the result directly to ``splitAndReduce``. """ from . import classify, match_utils, package_utils, shear_utils, utils, viz_utils From ca6e2488efd90036d86b31c861a52591b81ae43f Mon Sep 17 00:00:00 2001 From: Eric Charles Date: Wed, 14 Oct 2026 16:58:06 -0700 Subject: [PATCH 4/6] Add run_matching_pipeline.py combining StandarizeInputs + RunMatching notebooks Single self-contained script with a configuration block at the top for tract, file paths, shear parameters, and step-control flags. Steps: 1. Standardize Rubin MD / DESC MD / object / anacal catalogs 2. Rubin shear matching + joined object table 3. DESC shear matching + joined object table 4. WCS meta-matching across all catalogs Generated with AI Co-Authored-By: SLAC AI --- examples/rubin_dp2/run_matching_pipeline.py | 295 ++++++++++++++++++++ 1 file changed, 295 insertions(+) create mode 100644 examples/rubin_dp2/run_matching_pipeline.py diff --git a/examples/rubin_dp2/run_matching_pipeline.py b/examples/rubin_dp2/run_matching_pipeline.py new file mode 100644 index 0000000..40d72c0 --- /dev/null +++ b/examples/rubin_dp2/run_matching_pipeline.py @@ -0,0 +1,295 @@ +#!/usr/bin/env python +""" +HPMCM full matching pipeline for a single Rubin tract. + +Combines the actions of StandarizeInputs.ipynb and RunMatching.ipynb: + 1. Standardize inputs — split raw catalogs by shear type and reduce columns + 2. Rubin shear matching + 3. DESC shear matching + 4. WCS meta-matching across all object catalogs + +Edit the CONFIGURATION block below for your tract, then run: + + python run_matching_pipeline.py +""" +from __future__ import annotations + +import tables_io + +import hpmcm +from hpmcm import shear_utils, utils + +# ============================================================ +# CONFIGURATION — edit this block before running +# ============================================================ + +TRACT: int = 8006 + +# Tract centre (RA, Dec in degrees); used to build the shear-match WCS. +REF_DIR: tuple[float, float] = (340.5, -9.6694214876) + +# ── Raw input files ─────────────────────────────────────────────────────────── +# Full path to the raw joined Rubin metadetect catalog (all shear types in one +# file, produced before split-by-type). +RUBIN_MD_RAW_FILE: str = f"{TRACT}/shear_object_{TRACT}.parq" + +# Full path to the raw joined DESC metadetect catalog. +DESC_MD_RAW_FILE: str = f"{TRACT}/mdet-run-{TRACT}.parq" + +# Full path to the Rubin object catalog. +OBJECT_RAW_FILE: str = f"{TRACT}/object_{TRACT}.parq" + +# Full path to the anacal catalog. +ANACAL_RAW_FILE: str = f"{TRACT}/anacal_merged_{TRACT}.parq" + +# ── Shear parameters ────────────────────────────────────────────────────────── +# Magnitude of the applied shear (|g|). Used to set the matcher's deshear +# parameter (-RUBIN_SHEAR / -DESC_SHEAR) so matching is done on de-sheared +# coordinates. +RUBIN_SHEAR: float = 0.01 +DESC_SHEAR: float = 0.01 + +# Shear catalog types to include in each matching run. +# Rubin MD produces five catalogs; DESC MD produces three. +RUBIN_SHEAR_NAMES: list[str] = ["ns", "1m", "1p", "2m", "2p"] +DESC_SHEAR_NAMES: list[str] = ["ns", "1m", "1p"] + +# ── Shear matching parameters ───────────────────────────────────────────────── +SHEAR_PIXEL_MATCH_SCALE: int = 2 +SHEAR_PIXEL_R2_CUT: float = 16.0 + +# ── WCS meta-matching parameters ────────────────────────────────────────────── +# Pixel scale in degrees per pixel (0.2 arcsec/pix = 0.2/3600 deg/pix). +WCS_PIXEL_SCALE: float = 0.2 / 3600.0 +# Side length of the matching region in pixels. +WCS_N_PIX_REGION: int = 30_000 +WCS_CELL_SIZE: int = 150 +WCS_CELL_BUFFER: int = 5 +WCS_CELL_MAX_OBJECT: int = 10_000 +WCS_PIXEL_MATCH_SCALE: int = 4 +WCS_PIXEL_R2_CUT: float = 16.0 + +# ── Step control ────────────────────────────────────────────────────────────── +# Set any flag to False to skip that step (e.g. if the output files already +# exist from a previous run). +REDO_RUBIN_MD: bool = True # split + reduce raw Rubin MD catalog +REDO_DESC_MD: bool = True # split + reduce raw DESC MD catalog +REDO_OBJECT: bool = True # reduce Rubin object catalog +REDO_ANACAL: bool = True # reduce anacal catalog + +REMATCH_RUBIN: bool = True # run Rubin shear matching and save +REMATCH_DESC: bool = True # run DESC shear matching and save +REMATCH_OBJECTS: bool = True # run WCS meta-matching and save + +BUILD_RUBIN_JOINED: bool = True # build joined Rubin object table +BUILD_DESC_JOINED: bool = True # build joined DESC object table + +# Set True to run the analysis on a single cell only (quick smoke-test). +TEST_ONLY: bool = False + +# ============================================================ +# END CONFIGURATION — no edits needed below this line +# ============================================================ + +# ── Derived intermediate file paths ─────────────────────────────────────────── +# These follow the naming conventions of splitRubinMDByTypeAndClean / +# splitDESCMDByTypeAndClean and the reduce* helpers. + +def _rubin_split_file(st: str) -> str: + """Per-type Rubin file written by splitRubinMDByTypeAndClean.""" + return RUBIN_MD_RAW_FILE.replace(".parq", f"_uncleaned_{TRACT}_{st}.parq") + +def _desc_split_file(st: str) -> str: + """Per-type DESC file written by splitDESCMDByTypeAndClean.""" + return DESC_MD_RAW_FILE.replace(".parq", f"_uncleaned_{TRACT}_{st}.parq") + +def _rubin_reduced_file(st: str) -> str: + """Column-reduced Rubin file written by reduceRubinMDTable.""" + return f"{TRACT}/shear_object_reduced_uncleaned_{TRACT}_{st}.parq" + +def _desc_reduced_file(st: str) -> str: + """Column-reduced DESC file written by reduceDESCMDTable.""" + return f"{TRACT}/mdet-run-reduced_uncleaned_{TRACT}_{st}.parq" + +RUBIN_SHEAR_FILES: list[str] = [_rubin_reduced_file(st) for st in RUBIN_SHEAR_NAMES] +RUBIN_SHEAR_IDS: list[int] = list(range(len(RUBIN_SHEAR_FILES))) + +DESC_SHEAR_FILES: list[str] = [_desc_reduced_file(st) for st in DESC_SHEAR_NAMES] +DESC_SHEAR_IDS: list[int] = list(range(len(DESC_SHEAR_FILES))) + +# Meta-matching input files: object catalog, anacal, and the two shear-match +# joined object tables (produced by this script in steps 2 and 3). +ALL_FILES: list[str] = [ + f"{TRACT}/object_reduced_{TRACT}.parq", + f"{TRACT}/anacal_reduced_{TRACT}.parq", + f"{TRACT}/desc_md_match_object_joined.pq", + f"{TRACT}/rubin_md_match_object_joined.pq", +] +ALL_IDS: list[int] = list(range(len(ALL_FILES))) + +OUR_GEOM = shear_utils.ShearCellGeometry(ref_dir=REF_DIR) +_WCS_REGION_DEG = WCS_N_PIX_REGION * WCS_PIXEL_SCALE + +_TEST_RANGE = range(100, 101) # single cell used when TEST_ONLY=True + +# ── Step 1a: Standardize Rubin MD catalog ───────────────────────────────────── + +if REDO_RUBIN_MD: + print("Standardizing Rubin MD catalog...") + shear_utils.splitRubinMDByTypeAndClean( + RUBIN_MD_RAW_FILE, TRACT, RUBIN_SHEAR, deshear=False + ) + for st in RUBIN_SHEAR_NAMES: + utils.reduceRubinMDTable(_rubin_split_file(st), _rubin_reduced_file(st)) + print(" Done.") + +# ── Step 1b: Standardize DESC MD catalog ────────────────────────────────────── + +if REDO_DESC_MD: + print("Standardizing DESC MD catalog...") + shear_utils.splitDESCMDByTypeAndClean( + DESC_MD_RAW_FILE, TRACT, DESC_SHEAR, deshear=False + ) + for st in DESC_SHEAR_NAMES: + utils.reduceDESCMDTable(_desc_split_file(st), _desc_reduced_file(st)) + print(" Done.") + +# ── Step 1c: Reduce object catalog ──────────────────────────────────────────── + +if REDO_OBJECT: + print("Standardizing object catalog...") + utils.reduceObjectTable(OBJECT_RAW_FILE, f"{TRACT}/object_reduced_{TRACT}.parq") + print(" Done.") + +# ── Step 1d: Reduce anacal catalog ──────────────────────────────────────────── + +if REDO_ANACAL: + print("Standardizing anacal catalog...") + utils.reduceAnacalTable(ANACAL_RAW_FILE, f"{TRACT}/anacal_reduced_{TRACT}.parq") + print(" Done.") + +# ── Step 2: Rubin shear matching ────────────────────────────────────────────── + +if REMATCH_RUBIN: + print("Running Rubin shear matching...") + rubin_shear_match = hpmcm.ShearMatch.createShearMatch( + geometry=OUR_GEOM, + pixel_match_scale=SHEAR_PIXEL_MATCH_SCALE, + pixel_r2_cut=SHEAR_PIXEL_R2_CUT, + deshear=-RUBIN_SHEAR, + shear_names=RUBIN_SHEAR_NAMES, + ) + rubin_shear_match.reduceData(RUBIN_SHEAR_FILES, RUBIN_SHEAR_IDS) + if TEST_ONLY: + rubin_shear_match.analysisLoop(_TEST_RANGE, _TEST_RANGE) + else: + rubin_shear_match.analysisLoop() + rubin_stats = rubin_shear_match.extractStats() + rubin_stats.update(**rubin_shear_match.extractShearStats()) + tables_io.write(rubin_stats, f"{TRACT}/rubin_md_match_") + rubin_shear_match.save(f"{TRACT}/rubin_md_match") + print(" Done.") + +# ── Step 2 (cont.): Build joined Rubin object table ─────────────────────────── + +if BUILD_RUBIN_JOINED: + print("Building joined Rubin object table...") + rubin_obj_stats = tables_io.read(f"{TRACT}/rubin_md_match_object_stats.pq") + rubin_obj_assoc = tables_io.read(f"{TRACT}/rubin_md_match_object_assoc.pq") + joined_rubin = hpmcm.output_tables.buildJoinedObjectTable( + rubin_obj_stats, + rubin_obj_assoc, + RUBIN_SHEAR_FILES, + RUBIN_SHEAR_IDS, + object_shear=f"{TRACT}/rubin_md_match_object_shear.pq", + ) + stats_cols = [f"{b}_gaussFlux" for b in "griz"] + stats_cols += [f"{b}_gaussFluxErr" for b in "griz"] + drop_cols = [] + for col_prefix in ["tract", "ra", "dec", "gauss_snr", "patch_x", "patch_y", + "cell_idx_x", "cell_idx_y", "shear", "meta_step"]: + drop_cols += [f"{col_prefix}_{i}" for i in RUBIN_SHEAR_IDS] + rd_joined_rubin = hpmcm.reduceJoinedTable( + joined_rubin, + catalog_ids=RUBIN_SHEAR_IDS, + drop_cols=drop_cols, + stats_cols=stats_cols, + ) + tables_io.write(rd_joined_rubin, f"{TRACT}/rubin_md_match_object_joined.pq") + print(" Done.") + +# ── Step 3: DESC shear matching ─────────────────────────────────────────────── + +if REMATCH_DESC: + print("Running DESC shear matching...") + desc_shear_match = hpmcm.ShearMatch.createShearMatch( + geometry=OUR_GEOM, + pixel_match_scale=SHEAR_PIXEL_MATCH_SCALE, + pixel_r2_cut=SHEAR_PIXEL_R2_CUT, + deshear=-DESC_SHEAR, + shear_names=DESC_SHEAR_NAMES, + ) + desc_shear_match.reduceData(DESC_SHEAR_FILES, DESC_SHEAR_IDS) + if TEST_ONLY: + desc_shear_match.analysisLoop(_TEST_RANGE, _TEST_RANGE) + else: + desc_shear_match.analysisLoop() + desc_stats = desc_shear_match.extractStats() + desc_stats.update(**desc_shear_match.extractShearStats()) + tables_io.write(desc_stats, f"{TRACT}/desc_md_match_") + desc_shear_match.save(f"{TRACT}/desc_md_match") + print(" Done.") + +# ── Step 3 (cont.): Build joined DESC object table ──────────────────────────── + +if BUILD_DESC_JOINED: + print("Building joined DESC object table...") + desc_obj_stats = tables_io.read(f"{TRACT}/desc_md_match_object_stats.pq") + desc_obj_assoc = tables_io.read(f"{TRACT}/desc_md_match_object_assoc.pq") + joined_desc = hpmcm.output_tables.buildJoinedObjectTable( + desc_obj_stats, + desc_obj_assoc, + DESC_SHEAR_FILES, + DESC_SHEAR_IDS, + object_shear=f"{TRACT}/desc_md_match_object_shear.pq", + ) + stats_cols = [f"flux_{b}" for b in "riz"] + stats_cols += [f"flux_err_{b}" for b in "riz"] + drop_cols = [] + for col_prefix in ["ra", "dec", "gauss_snr", "patch_x", "patch_y", + "cell_idx_x", "cell_idx_y", "shear", "meta_step"]: + drop_cols += [f"{col_prefix}_{i}" for i in DESC_SHEAR_IDS] + rd_joined_desc = hpmcm.reduceJoinedTable( + joined_desc, + catalog_ids=DESC_SHEAR_IDS, + drop_cols=drop_cols, + stats_cols=stats_cols, + ) + tables_io.write(rd_joined_desc, f"{TRACT}/desc_md_match_object_joined.pq") + print(" Done.") + +# ── Step 4: WCS meta-matching ───────────────────────────────────────────────── + +if REMATCH_OBJECTS: + print("Running WCS meta-matching...") + wcs_match = hpmcm.WcsMatch.create( + ref_dir=REF_DIR, + region_size=(_WCS_REGION_DEG, _WCS_REGION_DEG), + pixel_size=WCS_PIXEL_SCALE, + cell_size=WCS_CELL_SIZE, + cell_buffer=WCS_CELL_BUFFER, + cell_max_object=WCS_CELL_MAX_OBJECT, + pixel_match_scale=WCS_PIXEL_MATCH_SCALE, + pixel_r2_cut=WCS_PIXEL_R2_CUT, + ) + wcs_match.reduceData(ALL_FILES, ALL_IDS) + if TEST_ONLY: + wcs_match.analysisLoop(_TEST_RANGE, _TEST_RANGE) + else: + wcs_match.analysisLoop() + stats = wcs_match.extractStats() + tables_io.write(stats, f"{TRACT}/all_match_") + wcs_match.save(f"{TRACT}/all_match") + print(" Done.") + +print("Pipeline complete.") From 578449712d2ae6885704204fae4479f30255e861 Mon Sep 17 00:00:00 2001 From: Eric Charles Date: Wed, 14 Oct 2026 17:07:01 -0700 Subject: [PATCH 5/6] Fix hierarchical object processing guard; update notebook outputs Add n_unique == n_src check before sub_mask.all() in heirarchicalProcessObject so that early return only fires when all sources are uniquely matched AND within the pixel r2 cut, not just when distances are small. Generated with AI Co-Authored-By: SLAC AI --- examples/rubin_dp2/CheckDescShearMatch.ipynb | 318 +++++++++++++- examples/rubin_dp2/CheckObjectMatch.ipynb | 76 ++++ examples/rubin_dp2/CheckRubinShearMatch.ipynb | 397 +++++++++++++++++- examples/rubin_dp2/RunMatching.ipynb | 76 ++-- examples/rubin_dp2/StudyObjectMatches.ipynb | 250 ++++++++++- examples/rubin_dp2/StudyRubinShearMatch.ipynb | 29 ++ src/hpmcm/match_utils.py | 2 +- 7 files changed, 1109 insertions(+), 39 deletions(-) diff --git a/examples/rubin_dp2/CheckDescShearMatch.ipynb b/examples/rubin_dp2/CheckDescShearMatch.ipynb index b0e0c68..df245b6 100644 --- a/examples/rubin_dp2/CheckDescShearMatch.ipynb +++ b/examples/rubin_dp2/CheckDescShearMatch.ipynb @@ -92,7 +92,7 @@ "outputs": [], "source": [ "# Join source catalog data onto the object stats table\n", - "joined_shear_obj = hpmcm.output_tables.buildJoinedObjectTable(\n", + "joined_shear_obj_full = hpmcm.output_tables.buildJoinedObjectTable(\n", " shear_desc_obj_stats,\n", " shear_desc_obj_assoc,\n", " input_files,\n", @@ -100,6 +100,27 @@ ")" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "3a19b6e4-a27d-4951-8886-da894882afaf", + "metadata": {}, + "outputs": [], + "source": [ + "center_mask = np.bitwise_and(np.fabs(joined_shear_obj_full.x_cent - 75) < 75, np.fabs(joined_shear_obj_full.y_cent - 75) < 75)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "63ee241e-df13-4fe2-8f90-0ffb7a43ed62", + "metadata": {}, + "outputs": [], + "source": [ + "joined_shear_obj = joined_shear_obj_full[center_mask]\n", + "joined_shear_obj_edge = joined_shear_obj_full[~center_mask]" + ] + }, { "cell_type": "code", "execution_count": null, @@ -110,6 +131,26 @@ "_ = plt.scatter(joined_shear_obj.y_cell_coadd_0, 3600*(joined_shear_obj.ra - joined_shear_obj.ra_0), s=1)" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "d329f707-8237-48dc-85ed-428025676345", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(joined_shear_obj_edge.x_cent, joined_shear_obj_edge.y_cent, s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "856c9a8a-48c1-472f-a110-6cbb51d6a542", + "metadata": {}, + "outputs": [], + "source": [ + "joined_shear_obj.x_cent" + ] + }, { "cell_type": "code", "execution_count": null, @@ -138,6 +179,281 @@ "id": "77d9bf1c-19f1-45c6-999c-65477ba8bc54", "metadata": {}, "outputs": [], + "source": [ + "cut = np.bitwise_and(joined_shear_obj.n_sib_objects==1, joined_shear_obj.n_sib_srcs>3)\n", + "print(cut.sum())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa3b6481-41df-4bf0-ab59-5fb4d43e53af", + "metadata": {}, + "outputs": [], + "source": [ + "good_mask = joined_shear_obj.catalog_mask == 7\n", + "bad_mask = joined_shear_obj.catalog_mask != 7\n", + "good_sub = joined_shear_obj[good_mask]\n", + "bad_sub = joined_shear_obj[bad_mask]\n", + "missing_ns_mask = np.bitwise_and(joined_shear_obj.catalog_mask, 1) != 1\n", + "missing_1m_mask = np.bitwise_and(joined_shear_obj.catalog_mask, 2) != 2\n", + "missing_1p_mask = np.bitwise_and(joined_shear_obj.catalog_mask, 4) != 4\n", + "missing_ns = joined_shear_obj[missing_ns_mask]\n", + "missing_1m = joined_shear_obj[missing_1m_mask]\n", + "missing_1p = joined_shear_obj[missing_1p_mask]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25e2f86f-5ada-4512-8062-30cf9c09642c", + "metadata": {}, + "outputs": [], + "source": [ + "print('good', np.bitwise_and(joined_shear_obj.n_src==3, joined_shear_obj.n_unique==3).sum())\n", + "print('conf', (joined_shear_obj.n_src != joined_shear_obj.n_unique).sum())\n", + "print('missing 1', np.bitwise_and(joined_shear_obj.n_src==2, joined_shear_obj.n_unique==2).sum())\n", + "print('missing 2', np.bitwise_and(joined_shear_obj.n_src==1, joined_shear_obj.n_unique==1).sum())\n", + "print('split', np.bitwise_and(joined_shear_obj.n_sib_objects>1, joined_shear_obj.n_sib_srcs==3).sum())\n", + "print('weird', np.bitwise_and(joined_shear_obj.n_sib_objects>1, joined_shear_obj.n_sib_srcs>3).sum())\n", + "print('weird', np.bitwise_and(joined_shear_obj.n_sib_objects>1, joined_shear_obj.n_sib_srcs<3).sum())\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e4f04843-b6c6-4c5e-9625-df11b637ae05", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(missing_ns.ra, missing_ns.dec, s=1)\n", + "_ = plt.scatter(missing_1m.ra, missing_1m.dec, s=1)\n", + "_ = plt.scatter(missing_1p.ra, missing_1p.dec, s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "be8d70e6-8172-408d-8a70-7f684c67bdb0", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(good_sub.n_sib_srcs, bins=np.linspace(0.5, 9.5, 10), alpha=0.5)\n", + "_ = plt.hist(bad_sub.n_sib_srcs, bins=np.linspace(0.5, 9.5, 10), alpha=0.5)\n", + "\n", + "_ = plt.yscale('log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "646d7c61-0f0f-4271-8e55-88e3d10b1a52", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(missing_ns.x_cent, bins=np.linspace(-5, 155, 16))\n", + "_ = plt.hist(missing_1m.x_cent, bins=np.linspace(-5, 155, 16))\n", + "_ = plt.hist(missing_1p.x_cent, bins=np.linspace(-5, 155, 16))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6fec2797-2fa6-4f23-bd22-f04f165dea58", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(good_sub.ra, good_sub.dec, s=1)\n", + "_ = plt.scatter(bad_sub.ra, bad_sub.dec, s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ce55fe23-9153-4c47-829d-2efaf8c4a2a6", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(good_sub.x_cent[::5], good_sub.y_cent[::5], s=1)\n", + "_ = plt.scatter(bad_sub.x_cent[::5], bad_sub.y_cent[::5], s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "95c4085a-79f8-45c1-acbd-09f7bd3e50ba", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist2d(3600*(good_sub.ra-good_sub.ra_0), 3600*(good_sub.dec-good_sub.dec_0),\n", + " bins=(np.linspace(-0.4, 0.4, 81), np.linspace(-0.4, 0.4, 81)), norm='log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e92fc228-1082-4b2e-b01b-1bbc23fabb26", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(good_sub.dist_rms, bins=np.linspace(0, 0.4, 201))\n", + "_ = plt.yscale('log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43bd12e8-3579-4213-86c2-8e4510e3cf16", + "metadata": {}, + "outputs": [], + "source": [ + "dist_0 = np.sqrt(good_sub.x_pix - good_sub.x_pix_0)**2 + (good_sub.y_pix - good_sub.y_pix_0)**2\n", + "dist_1 = np.sqrt(good_sub.x_pix - good_sub.x_pix_1)**2 + (good_sub.y_pix - good_sub.y_pix_1)**2\n", + "dist_2 = np.sqrt(good_sub.x_pix - good_sub.x_pix_2)**2 + (good_sub.y_pix - good_sub.y_pix_2)**2\n", + "\n", + "dist_0_b = np.sqrt(bad_sub.x_pix - bad_sub.x_pix_0)**2 + (bad_sub.y_pix - bad_sub.y_pix_0)**2\n", + "dist_1_b = np.sqrt(bad_sub.x_pix - bad_sub.x_pix_1)**2 + (bad_sub.y_pix - bad_sub.y_pix_1)**2\n", + "dist_2_b = np.sqrt(bad_sub.x_pix - bad_sub.x_pix_2)**2 + (bad_sub.y_pix - bad_sub.y_pix_2)**2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8bf2602-7e09-447d-b3a2-4f2f59c8858c", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(good_sub.y_cent, good_sub.x_pix - good_sub.x_pix_0, s=1)\n", + "_ = plt.scatter(good_sub.y_cent, good_sub.x_pix - good_sub.x_pix_1, s=1)\n", + "_ = plt.scatter(good_sub.y_cent, good_sub.x_pix - good_sub.x_pix_2, s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1b43922c-2582-4765-bef6-e458939401a5", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(good_sub.x_cent, good_sub.x_pix - good_sub.x_pix_0, s=1)\n", + "_ = plt.scatter(good_sub.x_cent, good_sub.x_pix - good_sub.x_pix_1, s=1)\n", + "_ = plt.scatter(good_sub.x_cent, good_sub.x_pix - good_sub.x_pix_2, s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db4c9ebb-8e36-424d-bc40-6ff5fd266749", + "metadata": {}, + "outputs": [], + "source": [ + "#_ = plt.hist2d(good_sub.x_cent, good_sub.x_pix - good_sub.x_pix_0, bins=(np.linspace(0, 150, 151), np.linspace(-2, 2, 101)), norm='log')\n", + "_ = plt.hist2d(good_sub.x_cent, good_sub.x_pix_0 - good_sub.x_pix_1, bins=(np.linspace(0, 150, 151), np.linspace(-2, 2, 101)), norm='log')\n", + "#_ = plt.hist2d(good_sub.x_cent, good_sub.x_pix_1 - good_sub.x_pix_2, bins=(np.linspace(0, 150, 151), np.linspace(-2, 2, 101)), norm='log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e3b32fc8-2ab4-43f7-9187-271659d2155c", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(good_sub.x_cent, good_sub.y_pix - good_sub.y_pix_0, s=1)\n", + "_ = plt.scatter(good_sub.x_cent, good_sub.y_pix - good_sub.y_pix_1, s=1)\n", + "_ = plt.scatter(good_sub.x_cent, good_sub.y_pix - good_sub.y_pix_2, s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cb7be02f-de30-4b53-a035-b899f8c4172d", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(dist_0, bins=np.linspace(0, 1.4, 201), alpha=0.2)\n", + "_ = plt.hist(dist_1, bins=np.linspace(0, 1.4, 201), alpha=0.2)\n", + "_ = plt.hist(dist_2, bins=np.linspace(0, 1.4, 201), alpha=0.2)\n", + "_ = plt.yscale('log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b37366ad-b065-48f2-b082-90434d2ee8f4", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(dist_0_b, bins=np.linspace(0, 1.4, 201), alpha=0.2)\n", + "_ = plt.hist(dist_1_b, bins=np.linspace(0, 1.4, 201), alpha=0.2)\n", + "_ = plt.hist(dist_2_b, bins=np.linspace(0, 1.4, 201), alpha=0.2)\n", + "_ = plt.yscale('log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7133154a-af09-4226-8ca3-d78b1bd84669", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist2d(3600*(bad_sub.ra-bad_sub.ra_0), 3600*(bad_sub.dec-bad_sub.dec_0),\n", + " bins=(np.linspace(-0.4, 0.4, 81), np.linspace(-0.4, 0.4, 81)), norm='log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b7d2a5c5-92ea-467d-801c-45044b9eab0e", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(shear_desc_ns.x_cell_coadd, bins=np.linspace(-5, 155, 16))\n", + "_ = plt.hist(shear_desc_1m.x_cell_coadd, bins=np.linspace(-5, 155, 16))\n", + "_ = plt.hist(shear_desc_1p.x_cell_coadd, bins=np.linspace(-5, 155, 16))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25b31837-479c-4793-8f4e-f19e1c757265", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(joined_shear_obj.x_cent, bins=np.linspace(-5, 155, 161))\n", + "_ = plt.hist(good_sub.x_cent, bins=np.linspace(-5, 155, 161))\n", + "_ = plt.hist(bad_sub.x_cent, bins=np.linspace(-5, 155, 161))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3de64923-7183-41ea-9e49-b787081cd036", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(good_sub.y_cent, bins=np.linspace(-5, 155, 161))\n", + "_ = plt.hist(bad_sub.y_cent, bins=np.linspace(-5, 155, 161))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d34f1dc2-f525-41da-97df-b0847fd95f88", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(joined_shear_obj.catalog_mask, bins=np.linspace(0.5, 7.5, 8), alpha=0.2)\n", + "_ = plt.hist(joined_shear_obj_edge.catalog_mask, bins=np.linspace(0.5, 7.5, 8), alpha=0.2)\n", + "_ = plt.yscale('log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c6c3eef-0314-48b6-bb06-a75bab3d2ab4", + "metadata": {}, + "outputs": [], "source": [] } ], diff --git a/examples/rubin_dp2/CheckObjectMatch.ipynb b/examples/rubin_dp2/CheckObjectMatch.ipynb index 07c65ab..f7f5b74 100644 --- a/examples/rubin_dp2/CheckObjectMatch.ipynb +++ b/examples/rubin_dp2/CheckObjectMatch.ipynb @@ -124,6 +124,82 @@ "id": "5e1e3368-79da-4235-968e-7b03587ec496", "metadata": {}, "outputs": [], + "source": [ + "cut = np.bitwise_and(joined_obj.n_sib_objects==1, joined_obj.n_sib_srcs>4)\n", + "print(cut.sum())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1e33d58d-2449-4375-947f-090496c3c7c0", + "metadata": {}, + "outputs": [], + "source": [ + "good_mask = joined_obj.catalog_mask == 15\n", + "bad_mask = joined_obj.catalog_mask != 15\n", + "good_sub = joined_obj[good_mask]\n", + "bad_sub = joined_obj[bad_mask]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36d436aa-0979-4ec8-8e38-0c9f3312277f", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(bad_sub.ra[::5], bad_sub.dec[::5], s=1)\n", + "_ = plt.scatter(good_sub.ra[::5], good_sub.dec[::5], s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d1cf6846-64b0-4f2c-a4b1-85174067c94c", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(bad_sub.x_cent[::5], bad_sub.y_cent[::5], s=1)\n", + "_ = plt.scatter(good_sub.x_cent[::5], good_sub.y_cent[::5], s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d9829e52-a386-4070-aaa2-9a7e18de47cb", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist2d(3600*(good_sub.ra-good_sub.ra_0), 3600*(good_sub.dec-good_sub.dec_0),\n", + " bins=(np.linspace(-0.8, 0.8, 81), np.linspace(-0.8, 0.8, 81)), norm='log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b0402d99-8b0c-4b52-b219-9688c7b36457", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(bad_sub.x_cent, bins=np.linspace(-5, 155, 161))\n", + "_ = plt.hist(good_sub.x_cent, bins=np.linspace(-5, 155, 161))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "254cab26-05c1-4f15-8ac1-7f66a8bd9d17", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "35f8b73a-63cb-40fa-8ded-85f738b4021a", + "metadata": {}, + "outputs": [], "source": [] } ], diff --git a/examples/rubin_dp2/CheckRubinShearMatch.ipynb b/examples/rubin_dp2/CheckRubinShearMatch.ipynb index eb83104..d35bb46 100644 --- a/examples/rubin_dp2/CheckRubinShearMatch.ipynb +++ b/examples/rubin_dp2/CheckRubinShearMatch.ipynb @@ -59,6 +59,26 @@ "catalog_ids = [0, 1, 2, 3, 4]" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "d7ff674d-b4a3-4f2a-b192-c8ddf5b026c3", + "metadata": {}, + "outputs": [], + "source": [ + "(shear_rubin_obj_stats.n_sib_objects == 2).sum()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cd9a5d68-07da-4d01-9ac4-e75b8f868b01", + "metadata": {}, + "outputs": [], + "source": [ + "np.bincount(shear_rubin_obj_stats.n_sib_srcs)" + ] + }, { "cell_type": "markdown", "id": "8b62b38e", @@ -94,7 +114,7 @@ "metadata": {}, "outputs": [], "source": [ - "joined_shear_obj = hpmcm.output_tables.buildJoinedObjectTable(\n", + "joined_shear_obj_full = hpmcm.output_tables.buildJoinedObjectTable(\n", " shear_rubin_obj_stats,\n", " shear_rubin_obj_assoc,\n", " input_files,\n", @@ -102,6 +122,49 @@ ")" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "0aade12d-4d75-4347-b859-c31bd3f19e5f", + "metadata": {}, + "outputs": [], + "source": [ + "(joined_shear_obj_full.n_src > joined_shear_obj_full.n_unique).sum()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b0bfc1f7-d092-4d14-acf7-fc33613a1c4f", + "metadata": {}, + "outputs": [], + "source": [ + "center_mask = np.bitwise_and(np.fabs(joined_shear_obj_full.x_cent - 75) < 75, np.fabs(joined_shear_obj_full.y_cent - 75) < 75)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cfaeb55c-7742-4e70-bb60-1b0558a757ad", + "metadata": {}, + "outputs": [], + "source": [ + "joined_shear_obj = joined_shear_obj_full[center_mask]\n", + "joined_shear_obj_edge = joined_shear_obj_full[~center_mask]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4fdb3a24-a83b-49be-9e9c-33fb395465c8", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(joined_shear_obj.catalog_mask, bins=np.linspace(0.5, 31.5, 32))\n", + "_ = plt.hist(joined_shear_obj_edge.catalog_mask, bins=np.linspace(0.5, 31.5, 32))\n", + "_ = plt.yscale('log')" + ] + }, { "cell_type": "code", "execution_count": null, @@ -140,6 +203,338 @@ "id": "089156b6-d083-41f2-8d11-3b8409d06d3e", "metadata": {}, "outputs": [], + "source": [ + "cut = np.bitwise_and(joined_shear_obj.n_sib_objects==1, joined_shear_obj.n_sib_srcs>5)\n", + "print(cut.sum())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cfd9a011-b6d4-4716-b71e-d40f1a6c8387", + "metadata": {}, + "outputs": [], + "source": [ + "print('good', np.bitwise_and(joined_shear_obj.n_src==5, joined_shear_obj.n_unique==5).sum())\n", + "print('conf', (joined_shear_obj.n_src != joined_shear_obj.n_unique).sum())\n", + "print('missing 1', np.bitwise_and(joined_shear_obj.n_src==4, joined_shear_obj.n_unique==4).sum())\n", + "print('missing 2', np.bitwise_and(joined_shear_obj.n_src==3, joined_shear_obj.n_unique==3).sum())\n", + "print('missing 3', np.bitwise_and(joined_shear_obj.n_src==2, joined_shear_obj.n_unique==2).sum())\n", + "print('missing 4', np.bitwise_and(joined_shear_obj.n_src==1, joined_shear_obj.n_unique==1).sum())\n", + "print('split', np.bitwise_and(joined_shear_obj.n_sib_objects>1, joined_shear_obj.n_sib_srcs==5).sum())\n", + "print('weird', np.bitwise_and(joined_shear_obj.n_sib_objects>1, joined_shear_obj.n_sib_srcs>5).sum())\n", + "print('weird', np.bitwise_and(joined_shear_obj.n_sib_objects>1, joined_shear_obj.n_sib_srcs<5).sum())\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "147f3e4b-d539-41a4-a04d-4fc5ba51a6af", + "metadata": {}, + "outputs": [], + "source": [ + "good_mask = joined_shear_obj.catalog_mask == 31\n", + "bad_mask = joined_shear_obj.catalog_mask != 31\n", + "good_sub = joined_shear_obj[good_mask]\n", + "bad_sub = joined_shear_obj[bad_mask]\n", + "\n", + "missing_1_mask = np.bitwise_and(joined_shear_obj.n_src==4, joined_shear_obj.n_unique==4)\n", + "missing_2_mask = np.bitwise_and(joined_shear_obj.n_src==3, joined_shear_obj.n_unique==3)\n", + "missing_3_mask = np.bitwise_and(joined_shear_obj.n_src==2, joined_shear_obj.n_unique==2)\n", + "missing_4_mask = np.bitwise_and(joined_shear_obj.n_src==1, joined_shear_obj.n_unique==1)\n", + "\n", + "missing_1 = joined_shear_obj[missing_1_mask]\n", + "missing_2 = joined_shear_obj[missing_2_mask]\n", + "missing_3 = joined_shear_obj[missing_3_mask]\n", + "missing_4 = joined_shear_obj[missing_4_mask]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bc697d06-60cd-4383-bc0f-da1f903b9f30", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2475ad1-9a19-4d2e-89d5-8e710cc0ef61", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(good_sub.ra, good_sub.dec, s=1)\n", + "_ = plt.scatter(bad_sub.ra, bad_sub.dec, s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c53016c7-a033-4f08-b92a-347a398a915d", + "metadata": {}, + "outputs": [], + "source": [ + "#_ = plt.scatter(good_sub.x_cent[::5], good_sub.y_cent[::5], s=1)\n", + "#_ = plt.scatter(bad_sub.x_cent[::5], bad_sub.y_cent[::5], s=1)\n", + "#_ = plt.scatter(missing_1.x_cent, missing_1.y_cent, s=1)\n", + "_ = plt.scatter(missing_2.x_cent, missing_2.y_cent, s=1)\n", + "_ = plt.scatter(missing_3.x_cent, missing_3.y_cent, s=1)\n", + "#_ = plt.scatter(missing_4.x_cent, missing_4.y_cent, s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3e076b70-8301-4df8-a7f4-8bf6d69c3a79", + "metadata": {}, + "outputs": [], + "source": [ + "missing_2_217 = missing_2[missing_2.cell_idx == 217]\n", + "missing_3_217 = missing_3[missing_3.cell_idx == 217]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b49898f4-b532-4661-a065-eee3c252bccb", + "metadata": {}, + "outputs": [], + "source": [ + "joined_shear_obj[joined_shear_obj.cell_idx == 217]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c07de349-b5ba-4cef-974c-b2619cc51b4f", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(joined_shear_obj.cell_idx_x_0, joined_shear_obj.cell_idx_y_0, s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21a58f53-99b2-4d71-966f-98320849a9b4", + "metadata": {}, + "outputs": [], + "source": [ + "(joined_shear_obj.cell_idx - (201*joined_shear_obj.cell_idx_x_0 + joined_shear_obj.cell_idx_y_0)).max()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47a4ec04-2d4f-4e5b-be1b-be9e6be8fec1", + "metadata": {}, + "outputs": [], + "source": [ + "((joined_shear_obj.cell_idx // 201 ) -1) % 20" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b2e027b8-2a02-4065-9a7c-a6ac7ef40172", + "metadata": {}, + "outputs": [], + "source": [ + "missing_2_217" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd0b8741-79a6-4c46-83fc-74705551d03e", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d28165f-f064-420a-a95e-c056bfadf984", + "metadata": {}, + "outputs": [], + "source": [ + "missing_3_217" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "286c7ef1-b478-4cfa-8244-cf244f755a8c", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d306816f-5f04-46fb-a7be-c35c27b5a959", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(missing_1.ra, missing_1.dec, s=1)\n", + "_ = plt.scatter(missing_4.ra, missing_4.dec, s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0809bb22-b2c9-4e37-a66b-a7dd25e218e6", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist2d(3600*(good_sub.ra-good_sub.ra_1), 3600*(good_sub.dec-good_sub.dec_1),\n", + " bins=(np.linspace(-0.4, 0.4, 81), np.linspace(-0.4, 0.4, 81)), norm='log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2c389f2-f90d-4a55-b336-b389c37bb1bc", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist((bad_sub.x_cent-75)**2 + (bad_sub.y_cent-75)**2, bins=np.linspace(-5, 2*75*75, 16))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4fb0e7e9-a52e-43f8-b383-6896ef11b858", + "metadata": {}, + "outputs": [], + "source": [ + "#_ = plt.hist(joined_shear_obj.x_cent, bins=np.linspace(5, 145, 15))\n", + "_ = plt.hist(missing_1.x_cent, bins=np.linspace(5, 145, 15), alpha=0.2, label='missing 1')\n", + "_ = plt.hist(missing_2.x_cent, bins=np.linspace(5, 145, 15), alpha=0.2, label='missing 2')\n", + "_ = plt.hist(missing_3.x_cent, bins=np.linspace(5, 145, 15), alpha=0.2, label='missing 3')\n", + "_ = plt.hist(missing_4.x_cent, bins=np.linspace(5, 145, 15), alpha=0.2, label='missing 4')\n", + "#_ = plt.yscale('log')\n", + "_ = plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b870f204-875a-4b1e-854b-371c778569eb", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14807fec-6478-4b81-a65d-9e6f2c908a05", + "metadata": {}, + "outputs": [], + "source": [ + "the_hist, the_edges = np.histogram(shear_rubin_ns.x_cell_coadd, bins=np.linspace(0, 150, 16))\n", + "the_hist += np.histogram(shear_rubin_1m.x_cell_coadd, bins=np.linspace(0, 150, 16))[0]\n", + "the_hist += np.histogram(shear_rubin_1p.x_cell_coadd, bins=np.linspace(0, 150, 16))[0]\n", + "the_hist += np.histogram(shear_rubin_2m.x_cell_coadd, bins=np.linspace(0, 150, 16))[0]\n", + "the_hist += np.histogram(shear_rubin_2p.x_cell_coadd, bins=np.linspace(0, 150, 16))[0]\n", + "_ = plt.stairs(the_hist, the_edges)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "85d437ae-43fa-4647-8712-fb8a47510bba", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "74298cfe-9bbb-43b5-ab78-87180ab989f2", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd51076f-5ca8-4558-bd50-d3ed747e2ea2", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(good_sub.x_cent, bins=np.linspace(-5, 155, 161))\n", + "_ = plt.hist(bad_sub.x_cent, bins=np.linspace(-5, 155, 161))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "33fb39e8-ff21-4bd8-b3e6-a9d5d50bc433", + "metadata": {}, + "outputs": [], + "source": [ + "bad_sub_single_mask = bad_sub.n_sib_objects == 1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b700b851-9096-4416-9ee9-34d5761b9fee", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(bad_sub.x_cent[bad_sub_single_mask], bins=np.linspace(-5, 155, 161))\n", + "_ = plt.hist(bad_sub.x_cent[~bad_sub_single_mask], bins=np.linspace(-5, 155, 161))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "66f405e8-df31-4186-94b5-c5917c785474", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(good_sub.y_cent, bins=np.linspace(-5, 15, 161))\n", + "_ = plt.hist(bad_sub.y_cent, bins=np.linspace(-5, 15, 161))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02218696-2595-4ceb-8c2b-c28ab3e02fc5", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist2d(good_sub.x_pix-good_sub.x_pix_0, good_sub.y_pix-good_sub.y_pix_0,\n", + " bins=(np.linspace(-4.0, 4.0, 81), np.linspace(-4.0, 4.0, 81)), norm='log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8f3c6afe-be71-4f5a-aa7d-0065766faf2d", + "metadata": {}, + "outputs": [], + "source": [ + "good_sub.x_cent-good_sub.x_cell_coadd_0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0798481b-ee26-458d-b58f-e4d54f04b6dc", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "002ccfb8-5023-48ec-8f5c-80395e44467b", + "metadata": {}, + "outputs": [], "source": [] } ], diff --git a/examples/rubin_dp2/RunMatching.ipynb b/examples/rubin_dp2/RunMatching.ipynb index 30c1b8a..8693f05 100644 --- a/examples/rubin_dp2/RunMatching.ipynb +++ b/examples/rubin_dp2/RunMatching.ipynb @@ -35,7 +35,9 @@ "tract = 8006\n", "# This will need to change with tract\n", "ref_dir=(340.5000000000, -9.6694214876)\n", - "our_geom = hpmcm.shear_utils.ShearCellGeometry(ref_dir=ref_dir)" + "our_geom = hpmcm.shear_utils.ShearCellGeometry(\n", + " ref_dir=ref_dir,\n", + ")" ] }, { @@ -92,13 +94,15 @@ "outputs": [], "source": [ "# Allow for re-running various matchers\n", - "rematch_rubin_shear = False\n", - "rematch_desc_shear = False\n", + "rematch_rubin_shear = True\n", + "rematch_desc_shear = True\n", "rematch_objects = True\n", "\n", "if rematch_rubin_shear:\n", " rubin_shear_match = hpmcm.shear_match.ShearMatch.createShearMatch(\n", " geometry=our_geom,\n", + " pixel_match_scale=2,\n", + " pixel_r2_cut=16,\n", " deshear=-0.01,\n", " shear_names=['ns', '1m', '1p', '2m', '2p'], \n", " )\n", @@ -117,6 +121,8 @@ " desc_shear_match = hpmcm.shear_match.ShearMatch.createShearMatch(\n", " geometry=our_geom,\n", " deshear=-0.01,\n", + " pixel_match_scale=2,\n", + " pixel_r2_cut=16, \n", " shear_names=['ns', '1m', '1p'],\n", " )\n", " desc_shear_match.reduceData(desc_shear_match_files, desc_shear_match_ids)\n", @@ -211,15 +217,16 @@ "outputs": [], "source": [ "# Draw a cluster from the rubin shear matching\n", - "cell_data = rubin_shear_match.cell_dict[20200]\n", - "cluster_id = list(cell_data.cluster_dict.keys())[1]\n", - "cluster_data = cell_data.cluster_dict[cluster_id]\n", - "img = cell_data.countsMap()\n", - "_ = hpmcm.viz_utils.showCluster(\n", - " img,\n", - " cluster_data,\n", - " cell_data,\n", - ")" + "if rematch_rubin_shear:\n", + " cell_data = rubin_shear_match.cell_dict[20200]\n", + " cluster_id = list(cell_data.cluster_dict.keys())[1]\n", + " cluster_data = cell_data.cluster_dict[cluster_id]\n", + " img = cell_data.countsMap()\n", + " _ = hpmcm.viz_utils.showCluster(\n", + " img,\n", + " cluster_data,\n", + " cell_data,\n", + " )" ] }, { @@ -230,15 +237,16 @@ "outputs": [], "source": [ "# Draw a cluster from the desc shear matching\n", - "cell_data = desc_shear_match.cell_dict[20200]\n", - "cluster_id = list(cell_data.cluster_dict.keys())[1]\n", - "cluster_data = cell_data.cluster_dict[cluster_id]\n", - "img = cell_data.countsMap()\n", - "_ = hpmcm.viz_utils.showCluster(\n", - " img,\n", - " cluster_data,\n", - " cell_data,\n", - ")" + "if rematch_desc_shear:\n", + " cell_data = desc_shear_match.cell_dict[20200]\n", + " cluster_id = list(cell_data.cluster_dict.keys())[1]\n", + " cluster_data = cell_data.cluster_dict[cluster_id]\n", + " img = cell_data.countsMap()\n", + " _ = hpmcm.viz_utils.showCluster(\n", + " img,\n", + " cluster_data,\n", + " cell_data,\n", + " )" ] }, { @@ -251,14 +259,17 @@ "if rematch_objects:\n", " pixel_scale = 0.2\n", " n_pix_region = 30000\n", + " cell_size = 150\n", " region_size_deg = n_pix_region*pixel_scale/3600\n", " wcs_match = hpmcm.WcsMatch.create(\n", " ref_dir=ref_dir,\n", " region_size=(region_size_deg, region_size_deg),\n", " pixel_size=pixel_scale/3600,\n", - " cell_size=150,\n", + " cell_size=cell_size,\n", " cell_buffer=5,\n", " cell_max_object=10000,\n", + " pixel_match_scale=4,\n", + " pixel_r2_cut=16, \n", " )\n", " \n", " wcs_match.reduceData(all_files, all_ids)\n", @@ -279,21 +290,22 @@ "outputs": [], "source": [ "# Draw a cluster from the meta matching\n", - "cell_data = wcs_match.cell_dict[20200]\n", - "cluster_id = list(cell_data.cluster_dict.keys())[0]\n", - "cluster_data = cell_data.cluster_dict[cluster_id]\n", - "img = cell_data.countsMap()\n", - "_ = hpmcm.viz_utils.showCluster(\n", - " img,\n", - " cluster_data,\n", - " cell_data,\n", - ")" + "if rematch_objects:\n", + " cell_data = wcs_match.cell_dict[20200]\n", + " cluster_id = list(cell_data.cluster_dict.keys())[0]\n", + " cluster_data = cell_data.cluster_dict[cluster_id]\n", + " img = cell_data.countsMap()\n", + " _ = hpmcm.viz_utils.showCluster(\n", + " img,\n", + " cluster_data,\n", + " cell_data,\n", + " )" ] }, { "cell_type": "code", "execution_count": null, - "id": "f588b1fc-64ab-4d4d-8d6c-d4b8e6d63e3e", + "id": "ee7ec80d-7681-44b5-92bf-5f5a4a82e6fa", "metadata": {}, "outputs": [], "source": [] diff --git a/examples/rubin_dp2/StudyObjectMatches.ipynb b/examples/rubin_dp2/StudyObjectMatches.ipynb index efef92b..076715f 100644 --- a/examples/rubin_dp2/StudyObjectMatches.ipynb +++ b/examples/rubin_dp2/StudyObjectMatches.ipynb @@ -60,8 +60,8 @@ "all_files = [\n", " f'{tract}/object_reduced_{tract}.parq',\n", " f'{tract}/anacal_reduced_{tract}.parq',\n", - " f'{tract}/desc_md_match_object_stats.pq',\n", - " f'{tract}/rubin_md_match_object_stats.pq',\n", + " f'{tract}/desc_md_match_object_joined.pq',\n", + " f'{tract}/rubin_md_match_object_joined.pq',\n", "]\n", "all_ids = list(range(len(all_files)))\n" ] @@ -82,6 +82,42 @@ ")" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "87917500-1d90-4281-964a-168fec2fb8da", + "metadata": {}, + "outputs": [], + "source": [ + "cut = np.bitwise_and(joined_shear_obj.n_sib_objects==1, joined_shear_obj.n_sib_srcs>4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19eb9204-60f6-4c17-af6d-3950a503c444", + "metadata": {}, + "outputs": [], + "source": [ + "a = joined_shear_obj[cut].iloc[0]\n", + "b = joined_shear_obj[cut].iloc[1]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b15571a3-726c-427e-8c04-f2cf95a42045", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "for k, x in a.items():\n", + " y = b[k]\n", + " if np.isfinite(x) and x != y:\n", + " print(k, x, y)" + ] + }, { "cell_type": "code", "execution_count": null, @@ -148,7 +184,164 @@ { "cell_type": "code", "execution_count": null, - "id": "533cb76e-aaf8-49d7-a11f-4d1740645cd9", + "id": "0d97ae89-8b30-4fae-8ad1-ad1895f2733a", + "metadata": {}, + "outputs": [], + "source": [ + "good_sub.r_gaapPsfFlux_0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "573d6e2f-8350-4b54-bc34-b3291db98004", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(bad_sub.r_gaapPsfFlux_0, bins=np.logspace(0, 6, 121))\n", + "_ = plt.hist(good_sub.r_gaapPsfFlux_0, bins=np.logspace(0, 6, 121))\n", + "_ = plt.xscale('log')\n", + "_ = plt.yscale('log')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22641fe8-3d48-4ed4-a325-f7a2e88d8a46", + "metadata": {}, + "outputs": [], + "source": [ + "def flux_to_mags(flux):\n", + " vals = -2.5 * np.log10(flux) + 31.4\n", + " return np.where(np.isfinite(vals), vals, np.nan)\n", + "\n", + "#_ = plt.scatter(flux_to_mags(good_sub.g_gaapPsfFlux_0), flux_to_mags(good_sub.g_gaapPsfFlux_0)-flux_to_mags(good_sub.flux_g_mean_2), s=1)\n", + "_ = plt.scatter(flux_to_mags(good_sub.r_gaapPsfFlux_0), flux_to_mags(good_sub.r_gaapPsfFlux_0)-flux_to_mags(good_sub.flux_r_mean_2), s=1)\n", + "_ = plt.scatter(flux_to_mags(good_sub.i_gaapPsfFlux_0), flux_to_mags(good_sub.i_gaapPsfFlux_0)-flux_to_mags(good_sub.flux_i_mean_2), s=1)\n", + "_ = plt.scatter(flux_to_mags(good_sub.z_gaapPsfFlux_0), flux_to_mags(good_sub.z_gaapPsfFlux_0)-flux_to_mags(good_sub.flux_z_mean_2), s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d89d187f-a076-4972-8edc-f34c73d2a554", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(flux_to_mags(good_sub.r_gaapPsfFlux_0)-flux_to_mags(good_sub.flux_r_mean_2), flux_to_mags(good_sub.i_gaapPsfFlux_0)-flux_to_mags(good_sub.flux_i_mean_2), s=1)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1a7c47d6-6f38-4081-90ee-81d9606b4b7f", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(bad_sub.x_pix, bad_sub.y_pix, s=1)\n", + "_ = plt.scatter(good_sub.x_pix, good_sub.y_pix, s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d4eeaab1-41b0-4e24-a1be-e503903149a6", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(bad_sub.r_gaapPsfFlux_0, bad_sub.i_gaapPsfFlux_0, s=1)\n", + "_ = plt.scatter(good_sub.r_gaapPsfFlux_0, good_sub.i_gaapPsfFlux_0, s=1)\n", + "_ = plt.xscale('log')\n", + "_ = plt.yscale('log')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05137ac3-dff7-4727-a428-6b50583c0ea8", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(3600*(good_sub.ra_0-good_sub.ra_3), 3600*(good_sub.dec_0-good_sub.dec_3), s=1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fb7ed887-c467-4c36-aa38-912e6b213a98", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.scatter(bad_sub.n_sib_objects, bad_sub.n_sib_srcs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd5a30f6-8ad8-425a-9d6d-d00125fb8376", + "metadata": {}, + "outputs": [], + "source": [ + "cut = np.bitwise_and(good_sub.n_sib_objects==1, good_sub.n_sib_srcs>4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16884f47-1627-465f-b5bc-d56112a81961", + "metadata": {}, + "outputs": [], + "source": [ + "good_sub[cut]." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b25c5e40-67b8-4b2e-9e28-94049ce5314b", + "metadata": {}, + "outputs": [], + "source": [ + "np.bincount(bad_sub.n_sib_objects.astype(int))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "87b28bdf-9862-4adb-a669-1932ed3a3085", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(good_sub.delta_g_1_1_3, bins=np.linspace(-1., 1, 201))\n", + "_ = plt.hist(good_sub.delta_g_2_2_3, bins=np.linspace(-1., 1, 201))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "982a0d95-5595-4aff-b31e-d96e94cad5cc", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "_ = plt.hist2d(good_sub.delta_g_1_1_2, good_sub.delta_g_1_1_3, bins=(np.linspace(-0.2, 0.2, 201), np.linspace(-0.2, 0.2, 201)), norm='log')\n", + "#_ = plt.hist2d(bad_sub.delta_g_1_1_2, bad_sub.delta_g_1_1_3, bins=(np.linspace(-0.2, 0.2, 201), np.linspace(-0.2, 0.2, 201)), norm='log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aefd2106-dd00-4d1e-a2c7-985f5ef3c8f5", + "metadata": {}, + "outputs": [], + "source": [ + "list(good_sub.columns)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "436a1f37-dc8b-4a93-b514-2f98d81780a3", "metadata": {}, "outputs": [], "source": [] @@ -156,7 +349,56 @@ { "cell_type": "code", "execution_count": null, - "id": "add1db1c-e59e-434e-a340-ee0ac8fd5b4b", + "id": "56ee292b-04a4-4dc0-8839-758c4b6661ae", + "metadata": {}, + "outputs": [], + "source": [ + "assoc_matrix = np.zeros((4,4))\n", + "for i in range(4):\n", + " for j in range(4):\n", + " check = np.bitwise_or(1 << i, 1 << j)\n", + " assoc_matrix[i, j] = (np.bitwise_and(joined_shear_obj.catalog_mask, check) == check).sum()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e06f22ef-45db-4277-8323-6968ae021619", + "metadata": {}, + "outputs": [], + "source": [ + "rubin_not_object_mask = np.bitwise_and(joined_shear_obj.catalog_mask, 9) == 8\n", + "rubin_and_object_mask = np.bitwise_and(joined_shear_obj.catalog_mask, 9) == 9\n", + "rubin_not_object = joined_shear_obj[rubin_not_object_mask]\n", + "rubin_and_object = joined_shear_obj[rubin_and_object_mask]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2845af38-b3c0-4f2f-b14e-d5e14df7b97f", + "metadata": {}, + "outputs": [], + "source": [ + "np.bincount(rubin_and_object.catalog_mask_3.astype(int))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bede5b8f-5c04-4dd9-afd5-a06c0d819a22", + "metadata": {}, + "outputs": [], + "source": [ + "_ = plt.hist(rubin_and_object.catalog_mask_3, bins=np.linspace(0.5, 31.5, 32), alpha=0.1)\n", + "_ = plt.hist(rubin_not_object.catalog_mask_3, bins=np.linspace(0.5, 31.5, 32), alpha=0.1)\n", + "_ = plt.yscale('log')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a10d2769-15b6-4311-b69a-d58bb4c9e2e1", "metadata": {}, "outputs": [], "source": [] diff --git a/examples/rubin_dp2/StudyRubinShearMatch.ipynb b/examples/rubin_dp2/StudyRubinShearMatch.ipynb index 7fc07df..d8bba39 100644 --- a/examples/rubin_dp2/StudyRubinShearMatch.ipynb +++ b/examples/rubin_dp2/StudyRubinShearMatch.ipynb @@ -126,6 +126,27 @@ "bad_sub = joined_shear_obj[bad_mask]\n" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "a01c2a34-8403-41ab-86a1-7db6294675e7", + "metadata": {}, + "outputs": [], + "source": [ + "cut = np.bitwise_and(joined_shear_obj.n_sib_objects==1, joined_shear_obj.n_sib_srcs>5)\n", + "print(cut.sum())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "83583473-4ac7-4289-b038-8a52ddbd524f", + "metadata": {}, + "outputs": [], + "source": [ + "joined_shear_obj[cut].recurse" + ] + }, { "cell_type": "code", "execution_count": null, @@ -177,6 +198,14 @@ "metadata": {}, "outputs": [], "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d44b9009-15b3-4de4-96e5-94761fcddb0f", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/src/hpmcm/match_utils.py b/src/hpmcm/match_utils.py index 44a41ee..62a2164 100644 --- a/src/hpmcm/match_utils.py +++ b/src/hpmcm/match_utils.py @@ -72,7 +72,7 @@ def heirarchicalProcessObject( obj_data.snr_rms = np.std(obj_data.data.snr.values) sub_mask = obj_data.dist_2 < pixel_r2_cut - if sub_mask.all(): + if (obj_data.n_unique == obj_data.n_src) and sub_mask.all(): return if recurse >= cell_data.matcher.max_sub_division: From 9b19e95e70b047cdeb180f30045b8f682e2326c9 Mon Sep 17 00:00:00 2001 From: Eric Charles Date: Thu, 15 Oct 2026 13:18:52 -0700 Subject: [PATCH 6/6] Split reduceObjectTable and reduceAnacalTable into DF + file-I/O layers Add reduceObjectDF and reduceAnacalDF that operate on DataFrames and return DataFrames. The existing Table wrappers now delegate to these, keeping file I/O separate from the transformation logic. Generated with AI Co-Authored-By: SLAC AI --- src/hpmcm/utils.py | 66 ++++++++++++++++++++++++++++------------------ 1 file changed, 41 insertions(+), 25 deletions(-) diff --git a/src/hpmcm/utils.py b/src/hpmcm/utils.py index b172fe6..dd4ed96 100644 --- a/src/hpmcm/utils.py +++ b/src/hpmcm/utils.py @@ -211,23 +211,11 @@ def associateSourcesToFootprints( return [findClusterIds(df, cluster_key, pixel_match_scale) for df in data] -def reduceObjectTable( - basefile: str, - outfile: str, +def reduceObjectDF( + t: pandas.DataFrame, extra_cols: list[str] | None = None, -) -> None: # pragma: no cover - """Reduce an object table to just the colums needed for matching - - Parameters - ---------- - basefile: - Original file name - - outfile: - Output file name - - extra_cols: - Extra columns to copy +) -> pandas.DataFrame: + """Reduce an object DataFrame to just the columns needed for matching. Notes ----- @@ -254,16 +242,13 @@ def reduceObjectTable( +-----------------------+---------------------------------------------------------------+ """ - t = tables_io.read(basefile) cols = ["tract", "patch", "coord_ra", "coord_dec", "objectId"] cols += [f"{band}_gaapPsfFlux" for band in "ugrizy"] cols += [f"{band}_gaapPsfFluxErr" for band in "ugrizy"] - if extra_cols is not None: cols += extra_cols tout = t[cols].copy(deep=True) - tout["ra"] = tout["coord_ra"] tout["dec"] = tout["coord_dec"] @@ -276,16 +261,15 @@ def reduceObjectTable( ) tout["snr"] = snr tout["id"] = tout["objectId"] + return tout - tout.to_parquet(outfile) - -def reduceAnacalTable( +def reduceObjectTable( basefile: str, outfile: str, extra_cols: list[str] | None = None, ) -> None: # pragma: no cover - """Reduce an object table to just the colums needed for matching + """Read an object catalog, reduce it, and write the result. Parameters ---------- @@ -297,6 +281,17 @@ def reduceAnacalTable( extra_cols: Extra columns to copy + """ + t = tables_io.read(basefile) + tout = reduceObjectDF(t, extra_cols) + tout.to_parquet(outfile) + + +def reduceAnacalDF( + t: pandas.DataFrame, + extra_cols: list[str] | None = None, +) -> pandas.DataFrame: + """Reduce an anacal DataFrame to just the columns needed for matching. Notes ----- @@ -323,7 +318,6 @@ def reduceAnacalTable( +-----------------------------+---------------------------------------------------------------+ """ - t = tables_io.read(basefile) cols = ["tract_id", "patch_x", "patch_y", "ra", "dec", "object_id"] cols += [f"lsst_{band}_flux_gauss2" for band in "riz"] cols += [f"lsst_{band}_flux_gauss2_err" for band in "riz"] @@ -331,7 +325,6 @@ def reduceAnacalTable( cols += extra_cols tout = t[cols].copy(deep=True) - tout["tract"] = tout["tract_id"] snr = np.zeros(len(tout["lsst_i_flux_gauss2_err"])) @@ -344,6 +337,29 @@ def reduceAnacalTable( tout["snr"] = snr tout["id"] = tout["object_id"] tout.drop(columns=["object_id"]) + return tout + + +def reduceAnacalTable( + basefile: str, + outfile: str, + extra_cols: list[str] | None = None, +) -> None: # pragma: no cover + """Read an anacal catalog, reduce it, and write the result. + + Parameters + ---------- + basefile: + Original file name + + outfile: + Output file name + + extra_cols: + Extra columns to copy + """ + t = tables_io.read(basefile) + tout = reduceAnacalDF(t, extra_cols) tout.to_parquet(outfile)