diff --git a/docs/non-developer-mode-example/des_y1_3x2pt.py b/docs/non-developer-mode-example/des_y1_3x2pt.py index 58e0772ea..4dcfc1110 100755 --- a/docs/non-developer-mode-example/des_y1_3x2pt.py +++ b/docs/non-developer-mode-example/des_y1_3x2pt.py @@ -2,12 +2,12 @@ import os -import firecrown.likelihood.gauss_family.statistic.source.weak_lensing as wl -import firecrown.likelihood.gauss_family.statistic.source.number_counts as nc +import firecrown.likelihood.weak_lensing as wl +import firecrown.likelihood.number_counts as nc -from firecrown.likelihood.gauss_family.statistic.two_point import TwoPoint +from firecrown.likelihood.two_point import TwoPoint -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian +from firecrown.likelihood.gaussian import ConstGaussian import sacc diff --git a/examples/cluster_number_counts/cluster_redshift_richness.py b/examples/cluster_number_counts/cluster_redshift_richness.py index 15e0799c6..a0ad5c2ae 100644 --- a/examples/cluster_number_counts/cluster_redshift_richness.py +++ b/examples/cluster_number_counts/cluster_redshift_richness.py @@ -5,8 +5,8 @@ import pyccl as ccl import sacc -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian -from firecrown.likelihood.gauss_family.statistic.binned_cluster_number_counts import ( +from firecrown.likelihood.gaussian import ConstGaussian +from firecrown.likelihood.binned_cluster_number_counts import ( BinnedClusterNumberCounts, ) from firecrown.likelihood.likelihood import Likelihood, NamedParameters diff --git a/examples/cosmicshear/cosmicshear.py b/examples/cosmicshear/cosmicshear.py index dead387fb..c0b8d9337 100644 --- a/examples/cosmicshear/cosmicshear.py +++ b/examples/cosmicshear/cosmicshear.py @@ -4,9 +4,9 @@ import sacc -import firecrown.likelihood.gauss_family.statistic.source.weak_lensing as wl -from firecrown.likelihood.gauss_family.statistic.two_point import TwoPoint -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian +import firecrown.likelihood.weak_lensing as wl +from firecrown.likelihood.two_point import TwoPoint +from firecrown.likelihood.gaussian import ConstGaussian def build_likelihood(_): diff --git a/examples/des_y1_3x2pt/des_y1_3x2pt.py b/examples/des_y1_3x2pt/des_y1_3x2pt.py index 8dd967f04..d51e79b24 100755 --- a/examples/des_y1_3x2pt/des_y1_3x2pt.py +++ b/examples/des_y1_3x2pt/des_y1_3x2pt.py @@ -3,10 +3,10 @@ import os import sacc -import firecrown.likelihood.gauss_family.statistic.source.weak_lensing as wl -import firecrown.likelihood.gauss_family.statistic.source.number_counts as nc -from firecrown.likelihood.gauss_family.statistic.two_point import TwoPoint -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian +import firecrown.likelihood.weak_lensing as wl +import firecrown.likelihood.number_counts as nc +from firecrown.likelihood.two_point import TwoPoint +from firecrown.likelihood.gaussian import ConstGaussian # The likelihood used for DES Y1 3x2pt analysis is a Gaussian likelihood, which diff --git a/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py b/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py index ac317e04c..21b7532e6 100755 --- a/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py +++ b/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py @@ -9,14 +9,14 @@ import pyccl as ccl import pyccl.nl_pt -import firecrown.likelihood.gauss_family.statistic.source.weak_lensing as wl -import firecrown.likelihood.gauss_family.statistic.source.number_counts as nc -from firecrown.likelihood.gauss_family.statistic.two_point import ( +import firecrown.likelihood.weak_lensing as wl +import firecrown.likelihood.number_counts as nc +from firecrown.likelihood.two_point import ( TwoPoint, TracerNames, TRACER_NAMES_TOTAL, ) -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian +from firecrown.likelihood.gaussian import ConstGaussian from firecrown.parameters import ParamsMap from firecrown.modeling_tools import ModelingTools from firecrown.likelihood.likelihood import Likelihood diff --git a/examples/des_y1_3x2pt/des_y1_cosmic_shear_TATT.py b/examples/des_y1_3x2pt/des_y1_cosmic_shear_TATT.py index eb836e867..640726aae 100644 --- a/examples/des_y1_3x2pt/des_y1_cosmic_shear_TATT.py +++ b/examples/des_y1_3x2pt/des_y1_cosmic_shear_TATT.py @@ -5,9 +5,9 @@ import pyccl as ccl import pyccl.nl_pt -import firecrown.likelihood.gauss_family.statistic.source.weak_lensing as wl -from firecrown.likelihood.gauss_family.statistic.two_point import TwoPoint -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian +import firecrown.likelihood.weak_lensing as wl +from firecrown.likelihood.two_point import TwoPoint +from firecrown.likelihood.gaussian import ConstGaussian from firecrown.parameters import ParamsMap from firecrown.modeling_tools import ModelingTools from firecrown.likelihood.likelihood import Likelihood diff --git a/examples/des_y1_3x2pt/des_y1_cosmic_shear_pk_modifier.py b/examples/des_y1_3x2pt/des_y1_cosmic_shear_pk_modifier.py index 9afa95856..6166b1fcf 100644 --- a/examples/des_y1_3x2pt/des_y1_cosmic_shear_pk_modifier.py +++ b/examples/des_y1_3x2pt/des_y1_cosmic_shear_pk_modifier.py @@ -9,9 +9,9 @@ import pyccl as ccl import pyccl.nl_pt -import firecrown.likelihood.gauss_family.statistic.source.weak_lensing as wl -from firecrown.likelihood.gauss_family.statistic.two_point import TwoPoint -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian +import firecrown.likelihood.weak_lensing as wl +from firecrown.likelihood.two_point import TwoPoint +from firecrown.likelihood.gaussian import ConstGaussian from firecrown.parameters import ParamsMap, create from firecrown.modeling_tools import ModelingTools, PowerspectrumModifier from firecrown.likelihood.likelihood import Likelihood diff --git a/examples/srd_sn/sn_srd.py b/examples/srd_sn/sn_srd.py index 67a0507ea..3a59fd46e 100644 --- a/examples/srd_sn/sn_srd.py +++ b/examples/srd_sn/sn_srd.py @@ -1,8 +1,8 @@ """Demonstration of the use of the :class:`Supernova` statistics object.""" import sacc -import firecrown.likelihood.gauss_family.statistic.supernova as sn -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian +import firecrown.likelihood.supernova as sn +from firecrown.likelihood.gaussian import ConstGaussian from firecrown.likelihood.likelihood import NamedParameters diff --git a/firecrown/connector/cosmosis/likelihood.py b/firecrown/connector/cosmosis/likelihood.py index f5084b97d..e58823b42 100644 --- a/firecrown/connector/cosmosis/likelihood.py +++ b/firecrown/connector/cosmosis/likelihood.py @@ -14,8 +14,8 @@ import pyccl as ccl from firecrown.connector.mapping import mapping_builder, MappingCosmoSIS -from firecrown.likelihood.gauss_family.gauss_family import GaussFamily -from firecrown.likelihood.gauss_family.statistic.two_point import TwoPoint +from firecrown.likelihood.gaussfamily import GaussFamily +from firecrown.likelihood.two_point import TwoPoint from firecrown.likelihood.likelihood import load_likelihood, Likelihood, NamedParameters from firecrown.parameters import ParamsMap from firecrown.updatable import MissingSamplerParameterError diff --git a/firecrown/connector/numcosmo/numcosmo.py b/firecrown/connector/numcosmo/numcosmo.py index 4562ea11f..f84557576 100644 --- a/firecrown/connector/numcosmo/numcosmo.py +++ b/firecrown/connector/numcosmo/numcosmo.py @@ -13,7 +13,7 @@ from firecrown.likelihood.likelihood import load_likelihood from firecrown.likelihood.likelihood import Likelihood from firecrown.likelihood.likelihood import NamedParameters -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian +from firecrown.likelihood.gaussian import ConstGaussian from firecrown.parameters import ParamsMap from firecrown.connector.mapping import Mapping, build_ccl_background_dict from firecrown.modeling_tools import ModelingTools diff --git a/firecrown/likelihood/binned_cluster_number_counts.py b/firecrown/likelihood/binned_cluster_number_counts.py new file mode 100644 index 000000000..2732134b5 --- /dev/null +++ b/firecrown/likelihood/binned_cluster_number_counts.py @@ -0,0 +1,150 @@ +"""This module holds classes needed to predict the binned cluster number counts. + +The binned cluster number counts statistic predicts the number of galaxy +clusters within a single redshift and mass bin. +""" + +from __future__ import annotations + +import numpy as np +import sacc + +# firecrown is needed for backward compatibility; remove support for deprecated +# directory structure is removed. +import firecrown # pylint: disable=unused-import # noqa: F401 +from firecrown.likelihood.source import SourceSystematic +from firecrown.likelihood.statistic import ( + DataVector, + Statistic, + TheoryVector, +) +from firecrown.modeling_tools import ModelingTools +from firecrown.models.cluster.abundance_data import AbundanceData +from firecrown.models.cluster.binning import SaccBin +from firecrown.models.cluster.properties import ClusterProperty +from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe + + +class BinnedClusterNumberCounts(Statistic): + """The Binned Cluster Number Counts statistic. + + This class will make a prediction for the number of clusters in a z, mass bin + and compare that prediction to the data provided in the sacc file. + """ + + def __init__( + self, + cluster_properties: ClusterProperty, + survey_name: str, + cluster_recipe: ClusterRecipe, + systematics: None | list[SourceSystematic] = None, + ): + super().__init__() + self.systematics = systematics or [] + self.theory_vector: None | TheoryVector = None + self.cluster_properties = cluster_properties + self.survey_name = survey_name + self.cluster_recipe = cluster_recipe + self.data_vector = DataVector.from_list([]) + self.sky_area = 0.0 + self.bins: list[SaccBin] = [] + + def read(self, sacc_data: sacc.Sacc) -> None: + """Read the data for this statistic and mark it as ready for use.""" + # Build the data vector and indices needed for the likelihood + if self.cluster_properties == ClusterProperty.NONE: + raise ValueError("You must specify at least one cluster property.") + + sacc_adapter = AbundanceData(sacc_data) + self.sky_area = sacc_adapter.get_survey_tracer(self.survey_name).sky_area + + data, indices = sacc_adapter.get_observed_data_and_indices_by_survey( + self.survey_name, self.cluster_properties + ) + self.data_vector = DataVector.from_list(data) + self.sacc_indices = np.array(indices) + + self.bins = sacc_adapter.get_bin_edges( + self.survey_name, self.cluster_properties + ) + for bin_edge in self.bins: + if bin_edge.dimension != self.bins[0].dimension: + raise ValueError( + "The cluster number counts statistic requires all bins to be the " + "same dimension." + ) + + super().read(sacc_data) + + def get_data_vector(self) -> DataVector: + """Gets the statistic data vector.""" + assert self.data_vector is not None + return self.data_vector + + def _compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: + """Compute a statistic from sources, concrete implementation.""" + assert tools.cluster_abundance is not None + + theory_vector_list: list[float] = [] + cluster_counts = [] + + cluster_counts = self.get_binned_cluster_counts(tools) + + for cl_property in ClusterProperty: + include_prop = cl_property & self.cluster_properties + if not include_prop: + continue + + if cl_property == ClusterProperty.COUNTS: + theory_vector_list += cluster_counts + continue + + theory_vector_list += self.get_binned_cluster_property( + tools, cluster_counts, cl_property + ) + + return TheoryVector.from_list(theory_vector_list) + + def get_binned_cluster_property( + self, + tools: ModelingTools, + cluster_counts: list[float], + cluster_properties: ClusterProperty, + ) -> list[float]: + """Computes the mean mass of clusters in each bin. + + Using the data from the sacc file, this function evaluates the likelihood for + a single point of the parameter space, and returns the predicted mean mass of + the clusters in each bin. + """ + assert tools.cluster_abundance is not None + + mean_values = [] + for this_bin, counts in zip(self.bins, cluster_counts): + total_observable = self.cluster_recipe.evaluate_theory_prediction( + tools.cluster_abundance, this_bin, self.sky_area, cluster_properties + ) + cluster_counts.append(counts) + + mean_observable = total_observable / counts + mean_values.append(mean_observable) + + return mean_values + + def get_binned_cluster_counts(self, tools: ModelingTools) -> list[float]: + """Computes the number of clusters in each bin. + + Using the data from the sacc file, this function evaluates the likelihood for + a single point of the parameter space, and returns the predicted number of + clusters in each bin. + """ + assert tools.cluster_abundance is not None + + cluster_counts = [] + for this_bin in self.bins: + counts = self.cluster_recipe.evaluate_theory_prediction( + tools.cluster_abundance, this_bin, self.sky_area + ) + cluster_counts.append(counts) + + return cluster_counts diff --git a/firecrown/likelihood/gauss_family/__init__.py b/firecrown/likelihood/gauss_family/__init__.py index aad814012..9c0fa90a1 100644 --- a/firecrown/likelihood/gauss_family/__init__.py +++ b/firecrown/likelihood/gauss_family/__init__.py @@ -1,10 +1 @@ -""" - -Gaussian Family Package -======================= - -This package provide all necessary tools for Gaussian family likelihoods. - -""" - # flake8: noqa diff --git a/firecrown/likelihood/gauss_family/gauss_family.py b/firecrown/likelihood/gauss_family/gauss_family.py index 27ef563dd..347f6cd84 100644 --- a/firecrown/likelihood/gauss_family/gauss_family.py +++ b/firecrown/likelihood/gauss_family/gauss_family.py @@ -1,391 +1,20 @@ -"""Gaussian Family Module. +"""Deprecated module with classes related to the Gaussian family of statistics.""" -Some notes. -""" +# flake8: noqa -from __future__ import annotations - -from enum import Enum -from functools import wraps -from typing import Sequence, Callable, TypeVar -from typing import final import warnings -from typing_extensions import ParamSpec -import numpy as np -import numpy.typing as npt -import scipy.linalg - -import sacc - -from firecrown.parameters import ParamsMap -from firecrown.likelihood.likelihood import Likelihood -from firecrown.modeling_tools import ModelingTools -from firecrown.updatable import UpdatableCollection -from firecrown.likelihood.gauss_family.statistic.statistic import ( - Statistic, - GuardedStatistic, -) -from firecrown.utils import save_to_sacc - - -class State(Enum): - """The states used in GaussFamily.""" - - INITIALIZED = 1 - READY = 2 - UPDATED = 3 - COMPUTED = 4 - - -T = TypeVar("T") -P = ParamSpec("P") - - -# See https://peps.python.org/pep-0612/ and -# https://stackoverflow.com/questions/66408662/type-annotations-for-decorators -# for how to specify the types of *args and **kwargs, and the return type of -# the method being decorated. - - -# Beware -def enforce_states( - *, - initial: State | list[State], - terminal: None | State = None, - failure_message: str, -) -> Callable[[Callable[P, T]], Callable[P, T]]: - """This decorator wraps a method, and enforces state machine behavior. - - If the object is not in one of the states in initial, an - AssertionError is raised with the given failure_message. - If terminal is None the state of the object is not modified. - If terminal is not None and the call to the wrapped method returns - normally the state of the object is set to terminal. - """ - initials: list[State] - if isinstance(initial, list): - initials = initial - else: - initials = [initial] - - def decorator_enforce_states(func: Callable[P, T]) -> Callable[P, T]: - """Part of the decorator which is the closure. - - This closure is what actually contains the values of initials, terminal, and - failure_message. - """ - - @wraps(func) - def wrapper_repeat(*args: P.args, **kwargs: P.kwargs) -> T: - """Part of the decorator which is the actual wrapped method. - - It is responsible for confirming a correct initial state, and - establishing the correct final state if the wrapped method - succeeds. - """ - # The odd use of args[0] instead of self seems to be the only way - # to have both the Python runtime and mypy agree on what is being - # passed to the method, and to allow access to the attribute - # 'state'. Recall that the syntax: - # o.foo() - # calls a bound function object accessible as o.foo; this bound - # function object calls the function foo() passing 'o' as the - # first argument, self. - assert isinstance(args[0], GaussFamily) - assert args[0].state in initials, failure_message - value = func(*args, **kwargs) - if terminal is not None: - args[0].state = terminal - return value - - return wrapper_repeat - - return decorator_enforce_states - - -class GaussFamily(Likelihood): - """GaussFamily is the base class for likelihoods based on a chi-squared calculation. - - It provides an implementation of Likelihood.compute_chisq. Derived classes must - implement the abstract method compute_loglike, which is inherited from Likelihood. - - GaussFamily (and all classes that inherit from it) must abide by the the - following rules regarding the order of calling of methods. - - 1. after a new object is created, :meth:`read` must be called before any - other method in the interfaqce. - 2. after :meth:`read` has been called it is legal to call - :meth:`get_data_vector`, or to call :meth:`update`. - 3. after :meth:`update` is called it is then legal to call - :meth:`calculate_loglike` or :meth:`get_data_vector`, or to reset - the object (returning to the pre-update state) by calling - :meth:`reset`. It is also legal to call :meth:`compute_theory_vector`. - 4. after :meth:`compute_theory_vector` is called it is legal to call - :meth:`get_theory_vector` to retrieve the already-calculated theory - vector. - - This state machine behavior is enforced through the use of the decorator - :meth:`enforce_states`, above. - """ - - def __init__( - self, - statistics: Sequence[Statistic], - ): - """Initialize the base class parts of a GaussFamily object.""" - super().__init__() - self.state: State = State.INITIALIZED - if len(statistics) == 0: - raise ValueError("GaussFamily requires at least one statistic") - - for i, s in enumerate(statistics): - if not isinstance(s, Statistic): - raise ValueError( - f"statistics[{i}] is not an instance of Statistic: {s}" - f"it is a {type(s)} instead." - ) - - self.statistics: UpdatableCollection[GuardedStatistic] = UpdatableCollection( - GuardedStatistic(s) for s in statistics - ) - self.cov: None | npt.NDArray[np.float64] = None - self.cholesky: None | npt.NDArray[np.float64] = None - self.inv_cov: None | npt.NDArray[np.float64] = None - self.cov_index_map: None | dict[int, int] = None - self.theory_vector: None | npt.NDArray[np.double] = None - self.data_vector: None | npt.NDArray[np.double] = None +import firecrown.likelihood.gaussfamily - @enforce_states( - initial=State.READY, - terminal=State.UPDATED, - failure_message="read() must be called before update()", - ) - def _update(self, _: ParamsMap) -> None: - """Handle the state resetting required by :class:`GaussFamily` likelihoods. +# pylint: disable=unused-import,unused-wildcard-import,wildcard-import +from firecrown.likelihood.gaussfamily import * - Any derived class that needs to implement :meth:`_update` - for its own reasons must be sure to do what this does: check the state - at the start of the method, and change the state at the end of the - method. - """ +# pylint: enable=unused-import,unused-wildcard-import,wildcard-import - @enforce_states( - initial=[State.UPDATED, State.COMPUTED], - terminal=State.READY, - failure_message="update() must be called before reset()", - ) - def _reset(self) -> None: - """Handle the state resetting required by :class:`GaussFamily` likelihoods. +assert not hasattr(firecrown.likelihood.gaussfamily, "__all__") - Any derived class that needs to implement :meth:`reset` - for its own reasons must be sure to do what this does: check the state - at the start of the method, and change the state at the end of the - method. - """ - self.theory_vector = None - - @enforce_states( - initial=State.INITIALIZED, - terminal=State.READY, - failure_message="read() must only be called once", - ) - def read(self, sacc_data: sacc.Sacc) -> None: - """Read the covariance matrix for this likelihood from the SACC file.""" - if sacc_data.covariance is None: - msg = ( - f"The {type(self).__name__} likelihood requires a covariance, " - f"but the SACC data object being read does not have one." - ) - raise RuntimeError(msg) - - covariance = sacc_data.covariance.dense - - indices_list = [] - data_vector_list = [] - for stat in self.statistics: - stat.read(sacc_data) - if stat.statistic.sacc_indices is None: - raise RuntimeError( - f"The statistic {stat.statistic} has no sacc_indices." - ) - indices_list.append(stat.statistic.sacc_indices.copy()) - data_vector_list.append(stat.statistic.get_data_vector()) - - indices = np.concatenate(indices_list) - data_vector = np.concatenate(data_vector_list) - cov = np.zeros((len(indices), len(indices))) - - for new_i, old_i in enumerate(indices): - for new_j, old_j in enumerate(indices): - cov[new_i, new_j] = covariance[old_i, old_j] - - self.data_vector = data_vector - self.cov_index_map = {old_i: new_i for new_i, old_i in enumerate(indices)} - self.cov = cov - self.cholesky = scipy.linalg.cholesky(self.cov, lower=True) - self.inv_cov = np.linalg.inv(cov) - - @final - @enforce_states( - initial=[State.READY, State.UPDATED, State.COMPUTED], - failure_message="read() must be called before get_cov()", - ) - def get_cov( - self, statistic: Statistic | list[Statistic] | None = None - ) -> npt.NDArray[np.float64]: - """Gets the current covariance matrix. - - :param statistic: The statistic for which the sub-covariance matrix - should be returned. If not specified, return the covariance of all - statistics. - """ - assert self.cov is not None - if statistic is None: - return self.cov - - assert self.cov_index_map is not None - if isinstance(statistic, Statistic): - statistic_list = [statistic] - else: - statistic_list = statistic - indices: list[int] = [] - for stat in statistic_list: - assert stat.sacc_indices is not None - temp = [self.cov_index_map[idx] for idx in stat.sacc_indices] - indices += temp - ixgrid = np.ix_(indices, indices) - - return self.cov[ixgrid] - - @final - @enforce_states( - initial=[State.READY, State.UPDATED, State.COMPUTED], - failure_message="read() must be called before get_data_vector()", - ) - def get_data_vector(self) -> npt.NDArray[np.float64]: - """Get the data vector from all statistics in the right order.""" - assert self.data_vector is not None - return self.data_vector - - @final - @enforce_states( - initial=[State.UPDATED, State.COMPUTED], - terminal=State.COMPUTED, - failure_message="update() must be called before compute_theory_vector()", - ) - def compute_theory_vector(self, tools: ModelingTools) -> npt.NDArray[np.float64]: - """Computes the theory vector using the current instance of pyccl.Cosmology. - - :param tools: Current ModelingTools object - """ - theory_vector_list: list[npt.NDArray[np.float64]] = [ - stat.compute_theory_vector(tools) for stat in self.statistics - ] - self.theory_vector = np.concatenate(theory_vector_list) - return self.theory_vector - - @final - @enforce_states( - initial=State.COMPUTED, - failure_message="compute_theory_vector() must be called before " - "get_theory_vector()", - ) - def get_theory_vector(self) -> npt.NDArray[np.float64]: - """Get the theory vector from all statistics in the right order.""" - assert ( - self.theory_vector is not None - ), "theory_vector is None after compute_theory_vector() has been called" - return self.theory_vector - - @final - @enforce_states( - initial=State.UPDATED, - failure_message="update() must be called before compute()", - ) - def compute( - self, tools: ModelingTools - ) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: - """Calculate and return both the data and theory vectors.""" - warnings.simplefilter("always", DeprecationWarning) - warnings.warn( - "The use of the `compute` method on Statistic is deprecated." - "The Statistic objects should implement `get_data` and " - "`compute_theory_vector` instead.", - category=DeprecationWarning, - ) - return self.get_data_vector(), self.compute_theory_vector(tools) - - @final - @enforce_states( - initial=[State.UPDATED, State.COMPUTED], - terminal=State.COMPUTED, - failure_message="update() must be called before compute_chisq()", - ) - def compute_chisq(self, tools: ModelingTools) -> float: - """Calculate and return the chi-squared for the given cosmology.""" - theory_vector: npt.NDArray[np.float64] - data_vector: npt.NDArray[np.float64] - residuals: npt.NDArray[np.float64] - try: - theory_vector = self.compute_theory_vector(tools) - data_vector = self.get_data_vector() - except NotImplementedError: - data_vector, theory_vector = self.compute(tools) - - assert len(data_vector) == len(theory_vector) - residuals = data_vector - theory_vector - - x = scipy.linalg.solve_triangular(self.cholesky, residuals, lower=True) - chisq = np.dot(x, x) - return chisq - - @enforce_states( - initial=[State.READY, State.UPDATED, State.COMPUTED], - failure_message="read() must be called before get_sacc_indices()", - ) - def get_sacc_indices( - self, statistic: Statistic | list[Statistic] | None = None - ) -> npt.NDArray[np.int64]: - """Get the SACC indices of the statistic or list of statistics. - - If no statistic is given, get the indices of all statistics of the likelihood. - """ - if statistic is None: - statistic = [stat.statistic for stat in self.statistics] - if isinstance(statistic, Statistic): - statistic = [statistic] - - sacc_indices_list = [] - for stat in statistic: - assert stat.sacc_indices is not None - sacc_indices_list.append(stat.sacc_indices.copy()) - - sacc_indices = np.concatenate(sacc_indices_list) - - return sacc_indices - - @enforce_states( - initial=State.COMPUTED, - failure_message="compute_theory_vector() must be called before " - "make_realization()", - ) - def make_realization( - self, sacc_data: sacc.Sacc, add_noise: bool = True, strict: bool = True - ) -> sacc.Sacc: - """Create a new realization of the model.""" - sacc_indices = self.get_sacc_indices() - - if add_noise: - new_data_vector = self.make_realization_vector() - else: - new_data_vector = self.get_theory_vector() - - new_sacc = save_to_sacc( - sacc_data=sacc_data, - data_vector=new_data_vector, - indices=sacc_indices, - strict=strict, - ) - - return new_sacc +warnings.warn( + "This module is deprecated. Use firecrown.likelihood.gaussfamily instead.", + DeprecationWarning, + stacklevel=2, +) diff --git a/firecrown/likelihood/gauss_family/gaussian.py b/firecrown/likelihood/gauss_family/gaussian.py index 5858db659..d8a1cc18b 100644 --- a/firecrown/likelihood/gauss_family/gaussian.py +++ b/firecrown/likelihood/gauss_family/gaussian.py @@ -1,25 +1,20 @@ -"""Provides GaussFamily concrete types.""" +"""Deprecated module with classes related to the Gaussian likelihood.""" -from __future__ import annotations -import numpy as np +# flake8: noqa -from firecrown.likelihood.gauss_family.gauss_family import GaussFamily -from firecrown.modeling_tools import ModelingTools +import warnings +import firecrown.likelihood.gaussian -class ConstGaussian(GaussFamily): - """A Gaussian log-likelihood with a constant covariance matrix.""" +# pylint: disable=unused-import,unused-wildcard-import,wildcard-import +from firecrown.likelihood.gaussian import * - def compute_loglike(self, tools: ModelingTools): - """Compute the log-likelihood.""" - return -0.5 * self.compute_chisq(tools) +# pylint: enable=unused-import,unused-wildcard-import,wildcard-import - def make_realization_vector(self) -> np.ndarray: - """Create a new realization of the model.""" - theory_vector = self.get_theory_vector() - assert self.cholesky is not None - new_data_vector = theory_vector + np.dot( - self.cholesky, np.random.randn(len(theory_vector)) - ) +assert not hasattr(firecrown.likelihood.gaussian, "__all__") - return new_data_vector +warnings.warn( + "This module is deprecated. Use firecrown.likelihood.gaussian instead.", + DeprecationWarning, + stacklevel=2, +) diff --git a/firecrown/likelihood/gauss_family/statistic/__init__.py b/firecrown/likelihood/gauss_family/statistic/__init__.py index a9fb3bdaa..9c0fa90a1 100644 --- a/firecrown/likelihood/gauss_family/statistic/__init__.py +++ b/firecrown/likelihood/gauss_family/statistic/__init__.py @@ -1,5 +1 @@ -"""This package provides statistics objects to be used by Gaussian Family -likelihoods. -""" - # flake8: noqa diff --git a/firecrown/likelihood/gauss_family/statistic/binned_cluster_number_counts.py b/firecrown/likelihood/gauss_family/statistic/binned_cluster_number_counts.py index 00aaa9259..32192ce02 100644 --- a/firecrown/likelihood/gauss_family/statistic/binned_cluster_number_counts.py +++ b/firecrown/likelihood/gauss_family/statistic/binned_cluster_number_counts.py @@ -1,147 +1,20 @@ -"""This module holds classes needed to predict the binned cluster number counts. +"""Deprecated module with classes to predict binned cluster number counts.""" -The binned cluster number counts statistic predicts the number of galaxy -clusters within a single redshift and mass bin. -""" +# flake8: noqa -from __future__ import annotations +import warnings -import numpy as np -import sacc +import firecrown.likelihood.binned_cluster_number_counts -from firecrown.likelihood.gauss_family.statistic.source.source import SourceSystematic -from firecrown.likelihood.gauss_family.statistic.statistic import ( - DataVector, - Statistic, - TheoryVector, -) -from firecrown.modeling_tools import ModelingTools -from firecrown.models.cluster.abundance_data import AbundanceData -from firecrown.models.cluster.binning import SaccBin -from firecrown.models.cluster.properties import ClusterProperty -from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe - - -class BinnedClusterNumberCounts(Statistic): - """The Binned Cluster Number Counts statistic. - - This class will make a prediction for the number of clusters in a z, mass bin - and compare that prediction to the data provided in the sacc file. - """ - - def __init__( - self, - cluster_properties: ClusterProperty, - survey_name: str, - cluster_recipe: ClusterRecipe, - systematics: None | list[SourceSystematic] = None, - ): - super().__init__() - self.systematics = systematics or [] - self.theory_vector: None | TheoryVector = None - self.cluster_properties = cluster_properties - self.survey_name = survey_name - self.cluster_recipe = cluster_recipe - self.data_vector = DataVector.from_list([]) - self.sky_area = 0.0 - self.bins: list[SaccBin] = [] - - def read(self, sacc_data: sacc.Sacc) -> None: - """Read the data for this statistic and mark it as ready for use.""" - # Build the data vector and indices needed for the likelihood - if self.cluster_properties == ClusterProperty.NONE: - raise ValueError("You must specify at least one cluster property.") - - sacc_adapter = AbundanceData(sacc_data) - self.sky_area = sacc_adapter.get_survey_tracer(self.survey_name).sky_area - - data, indices = sacc_adapter.get_observed_data_and_indices_by_survey( - self.survey_name, self.cluster_properties - ) - self.data_vector = DataVector.from_list(data) - self.sacc_indices = np.array(indices) - - self.bins = sacc_adapter.get_bin_edges( - self.survey_name, self.cluster_properties - ) - for bin_edge in self.bins: - if bin_edge.dimension != self.bins[0].dimension: - raise ValueError( - "The cluster number counts statistic requires all bins to be the " - "same dimension." - ) - - super().read(sacc_data) - - def get_data_vector(self) -> DataVector: - """Gets the statistic data vector.""" - assert self.data_vector is not None - return self.data_vector - - def _compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: - """Compute a statistic from sources, concrete implementation.""" - assert tools.cluster_abundance is not None - - theory_vector_list: list[float] = [] - cluster_counts = [] +# pylint: disable=unused-import,unused-wildcard-import,wildcard-import +from firecrown.likelihood.binned_cluster_number_counts import * - cluster_counts = self.get_binned_cluster_counts(tools) +# pylint: enable=unused-import,unused-wildcard-import,wildcard-import - for cl_property in ClusterProperty: - include_prop = cl_property & self.cluster_properties - if not include_prop: - continue +assert not hasattr(firecrown.likelihood.binned_cluster_number_counts, "__all__") - if cl_property == ClusterProperty.COUNTS: - theory_vector_list += cluster_counts - continue - - theory_vector_list += self.get_binned_cluster_property( - tools, cluster_counts, cl_property - ) - - return TheoryVector.from_list(theory_vector_list) - - def get_binned_cluster_property( - self, - tools: ModelingTools, - cluster_counts: list[float], - cluster_properties: ClusterProperty, - ) -> list[float]: - """Computes the mean mass of clusters in each bin. - - Using the data from the sacc file, this function evaluates the likelihood for - a single point of the parameter space, and returns the predicted mean mass of - the clusters in each bin. - """ - assert tools.cluster_abundance is not None - - mean_values = [] - for this_bin, counts in zip(self.bins, cluster_counts): - total_observable = self.cluster_recipe.evaluate_theory_prediction( - tools.cluster_abundance, this_bin, self.sky_area, cluster_properties - ) - cluster_counts.append(counts) - - mean_observable = total_observable / counts - mean_values.append(mean_observable) - - return mean_values - - def get_binned_cluster_counts(self, tools: ModelingTools) -> list[float]: - """Computes the number of clusters in each bin. - - Using the data from the sacc file, this function evaluates the likelihood for - a single point of the parameter space, and returns the predicted number of - clusters in each bin. - """ - assert tools.cluster_abundance is not None - - cluster_counts = [] - for this_bin in self.bins: - counts = self.cluster_recipe.evaluate_theory_prediction( - tools.cluster_abundance, this_bin, self.sky_area - ) - cluster_counts.append(counts) - - return cluster_counts +warnings.warn( + "This module is deprecated. Use firecrown.likelihood.binned_cluster_number_counts instead.", + DeprecationWarning, + stacklevel=2, +) diff --git a/firecrown/likelihood/gauss_family/statistic/source/__init__.py b/firecrown/likelihood/gauss_family/statistic/source/__init__.py index b57beda53..9c0fa90a1 100644 --- a/firecrown/likelihood/gauss_family/statistic/source/__init__.py +++ b/firecrown/likelihood/gauss_family/statistic/source/__init__.py @@ -1,5 +1 @@ -"""This package contains the base class :class:`Source` for :class:`TwoPoint` -sources and implementations. -""" - # flake8: noqa diff --git a/firecrown/likelihood/gauss_family/statistic/source/number_counts.py b/firecrown/likelihood/gauss_family/statistic/source/number_counts.py index 8979ad90b..51a2dc573 100644 --- a/firecrown/likelihood/gauss_family/statistic/source/number_counts.py +++ b/firecrown/likelihood/gauss_family/statistic/source/number_counts.py @@ -1,546 +1,20 @@ -"""Number counts source and systematics.""" +"""Deprecated module for number counts likelihoods.""" -from __future__ import annotations +# flake8: noqa -from abc import abstractmethod -from dataclasses import dataclass, replace -from typing import Sequence, final +import warnings -import numpy as np -import numpy.typing as npt -import pyccl +import firecrown.likelihood.number_counts -from firecrown import parameters -from firecrown.likelihood.gauss_family.statistic.source.source import ( - SourceGalaxy, - SourceGalaxyArgs, - SourceGalaxyPhotoZShift, - SourceGalaxySelectField, - SourceGalaxySystematic, - Tracer, -) -from firecrown.metadata.two_point import InferredGalaxyZDist -from firecrown.modeling_tools import ModelingTools -from firecrown.parameters import DerivedParameter, DerivedParameterCollection, ParamsMap -from firecrown.updatable import UpdatableCollection - -__all__ = ["NumberCounts"] - - -@dataclass(frozen=True) -class NumberCountsArgs(SourceGalaxyArgs): - """Class for number counts tracer builder argument.""" - - bias: None | npt.NDArray[np.float64] = None - mag_bias: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None - has_pt: bool = False - has_hm: bool = False - b_2: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None - b_s: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None - - -class NumberCountsSystematic(SourceGalaxySystematic[NumberCountsArgs]): - """Abstract base class for systematics for Number Counts sources. - - Derived classes must implement :python`apply` with the correct signature. - """ - - @abstractmethod - def apply( - self, tools: ModelingTools, tracer_arg: NumberCountsArgs - ) -> NumberCountsArgs: - """Apply method to include systematics in the tracer_arg.""" - - -class PhotoZShift(SourceGalaxyPhotoZShift[NumberCountsArgs]): - """Photo-z shift systematic.""" - - -class SelectField(SourceGalaxySelectField[NumberCountsArgs]): - """Systematic to select 3D field.""" - - -LINEAR_BIAS_DEFAULT_ALPHAZ = 0.0 -LINEAR_BIAS_DEFAULT_ALPHAG = 1.0 -LINEAR_BIAS_DEFAULT_Z_PIV = 0.5 - - -class LinearBiasSystematic(NumberCountsSystematic): - """Linear bias systematic. - - This systematic adds a linear bias model which varies with redshift and - the growth function. - - The following parameters are special Updatable parameters, which means that - they can be updated by the sampler, sacc_tracer is going to be used as a - prefix for the parameters: - - :ivar alphaz: the redshift exponent of the bias. - :ivar alphag: the growth function exponent of the bias. - :ivar z_piv: the pivot redshift of the bias. - """ - - def __init__(self, sacc_tracer: str): - """Initialize the LinearBiasSystematic. - - :param sacc_tracer: the name of the tracer in the SACC file. This is used - as a prefix for its parameters. - - """ - super().__init__(parameter_prefix=sacc_tracer) - - self.alphaz = parameters.register_new_updatable_parameter( - default_value=LINEAR_BIAS_DEFAULT_ALPHAZ - ) - self.alphag = parameters.register_new_updatable_parameter( - default_value=LINEAR_BIAS_DEFAULT_ALPHAG - ) - self.z_piv = parameters.register_new_updatable_parameter( - default_value=LINEAR_BIAS_DEFAULT_Z_PIV - ) - - def apply( - self, tools: ModelingTools, tracer_arg: NumberCountsArgs - ) -> NumberCountsArgs: - """Apply a linear bias systematic. - - Parameters - ---------- - cosmo : Cosmology - A Cosmology object. - tracer_arg : NumberCountsArgs - The source to which apply the shear bias. - """ - ccl_cosmo = tools.get_ccl_cosmology() - pref = ((1.0 + tracer_arg.z) / (1.0 + self.z_piv)) ** self.alphaz - pref *= ( - pyccl.growth_factor(ccl_cosmo, 1.0 / (1.0 + tracer_arg.z)) ** self.alphag - ) - - if tracer_arg.bias is None: - bias = np.ones_like(tracer_arg.z) - else: - bias = tracer_arg.bias - bias = bias * pref - - return replace( - tracer_arg, - bias=bias, - ) - - -PT_NON_LINEAR_BIAS_DEFAULT_B_2 = 1.0 -PT_NON_LINEAR_BIAS_DEFAULT_B_S = 1.0 - - -class PTNonLinearBiasSystematic(NumberCountsSystematic): - """Non-linear bias systematic. - - This systematic adds a linear bias model which varies with redshift and - the growth function. - - The following parameters are special Updatable parameters, which means that - they can be updated by the sampler, sacc_tracer is going to be used as a - prefix for the parameters: - - :ivar b_2: the quadratic bias. - :ivar b_s: the stochastic bias. - """ - - def __init__(self, sacc_tracer: str): - """Initialize the PTNonLinearBiasSystematic. - - :param sacc_tracer: the name of the tracer in the SACC file. This is used - as a prefix for its parameters. - - """ - super().__init__(parameter_prefix=sacc_tracer) - - self.b_2 = parameters.register_new_updatable_parameter( - default_value=PT_NON_LINEAR_BIAS_DEFAULT_B_2 - ) - self.b_s = parameters.register_new_updatable_parameter( - default_value=PT_NON_LINEAR_BIAS_DEFAULT_B_S - ) - - def apply( - self, tools: ModelingTools, tracer_arg: NumberCountsArgs - ) -> NumberCountsArgs: - z = tracer_arg.z - b_2_z = self.b_2 * np.ones_like(z) - b_s_z = self.b_s * np.ones_like(z) - # b_1 uses the "bias" field - return replace( - tracer_arg, - has_pt=True, - b_2=(z, b_2_z), - b_s=(z, b_s_z), - ) - - -class MagnificationBiasSystematic(NumberCountsSystematic): - """Magnification bias systematic. - - This systematic adds a magnification bias model for galaxy number contrast - following Joachimi & Bridle (2010), arXiv:0911.2454. - - The following parameters are special Updatable parameters, which means that - they can be updated by the sampler, sacc_tracer is going to be used as a - prefix for the parameters: - - :ivar r_lim: the limiting magnitude. - :ivar sig_c: the intrinsic dispersion of the source redshift distribution. - :ivar eta: the slope of the luminosity function. - :ivar z_c: the characteristic redshift of the source distribution. - :ivar z_m: the slope of the source redshift distribution. - """ - - def __init__(self, sacc_tracer: str): - """Initialize the MagnificationBiasSystematic. - - :param sacc_tracer: the name of the tracer in the SACC file. This is used - as a prefix for its parameters. - """ - super().__init__(parameter_prefix=sacc_tracer) - - self.r_lim = parameters.register_new_updatable_parameter() - self.sig_c = parameters.register_new_updatable_parameter() - self.eta = parameters.register_new_updatable_parameter() - self.z_c = parameters.register_new_updatable_parameter() - self.z_m = parameters.register_new_updatable_parameter() - - def apply( - self, tools: ModelingTools, tracer_arg: NumberCountsArgs - ) -> NumberCountsArgs: - """Apply a magnification bias systematic. - - :param tools: a ModelingTools object - :param tracer_arg: a NumberCountsArgs object - - :return: a NumberCountsArgs object - """ - z_bar = self.z_c + self.z_m * (self.r_lim - 24.0) - # The slope of log(n_tot(z,r_lim)) with respect to r_lim - # where n_tot(z,r_lim) is the luminosity function after using fit (C.1) - s = ( - self.eta / self.r_lim - - 3.0 * self.z_m / z_bar - + 1.5 * self.z_m * np.power(tracer_arg.z / z_bar, 1.5) / z_bar - ) - - if tracer_arg.mag_bias is None: - mag_bias = np.ones_like(tracer_arg.z) - else: - mag_bias = tracer_arg.mag_bias[1] - mag_bias = mag_bias * s / np.log(10) - - return replace( - tracer_arg, - mag_bias=(tracer_arg.z, mag_bias), - ) - - -CONSTANT_MAGNIFICATION_BIAS_DEFAULT_MAG_BIAS = 1.0 - - -class ConstantMagnificationBiasSystematic(NumberCountsSystematic): - """Simple constant magnification bias systematic. +# pylint: disable=unused-import,unused-wildcard-import,wildcard-import +from firecrown.likelihood.number_counts import * - This systematic adds a constant magnification bias model for galaxy number - contrast. +# pylint: enable=unused-import,unused-wildcard-import,wildcard-import - The following parameters are special Updatable parameters, which means that - they can be updated by the sampler, sacc_tracer is going to be used as a - prefix for the parameters: +assert not hasattr(firecrown.likelihood.number_counts, "__all__") - :ivar mag_bias: the magnification bias. - """ - - def __init__(self, sacc_tracer: str): - """Initialize the ConstantMagnificationBiasSystematic. - - :param sacc_tracer: the name of the tracer in the SACC file. This is used - as a prefix for its parameters. - """ - super().__init__(parameter_prefix=sacc_tracer) - - self.mag_bias = parameters.register_new_updatable_parameter( - default_value=CONSTANT_MAGNIFICATION_BIAS_DEFAULT_MAG_BIAS - ) - - def apply( - self, tools: ModelingTools, tracer_arg: NumberCountsArgs - ) -> NumberCountsArgs: - return replace( - tracer_arg, - mag_bias=(tracer_arg.z, np.ones_like(tracer_arg.z) * self.mag_bias), - ) - - -NUMBER_COUNTS_DEFAULT_BIAS = 1.5 - - -class NumberCounts(SourceGalaxy[NumberCountsArgs]): - """Source class for number counts.""" - - def __init__( - self, - *, - sacc_tracer: str, - has_rsd: bool = False, - derived_scale: bool = False, - scale: float = 1.0, - systematics: None | Sequence[SourceGalaxySystematic[NumberCountsArgs]] = None, - ): - """Initialize the NumberCounts object. - - :param sacc_tracer: the name of the tracer in the SACC file. This is used - as a prefix for its parameters. - :param has_rsd: whether to include RSD in the tracer. - :param derived_scale: whether to include a derived parameter for the scale - of the tracer. - :param scale: the initial scale of the tracer. - :param systematics: a list of systematics to apply to the tracer. - """ - super().__init__(sacc_tracer=sacc_tracer, systematics=systematics) - - self.sacc_tracer = sacc_tracer - self.has_rsd = has_rsd - self.derived_scale = derived_scale - - self.bias = parameters.register_new_updatable_parameter( - default_value=NUMBER_COUNTS_DEFAULT_BIAS - ) - self.systematics: UpdatableCollection[ - SourceGalaxySystematic[NumberCountsArgs] - ] = UpdatableCollection(systematics) - self.scale = scale - self.current_tracer_args: None | NumberCountsArgs = None - self.tracer_args: NumberCountsArgs - - @classmethod - def create_ready( - cls, - inferred_zdist: InferredGalaxyZDist, - has_rsd: bool = False, - derived_scale: bool = False, - scale: float = 1.0, - systematics: None | list[SourceGalaxySystematic[NumberCountsArgs]] = None, - ) -> NumberCounts: - """Create a NumberCounts object with the given tracer name and scale.""" - obj = cls( - sacc_tracer=inferred_zdist.bin_name, - systematics=systematics, - has_rsd=has_rsd, - derived_scale=derived_scale, - scale=scale, - ) - obj.tracer_args = NumberCountsArgs( - scale=obj.scale, - z=inferred_zdist.z, - dndz=inferred_zdist.dndz, - bias=None, - mag_bias=None, - ) - - return obj - - @final - def _update_source(self, params: ParamsMap): - """Perform any updates necessary after the parameters have being updated. - - This implementation must update all contained Updatable instances. - """ - self.systematics.update(params) - - @final - def _get_derived_parameters(self) -> DerivedParameterCollection: - if self.derived_scale: - assert self.current_tracer_args is not None - derived_scale = DerivedParameter( - "TwoPoint", - f"NumberCountsScale_{self.sacc_tracer}", - self.current_tracer_args.scale, - ) - derived_parameters = DerivedParameterCollection([derived_scale]) - else: - derived_parameters = DerivedParameterCollection([]) - - return derived_parameters - - def _read(self, sacc_data): - """Read the data for this source from the SACC file. - - Parameters - ---------- - sacc_data : sacc.Sacc - The data in the sacc format. - """ - self.tracer_args = NumberCountsArgs( - scale=self.scale, - z=np.array([]), - dndz=np.array([]), - bias=None, - mag_bias=None, - ) - super()._read(sacc_data) - - def create_tracers( - self, tools: ModelingTools - ) -> tuple[list[Tracer], NumberCountsArgs]: - """Create the tracers for this source.""" - tracer_args = self.tracer_args - tracer_args = replace(tracer_args, bias=self.bias * np.ones_like(tracer_args.z)) - - ccl_cosmo = tools.get_ccl_cosmology() - for systematic in self.systematics: - tracer_args = systematic.apply(tools, tracer_args) - - tracers = [] - - if not tracer_args.has_pt or tracer_args.mag_bias is not None or self.has_rsd: - # Create a normal pyccl.NumberCounts tracer if there's no PT, or - # in case there's magnification or RSD. - tracer_names = [] - if tracer_args.has_pt: - # use PT for galaxy bias - bias = None - else: - bias = (tracer_args.z, tracer_args.bias) - tracer_names += ["galaxies"] - if tracer_args.mag_bias is not None: - tracer_names += ["magnification"] - if self.has_rsd: - tracer_names += ["rsd"] - - ccl_mag_tracer = pyccl.NumberCountsTracer( - ccl_cosmo, - has_rsd=self.has_rsd, - dndz=(tracer_args.z, tracer_args.dndz), - bias=bias, - mag_bias=tracer_args.mag_bias, - ) - - tracers.append( - Tracer( - ccl_mag_tracer, - tracer_name="+".join(tracer_names), - field=tracer_args.field, - ) - ) - if tracer_args.has_pt: - nc_pt_tracer = pyccl.nl_pt.PTNumberCountsTracer( - b1=(tracer_args.z, tracer_args.bias), - b2=tracer_args.b_2, - bs=tracer_args.b_s, - ) - - ccl_nc_dummy_tracer = pyccl.NumberCountsTracer( - ccl_cosmo, - has_rsd=False, - dndz=(tracer_args.z, tracer_args.dndz), - bias=(tracer_args.z, np.ones_like(tracer_args.z)), - ) - nc_pt_tracer = Tracer( - ccl_nc_dummy_tracer, tracer_name="galaxies", pt_tracer=nc_pt_tracer - ) - tracers.append(nc_pt_tracer) - - self.current_tracer_args = tracer_args - - return tracers, tracer_args - - def get_scale(self): - """Return the scale for this source.""" - assert self.current_tracer_args - return self.current_tracer_args.scale - - -class NumberCountsSystematicFactory: - """Factory class for NumberCountsSystematic objects.""" - - @abstractmethod - def create( - self, inferred_zdist: InferredGalaxyZDist - ) -> SourceGalaxySystematic[NumberCountsArgs]: - """Create a NumberCountsSystematic object with the given tracer name.""" - - -class PhotoZShiftFactory(NumberCountsSystematicFactory): - """Factory class for PhotoZShift objects.""" - - def create(self, inferred_zdist: InferredGalaxyZDist) -> PhotoZShift: - """Create a PhotoZShift object with the given tracer name.""" - return PhotoZShift(inferred_zdist.bin_name) - - -class LinearBiasSystematicFactory(NumberCountsSystematicFactory): - """Factory class for LinearBiasSystematic objects.""" - - def create(self, inferred_zdist: InferredGalaxyZDist) -> LinearBiasSystematic: - """Create a LinearBiasSystematic object with the given tracer name.""" - return LinearBiasSystematic(inferred_zdist.bin_name) - - -class PTNonLinearBiasSystematicFactory(NumberCountsSystematicFactory): - """Factory class for PTNonLinearBiasSystematic objects.""" - - def create(self, inferred_zdist: InferredGalaxyZDist) -> PTNonLinearBiasSystematic: - """Create a PTNonLinearBiasSystematic object with the given tracer name.""" - return PTNonLinearBiasSystematic(inferred_zdist.bin_name) - - -class MagnificationBiasSystematicFactory(NumberCountsSystematicFactory): - """Factory class for MagnificationBiasSystematic objects.""" - - def create( - self, inferred_zdist: InferredGalaxyZDist - ) -> MagnificationBiasSystematic: - """Create a MagnificationBiasSystematic object with the given tracer name.""" - return MagnificationBiasSystematic(inferred_zdist.bin_name) - - -class ConstantMagnificationBiasSystematicFactory(NumberCountsSystematicFactory): - """Factory class for ConstantMagnificationBiasSystematic objects.""" - - def create( - self, inferred_zdist: InferredGalaxyZDist - ) -> ConstantMagnificationBiasSystematic: - """Create a ConstantMagnificationBiasSystematic object. - - Use the inferred_zdist to create the systematic. - """ - return ConstantMagnificationBiasSystematic(inferred_zdist.bin_name) - - -class NumberCountsFactory: - """Factory class for NumberCounts objects.""" - - def __init__( - self, - per_bin_systematics: list[NumberCountsSystematicFactory], - global_systematics: Sequence[NumberCountsSystematic], - ) -> None: - """Initialize the NumberCountsFactory.""" - self.per_bin_systematics: list[NumberCountsSystematicFactory] = ( - per_bin_systematics - ) - self.global_systematics: Sequence[NumberCountsSystematic] = global_systematics - self.cache: dict[int, NumberCounts] = {} - - def create(self, inferred_zdist: InferredGalaxyZDist) -> NumberCounts: - """Create a NumberCounts object with the given tracer name and scale.""" - inferred_zdist_id = id(inferred_zdist) - if inferred_zdist_id in self.cache: - return self.cache[inferred_zdist_id] - - systematics: list[SourceGalaxySystematic[NumberCountsArgs]] = [ - systematic_factory.create(inferred_zdist) - for systematic_factory in self.per_bin_systematics - ] - systematics.extend(self.global_systematics) - - nc = NumberCounts.create_ready(inferred_zdist, systematics=systematics) - self.cache[inferred_zdist_id] = nc - - return nc +warnings.warn( + "This module is deprecated. Use firecrown.likelihood.number_counts instead.", + DeprecationWarning, + stacklevel=2, +) diff --git a/firecrown/likelihood/gauss_family/statistic/source/source.py b/firecrown/likelihood/gauss_family/statistic/source/source.py index 2de1050dd..9f658f20c 100644 --- a/firecrown/likelihood/gauss_family/statistic/source/source.py +++ b/firecrown/likelihood/gauss_family/statistic/source/source.py @@ -1,324 +1,21 @@ -"""Abstract base classes for TwoPoint Statistics sources.""" +"""Deprecated module for two-point statistic sources.""" -from __future__ import annotations +# flake8: noqa -from abc import abstractmethod -from dataclasses import dataclass, replace -from typing import Generic, Sequence, TypeVar, final +import warnings -import numpy as np -import numpy.typing as npt -import pyccl -import pyccl.nl_pt -import sacc -from scipy.interpolate import Akima1DInterpolator +import firecrown.likelihood.source -from firecrown import parameters -from firecrown.modeling_tools import ModelingTools -from firecrown.parameters import ParamsMap -from firecrown.updatable import Updatable, UpdatableCollection +# pylint: disable=unused-import,unused-wildcard-import,wildcard-import +from firecrown.likelihood.source import * +# pylint: enable=unused-import,unused-wildcard-import,wildcard-import -class SourceSystematic(Updatable): - """An abstract systematic class (e.g., shear biases, photo-z shifts, etc.). - This class currently has no methods at all, because the argument types for - the `apply` method of different subclasses are different. - """ +assert not hasattr(firecrown.likelihood.source, "__all__") - def read(self, sacc_data: sacc.Sacc): - """Call to allow this object to read from the appropriate sacc data.""" - - -class Source(Updatable): - """An abstract source class (e.g., a sample of lenses). - - Parameters - ---------- - systematics : list of str, optional - A list of the source-level systematics to apply to the source. The - default of `None` implies no systematics. - """ - - systematics: Sequence[SourceSystematic] - cosmo_hash: None | int - tracers: Sequence[Tracer] - - def __init__(self, sacc_tracer: str) -> None: - """Create a Source object that uses the named tracer. - - :param sacc_tracer: the name of the tracer in the SACC file. This is used - as a prefix for its parameters. - """ - super().__init__(parameter_prefix=sacc_tracer) - self.sacc_tracer = sacc_tracer - - @final - def read(self, sacc_data: sacc.Sacc): - """Read the data for this source from the SACC file.""" - if hasattr(self, "systematics"): - for systematic in self.systematics: - systematic.read(sacc_data) - self._read(sacc_data) - - @abstractmethod - def _read(self, sacc_data: sacc.Sacc): - """Abstract method to read the data for this source from the SACC file.""" - - def _update_source(self, params: ParamsMap): - """Method to update the source from the given ParamsMap. - - Any subclass that needs to do more than update its contained :class:`Updatable` - objects should implement this method. - """ - - @final - def _update(self, params: ParamsMap): - """Implementation of Updatable interface method `_update`. - - This clears the current hash and tracer, and calls the abstract method - `_update_source`, which must be implemented in all subclasses. - """ - self.cosmo_hash = None - self.tracers = [] - self._update_source(params) - - @abstractmethod - def get_scale(self) -> float: - """Abstract method to return the scales for this `Source`.""" - - @abstractmethod - def create_tracers(self, tools: ModelingTools): - """Create tracers for this `Source`, for the given cosmology.""" - - @final - def get_tracers(self, tools: ModelingTools) -> Sequence[Tracer]: - """Return the tracer for the given cosmology. - - This method caches its result, so if called a second time with the same - cosmology, no calculation needs to be done. - """ - ccl_cosmo = tools.get_ccl_cosmology() - - cur_hash = hash(ccl_cosmo) - if hasattr(self, "cosmo_hash") and self.cosmo_hash == cur_hash: - return self.tracers - - self.tracers, _ = self.create_tracers(tools) - self.cosmo_hash = cur_hash - return self.tracers - - -class Tracer: - """Extending the pyccl.Tracer object with additional information. - - Bundles together a pyccl.Tracer object with optional information about the - underlying 3D field, a pyccl.nl_pt.PTTracer, and halo profiles. - """ - - @staticmethod - def determine_field_name(field: None | str, tracer: None | str) -> str: - """Gets a field name for a tracer. - - This function encapsulates the policy for determining the value to be - assigned to the :attr:`field` attribute of a :class:`Tracer`. - - It is a static method only to keep it grouped with the class for which it is - defining the initialization policy. - """ - if field is not None: - return field - if tracer is not None: - return tracer - return "delta_matter" - - def __init__( - self, - tracer: pyccl.Tracer, - tracer_name: None | str = None, - field: None | str = None, - pt_tracer: None | pyccl.nl_pt.PTTracer = None, - halo_profile: None | pyccl.halos.HaloProfile = None, - halo_2pt: None | pyccl.halos.Profile2pt = None, - ): - """Initialize a new Tracer based on the pyccl.Tracer which must not be None. - - Note that the :class:`pyccl.Tracer` is not copied; we store a reference to the - original tracer. Be careful not to accidentally share :class:`pyccl.Tracer`s. - - If no tracer_name is supplied, then the tracer_name is set to the name of the - :class:`pyccl.Tracer` class that was used. - - If no `field` is given, then the attribute :attr:`field` is set to either - (1) the tracer_name, if one was given, or (2) 'delta_matter'. - """ - assert tracer is not None - self.ccl_tracer = tracer - self.tracer_name: str = tracer_name or tracer.__class__.__name__ - self.field = Tracer.determine_field_name(field, tracer_name) - self.pt_tracer = pt_tracer - self.halo_profile = halo_profile - self.halo_2pt = halo_2pt - - @property - def has_pt(self) -> bool: - """Return True if we have a pt_tracer, and False if not.""" - return self.pt_tracer is not None - - @property - def has_hm(self) -> bool: - """Return True if we have a halo_profile, and False if not.""" - return self.halo_profile is not None - - -# Sources of galaxy distributions - - -@dataclass(frozen=True) -class SourceGalaxyArgs: - """Class for galaxy based sources arguments.""" - - z: npt.NDArray[np.float64] - dndz: npt.NDArray[np.float64] - - scale: float = 1.0 - - field: str = "delta_matter" - - -_SourceGalaxyArgsT = TypeVar("_SourceGalaxyArgsT", bound=SourceGalaxyArgs) - - -class SourceGalaxySystematic(SourceSystematic, Generic[_SourceGalaxyArgsT]): - """Abstract base class for all galaxy based source systematics.""" - - @abstractmethod - def apply( - self, tools: ModelingTools, tracer_arg: _SourceGalaxyArgsT - ) -> _SourceGalaxyArgsT: - """Apply method to include systematics in the tracer_arg.""" - - -_SourceGalaxySystematicT = TypeVar( - "_SourceGalaxySystematicT", bound=SourceGalaxySystematic +warnings.warn( + "This module is deprecated. Use firecrown.likelihood.source instead.", + DeprecationWarning, + stacklevel=2, ) - - -SOURCE_GALAXY_SYSTEMATIC_DEFAULT_DELTA_Z = 0.0 - - -class SourceGalaxyPhotoZShift( - SourceGalaxySystematic[_SourceGalaxyArgsT], Generic[_SourceGalaxyArgsT] -): - """A photo-z shift bias. - - This systematic shifts the photo-z distribution by some amount `delta_z`. - - The following parameters are special Updatable parameters, which means that - they can be updated by the sampler, sacc_tracer is going to be used as a - prefix for the parameters: - - :ivar delta_z: the photo-z shift. - """ - - def __init__(self, sacc_tracer: str) -> None: - """Create a PhotoZShift object, using the specified tracer name. - - :param sacc_tracer: the name of the tracer in the SACC file. This is used - as a prefix for its parameters. - """ - super().__init__(parameter_prefix=sacc_tracer) - - self.delta_z = parameters.register_new_updatable_parameter( - default_value=SOURCE_GALAXY_SYSTEMATIC_DEFAULT_DELTA_Z - ) - - def apply(self, tools: ModelingTools, tracer_arg: _SourceGalaxyArgsT): - """Apply a shift to the photo-z distribution of a source.""" - dndz_interp = Akima1DInterpolator(tracer_arg.z, tracer_arg.dndz) - - dndz = dndz_interp(tracer_arg.z - self.delta_z, extrapolate=False) - dndz[np.isnan(dndz)] = 0.0 - - return replace( - tracer_arg, - dndz=dndz, - ) - - -class SourceGalaxySelectField( - SourceGalaxySystematic[_SourceGalaxyArgsT], Generic[_SourceGalaxyArgsT] -): - """The source galaxy select field systematic. - - A systematic that allows specifying the 3D field that will be used - to select the 3D power spectrum when computing the angular power - spectrum. - """ - - def __init__(self, field: str = "delta_matter"): - """Specify which 3D field should be used when computing angular power spectra. - - :param field: the name of the 3D field that is associated to the tracer. - Default: `"delta_matter"` - """ - super().__init__() - self.field = field - - def apply( - self, tools: ModelingTools, tracer_arg: _SourceGalaxyArgsT - ) -> _SourceGalaxyArgsT: - """Apply method to include systematics in the tracer_arg.""" - return replace(tracer_arg, field=self.field) - - -class SourceGalaxy(Source, Generic[_SourceGalaxyArgsT]): - """Source class for galaxy based sources.""" - - def __init__( - self, - *, - sacc_tracer: str, - systematics: None | Sequence[SourceGalaxySystematic] = None, - ): - """Initialize the SourceGalaxy object. - - :param sacc_tracer: the name of the tracer in the SACC file. This is used - as a prefix for its parameters. - - """ - super().__init__(sacc_tracer) - - self.sacc_tracer = sacc_tracer - self.current_tracer_args: None | _SourceGalaxyArgsT = None - self.systematics: UpdatableCollection[SourceGalaxySystematic] = ( - UpdatableCollection(systematics) - ) - self.tracer_args: _SourceGalaxyArgsT - - def _read(self, sacc_data: sacc.Sacc): - """Read the galaxy redshift distribution model from a sacc file. - - All derived classes must call this method in their own `_read` method - after they have read their own data and initialized their tracer_args. - """ - try: - tracer_args = self.tracer_args - except AttributeError as exc: - raise RuntimeError( - "Must initialize tracer_args before calling _read on SourceGalaxy" - ) from exc - - tracer = sacc_data.get_tracer(self.sacc_tracer) - - z = getattr(tracer, "z").copy().flatten() - nz = getattr(tracer, "nz").copy().flatten() - indices = np.argsort(z) - z = z[indices] - nz = nz[indices] - - self.tracer_args = replace( - tracer_args, - z=z, - dndz=nz, - ) diff --git a/firecrown/likelihood/gauss_family/statistic/source/weak_lensing.py b/firecrown/likelihood/gauss_family/statistic/source/weak_lensing.py index 529467639..1be4e14e1 100644 --- a/firecrown/likelihood/gauss_family/statistic/source/weak_lensing.py +++ b/firecrown/likelihood/gauss_family/statistic/source/weak_lensing.py @@ -1,405 +1,20 @@ -"""Weak lensing source and systematics.""" +"""Deprecated module for weak lensing likelihoods.""" -from __future__ import annotations +# flake8: noqa -from abc import abstractmethod -from dataclasses import dataclass, replace -from typing import Sequence, final +import warnings -import numpy as np -import numpy.typing as npt -import pyccl -import pyccl.nl_pt -import sacc - -from firecrown import parameters -from firecrown.likelihood.gauss_family.statistic.source.source import ( - SourceGalaxy, - SourceGalaxyArgs, - SourceGalaxyPhotoZShift, - SourceGalaxySelectField, - SourceGalaxySystematic, - Tracer, -) -from firecrown.metadata.two_point import InferredGalaxyZDist -from firecrown.modeling_tools import ModelingTools -from firecrown.parameters import ParamsMap - -__all__ = ["WeakLensing"] - - -@dataclass(frozen=True) -class WeakLensingArgs(SourceGalaxyArgs): - """Class for weak lensing tracer builder argument.""" - - ia_bias: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None - - has_pt: bool = False - has_hm: bool = False - - ia_pt_c_1: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None - ia_pt_c_d: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None - ia_pt_c_2: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None - - -class WeakLensingSystematic(SourceGalaxySystematic[WeakLensingArgs]): - """Abstract base class for all weak lensing systematics.""" - - @abstractmethod - def apply( - self, tools: ModelingTools, tracer_arg: WeakLensingArgs - ) -> WeakLensingArgs: - """Apply method to include systematics in the tracer_arg.""" - - -class PhotoZShift(SourceGalaxyPhotoZShift[WeakLensingArgs]): - """Photo-z shift systematic.""" - - -class SelectField(SourceGalaxySelectField[WeakLensingArgs]): - """Systematic to select 3D field.""" - - -MULTIPLICATIVE_SHEAR_BIAS_DEFAULT_BIAS = 1.0 - - -class MultiplicativeShearBias(WeakLensingSystematic): - """Multiplicative shear bias systematic. - - This systematic adjusts the `scale_` of a source by `(1 + m)`. - - The following parameters are special Updatable parameters, which means that - they can be updated by the sampler, sacc_tracer is going to be used as a - prefix for the parameters: - - :ivar mult_bias: the multiplicative shear bias parameter. - """ - - def __init__(self, sacc_tracer: str) -> None: - """Create a MultiplicativeShearBias object that uses the named tracer. - - :param sacc_tracer: the name of the tracer in the SACC file. This is used - as a prefix for its parameters. - """ - super().__init__(parameter_prefix=sacc_tracer) - - self.mult_bias = parameters.register_new_updatable_parameter( - default_value=MULTIPLICATIVE_SHEAR_BIAS_DEFAULT_BIAS - ) - - def apply( - self, tools: ModelingTools, tracer_arg: WeakLensingArgs - ) -> WeakLensingArgs: - """Apply multiplicative shear bias to a source. - - The `scale_` of the source is multiplied by `(1 + m)`. - - :param tools: A ModelingTools object. - :param tracer_arg: The WeakLensingArgs to which apply the shear bias. - - :returns: A new WeakLensingArgs object with the shear bias applied. - """ - return replace( - tracer_arg, - scale=tracer_arg.scale * (1.0 + self.mult_bias), - ) - - -LINEAR_ALIGNMENT_DEFAULT_IA_BIAS = 0.5 -LINEAR_ALIGNMENT_DEFAULT_ALPHAZ = 0.0 -LINEAR_ALIGNMENT_DEFAULT_ALPHAG = 1.0 -LINEAR_ALIGNMENT_DEFAULT_Z_PIV = 0.5 - - -class LinearAlignmentSystematic(WeakLensingSystematic): - """Linear alignment systematic. - - This systematic adds a linear intrinsic alignment model systematic - which varies with redshift and the growth function. - - The following parameters are special Updatable parameters, which means that - they can be updated by the sampler, sacc_tracer is going to be used as a - prefix for the parameters: - - :ivar ia_bias: the intrinsic alignment bias parameter. - :ivar alphaz: the redshift dependence of the intrinsic alignment bias. - :ivar alphag: the growth function dependence of the intrinsic alignment bias. - :ivar z_piv: the pivot redshift for the intrinsic alignment bias. - """ - - def __init__(self, sacc_tracer: None | str = None, alphag=1.0): - """Create a LinearAlignmentSystematic object, using the specified tracer name. - - :param sacc_tracer: the name of the tracer in the SACC file. This is used - as a prefix for its parameters. - - """ - super().__init__(parameter_prefix=sacc_tracer) - - self.ia_bias = parameters.register_new_updatable_parameter( - default_value=LINEAR_ALIGNMENT_DEFAULT_IA_BIAS - ) - self.alphaz = parameters.register_new_updatable_parameter( - default_value=LINEAR_ALIGNMENT_DEFAULT_ALPHAZ - ) - self.alphag = parameters.register_new_updatable_parameter( - alphag, default_value=LINEAR_ALIGNMENT_DEFAULT_ALPHAG - ) - self.z_piv = parameters.register_new_updatable_parameter( - default_value=LINEAR_ALIGNMENT_DEFAULT_Z_PIV - ) - - def apply( - self, tools: ModelingTools, tracer_arg: WeakLensingArgs - ) -> WeakLensingArgs: - """Return a new linear alignment systematic. - - This choice is based on the given tracer_arg, in the context of the given - cosmology. - """ - ccl_cosmo = tools.get_ccl_cosmology() - - pref = ((1.0 + tracer_arg.z) / (1.0 + self.z_piv)) ** self.alphaz - pref *= pyccl.growth_factor(ccl_cosmo, 1.0 / (1.0 + tracer_arg.z)) ** ( - self.alphag - 1.0 - ) - - ia_bias_array = pref * self.ia_bias - - return replace( - tracer_arg, - ia_bias=(tracer_arg.z, ia_bias_array), - ) - - -TATT_ALIGNMENT_DEFAULT_IA_A_1 = 1.0 -TATT_ALIGNMENT_DEFAULT_IA_A_2 = 0.5 -TATT_ALIGNMENT_DEFAULT_IA_A_D = 0.5 +import firecrown.likelihood.weak_lensing +# pylint: disable=unused-import,unused-wildcard-import,wildcard-import +from firecrown.likelihood.weak_lensing import * -class TattAlignmentSystematic(WeakLensingSystematic): - """TATT alignment systematic. +# pylint: enable=unused-import,unused-wildcard-import,wildcard-import - This systematic adds a TATT (nonlinear) intrinsic alignment model systematic. +assert not hasattr(firecrown.likelihood.weak_lensing, "__all__") - The following parameters are special Updatable parameters, which means that - they can be updated by the sampler, sacc_tracer is going to be used as a - prefix for the parameters: - - :ivar ia_a_1: the amplitude of the linear alignment model. - :ivar ia_a_2: the amplitude of the quadratic alignment model. - :ivar ia_a_d: the amplitude of the density-dependent alignment model. - """ - - def __init__(self, sacc_tracer: None | str = None): - """Create a TattAlignmentSystematic object, using the specified tracer name. - - :param sacc_tracer: the name of the tracer in the SACC file. This is used - as a prefix for its parameters. - """ - super().__init__(parameter_prefix=sacc_tracer) - self.ia_a_1 = parameters.register_new_updatable_parameter( - default_value=TATT_ALIGNMENT_DEFAULT_IA_A_1 - ) - self.ia_a_2 = parameters.register_new_updatable_parameter( - default_value=TATT_ALIGNMENT_DEFAULT_IA_A_2 - ) - self.ia_a_d = parameters.register_new_updatable_parameter( - default_value=TATT_ALIGNMENT_DEFAULT_IA_A_D - ) - - def apply( - self, tools: ModelingTools, tracer_arg: WeakLensingArgs - ) -> WeakLensingArgs: - """Return a new linear alignment systematic. - - This choice is based on the given tracer_arg, in the context of the given - cosmology. - """ - ccl_cosmo = tools.get_ccl_cosmology() - z = tracer_arg.z - c_1, c_d, c_2 = pyccl.nl_pt.translate_IA_norm( - ccl_cosmo, - z=z, - a1=self.ia_a_1, - a1delta=self.ia_a_d, - a2=self.ia_a_2, - Om_m2_for_c2=False, - ) - - return replace( - tracer_arg, - has_pt=True, - ia_pt_c_1=(z, c_1), - ia_pt_c_d=(z, c_d), - ia_pt_c_2=(z, c_2), - ) - - -class WeakLensing(SourceGalaxy[WeakLensingArgs]): - """Source class for weak lensing.""" - - def __init__( - self, - *, - sacc_tracer: str, - scale: float = 1.0, - systematics: None | Sequence[SourceGalaxySystematic[WeakLensingArgs]] = None, - ): - """Initialize the WeakLensing object. - - :param sacc_tracer: the name of the tracer in the SACC file. This is used - as a prefix for its parameters. - :param scale: the scale of the source. This is used to scale the shear - power spectrum. - :param systematics: a list of WeakLensingSystematic objects to apply to - this source. - - """ - super().__init__(sacc_tracer=sacc_tracer, systematics=systematics) - - self.sacc_tracer = sacc_tracer - self.scale = scale - self.current_tracer_args: None | WeakLensingArgs = None - self.tracer_args: WeakLensingArgs - - @classmethod - def create_ready( - cls, - inferred_zdist: InferredGalaxyZDist, - systematics: None | list[SourceGalaxySystematic[WeakLensingArgs]] = None, - ) -> WeakLensing: - """Create a WeakLensing object with the given tracer name and scale.""" - obj = cls(sacc_tracer=inferred_zdist.bin_name, systematics=systematics) - obj.tracer_args = WeakLensingArgs( - scale=obj.scale, z=inferred_zdist.z, dndz=inferred_zdist.dndz, ia_bias=None - ) - - return obj - - @final - def _update_source(self, params: ParamsMap): - """Implementation of Source interface `_update_source`. - - This updates all the contained systematics. - """ - self.systematics.update(params) - - def _read(self, sacc_data: sacc.Sacc) -> None: - """Read the data for this source from the SACC file. - - This sets self.tracer_args, based on the data in `sacc_data` associated with - this object's `sacc_tracer` name. - """ - self.tracer_args = WeakLensingArgs( - scale=self.scale, z=np.array([]), dndz=np.array([]), ia_bias=None - ) - - super()._read(sacc_data) - - def create_tracers(self, tools: ModelingTools): - """Render a source by applying systematics.""" - ccl_cosmo = tools.get_ccl_cosmology() - tracer_args = self.tracer_args - - assert self.systematics is not None - for systematic in self.systematics: - tracer_args = systematic.apply(tools, tracer_args) - - ccl_wl_tracer = pyccl.WeakLensingTracer( - ccl_cosmo, - dndz=(tracer_args.z, tracer_args.dndz), - ia_bias=tracer_args.ia_bias, - ) - tracers = [Tracer(ccl_wl_tracer, tracer_name="shear", field=tracer_args.field)] - - if tracer_args.has_pt: - ia_pt_tracer = pyccl.nl_pt.PTIntrinsicAlignmentTracer( - c1=tracer_args.ia_pt_c_1, - cdelta=tracer_args.ia_pt_c_d, - c2=tracer_args.ia_pt_c_2, - ) - - ccl_wl_dummy_tracer = pyccl.WeakLensingTracer( - ccl_cosmo, - has_shear=False, - use_A_ia=False, - dndz=(tracer_args.z, tracer_args.dndz), - ia_bias=(tracer_args.z, np.ones_like(tracer_args.z)), - ) - ia_tracer = Tracer( - ccl_wl_dummy_tracer, tracer_name="intrinsic_pt", pt_tracer=ia_pt_tracer - ) - tracers.append(ia_tracer) - - self.current_tracer_args = tracer_args - - return tracers, tracer_args - - def get_scale(self): - """Returns the scales for this Source.""" - assert self.current_tracer_args - return self.current_tracer_args.scale - - -class WeakLensingSystematicFactory: - """Factory class for WeakLensingSystematic objects.""" - - @abstractmethod - def create( - self, inferred_zdist: InferredGalaxyZDist - ) -> SourceGalaxySystematic[WeakLensingArgs]: - """Create a WeakLensingSystematic object with the given tracer name.""" - - -class MultiplicativeShearBiasFactory(WeakLensingSystematicFactory): - """Factory class for MultiplicativeShearBias objects.""" - - def create(self, inferred_zdist: InferredGalaxyZDist) -> MultiplicativeShearBias: - return MultiplicativeShearBias(inferred_zdist.bin_name) - - -class TattAlignmentSystematicFactory(WeakLensingSystematicFactory): - """Factory class for TattAlignmentSystematic objects.""" - - def create(self, inferred_zdist: InferredGalaxyZDist) -> TattAlignmentSystematic: - return TattAlignmentSystematic(inferred_zdist.bin_name) - - -class PhotoZShiftFactory(WeakLensingSystematicFactory): - """Factory class for PhotoZShift objects.""" - - def create(self, inferred_zdist: InferredGalaxyZDist) -> PhotoZShift: - return PhotoZShift(inferred_zdist.bin_name) - - -class WeakLensingFactory: - """Factory class for WeakLensing objects.""" - - def __init__( - self, - per_bin_systematics: list[WeakLensingSystematicFactory], - global_systematics: Sequence[WeakLensingSystematic], - ) -> None: - self.per_bin_systematics: list[WeakLensingSystematicFactory] = ( - per_bin_systematics - ) - self.global_systematics: Sequence[WeakLensingSystematic] = global_systematics - self.cache: dict[int, WeakLensing] = {} - - def create(self, inferred_zdist: InferredGalaxyZDist) -> WeakLensing: - """Create a WeakLensing object with the given tracer name and scale.""" - inferred_zdist_id = id(inferred_zdist) - if inferred_zdist_id in self.cache: - return self.cache[inferred_zdist_id] - - systematics: list[SourceGalaxySystematic[WeakLensingArgs]] = [ - systematic_factory.create(inferred_zdist) - for systematic_factory in self.per_bin_systematics - ] - systematics.extend(self.global_systematics) - - wl = WeakLensing.create_ready(inferred_zdist, systematics) - self.cache[inferred_zdist_id] = wl - - return wl +warnings.warn( + "This module is deprecated. Use firecrown.likelihood.weak_lensing instead.", + DeprecationWarning, + stacklevel=2, +) diff --git a/firecrown/likelihood/gauss_family/statistic/statistic.py b/firecrown/likelihood/gauss_family/statistic/statistic.py index e951ff552..5c1c68221 100644 --- a/firecrown/likelihood/gauss_family/statistic/statistic.py +++ b/firecrown/likelihood/gauss_family/statistic/statistic.py @@ -1,292 +1,21 @@ -"""Gaussian Family Statistic Module. +""""Deprecated module with classes related to Statistic.""" -The Statistic class describing objects that implement methods to compute the -data and theory vectors for a :class:`GaussFamily` subclass. -""" +# flake8: noqa -from __future__ import annotations -from typing import final -from dataclasses import dataclass -from abc import abstractmethod import warnings -import numpy as np -import numpy.typing as npt -import sacc -import firecrown.parameters -from firecrown.parameters import DerivedParameterCollection, RequiredParameters -from firecrown.modeling_tools import ModelingTools -from firecrown.updatable import Updatable +import firecrown.likelihood.statistic +# pylint: disable=unused-import,unused-wildcard-import,wildcard-import +from firecrown.likelihood.statistic import * -class DataVector(npt.NDArray[np.float64]): - """This class wraps a np.ndarray that represents some observed data values.""" +# pylint: enable=unused-import,unused-wildcard-import,wildcard-import - @classmethod - def create(cls, vals: npt.NDArray[np.float64]) -> DataVector: - """Create a DataVector that wraps a copy of the given array vals.""" - return vals.view(cls) - @classmethod - def from_list(cls, vals: list[float]) -> DataVector: - """Create a DataVector from the given list of floats.""" - array = np.array(vals) - return cls.create(array) +assert not hasattr(firecrown.likelihood.statistic, "__all__") - -class TheoryVector(npt.NDArray[np.float64]): - """This class represents an observation predicted by some theory.""" - - @classmethod - def create(cls, vals: npt.NDArray[np.float64]) -> TheoryVector: - """Create a TheoryVector that wraps a copy of the given array vals.""" - return vals.view(cls) - - @classmethod - def from_list(cls, vals: list[float]) -> TheoryVector: - """Create a TheoryVector from the given list of floats.""" - array = np.array(vals) - return cls.create(array) - - -def residuals(data: DataVector, theory: TheoryVector) -> npt.NDArray[np.float64]: - """Return a bare np.ndarray with the difference between `data` and `theory`. - - This is to be preferred to using arithmetic on the vectors directly. - """ - assert isinstance(data, DataVector) - assert isinstance(theory, TheoryVector) - return (data - theory).view(np.ndarray) - - -@dataclass -class StatisticsResult: - """This is the type returned by the `compute` method of any `Statistic`.""" - - data: DataVector - theory: TheoryVector - - def __post_init__(self): - """Make sure the data and theory vectors are of the same shape.""" - assert self.data.shape == self.theory.shape - - def residuals(self) -> npt.NDArray[np.float64]: - """Return the residuals -- the difference between data and theory.""" - return self.data - self.theory - - def __iter__(self): - """Iterate through the data members. - - This is to allow automatic unpacking, as if the StatisticsResult were a tuple - of (data, theory). - - This method is a temporary measure to help code migrate to the newer, - safer interface for Statistic.compute(). - """ - warnings.warn( - "Iteration and tuple unpacking for StatisticsResult is " - "deprecated.\nPlease use the StatisticsResult class accessors" - ".data and .theory by name." - ) - yield self.data - yield self.theory - - -class StatisticUnreadError(RuntimeError): - """Error raised when accessing an un-read statistic. - - Run-time error indicating an attempt has been made to use a statistic - that has not had `read` called in it. - """ - - def __init__(self, stat: Statistic): - msg = ( - f"The statistic {stat} was used for calculation before `read` " - f"was called.\nIt may be that a likelihood factory function did not" - f"call `read` before returning the likelihood." - ) - super().__init__(msg) - self.statstic = stat - - -class Statistic(Updatable): - """The abstract base class for all physics-related statistics. - - Statistics read data from a SACC object as part of a multi-phase - initialization. They manage a :class:`DataVector` and, given a - :class:`ModelingTools` object, can compute a :class:`TheoryVector`. - - Statistics represent things like two-point functions and mass functions. - """ - - def __init__(self, parameter_prefix: None | str = None): - super().__init__(parameter_prefix=parameter_prefix) - self.sacc_indices: None | npt.NDArray[np.int64] - self.ready = False - self.computed_theory_vector = False - self.theory_vector: None | TheoryVector = None - - def read(self, _: sacc.Sacc) -> None: - """Read the data for this statistic and mark it as ready for use. - - Derived classes that override this function should make sure to call the base - class method using: - - .. code-block:: python - - super().read(sacc_data) - - as the last thing they do in `__init__`. - - Note that currently the argument is not used; it is present so that this - method will have the correct argument type for the override. - """ - self.ready = True - if len(self.get_data_vector()) == 0: - raise RuntimeError( - f"the statistic {self} has read a data vector " - f"of length 0; the length must be positive" - ) - - def _reset(self): - """Reset this statistic. - - All subclasses implementations must call super()._reset() - """ - self.computed_theory_vector = False - self.theory_vector = None - - @abstractmethod - def get_data_vector(self) -> DataVector: - """Gets the statistic data vector.""" - - @final - def compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: - """Compute a statistic from sources, applying any systematics.""" - if not self.is_updated(): - raise RuntimeError( - f"The statistic {self} has not been updated with parameters." - ) - self.theory_vector = self._compute_theory_vector(tools) - self.computed_theory_vector = True - - return self.theory_vector - - @abstractmethod - def _compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: - """Compute a statistic from sources, concrete implementation.""" - - def get_theory_vector(self) -> TheoryVector: - """Returns the last computed theory vector. - - Raises a RuntimeError if the vector has not been computed. - """ - if not self.computed_theory_vector: - raise RuntimeError( - f"The theory for statistic {self} has not been computed yet." - ) - assert self.theory_vector is not None, ( - "implementation error, " - "computed_theory_vector is True but theory_vector is None" - ) - return self.theory_vector - - -class GuardedStatistic(Updatable): - """An internal class used to maintain state on statistics. - - :class:`GuardedStatistic` is used by the framework to maintain and - validate the state of instances of classes derived from :class:`Statistic`. - """ - - def __init__(self, stat: Statistic): - """Initialize the GuardedStatistic to contain the given :class:`Statistic`.""" - super().__init__() - assert isinstance(stat, Statistic) - self.statistic = stat - - def read(self, sacc_data: sacc.Sacc) -> None: - """Read whatever data is needed from the given :class:`sacc.Sacc` object. - - After this function is called, the object should be prepared for the - calling of the methods :meth:`get_data_vector` and - :meth:`compute_theory_vector`. - """ - if self.statistic.ready: - raise RuntimeError("Firecrown has called read twice on a GuardedStatistic") - try: - self.statistic.read(sacc_data) - except TypeError as exc: - msg = ( - f"A statistic of type {type(self.statistic).__name__} has raised " - f"an exception during `read`.\nThe problem may be a malformed " - f"SACC data object." - ) - raise RuntimeError(msg) from exc - - def get_data_vector(self) -> DataVector: - """Return the contained :class:`Statistic`'s data vector. - - :class:`GuardedStatistic` ensures that :meth:`read` has been called. - first. - """ - if not self.statistic.ready: - raise StatisticUnreadError(self.statistic) - return self.statistic.get_data_vector() - - def compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: - """Return the contained :class:`Statistic`'s computed theory vector. - - :class:`GuardedStatistic` ensures that :meth:`read` has been called. - first. - """ - if not self.statistic.ready: - raise StatisticUnreadError(self.statistic) - return self.statistic.compute_theory_vector(tools) - - -class TrivialStatistic(Statistic): - """A minimal statistic only to be used for testing Gaussian likelihoods. - - It returns a :class:`DataVector` and :class:`TheoryVector` each of which is - three elements long. The SACC data provided to :meth:`TrivialStatistic.read` - must supply the necessary values. - """ - - def __init__(self) -> None: - """Initialize this statistic.""" - super().__init__() - # Data and theory will both be of length self.count - self.count = 3 - self.data_vector: None | DataVector = None - self.mean = firecrown.parameters.register_new_updatable_parameter( - default_value=0.0 - ) - self.computed_theory_vector = False - - def read(self, sacc_data: sacc.Sacc): - """Read the necessary items from the sacc data.""" - our_data = sacc_data.get_mean(data_type="count") - assert len(our_data) == self.count - self.data_vector = DataVector.from_list(our_data) - self.sacc_indices = np.arange(len(self.data_vector)) - super().read(sacc_data) - - @final - def _required_parameters(self) -> RequiredParameters: - """Return an empty RequiredParameters.""" - return RequiredParameters([]) - - @final - def _get_derived_parameters(self) -> DerivedParameterCollection: - """Return an empty DerivedParameterCollection.""" - return DerivedParameterCollection([]) - - def get_data_vector(self) -> DataVector: - """Return the data vector; raise exception if there is none.""" - assert self.data_vector is not None - return self.data_vector - - def _compute_theory_vector(self, _: ModelingTools) -> TheoryVector: - """Return a fixed theory vector.""" - return TheoryVector.from_list([self.mean] * self.count) +warnings.warn( + "This module is deprecated. Use firecrown.likelihood.statistic instead.", + DeprecationWarning, + stacklevel=2, +) diff --git a/firecrown/likelihood/gauss_family/statistic/supernova.py b/firecrown/likelihood/gauss_family/statistic/supernova.py index 83253f152..b18acd9b2 100644 --- a/firecrown/likelihood/gauss_family/statistic/supernova.py +++ b/firecrown/likelihood/gauss_family/statistic/supernova.py @@ -1,76 +1,21 @@ -"""Supernova statistic support.""" +"""Deprecated module with classes related to the Supernova statist.""" -from __future__ import annotations +# flake8: noqa -import numpy as np -import numpy.typing as npt -import pyccl -import sacc -from sacc.tracers import MiscTracer +import warnings -from firecrown import parameters -from firecrown.likelihood.gauss_family.statistic.statistic import ( - DataVector, - Statistic, - TheoryVector, -) -from firecrown.modeling_tools import ModelingTools - -SNIA_DEFAULT_M = -19.2 - - -class Supernova(Statistic): - """A supernova statistic. +import firecrown.likelihood.supernova - This statistic that applies an additive shift M to a supernova's distance - modulus. - """ +# pylint: disable=unused-import,unused-wildcard-import,wildcard-import +from firecrown.likelihood.supernova import * - def __init__(self, sacc_tracer: str) -> None: - """Initialize this statistic.""" - super().__init__(parameter_prefix=sacc_tracer) +# pylint: enable=unused-import,unused-wildcard-import,wildcard-import - self.sacc_tracer = sacc_tracer - self.data_vector: None | DataVector = None - self.a: None | npt.NDArray[np.float64] = None - self.M = parameters.register_new_updatable_parameter( - default_value=SNIA_DEFAULT_M - ) - def read(self, sacc_data: sacc.Sacc) -> None: - """Read the data for this statistic from the SACC file.""" - # We do not actually need the tracer, but we want to make sure the SACC - # data contains this tracer. - # TODO: remove the work-around when the new version of SACC supporting - # sacc.Sacc.has_tracer is available. - try: - tracer = sacc_data.get_tracer(self.sacc_tracer) - except KeyError as exc: - # Translate to the error type we want - raise ValueError( - f"The SACC file does not contain the MiscTracer {self.sacc_tracer}" - ) from exc - if not isinstance(tracer, MiscTracer): - raise ValueError( - f"The SACC tracer {self.sacc_tracer} is not a " f"MiscTracer" - ) +assert not hasattr(firecrown.likelihood.supernova, "__all__") - data_points = sacc_data.get_data_points( - data_type="supernova_distance_mu", tracers=(self.sacc_tracer,) - ) - z = np.array([dp.get_tag("z") for dp in data_points]) - self.a = 1.0 / (1.0 + z) - self.data_vector = DataVector.from_list([dp.value for dp in data_points]) - self.sacc_indices = np.arange(len(self.data_vector)) - super().read(sacc_data) - - def get_data_vector(self) -> DataVector: - """Return the data vector; raise exception if there is none.""" - assert self.data_vector is not None - return self.data_vector - - def _compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: - """Compute SNIa distance statistic using CCL.""" - ccl_cosmo = tools.get_ccl_cosmology() - prediction = self.M + pyccl.distance_modulus(ccl_cosmo, self.a) - return TheoryVector.create(prediction) +warnings.warn( + "This module is deprecated. Use firecrown.likelihood.supernova instead.", + DeprecationWarning, + stacklevel=2, +) diff --git a/firecrown/likelihood/gauss_family/statistic/two_point.py b/firecrown/likelihood/gauss_family/statistic/two_point.py index 1d050bc7a..50a82bc0d 100644 --- a/firecrown/likelihood/gauss_family/statistic/two_point.py +++ b/firecrown/likelihood/gauss_family/statistic/two_point.py @@ -1,776 +1,20 @@ -"""Two point statistic support.""" +"""Deprecated module with classes related to the TwoPoint statistic.""" -from __future__ import annotations +# flake8: noqa -import copy import warnings -from typing import Sequence, TypedDict -import numpy as np -import numpy.typing as npt -import pyccl -import pyccl.nl_pt -import sacc.windows -import scipy.interpolate +import firecrown.likelihood.two_point -from firecrown.likelihood.gauss_family.statistic.source.source import Source, Tracer -from firecrown.likelihood.gauss_family.statistic.source.weak_lensing import ( - WeakLensingFactory, - WeakLensing, -) -from firecrown.likelihood.gauss_family.statistic.source.number_counts import ( - NumberCountsFactory, - NumberCounts, -) -from firecrown.likelihood.gauss_family.statistic.statistic import ( - DataVector, - Statistic, - TheoryVector, -) -from firecrown.metadata.two_point import ( - GalaxyMeasuredType, - InferredGalaxyZDist, - TRACER_NAMES_TOTAL, - TracerNames, - TwoPointCells, - Window, - extract_window_function, -) -from firecrown.modeling_tools import ModelingTools -from firecrown.updatable import UpdatableCollection - -# only supported types are here, anything else will throw -# a value error -SACC_DATA_TYPE_TO_CCL_KIND = { - "galaxy_density_cl": "cl", - "galaxy_density_xi": "NN", - "galaxy_shearDensity_cl_e": "cl", - "galaxy_shearDensity_xi_t": "NG", - "galaxy_shear_cl_ee": "cl", - "galaxy_shear_xi_minus": "GG-", - "galaxy_shear_xi_plus": "GG+", - "cmbGalaxy_convergenceDensity_xi": "NN", - "cmbGalaxy_convergenceShear_xi_t": "NG", -} - -ELL_FOR_XI_DEFAULTS = {"minimum": 2, "midpoint": 50, "maximum": 60_000, "n_log": 200} - - -def _ell_for_xi( - *, minimum: int, midpoint: int, maximum: int, n_log: int -) -> npt.NDArray[np.int64]: - """Create an array of ells to sample the power spectrum. - - This is used for for real-space predictions. The result will contain - each integral value from min to mid. Starting from mid, and going up - to max, there will be n_log logarithmically spaced values. - - All values are rounded to the nearest integer. - """ - assert minimum >= 0 - assert minimum < midpoint - assert midpoint < maximum - lower_range = np.linspace(minimum, midpoint - 1, midpoint - minimum) - upper_range = np.logspace(np.log10(midpoint), np.log10(maximum), n_log) - concatenated = np.concatenate((lower_range, upper_range)) - # Round the results to the nearest integer values. - # N.B. the dtype of the result is np.dtype[float64] - return np.unique(np.around(concatenated)).astype(np.int64) - - -def _generate_ell_or_theta(*, minimum, maximum, n, binning="log"): - if binning == "log": - edges = np.logspace(np.log10(minimum), np.log10(maximum), n + 1) - return np.sqrt(edges[1:] * edges[:-1]) - edges = np.linspace(minimum, maximum, n + 1) - return (edges[1:] + edges[:-1]) / 2.0 - - -# @functools.lru_cache(maxsize=128) -def _cached_angular_cl(cosmo, tracers, ells, p_of_k_a=None): - return pyccl.angular_cl( - cosmo, tracers[0], tracers[1], np.array(ells), p_of_k_a=p_of_k_a - ) - - -def make_log_interpolator(x, y): - """Return a function object that does 1D spline interpolation. - - If all the y values are greater than 0, the function - interpolates log(y) as a function of log(x). - Otherwise, the function interpolates y as a function of log(x). - The resulting interpolater will not extrapolate; if called with - an out-of-range argument it will raise a ValueError. - """ - # TODO: There is no code in Firecrown, neither test nor example, that uses - # this in any way. - if np.all(y > 0): - # use log-log interpolation - intp = scipy.interpolate.InterpolatedUnivariateSpline( - np.log(x), np.log(y), ext=2 - ) - return lambda x_, intp=intp: np.exp(intp(np.log(x_))) - # only use log for x - intp = scipy.interpolate.InterpolatedUnivariateSpline(np.log(x), y, ext=2) - return lambda x_, intp=intp: intp(np.log(x_)) - - -def calculate_ells_for_interpolation(w: Window) -> npt.NDArray[np.int64]: - """See _ell_for_xi. - - This method mixes together: - 1. the default parameters in ELL_FOR_XI_DEFAULTS - 2. the first and last values in w. - - and then calls _ell_for_xi with those arguments, returning whatever it - returns. - """ - ell_config = { - **ELL_FOR_XI_DEFAULTS, - "maximum": w.ells[-1], - } - ell_config["minimum"] = max(ell_config["minimum"], w.ells[0]) - return _ell_for_xi(**ell_config) - - -class EllOrThetaConfig(TypedDict): - """A dictionary of options for generating the ell or theta. - - This dictionary contains the minimum, maximum and number of - bins to generate the ell or theta values at which to compute the statistics. - - :param minimum: The start of the binning. - :param maximum: The end of the binning. - :param n: The number of bins. - :param binning: Pass 'log' to get logarithmic spaced bins and 'lin' to get linearly - spaced bins. Default is 'log'. - - """ - - minimum: float - maximum: float - n: int - binning: str - - -def generate_ells_cells(ell_config: EllOrThetaConfig): - """Generate ells or theta values from the configuration dictionary.""" - ells = _generate_ell_or_theta(**ell_config) - Cells = np.zeros_like(ells) - - return ells, Cells - - -def generate_theta_xis(theta_config: EllOrThetaConfig): - """Generate theta and xi values from the configuration dictionary.""" - thetas = _generate_ell_or_theta(**theta_config) - xis = np.zeros_like(thetas) - - return thetas, xis - - -def apply_ells_min_max( - ells: npt.NDArray[np.int64], - Cells: npt.NDArray[np.float64], - indices: None | npt.NDArray[np.int64], - ell_min: None | int, - ell_max: None | int, -) -> tuple[ - npt.NDArray[np.int64], npt.NDArray[np.float64], None | npt.NDArray[np.int64] -]: - """Apply the minimum and maximum ell values to the ells and Cells.""" - if ell_min is not None: - locations = np.where(ells >= ell_min) - ells = ells[locations] - Cells = Cells[locations] - if indices is not None: - indices = indices[locations] - - if ell_max is not None: - locations = np.where(ells <= ell_max) - ells = ells[locations] - Cells = Cells[locations] - if indices is not None: - indices = indices[locations] - - return ells, Cells, indices - - -def apply_theta_min_max( - thetas: npt.NDArray[np.float64], - xis: npt.NDArray[np.float64], - indices: None | npt.NDArray[np.int64], - theta_min: None | float, - theta_max: None | float, -) -> tuple[ - npt.NDArray[np.float64], npt.NDArray[np.float64], None | npt.NDArray[np.int64] -]: - """Apply the minimum and maximum theta values to the thetas and xis.""" - if theta_min is not None: - locations = np.where(thetas >= theta_min) - thetas = thetas[locations] - xis = xis[locations] - if indices is not None: - indices = indices[locations] - - if theta_max is not None: - locations = np.where(thetas <= theta_max) - thetas = thetas[locations] - xis = xis[locations] - if indices is not None: - indices = indices[locations] - - return thetas, xis, indices - - -def use_source_factory( - inferred_galaxy_zdist: InferredGalaxyZDist, - wl_factory: WeakLensingFactory | None = None, - nc_factory: NumberCountsFactory | None = None, -) -> WeakLensing | NumberCounts: - """Apply the factory to the inferred galaxy redshift distribution.""" - source: WeakLensing | NumberCounts - match inferred_galaxy_zdist.measured_type: - case GalaxyMeasuredType.COUNTS: - assert nc_factory is not None - source = nc_factory.create(inferred_galaxy_zdist) - case GalaxyMeasuredType.SHEAR_E | GalaxyMeasuredType.SHEAR_T: - assert wl_factory is not None - source = wl_factory.create(inferred_galaxy_zdist) - case _: - raise ValueError( - f"Measured type {inferred_galaxy_zdist.measured_type} not supported!" - ) - return source - - -def create_sacc(metadata: list[TwoPointCells]): - """Fill the SACC file with the inferred galaxy redshift distributions.""" - sacc_data = sacc.Sacc() - - dv = [] - for cell in metadata: - for inferred_galaxy_zdist in (cell.XY.x, cell.XY.y): - if inferred_galaxy_zdist.bin_name not in sacc_data.tracers: - sacc_data.add_tracer( - "NZ", - inferred_galaxy_zdist.bin_name, - inferred_galaxy_zdist.z, - inferred_galaxy_zdist.dndz, - ) - cells = np.ones_like(cell.ells) - sacc_data.add_ell_cl( - cell.get_sacc_name(), - cell.XY.x.bin_name, - cell.XY.y.bin_name, - ell=cell.ells, - x=cells, - ) - dv.append(cells) - - delta_v = np.concatenate(dv, axis=0) - sacc_data.add_covariance(np.diag(np.ones_like(delta_v))) - - return sacc_data - - -class TwoPoint(Statistic): - """A two-point statistic. - - For example, shear correlation function, galaxy-shear correlation function, etc. +# pylint: disable=unused-import,unused-wildcard-import,wildcard-import +from firecrown.likelihood.two_point import * - Parameters - ---------- - sacc_data_type : str - The kind of two-point statistic. This must be a valid SACC data type that - maps to one of the CCL correlation function kinds or a power spectra. - Possible options are +# pylint: enable=unused-import,unused-wildcard-import,wildcard-import - - galaxy_density_cl : maps to 'cl' (a CCL angular power spectrum) - - galaxy_density_xi : maps to 'gg' (a CCL angular position corr. function) - - galaxy_shearDensity_cl_e : maps to 'cl' (a CCL angular power spectrum) - - galaxy_shearDensity_xi_t : maps to 'gl' (a CCL angular cross-correlation - between position and shear) - - galaxy_shear_cl_ee : maps to 'cl' (a CCL angular power spectrum) - - galaxy_shear_xi_minus : maps to 'l-' (a CCL angular shear corr. - function xi-) - - galaxy_shear_xi_plus : maps to 'l+' (a CCL angular shear corr. - function xi-) - - cmbGalaxy_convergenceDensity_xi : maps to 'gg' (a CCL angular position - corr. function) - - cmbGalaxy_convergenceShear_xi_t : maps to 'gl' (a CCL angular cross- - correlation between position and shear) +assert not hasattr(firecrown.likelihood.two_point, "__all__") - source0 : Source - The first sources needed to compute this statistic. - source1 : Source - The second sources needed to compute this statistic. - ell_or_theta : dict, optional - A dictionary of options for generating the ell or theta values at which - to compute the statistics. This option can be used to have firecrown - generate data without the corresponding 2pt data in the input SACC file. - The options are: - - - minimun : float - The start of the binning. - - maximun : float - The end of the binning. - - n : int - The number of bins. Note that the edges of the bins start - at `min` and end at `max`. The actual bin locations will be at the - (possibly geometric) midpoint of the bin. - - binning : str, optional - Pass 'log' to get logarithmic spaced bins and 'lin' - to get linearly spaced bins. Default is 'log'. - - ell_or_theta_min : float, optional - The minimum ell or theta value to keep. This minimum is applied after - the ell or theta values are read and/or generated. - ell_or_theta_max : float, optional - The maximum ell or theta value to keep. This maximum is applied after - the ell or theta values are read and/or generated. - ell_for_xi : dict, optional - A dictionary of options for making the ell values at which to compute - Cls for use in real-space integrations. The possible keys are: - - - minimum : int, optional - The minimum angular wavenumber to use for - real-space integrations. Default is 2. - - midpoint : int, optional - The midpoint angular wavenumber to use - for real-space integrations. The angular wavenumber samples are - linearly spaced at integers between `minimum` and `midpoint`. Default - is 50. - - maximum : int, optional - The maximum angular wavenumber to use for - real-space integrations. The angular wavenumber samples are - logarithmically spaced between `midpoint` and `maximum`. Default is - 60,000. - - n_log : int, optional - The number of logarithmically spaced angular - wavenumber samples between `mid` and `max`. Default is 200. - - Attributes - ---------- - ccl_kind : str - The CCL correlation function kind or 'cl' for power spectra corresponding - to the SACC data type. - sacc_tracers : 2-tuple of str - A tuple of the SACC tracer names for this 2pt statistic. Set after a - call to read. - ell_or_theta_ : npt.NDArray[np.float64] - The final array of ell/theta values for the statistic. Set after - `compute` is called. - measured_statistic_ : npt.NDArray[np.float64] - The measured value for the statistic. - predicted_statistic_ : npt.NDArray[np.float64] - The final prediction for the statistic. Set after `compute` is called. - - """ - - def __init__( - self, - sacc_data_type: str, - source0: Source, - source1: Source, - *, - ell_for_xi: None | dict[str, int] = None, - ell_or_theta: None | EllOrThetaConfig = None, - ell_or_theta_min: None | float | int = None, - ell_or_theta_max: None | float | int = None, - ) -> None: - super().__init__() - - assert isinstance(source0, Source) - assert isinstance(source1, Source) - - self.sacc_data_type: str = sacc_data_type - self.source0: Source = source0 - self.source1 = source1 - self.ell_for_xi_config: dict[str, int] = copy.deepcopy(ELL_FOR_XI_DEFAULTS) - if ell_for_xi is not None: - self.ell_for_xi_config.update(ell_for_xi) - # What is the difference between the following 3 instance variables? - # ell_or_theta - # _ell_or_theta - # ell_or_theta_ - self.ell_or_theta_config = ell_or_theta - self.ell_or_theta_min = ell_or_theta_min - self.ell_or_theta_max = ell_or_theta_max - self.window: None | Window = None - - self.data_vector: None | DataVector = None - self.theory_vector: None | TheoryVector = None - - self.sacc_tracers: TracerNames - self.ells: None | npt.NDArray[np.int64] = None - self.thetas: None | npt.NDArray[np.float64] = None - self.mean_ells: None | npt.NDArray[np.float64] = None - self.ells_for_xi: None | npt.NDArray[np.int64] = None - - self.cells: dict[TracerNames, npt.NDArray[np.float64]] = {} - if self.sacc_data_type in SACC_DATA_TYPE_TO_CCL_KIND: - self.ccl_kind = SACC_DATA_TYPE_TO_CCL_KIND[self.sacc_data_type] - else: - raise ValueError( - f"The SACC data type {sacc_data_type}'%s' is not " f"supported!" - ) - - @classmethod - def from_metadata_cells( - cls, - metadata: list[TwoPointCells], - wl_factory: WeakLensingFactory | None = None, - nc_factory: NumberCountsFactory | None = None, - ) -> UpdatableCollection[TwoPoint]: - """Create a TwoPoint statistic from a TwoPointCells metadata object.""" - sacc_data = create_sacc(metadata) - - two_point_list = [ - cls( - cell.get_sacc_name(), - use_source_factory( - cell.XY.x, wl_factory=wl_factory, nc_factory=nc_factory - ), - use_source_factory( - cell.XY.y, wl_factory=wl_factory, nc_factory=nc_factory - ), - ) - for cell in metadata - ] - - for two_point in two_point_list: - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - two_point.read(sacc_data) - - return UpdatableCollection(two_point_list) - - def read_ell_cells( - self, sacc_data_type: str, sacc_data: sacc.Sacc, tracers: TracerNames - ) -> ( - None - | tuple[npt.NDArray[np.int64], npt.NDArray[np.float64], npt.NDArray[np.int64]] - ): - """Read and return ell and Cell.""" - ells, Cells = sacc_data.get_ell_cl(sacc_data_type, *tracers, return_cov=False) - # As version 0.13 of sacc, the method get_ell_cl returns the - # ell values and the Cl values in arrays of the same length. - assert len(ells) == len(Cells) - common_length = len(ells) - sacc_indices = None - - if common_length == 0: - return None - sacc_indices = np.atleast_1d(sacc_data.indices(self.sacc_data_type, tracers)) - assert sacc_indices is not None # Needed for mypy - assert len(sacc_indices) == common_length - - return ells, Cells, sacc_indices - - def read_theta_xis( - self, sacc_data_type: str, sacc_data: sacc.Sacc, tracers: TracerNames - ) -> ( - None - | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.int64]] - ): - """Read and return theta and xi.""" - thetas, xis = sacc_data.get_theta_xi(sacc_data_type, *tracers, return_cov=False) - # As version 0.13 of sacc, the method get_theta_xi returns the - # theta values and the xi values in arrays of the same length. - assert len(thetas) == len(xis) - - common_length = len(thetas) - if common_length == 0: - return None - sacc_indices = np.atleast_1d(sacc_data.indices(self.sacc_data_type, tracers)) - assert sacc_indices is not None # Needed for mypy - assert len(sacc_indices) == common_length - return thetas, xis, sacc_indices - - def read(self, sacc_data: sacc.Sacc) -> None: - """Read the data for this statistic from the SACC file. - - :param sacc_data: The data in the sacc format. - """ - self.sacc_tracers = self.initialize_sources(sacc_data) - - if self.ccl_kind == "cl": - self.read_harmonic_space(sacc_data) - else: - self.read_real_space(sacc_data) - - super().read(sacc_data) - - def read_real_space(self, sacc_data: sacc.Sacc): - """Read the data for this statistic from the SACC file.""" - thetas_xis_indices = self.read_theta_xis( - self.sacc_data_type, sacc_data, self.sacc_tracers - ) - # We do not support window functions for real space statistics - if thetas_xis_indices is not None: - thetas, xis, sacc_indices = thetas_xis_indices - if self.ell_or_theta_config is not None: - # If we have data from our construction, and also have data in the - # SACC object, emit a warning and use the information read from the - # SACC object. - warnings.warn( - f"Tracers '{self.sacc_tracers}' have 2pt data and you have " - f"specified `theta` in the configuration. `theta` is being " - f"ignored!", - stacklevel=2, - ) - else: - if self.ell_or_theta_config is None: - # The SACC file has no data points, just a tracer, in this case we - # are building the statistic from scratch. In this case the user - # must have set the dictionary ell_or_theta, containing the - # minimum, maximum and number of bins to generate the ell values. - raise RuntimeError( - f"Tracers '{self.sacc_tracers}' for data type " - f"'{self.sacc_data_type}' " - "have no 2pt data in the SACC file and no input theta values " - "were given!" - ) - thetas, xis = generate_theta_xis(self.ell_or_theta_config) - sacc_indices = None - assert isinstance(self.ell_or_theta_min, (float, type(None))) - assert isinstance(self.ell_or_theta_max, (float, type(None))) - thetas, xis, sacc_indices = apply_theta_min_max( - thetas, xis, sacc_indices, self.ell_or_theta_min, self.ell_or_theta_max - ) - self.ells_for_xi = _ell_for_xi(**self.ell_for_xi_config) - self.thetas = thetas - self.sacc_indices = sacc_indices - self.data_vector = DataVector.create(xis) - - def read_harmonic_space(self, sacc_data: sacc.Sacc): - """Read the data for this statistic from the SACC file.""" - ells_cells_indices = self.read_ell_cells( - self.sacc_data_type, sacc_data, self.sacc_tracers - ) - if ells_cells_indices is not None: - ells, Cells, sacc_indices = ells_cells_indices - if self.ell_or_theta_config is not None: - # If we have data from our construction, and also have data in the - # SACC object, emit a warning and use the information read from the - # SACC object. - warnings.warn( - f"Tracers '{self.sacc_tracers}' have 2pt data and you have " - f"specified `ell` in the configuration. `ell` is being ignored!", - stacklevel=2, - ) - window = extract_window_function(sacc_data, sacc_indices) - if window is not None: - # When using a window function, we do not calculate all Cl's. - # For this reason we have a default set of ells that we use - # to compute Cl's, and we have a set of ells used for - # interpolation. - window.ells_for_interpolation = calculate_ells_for_interpolation(window) - - else: - if self.ell_or_theta_config is None: - # The SACC file has no data points, just a tracer, in this case we - # are building the statistic from scratch. In this case the user - # must have set the dictionary ell_or_theta, containing the - # minimum, maximum and number of bins to generate the ell values. - raise RuntimeError( - f"Tracers '{self.sacc_tracers}' for data type " - f"'{self.sacc_data_type}' " - "have no 2pt data in the SACC file and no input ell values " - "were given!" - ) - ells, Cells = generate_ells_cells(self.ell_or_theta_config) - sacc_indices = None - - # When generating the ells and Cells we do not have a window function - window = None - assert isinstance(self.ell_or_theta_min, (int, type(None))) - assert isinstance(self.ell_or_theta_max, (int, type(None))) - ells, Cells, sacc_indices = apply_ells_min_max( - ells, Cells, sacc_indices, self.ell_or_theta_min, self.ell_or_theta_max - ) - self.ells = ells - self.window = window - self.sacc_indices = sacc_indices - self.data_vector = DataVector.create(Cells) - - def initialize_sources(self, sacc_data: sacc.Sacc) -> TracerNames: - """Initialize this TwoPoint's sources, and return the tracer names.""" - self.source0.read(sacc_data) - if self.source0 is not self.source1: - self.source1.read(sacc_data) - assert self.source0.sacc_tracer is not None - assert self.source1.sacc_tracer is not None - tracers = (self.source0.sacc_tracer, self.source1.sacc_tracer) - return TracerNames(*tracers) - - def get_data_vector(self) -> DataVector: - """Return this statistic's data vector.""" - assert self.data_vector is not None - return self.data_vector - - def compute_theory_vector_real_space(self, tools: ModelingTools) -> TheoryVector: - """Compute a two-point statistic in real space. - - This method computes the two-point statistic in real space. It first computes - the Cl's in harmonic space and then translates them to real space using CCL. - """ - tracers0 = self.source0.get_tracers(tools) - tracers1 = self.source1.get_tracers(tools) - scale0 = self.source0.get_scale() - scale1 = self.source1.get_scale() - - assert self.ccl_kind != "cl" - assert self.thetas is not None - assert self.ells_for_xi is not None - - self.cells = {} - print(self.source0.parameter_prefix, self.source1.parameter_prefix) - cells_for_xi = self.compute_cells( - self.ells_for_xi, scale0, scale1, tools, tracers0, tracers1 - ) - - theory_vector = pyccl.correlation( - tools.get_ccl_cosmology(), - ell=self.ells_for_xi, - C_ell=cells_for_xi, - theta=self.thetas / 60, - type=self.ccl_kind, - ) - assert self.data_vector is not None - return TheoryVector.create(theory_vector) - - def compute_theory_vector_harmonic_space( - self, tools: ModelingTools - ) -> TheoryVector: - """Compute a two-point statistic in harmonic space. - - This method computes the two-point statistic in harmonic space. It computes - either the Cl's at the ells provided by the SACC file or the ells required - for the window function. - """ - tracers0 = self.source0.get_tracers(tools) - tracers1 = self.source1.get_tracers(tools) - scale0 = self.source0.get_scale() - scale1 = self.source1.get_scale() - - assert self.ccl_kind == "cl" - assert self.ells is not None - - if self.window is not None: - # If a window function is provided, we need to compute the Cl's - # for the ells used in the window function. To do this, we will - # first compute the Cl's for the ells used in the interpolation - # and then interpolate the results to the ells used in the window - # function. - assert self.window.ells_for_interpolation is not None - cells_for_interpolation = self.compute_cells( - self.window.ells_for_interpolation, - scale0, - scale1, - tools, - tracers0, - tracers1, - ) - - # TODO: There is no code in Firecrown, neither test nor example, - # that exercises a theory window function in any way. - cell_interpolator = make_log_interpolator( - self.window.ells_for_interpolation, cells_for_interpolation - ) - # Deal with ell=0 and ell=1 - cells_interpolated = np.zeros(self.window.ells.size) - cells_interpolated[2:] = cell_interpolator(self.window.ells[2:]) - - # Here we left multiply the computed Cl's by the window function - # to get the final Cl's. - theory_vector = np.einsum( - "lb, l -> b", - self.window.weights, - cells_interpolated, - ) - # We also compute the mean ell value associated with each bin. - self.mean_ells = np.einsum( - "lb, l -> b", self.window.weights, self.window.ells - ) - - assert self.data_vector is not None - return TheoryVector.create(theory_vector) - - # If we get here, we are working in harmonic space without a window function. - assert self.ells is not None - theory_vector = self.compute_cells( - self.ells, - scale0, - scale1, - tools, - tracers0, - tracers1, - ) - - assert self.data_vector is not None - return TheoryVector.create(theory_vector) - - def _compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: - """Compute a two-point statistic from sources.""" - if self.ccl_kind == "cl": - return self.compute_theory_vector_harmonic_space(tools) - - return self.compute_theory_vector_real_space(tools) - - def compute_cells( - self, - ells: npt.NDArray[np.int64], - scale0: float, - scale1: float, - tools: ModelingTools, - tracers0: Sequence[Tracer], - tracers1: Sequence[Tracer], - ) -> npt.NDArray[np.float64]: - """Compute the power spectrum for the given ells and tracers.""" - for tracer0 in tracers0: - for tracer1 in tracers1: - pk_name = f"{tracer0.field}:{tracer1.field}" - tn = TracerNames(tracer0.tracer_name, tracer1.tracer_name) - if tn in self.cells: - # Already computed this combination, skipping - continue - pk = self.calculate_pk(pk_name, tools, tracer0, tracer1) - - self.cells[tn] = ( - _cached_angular_cl( - tools.get_ccl_cosmology(), - (tracer0.ccl_tracer, tracer1.ccl_tracer), - tuple(ells.tolist()), - p_of_k_a=pk, - ) - * scale0 - * scale1 - ) - # Add up all the contributions to the cells - self.cells[TRACER_NAMES_TOTAL] = np.array(sum(self.cells.values())) - theory_vector = self.cells[TRACER_NAMES_TOTAL] - return theory_vector - - def calculate_pk( - self, pk_name: str, tools: ModelingTools, tracer0: Tracer, tracer1: Tracer - ): - """Return the power spectrum named by pk_name.""" - if tools.has_pk(pk_name): - # Use existing power spectrum - pk = tools.get_pk(pk_name) - elif tracer0.has_pt or tracer1.has_pt: - if not (tracer0.has_pt and tracer1.has_pt): - # Mixture of PT and non-PT tracers - # Create a dummy matter PT tracer for the non-PT part - matter_pt_tracer = pyccl.nl_pt.PTMatterTracer() - if not tracer0.has_pt: - tracer0.pt_tracer = matter_pt_tracer - else: - tracer1.pt_tracer = matter_pt_tracer - # Compute perturbation power spectrum - - pt_calculator = tools.get_pt_calculator() - pk = pt_calculator.get_biased_pk2d( - tracer1=tracer0.pt_tracer, - tracer2=tracer1.pt_tracer, - ) - elif tracer0.has_hm or tracer1.has_hm: - # Compute halo model power spectrum - raise NotImplementedError("Halo model power spectra not supported yet") - else: - raise ValueError(f"No power spectrum for {pk_name} can be found.") - return pk +warnings.warn( + "This module is deprecated. Use firecrown.likelihood.two_point instead.", + DeprecationWarning, + stacklevel=2, +) diff --git a/firecrown/likelihood/gauss_family/student_t.py b/firecrown/likelihood/gauss_family/student_t.py index 7d03ec6b8..b82065ccb 100644 --- a/firecrown/likelihood/gauss_family/student_t.py +++ b/firecrown/likelihood/gauss_family/student_t.py @@ -1,39 +1,20 @@ -"""The Student-t likelihood.""" +"""Deprecated module with classes related to the Student T distribution.""" -from __future__ import annotations +# flake8: noqa -import numpy as np +import warnings -from firecrown.likelihood.gauss_family.gauss_family import GaussFamily -from firecrown.modeling_tools import ModelingTools -from firecrown.likelihood.gauss_family.statistic.statistic import Statistic -from firecrown import parameters +import firecrown.likelihood.student_t +# pylint: disable=unused-import,unused-wildcard-import,wildcard-import +from firecrown.likelihood.student_t import * -class StudentT(GaussFamily): - r"""A T-distribution for the log-likelihood. +# pylint: enable=unused-import,unused-wildcard-import,wildcard-import - This distribution is appropriate when the covariance has been obtained - from a finite number of simulations. See Sellentin & Heavens - (2016; arXiv:1511.05969). As the number of simulations increases, the - T-distribution approaches a Gaussian. +assert not hasattr(firecrown.likelihood.student_t, "__all__") - :param statistics: list of statistics to build the theory and data vectors - :param nu: The Student-t $\nu$ parameter - """ - - def __init__( - self, - statistics: list[Statistic], - nu: None | float = None, - ): - super().__init__(statistics) - self.nu = parameters.register_new_updatable_parameter(nu) - - def compute_loglike(self, tools: ModelingTools): - """Compute the log-likelihood. - - :param cosmo: Current Cosmology object - """ - chi2 = self.compute_chisq(tools) - return -0.5 * self.nu * np.log(1.0 + chi2 / (self.nu - 1.0)) +warnings.warn( + "This module is deprecated. Use firecrown.likelihood.student_t instead.", + DeprecationWarning, + stacklevel=2, +) diff --git a/firecrown/likelihood/gaussfamily.py b/firecrown/likelihood/gaussfamily.py new file mode 100644 index 000000000..8880283c0 --- /dev/null +++ b/firecrown/likelihood/gaussfamily.py @@ -0,0 +1,391 @@ +"""Support for the family of Gaussian likelihood.""" + +from __future__ import annotations + +from enum import Enum +from functools import wraps +from typing import Sequence, Callable, TypeVar +from typing import final +import warnings + +from typing_extensions import ParamSpec +import numpy as np +import numpy.typing as npt +import scipy.linalg + +import sacc + +# firecrown is needed for backward compatibility; remove support for deprecated +# directory structure is removed. +import firecrown # pylint: disable=unused-import # noqa: F401 +from firecrown.parameters import ParamsMap +from firecrown.likelihood.likelihood import Likelihood +from firecrown.modeling_tools import ModelingTools +from firecrown.updatable import UpdatableCollection +from firecrown.likelihood.statistic import ( + Statistic, + GuardedStatistic, +) +from firecrown.utils import save_to_sacc + + +class State(Enum): + """The states used in GaussFamily.""" + + INITIALIZED = 1 + READY = 2 + UPDATED = 3 + COMPUTED = 4 + + +T = TypeVar("T") +P = ParamSpec("P") + + +# See https://peps.python.org/pep-0612/ and +# https://stackoverflow.com/questions/66408662/type-annotations-for-decorators +# for how to specify the types of *args and **kwargs, and the return type of +# the method being decorated. + + +# Beware +def enforce_states( + *, + initial: State | list[State], + terminal: None | State = None, + failure_message: str, +) -> Callable[[Callable[P, T]], Callable[P, T]]: + """This decorator wraps a method, and enforces state machine behavior. + + If the object is not in one of the states in initial, an + AssertionError is raised with the given failure_message. + If terminal is None the state of the object is not modified. + If terminal is not None and the call to the wrapped method returns + normally the state of the object is set to terminal. + """ + initials: list[State] + if isinstance(initial, list): + initials = initial + else: + initials = [initial] + + def decorator_enforce_states(func: Callable[P, T]) -> Callable[P, T]: + """Part of the decorator which is the closure. + + This closure is what actually contains the values of initials, terminal, and + failure_message. + """ + + @wraps(func) + def wrapper_repeat(*args: P.args, **kwargs: P.kwargs) -> T: + """Part of the decorator which is the actual wrapped method. + + It is responsible for confirming a correct initial state, and + establishing the correct final state if the wrapped method + succeeds. + """ + # The odd use of args[0] instead of self seems to be the only way + # to have both the Python runtime and mypy agree on what is being + # passed to the method, and to allow access to the attribute + # 'state'. Recall that the syntax: + # o.foo() + # calls a bound function object accessible as o.foo; this bound + # function object calls the function foo() passing 'o' as the + # first argument, self. + assert isinstance(args[0], GaussFamily) + assert args[0].state in initials, failure_message + value = func(*args, **kwargs) + if terminal is not None: + args[0].state = terminal + return value + + return wrapper_repeat + + return decorator_enforce_states + + +class GaussFamily(Likelihood): + """GaussFamily is the base class for likelihoods based on a chi-squared calculation. + + It provides an implementation of Likelihood.compute_chisq. Derived classes must + implement the abstract method compute_loglike, which is inherited from Likelihood. + + GaussFamily (and all classes that inherit from it) must abide by the the + following rules regarding the order of calling of methods. + + 1. after a new object is created, :meth:`read` must be called before any + other method in the interfaqce. + 2. after :meth:`read` has been called it is legal to call + :meth:`get_data_vector`, or to call :meth:`update`. + 3. after :meth:`update` is called it is then legal to call + :meth:`calculate_loglike` or :meth:`get_data_vector`, or to reset + the object (returning to the pre-update state) by calling + :meth:`reset`. It is also legal to call :meth:`compute_theory_vector`. + 4. after :meth:`compute_theory_vector` is called it is legal to call + :meth:`get_theory_vector` to retrieve the already-calculated theory + vector. + + This state machine behavior is enforced through the use of the decorator + :meth:`enforce_states`, above. + """ + + def __init__( + self, + statistics: Sequence[Statistic], + ): + """Initialize the base class parts of a GaussFamily object.""" + super().__init__() + self.state: State = State.INITIALIZED + if len(statistics) == 0: + raise ValueError("GaussFamily requires at least one statistic") + + for i, s in enumerate(statistics): + if not isinstance(s, Statistic): + raise ValueError( + f"statistics[{i}] is not an instance of Statistic: {s}" + f"it is a {type(s)} instead." + ) + + self.statistics: UpdatableCollection[GuardedStatistic] = UpdatableCollection( + GuardedStatistic(s) for s in statistics + ) + self.cov: None | npt.NDArray[np.float64] = None + self.cholesky: None | npt.NDArray[np.float64] = None + self.inv_cov: None | npt.NDArray[np.float64] = None + self.cov_index_map: None | dict[int, int] = None + self.theory_vector: None | npt.NDArray[np.double] = None + self.data_vector: None | npt.NDArray[np.double] = None + + @enforce_states( + initial=State.READY, + terminal=State.UPDATED, + failure_message="read() must be called before update()", + ) + def _update(self, _: ParamsMap) -> None: + """Handle the state resetting required by :class:`GaussFamily` likelihoods. + + Any derived class that needs to implement :meth:`_update` + for its own reasons must be sure to do what this does: check the state + at the start of the method, and change the state at the end of the + method. + """ + + @enforce_states( + initial=[State.UPDATED, State.COMPUTED], + terminal=State.READY, + failure_message="update() must be called before reset()", + ) + def _reset(self) -> None: + """Handle the state resetting required by :class:`GaussFamily` likelihoods. + + Any derived class that needs to implement :meth:`reset` + for its own reasons must be sure to do what this does: check the state + at the start of the method, and change the state at the end of the + method. + """ + self.theory_vector = None + + @enforce_states( + initial=State.INITIALIZED, + terminal=State.READY, + failure_message="read() must only be called once", + ) + def read(self, sacc_data: sacc.Sacc) -> None: + """Read the covariance matrix for this likelihood from the SACC file.""" + if sacc_data.covariance is None: + msg = ( + f"The {type(self).__name__} likelihood requires a covariance, " + f"but the SACC data object being read does not have one." + ) + raise RuntimeError(msg) + + covariance = sacc_data.covariance.dense + + indices_list = [] + data_vector_list = [] + for stat in self.statistics: + stat.read(sacc_data) + if stat.statistic.sacc_indices is None: + raise RuntimeError( + f"The statistic {stat.statistic} has no sacc_indices." + ) + indices_list.append(stat.statistic.sacc_indices.copy()) + data_vector_list.append(stat.statistic.get_data_vector()) + + indices = np.concatenate(indices_list) + data_vector = np.concatenate(data_vector_list) + cov = np.zeros((len(indices), len(indices))) + + for new_i, old_i in enumerate(indices): + for new_j, old_j in enumerate(indices): + cov[new_i, new_j] = covariance[old_i, old_j] + + self.data_vector = data_vector + self.cov_index_map = {old_i: new_i for new_i, old_i in enumerate(indices)} + self.cov = cov + self.cholesky = scipy.linalg.cholesky(self.cov, lower=True) + self.inv_cov = np.linalg.inv(cov) + + @final + @enforce_states( + initial=[State.READY, State.UPDATED, State.COMPUTED], + failure_message="read() must be called before get_cov()", + ) + def get_cov( + self, statistic: Statistic | list[Statistic] | None = None + ) -> npt.NDArray[np.float64]: + """Gets the current covariance matrix. + + :param statistic: The statistic for which the sub-covariance matrix + should be returned. If not specified, return the covariance of all + statistics. + """ + assert self.cov is not None + if statistic is None: + return self.cov + + assert self.cov_index_map is not None + if isinstance(statistic, Statistic): + statistic_list = [statistic] + else: + statistic_list = statistic + indices: list[int] = [] + for stat in statistic_list: + assert stat.sacc_indices is not None + temp = [self.cov_index_map[idx] for idx in stat.sacc_indices] + indices += temp + ixgrid = np.ix_(indices, indices) + + return self.cov[ixgrid] + + @final + @enforce_states( + initial=[State.READY, State.UPDATED, State.COMPUTED], + failure_message="read() must be called before get_data_vector()", + ) + def get_data_vector(self) -> npt.NDArray[np.float64]: + """Get the data vector from all statistics in the right order.""" + assert self.data_vector is not None + return self.data_vector + + @final + @enforce_states( + initial=[State.UPDATED, State.COMPUTED], + terminal=State.COMPUTED, + failure_message="update() must be called before compute_theory_vector()", + ) + def compute_theory_vector(self, tools: ModelingTools) -> npt.NDArray[np.float64]: + """Computes the theory vector using the current instance of pyccl.Cosmology. + + :param tools: Current ModelingTools object + """ + theory_vector_list: list[npt.NDArray[np.float64]] = [ + stat.compute_theory_vector(tools) for stat in self.statistics + ] + self.theory_vector = np.concatenate(theory_vector_list) + return self.theory_vector + + @final + @enforce_states( + initial=State.COMPUTED, + failure_message="compute_theory_vector() must be called before " + "get_theory_vector()", + ) + def get_theory_vector(self) -> npt.NDArray[np.float64]: + """Get the theory vector from all statistics in the right order.""" + assert ( + self.theory_vector is not None + ), "theory_vector is None after compute_theory_vector() has been called" + return self.theory_vector + + @final + @enforce_states( + initial=State.UPDATED, + failure_message="update() must be called before compute()", + ) + def compute( + self, tools: ModelingTools + ) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: + """Calculate and return both the data and theory vectors.""" + warnings.simplefilter("always", DeprecationWarning) + warnings.warn( + "The use of the `compute` method on Statistic is deprecated." + "The Statistic objects should implement `get_data` and " + "`compute_theory_vector` instead.", + category=DeprecationWarning, + ) + return self.get_data_vector(), self.compute_theory_vector(tools) + + @final + @enforce_states( + initial=[State.UPDATED, State.COMPUTED], + terminal=State.COMPUTED, + failure_message="update() must be called before compute_chisq()", + ) + def compute_chisq(self, tools: ModelingTools) -> float: + """Calculate and return the chi-squared for the given cosmology.""" + theory_vector: npt.NDArray[np.float64] + data_vector: npt.NDArray[np.float64] + residuals: npt.NDArray[np.float64] + try: + theory_vector = self.compute_theory_vector(tools) + data_vector = self.get_data_vector() + except NotImplementedError: + data_vector, theory_vector = self.compute(tools) + + assert len(data_vector) == len(theory_vector) + residuals = data_vector - theory_vector + + x = scipy.linalg.solve_triangular(self.cholesky, residuals, lower=True) + chisq = np.dot(x, x) + return chisq + + @enforce_states( + initial=[State.READY, State.UPDATED, State.COMPUTED], + failure_message="read() must be called before get_sacc_indices()", + ) + def get_sacc_indices( + self, statistic: Statistic | list[Statistic] | None = None + ) -> npt.NDArray[np.int64]: + """Get the SACC indices of the statistic or list of statistics. + + If no statistic is given, get the indices of all statistics of the likelihood. + """ + if statistic is None: + statistic = [stat.statistic for stat in self.statistics] + if isinstance(statistic, Statistic): + statistic = [statistic] + + sacc_indices_list = [] + for stat in statistic: + assert stat.sacc_indices is not None + sacc_indices_list.append(stat.sacc_indices.copy()) + + sacc_indices = np.concatenate(sacc_indices_list) + + return sacc_indices + + @enforce_states( + initial=State.COMPUTED, + failure_message="compute_theory_vector() must be called before " + "make_realization()", + ) + def make_realization( + self, sacc_data: sacc.Sacc, add_noise: bool = True, strict: bool = True + ) -> sacc.Sacc: + """Create a new realization of the model.""" + sacc_indices = self.get_sacc_indices() + + if add_noise: + new_data_vector = self.make_realization_vector() + else: + new_data_vector = self.get_theory_vector() + + new_sacc = save_to_sacc( + sacc_data=sacc_data, + data_vector=new_data_vector, + indices=sacc_indices, + strict=strict, + ) + + return new_sacc diff --git a/firecrown/likelihood/gaussian.py b/firecrown/likelihood/gaussian.py new file mode 100644 index 000000000..f503e6b6d --- /dev/null +++ b/firecrown/likelihood/gaussian.py @@ -0,0 +1,28 @@ +"""Provides GaussFamily concrete types.""" + +from __future__ import annotations +import numpy as np + +# firecrown is needed for backward compatibility; remove support for deprecated +# directory structure is removed. +import firecrown # pylint: disable=unused-import # noqa: F401 +from firecrown.likelihood.gaussfamily import GaussFamily +from firecrown.modeling_tools import ModelingTools + + +class ConstGaussian(GaussFamily): + """A Gaussian log-likelihood with a constant covariance matrix.""" + + def compute_loglike(self, tools: ModelingTools): + """Compute the log-likelihood.""" + return -0.5 * self.compute_chisq(tools) + + def make_realization_vector(self) -> np.ndarray: + """Create a new realization of the model.""" + theory_vector = self.get_theory_vector() + assert self.cholesky is not None + new_data_vector = theory_vector + np.dot( + self.cholesky, np.random.randn(len(theory_vector)) + ) + + return new_data_vector diff --git a/firecrown/likelihood/number_counts.py b/firecrown/likelihood/number_counts.py new file mode 100644 index 000000000..0294c36d1 --- /dev/null +++ b/firecrown/likelihood/number_counts.py @@ -0,0 +1,565 @@ +"""Number counts source and systematics.""" + +from __future__ import annotations + +from abc import abstractmethod +from dataclasses import dataclass, replace +from typing import Sequence, final + +import numpy as np +import numpy.typing as npt +import pyccl + +# firecrown is needed for backward compatibility; remove support for deprecated +# directory structure is removed. +import firecrown # pylint: disable=unused-import # noqa: F401 +from firecrown import parameters +from firecrown.likelihood.source import ( + SourceGalaxy, + SourceGalaxyArgs, + SourceGalaxyPhotoZShift, + SourceGalaxySelectField, + SourceGalaxySystematic, + Tracer, +) +from firecrown.metadata.two_point import InferredGalaxyZDist +from firecrown.modeling_tools import ModelingTools +from firecrown.parameters import DerivedParameter, DerivedParameterCollection, ParamsMap +from firecrown.updatable import UpdatableCollection + + +@dataclass(frozen=True) +class NumberCountsArgs(SourceGalaxyArgs): + """Class for number counts tracer builder argument.""" + + bias: None | npt.NDArray[np.float64] = None + mag_bias: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None + has_pt: bool = False + has_hm: bool = False + b_2: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None + b_s: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None + + +class NumberCountsSystematic(SourceGalaxySystematic[NumberCountsArgs]): + """Abstract base class for systematics for Number Counts sources. + + Derived classes must implement :python`apply` with the correct signature. + """ + + @abstractmethod + def apply( + self, tools: ModelingTools, tracer_arg: NumberCountsArgs + ) -> NumberCountsArgs: + """Apply method to include systematics in the tracer_arg.""" + + +class PhotoZShift(SourceGalaxyPhotoZShift[NumberCountsArgs]): + """Photo-z shift systematic.""" + + +class SelectField(SourceGalaxySelectField[NumberCountsArgs]): + """Systematic to select 3D field.""" + + +LINEAR_BIAS_DEFAULT_ALPHAZ = 0.0 +LINEAR_BIAS_DEFAULT_ALPHAG = 1.0 +LINEAR_BIAS_DEFAULT_Z_PIV = 0.5 + + +class LinearBiasSystematic(NumberCountsSystematic): + """Linear bias systematic. + + This systematic adds a linear bias model which varies with redshift and + the growth function. + + The following parameters are special Updatable parameters, which means that + they can be updated by the sampler, sacc_tracer is going to be used as a + prefix for the parameters: + + :ivar alphaz: the redshift exponent of the bias. + :ivar alphag: the growth function exponent of the bias. + :ivar z_piv: the pivot redshift of the bias. + """ + + def __init__(self, sacc_tracer: str): + """Initialize the LinearBiasSystematic. + + :param sacc_tracer: the name of the tracer in the SACC file. This is used + as a prefix for its parameters. + + """ + super().__init__(parameter_prefix=sacc_tracer) + + self.alphaz = parameters.register_new_updatable_parameter( + default_value=LINEAR_BIAS_DEFAULT_ALPHAZ + ) + self.alphag = parameters.register_new_updatable_parameter( + default_value=LINEAR_BIAS_DEFAULT_ALPHAG + ) + self.z_piv = parameters.register_new_updatable_parameter( + default_value=LINEAR_BIAS_DEFAULT_Z_PIV + ) + + def apply( + self, tools: ModelingTools, tracer_arg: NumberCountsArgs + ) -> NumberCountsArgs: + """Apply a linear bias systematic. + + Parameters + ---------- + cosmo : Cosmology + A Cosmology object. + tracer_arg : NumberCountsArgs + The source to which apply the shear bias. + """ + ccl_cosmo = tools.get_ccl_cosmology() + pref = ((1.0 + tracer_arg.z) / (1.0 + self.z_piv)) ** self.alphaz + pref *= ( + pyccl.growth_factor(ccl_cosmo, 1.0 / (1.0 + tracer_arg.z)) ** self.alphag + ) + + if tracer_arg.bias is None: + bias = np.ones_like(tracer_arg.z) + else: + bias = tracer_arg.bias + bias = bias * pref + + return replace( + tracer_arg, + bias=bias, + ) + + +PT_NON_LINEAR_BIAS_DEFAULT_B_2 = 1.0 +PT_NON_LINEAR_BIAS_DEFAULT_B_S = 1.0 + + +class PTNonLinearBiasSystematic(NumberCountsSystematic): + """Non-linear bias systematic. + + This systematic adds a linear bias model which varies with redshift and + the growth function. + + The following parameters are special Updatable parameters, which means that + they can be updated by the sampler, sacc_tracer is going to be used as a + prefix for the parameters: + + :ivar b_2: the quadratic bias. + :ivar b_s: the stochastic bias. + """ + + def __init__(self, sacc_tracer: str): + """Initialize the PTNonLinearBiasSystematic. + + :param sacc_tracer: the name of the tracer in the SACC file. This is used + as a prefix for its parameters. + + """ + super().__init__(parameter_prefix=sacc_tracer) + + self.b_2 = parameters.register_new_updatable_parameter( + default_value=PT_NON_LINEAR_BIAS_DEFAULT_B_2 + ) + self.b_s = parameters.register_new_updatable_parameter( + default_value=PT_NON_LINEAR_BIAS_DEFAULT_B_S + ) + + def apply( + self, tools: ModelingTools, tracer_arg: NumberCountsArgs + ) -> NumberCountsArgs: + """Apply a non-linear bias systematic. + + :param tools: currently unused, but required by interface + :param tracer_arg: a NumberCountsArgs object with values to be updated + + :return: the updated NumberCountsArgs object + """ + z = tracer_arg.z + b_2_z = self.b_2 * np.ones_like(z) + b_s_z = self.b_s * np.ones_like(z) + # b_1 uses the "bias" field + return replace( + tracer_arg, + has_pt=True, + b_2=(z, b_2_z), + b_s=(z, b_s_z), + ) + + +class MagnificationBiasSystematic(NumberCountsSystematic): + """Magnification bias systematic. + + This systematic adds a magnification bias model for galaxy number contrast + following Joachimi & Bridle (2010), arXiv:0911.2454. + + The following parameters are special Updatable parameters, which means that + they can be updated by the sampler, sacc_tracer is going to be used as a + prefix for the parameters: + + :ivar r_lim: the limiting magnitude. + :ivar sig_c: the intrinsic dispersion of the source redshift distribution. + :ivar eta: the slope of the luminosity function. + :ivar z_c: the characteristic redshift of the source distribution. + :ivar z_m: the slope of the source redshift distribution. + """ + + def __init__(self, sacc_tracer: str): + """Initialize the MagnificationBiasSystematic. + + :param sacc_tracer: the name of the tracer in the SACC file. This is used + as a prefix for its parameters. + """ + super().__init__(parameter_prefix=sacc_tracer) + + self.r_lim = parameters.register_new_updatable_parameter() + self.sig_c = parameters.register_new_updatable_parameter() + self.eta = parameters.register_new_updatable_parameter() + self.z_c = parameters.register_new_updatable_parameter() + self.z_m = parameters.register_new_updatable_parameter() + + def apply( + self, tools: ModelingTools, tracer_arg: NumberCountsArgs + ) -> NumberCountsArgs: + """Apply a magnification bias systematic. + + :param tools: a ModelingTools object + :param tracer_arg: a NumberCountsArgs object + + :return: a NumberCountsArgs object + """ + z_bar = self.z_c + self.z_m * (self.r_lim - 24.0) + # The slope of log(n_tot(z,r_lim)) with respect to r_lim + # where n_tot(z,r_lim) is the luminosity function after using fit (C.1) + s = ( + self.eta / self.r_lim + - 3.0 * self.z_m / z_bar + + 1.5 * self.z_m * np.power(tracer_arg.z / z_bar, 1.5) / z_bar + ) + + if tracer_arg.mag_bias is None: + mag_bias = np.ones_like(tracer_arg.z) + else: + mag_bias = tracer_arg.mag_bias[1] + mag_bias = mag_bias * s / np.log(10) + + return replace( + tracer_arg, + mag_bias=(tracer_arg.z, mag_bias), + ) + + +CONSTANT_MAGNIFICATION_BIAS_DEFAULT_MAG_BIAS = 1.0 + + +class ConstantMagnificationBiasSystematic(NumberCountsSystematic): + """Simple constant magnification bias systematic. + + This systematic adds a constant magnification bias model for galaxy number + contrast. + + The following parameters are special Updatable parameters, which means that + they can be updated by the sampler, sacc_tracer is going to be used as a + prefix for the parameters: + + :ivar mag_bias: the magnification bias. + """ + + def __init__(self, sacc_tracer: str): + """Initialize the ConstantMagnificationBiasSystematic. + + :param sacc_tracer: the name of the tracer in the SACC file. This is used + as a prefix for its parameters. + """ + super().__init__(parameter_prefix=sacc_tracer) + + self.mag_bias = parameters.register_new_updatable_parameter( + default_value=CONSTANT_MAGNIFICATION_BIAS_DEFAULT_MAG_BIAS + ) + + def apply( + self, tools: ModelingTools, tracer_arg: NumberCountsArgs + ) -> NumberCountsArgs: + """Apply a constant magnification bias systematic. + + :param tools: currently unused, but required by interface + :param tracer_arg: a NumberCountsArgs object with values to be updated + + :return: the updated NumberCountsArgs object + """ + return replace( + tracer_arg, + mag_bias=(tracer_arg.z, np.ones_like(tracer_arg.z) * self.mag_bias), + ) + + +NUMBER_COUNTS_DEFAULT_BIAS = 1.5 + + +class NumberCounts(SourceGalaxy[NumberCountsArgs]): + """Source class for number counts.""" + + def __init__( + self, + *, + sacc_tracer: str, + has_rsd: bool = False, + derived_scale: bool = False, + scale: float = 1.0, + systematics: None | Sequence[SourceGalaxySystematic[NumberCountsArgs]] = None, + ): + """Initialize the NumberCounts object. + + :param sacc_tracer: the name of the tracer in the SACC file. This is used + as a prefix for its parameters. + :param has_rsd: whether to include RSD in the tracer. + :param derived_scale: whether to include a derived parameter for the scale + of the tracer. + :param scale: the initial scale of the tracer. + :param systematics: a list of systematics to apply to the tracer. + """ + super().__init__(sacc_tracer=sacc_tracer, systematics=systematics) + + self.sacc_tracer = sacc_tracer + self.has_rsd = has_rsd + self.derived_scale = derived_scale + + self.bias = parameters.register_new_updatable_parameter( + default_value=NUMBER_COUNTS_DEFAULT_BIAS + ) + self.systematics: UpdatableCollection[ + SourceGalaxySystematic[NumberCountsArgs] + ] = UpdatableCollection(systematics) + self.scale = scale + self.current_tracer_args: None | NumberCountsArgs = None + self.tracer_args: NumberCountsArgs + + @classmethod + def create_ready( + cls, + inferred_zdist: InferredGalaxyZDist, + has_rsd: bool = False, + derived_scale: bool = False, + scale: float = 1.0, + systematics: None | list[SourceGalaxySystematic[NumberCountsArgs]] = None, + ) -> NumberCounts: + """Create a NumberCounts object with the given tracer name and scale.""" + obj = cls( + sacc_tracer=inferred_zdist.bin_name, + systematics=systematics, + has_rsd=has_rsd, + derived_scale=derived_scale, + scale=scale, + ) + # pylint: disable=unexpected-keyword-arg + obj.tracer_args = NumberCountsArgs( + scale=obj.scale, + z=inferred_zdist.z, + dndz=inferred_zdist.dndz, + bias=None, + mag_bias=None, + ) + # pylint: enable=unexpected-keyword-arg + + return obj + + @final + def _update_source(self, params: ParamsMap): + """Perform any updates necessary after the parameters have being updated. + + This implementation must update all contained Updatable instances. + """ + self.systematics.update(params) + + @final + def _get_derived_parameters(self) -> DerivedParameterCollection: + if self.derived_scale: + assert self.current_tracer_args is not None + derived_scale = DerivedParameter( + "TwoPoint", + f"NumberCountsScale_{self.sacc_tracer}", + self.current_tracer_args.scale, + ) + derived_parameters = DerivedParameterCollection([derived_scale]) + else: + derived_parameters = DerivedParameterCollection([]) + + return derived_parameters + + def _read(self, sacc_data): + """Read the data for this source from the SACC file. + + Parameters + ---------- + sacc_data : sacc.Sacc + The data in the sacc format. + """ + # pylint: disable=unexpected-keyword-arg + self.tracer_args = NumberCountsArgs( + scale=self.scale, + z=np.array([]), + dndz=np.array([]), + bias=None, + mag_bias=None, + ) + # pylint: enable=unexpected-keyword-arg + super()._read(sacc_data) + + def create_tracers( + self, tools: ModelingTools + ) -> tuple[list[Tracer], NumberCountsArgs]: + """Create the tracers for this source.""" + tracer_args = self.tracer_args + tracer_args = replace(tracer_args, bias=self.bias * np.ones_like(tracer_args.z)) + + ccl_cosmo = tools.get_ccl_cosmology() + for systematic in self.systematics: + tracer_args = systematic.apply(tools, tracer_args) + + tracers = [] + + if not tracer_args.has_pt or tracer_args.mag_bias is not None or self.has_rsd: + # Create a normal pyccl.NumberCounts tracer if there's no PT, or + # in case there's magnification or RSD. + tracer_names = [] + if tracer_args.has_pt: + # use PT for galaxy bias + bias = None + else: + bias = (tracer_args.z, tracer_args.bias) + tracer_names += ["galaxies"] + if tracer_args.mag_bias is not None: + tracer_names += ["magnification"] + if self.has_rsd: + tracer_names += ["rsd"] + + ccl_mag_tracer = pyccl.NumberCountsTracer( + ccl_cosmo, + has_rsd=self.has_rsd, + dndz=(tracer_args.z, tracer_args.dndz), + bias=bias, + mag_bias=tracer_args.mag_bias, + ) + + tracers.append( + Tracer( + ccl_mag_tracer, + tracer_name="+".join(tracer_names), + field=tracer_args.field, + ) + ) + if tracer_args.has_pt: + nc_pt_tracer = pyccl.nl_pt.PTNumberCountsTracer( + b1=(tracer_args.z, tracer_args.bias), + b2=tracer_args.b_2, + bs=tracer_args.b_s, + ) + + ccl_nc_dummy_tracer = pyccl.NumberCountsTracer( + ccl_cosmo, + has_rsd=False, + dndz=(tracer_args.z, tracer_args.dndz), + bias=(tracer_args.z, np.ones_like(tracer_args.z)), + ) + nc_pt_tracer = Tracer( + ccl_nc_dummy_tracer, tracer_name="galaxies", pt_tracer=nc_pt_tracer + ) + tracers.append(nc_pt_tracer) + + self.current_tracer_args = tracer_args + + return tracers, tracer_args + + def get_scale(self): + """Return the scale for this source.""" + assert self.current_tracer_args + return self.current_tracer_args.scale + + +class NumberCountsSystematicFactory: + """Factory class for NumberCountsSystematic objects.""" + + @abstractmethod + def create( + self, inferred_zdist: InferredGalaxyZDist + ) -> SourceGalaxySystematic[NumberCountsArgs]: + """Create a NumberCountsSystematic object with the given tracer name.""" + + +class PhotoZShiftFactory(NumberCountsSystematicFactory): + """Factory class for PhotoZShift objects.""" + + def create(self, inferred_zdist: InferredGalaxyZDist) -> PhotoZShift: + """Create a PhotoZShift object with the given tracer name.""" + return PhotoZShift(inferred_zdist.bin_name) + + +class LinearBiasSystematicFactory(NumberCountsSystematicFactory): + """Factory class for LinearBiasSystematic objects.""" + + def create(self, inferred_zdist: InferredGalaxyZDist) -> LinearBiasSystematic: + """Create a LinearBiasSystematic object with the given tracer name.""" + return LinearBiasSystematic(inferred_zdist.bin_name) + + +class PTNonLinearBiasSystematicFactory(NumberCountsSystematicFactory): + """Factory class for PTNonLinearBiasSystematic objects.""" + + def create(self, inferred_zdist: InferredGalaxyZDist) -> PTNonLinearBiasSystematic: + """Create a PTNonLinearBiasSystematic object with the given tracer name.""" + return PTNonLinearBiasSystematic(inferred_zdist.bin_name) + + +class MagnificationBiasSystematicFactory(NumberCountsSystematicFactory): + """Factory class for MagnificationBiasSystematic objects.""" + + def create( + self, inferred_zdist: InferredGalaxyZDist + ) -> MagnificationBiasSystematic: + """Create a MagnificationBiasSystematic object with the given tracer name.""" + return MagnificationBiasSystematic(inferred_zdist.bin_name) + + +class ConstantMagnificationBiasSystematicFactory(NumberCountsSystematicFactory): + """Factory class for ConstantMagnificationBiasSystematic objects.""" + + def create( + self, inferred_zdist: InferredGalaxyZDist + ) -> ConstantMagnificationBiasSystematic: + """Create a ConstantMagnificationBiasSystematic object. + + Use the inferred_zdist to create the systematic. + """ + return ConstantMagnificationBiasSystematic(inferred_zdist.bin_name) + + +class NumberCountsFactory: + """Factory class for NumberCounts objects.""" + + def __init__( + self, + per_bin_systematics: list[NumberCountsSystematicFactory], + global_systematics: Sequence[NumberCountsSystematic], + ) -> None: + """Initialize the NumberCountsFactory.""" + self.per_bin_systematics: list[NumberCountsSystematicFactory] = ( + per_bin_systematics + ) + self.global_systematics: Sequence[NumberCountsSystematic] = global_systematics + self.cache: dict[int, NumberCounts] = {} + + def create(self, inferred_zdist: InferredGalaxyZDist) -> NumberCounts: + """Create a NumberCounts object with the given tracer name and scale.""" + inferred_zdist_id = id(inferred_zdist) + if inferred_zdist_id in self.cache: + return self.cache[inferred_zdist_id] + + systematics: list[SourceGalaxySystematic[NumberCountsArgs]] = [ + systematic_factory.create(inferred_zdist) + for systematic_factory in self.per_bin_systematics + ] + systematics.extend(self.global_systematics) + + nc = NumberCounts.create_ready(inferred_zdist, systematics=systematics) + self.cache[inferred_zdist_id] = nc + + return nc diff --git a/firecrown/likelihood/source.py b/firecrown/likelihood/source.py new file mode 100644 index 000000000..8b6c491dd --- /dev/null +++ b/firecrown/likelihood/source.py @@ -0,0 +1,327 @@ +"""Abstract base classes for TwoPoint Statistics sources.""" + +from __future__ import annotations + +from abc import abstractmethod +from dataclasses import dataclass, replace +from typing import Generic, Sequence, TypeVar, final + +import numpy as np +import numpy.typing as npt +import pyccl +import pyccl.nl_pt +import sacc +from scipy.interpolate import Akima1DInterpolator + +# firecrown is needed for backward compatibility; remove support for deprecated +# directory structure is removed. +import firecrown # pylint: disable=unused-import # noqa: F401 +from firecrown import parameters +from firecrown.modeling_tools import ModelingTools +from firecrown.parameters import ParamsMap +from firecrown.updatable import Updatable, UpdatableCollection + + +class SourceSystematic(Updatable): + """An abstract systematic class (e.g., shear biases, photo-z shifts, etc.). + + This class currently has no methods at all, because the argument types for + the `apply` method of different subclasses are different. + """ + + def read(self, sacc_data: sacc.Sacc): + """Call to allow this object to read from the appropriate sacc data.""" + + +class Source(Updatable): + """An abstract source class (e.g., a sample of lenses). + + Parameters + ---------- + systematics : list of str, optional + A list of the source-level systematics to apply to the source. The + default of `None` implies no systematics. + """ + + systematics: Sequence[SourceSystematic] + cosmo_hash: None | int + tracers: Sequence[Tracer] + + def __init__(self, sacc_tracer: str) -> None: + """Create a Source object that uses the named tracer. + + :param sacc_tracer: the name of the tracer in the SACC file. This is used + as a prefix for its parameters. + """ + super().__init__(parameter_prefix=sacc_tracer) + self.sacc_tracer = sacc_tracer + + @final + def read(self, sacc_data: sacc.Sacc): + """Read the data for this source from the SACC file.""" + if hasattr(self, "systematics"): + for systematic in self.systematics: + systematic.read(sacc_data) + self._read(sacc_data) + + @abstractmethod + def _read(self, sacc_data: sacc.Sacc): + """Abstract method to read the data for this source from the SACC file.""" + + def _update_source(self, params: ParamsMap): + """Method to update the source from the given ParamsMap. + + Any subclass that needs to do more than update its contained :class:`Updatable` + objects should implement this method. + """ + + @final + def _update(self, params: ParamsMap): + """Implementation of Updatable interface method `_update`. + + This clears the current hash and tracer, and calls the abstract method + `_update_source`, which must be implemented in all subclasses. + """ + self.cosmo_hash = None + self.tracers = [] + self._update_source(params) + + @abstractmethod + def get_scale(self) -> float: + """Abstract method to return the scales for this `Source`.""" + + @abstractmethod + def create_tracers(self, tools: ModelingTools): + """Create tracers for this `Source`, for the given cosmology.""" + + @final + def get_tracers(self, tools: ModelingTools) -> Sequence[Tracer]: + """Return the tracer for the given cosmology. + + This method caches its result, so if called a second time with the same + cosmology, no calculation needs to be done. + """ + ccl_cosmo = tools.get_ccl_cosmology() + + cur_hash = hash(ccl_cosmo) + if hasattr(self, "cosmo_hash") and self.cosmo_hash == cur_hash: + return self.tracers + + self.tracers, _ = self.create_tracers(tools) + self.cosmo_hash = cur_hash + return self.tracers + + +class Tracer: + """Extending the pyccl.Tracer object with additional information. + + Bundles together a pyccl.Tracer object with optional information about the + underlying 3D field, a pyccl.nl_pt.PTTracer, and halo profiles. + """ + + @staticmethod + def determine_field_name(field: None | str, tracer: None | str) -> str: + """Gets a field name for a tracer. + + This function encapsulates the policy for determining the value to be + assigned to the :attr:`field` attribute of a :class:`Tracer`. + + It is a static method only to keep it grouped with the class for which it is + defining the initialization policy. + """ + if field is not None: + return field + if tracer is not None: + return tracer + return "delta_matter" + + def __init__( + self, + tracer: pyccl.Tracer, + tracer_name: None | str = None, + field: None | str = None, + pt_tracer: None | pyccl.nl_pt.PTTracer = None, + halo_profile: None | pyccl.halos.HaloProfile = None, + halo_2pt: None | pyccl.halos.Profile2pt = None, + ): + """Initialize a new Tracer based on the pyccl.Tracer which must not be None. + + Note that the :class:`pyccl.Tracer` is not copied; we store a reference to the + original tracer. Be careful not to accidentally share :class:`pyccl.Tracer`s. + + If no tracer_name is supplied, then the tracer_name is set to the name of the + :class:`pyccl.Tracer` class that was used. + + If no `field` is given, then the attribute :attr:`field` is set to either + (1) the tracer_name, if one was given, or (2) 'delta_matter'. + """ + assert tracer is not None + self.ccl_tracer = tracer + self.tracer_name: str = tracer_name or tracer.__class__.__name__ + self.field = Tracer.determine_field_name(field, tracer_name) + self.pt_tracer = pt_tracer + self.halo_profile = halo_profile + self.halo_2pt = halo_2pt + + @property + def has_pt(self) -> bool: + """Return True if we have a pt_tracer, and False if not.""" + return self.pt_tracer is not None + + @property + def has_hm(self) -> bool: + """Return True if we have a halo_profile, and False if not.""" + return self.halo_profile is not None + + +# Sources of galaxy distributions + + +@dataclass(frozen=True) +class SourceGalaxyArgs: + """Class for galaxy based sources arguments.""" + + z: npt.NDArray[np.float64] + dndz: npt.NDArray[np.float64] + + scale: float = 1.0 + + field: str = "delta_matter" + + +_SourceGalaxyArgsT = TypeVar("_SourceGalaxyArgsT", bound=SourceGalaxyArgs) + + +class SourceGalaxySystematic(SourceSystematic, Generic[_SourceGalaxyArgsT]): + """Abstract base class for all galaxy based source systematics.""" + + @abstractmethod + def apply( + self, tools: ModelingTools, tracer_arg: _SourceGalaxyArgsT + ) -> _SourceGalaxyArgsT: + """Apply method to include systematics in the tracer_arg.""" + + +_SourceGalaxySystematicT = TypeVar( + "_SourceGalaxySystematicT", bound=SourceGalaxySystematic +) + + +SOURCE_GALAXY_SYSTEMATIC_DEFAULT_DELTA_Z = 0.0 + + +class SourceGalaxyPhotoZShift( + SourceGalaxySystematic[_SourceGalaxyArgsT], Generic[_SourceGalaxyArgsT] +): + """A photo-z shift bias. + + This systematic shifts the photo-z distribution by some amount `delta_z`. + + The following parameters are special Updatable parameters, which means that + they can be updated by the sampler, sacc_tracer is going to be used as a + prefix for the parameters: + + :ivar delta_z: the photo-z shift. + """ + + def __init__(self, sacc_tracer: str) -> None: + """Create a PhotoZShift object, using the specified tracer name. + + :param sacc_tracer: the name of the tracer in the SACC file. This is used + as a prefix for its parameters. + """ + super().__init__(parameter_prefix=sacc_tracer) + + self.delta_z = parameters.register_new_updatable_parameter( + default_value=SOURCE_GALAXY_SYSTEMATIC_DEFAULT_DELTA_Z + ) + + def apply(self, tools: ModelingTools, tracer_arg: _SourceGalaxyArgsT): + """Apply a shift to the photo-z distribution of a source.""" + dndz_interp = Akima1DInterpolator(tracer_arg.z, tracer_arg.dndz) + + dndz = dndz_interp(tracer_arg.z - self.delta_z, extrapolate=False) + dndz[np.isnan(dndz)] = 0.0 + + return replace( + tracer_arg, + dndz=dndz, + ) + + +class SourceGalaxySelectField( + SourceGalaxySystematic[_SourceGalaxyArgsT], Generic[_SourceGalaxyArgsT] +): + """The source galaxy select field systematic. + + A systematic that allows specifying the 3D field that will be used + to select the 3D power spectrum when computing the angular power + spectrum. + """ + + def __init__(self, field: str = "delta_matter"): + """Specify which 3D field should be used when computing angular power spectra. + + :param field: the name of the 3D field that is associated to the tracer. + Default: `"delta_matter"` + """ + super().__init__() + self.field = field + + def apply( + self, tools: ModelingTools, tracer_arg: _SourceGalaxyArgsT + ) -> _SourceGalaxyArgsT: + """Apply method to include systematics in the tracer_arg.""" + return replace(tracer_arg, field=self.field) + + +class SourceGalaxy(Source, Generic[_SourceGalaxyArgsT]): + """Source class for galaxy based sources.""" + + def __init__( + self, + *, + sacc_tracer: str, + systematics: None | Sequence[SourceGalaxySystematic] = None, + ): + """Initialize the SourceGalaxy object. + + :param sacc_tracer: the name of the tracer in the SACC file. This is used + as a prefix for its parameters. + + """ + super().__init__(sacc_tracer) + + self.sacc_tracer = sacc_tracer + self.current_tracer_args: None | _SourceGalaxyArgsT = None + self.systematics: UpdatableCollection[SourceGalaxySystematic] = ( + UpdatableCollection(systematics) + ) + self.tracer_args: _SourceGalaxyArgsT + + def _read(self, sacc_data: sacc.Sacc): + """Read the galaxy redshift distribution model from a sacc file. + + All derived classes must call this method in their own `_read` method + after they have read their own data and initialized their tracer_args. + """ + try: + tracer_args = self.tracer_args + except AttributeError as exc: + raise RuntimeError( + "Must initialize tracer_args before calling _read on SourceGalaxy" + ) from exc + + tracer = sacc_data.get_tracer(self.sacc_tracer) + + z = getattr(tracer, "z").copy().flatten() + nz = getattr(tracer, "nz").copy().flatten() + indices = np.argsort(z) + z = z[indices] + nz = nz[indices] + + self.tracer_args = replace( + tracer_args, + z=z, + dndz=nz, + ) diff --git a/firecrown/likelihood/statistic.py b/firecrown/likelihood/statistic.py new file mode 100644 index 000000000..e951ff552 --- /dev/null +++ b/firecrown/likelihood/statistic.py @@ -0,0 +1,292 @@ +"""Gaussian Family Statistic Module. + +The Statistic class describing objects that implement methods to compute the +data and theory vectors for a :class:`GaussFamily` subclass. +""" + +from __future__ import annotations +from typing import final +from dataclasses import dataclass +from abc import abstractmethod +import warnings +import numpy as np +import numpy.typing as npt +import sacc + +import firecrown.parameters +from firecrown.parameters import DerivedParameterCollection, RequiredParameters +from firecrown.modeling_tools import ModelingTools +from firecrown.updatable import Updatable + + +class DataVector(npt.NDArray[np.float64]): + """This class wraps a np.ndarray that represents some observed data values.""" + + @classmethod + def create(cls, vals: npt.NDArray[np.float64]) -> DataVector: + """Create a DataVector that wraps a copy of the given array vals.""" + return vals.view(cls) + + @classmethod + def from_list(cls, vals: list[float]) -> DataVector: + """Create a DataVector from the given list of floats.""" + array = np.array(vals) + return cls.create(array) + + +class TheoryVector(npt.NDArray[np.float64]): + """This class represents an observation predicted by some theory.""" + + @classmethod + def create(cls, vals: npt.NDArray[np.float64]) -> TheoryVector: + """Create a TheoryVector that wraps a copy of the given array vals.""" + return vals.view(cls) + + @classmethod + def from_list(cls, vals: list[float]) -> TheoryVector: + """Create a TheoryVector from the given list of floats.""" + array = np.array(vals) + return cls.create(array) + + +def residuals(data: DataVector, theory: TheoryVector) -> npt.NDArray[np.float64]: + """Return a bare np.ndarray with the difference between `data` and `theory`. + + This is to be preferred to using arithmetic on the vectors directly. + """ + assert isinstance(data, DataVector) + assert isinstance(theory, TheoryVector) + return (data - theory).view(np.ndarray) + + +@dataclass +class StatisticsResult: + """This is the type returned by the `compute` method of any `Statistic`.""" + + data: DataVector + theory: TheoryVector + + def __post_init__(self): + """Make sure the data and theory vectors are of the same shape.""" + assert self.data.shape == self.theory.shape + + def residuals(self) -> npt.NDArray[np.float64]: + """Return the residuals -- the difference between data and theory.""" + return self.data - self.theory + + def __iter__(self): + """Iterate through the data members. + + This is to allow automatic unpacking, as if the StatisticsResult were a tuple + of (data, theory). + + This method is a temporary measure to help code migrate to the newer, + safer interface for Statistic.compute(). + """ + warnings.warn( + "Iteration and tuple unpacking for StatisticsResult is " + "deprecated.\nPlease use the StatisticsResult class accessors" + ".data and .theory by name." + ) + yield self.data + yield self.theory + + +class StatisticUnreadError(RuntimeError): + """Error raised when accessing an un-read statistic. + + Run-time error indicating an attempt has been made to use a statistic + that has not had `read` called in it. + """ + + def __init__(self, stat: Statistic): + msg = ( + f"The statistic {stat} was used for calculation before `read` " + f"was called.\nIt may be that a likelihood factory function did not" + f"call `read` before returning the likelihood." + ) + super().__init__(msg) + self.statstic = stat + + +class Statistic(Updatable): + """The abstract base class for all physics-related statistics. + + Statistics read data from a SACC object as part of a multi-phase + initialization. They manage a :class:`DataVector` and, given a + :class:`ModelingTools` object, can compute a :class:`TheoryVector`. + + Statistics represent things like two-point functions and mass functions. + """ + + def __init__(self, parameter_prefix: None | str = None): + super().__init__(parameter_prefix=parameter_prefix) + self.sacc_indices: None | npt.NDArray[np.int64] + self.ready = False + self.computed_theory_vector = False + self.theory_vector: None | TheoryVector = None + + def read(self, _: sacc.Sacc) -> None: + """Read the data for this statistic and mark it as ready for use. + + Derived classes that override this function should make sure to call the base + class method using: + + .. code-block:: python + + super().read(sacc_data) + + as the last thing they do in `__init__`. + + Note that currently the argument is not used; it is present so that this + method will have the correct argument type for the override. + """ + self.ready = True + if len(self.get_data_vector()) == 0: + raise RuntimeError( + f"the statistic {self} has read a data vector " + f"of length 0; the length must be positive" + ) + + def _reset(self): + """Reset this statistic. + + All subclasses implementations must call super()._reset() + """ + self.computed_theory_vector = False + self.theory_vector = None + + @abstractmethod + def get_data_vector(self) -> DataVector: + """Gets the statistic data vector.""" + + @final + def compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: + """Compute a statistic from sources, applying any systematics.""" + if not self.is_updated(): + raise RuntimeError( + f"The statistic {self} has not been updated with parameters." + ) + self.theory_vector = self._compute_theory_vector(tools) + self.computed_theory_vector = True + + return self.theory_vector + + @abstractmethod + def _compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: + """Compute a statistic from sources, concrete implementation.""" + + def get_theory_vector(self) -> TheoryVector: + """Returns the last computed theory vector. + + Raises a RuntimeError if the vector has not been computed. + """ + if not self.computed_theory_vector: + raise RuntimeError( + f"The theory for statistic {self} has not been computed yet." + ) + assert self.theory_vector is not None, ( + "implementation error, " + "computed_theory_vector is True but theory_vector is None" + ) + return self.theory_vector + + +class GuardedStatistic(Updatable): + """An internal class used to maintain state on statistics. + + :class:`GuardedStatistic` is used by the framework to maintain and + validate the state of instances of classes derived from :class:`Statistic`. + """ + + def __init__(self, stat: Statistic): + """Initialize the GuardedStatistic to contain the given :class:`Statistic`.""" + super().__init__() + assert isinstance(stat, Statistic) + self.statistic = stat + + def read(self, sacc_data: sacc.Sacc) -> None: + """Read whatever data is needed from the given :class:`sacc.Sacc` object. + + After this function is called, the object should be prepared for the + calling of the methods :meth:`get_data_vector` and + :meth:`compute_theory_vector`. + """ + if self.statistic.ready: + raise RuntimeError("Firecrown has called read twice on a GuardedStatistic") + try: + self.statistic.read(sacc_data) + except TypeError as exc: + msg = ( + f"A statistic of type {type(self.statistic).__name__} has raised " + f"an exception during `read`.\nThe problem may be a malformed " + f"SACC data object." + ) + raise RuntimeError(msg) from exc + + def get_data_vector(self) -> DataVector: + """Return the contained :class:`Statistic`'s data vector. + + :class:`GuardedStatistic` ensures that :meth:`read` has been called. + first. + """ + if not self.statistic.ready: + raise StatisticUnreadError(self.statistic) + return self.statistic.get_data_vector() + + def compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: + """Return the contained :class:`Statistic`'s computed theory vector. + + :class:`GuardedStatistic` ensures that :meth:`read` has been called. + first. + """ + if not self.statistic.ready: + raise StatisticUnreadError(self.statistic) + return self.statistic.compute_theory_vector(tools) + + +class TrivialStatistic(Statistic): + """A minimal statistic only to be used for testing Gaussian likelihoods. + + It returns a :class:`DataVector` and :class:`TheoryVector` each of which is + three elements long. The SACC data provided to :meth:`TrivialStatistic.read` + must supply the necessary values. + """ + + def __init__(self) -> None: + """Initialize this statistic.""" + super().__init__() + # Data and theory will both be of length self.count + self.count = 3 + self.data_vector: None | DataVector = None + self.mean = firecrown.parameters.register_new_updatable_parameter( + default_value=0.0 + ) + self.computed_theory_vector = False + + def read(self, sacc_data: sacc.Sacc): + """Read the necessary items from the sacc data.""" + our_data = sacc_data.get_mean(data_type="count") + assert len(our_data) == self.count + self.data_vector = DataVector.from_list(our_data) + self.sacc_indices = np.arange(len(self.data_vector)) + super().read(sacc_data) + + @final + def _required_parameters(self) -> RequiredParameters: + """Return an empty RequiredParameters.""" + return RequiredParameters([]) + + @final + def _get_derived_parameters(self) -> DerivedParameterCollection: + """Return an empty DerivedParameterCollection.""" + return DerivedParameterCollection([]) + + def get_data_vector(self) -> DataVector: + """Return the data vector; raise exception if there is none.""" + assert self.data_vector is not None + return self.data_vector + + def _compute_theory_vector(self, _: ModelingTools) -> TheoryVector: + """Return a fixed theory vector.""" + return TheoryVector.from_list([self.mean] * self.count) diff --git a/firecrown/likelihood/student_t.py b/firecrown/likelihood/student_t.py new file mode 100644 index 000000000..aea79d450 --- /dev/null +++ b/firecrown/likelihood/student_t.py @@ -0,0 +1,42 @@ +"""The Student-t likelihood.""" + +from __future__ import annotations + +import numpy as np + +# firecrown is needed for backward compatibility; remove support for deprecated +# directory structure is removed. +import firecrown # pylint: disable=unused-import # noqa: F401 +from firecrown.likelihood.gaussfamily import GaussFamily +from firecrown.modeling_tools import ModelingTools +from firecrown.likelihood.statistic import Statistic +from firecrown import parameters + + +class StudentT(GaussFamily): + r"""A T-distribution for the log-likelihood. + + This distribution is appropriate when the covariance has been obtained + from a finite number of simulations. See Sellentin & Heavens + (2016; arXiv:1511.05969). As the number of simulations increases, the + T-distribution approaches a Gaussian. + + :param statistics: list of statistics to build the theory and data vectors + :param nu: The Student-t $\nu$ parameter + """ + + def __init__( + self, + statistics: list[Statistic], + nu: None | float = None, + ): + super().__init__(statistics) + self.nu = parameters.register_new_updatable_parameter(nu) + + def compute_loglike(self, tools: ModelingTools): + """Compute the log-likelihood. + + :param cosmo: Current Cosmology object + """ + chi2 = self.compute_chisq(tools) + return -0.5 * self.nu * np.log(1.0 + chi2 / (self.nu - 1.0)) diff --git a/firecrown/likelihood/supernova.py b/firecrown/likelihood/supernova.py new file mode 100644 index 000000000..40ee63e34 --- /dev/null +++ b/firecrown/likelihood/supernova.py @@ -0,0 +1,79 @@ +"""Supernova statistic support.""" + +from __future__ import annotations + +import numpy as np +import numpy.typing as npt +import pyccl +import sacc +from sacc.tracers import MiscTracer + +# firecrown is needed for backward compatibility; remove support for deprecated +# directory structure is removed. +import firecrown # pylint: disable=unused-import # noqa: F401 +from firecrown import parameters +from firecrown.likelihood.statistic import ( + DataVector, + Statistic, + TheoryVector, +) +from firecrown.modeling_tools import ModelingTools + +SNIA_DEFAULT_M = -19.2 + + +class Supernova(Statistic): + """A supernova statistic. + + This statistic that applies an additive shift M to a supernova's distance + modulus. + """ + + def __init__(self, sacc_tracer: str) -> None: + """Initialize this statistic.""" + super().__init__(parameter_prefix=sacc_tracer) + + self.sacc_tracer = sacc_tracer + self.data_vector: None | DataVector = None + self.a: None | npt.NDArray[np.float64] = None + self.M = parameters.register_new_updatable_parameter( + default_value=SNIA_DEFAULT_M + ) + + def read(self, sacc_data: sacc.Sacc) -> None: + """Read the data for this statistic from the SACC file.""" + # We do not actually need the tracer, but we want to make sure the SACC + # data contains this tracer. + # TODO: remove the work-around when the new version of SACC supporting + # sacc.Sacc.has_tracer is available. + try: + tracer = sacc_data.get_tracer(self.sacc_tracer) + except KeyError as exc: + # Translate to the error type we want + raise ValueError( + f"The SACC file does not contain the MiscTracer {self.sacc_tracer}" + ) from exc + if not isinstance(tracer, MiscTracer): + raise ValueError( + f"The SACC tracer {self.sacc_tracer} is not a " f"MiscTracer" + ) + + data_points = sacc_data.get_data_points( + data_type="supernova_distance_mu", tracers=(self.sacc_tracer,) + ) + z = np.array([dp.get_tag("z") for dp in data_points]) + self.a = 1.0 / (1.0 + z) + self.data_vector = DataVector.from_list([dp.value for dp in data_points]) + self.sacc_indices = np.arange(len(self.data_vector)) + super().read(sacc_data) + + def get_data_vector(self) -> DataVector: + """Return the data vector; raise exception if there is none.""" + assert self.data_vector is not None + return self.data_vector + + def _compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: + """Compute SNIa distance statistic using CCL.""" + ccl_cosmo = tools.get_ccl_cosmology() + prediction = self.M + pyccl.distance_modulus(ccl_cosmo, self.a) + return TheoryVector.create(prediction) diff --git a/firecrown/likelihood/two_point.py b/firecrown/likelihood/two_point.py new file mode 100644 index 000000000..c7266115e --- /dev/null +++ b/firecrown/likelihood/two_point.py @@ -0,0 +1,779 @@ +"""Two point statistic support.""" + +from __future__ import annotations + +import copy +import warnings +from typing import Sequence, TypedDict + +import numpy as np +import numpy.typing as npt +import pyccl +import pyccl.nl_pt +import sacc.windows +import scipy.interpolate + +# firecrown is needed for backward compatibility; remove support for deprecated +# directory structure is removed. +import firecrown # pylint: disable=unused-import # noqa: F401 +from firecrown.likelihood.source import Source, Tracer +from firecrown.likelihood.weak_lensing import ( + WeakLensingFactory, + WeakLensing, +) +from firecrown.likelihood.number_counts import ( + NumberCountsFactory, + NumberCounts, +) +from firecrown.likelihood.statistic import ( + DataVector, + Statistic, + TheoryVector, +) +from firecrown.metadata.two_point import ( + GalaxyMeasuredType, + InferredGalaxyZDist, + TRACER_NAMES_TOTAL, + TracerNames, + TwoPointCells, + Window, + extract_window_function, +) +from firecrown.modeling_tools import ModelingTools +from firecrown.updatable import UpdatableCollection + +# only supported types are here, anything else will throw +# a value error +SACC_DATA_TYPE_TO_CCL_KIND = { + "galaxy_density_cl": "cl", + "galaxy_density_xi": "NN", + "galaxy_shearDensity_cl_e": "cl", + "galaxy_shearDensity_xi_t": "NG", + "galaxy_shear_cl_ee": "cl", + "galaxy_shear_xi_minus": "GG-", + "galaxy_shear_xi_plus": "GG+", + "cmbGalaxy_convergenceDensity_xi": "NN", + "cmbGalaxy_convergenceShear_xi_t": "NG", +} + +ELL_FOR_XI_DEFAULTS = {"minimum": 2, "midpoint": 50, "maximum": 60_000, "n_log": 200} + + +def _ell_for_xi( + *, minimum: int, midpoint: int, maximum: int, n_log: int +) -> npt.NDArray[np.int64]: + """Create an array of ells to sample the power spectrum. + + This is used for for real-space predictions. The result will contain + each integral value from min to mid. Starting from mid, and going up + to max, there will be n_log logarithmically spaced values. + + All values are rounded to the nearest integer. + """ + assert minimum >= 0 + assert minimum < midpoint + assert midpoint < maximum + lower_range = np.linspace(minimum, midpoint - 1, midpoint - minimum) + upper_range = np.logspace(np.log10(midpoint), np.log10(maximum), n_log) + concatenated = np.concatenate((lower_range, upper_range)) + # Round the results to the nearest integer values. + # N.B. the dtype of the result is np.dtype[float64] + return np.unique(np.around(concatenated)).astype(np.int64) + + +def _generate_ell_or_theta(*, minimum, maximum, n, binning="log"): + if binning == "log": + edges = np.logspace(np.log10(minimum), np.log10(maximum), n + 1) + return np.sqrt(edges[1:] * edges[:-1]) + edges = np.linspace(minimum, maximum, n + 1) + return (edges[1:] + edges[:-1]) / 2.0 + + +# @functools.lru_cache(maxsize=128) +def _cached_angular_cl(cosmo, tracers, ells, p_of_k_a=None): + return pyccl.angular_cl( + cosmo, tracers[0], tracers[1], np.array(ells), p_of_k_a=p_of_k_a + ) + + +def make_log_interpolator(x, y): + """Return a function object that does 1D spline interpolation. + + If all the y values are greater than 0, the function + interpolates log(y) as a function of log(x). + Otherwise, the function interpolates y as a function of log(x). + The resulting interpolater will not extrapolate; if called with + an out-of-range argument it will raise a ValueError. + """ + # TODO: There is no code in Firecrown, neither test nor example, that uses + # this in any way. + if np.all(y > 0): + # use log-log interpolation + intp = scipy.interpolate.InterpolatedUnivariateSpline( + np.log(x), np.log(y), ext=2 + ) + return lambda x_, intp=intp: np.exp(intp(np.log(x_))) + # only use log for x + intp = scipy.interpolate.InterpolatedUnivariateSpline(np.log(x), y, ext=2) + return lambda x_, intp=intp: intp(np.log(x_)) + + +def calculate_ells_for_interpolation(w: Window) -> npt.NDArray[np.int64]: + """See _ell_for_xi. + + This method mixes together: + 1. the default parameters in ELL_FOR_XI_DEFAULTS + 2. the first and last values in w. + + and then calls _ell_for_xi with those arguments, returning whatever it + returns. + """ + ell_config = { + **ELL_FOR_XI_DEFAULTS, + "maximum": w.ells[-1], + } + ell_config["minimum"] = max(ell_config["minimum"], w.ells[0]) + return _ell_for_xi(**ell_config) + + +class EllOrThetaConfig(TypedDict): + """A dictionary of options for generating the ell or theta. + + This dictionary contains the minimum, maximum and number of + bins to generate the ell or theta values at which to compute the statistics. + + :param minimum: The start of the binning. + :param maximum: The end of the binning. + :param n: The number of bins. + :param binning: Pass 'log' to get logarithmic spaced bins and 'lin' to get linearly + spaced bins. Default is 'log'. + + """ + + minimum: float + maximum: float + n: int + binning: str + + +def generate_ells_cells(ell_config: EllOrThetaConfig): + """Generate ells or theta values from the configuration dictionary.""" + ells = _generate_ell_or_theta(**ell_config) + Cells = np.zeros_like(ells) + + return ells, Cells + + +def generate_theta_xis(theta_config: EllOrThetaConfig): + """Generate theta and xi values from the configuration dictionary.""" + thetas = _generate_ell_or_theta(**theta_config) + xis = np.zeros_like(thetas) + + return thetas, xis + + +def apply_ells_min_max( + ells: npt.NDArray[np.int64], + Cells: npt.NDArray[np.float64], + indices: None | npt.NDArray[np.int64], + ell_min: None | int, + ell_max: None | int, +) -> tuple[ + npt.NDArray[np.int64], npt.NDArray[np.float64], None | npt.NDArray[np.int64] +]: + """Apply the minimum and maximum ell values to the ells and Cells.""" + if ell_min is not None: + locations = np.where(ells >= ell_min) + ells = ells[locations] + Cells = Cells[locations] + if indices is not None: + indices = indices[locations] + + if ell_max is not None: + locations = np.where(ells <= ell_max) + ells = ells[locations] + Cells = Cells[locations] + if indices is not None: + indices = indices[locations] + + return ells, Cells, indices + + +def apply_theta_min_max( + thetas: npt.NDArray[np.float64], + xis: npt.NDArray[np.float64], + indices: None | npt.NDArray[np.int64], + theta_min: None | float, + theta_max: None | float, +) -> tuple[ + npt.NDArray[np.float64], npt.NDArray[np.float64], None | npt.NDArray[np.int64] +]: + """Apply the minimum and maximum theta values to the thetas and xis.""" + if theta_min is not None: + locations = np.where(thetas >= theta_min) + thetas = thetas[locations] + xis = xis[locations] + if indices is not None: + indices = indices[locations] + + if theta_max is not None: + locations = np.where(thetas <= theta_max) + thetas = thetas[locations] + xis = xis[locations] + if indices is not None: + indices = indices[locations] + + return thetas, xis, indices + + +def use_source_factory( + inferred_galaxy_zdist: InferredGalaxyZDist, + wl_factory: WeakLensingFactory | None = None, + nc_factory: NumberCountsFactory | None = None, +) -> WeakLensing | NumberCounts: + """Apply the factory to the inferred galaxy redshift distribution.""" + source: WeakLensing | NumberCounts + match inferred_galaxy_zdist.measured_type: + case GalaxyMeasuredType.COUNTS: + assert nc_factory is not None + source = nc_factory.create(inferred_galaxy_zdist) + case GalaxyMeasuredType.SHEAR_E | GalaxyMeasuredType.SHEAR_T: + assert wl_factory is not None + source = wl_factory.create(inferred_galaxy_zdist) + case _: + raise ValueError( + f"Measured type {inferred_galaxy_zdist.measured_type} not supported!" + ) + return source + + +def create_sacc(metadata: list[TwoPointCells]): + """Fill the SACC file with the inferred galaxy redshift distributions.""" + sacc_data = sacc.Sacc() + + dv = [] + for cell in metadata: + for inferred_galaxy_zdist in (cell.XY.x, cell.XY.y): + if inferred_galaxy_zdist.bin_name not in sacc_data.tracers: + sacc_data.add_tracer( + "NZ", + inferred_galaxy_zdist.bin_name, + inferred_galaxy_zdist.z, + inferred_galaxy_zdist.dndz, + ) + cells = np.ones_like(cell.ells) + sacc_data.add_ell_cl( + cell.get_sacc_name(), + cell.XY.x.bin_name, + cell.XY.y.bin_name, + ell=cell.ells, + x=cells, + ) + dv.append(cells) + + delta_v = np.concatenate(dv, axis=0) + sacc_data.add_covariance(np.diag(np.ones_like(delta_v))) + + return sacc_data + + +class TwoPoint(Statistic): + """A two-point statistic. + + For example, shear correlation function, galaxy-shear correlation function, etc. + + Parameters + ---------- + sacc_data_type : str + The kind of two-point statistic. This must be a valid SACC data type that + maps to one of the CCL correlation function kinds or a power spectra. + Possible options are + + - galaxy_density_cl : maps to 'cl' (a CCL angular power spectrum) + - galaxy_density_xi : maps to 'gg' (a CCL angular position corr. function) + - galaxy_shearDensity_cl_e : maps to 'cl' (a CCL angular power spectrum) + - galaxy_shearDensity_xi_t : maps to 'gl' (a CCL angular cross-correlation + between position and shear) + - galaxy_shear_cl_ee : maps to 'cl' (a CCL angular power spectrum) + - galaxy_shear_xi_minus : maps to 'l-' (a CCL angular shear corr. + function xi-) + - galaxy_shear_xi_plus : maps to 'l+' (a CCL angular shear corr. + function xi-) + - cmbGalaxy_convergenceDensity_xi : maps to 'gg' (a CCL angular position + corr. function) + - cmbGalaxy_convergenceShear_xi_t : maps to 'gl' (a CCL angular cross- + correlation between position and shear) + + source0 : Source + The first sources needed to compute this statistic. + source1 : Source + The second sources needed to compute this statistic. + ell_or_theta : dict, optional + A dictionary of options for generating the ell or theta values at which + to compute the statistics. This option can be used to have firecrown + generate data without the corresponding 2pt data in the input SACC file. + The options are: + + - minimun : float - The start of the binning. + - maximun : float - The end of the binning. + - n : int - The number of bins. Note that the edges of the bins start + at `min` and end at `max`. The actual bin locations will be at the + (possibly geometric) midpoint of the bin. + - binning : str, optional - Pass 'log' to get logarithmic spaced bins and 'lin' + to get linearly spaced bins. Default is 'log'. + + ell_or_theta_min : float, optional + The minimum ell or theta value to keep. This minimum is applied after + the ell or theta values are read and/or generated. + ell_or_theta_max : float, optional + The maximum ell or theta value to keep. This maximum is applied after + the ell or theta values are read and/or generated. + ell_for_xi : dict, optional + A dictionary of options for making the ell values at which to compute + Cls for use in real-space integrations. The possible keys are: + + - minimum : int, optional - The minimum angular wavenumber to use for + real-space integrations. Default is 2. + - midpoint : int, optional - The midpoint angular wavenumber to use + for real-space integrations. The angular wavenumber samples are + linearly spaced at integers between `minimum` and `midpoint`. Default + is 50. + - maximum : int, optional - The maximum angular wavenumber to use for + real-space integrations. The angular wavenumber samples are + logarithmically spaced between `midpoint` and `maximum`. Default is + 60,000. + - n_log : int, optional - The number of logarithmically spaced angular + wavenumber samples between `mid` and `max`. Default is 200. + + Attributes + ---------- + ccl_kind : str + The CCL correlation function kind or 'cl' for power spectra corresponding + to the SACC data type. + sacc_tracers : 2-tuple of str + A tuple of the SACC tracer names for this 2pt statistic. Set after a + call to read. + ell_or_theta_ : npt.NDArray[np.float64] + The final array of ell/theta values for the statistic. Set after + `compute` is called. + measured_statistic_ : npt.NDArray[np.float64] + The measured value for the statistic. + predicted_statistic_ : npt.NDArray[np.float64] + The final prediction for the statistic. Set after `compute` is called. + + """ + + def __init__( + self, + sacc_data_type: str, + source0: Source, + source1: Source, + *, + ell_for_xi: None | dict[str, int] = None, + ell_or_theta: None | EllOrThetaConfig = None, + ell_or_theta_min: None | float | int = None, + ell_or_theta_max: None | float | int = None, + ) -> None: + super().__init__() + + assert isinstance(source0, Source) + assert isinstance(source1, Source) + + self.sacc_data_type: str = sacc_data_type + self.source0: Source = source0 + self.source1 = source1 + self.ell_for_xi_config: dict[str, int] = copy.deepcopy(ELL_FOR_XI_DEFAULTS) + if ell_for_xi is not None: + self.ell_for_xi_config.update(ell_for_xi) + # What is the difference between the following 3 instance variables? + # ell_or_theta + # _ell_or_theta + # ell_or_theta_ + self.ell_or_theta_config = ell_or_theta + self.ell_or_theta_min = ell_or_theta_min + self.ell_or_theta_max = ell_or_theta_max + self.window: None | Window = None + + self.data_vector: None | DataVector = None + self.theory_vector: None | TheoryVector = None + + self.sacc_tracers: TracerNames + self.ells: None | npt.NDArray[np.int64] = None + self.thetas: None | npt.NDArray[np.float64] = None + self.mean_ells: None | npt.NDArray[np.float64] = None + self.ells_for_xi: None | npt.NDArray[np.int64] = None + + self.cells: dict[TracerNames, npt.NDArray[np.float64]] = {} + if self.sacc_data_type in SACC_DATA_TYPE_TO_CCL_KIND: + self.ccl_kind = SACC_DATA_TYPE_TO_CCL_KIND[self.sacc_data_type] + else: + raise ValueError( + f"The SACC data type {sacc_data_type}'%s' is not " f"supported!" + ) + + @classmethod + def from_metadata_cells( + cls, + metadata: list[TwoPointCells], + wl_factory: WeakLensingFactory | None = None, + nc_factory: NumberCountsFactory | None = None, + ) -> UpdatableCollection[TwoPoint]: + """Create a TwoPoint statistic from a TwoPointCells metadata object.""" + sacc_data = create_sacc(metadata) + + two_point_list = [ + cls( + cell.get_sacc_name(), + use_source_factory( + cell.XY.x, wl_factory=wl_factory, nc_factory=nc_factory + ), + use_source_factory( + cell.XY.y, wl_factory=wl_factory, nc_factory=nc_factory + ), + ) + for cell in metadata + ] + + for two_point in two_point_list: + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + two_point.read(sacc_data) + + return UpdatableCollection(two_point_list) + + def read_ell_cells( + self, sacc_data_type: str, sacc_data: sacc.Sacc, tracers: TracerNames + ) -> ( + None + | tuple[npt.NDArray[np.int64], npt.NDArray[np.float64], npt.NDArray[np.int64]] + ): + """Read and return ell and Cell.""" + ells, Cells = sacc_data.get_ell_cl(sacc_data_type, *tracers, return_cov=False) + # As version 0.13 of sacc, the method get_ell_cl returns the + # ell values and the Cl values in arrays of the same length. + assert len(ells) == len(Cells) + common_length = len(ells) + sacc_indices = None + + if common_length == 0: + return None + sacc_indices = np.atleast_1d(sacc_data.indices(self.sacc_data_type, tracers)) + assert sacc_indices is not None # Needed for mypy + assert len(sacc_indices) == common_length + + return ells, Cells, sacc_indices + + def read_theta_xis( + self, sacc_data_type: str, sacc_data: sacc.Sacc, tracers: TracerNames + ) -> ( + None + | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.int64]] + ): + """Read and return theta and xi.""" + thetas, xis = sacc_data.get_theta_xi(sacc_data_type, *tracers, return_cov=False) + # As version 0.13 of sacc, the method get_theta_xi returns the + # theta values and the xi values in arrays of the same length. + assert len(thetas) == len(xis) + + common_length = len(thetas) + if common_length == 0: + return None + sacc_indices = np.atleast_1d(sacc_data.indices(self.sacc_data_type, tracers)) + assert sacc_indices is not None # Needed for mypy + assert len(sacc_indices) == common_length + return thetas, xis, sacc_indices + + def read(self, sacc_data: sacc.Sacc) -> None: + """Read the data for this statistic from the SACC file. + + :param sacc_data: The data in the sacc format. + """ + self.sacc_tracers = self.initialize_sources(sacc_data) + + if self.ccl_kind == "cl": + self.read_harmonic_space(sacc_data) + else: + self.read_real_space(sacc_data) + + super().read(sacc_data) + + def read_real_space(self, sacc_data: sacc.Sacc): + """Read the data for this statistic from the SACC file.""" + thetas_xis_indices = self.read_theta_xis( + self.sacc_data_type, sacc_data, self.sacc_tracers + ) + # We do not support window functions for real space statistics + if thetas_xis_indices is not None: + thetas, xis, sacc_indices = thetas_xis_indices + if self.ell_or_theta_config is not None: + # If we have data from our construction, and also have data in the + # SACC object, emit a warning and use the information read from the + # SACC object. + warnings.warn( + f"Tracers '{self.sacc_tracers}' have 2pt data and you have " + f"specified `theta` in the configuration. `theta` is being " + f"ignored!", + stacklevel=2, + ) + else: + if self.ell_or_theta_config is None: + # The SACC file has no data points, just a tracer, in this case we + # are building the statistic from scratch. In this case the user + # must have set the dictionary ell_or_theta, containing the + # minimum, maximum and number of bins to generate the ell values. + raise RuntimeError( + f"Tracers '{self.sacc_tracers}' for data type " + f"'{self.sacc_data_type}' " + "have no 2pt data in the SACC file and no input theta values " + "were given!" + ) + thetas, xis = generate_theta_xis(self.ell_or_theta_config) + sacc_indices = None + assert isinstance(self.ell_or_theta_min, (float, type(None))) + assert isinstance(self.ell_or_theta_max, (float, type(None))) + thetas, xis, sacc_indices = apply_theta_min_max( + thetas, xis, sacc_indices, self.ell_or_theta_min, self.ell_or_theta_max + ) + self.ells_for_xi = _ell_for_xi(**self.ell_for_xi_config) + self.thetas = thetas + self.sacc_indices = sacc_indices + self.data_vector = DataVector.create(xis) + + def read_harmonic_space(self, sacc_data: sacc.Sacc): + """Read the data for this statistic from the SACC file.""" + ells_cells_indices = self.read_ell_cells( + self.sacc_data_type, sacc_data, self.sacc_tracers + ) + if ells_cells_indices is not None: + ells, Cells, sacc_indices = ells_cells_indices + if self.ell_or_theta_config is not None: + # If we have data from our construction, and also have data in the + # SACC object, emit a warning and use the information read from the + # SACC object. + warnings.warn( + f"Tracers '{self.sacc_tracers}' have 2pt data and you have " + f"specified `ell` in the configuration. `ell` is being ignored!", + stacklevel=2, + ) + window = extract_window_function(sacc_data, sacc_indices) + if window is not None: + # When using a window function, we do not calculate all Cl's. + # For this reason we have a default set of ells that we use + # to compute Cl's, and we have a set of ells used for + # interpolation. + window.ells_for_interpolation = calculate_ells_for_interpolation(window) + + else: + if self.ell_or_theta_config is None: + # The SACC file has no data points, just a tracer, in this case we + # are building the statistic from scratch. In this case the user + # must have set the dictionary ell_or_theta, containing the + # minimum, maximum and number of bins to generate the ell values. + raise RuntimeError( + f"Tracers '{self.sacc_tracers}' for data type " + f"'{self.sacc_data_type}' " + "have no 2pt data in the SACC file and no input ell values " + "were given!" + ) + ells, Cells = generate_ells_cells(self.ell_or_theta_config) + sacc_indices = None + + # When generating the ells and Cells we do not have a window function + window = None + assert isinstance(self.ell_or_theta_min, (int, type(None))) + assert isinstance(self.ell_or_theta_max, (int, type(None))) + ells, Cells, sacc_indices = apply_ells_min_max( + ells, Cells, sacc_indices, self.ell_or_theta_min, self.ell_or_theta_max + ) + self.ells = ells + self.window = window + self.sacc_indices = sacc_indices + self.data_vector = DataVector.create(Cells) + + def initialize_sources(self, sacc_data: sacc.Sacc) -> TracerNames: + """Initialize this TwoPoint's sources, and return the tracer names.""" + self.source0.read(sacc_data) + if self.source0 is not self.source1: + self.source1.read(sacc_data) + assert self.source0.sacc_tracer is not None + assert self.source1.sacc_tracer is not None + tracers = (self.source0.sacc_tracer, self.source1.sacc_tracer) + return TracerNames(*tracers) + + def get_data_vector(self) -> DataVector: + """Return this statistic's data vector.""" + assert self.data_vector is not None + return self.data_vector + + def compute_theory_vector_real_space(self, tools: ModelingTools) -> TheoryVector: + """Compute a two-point statistic in real space. + + This method computes the two-point statistic in real space. It first computes + the Cl's in harmonic space and then translates them to real space using CCL. + """ + tracers0 = self.source0.get_tracers(tools) + tracers1 = self.source1.get_tracers(tools) + scale0 = self.source0.get_scale() + scale1 = self.source1.get_scale() + + assert self.ccl_kind != "cl" + assert self.thetas is not None + assert self.ells_for_xi is not None + + self.cells = {} + print(self.source0.parameter_prefix, self.source1.parameter_prefix) + cells_for_xi = self.compute_cells( + self.ells_for_xi, scale0, scale1, tools, tracers0, tracers1 + ) + + theory_vector = pyccl.correlation( + tools.get_ccl_cosmology(), + ell=self.ells_for_xi, + C_ell=cells_for_xi, + theta=self.thetas / 60, + type=self.ccl_kind, + ) + assert self.data_vector is not None + return TheoryVector.create(theory_vector) + + def compute_theory_vector_harmonic_space( + self, tools: ModelingTools + ) -> TheoryVector: + """Compute a two-point statistic in harmonic space. + + This method computes the two-point statistic in harmonic space. It computes + either the Cl's at the ells provided by the SACC file or the ells required + for the window function. + """ + tracers0 = self.source0.get_tracers(tools) + tracers1 = self.source1.get_tracers(tools) + scale0 = self.source0.get_scale() + scale1 = self.source1.get_scale() + + assert self.ccl_kind == "cl" + assert self.ells is not None + + if self.window is not None: + # If a window function is provided, we need to compute the Cl's + # for the ells used in the window function. To do this, we will + # first compute the Cl's for the ells used in the interpolation + # and then interpolate the results to the ells used in the window + # function. + assert self.window.ells_for_interpolation is not None + cells_for_interpolation = self.compute_cells( + self.window.ells_for_interpolation, + scale0, + scale1, + tools, + tracers0, + tracers1, + ) + + # TODO: There is no code in Firecrown, neither test nor example, + # that exercises a theory window function in any way. + cell_interpolator = make_log_interpolator( + self.window.ells_for_interpolation, cells_for_interpolation + ) + # Deal with ell=0 and ell=1 + cells_interpolated = np.zeros(self.window.ells.size) + cells_interpolated[2:] = cell_interpolator(self.window.ells[2:]) + + # Here we left multiply the computed Cl's by the window function + # to get the final Cl's. + theory_vector = np.einsum( + "lb, l -> b", + self.window.weights, + cells_interpolated, + ) + # We also compute the mean ell value associated with each bin. + self.mean_ells = np.einsum( + "lb, l -> b", self.window.weights, self.window.ells + ) + + assert self.data_vector is not None + return TheoryVector.create(theory_vector) + + # If we get here, we are working in harmonic space without a window function. + assert self.ells is not None + theory_vector = self.compute_cells( + self.ells, + scale0, + scale1, + tools, + tracers0, + tracers1, + ) + + assert self.data_vector is not None + return TheoryVector.create(theory_vector) + + def _compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: + """Compute a two-point statistic from sources.""" + if self.ccl_kind == "cl": + return self.compute_theory_vector_harmonic_space(tools) + + return self.compute_theory_vector_real_space(tools) + + def compute_cells( + self, + ells: npt.NDArray[np.int64], + scale0: float, + scale1: float, + tools: ModelingTools, + tracers0: Sequence[Tracer], + tracers1: Sequence[Tracer], + ) -> npt.NDArray[np.float64]: + """Compute the power spectrum for the given ells and tracers.""" + for tracer0 in tracers0: + for tracer1 in tracers1: + pk_name = f"{tracer0.field}:{tracer1.field}" + tn = TracerNames(tracer0.tracer_name, tracer1.tracer_name) + if tn in self.cells: + # Already computed this combination, skipping + continue + pk = self.calculate_pk(pk_name, tools, tracer0, tracer1) + + self.cells[tn] = ( + _cached_angular_cl( + tools.get_ccl_cosmology(), + (tracer0.ccl_tracer, tracer1.ccl_tracer), + tuple(ells.tolist()), + p_of_k_a=pk, + ) + * scale0 + * scale1 + ) + # Add up all the contributions to the cells + self.cells[TRACER_NAMES_TOTAL] = np.array(sum(self.cells.values())) + theory_vector = self.cells[TRACER_NAMES_TOTAL] + return theory_vector + + def calculate_pk( + self, pk_name: str, tools: ModelingTools, tracer0: Tracer, tracer1: Tracer + ): + """Return the power spectrum named by pk_name.""" + if tools.has_pk(pk_name): + # Use existing power spectrum + pk = tools.get_pk(pk_name) + elif tracer0.has_pt or tracer1.has_pt: + if not (tracer0.has_pt and tracer1.has_pt): + # Mixture of PT and non-PT tracers + # Create a dummy matter PT tracer for the non-PT part + matter_pt_tracer = pyccl.nl_pt.PTMatterTracer() + if not tracer0.has_pt: + tracer0.pt_tracer = matter_pt_tracer + else: + tracer1.pt_tracer = matter_pt_tracer + # Compute perturbation power spectrum + + pt_calculator = tools.get_pt_calculator() + pk = pt_calculator.get_biased_pk2d( + tracer1=tracer0.pt_tracer, + tracer2=tracer1.pt_tracer, + ) + elif tracer0.has_hm or tracer1.has_hm: + # Compute halo model power spectrum + raise NotImplementedError("Halo model power spectra not supported yet") + else: + raise ValueError(f"No power spectrum for {pk_name} can be found.") + return pk diff --git a/firecrown/likelihood/weak_lensing.py b/firecrown/likelihood/weak_lensing.py new file mode 100644 index 000000000..f312c6a48 --- /dev/null +++ b/firecrown/likelihood/weak_lensing.py @@ -0,0 +1,428 @@ +"""Weak lensing source and systematics.""" + +from __future__ import annotations + +from abc import abstractmethod +from dataclasses import dataclass, replace +from typing import Sequence, final + +import numpy as np +import numpy.typing as npt +import pyccl +import pyccl.nl_pt +import sacc + +# firecrown is needed for backward compatibility; remove support for deprecated +# directory structure is removed. +import firecrown # pylint: disable=unused-import # noqa: F401 +from firecrown import parameters +from firecrown.likelihood.source import ( + SourceGalaxy, + SourceGalaxyArgs, + SourceGalaxyPhotoZShift, + SourceGalaxySelectField, + SourceGalaxySystematic, + Tracer, +) +from firecrown.metadata.two_point import InferredGalaxyZDist +from firecrown.modeling_tools import ModelingTools +from firecrown.parameters import ParamsMap + + +@dataclass(frozen=True) +class WeakLensingArgs(SourceGalaxyArgs): + """Class for weak lensing tracer builder argument.""" + + ia_bias: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None + + has_pt: bool = False + has_hm: bool = False + + ia_pt_c_1: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None + ia_pt_c_d: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None + ia_pt_c_2: None | tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]] = None + + +class WeakLensingSystematic(SourceGalaxySystematic[WeakLensingArgs]): + """Abstract base class for all weak lensing systematics.""" + + @abstractmethod + def apply( + self, tools: ModelingTools, tracer_arg: WeakLensingArgs + ) -> WeakLensingArgs: + """Apply method to include systematics in the tracer_arg.""" + + +class PhotoZShift(SourceGalaxyPhotoZShift[WeakLensingArgs]): + """Photo-z shift systematic.""" + + +class SelectField(SourceGalaxySelectField[WeakLensingArgs]): + """Systematic to select 3D field.""" + + +MULTIPLICATIVE_SHEAR_BIAS_DEFAULT_BIAS = 1.0 + + +class MultiplicativeShearBias(WeakLensingSystematic): + """Multiplicative shear bias systematic. + + This systematic adjusts the `scale_` of a source by `(1 + m)`. + + The following parameters are special Updatable parameters, which means that + they can be updated by the sampler, sacc_tracer is going to be used as a + prefix for the parameters: + + :ivar mult_bias: the multiplicative shear bias parameter. + """ + + def __init__(self, sacc_tracer: str) -> None: + """Create a MultiplicativeShearBias object that uses the named tracer. + + :param sacc_tracer: the name of the tracer in the SACC file. This is used + as a prefix for its parameters. + """ + super().__init__(parameter_prefix=sacc_tracer) + + self.mult_bias = parameters.register_new_updatable_parameter( + default_value=MULTIPLICATIVE_SHEAR_BIAS_DEFAULT_BIAS + ) + + def apply( + self, tools: ModelingTools, tracer_arg: WeakLensingArgs + ) -> WeakLensingArgs: + """Apply multiplicative shear bias to a source. + + The `scale_` of the source is multiplied by `(1 + m)`. + + :param tools: A ModelingTools object. + :param tracer_arg: The WeakLensingArgs to which apply the shear bias. + + :returns: A new WeakLensingArgs object with the shear bias applied. + """ + return replace( + tracer_arg, + scale=tracer_arg.scale * (1.0 + self.mult_bias), + ) + + +LINEAR_ALIGNMENT_DEFAULT_IA_BIAS = 0.5 +LINEAR_ALIGNMENT_DEFAULT_ALPHAZ = 0.0 +LINEAR_ALIGNMENT_DEFAULT_ALPHAG = 1.0 +LINEAR_ALIGNMENT_DEFAULT_Z_PIV = 0.5 + + +class LinearAlignmentSystematic(WeakLensingSystematic): + """Linear alignment systematic. + + This systematic adds a linear intrinsic alignment model systematic + which varies with redshift and the growth function. + + The following parameters are special Updatable parameters, which means that + they can be updated by the sampler, sacc_tracer is going to be used as a + prefix for the parameters: + + :ivar ia_bias: the intrinsic alignment bias parameter. + :ivar alphaz: the redshift dependence of the intrinsic alignment bias. + :ivar alphag: the growth function dependence of the intrinsic alignment bias. + :ivar z_piv: the pivot redshift for the intrinsic alignment bias. + """ + + def __init__(self, sacc_tracer: None | str = None, alphag=1.0): + """Create a LinearAlignmentSystematic object, using the specified tracer name. + + :param sacc_tracer: the name of the tracer in the SACC file. This is used + as a prefix for its parameters. + + """ + super().__init__(parameter_prefix=sacc_tracer) + + self.ia_bias = parameters.register_new_updatable_parameter( + default_value=LINEAR_ALIGNMENT_DEFAULT_IA_BIAS + ) + self.alphaz = parameters.register_new_updatable_parameter( + default_value=LINEAR_ALIGNMENT_DEFAULT_ALPHAZ + ) + self.alphag = parameters.register_new_updatable_parameter( + alphag, default_value=LINEAR_ALIGNMENT_DEFAULT_ALPHAG + ) + self.z_piv = parameters.register_new_updatable_parameter( + default_value=LINEAR_ALIGNMENT_DEFAULT_Z_PIV + ) + + def apply( + self, tools: ModelingTools, tracer_arg: WeakLensingArgs + ) -> WeakLensingArgs: + """Return a new linear alignment systematic. + + This choice is based on the given tracer_arg, in the context of the given + cosmology. + """ + ccl_cosmo = tools.get_ccl_cosmology() + + pref = ((1.0 + tracer_arg.z) / (1.0 + self.z_piv)) ** self.alphaz + pref *= pyccl.growth_factor(ccl_cosmo, 1.0 / (1.0 + tracer_arg.z)) ** ( + self.alphag - 1.0 + ) + + ia_bias_array = pref * self.ia_bias + + return replace( + tracer_arg, + ia_bias=(tracer_arg.z, ia_bias_array), + ) + + +TATT_ALIGNMENT_DEFAULT_IA_A_1 = 1.0 +TATT_ALIGNMENT_DEFAULT_IA_A_2 = 0.5 +TATT_ALIGNMENT_DEFAULT_IA_A_D = 0.5 + + +class TattAlignmentSystematic(WeakLensingSystematic): + """TATT alignment systematic. + + This systematic adds a TATT (nonlinear) intrinsic alignment model systematic. + + The following parameters are special Updatable parameters, which means that + they can be updated by the sampler, sacc_tracer is going to be used as a + prefix for the parameters: + + :ivar ia_a_1: the amplitude of the linear alignment model. + :ivar ia_a_2: the amplitude of the quadratic alignment model. + :ivar ia_a_d: the amplitude of the density-dependent alignment model. + """ + + def __init__(self, sacc_tracer: None | str = None): + """Create a TattAlignmentSystematic object, using the specified tracer name. + + :param sacc_tracer: the name of the tracer in the SACC file. This is used + as a prefix for its parameters. + """ + super().__init__(parameter_prefix=sacc_tracer) + self.ia_a_1 = parameters.register_new_updatable_parameter( + default_value=TATT_ALIGNMENT_DEFAULT_IA_A_1 + ) + self.ia_a_2 = parameters.register_new_updatable_parameter( + default_value=TATT_ALIGNMENT_DEFAULT_IA_A_2 + ) + self.ia_a_d = parameters.register_new_updatable_parameter( + default_value=TATT_ALIGNMENT_DEFAULT_IA_A_D + ) + + def apply( + self, tools: ModelingTools, tracer_arg: WeakLensingArgs + ) -> WeakLensingArgs: + """Return a new linear alignment systematic. + + This choice is based on the given tracer_arg, in the context of the given + cosmology. + """ + ccl_cosmo = tools.get_ccl_cosmology() + z = tracer_arg.z + c_1, c_d, c_2 = pyccl.nl_pt.translate_IA_norm( + ccl_cosmo, + z=z, + a1=self.ia_a_1, + a1delta=self.ia_a_d, + a2=self.ia_a_2, + Om_m2_for_c2=False, + ) + + return replace( + tracer_arg, + has_pt=True, + ia_pt_c_1=(z, c_1), + ia_pt_c_d=(z, c_d), + ia_pt_c_2=(z, c_2), + ) + + +class WeakLensing(SourceGalaxy[WeakLensingArgs]): + """Source class for weak lensing.""" + + def __init__( + self, + *, + sacc_tracer: str, + scale: float = 1.0, + systematics: None | Sequence[SourceGalaxySystematic[WeakLensingArgs]] = None, + ): + """Initialize the WeakLensing object. + + :param sacc_tracer: the name of the tracer in the SACC file. This is used + as a prefix for its parameters. + :param scale: the scale of the source. This is used to scale the shear + power spectrum. + :param systematics: a list of WeakLensingSystematic objects to apply to + this source. + + """ + super().__init__(sacc_tracer=sacc_tracer, systematics=systematics) + + self.sacc_tracer = sacc_tracer + self.scale = scale + self.current_tracer_args: None | WeakLensingArgs = None + self.tracer_args: WeakLensingArgs + + @classmethod + def create_ready( + cls, + inferred_zdist: InferredGalaxyZDist, + systematics: None | list[SourceGalaxySystematic[WeakLensingArgs]] = None, + ) -> WeakLensing: + """Create a WeakLensing object with the given tracer name and scale.""" + obj = cls(sacc_tracer=inferred_zdist.bin_name, systematics=systematics) + # pylint: disable=unexpected-keyword-arg + obj.tracer_args = WeakLensingArgs( + scale=obj.scale, z=inferred_zdist.z, dndz=inferred_zdist.dndz, ia_bias=None + ) + # pylint: enable=unexpected-keyword-arg + + return obj + + @final + def _update_source(self, params: ParamsMap): + """Implementation of Source interface `_update_source`. + + This updates all the contained systematics. + """ + self.systematics.update(params) + + def _read(self, sacc_data: sacc.Sacc) -> None: + """Read the data for this source from the SACC file. + + This sets self.tracer_args, based on the data in `sacc_data` associated with + this object's `sacc_tracer` name. + """ + # pylint: disable=unexpected-keyword-arg + self.tracer_args = WeakLensingArgs( + scale=self.scale, z=np.array([]), dndz=np.array([]), ia_bias=None + ) + # pylint: enable=unexpected-keyword-arg + + super()._read(sacc_data) + + def create_tracers(self, tools: ModelingTools): + """Render a source by applying systematics.""" + ccl_cosmo = tools.get_ccl_cosmology() + tracer_args = self.tracer_args + + assert self.systematics is not None + for systematic in self.systematics: + tracer_args = systematic.apply(tools, tracer_args) + + ccl_wl_tracer = pyccl.WeakLensingTracer( + ccl_cosmo, + dndz=(tracer_args.z, tracer_args.dndz), + ia_bias=tracer_args.ia_bias, + ) + tracers = [Tracer(ccl_wl_tracer, tracer_name="shear", field=tracer_args.field)] + + if tracer_args.has_pt: + ia_pt_tracer = pyccl.nl_pt.PTIntrinsicAlignmentTracer( + c1=tracer_args.ia_pt_c_1, + cdelta=tracer_args.ia_pt_c_d, + c2=tracer_args.ia_pt_c_2, + ) + + ccl_wl_dummy_tracer = pyccl.WeakLensingTracer( + ccl_cosmo, + has_shear=False, + use_A_ia=False, + dndz=(tracer_args.z, tracer_args.dndz), + ia_bias=(tracer_args.z, np.ones_like(tracer_args.z)), + ) + ia_tracer = Tracer( + ccl_wl_dummy_tracer, tracer_name="intrinsic_pt", pt_tracer=ia_pt_tracer + ) + tracers.append(ia_tracer) + + self.current_tracer_args = tracer_args + + return tracers, tracer_args + + def get_scale(self): + """Returns the scales for this Source.""" + assert self.current_tracer_args + return self.current_tracer_args.scale + + +class WeakLensingSystematicFactory: + """Factory class for WeakLensingSystematic objects.""" + + @abstractmethod + def create( + self, inferred_zdist: InferredGalaxyZDist + ) -> SourceGalaxySystematic[WeakLensingArgs]: + """Create a WeakLensingSystematic object with the given tracer name.""" + + +class MultiplicativeShearBiasFactory(WeakLensingSystematicFactory): + """Factory class for MultiplicativeShearBias objects.""" + + def create(self, inferred_zdist: InferredGalaxyZDist) -> MultiplicativeShearBias: + """Create a MultiplicativeShearBias object. + + :param inferred_zdist: The inferred galaxy redshift distribution for + the created MultiplicativeShearBias object. + :return: The created MultiplicativeShearBias object. + """ + return MultiplicativeShearBias(inferred_zdist.bin_name) + + +class TattAlignmentSystematicFactory(WeakLensingSystematicFactory): + """Factory class for TattAlignmentSystematic objects.""" + + def create(self, inferred_zdist: InferredGalaxyZDist) -> TattAlignmentSystematic: + """Create a TattAlignmentSystematic object. + + :param inferred_zdist: The inferred galaxy redshift distribution for + the created TattAlignmentSystematic object. + :return: The created TattAlignmentSystematic object. + """ + return TattAlignmentSystematic(inferred_zdist.bin_name) + + +class PhotoZShiftFactory(WeakLensingSystematicFactory): + """Factory class for PhotoZShift objects.""" + + def create(self, inferred_zdist: InferredGalaxyZDist) -> PhotoZShift: + """Create a PhotoZShift object. + + :param inferred_zdist: The inferred galaxy redshift distribution for + the created PhotoZShift object. + :return: The created PhotoZShift object. + """ + return PhotoZShift(inferred_zdist.bin_name) + + +class WeakLensingFactory: + """Factory class for WeakLensing objects.""" + + def __init__( + self, + per_bin_systematics: list[WeakLensingSystematicFactory], + global_systematics: Sequence[WeakLensingSystematic], + ) -> None: + self.per_bin_systematics: list[WeakLensingSystematicFactory] = ( + per_bin_systematics + ) + self.global_systematics: Sequence[WeakLensingSystematic] = global_systematics + self.cache: dict[int, WeakLensing] = {} + + def create(self, inferred_zdist: InferredGalaxyZDist) -> WeakLensing: + """Create a WeakLensing object with the given tracer name and scale.""" + inferred_zdist_id = id(inferred_zdist) + if inferred_zdist_id in self.cache: + return self.cache[inferred_zdist_id] + + systematics: list[SourceGalaxySystematic[WeakLensingArgs]] = [ + systematic_factory.create(inferred_zdist) + for systematic_factory in self.per_bin_systematics + ] + systematics.extend(self.global_systematics) + + wl = WeakLensing.create_ready(inferred_zdist, systematics) + self.cache[inferred_zdist_id] = wl + + return wl diff --git a/tests/conftest.py b/tests/conftest.py index fd0af10ad..d74b2ea19 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -14,7 +14,7 @@ import pyccl from firecrown.utils import upper_triangle_indices -from firecrown.likelihood.gauss_family.statistic.statistic import TrivialStatistic +from firecrown.likelihood.statistic import TrivialStatistic from firecrown.parameters import ParamsMap from firecrown.connector.mapping import MappingCosmoSIS, mapping_builder from firecrown.modeling_tools import ModelingTools diff --git a/tests/connector/numcosmo/test_numcosmo_mapping.py b/tests/connector/numcosmo/test_numcosmo_mapping.py index 655d4c4c4..a116ede25 100644 --- a/tests/connector/numcosmo/test_numcosmo_mapping.py +++ b/tests/connector/numcosmo/test_numcosmo_mapping.py @@ -6,7 +6,7 @@ from numcosmo_py import Ncm, Nc, GObject from firecrown.likelihood.likelihood import Likelihood, NamedParameters -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian +from firecrown.likelihood.gaussian import ConstGaussian from firecrown.modeling_tools import ModelingTools from firecrown.connector.numcosmo.numcosmo import ( NumCosmoData, diff --git a/tests/likelihood/gauss_family/lkscript_const_gaussian.py b/tests/likelihood/gauss_family/lkscript_const_gaussian.py index 99d592f67..c58be271b 100644 --- a/tests/likelihood/gauss_family/lkscript_const_gaussian.py +++ b/tests/likelihood/gauss_family/lkscript_const_gaussian.py @@ -5,8 +5,8 @@ import sacc -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian -from firecrown.likelihood.gauss_family.statistic.supernova import Supernova +from firecrown.likelihood.gaussian import ConstGaussian +from firecrown.likelihood.supernova import Supernova def build_likelihood(_): diff --git a/tests/likelihood/gauss_family/lkscript_two_point.py b/tests/likelihood/gauss_family/lkscript_two_point.py index fc9e80f80..b8ec49e34 100644 --- a/tests/likelihood/gauss_family/lkscript_two_point.py +++ b/tests/likelihood/gauss_family/lkscript_two_point.py @@ -7,9 +7,9 @@ import sacc from firecrown.likelihood.likelihood import NamedParameters -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian -from firecrown.likelihood.gauss_family.statistic.two_point import TwoPoint -from firecrown.likelihood.gauss_family.statistic.source.number_counts import ( +from firecrown.likelihood.gaussian import ConstGaussian +from firecrown.likelihood.two_point import TwoPoint +from firecrown.likelihood.number_counts import ( NumberCounts, ) diff --git a/tests/likelihood/gauss_family/statistic/source/test_source.py b/tests/likelihood/gauss_family/statistic/source/test_source.py index 772134b08..841565ab1 100644 --- a/tests/likelihood/gauss_family/statistic/source/test_source.py +++ b/tests/likelihood/gauss_family/statistic/source/test_source.py @@ -1,5 +1,5 @@ """ -Tests for the module firecrown.likelihood.gauss_family.statistic.source.source. +Tests for the module firecrown.likelihood.statistic.source. """ import pytest @@ -10,14 +10,14 @@ import sacc from firecrown.modeling_tools import ModelingTools -from firecrown.likelihood.gauss_family.statistic.source.source import ( +from firecrown.likelihood.source import ( Tracer, SourceGalaxy, SourceGalaxyArgs, SourceGalaxySelectField, ) -import firecrown.likelihood.gauss_family.statistic.source.number_counts as nc -import firecrown.likelihood.gauss_family.statistic.source.weak_lensing as wl +import firecrown.likelihood.number_counts as nc +import firecrown.likelihood.weak_lensing as wl from firecrown.metadata.two_point import ( extract_all_tracers, InferredGalaxyZDist, diff --git a/tests/likelihood/gauss_family/statistic/test_binned_cluster_number_counts.py b/tests/likelihood/gauss_family/statistic/test_binned_cluster_number_counts.py index 212d1aacb..9597320e1 100644 --- a/tests/likelihood/gauss_family/statistic/test_binned_cluster_number_counts.py +++ b/tests/likelihood/gauss_family/statistic/test_binned_cluster_number_counts.py @@ -5,12 +5,12 @@ import pytest import pyccl from firecrown.models.cluster.recipes.cluster_recipe import ClusterRecipe -from firecrown.likelihood.gauss_family.statistic.source.source import SourceSystematic +from firecrown.likelihood.source import SourceSystematic from firecrown.modeling_tools import ModelingTools from firecrown.models.cluster.properties import ClusterProperty from firecrown.parameters import ParamsMap from firecrown.models.cluster.abundance import ClusterAbundance -from firecrown.likelihood.gauss_family.statistic.binned_cluster_number_counts import ( +from firecrown.likelihood.binned_cluster_number_counts import ( BinnedClusterNumberCounts, ) diff --git a/tests/likelihood/gauss_family/statistic/test_statistic.py b/tests/likelihood/gauss_family/statistic/test_statistic.py index 8b453342d..0738ab303 100644 --- a/tests/likelihood/gauss_family/statistic/test_statistic.py +++ b/tests/likelihood/gauss_family/statistic/test_statistic.py @@ -1,5 +1,5 @@ """ -Tests for the module firecrown.likelihood.gauss_family.statistic.statistic. +Tests for the module firecrown.likelihood.statistic. """ import numpy as np @@ -7,7 +7,7 @@ import pytest import sacc -import firecrown.likelihood.gauss_family.statistic.statistic as stat +import firecrown.likelihood.statistic as stat from firecrown.modeling_tools import ModelingTools from firecrown.parameters import ParamsMap diff --git a/tests/likelihood/gauss_family/statistic/test_supernova.py b/tests/likelihood/gauss_family/statistic/test_supernova.py index fe6cf7e19..4a433bfc8 100644 --- a/tests/likelihood/gauss_family/statistic/test_supernova.py +++ b/tests/likelihood/gauss_family/statistic/test_supernova.py @@ -6,7 +6,7 @@ import sacc import pyccl -from firecrown.likelihood.gauss_family.statistic.supernova import Supernova +from firecrown.likelihood.supernova import Supernova from firecrown.modeling_tools import ModelingTools from firecrown.parameters import ParamsMap diff --git a/tests/likelihood/gauss_family/statistic/test_two_point.py b/tests/likelihood/gauss_family/statistic/test_two_point.py index 542c717da..6931dca1c 100644 --- a/tests/likelihood/gauss_family/statistic/test_two_point.py +++ b/tests/likelihood/gauss_family/statistic/test_two_point.py @@ -11,13 +11,13 @@ from firecrown.modeling_tools import ModelingTools from firecrown.parameters import ParamsMap -from firecrown.likelihood.gauss_family.statistic.source.number_counts import ( +from firecrown.likelihood.number_counts import ( NumberCounts, ) -from firecrown.likelihood.gauss_family.statistic.source.weak_lensing import ( +from firecrown.likelihood.weak_lensing import ( WeakLensing, ) -from firecrown.likelihood.gauss_family.statistic.two_point import ( +from firecrown.likelihood.two_point import ( _ell_for_xi, TwoPoint, TracerNames, diff --git a/tests/likelihood/gauss_family/test_const_gaussian.py b/tests/likelihood/gauss_family/test_const_gaussian.py index d6a270afe..39e2b0612 100644 --- a/tests/likelihood/gauss_family/test_const_gaussian.py +++ b/tests/likelihood/gauss_family/test_const_gaussian.py @@ -12,9 +12,9 @@ import sacc import firecrown.parameters -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian -from firecrown.likelihood.gauss_family.gauss_family import Statistic -from firecrown.likelihood.gauss_family.statistic.statistic import TrivialStatistic +from firecrown.likelihood.gaussian import ConstGaussian +from firecrown.likelihood.gaussfamily import Statistic +from firecrown.likelihood.statistic import TrivialStatistic from firecrown.modeling_tools import ModelingTools from firecrown.parameters import ( RequiredParameters, diff --git a/tests/likelihood/gauss_family/test_student_t.py b/tests/likelihood/gauss_family/test_student_t.py index b4c9cf74b..d604ab430 100644 --- a/tests/likelihood/gauss_family/test_student_t.py +++ b/tests/likelihood/gauss_family/test_student_t.py @@ -7,8 +7,8 @@ import sacc import firecrown.parameters -from firecrown.likelihood.gauss_family.student_t import StudentT -from firecrown.likelihood.gauss_family.gauss_family import Statistic +from firecrown.likelihood.student_t import StudentT +from firecrown.likelihood.gaussfamily import Statistic from firecrown.modeling_tools import ModelingTools from firecrown.parameters import ( RequiredParameters, diff --git a/tests/metadata/test_metadata_two_point.py b/tests/metadata/test_metadata_two_point.py index 2bd27a8e0..3c277c03a 100644 --- a/tests/metadata/test_metadata_two_point.py +++ b/tests/metadata/test_metadata_two_point.py @@ -28,10 +28,10 @@ type_to_sacc_string_real as real, Window, ) -from firecrown.likelihood.gauss_family.statistic.source.source import SourceGalaxy -import firecrown.likelihood.gauss_family.statistic.source.weak_lensing as wl -import firecrown.likelihood.gauss_family.statistic.source.number_counts as nc -from firecrown.likelihood.gauss_family.statistic.two_point import TwoPoint +from firecrown.likelihood.source import SourceGalaxy +import firecrown.likelihood.weak_lensing as wl +import firecrown.likelihood.number_counts as nc +from firecrown.likelihood.two_point import TwoPoint @pytest.fixture( diff --git a/tests/test_bug398.py b/tests/test_bug398.py index cfd44adc6..bf99b6fbc 100644 --- a/tests/test_bug398.py +++ b/tests/test_bug398.py @@ -6,9 +6,9 @@ from numpy.testing import assert_allclose -import firecrown.likelihood.gauss_family.statistic.source.weak_lensing as wl -from firecrown.likelihood.gauss_family.statistic.two_point import TwoPoint -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian +import firecrown.likelihood.weak_lensing as wl +from firecrown.likelihood.two_point import TwoPoint +from firecrown.likelihood.gaussian import ConstGaussian from firecrown.modeling_tools import ModelingTools from firecrown.likelihood.likelihood import Likelihood, NamedParameters from firecrown.parameters import ParamsMap diff --git a/tests/test_bug_413.py b/tests/test_bug_413.py index 8aea0bfa0..1091a7fcc 100644 --- a/tests/test_bug_413.py +++ b/tests/test_bug_413.py @@ -4,11 +4,11 @@ import pyccl import sacc -from firecrown.likelihood.gauss_family.statistic.source.number_counts import ( +from firecrown.likelihood.number_counts import ( NumberCounts, PTNonLinearBiasSystematic, ) -from firecrown.likelihood.gauss_family.statistic.two_point import TwoPoint +from firecrown.likelihood.two_point import TwoPoint from firecrown.parameters import ParamsMap from firecrown.modeling_tools import ModelingTools diff --git a/tests/test_old_imports.py b/tests/test_old_imports.py new file mode 100644 index 000000000..efd8de392 --- /dev/null +++ b/tests/test_old_imports.py @@ -0,0 +1,86 @@ +""" +Tests to ensure that oldmodule imports still work, and yield modules that contain the +same names as the newmodule imports. +""" + +# pylint: disable=import-outside-toplevel + +import re +import types + +HIDDEN_NAME_PATTERN = re.compile(r"_\w+") + + +def diff_module_names( + oldmodule: types.ModuleType, newmodule: types.ModuleType +) -> set[str]: + """Return a set containing the names specified by the given module.""" + old_names = set(x for x in dir(oldmodule) if not HIDDEN_NAME_PATTERN.match(x)) + new_names = set(x for x in dir(newmodule) if not HIDDEN_NAME_PATTERN.match(x)) + + return old_names ^ new_names + + +def test_import_binned_cluster_number_counts(): + import firecrown.likelihood.gauss_family.statistic.binned_cluster_number_counts as oldmodule # pylint: disable=line-too-long # noqa: E501 + import firecrown.likelihood.binned_cluster_number_counts as newmodule + + assert diff_module_names(oldmodule, newmodule) == set(["warnings"]) + + +def test_import_gaussian(): + import firecrown.likelihood.gauss_family.gaussian as oldmodule + import firecrown.likelihood.gaussian as newmodule + + assert diff_module_names(oldmodule, newmodule) == set(["warnings"]) + + +def test_import_number_counts(): + import firecrown.likelihood.gauss_family.statistic.source.number_counts as oldmodule + import firecrown.likelihood.number_counts as newmodule + + assert diff_module_names(oldmodule, newmodule) == set(["warnings"]) + + +def test_import_source(): + import firecrown.likelihood.gauss_family.statistic.source.source as oldmodule + import firecrown.likelihood.source as newmodule + + assert diff_module_names(oldmodule, newmodule) == set(["warnings"]) + + +def test_import_statistic(): + import firecrown.likelihood.gauss_family.statistic.statistic as oldmodule + import firecrown.likelihood.statistic as newmodule + + # The statistic module contains warnings. + assert diff_module_names(oldmodule, newmodule) == set() + + +def test_import_student_t(): + import firecrown.likelihood.gauss_family.student_t as oldmodule + import firecrown.likelihood.student_t as newmodule + + assert diff_module_names(oldmodule, newmodule) == set(["warnings"]) + + +def test_import_two_point(): + import firecrown.likelihood.gauss_family.statistic.two_point as oldmodule + import firecrown.likelihood.two_point as newmodule + + # The two_point module contains warnings. + assert diff_module_names(oldmodule, newmodule) == set() + + +def test_import_supernova(): + import firecrown.likelihood.gauss_family.statistic.supernova as oldmodule + import firecrown.likelihood.supernova as newmodule + + assert diff_module_names(oldmodule, newmodule) == set(["warnings"]) + + +def test_import_weak_lensing(): + import firecrown.likelihood.gauss_family.statistic.source.weak_lensing as oldmodule + import firecrown.likelihood.weak_lensing as newmodule + + assert diff_module_names(oldmodule, newmodule) == set(["warnings"]) diff --git a/tests/test_pt_systematics.py b/tests/test_pt_systematics.py index 0ed56a854..727005c45 100644 --- a/tests/test_pt_systematics.py +++ b/tests/test_pt_systematics.py @@ -12,14 +12,14 @@ import pyccl.nl_pt as pt import sacc -import firecrown.likelihood.gauss_family.statistic.source.weak_lensing as wl -import firecrown.likelihood.gauss_family.statistic.source.number_counts as nc -from firecrown.likelihood.gauss_family.statistic.two_point import ( +import firecrown.likelihood.weak_lensing as wl +import firecrown.likelihood.number_counts as nc +from firecrown.likelihood.two_point import ( TwoPoint, TracerNames, TRACER_NAMES_TOTAL, ) -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian +from firecrown.likelihood.gaussian import ConstGaussian from firecrown.modeling_tools import ModelingTools from firecrown.parameters import ParamsMap diff --git a/tutorial/_building_blocks_for_the_gaussfamily_likelihoods.qmd b/tutorial/_building_blocks_for_the_gaussfamily_likelihoods.qmd index 8e82000e6..b3279dc0f 100644 --- a/tutorial/_building_blocks_for_the_gaussfamily_likelihoods.qmd +++ b/tutorial/_building_blocks_for_the_gaussfamily_likelihoods.qmd @@ -14,7 +14,7 @@ A given `SACC` object may contain information from observations in many bins, bu All `GaussFamily` likelihoods have an implementation of the `read` method that reads data covariance information from the provided `sacc.SACC` object. These likelihoods use the indices from all of its (possibly many) `Statistic`s to build the covariance matrix for the likelihood. -The class `firecrown.likelihood.gauss_family.two_point.TwoPoint` is a statistic that represents a two-point function. +The class `firecrown.likelihood.two_point.TwoPoint` is a statistic that represents a two-point function. A `TwoPoint` object has two `Source`s, each of which is associated with one or more tracer names. To calculate an autocorrelation, use the same `Source` twice. Each `Source` will produce one or more [`pyccl.Tracer`s](https://ccl.readthedocs.io/en/latest/api/pyccl.tracers.html#pyccl.tracers.Tracer).[^tracer] diff --git a/tutorial/_code_development_hygiene.qmd b/tutorial/_code_development_hygiene.qmd index 1be9dc5f8..45e305ef9 100644 --- a/tutorial/_code_development_hygiene.qmd +++ b/tutorial/_code_development_hygiene.qmd @@ -52,15 +52,13 @@ The following is the set of commands using these tools that are used by the CI s Since a pull request that fails any of these will be automatically rejected by the CI system, we strongly recommend running them before committing and pushing your code. Note that we have not yet completed the cleanup of the whole Firecrown repository, and so we do not yet apply `pylint` to all of the code. We strongly recommend that any new code you write *should* be checked with `pylint` before it is committed to the repository. -We are actively working toward full coverage of the code, and will activate additional checking in the CI system as this work progresses. ```{.bash} -black --check firecrown tests examples -flake8 firecrown tests examples -mypy firecrown tests examples -pylint --rcfile tests/pylintrc tests +black firecrown/ tests/ examples/ +flake8 firecrown/ tests/ examples/ +pyline firecrown pylint --rcfile firecrown/models/pylintrc firecrown/models -pylint firecrown/connector firecrown/*.py firecrown/likelihood/*.py \ - firecrown/likelihood/gauss_family/*.py -python -m pytest -v tests +pylint --rcfile tests/pylintrc tests +mypy -p firecrown -p tests -p examples +python -m pytest --runslow -vv --integration tests ``` diff --git a/tutorial/_high_level_firecrown_classes.qmd b/tutorial/_high_level_firecrown_classes.qmd index 4e8ca3b22..c6db16b89 100644 --- a/tutorial/_high_level_firecrown_classes.qmd +++ b/tutorial/_high_level_firecrown_classes.qmd @@ -10,8 +10,8 @@ A `ModelingTools` object associates a cosmology with a set of objects representi Each of these may be used more than once in the evaluation of the likelihood. This is why they are gathered together in one location: to help assure that different parts of a likelihood calculation that require the same theoretical calculation get the identical theoretical calculation for a given cosmology. -Moreover, we define for `ModelingTools` an abstract class [^abstract-class] for each additional tool that can be used in the likelihood calculation. -Thus, the same tool can have different implementations, and the user can choose which one to use. +Moreover, we define for `ModelingTools` an abstract class [^abstract-class] for each additional tool that can be used in the likelihood calculation. +Thus, the same tool can have different implementations, and the user can choose which one to use. This is intended to partially address the issue of systematic effects, as we can have different implementations of the same tool, each one representing a different systematic effect. For example, we can have different implementations of the halo model, each one including a different effects. This is in contrast to the current implementation, where we would have a single halo model that needs to have its results modified by `Systematic` objects. @@ -31,7 +31,6 @@ calculate_loglike(tools: ModelingTools) -> float The method `read` reads the necessary data (data vectors and covariances) from the provided `sacc.SACC` object. This specifies the data for which we are calculating the likelihood. The method `calculate_loglike` return the (natural) logarithm of the likelihood for the data given the cosmology and models in `tools`. -Gaussian-related likelihoods are subclasses of `firecrown.likelihood.gauss_family.GaussFamily`. +Gaussian-related likelihoods are subclasses of `firecrown.likelihood.gaussfamily.GaussFamily`. Currently-implemented subclasses include `ConstGaussian` and `StudentT`. `ConstGaussian` assumes a Gaussian distribution in which the covariance of the data is constant. - diff --git a/tutorial/development_example.qmd b/tutorial/development_example.qmd index 380e3bd11..a415bef90 100644 --- a/tutorial/development_example.qmd +++ b/tutorial/development_example.qmd @@ -55,8 +55,8 @@ In Firecrown, this means writing a new *likelihood factory function* -- the func ```{python} #| eval: false import os.path -import firecrown.likelihood.gauss_family.statistic.supernova as sn -from firecrown.likelihood.gauss_family.gaussian import ConstGaussian +import firecrown.likelihood.supernova as sn +from firecrown.likelihood.gaussian import ConstGaussian import sacc SACC_FILENAME = os.path.expandvars( \ diff --git a/tutorial/two_point_generators.qmd b/tutorial/two_point_generators.qmd index aca52b308..55ca8be38 100644 --- a/tutorial/two_point_generators.qmd +++ b/tutorial/two_point_generators.qmd @@ -14,7 +14,7 @@ This document outlines how to use `InferredGalaxyZDist` objects to derive the tw Two-point statistics are widely used in cosmological analyses. They provide a statistical summary of the distribution, in either real of harmonic space, of galaxies and the matter density field, or of other observables. -In Firecrown, two-point statistics are represented by the `TwoPoint` class, in the module `firecrown.likelihood.gauss_family.statistic.two_point`. +In Firecrown, two-point statistics are represented by the `TwoPoint` class, in the module `firecrown.likelihood.two_point`. ## Generating the LSST Year 1 Redshift Distribution Bins @@ -104,7 +104,7 @@ These are represented by the classes `LinearAlignmentSystematic`, `Multiplicati Firecrown contains several factories for generating the objects that represent these systematics. The purpose of the factories is to ensure that each of the systematics created by a given factory are consistent. -[^1]: The classes `LinearAlignmentSystematic`, `MultiplicativeShearBias`, and `PhotoZShift` are all found in the module `firecrown.likelihood.gauss_family.statistic.source.weak_lensing`. +[^1]: The classes `LinearAlignmentSystematic`, `MultiplicativeShearBias`, and `PhotoZShift` are all found in the module `firecrown.likelihood.weak_lensing`. In order to generate the two-point statistics we need to define the factories for the weak lensing and number counts systematics. These factories are responsible for generating the systematics that will be applied to the two-point statistics. @@ -112,9 +112,9 @@ The code below defines the factories for the weak lensing and number counts syst ```{python} -import firecrown.likelihood.gauss_family.statistic.source.weak_lensing as wl -import firecrown.likelihood.gauss_family.statistic.source.number_counts as nc -import firecrown.likelihood.gauss_family.statistic.two_point as tp +import firecrown.likelihood.weak_lensing as wl +import firecrown.likelihood.number_counts as nc +import firecrown.likelihood.two_point as tp import warnings # WeakLensing systematics -- global