Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 4 additions & 28 deletions src/hpmcm/cell.py
Original file line number Diff line number Diff line change
Expand Up @@ -232,7 +232,7 @@ def _newObject(
return ObjectData(cluster, object_id, mask)

def _emptyCountsMaps(self) -> np.ndarray:
to_fill = np.zeros(np.ceil(self.n_pix).astype(int))
to_fill = np.zeros(np.ceil(self.n_pix / self.pixel_match_scale).astype(int))
return to_fill

def _singleCatalogCountsMap(
Expand All @@ -242,6 +242,7 @@ def _singleCatalogCountsMap(
df,
n_pix=self.n_pix,
weight_name=weight_name,
pixel_match_scale=self.pixel_match_scale,
)

def _buildClusterData(
Expand All @@ -250,7 +251,7 @@ def _buildClusterData(
return ClusterData(i_cluster, footprint, sources, pixel_match_scale=self.pixel_match_scale)

def _getFootprints(self, counts_map: np.ndarray) -> dict:
return utils.getFootprints(counts_map, buf=self.buf)
return utils.getFootprints(counts_map, buf=self.buf, pixel_match_scale=self.pixel_match_scale)

def _associateSourcesToFootprints(
self,
Expand All @@ -260,6 +261,7 @@ def _associateSourcesToFootprints(
return utils.associateSourcesToFootprints(
data,
cluster_key,
pixel_match_scale=self.pixel_match_scale,
)

def getRaDec(
Expand Down Expand Up @@ -308,21 +310,6 @@ def _newObject(
) -> ObjectData:
return ShearObjectData(cluster, object_id, mask)

def _emptyCountsMaps(self) -> np.ndarray:
pixel_match_scale = self.pixel_match_scale
to_fill = np.zeros(np.ceil(self.n_pix / pixel_match_scale).astype(int))
return to_fill

def _singleCatalogCountsMap(
self, df: pandas.DataFrame, weight_name: str | None = None
) -> np.ndarray:
return utils.fillCountsMapFromDf(
df,
n_pix=self.n_pix,
weight_name=weight_name,
pixel_match_scale=self.pixel_match_scale,
)

def _buildClusterData(
self, i_cluster: int, footprint: Footprint, sources: np.ndarray
) -> ClusterData:
Expand All @@ -339,14 +326,3 @@ def _getFootprints(self, counts_map: np.ndarray) -> dict:
counts_map, buf=0, pixel_match_scale=self.pixel_match_scale
)

def _associateSourcesToFootprints(
self,
data: list[pandas.DataFrame],
cluster_key: np.ndarray,
) -> list[np.ndarray]:
return utils.associateSourcesToFootprints(
data,
cluster_key,
pixel_match_scale=self.pixel_match_scale,
)

81 changes: 62 additions & 19 deletions src/hpmcm/match.py
Original file line number Diff line number Diff line change
Expand Up @@ -291,9 +291,40 @@ def analysisLoop(
sys.stdout.write(" Done!\n")
sys.stdout.flush()

def extractStats(self) -> dict[str, pandas.DataFrame]:
def _getCentralIds(
self, cell_data: CellData
) -> tuple[frozenset[int], frozenset[int]]:
"""Return (object_ids, cluster_ids) whose centroids lie in the inner cell region.

The inner region is ``[0, cell_data.size[i])`` in x_cell_coadd coordinates,
which corresponds to ``[cell_data.buf, cell_data.buf + cell_data.size[i])``
in the internal x_cell convention stored on each object/cluster.
"""
lo_x, hi_x = float(cell_data.buf), float(cell_data.buf + cell_data.size[0])
lo_y, hi_y = float(cell_data.buf), float(cell_data.buf + cell_data.size[1])
obj_ids = frozenset(
oid
for oid, obj in cell_data.object_dict.items()
if lo_x <= obj.x_cent < hi_x and lo_y <= obj.y_cent < hi_y
)
clust_ids = frozenset(
cid
for cid, cl in cell_data.cluster_dict.items()
if lo_x <= cl.x_cent < hi_x and lo_y <= cl.y_cent < hi_y
)
return obj_ids, clust_ids

def extractStats(self, central_only: bool = True) -> dict[str, pandas.DataFrame]:
"""Extracts cluster statisistics

Parameters
----------
central_only:
When ``True`` (default), exclude objects and clusters whose centroid
(in x_cell_coadd coordinates) falls outside the inner cell region
``[0, cell_size)``. This avoids double-counting edge objects that
are also captured by adjacent cells.

Returns
-------
Dict with keys:
Expand All @@ -314,18 +345,23 @@ def extractStats(self) -> dict[str, pandas.DataFrame]:
if i_cell not in self.cell_dict:
continue
cell_data = self.cell_dict[i_cell]
cluster_assoc_tables.append(
output_tables.ClusterAssocTable.buildFromCellData(cell_data).data,
)
object_assoc_tables.append(
output_tables.ObjectAssocTable.buildFromCellData(cell_data).data,
)
cluster_stats_tables.append(
output_tables.ClusterStatsTable.buildFromCellData(cell_data).data,
)
object_stats_tables.append(
output_tables.ObjectStatsTable.buildFromCellData(cell_data).data,
)

ca = output_tables.ClusterAssocTable.buildFromCellData(cell_data).data
oa = output_tables.ObjectAssocTable.buildFromCellData(cell_data).data
cs = output_tables.ClusterStatsTable.buildFromCellData(cell_data).data
os_ = output_tables.ObjectStatsTable.buildFromCellData(cell_data).data

if central_only:
c_obj, c_clust = self._getCentralIds(cell_data)
ca = ca[ca["cluster_id"].isin(c_clust)]
oa = oa[oa["object_id"].isin(c_obj)]
cs = cs[cs["cluster_id"].isin(c_clust)]
os_ = os_[os_["object_id"].isin(c_obj)]

cluster_assoc_tables.append(ca)
object_assoc_tables.append(oa)
cluster_stats_tables.append(cs)
object_stats_tables.append(os_)
if ix == 0:
pass
elif ix % 10 == 0:
Expand Down Expand Up @@ -370,7 +406,7 @@ def _geometryDict(self) -> dict[str, Any]:
"catalog_id_map": {str(k): int(v) for k, v in self.catalog_id_map.items()},
}

def save(self, save_dir: str | Path) -> None:
def save(self, save_dir: str | Path, central_only: bool = True) -> None:
"""Persist the full match state to a directory.

Saves geometry, association/stats tables, and one parquet per catalog
Expand All @@ -381,14 +417,18 @@ def save(self, save_dir: str | Path) -> None:
----------
save_dir:
Directory to write into (created if absent).
central_only:
Passed through to :meth:`extractStats`; when ``True`` (default),
only objects and clusters whose centroid lies in the inner cell
region are written.
"""
save_dir = Path(save_dir)
save_dir.mkdir(parents=True, exist_ok=True)

with open(save_dir / "geometry.json", "w") as fh:
json.dump(_to_json_safe(self._geometryDict()), fh, indent=2)

stats = self.extractStats()
stats = self.extractStats(central_only=central_only)
for name, df in stats.items():
df.to_parquet(save_dir / f"{name}.parquet")

Expand Down Expand Up @@ -560,8 +600,10 @@ def _reconstructCells(
cluster = cell._buildClusterData(cluster_id, fp, sources)
cluster.extract(cell)

cluster.x_cent = float(sr["x_cent"])
cluster.y_cent = float(sr["y_cent"])
# x_cent saved in x_cell_coadd convention (0 at inner edge);
# restore to x_cell convention (0 at outer edge) for internal use
cluster.x_cent = float(sr["x_cent"]) + cell.buf
cluster.y_cent = float(sr["y_cent"]) + cell.buf
cluster.dist_2 = (ca_rows["distance"].values / pix_to_arcsec) ** 2
cluster.rms_dist = float(sr["dist_rms"]) / pix_to_arcsec
cluster.snr_mean = float(sr["snr"])
Expand All @@ -575,8 +617,9 @@ def _reconstructCells(
obj = cell._newObject(cluster, obj_id, mask)
or_ = os_idx.loc[obj_id]
obj.dist_2 = (obj_rows["distance"].values / pix_to_arcsec) ** 2
obj.x_cent = float(or_["x_cent"])
obj.y_cent = float(or_["y_cent"])
# x_cent saved in x_cell_coadd convention; restore to x_cell
obj.x_cent = float(or_["x_cent"]) + cell.buf
obj.y_cent = float(or_["y_cent"]) + cell.buf
obj.rms_dist = float(or_["dist_rms"]) / pix_to_arcsec
obj.snr_mean = float(or_["snr"])
obj.snr_rms = float(or_["snr_rms"])
Expand Down
64 changes: 55 additions & 9 deletions src/hpmcm/output_tables.py
Original file line number Diff line number Diff line change
Expand Up @@ -98,8 +98,8 @@ class ObjectStatsTable(TableInterface):
),
ra=TableColumnInfo(float, "RA of object centroid"),
dec=TableColumnInfo(float, "DEC of object centroid"),
x_cent=TableColumnInfo(float, "X-value of object centroid in cell pixels"),
y_cent=TableColumnInfo(float, "Y-value of object centroid in cell pixels"),
x_cent=TableColumnInfo(float, "X-value of object centroid in x_cell_coadd coords (0 at inner left edge)"),
y_cent=TableColumnInfo(float, "Y-value of object centroid in x_cell_coadd coords (0 at inner left edge)"),
x_pix=TableColumnInfo(float, "X-value of object centroid in global WCS pixels"),
y_pix=TableColumnInfo(float, "Y-value of object centroid in global WCS pixels"),
snr=TableColumnInfo(float, "Mean signal-to-noise ratio"),
Expand Down Expand Up @@ -163,6 +163,9 @@ def buildFromCellData(cell_data: CellData) -> ObjectStatsTable:
dist_rms *= cell_data.matcher.pixToArcsec()
x_pix = x_cents + cell_data.min_pix[0]
y_pix = y_cents + cell_data.min_pix[1]
# Convert centroids from x_cell (0 at outer edge) to x_cell_coadd (0 at inner edge)
x_cents_coadd = x_cents - cell_data.buf
y_cents_coadd = y_cents - cell_data.buf

return ObjectStatsTable(
cluster_id=cluster_ids,
Expand All @@ -172,8 +175,8 @@ def buildFromCellData(cell_data: CellData) -> ObjectStatsTable:
dist_rms=dist_rms,
ra=ra,
dec=dec,
x_cent=x_cents,
y_cent=y_cents,
x_cent=x_cents_coadd,
y_cent=y_cents_coadd,
x_pix=x_pix,
y_pix=y_pix,
snr=snrs,
Expand Down Expand Up @@ -260,8 +263,8 @@ class ClusterStatsTable(TableInterface):
),
ra=TableColumnInfo(float, "RA of cluster centroid"),
dec=TableColumnInfo(float, "DEC of cluster centroid"),
x_cent=TableColumnInfo(float, "X-value of cluster centroid in cell pixels"),
y_cent=TableColumnInfo(float, "Y-value of cluster centroid in cell pixels"),
x_cent=TableColumnInfo(float, "X-value of cluster centroid in x_cell_coadd coords (0 at inner left edge)"),
y_cent=TableColumnInfo(float, "Y-value of cluster centroid in x_cell_coadd coords (0 at inner left edge)"),
x_pix=TableColumnInfo(float, "X-value of cluster centroid in global WCS pixels"),
y_pix=TableColumnInfo(float, "Y-value of cluster centroid in global WCS pixels"),
snr=TableColumnInfo(float, "Mean signal-to-noise ratio"),
Expand Down Expand Up @@ -327,6 +330,9 @@ def buildFromCellData(cell_data: CellData) -> ClusterStatsTable:
dist_rms *= cell_data.matcher.pixToArcsec()
x_pix = x_cents + cell_data.min_pix[0]
y_pix = y_cents + cell_data.min_pix[1]
# Convert centroids from x_cell (0 at outer edge) to x_cell_coadd (0 at inner edge)
x_cents_coadd = x_cents - cell_data.buf
y_cents_coadd = y_cents - cell_data.buf

return ClusterStatsTable(
cluster_id=cluster_ids,
Expand All @@ -336,8 +342,8 @@ def buildFromCellData(cell_data: CellData) -> ClusterStatsTable:
dist_rms=dist_rms,
ra=ra,
dec=dec,
x_cent=x_cents,
y_cent=y_cents,
x_cent=x_cents_coadd,
y_cent=y_cents_coadd,
x_pix=x_pix,
y_pix=y_pix,
snr=snrs,
Expand Down Expand Up @@ -464,6 +470,8 @@ def buildJoinedObjectTable(
catalog_ids: list[int],
source_cols: SourceColsType = None,
object_shear: str | Path | None = None,
central_only: bool = True,
cell_size: float | None = None,
) -> pandas.DataFrame:
"""Build a wide joined table from an ObjectStatsTable and source catalogs.

Expand Down Expand Up @@ -500,13 +508,33 @@ def buildJoinedObjectTable(
shear statistics are joined on ``object_id``. The ``cluster_id``
column is dropped before joining because it is already present in
``object_stats``.
central_only:
When ``True`` (default), rows whose centroid (``x_cent``, ``y_cent``
in x_cell_coadd coordinates) has a negative x or y component are
dropped before joining. When ``cell_size`` is also given, objects
whose centroid is >= ``cell_size`` in either axis are also removed.
This mirrors the filtering applied by
:meth:`~hpmcm.match.Match.extractStats` and is a no-op when the
input already comes from that method with ``central_only=True``.
cell_size:
Inner cell size in pixels. When provided together with
``central_only=True``, objects with ``x_cent >= cell_size`` or
``y_cent >= cell_size`` are also removed. Pass ``matcher.cell_size``
when calling with unfiltered stats.

Returns
-------
DataFrame with one row per ``object_id``. Columns from each catalog are
renamed ``{col}_{catalog_id}``. Objects that have no source in a
given catalog will have ``NaN`` for that catalog's columns.
"""
if central_only:
mask = (object_stats["x_cent"] >= 0) & (object_stats["y_cent"] >= 0)
if cell_size is not None:
mask &= (object_stats["x_cent"] < cell_size) & (object_stats["y_cent"] < cell_size)
object_stats = object_stats[mask]
object_assoc = object_assoc[object_assoc["object_id"].isin(object_stats["object_id"])]

catalog_file_map = dict(zip(catalog_ids, input_files))
base = object_stats.set_index("object_id")

Expand Down Expand Up @@ -544,13 +572,14 @@ def buildJoinedObjectTable(


def buildJoinedClusterTable(

cluster_stats: pandas.DataFrame,
cluster_assoc: pandas.DataFrame,
input_files: list[str],
catalog_ids: list[int],
source_cols: SourceColsType = None,
cluster_shear: str | Path | None = None,
central_only: bool = True,
cell_size: float | None = None,
) -> pandas.DataFrame:
"""Build a wide joined table from a ClusterStatsTable and source catalogs.

Expand Down Expand Up @@ -582,13 +611,30 @@ def buildJoinedClusterTable(
cluster_shear:
Optional path to a ``ClusterShearTable`` parquet file. When provided,
shear statistics are joined on ``cluster_id``.
central_only:
When ``True`` (default), clusters whose centroid (``x_cent``,
``y_cent`` in x_cell_coadd coordinates) has a negative component are
dropped before joining. When ``cell_size`` is also given, clusters
with centroid >= ``cell_size`` are also removed.
cell_size:
Inner cell size in pixels. When provided together with
``central_only=True``, clusters with ``x_cent >= cell_size`` or
``y_cent >= cell_size`` are also removed. Pass ``matcher.cell_size``
when calling with unfiltered stats.

Returns
-------
DataFrame with one row per ``cluster_id``. Columns from each catalog are
renamed ``{col}_{catalog_id}``. Clusters that have no source in a
given catalog will have ``NaN`` for that catalog's columns.
"""
if central_only:
mask = (cluster_stats["x_cent"] >= 0) & (cluster_stats["y_cent"] >= 0)
if cell_size is not None:
mask &= (cluster_stats["x_cent"] < cell_size) & (cluster_stats["y_cent"] < cell_size)
cluster_stats = cluster_stats[mask]
cluster_assoc = cluster_assoc[cluster_assoc["cluster_id"].isin(cluster_stats["cluster_id"])]

catalog_file_map = dict(zip(catalog_ids, input_files))
base = cluster_stats.set_index("cluster_id")

Expand Down
26 changes: 19 additions & 7 deletions src/hpmcm/shear_match.py
Original file line number Diff line number Diff line change
Expand Up @@ -270,9 +270,16 @@ def _buildCellData(
) -> CellData:
return ShearCellData(self, id_offset, corner, size, idx, self.cell_buffer)

def extractShearStats(self) -> dict[str, pandas.DataFrame]:
def extractShearStats(self, central_only: bool = True) -> dict[str, pandas.DataFrame]:
"""Extract shear stats

Parameters
----------
central_only:
When ``True`` (default), exclude objects and clusters whose centroid
falls outside the inner cell region. Mirrors the behaviour of
:meth:`~hpmcm.match.Match.extractStats`.

Returns
-------
Dict with keys:
Expand All @@ -289,12 +296,17 @@ def extractShearStats(self) -> dict[str, pandas.DataFrame]:
continue
cell_data = self.cell_dict[i_cell]
assert isinstance(cell_data, ShearCellData)
cluster_shear_stats_tables.append(
output_tables.ShearTable.buildClusterShearStats(cell_data).data
)
object_shear_stats_tables.append(
output_tables.ShearTable.buildObjectShearStats(cell_data).data
)

cs = output_tables.ShearTable.buildClusterShearStats(cell_data).data
os_ = output_tables.ShearTable.buildObjectShearStats(cell_data).data

if central_only:
c_obj, c_clust = self._getCentralIds(cell_data)
cs = cs[cs["cluster_id"].isin(c_clust)]
os_ = os_[os_["object_id"].isin(c_obj)]

cluster_shear_stats_tables.append(cs)
object_shear_stats_tables.append(os_)

return {
"cluster_shear": pandas.concat(cluster_shear_stats_tables),
Expand Down
Loading
Loading