Adding support for TwoPointMeasurement filters. (#479)

* Adding support for TwoPointMeasurement filters.
* Adding tests for filters.
* Removed unused import.
* Using two-point pair as filter specification.
* Adding support for serializing and deserializing filters.
* Adding tests for factories.
* Added check for never reached branch.
* TwoPointFilter documentation first draft
* Simplify logic
* If _path is set do not search current directory
* Release cython version restriction
* Do not load our duplicate_code plugin
* Update finding of SACC files in some tests
* Update version tag
* Refactor for improved test coverage and fix missing error case
* Complete branch coverage
* Improving tutorial.
* Correct serialization for TwoPointFactory.
* Added test for TwoPointFactory serialization.

---------

Co-authored-by: paulrogozenski <[email protected]>
Co-authored-by: Marc Paterno <[email protected]>

Author: 3 people

docs/conf.py

@@ -23,7 +23,7 @@
 author = "LSST DESC Firecrown Contributors"
 
 # The full version, including alpha/beta/rc tags
-release = "1.8.0"
+release = "1.9.0a0"
 
 
 # -- General configuration ---------------------------------------------------

environment.yml

@@ -9,7 +9,7 @@ dependencies:
   - cosmosis >= 3.0
   - cosmosis-build-standard-library
   - coverage
-  - cython < 3.0.0
+  - cython
   - dill
   - fitsio
   - flake8

fctools/tracer.py

@@ -24,7 +24,7 @@
 
 # some global context to be used in the tracing. We are relying on
 # 'trace_call' to act as a closure that captures these names.
-tracefile = None  # the file used for logging
+tracefile: TextIO | None = None  # the file used for logging
 level = 0  # the call nesting level
 entry = 0  # sequential entry number for each record
 

firecrown/data_functions.py

@@ -4,20 +4,35 @@
 """
 
 import hashlib
-from typing import Callable, Sequence
+from typing import Callable, Sequence, Annotated
+from typing_extensions import assert_never
+
+import sacc
+from pydantic import (
+    BaseModel,
+    BeforeValidator,
+    ConfigDict,
+    Field,
+    model_validator,
+    PrivateAttr,
+    field_serializer,
+)
 import numpy as np
 import numpy.typing as npt
-import sacc
+
 from firecrown.metadata_types import (
     TwoPointHarmonic,
     TwoPointReal,
+    Measurement,
 )
 from firecrown.metadata_functions import (
     extract_all_tracers_inferred_galaxy_zdists,
     extract_window_function,
     extract_all_harmonic_metadata_indices,
     extract_all_real_metadata_indices,
     make_two_point_xy,
+    make_measurement,
+    make_measurement_dict,
 )
 from firecrown.data_types import TwoPointMeasurement
 
@@ -222,3 +237,275 @@ def check_two_point_consistence_real(
 ) -> None:
     """Check the indices of the real-space two-point functions."""
     check_consistence(two_point_reals, lambda m: m.is_real(), "TwoPointReal")
+
+
+class TwoPointTracerSpec(BaseModel):
+    """Class defining a tracer bin specification."""
+
+    model_config = ConfigDict(extra="forbid", frozen=True)
+
+    name: Annotated[str, Field(description="The name of the tracer bin.")]
+    measurement: Annotated[
+        Measurement,
+        Field(description="The measurement of the tracer bin."),
+        BeforeValidator(make_measurement),
+    ]
+
+    @field_serializer("measurement")
+    @classmethod
+    def serialize_measurement(cls, value: Measurement) -> dict[str, str]:
+        """Serialize the Measurement."""
+        return make_measurement_dict(value)
+
+
+def make_interval_from_list(
+    values: list[float] | tuple[float, float],
+) -> tuple[float, float]:
+    """Create an interval from a list of values."""
+    if isinstance(values, list):
+        if len(values) != 2:
+            raise ValueError("The list should have two values.")
+        if not all(isinstance(v, float) for v in values):
+            raise ValueError("The list should have two float values.")
+
+        return (values[0], values[1])
+    if isinstance(values, tuple):
+        return values
+
+    raise ValueError("The values should be a list or a tuple.")
+
+
+class TwoPointBinFilter(BaseModel):
+    """Class defining a filter for a bin."""
+
+    model_config = ConfigDict(extra="forbid", frozen=True)
+
+    spec: Annotated[
+        list[TwoPointTracerSpec],
+        Field(
+            description="The two-point bin specification.",
+        ),
+    ]
+    interval: Annotated[
+        tuple[float, float],
+        BeforeValidator(make_interval_from_list),
+        Field(description="The range of the bin to filter."),
+    ]
+
+    @model_validator(mode="after")
+    def check_bin_filter(self) -> "TwoPointBinFilter":
+        """Check the bin filter."""
+        if self.interval[0] >= self.interval[1]:
+            raise ValueError("The bin filter should be a valid range.")
+        if not 1 <= len(self.spec) <= 2:
+            raise ValueError("The bin_spec must contain one or two elements.")
+        return self
+
+    @field_serializer("interval")
+    @classmethod
+    def serialize_interval(cls, value: tuple[float, float]) -> list[float]:
+        """Serialize the Measurement."""
+        return list(value)
+
+    @classmethod
+    def from_args(
+        cls,
+        name1: str,
+        measurement1: Measurement,
+        name2: str,
+        measurement2: Measurement,
+        lower: float,
+        upper: float,
+    ) -> "TwoPointBinFilter":
+        """Create a TwoPointBinFilter from the arguments."""
+        return cls(
+            spec=[
+                TwoPointTracerSpec(name=name1, measurement=measurement1),
+                TwoPointTracerSpec(name=name2, measurement=measurement2),
+            ],
+            interval=(lower, upper),
+        )
+
+    @classmethod
+    def from_args_auto(
+        cls, name: str, measurement: Measurement, lower: float, upper: float
+    ) -> "TwoPointBinFilter":
+        """Create a TwoPointBinFilter from the arguments."""
+        return cls(
+            spec=[
+                TwoPointTracerSpec(name=name, measurement=measurement),
+            ],
+            interval=(lower, upper),
+        )
+
+
+BinSpec = frozenset[TwoPointTracerSpec]
+
+
+def bin_spec_from_metadata(metadata: TwoPointReal | TwoPointHarmonic) -> BinSpec:
+    """Return the bin spec from the metadata."""
+    return frozenset(
+        (
+            TwoPointTracerSpec(
+                name=metadata.XY.x.bin_name,
+                measurement=metadata.XY.x_measurement,
+            ),
+            TwoPointTracerSpec(
+                name=metadata.XY.y.bin_name,
+                measurement=metadata.XY.y_measurement,
+            ),
+        )
+    )
+
+
+class TwoPointBinFilterCollection(BaseModel):
+    """Class defining a collection of bin filters."""
+
+    model_config = ConfigDict(extra="forbid", frozen=True)
+
+    require_filter_for_all: bool = Field(
+        default=False,
+        description="If True, all bins should match a filter.",
+    )
+    allow_empty: bool = Field(
+        default=False,
+        description=(
+            "When true, objects with no elements remaining after applying "
+            "the filter will be ignored rather than treated as an error."
+        ),
+    )
+    filters: list[TwoPointBinFilter] = Field(
+        description="The list of bin filters.",
+    )
+
+    _bin_filter_dict: dict[BinSpec, tuple[float, float]] = PrivateAttr()
+
+    @model_validator(mode="after")
+    def check_bin_filters(self) -> "TwoPointBinFilterCollection":
+        """Check the bin filters."""
+        bin_specs = set()
+        for bin_filter in self.filters:
+            bin_spec = frozenset(bin_filter.spec)
+            if bin_spec in bin_specs:
+                raise ValueError(
+                    f"The bin name {bin_filter.spec} is repeated "
+                    f"in the bin filters."
+                )
+            bin_specs.add(bin_spec)
+
+        self._bin_filter_dict = {
+            frozenset(bin_filter.spec): bin_filter.interval
+            for bin_filter in self.filters
+        }
+        return self
+
+    @property
+    def bin_filter_dict(self) -> dict[BinSpec, tuple[float, float]]:
+        """Return the bin filter dictionary."""
+        return self._bin_filter_dict
+
+    def filter_match(self, tpm: TwoPointMeasurement) -> bool:
+        """Check if the TwoPointMeasurement matches the filter."""
+        bin_spec_key = bin_spec_from_metadata(tpm.metadata)
+        return bin_spec_key in self._bin_filter_dict
+
+    def run_bin_filter(
+        self,
+        bin_filter: tuple[float, float],
+        vals: npt.NDArray[np.float64] | npt.NDArray[np.int64],
+    ) -> npt.NDArray[np.bool_]:
+        """Run the filter merge."""
+        return (vals >= bin_filter[0]) & (vals <= bin_filter[1])
+
+    def apply_filter_single(
+        self, tpm: TwoPointMeasurement
+    ) -> tuple[npt.NDArray[np.bool_], npt.NDArray[np.bool_]]:
+        """Apply the filter to a single TwoPointMeasurement."""
+        assert self.filter_match(tpm)
+        bin_spec_key = bin_spec_from_metadata(tpm.metadata)
+        bin_filter = self._bin_filter_dict[bin_spec_key]
+        if tpm.is_real():
+            assert isinstance(tpm.metadata, TwoPointReal)
+            match_elements = self.run_bin_filter(bin_filter, tpm.metadata.thetas)
+            return match_elements, match_elements
+
+        assert isinstance(tpm.metadata, TwoPointHarmonic)
+        match_elements = self.run_bin_filter(bin_filter, tpm.metadata.ells)
+        match_obs = match_elements
+        if tpm.metadata.window is not None:
+            # The window function is represented by a matrix where each column
+            # corresponds to the weights for the ell values of each observation. We
+            # need to ensure that the window function is filtered correctly. To do this,
+            # we will check each column of the matrix and verify that all non-zero
+            # elements are within the filtered set. If any non-zero element falls
+            # outside the filtered set, the match_elements will be set to False for that
+            # observation.
+            non_zero_window = tpm.metadata.window > 0
+            match_obs = (
+                np.all(
+                    (non_zero_window & match_elements[:, None]) == non_zero_window,
+                    axis=0,
+                )
+                .ravel()
+                .astype(np.bool_)
+            )
+
+        return match_elements, match_obs
+
+    def __call__(
+        self, tpms: Sequence[TwoPointMeasurement]
+    ) -> list[TwoPointMeasurement]:
+        """Filter the two-point measurements."""
+        result = []
+
+        for tpm in tpms:
+            if not self.filter_match(tpm):
+                if not self.require_filter_for_all:
+                    result.append(tpm)
+                    continue
+                raise ValueError(f"The bin name {tpm.metadata} does not have a filter.")
+
+            match_elements, match_obs = self.apply_filter_single(tpm)
+            if not match_obs.any():
+                if not self.allow_empty:
+                    # If empty results are not allowed, we raise an error
+                    raise ValueError(
+                        f"The TwoPointMeasurement {tpm.metadata} does not "
+                        f"have any elements matching the filter."
+                    )
+                # If the filter is empty, we skip this measurement
+                continue
+
+            assert isinstance(tpm.metadata, (TwoPointReal, TwoPointHarmonic))
+            new_metadata: TwoPointReal | TwoPointHarmonic
+            match tpm.metadata:
+                case TwoPointReal():
+                    new_metadata = TwoPointReal(
+                        XY=tpm.metadata.XY,
+                        thetas=tpm.metadata.thetas[match_elements],
+                    )
+                case TwoPointHarmonic():
+                    # If the window function is not None, we need to filter it as well
+                    # and update the metadata accordingly.
+                    new_metadata = TwoPointHarmonic(
+                        XY=tpm.metadata.XY,
+                        window=(
+                            tpm.metadata.window[:, match_obs][match_elements, :]
+                            if tpm.metadata.window is not None
+                            else None
+                        ),
+                        ells=tpm.metadata.ells[match_elements],
+                    )
+                case _ as unreachable:
+                    assert_never(unreachable)
+
+            result.append(
+                TwoPointMeasurement(
+                    data=tpm.data[match_obs],
+                    indices=tpm.indices[match_obs],
+                    covariance_name=tpm.covariance_name,
+                    metadata=new_metadata,
+                )
+            )
+
+        return result

firecrown/generators/inferred_galaxy_zdist.py

@@ -13,15 +13,12 @@
 
 from numcosmo_py import Ncm
 
-from firecrown.metadata_types import (
-    InferredGalaxyZDist,
-    ALL_MEASUREMENT_TYPES,
+from firecrown.metadata_types import InferredGalaxyZDist, Galaxies
+from firecrown.metadata_functions import (
+    Measurement,
+    make_measurements,
     make_measurements_dict,
-    Galaxies,
-    CMB,
-    Clusters,
 )
-from firecrown.metadata_functions import Measurement
 
 
 BinsType = TypedDict("BinsType", {"edges": npt.NDArray, "sigma_z": float})
@@ -446,41 +443,6 @@ def generate(self) -> npt.NDArray:
 Grid1D = LinearGrid1D | RawGrid1D
 
 
-def make_measurements(
-    value: set[Measurement] | list[dict[str, Any]],
-) -> set[Measurement]:
-    """Create a Measurement object from a dictionary."""
-    if isinstance(value, set) and all(
-        isinstance(v, ALL_MEASUREMENT_TYPES) for v in value
-    ):
-        return value
-
-    measurements: set[Measurement] = set()
-    for measurement_dict in value:
-        if not isinstance(measurement_dict, dict):
-            raise ValueError(f"Invalid Measurement: {value} is not a dictionary")
-
-        if "subject" not in measurement_dict:
-            raise ValueError(
-                "Invalid Measurement: dictionary does not contain 'subject'"
-            )
-
-        subject = measurement_dict["subject"]
-
-        match subject:
-            case "Galaxies":
-                measurements.update({Galaxies[measurement_dict["property"]]})
-            case "CMB":
-                measurements.update({CMB[measurement_dict["property"]]})
-            case "Clusters":
-                measurements.update({Clusters[measurement_dict["property"]]})
-            case _:
-                raise ValueError(
-                    f"Invalid Measurement: subject: '{subject}' is not recognized"
-                )
-    return measurements
-
-
 class ZDistLSSTSRDBin(BaseModel):
     """LSST Inferred galaxy redshift distributions in bins."""
 

firecrown/likelihood/factories.py

@@ -21,7 +21,7 @@
 
 from typing_extensions import assert_never
 import yaml
-from pydantic import BaseModel, ConfigDict, BeforeValidator
+from pydantic import BaseModel, ConfigDict, BeforeValidator, Field, field_serializer
 
 import sacc
 from firecrown.likelihood.likelihood import Likelihood, NamedParameters
@@ -34,6 +34,7 @@
     extract_all_harmonic_data,
     check_two_point_consistence_real,
     check_two_point_consistence_harmonic,
+    TwoPointBinFilterCollection,
 )
 from firecrown.modeling_tools import ModelingTools
 from firecrown.ccl_factory import CCLFactory
@@ -56,10 +57,12 @@ def _generate_next_value_(name, _start, _count, _last_values):
     HARMONIC = auto()
 
 
-def _validate_correlation_space(value):
-    if isinstance(value, str):
+def _validate_correlation_space(value: TwoPointCorrelationSpace | str):
+    if not isinstance(value, TwoPointCorrelationSpace) and isinstance(value, str):
         try:
-            return TwoPointCorrelationSpace(value)  # Convert from string to Enum
+            return TwoPointCorrelationSpace(
+                value.lower()
+            )  # Convert from string to Enum
         except ValueError as exc:
             raise ValueError(
                 f"Invalid value for TwoPointCorrelationSpace: {value}"
@@ -73,42 +76,60 @@ class TwoPointFactory(BaseModel):
     model_config = ConfigDict(extra="forbid", frozen=True)
 
     correlation_space: Annotated[
-        TwoPointCorrelationSpace, BeforeValidator(_validate_correlation_space)
+        TwoPointCorrelationSpace,
+        BeforeValidator(_validate_correlation_space),
+        Field(description="The two-point correlation space."),
     ]
     weak_lensing_factory: WeakLensingFactory
     number_counts_factory: NumberCountsFactory
 
     def model_post_init(self, __context) -> None:
         """Initialize the WeakLensingFactory object."""
 
+    @field_serializer("correlation_space")
+    @classmethod
+    def serialize_correlation_space(cls, value: TwoPointCorrelationSpace) -> str:
+        """Serialize the amplitude parameter."""
+        return value.name
+
 
 class DataSourceSacc(BaseModel):
     """Model for the data source in a likelihood configuration."""
 
     sacc_data_file: str
+    filters: TwoPointBinFilterCollection | None = None
     _path: Path | None = None
 
     def set_path(self, path: Path) -> None:
         """Set the path for the data source."""
         self._path = path
 
-    def get_sacc_data(self) -> sacc.Sacc:
-        """Load the SACC data file."""
+    def get_filepath(self) -> Path:
+        """Return the filename of the data source.
+
+        Raises a FileNotFoundError if the file does not exist.
+        :return: The filename
+        """
         sacc_data_path = Path(self.sacc_data_file)
         # If sacc_data_file is absolute, use it directly
-        if sacc_data_path.is_absolute():
-            return sacc.Sacc.load_fits(self.sacc_data_file)
+        if sacc_data_path.is_absolute() and sacc_data_path.exists():
+            return Path(self.sacc_data_file)
         # If path is set, use it to find the file
         if self._path is not None:
             full_sacc_data_path = self._path / sacc_data_path
             if full_sacc_data_path.exists():
-                return sacc.Sacc.load_fits(full_sacc_data_path)
+                return full_sacc_data_path
         # If path is not set, use the current directory
-        if sacc_data_path.exists():
-            return sacc.Sacc.load_fits(sacc_data_path)
+        elif sacc_data_path.exists():
+            return sacc_data_path
         # If the file does not exist, raise an error
         raise FileNotFoundError(f"File {sacc_data_path} does not exist")
 
+    def get_sacc_data(self) -> sacc.Sacc:
+        """Load the SACC data file."""
+        filename = self.get_filepath()
+        return sacc.Sacc.load_fits(filename)
+
 
 def ensure_path(file: str | Path) -> Path:
     """Ensure the file path is a Path object."""
@@ -130,6 +151,8 @@ class TwoPointExperiment(BaseModel):
 
     def model_post_init(self, __context) -> None:
         """Initialize the TwoPointExperiment object."""
+        if self.ccl_factory is None:
+            self.ccl_factory = CCLFactory()
 
     @classmethod
     def load_from_yaml(cls, file: str | Path) -> "TwoPointExperiment":
@@ -144,6 +167,32 @@ def load_from_yaml(cls, file: str | Path) -> "TwoPointExperiment":
         tpe.data_source.set_path(filepath.parent)
         return tpe
 
+    def make_likelihood(self) -> Likelihood:
+        """Create a likelihood object for two-point statistics from a SACC file."""
+        # Load the SACC file
+        sacc_data = self.data_source.get_sacc_data()
+
+        likelihood: None | Likelihood = None
+        match self.two_point_factory.correlation_space:
+            case TwoPointCorrelationSpace.REAL:
+                likelihood = _build_two_point_likelihood_real(
+                    sacc_data,
+                    self.two_point_factory.weak_lensing_factory,
+                    self.two_point_factory.number_counts_factory,
+                    filters=self.data_source.filters,
+                )
+            case TwoPointCorrelationSpace.HARMONIC:
+                likelihood = _build_two_point_likelihood_harmonic(
+                    sacc_data,
+                    self.two_point_factory.weak_lensing_factory,
+                    self.two_point_factory.number_counts_factory,
+                    filters=self.data_source.filters,
+                )
+            case _ as unreachable:
+                assert_never(unreachable)
+        assert likelihood is not None
+        return likelihood
+
 
 def build_two_point_likelihood(
     build_parameters: NamedParameters,
@@ -165,24 +214,7 @@ def build_two_point_likelihood(
     exp = TwoPointExperiment.load_from_yaml(likelihood_config_file)
     modeling_tools = ModelingTools(ccl_factory=exp.ccl_factory)
 
-    # Load the SACC file
-    sacc_data = exp.data_source.get_sacc_data()
-
-    match exp.two_point_factory.correlation_space:
-        case TwoPointCorrelationSpace.REAL:
-            likelihood = _build_two_point_likelihood_real(
-                sacc_data,
-                exp.two_point_factory.weak_lensing_factory,
-                exp.two_point_factory.number_counts_factory,
-            )
-        case TwoPointCorrelationSpace.HARMONIC:
-            likelihood = _build_two_point_likelihood_harmonic(
-                sacc_data,
-                exp.two_point_factory.weak_lensing_factory,
-                exp.two_point_factory.number_counts_factory,
-            )
-        case _ as unreachable:
-            assert_never(unreachable)
+    likelihood = exp.make_likelihood()
 
     return likelihood, modeling_tools
 
@@ -191,6 +223,7 @@ def _build_two_point_likelihood_harmonic(
     sacc_data: sacc.Sacc,
     wl_factory: WeakLensingFactory,
     nc_factory: NumberCountsFactory,
+    filters: TwoPointBinFilterCollection | None = None,
 ):
     """
     Build a likelihood object for two-point statistics in harmonic space.
@@ -211,8 +244,9 @@ def _build_two_point_likelihood_harmonic(
         raise ValueError(
             "No two-point measurements in harmonic space found in the SACC file."
         )
-
     check_two_point_consistence_harmonic(tpms)
+    if filters is not None:
+        tpms = filters(tpms)
 
     two_points = TwoPoint.from_measurement(
         tpms, wl_factory=wl_factory, nc_factory=nc_factory
@@ -227,6 +261,7 @@ def _build_two_point_likelihood_real(
     sacc_data: sacc.Sacc,
     wl_factory: WeakLensingFactory,
     nc_factory: NumberCountsFactory,
+    filters: TwoPointBinFilterCollection | None = None,
 ):
     """
     Build a likelihood object for two-point statistics in real space.
@@ -248,6 +283,8 @@ def _build_two_point_likelihood_real(
             "No two-point measurements in real space found in the SACC file."
         )
     check_two_point_consistence_real(tpms)
+    if filters is not None:
+        tpms = filters(tpms)
 
     two_points = TwoPoint.from_measurement(
         tpms, wl_factory=wl_factory, nc_factory=nc_factory

firecrown/metadata_functions.py

@@ -5,7 +5,7 @@
 """
 
 from itertools import combinations_with_replacement, product
-from typing import TypedDict
+from typing import TypedDict, Any
 
 import numpy as np
 import numpy.typing as npt
@@ -25,6 +25,10 @@
     measurement_is_compatible,
     GALAXY_LENS_TYPES,
     GALAXY_SOURCE_TYPES,
+    ALL_MEASUREMENT_TYPES,
+    Galaxies,
+    CMB,
+    Clusters,
 )
 
 # TwoPointRealIndex is a type used to create intermediate objects when reading SACC
@@ -48,6 +52,63 @@
 )
 
 
+def make_measurement(value: Measurement | dict[str, Any]) -> Measurement:
+    """Create a Measurement object from a dictionary."""
+    if isinstance(value, ALL_MEASUREMENT_TYPES):
+        return value
+
+    if not isinstance(value, dict):
+        raise ValueError(f"Invalid Measurement: {value} is not a dictionary")
+
+    if "subject" not in value:
+        raise ValueError("Invalid Measurement: dictionary does not contain 'subject'")
+
+    subject = value["subject"]
+
+    match subject:
+        case "Galaxies":
+            return Galaxies[value["property"]]
+        case "CMB":
+            return CMB[value["property"]]
+        case "Clusters":
+            return Clusters[value["property"]]
+        case _:
+            raise ValueError(
+                f"Invalid Measurement: subject: '{subject}' is not recognized"
+            )
+
+
+def make_measurements(
+    value: set[Measurement] | list[dict[str, Any]],
+) -> set[Measurement]:
+    """Create a Measurement object from a dictionary."""
+    if isinstance(value, set) and all(
+        isinstance(v, ALL_MEASUREMENT_TYPES) for v in value
+    ):
+        return value
+
+    measurements: set[Measurement] = set()
+    for measurement_dict in value:
+        measurements.update([make_measurement(measurement_dict)])
+    return measurements
+
+
+def make_measurement_dict(value: Measurement) -> dict[str, str]:
+    """Create a dictionary from a Measurement object.
+
+    :param value: the measurement to turn into a dictionary
+    """
+    return {"subject": type(value).__name__, "property": value.name}
+
+
+def make_measurements_dict(value: set[Measurement]) -> list[dict[str, str]]:
+    """Create a dictionary from a Measurement object.
+
+    :param value: the measurement to turn into a dictionary
+    """
+    return [make_measurement_dict(measurement) for measurement in value]
+
+
 def _extract_all_candidate_measurement_types(
     data_points: list[sacc.DataPoint],
     include_maybe_types: bool = False,

firecrown/metadata_types.py

@@ -312,17 +312,6 @@ def __eq__(self, other):
         )
 
 
-def make_measurements_dict(value: set[Measurement]) -> list[dict[str, str]]:
-    """Create a dictionary from a Measurement object.
-
-    :param value: the measurement to turn into a dictionary
-    """
-    return [
-        {"subject": type(measurement).__name__, "property": measurement.name}
-        for measurement in value
-    ]
-
-
 def measurement_is_compatible(a: Measurement, b: Measurement) -> bool:
     """Check if two Measurement are compatible.
 

firecrown/utils.py

@@ -49,7 +49,9 @@ def base_model_from_yaml(cls: type, yaml_str: str):
 
 def base_model_to_yaml(model: BaseModel) -> str:
     """Convert a base model to a yaml string."""
-    return yaml.dump(model.model_dump(), default_flow_style=False, sort_keys=False)
+    return yaml.dump(
+        model.model_dump(), default_flow_style=None, sort_keys=False, width=80
+    )
 
 
 def upper_triangle_indices(n: int) -> Generator[tuple[int, int], None, None]:

firecrown/version.py

@@ -7,6 +7,6 @@
 """
 
 FIRECROWN_MAJOR = 1
-FIRECROWN_MINOR = 8
-FIRECROWN_PATCH = 0
+FIRECROWN_MINOR = 9
+FIRECROWN_PATCH = "0a0"
 __version__ = f"{FIRECROWN_MAJOR}.{FIRECROWN_MINOR}.{FIRECROWN_PATCH}"

pylintrc

@@ -12,9 +12,6 @@ py-version=3.10
 # Discover python modules and packages in the file system subtree.
 recursive=yes
 
-# Add custom pylint plugins
-load-plugins=pylint_plugins.duplicate_code
-
 [MESSAGES CONTROL]
 
 # Enable the message, report, category or checker with the given id(s). You can

tests/conftest.py

@@ -206,6 +206,25 @@ def make_harmonic_bin_2(request) -> InferredGalaxyZDist:
     return x
 
 
+@pytest.fixture(
+    name="all_harmonic_bins",
+)
+def make_all_harmonic_bins() -> list[InferredGalaxyZDist]:
+    """Generate a list of InferredGalaxyZDist objects with 5 bins."""
+    z = np.linspace(0.0, 1.0, 256)
+    dndzs = [
+        np.exp(-0.5 * (z - 0.5) ** 2 / 0.05**2) / (np.sqrt(2 * np.pi) * 0.05),
+        np.exp(-0.5 * (z - 0.6) ** 2 / 0.05**2) / (np.sqrt(2 * np.pi) * 0.05),
+    ]
+    return [
+        InferredGalaxyZDist(
+            bin_name=f"bin_{i + 1}", z=z, dndz=dndzs[i], measurements={m}
+        )
+        for i in range(2)
+        for m in [Galaxies.COUNTS, Galaxies.SHEAR_E]
+    ]
+
+
 @pytest.fixture(
     name="real_bin_1",
     params=[
@@ -248,6 +267,28 @@ def make_real_bin_2(request) -> InferredGalaxyZDist:
     return x
 
 
+@pytest.fixture(
+    name="all_real_bins",
+)
+def make_all_real_bins() -> list[InferredGalaxyZDist]:
+    """Generate a list of InferredGalaxyZDist objects with 5 bins."""
+    return [
+        InferredGalaxyZDist(
+            bin_name=f"bin_{i + 1}",
+            z=np.linspace(0, 1, 5),
+            dndz=np.array([0.1, 0.5, 0.2, 0.3, 0.4]),
+            measurements={m},
+        )
+        for i in range(2)
+        for m in [
+            Galaxies.COUNTS,
+            Galaxies.SHEAR_T,
+            Galaxies.SHEAR_MINUS,
+            Galaxies.SHEAR_PLUS,
+        ]
+    ]
+
+
 @pytest.fixture(name="window_1")
 def make_window_1() -> tuple[npt.NDArray[np.int64], npt.NDArray[np.float64]]:
     """Generate a Window object with 100 ells."""

tests/generators/test_inferred_galaxy_zdist.py

@@ -467,7 +467,7 @@ def test_make_measurement_from_dictionary():
         _ = make_measurements([{}])
 
     with pytest.raises(
-        ValueError, match=re.escape(r"Invalid Measurement: {3} is not a dictionary")
+        ValueError, match=re.escape(r"Invalid Measurement: 3 is not a dictionary")
     ):
         _ = make_measurements({3})  # type: ignore
 

tests/likelihood/test_factories.py

@@ -40,7 +40,16 @@ def test_two_point_cells_with_data(harmonic_two_point_xy: TwoPointXY):
 
 def test_two_point_two_point_cwindow_with_data(harmonic_two_point_xy: TwoPointXY):
     ells = np.array(np.linspace(0, 100, 100), dtype=np.int64)
-    weights = np.ones(400).reshape(-1, 4)
+    weights = np.zeros((100, 4), dtype=np.float64)
+    # Create a window with 4 bins, each containing 25 elements with a weight of 1.0.
+    # The bins are defined as follows:
+    # - Bin 1: Elements 0 to 24
+    # - Bin 2: Elements 25 to 49
+    # - Bin 3: Elements 50 to 74
+    # - Bin 4: Elements 75 to 99
+    rows = np.arange(100)
+    cols = rows // 25
+    weights[rows, cols] = 1.0
 
     ells = np.array(np.linspace(0, 100, 100), dtype=np.int64)
     data = np.array(np.zeros(4) + 1.1, dtype=np.float64)

tests/metadata/test_data_functions.py

@@ -8,4 +8,4 @@
 
 
 def test_version():
-    assert firecrown.__version__ == "1.8.0"
+    assert firecrown.__version__ == "1.9.0a0"

tests/metadata/test_metadata_two_point_measurement.py

@@ -58,7 +58,7 @@ format:
 
 reference-location: margin
 citation-location: margin
-subtitle: "version 1.8.0"
+subtitle: "version 1.9.0a0"
 authors:
   - Marc Paterno
   - Sandro Vitenti

tests/test_version.py

@@ -1,6 +1,6 @@
 ---
 title: An Introduction to Firecrown
-subtitle: "Version 1.8.0"
+subtitle: "Version 1.9.0a0"
 authors:
   - Marc Paterno
   - Sandro Vitenti

tutorial/_quarto.yml

@@ -37,7 +37,7 @@ sacc_data = sacc.Sacc.load_fits("../examples/des_y1_3x2pt/sacc_data.fits")
 all_meta = extract_all_real_metadata_indices(sacc_data)
 ```
 
-The metadata can be seem below:
+The metadata can be seen below:
 ```{python}
 # | code-fold: true
 import yaml
@@ -194,6 +194,10 @@ from firecrown.modeling_tools import ModelingTools
 from firecrown.ccl_factory import CCLFactory
 from firecrown.updatable import get_default_params_map
 from firecrown.parameters import ParamsMap
+from firecrown.utils import base_model_to_yaml
+from firecrown.data_functions import TwoPointBinFilterCollection, TwoPointBinFilter
+from firecrown.metadata_types import Galaxies
+
 
 tools = ModelingTools(ccl_factory=CCLFactory(require_nonlinear_pk=True))
 params = get_default_params_map(tools, likelihood)
@@ -211,3 +215,157 @@ print(f"Loglike from metadata only: {likelihood.compute_loglike(tools)}")
 print(f"Loglike from ready state: {likelihood_ready.compute_loglike(tools)}")
 ```
 
+## Filtering Data: Scale-cuts
+
+Real analyses use only a subset of the measured two-points statistics, where the utilized data is typically limited my the accuracy of the models used to fit the data. 
+It is then useful to define the physical scales (corresponding to the data) that should be analyzed in a given likelihood evaluation of two-point statistics. 
+Firecrown can implement this feature though its factories, notably by defining a `TwoPointBinFilterCollection` object. 
+This object is a collection of `TwoPointBinFilter` objects, which define the valid data analysis range for a given combination of two-point tracers. 
+For instance, we can define the filtered range of galaxy clustering auto-correlations as follows:
+
+```{python}
+tp_collection = TwoPointBinFilterCollection(
+    filters=[
+        TwoPointBinFilter.from_args(
+            name1=f"lens{i}",
+            measurement1=Galaxies.COUNTS,
+            name2=f"lens{i}",
+            measurement2=Galaxies.COUNTS,
+            lower=2,
+            upper=300,
+        )
+        for i in range(5)
+    ],
+    require_filter_for_all=True,
+    allow_empty=True,
+)
+Markdown(f"```yaml\n{base_model_to_yaml(tp_collection)}\n```")
+```
+
+Equivalently, we may reduce the complexity of the code slightly and specify the use of auto-correlations only:
+
+```{python}
+tp_collection = TwoPointBinFilterCollection(
+                filters=[
+                    TwoPointBinFilter.from_args_auto(
+                        name=f"lens{i}",
+                        measurement=Galaxies.COUNTS,
+                        lower=2,
+                        upper=300,
+                    )
+                    for i in range(5)
+                ],
+                require_filter_for_all=True,
+                allow_empty=True,
+)
+Markdown(f"```yaml\n{base_model_to_yaml(tp_collection)}\n```")
+```
+
+One may alternatively define the tracers directly (instead of from arguments) as `TwoPointTracerSpec` objects. 
+
+A `TwoPointExperiment` object is able to keep track of the relevant `Factory` instances to generate the two-point configurations of the analysis (either in configuration or harmonic space) and the scale-cut/data filtering choices to evaluate a defined likelihood. 
+The interpretation of the filtered lower and upper limits of the data depend on the definition of the `TwoPointExperiment` factories in either configuration or harmonic space.
+
+With this formalism, we are able to evaluate the likelihood exactly as the previous section by defining filters to be very wide. 
+Alternatively, by setting a restrictively small filtered range, we can remove data from the analysis and do so in the example below by filtering-out all galaxy clustering data. 
+
+```{python}
+from firecrown.likelihood.factories import (
+    DataSourceSacc,
+    TwoPointCorrelationSpace,
+    TwoPointExperiment,
+    TwoPointFactory,
+)
+
+tpf = TwoPointFactory(
+    correlation_space=TwoPointCorrelationSpace.REAL,
+    weak_lensing_factory=weak_lensing_factory,
+    number_counts_factory=number_counts_factory,
+)
+
+two_point_experiment = TwoPointExperiment(
+    two_point_factory=tpf,
+    data_source=DataSourceSacc(
+        sacc_data_file="../examples/des_y1_3x2pt/sacc_data.fits",
+        filters=TwoPointBinFilterCollection(
+            require_filter_for_all=False,
+            allow_empty=True,
+            filters=[
+                TwoPointBinFilter.from_args_auto(
+                    name=f"lens{i}",
+                    measurement=Galaxies.COUNTS,
+                    lower=0.5,
+                    upper=300,
+                )
+                for i in range(5)
+            ],
+        ),
+    ),
+)
+
+two_point_experiment_filtered = TwoPointExperiment(
+    two_point_factory=tpf,
+    data_source=DataSourceSacc(
+        sacc_data_file="../examples/des_y1_3x2pt/sacc_data.fits",
+        filters=TwoPointBinFilterCollection(
+            require_filter_for_all=False,
+            allow_empty=True,
+            filters=[
+                TwoPointBinFilter.from_args_auto(
+                    name=f"lens{i}",
+                    measurement=Galaxies.COUNTS,
+                    lower=2999,
+                    upper=3000,
+                )
+                for i in range(5)
+            ],
+        ),
+    ),
+)
+```
+
+The `TwoPointExperiment` objects can also be used to create likelihoods in the ready state. 
+Additionally, they can be serialized into a yaml file, making it easier to share specific analysis choices with other users and collaborators.
+
+The `yaml` below shows the first experiment.
+```{python}
+# | code-fold: true
+Markdown(f"```yaml\n{base_model_to_yaml(two_point_experiment)}\n```")
+```
+
+The `yaml` below shows the second experiment.
+```{python}
+# | code-fold: true
+Markdown(f"```yaml\n{base_model_to_yaml(two_point_experiment_filtered)}\n```")
+```
+
+Next, we can create likelihoods from the `TwoPointExperiment` objects and compare the loglike values.
+
+```{python}
+likelihood_tpe = two_point_experiment.make_likelihood()
+
+params = get_default_params_map(tools, likelihood_tpe)
+
+tools = ModelingTools()
+tools.update(params)
+tools.prepare()
+likelihood_tpe.update(params)
+
+likelihood_tpe_filtered = two_point_experiment_filtered.make_likelihood()
+
+params = get_default_params_map(tools, likelihood_tpe_filtered)
+
+tools = ModelingTools()
+tools.update(params)
+tools.prepare()
+likelihood_tpe_filtered.update(params)
+
+```
+
+```{python}
+# | code-fold: true
+print(f"Loglike from metadata only: {likelihood.compute_loglike(tools)}")
+print(f"Loglike from ready state: {likelihood_ready.compute_loglike(tools)}")
+print(f"Loglike from TwoPointExperiment: {likelihood_tpe.compute_loglike(tools)}")
+print(f"Loglike from filtered TwoPointExperiment: {likelihood_tpe_filtered.compute_loglike(tools)}")
+```

tutorial/introduction_to_firecrown.qmd

@@ -93,7 +93,7 @@ wl_photoz = wl.PhotoZShiftFactory()
 wl_mult_bias = wl.MultiplicativeShearBiasFactory()
 
 # NumberCounts systematics -- global
-# As for Firecrown 1.8.0, we do not have any global systematics for number counts
+# As for Firecrown 1.9.0a0, we do not have any global systematics for number counts
 # NumberCounts systematics -- per-bin
 nc_photoz = nc.PhotoZShiftFactory()