From d39ba1f3460cf2c4b1acd87e9d713cca5bcf055d Mon Sep 17 00:00:00 2001 From: Tilman Troester Date: Tue, 30 Jan 2024 20:10:17 +0100 Subject: [PATCH] Save predictions to sacc (#349) * Improve type usage, now UpdatableCollection allows the user to specify the type of its contents. * Reorganized theory_vector and data_vector set/get and computation. * Combining new state machine and new methods. * Fixed reset calls in cosmosis connector, must be called after all possible gets in the end of execute(). * Cleaning all quantities computed after updated. * Testing new methods of GaussFamily (and older not tested ones). * Adding optional noise to the realizations and testing it. * Updated documentation. * More tests for Statistics. * Added make_realization_vector to compute a new realization vector. * Added make_realization to generate a new sacc with or without noise. * Simplify checking of RE match * Add COMPUTED state to GaussFamily * Address failure to test line gauss_family:202 * Require 100% coverage on changed lines (testing codecov configuration) * Remove needless call to super().__init__ * Support getting covariance for list of statistics --------- Co-authored-by: Sandro Dias Pinto Vitenti Co-authored-by: Marc Paterno --- codecov.yml | 5 +- examples/des_y1_3x2pt/des_y1_3x2pt_PT.py | 9 +- .../des_y1_3x2pt/des_y1_cosmic_shear_TATT.py | 1 + .../des_y1_cosmic_shear_pk_modifier.py | 3 +- firecrown/connector/cosmosis/likelihood.py | 46 +-- .../likelihood/gauss_family/gauss_family.py | 135 +++++++- firecrown/likelihood/gauss_family/gaussian.py | 10 + .../statistic/binned_cluster_number_counts.py | 2 +- .../statistic/source/number_counts.py | 16 +- .../gauss_family/statistic/source/source.py | 6 +- .../statistic/source/weak_lensing.py | 6 +- .../gauss_family/statistic/statistic.py | 45 ++- .../gauss_family/statistic/supernova.py | 2 +- .../gauss_family/statistic/two_point.py | 18 +- firecrown/likelihood/likelihood.py | 28 +- firecrown/modeling_tools.py | 2 +- firecrown/updatable.py | 29 +- .../statistic/test_cluster_number_counts.py | 3 + .../gauss_family/statistic/test_statistic.py | 70 ++++ .../gauss_family/test_const_gaussian.py | 322 +++++++++++++++++- tests/test_modeling_tools.py | 2 +- tests/test_pt_systematics.py | 4 + tests/test_updatable.py | 25 +- 23 files changed, 677 insertions(+), 112 deletions(-) diff --git a/codecov.yml b/codecov.yml index c9d3d1230..c02ea7cc7 100644 --- a/codecov.yml +++ b/codecov.yml @@ -31,7 +31,7 @@ coverage: patch: default: # basic - target: auto + target: 100% threshold: 0% # advanced @@ -63,4 +63,5 @@ coverage: removed_code_behavior: fully_covered_patch #off, removals_only, adjust_base", github_checks: - annotations: true #,false \ No newline at end of file + annotations: true #,false + diff --git a/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py b/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py index 68fad9d15..5c8386084 100755 --- a/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py +++ b/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py @@ -43,6 +43,9 @@ class CclSetup: @dataclass class CElls: + """A package of related C_ell values, to reduce the number of variables + used in the :meth:`run_likelihood` method.""" + GG: np.ndarray GI: np.ndarray II: np.ndarray @@ -234,6 +237,7 @@ def run_likelihood() -> None: assert likelihood.cov is not None stat0 = likelihood.statistics[0].statistic + assert isinstance(stat0, TwoPoint) # x = likelihood.statistics[0].ell_or_theta_ # y_data = likelihood.statistics[0].measured_statistic_ @@ -243,11 +247,12 @@ def run_likelihood() -> None: print(list(stat0.cells.keys())) - stat2 = likelihood.statistics[2].statistic + stat2 = likelihood.statistics[2].statistic # pylint: disable=no-member assert isinstance(stat2, TwoPoint) print(list(stat2.cells.keys())) - stat3 = likelihood.statistics[3].statistic + stat3 = likelihood.statistics[3].statistic # pylint: disable=no-member + assert isinstance(stat3, TwoPoint) print(list(stat3.cells.keys())) plot_predicted_and_measured_statistics( 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 e918a2676..b7973b1fd 100644 --- a/examples/des_y1_3x2pt/des_y1_cosmic_shear_TATT.py +++ b/examples/des_y1_3x2pt/des_y1_cosmic_shear_TATT.py @@ -158,6 +158,7 @@ def run_likelihood() -> None: print(f"Log-like = {log_like:.1f}") # Plot the predicted and measured statistic + assert isinstance(likelihood, ConstGaussian) two_point_0 = likelihood.statistics[0].statistic assert isinstance(two_point_0, TwoPoint) 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 b739f122d..710e7cf70 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 @@ -156,10 +156,9 @@ def run_likelihood() -> None: print(f"Log-like = {log_like:.1f}") # Plot the predicted and measured statistic + assert isinstance(likelihood, ConstGaussian) two_point_0 = likelihood.statistics[0].statistic assert isinstance(two_point_0, TwoPoint) - - assert isinstance(likelihood, ConstGaussian) assert likelihood.cov is not None # Predict CCL Cl diff --git a/firecrown/connector/cosmosis/likelihood.py b/firecrown/connector/cosmosis/likelihood.py index 3a027c5e2..00e1ef91c 100644 --- a/firecrown/connector/cosmosis/likelihood.py +++ b/firecrown/connector/cosmosis/likelihood.py @@ -41,9 +41,6 @@ class FirecrownLikelihood: :param config: current CosmoSIS datablock """ - likelihood: Likelihood - map: MappingCosmoSIS - def __init__(self, config: cosmosis.datablock): """Create the FirecrownLikelihood object from the given configuration.""" likelihood_source = config.get_string(option_section, "likelihood_source", "") @@ -61,6 +58,7 @@ def __init__(self, config: cosmosis.datablock): self.firecrown_module_name = option_section self.sampling_sections = sections + self.likelihood: Likelihood try: self.likelihood, self.tools = load_likelihood( likelihood_source, build_parameters @@ -70,7 +68,7 @@ def __init__(self, config: cosmosis.datablock): print(f"The Firecrown likelihood needs a required parameter: {err}") print("*" * 30) raise - self.map = mapping_builder( + self.map: MappingCosmoSIS = mapping_builder( input_style="CosmoSIS", require_nonlinear_pk=require_nonlinear_pk ) @@ -127,23 +125,31 @@ def execute(self, sample: cosmosis.datablock) -> int: for section, name, val in derived_params_collection: sample.put(section, name, val) - self.likelihood.reset() - self.tools.reset() + if not isinstance(self.likelihood, GaussFamily): + self.likelihood.reset() + self.tools.reset() + return 0 - # Save concatenated data vector and inverse covariance to enable support + # If we get here, we have a GaussFamily likelihood, and we need to + # save concatenated data vector and inverse covariance to enable support # for the CosmoSIS Fisher sampler. This can only work for likelihoods # that have these quantities. Currently, this is only GaussFamily. - if isinstance(self.likelihood, GaussFamily): - sample.put( - "data_vector", "firecrown_theory", self.likelihood.predicted_data_vector - ) - sample.put( - "data_vector", "firecrown_data", self.likelihood.measured_data_vector - ) - sample.put( - "data_vector", "firecrown_inverse_covariance", self.likelihood.inv_cov - ) + sample.put( + "data_vector", + "firecrown_theory", + self.likelihood.get_theory_vector(), + ) + sample.put( + "data_vector", + "firecrown_data", + self.likelihood.get_data_vector(), + ) + sample.put( + "data_vector", + "firecrown_inverse_covariance", + self.likelihood.inv_cov, + ) # Write out theory and data vectors to the data block the ease # debugging. @@ -164,14 +170,16 @@ def execute(self, sample: cosmosis.datablock) -> int: sample.put( "data_vector", f"theory_{stat.sacc_data_type}_{tracer}", - stat.predicted_statistic_, + stat.get_theory_vector(), ) sample.put( "data_vector", f"data_{stat.sacc_data_type}_{tracer}", - stat.measured_statistic_, + stat.get_data_vector(), ) + self.likelihood.reset() + self.tools.reset() return 0 def form_error_message(self, exc: MissingSamplerParameterError) -> str: diff --git a/firecrown/likelihood/gauss_family/gauss_family.py b/firecrown/likelihood/gauss_family/gauss_family.py index 84179a0e0..9765ae5eb 100644 --- a/firecrown/likelihood/gauss_family/gauss_family.py +++ b/firecrown/likelihood/gauss_family/gauss_family.py @@ -8,9 +8,10 @@ """ from __future__ import annotations + from enum import Enum from functools import wraps -from typing import List, Optional, Tuple, Sequence, Callable, Union, TypeVar +from typing import List, Optional, Tuple, Sequence, Callable, Union, TypeVar, Dict from typing import final import warnings @@ -34,6 +35,7 @@ class State(Enum): INITIALIZED = 1 READY = 2 UPDATED = 3 + COMPUTED = 4 T = TypeVar("T") @@ -114,7 +116,10 @@ class GaussFamily(Likelihood): 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`. + :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. @@ -129,12 +134,23 @@ def __init__( self.state: State = State.INITIALIZED if len(statistics) == 0: raise ValueError("GaussFamily requires at least one statistic") - self.statistics: UpdatableCollection = UpdatableCollection( + + 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: Optional[npt.NDArray[np.float64]] = None self.cholesky: Optional[npt.NDArray[np.float64]] = None self.inv_cov: Optional[npt.NDArray[np.float64]] = None + self.cov_index_map: Optional[Dict[int, int]] = None + self.theory_vector: Optional[npt.NDArray[np.double]] = None + self.data_vector: Optional[npt.NDArray[np.double]] = None @enforce_states( initial=State.READY, @@ -149,7 +165,7 @@ def _update(self, _: ParamsMap) -> None: method.""" @enforce_states( - initial=State.UPDATED, + initial=[State.UPDATED, State.COMPUTED], terminal=State.READY, failure_message="update() must be called before reset()", ) @@ -159,6 +175,7 @@ def _reset(self) -> None: 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, @@ -176,17 +193,28 @@ def read(self, sacc_data: sacc.Sacc) -> None: 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_list = [s.statistic.sacc_indices.copy() for s in self.statistics] 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) @@ -196,27 +224,48 @@ def read(self, sacc_data: sacc.Sacc) -> None: failure_message="read() must be called before get_cov()", ) @final - def get_cov(self) -> npt.NDArray[np.float64]: - """Gets the current covariance matrix.""" + def get_cov( + self, statistic: Union[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 return. If not specified, return the covariance of all + statistics. + """ assert self.cov is not None - return self.cov + 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], + 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 and concatenate in the right order.""" - data_vector_list: List[npt.NDArray[np.float64]] = [ - stat.get_data_vector() for stat in self.statistics - ] - return np.concatenate(data_vector_list) + assert self.data_vector is not None + return self.data_vector @final @enforce_states( initial=State.UPDATED, + terminal=State.COMPUTED, failure_message="update() must be called before compute_theory_vector()", ) def compute_theory_vector(self, tools: ModelingTools) -> npt.NDArray[np.float64]: @@ -227,7 +276,22 @@ def compute_theory_vector(self, tools: ModelingTools) -> npt.NDArray[np.float64] theory_vector_list: List[npt.NDArray[np.float64]] = [ stat.compute_theory_vector(tools) for stat in self.statistics ] - return np.concatenate(theory_vector_list) + 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 and concatenate 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( @@ -249,7 +313,8 @@ def compute( @final @enforce_states( - initial=State.UPDATED, + initial=[State.UPDATED, State.COMPUTED], + terminal=State.COMPUTED, failure_message="update() must be called before compute_chisq()", ) def compute_chisq(self, tools: ModelingTools) -> float: @@ -266,9 +331,43 @@ def compute_chisq(self, tools: ModelingTools) -> float: assert len(data_vector) == len(theory_vector) residuals = data_vector - theory_vector - self.predicted_data_vector: npt.NDArray[np.float64] = theory_vector - self.measured_data_vector: npt.NDArray[np.float64] = data_vector - x = scipy.linalg.solve_triangular(self.cholesky, residuals, lower=True) chisq = np.dot(x, x) return chisq + + @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: + new_sacc = sacc_data.copy() + + sacc_indices_list = [] + for stat in self.statistics: + assert stat.statistic.sacc_indices is not None + sacc_indices_list.append(stat.statistic.sacc_indices.copy()) + + sacc_indices = np.concatenate(sacc_indices_list) + + if add_noise: + new_data_vector = self.make_realization_vector() + else: + new_data_vector = self.get_theory_vector() + + assert len(sacc_indices) == len(new_data_vector) + + if strict: + if set(sacc_indices.tolist()) != set(sacc_data.indices()): + raise RuntimeError( + "The predicted data does not cover all the data in the " + "sacc object. To write only the calculated predictions, " + "set strict=False." + ) + + for prediction_idx, sacc_idx in enumerate(sacc_indices): + new_sacc.data[sacc_idx].value = new_data_vector[prediction_idx] + + return new_sacc diff --git a/firecrown/likelihood/gauss_family/gaussian.py b/firecrown/likelihood/gauss_family/gaussian.py index af25bab48..852807677 100644 --- a/firecrown/likelihood/gauss_family/gaussian.py +++ b/firecrown/likelihood/gauss_family/gaussian.py @@ -3,6 +3,7 @@ """ from __future__ import annotations +import numpy as np from .gauss_family import GaussFamily from ...modeling_tools import ModelingTools @@ -15,3 +16,12 @@ def compute_loglike(self, tools: ModelingTools): """Compute the log-likelihood.""" return -0.5 * self.compute_chisq(tools) + + def make_realization_vector(self) -> np.ndarray: + 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/gauss_family/statistic/binned_cluster_number_counts.py b/firecrown/likelihood/gauss_family/statistic/binned_cluster_number_counts.py index 234657d28..f3791cac2 100644 --- a/firecrown/likelihood/gauss_family/statistic/binned_cluster_number_counts.py +++ b/firecrown/likelihood/gauss_family/statistic/binned_cluster_number_counts.py @@ -78,7 +78,7 @@ def get_data_vector(self) -> DataVector: assert self.data_vector is not None return self.data_vector - def compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: + def _compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: assert tools.cluster_abundance is not None theory_vector_list: List[float] = [] diff --git a/firecrown/likelihood/gauss_family/statistic/source/number_counts.py b/firecrown/likelihood/gauss_family/statistic/source/number_counts.py index f751c8435..fd17b5df4 100644 --- a/firecrown/likelihood/gauss_family/statistic/source/number_counts.py +++ b/firecrown/likelihood/gauss_family/statistic/source/number_counts.py @@ -3,7 +3,7 @@ """ from __future__ import annotations -from typing import List, Union, Tuple, Optional, final +from typing import List, Tuple, Optional, final from dataclasses import dataclass, replace from abc import abstractmethod @@ -263,9 +263,6 @@ def apply( class NumberCounts(SourceGalaxy[NumberCountsArgs]): """Source class for number counts.""" - systematics: UpdatableCollection - tracer_args: NumberCountsArgs - def __init__( self, *, @@ -273,11 +270,7 @@ def __init__( has_rsd: bool = False, derived_scale: bool = False, scale: float = 1.0, - systematics: Optional[ - List[ - Union[NumberCountsSystematic, SourceGalaxySystematic[NumberCountsArgs]] - ] - ] = None, + systematics: Optional[List[SourceGalaxySystematic[NumberCountsArgs]]] = None, ): """Initialize the NumberCounts object. @@ -296,9 +289,12 @@ def __init__( self.derived_scale = derived_scale self.bias = parameters.register_new_updatable_parameter() - self.systematics = UpdatableCollection(systematics) + self.systematics: UpdatableCollection[ + SourceGalaxySystematic[NumberCountsArgs] + ] = UpdatableCollection(systematics) self.scale = scale self.current_tracer_args: Optional[NumberCountsArgs] = None + self.tracer_args: NumberCountsArgs @final def _update_source(self, params: ParamsMap): diff --git a/firecrown/likelihood/gauss_family/statistic/source/source.py b/firecrown/likelihood/gauss_family/statistic/source/source.py index 324c95f28..a26dac3a1 100644 --- a/firecrown/likelihood/gauss_family/statistic/source/source.py +++ b/firecrown/likelihood/gauss_family/statistic/source/source.py @@ -213,7 +213,7 @@ class SourceGalaxyPhotoZShift( :ivar delta_z: the photo-z shift. """ - def __init__(self, sacc_tracer: str): + 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 @@ -280,7 +280,9 @@ def __init__( self.sacc_tracer = sacc_tracer self.current_tracer_args: Optional[_SourceGalaxyArgsT] = None - self.systematics: UpdatableCollection = UpdatableCollection(systematics) + self.systematics: UpdatableCollection[SourceGalaxySystematic] = ( + UpdatableCollection(systematics) + ) self.tracer_args: _SourceGalaxyArgsT def _read(self, sacc_data: sacc.Sacc): diff --git a/firecrown/likelihood/gauss_family/statistic/source/weak_lensing.py b/firecrown/likelihood/gauss_family/statistic/source/weak_lensing.py index 5d225cf65..c2ee9efba 100644 --- a/firecrown/likelihood/gauss_family/statistic/source/weak_lensing.py +++ b/firecrown/likelihood/gauss_family/statistic/source/weak_lensing.py @@ -3,7 +3,7 @@ """ from __future__ import annotations -from typing import List, Tuple, Optional, Union, final +from typing import List, Tuple, Optional, final from dataclasses import dataclass, replace from abc import abstractmethod @@ -214,9 +214,7 @@ def __init__( *, sacc_tracer: str, scale: float = 1.0, - systematics: Optional[ - List[Union[WeakLensingSystematic, SourceGalaxySystematic[WeakLensingArgs]]] - ] = None, + systematics: Optional[List[SourceGalaxySystematic[WeakLensingArgs]]] = None, ): """Initialize the WeakLensing object. diff --git a/firecrown/likelihood/gauss_family/statistic/statistic.py b/firecrown/likelihood/gauss_family/statistic/statistic.py index 7ec4fe59d..9552f2dc4 100644 --- a/firecrown/likelihood/gauss_family/statistic/statistic.py +++ b/firecrown/likelihood/gauss_family/statistic/statistic.py @@ -120,6 +120,8 @@ def __init__(self, parameter_prefix: Optional[str] = None): super().__init__(parameter_prefix=parameter_prefix) self.sacc_indices: Optional[npt.NDArray[np.int64]] self.ready = False + self.computed_theory_vector = False + self.theory_vector: Optional[TheoryVector] = None def read(self, _: sacc.Sacc) -> None: """Read the data for this statistic from the SACC data, and mark it @@ -138,13 +140,45 @@ def read(self, _: sacc.Sacc) -> None: f"of length 0; the length must be positive" ) + def _reset(self): + """Reset this statistic, 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.""" - @abstractmethod + @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): @@ -156,6 +190,7 @@ 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: @@ -223,11 +258,6 @@ def read(self, sacc_data: sacc.Sacc): self.sacc_indices = np.arange(len(self.data_vector)) super().read(sacc_data) - @final - def _reset(self): - """Reset this statistic.""" - self.computed_theory_vector = False - @final def _required_parameters(self) -> RequiredParameters: """Return an empty RequiredParameters.""" @@ -243,7 +273,6 @@ def get_data_vector(self) -> DataVector: assert self.data_vector is not None return self.data_vector - def compute_theory_vector(self, _: ModelingTools) -> TheoryVector: + def _compute_theory_vector(self, _: ModelingTools) -> TheoryVector: """Return a fixed theory vector.""" - self.computed_theory_vector = True return TheoryVector.from_list([self.mean] * self.count) diff --git a/firecrown/likelihood/gauss_family/statistic/supernova.py b/firecrown/likelihood/gauss_family/statistic/supernova.py index 0c1a9e8aa..63ce81261 100644 --- a/firecrown/likelihood/gauss_family/statistic/supernova.py +++ b/firecrown/likelihood/gauss_family/statistic/supernova.py @@ -62,7 +62,7 @@ def get_data_vector(self) -> DataVector: assert self.data_vector is not None return self.data_vector - def compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: + def _compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: """Compute SNIa distance statistic using CCL.""" ccl_cosmo = tools.get_ccl_cosmology() diff --git a/firecrown/likelihood/gauss_family/statistic/two_point.py b/firecrown/likelihood/gauss_family/statistic/two_point.py index 44005edd8..74c99e444 100644 --- a/firecrown/likelihood/gauss_family/statistic/two_point.py +++ b/firecrown/likelihood/gauss_family/statistic/two_point.py @@ -2,7 +2,7 @@ """ from __future__ import annotations -from typing import Dict, Tuple, Optional, final, Union +from typing import Dict, Tuple, Optional, Union import copy import functools import warnings @@ -207,8 +207,6 @@ def __init__( self.data_vector: Optional[DataVector] = None self.theory_vector: Optional[TheoryVector] = None self._ell_or_theta: Optional[npt.NDArray[np.float64]] = None - self.predicted_statistic_: Optional[TheoryVector] = None - self.measured_statistic_: Optional[DataVector] = None self.ell_or_theta_: Optional[npt.NDArray[np.float64]] = None self.sacc_tracers: Tuple[str, str] @@ -222,13 +220,6 @@ def __init__( f"The SACC data type {sacc_data_type}'%s' is not " f"supported!" ) - @final - def _reset(self) -> None: - """Prepared to be called again for a new cosmology.""" - # TODO: Why is self.predicted_statistic_ not re-set to None here? - # If we do that, then the CosmoSIS module fails -- because this data - # is accessed from that code. - def read(self, sacc_data: sacc.Sacc) -> None: """Read the data for this statistic from the SACC file. @@ -298,7 +289,7 @@ def read(self, sacc_data: sacc.Sacc) -> None: # I don't think we need these copies, but being safe here. self._ell_or_theta = _ell_or_theta.copy() self.data_vector = DataVector.create(_stat) - self.measured_statistic_ = self.data_vector + self.data_vector = self.data_vector self.sacc_tracers = tracers super().read(sacc_data) @@ -331,7 +322,7 @@ def get_data_vector(self) -> DataVector: assert self.data_vector is not None return self.data_vector - def compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: + def _compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: """Compute a two-point statistic from sources.""" assert self._ell_or_theta is not None @@ -410,10 +401,7 @@ def compute_theory_vector(self, tools: ModelingTools) -> TheoryVector: "lb, l -> b", self.theory_window_function.weight, ell ) - self.predicted_statistic_ = TheoryVector.create(theory_vector) - assert self.data_vector is not None - return TheoryVector.create(theory_vector) def calculate_pk( diff --git a/firecrown/likelihood/likelihood.py b/firecrown/likelihood/likelihood.py index f185ce9a9..5d2c485c0 100644 --- a/firecrown/likelihood/likelihood.py +++ b/firecrown/likelihood/likelihood.py @@ -65,7 +65,7 @@ import sacc -from ..updatable import Updatable, UpdatableCollection +from ..updatable import Updatable from ..modeling_tools import ModelingTools @@ -83,15 +83,31 @@ def __init__(self, parameter_prefix: Optional[str] = None) -> None: """Default initialization for a base Likelihood object.""" super().__init__(parameter_prefix=parameter_prefix) - self.predicted_data_vector: Optional[npt.NDArray[np.double]] = None - self.measured_data_vector: Optional[npt.NDArray[np.double]] = None - self.inv_cov: Optional[npt.NDArray[np.double]] = None - self.statistics: UpdatableCollection = UpdatableCollection() - @abstractmethod def read(self, sacc_data: sacc.Sacc) -> None: """Read the covariance matrix for this likelihood from the SACC file.""" + def make_realization_vector(self) -> npt.NDArray[np.float64]: + """Create a new realization of the model using the previously computed + theory vector and covariance matrix. + """ + raise NotImplementedError( + "This class does not implement make_realization_vector." + ) + + def make_realization( + self, sacc_data: sacc.Sacc, add_noise: bool = True, strict: bool = True + ) -> sacc.Sacc: + """Create a new realization of the model using the previously computed + theory vector and covariance matrix. + + :param sacc_data: The SACC data object containing the covariance matrix + :param add_noise: If True, add noise to the realization. If False, return + only the theory vector. + :param strict: If True, check that the indices of the realization cover + all the indices of the SACC data object. + """ + @abstractmethod def compute_loglike(self, tools: ModelingTools) -> float: """Compute the log-likelihood of generic CCL data.""" diff --git a/firecrown/modeling_tools.py b/firecrown/modeling_tools.py index bee5af689..6167d8599 100644 --- a/firecrown/modeling_tools.py +++ b/firecrown/modeling_tools.py @@ -77,7 +77,7 @@ def prepare(self, ccl_cosmo: pyccl.Cosmology) -> None: """ if not self.is_updated(): - raise RuntimeError("ModelingTools has not been prepared") + raise RuntimeError("ModelingTools has not been updated.") if self._prepared: raise RuntimeError("ModelingTools has already been prepared") diff --git a/firecrown/updatable.py b/firecrown/updatable.py index 33004f675..18c06eed1 100644 --- a/firecrown/updatable.py +++ b/firecrown/updatable.py @@ -14,7 +14,18 @@ """ from __future__ import annotations -from typing import final, Dict, Optional, Any, List, Union, Iterable +from typing import ( + Any, + cast, + Dict, + final, + Generic, + Iterable, + List, + Optional, + TypeVar, + Union, +) from abc import ABC from collections import UserList from .parameters import ( @@ -299,7 +310,10 @@ def _get_derived_parameters(self) -> DerivedParameterCollection: return DerivedParameterCollection([]) -class UpdatableCollection(UserList[Any]): +T = TypeVar("T", bound=Updatable) + + +class UpdatableCollection(UserList[T], Generic[T]): """UpdatableCollection is a list of Updatable objects and is itself supports :meth:`update` and :meth:`reset` (although it does not inherit from :class:`Updatable`). @@ -309,7 +323,7 @@ class UpdatableCollection(UserList[Any]): updated. """ - def __init__(self, iterable: Optional[Iterable[Any]] = None) -> None: + def __init__(self, iterable: Optional[Iterable[T]] = None) -> None: """Initialize the UpdatableCollection from the supplied iterable. If the iterable contains any object that is not Updatable, a TypeError @@ -378,7 +392,7 @@ def get_derived_parameters(self) -> Optional[DerivedParameterCollection]: return derived_parameters return None - def append(self, item: Updatable) -> None: + def append(self, item: T) -> None: """Append the given item to self. If the item is not Updatable a TypeError is raised. @@ -391,10 +405,11 @@ def append(self, item: Updatable) -> None: ) super().append(item) - def __setitem__(self, key, value) -> None: + def __setitem__(self, key, value): """Set self[key] to value; raise TypeError if Value is not Updatable.""" if not isinstance(value, Updatable): raise TypeError( - "Values inserted into an UpdatableCollection must be Updatable" + "Only updatable items can be appended to an UpdatableCollection" ) - super().__setitem__(key, value) + + super().__setitem__(key, cast(T, value)) diff --git a/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py b/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py index aaf3f3c9c..212d1aacb 100644 --- a/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py +++ b/tests/likelihood/gauss_family/statistic/test_cluster_number_counts.py @@ -103,12 +103,14 @@ def test_compute_theory_vector(cluster_sacc_data: sacc.Sacc): bnc = BinnedClusterNumberCounts(ClusterProperty.COUNTS, "my_survey", recipe) bnc.read(cluster_sacc_data) + bnc.update(params) tv = bnc.compute_theory_vector(tools) assert tv is not None assert len(tv) == 2 bnc = BinnedClusterNumberCounts(ClusterProperty.MASS, "my_survey", recipe) bnc.read(cluster_sacc_data) + bnc.update(params) tv = bnc.compute_theory_vector(tools) assert tv is not None assert len(tv) == 2 @@ -117,6 +119,7 @@ def test_compute_theory_vector(cluster_sacc_data: sacc.Sacc): (ClusterProperty.COUNTS | ClusterProperty.MASS), "my_survey", recipe ) bnc.read(cluster_sacc_data) + bnc.update(params) tv = bnc.compute_theory_vector(tools) assert tv is not None assert len(tv) == 4 diff --git a/tests/likelihood/gauss_family/statistic/test_statistic.py b/tests/likelihood/gauss_family/statistic/test_statistic.py index 904954e0d..4f10324c8 100644 --- a/tests/likelihood/gauss_family/statistic/test_statistic.py +++ b/tests/likelihood/gauss_family/statistic/test_statistic.py @@ -4,10 +4,13 @@ from typing import List import numpy as np +from numpy.testing import assert_allclose import pytest import sacc import firecrown.likelihood.gauss_family.statistic.statistic as stat +from firecrown.modeling_tools import ModelingTools +from firecrown.parameters import ParamsMap VECTOR_CLASSES = (stat.TheoryVector, stat.DataVector) @@ -105,3 +108,70 @@ def test_guarded_statistic_get_data_before_read(trivial_stats): ): g = stat.GuardedStatistic(s) _ = g.get_data_vector() + + +def test_statistic_get_data_vector_before_read(): + s = stat.TrivialStatistic() + with pytest.raises(AssertionError): + _ = s.get_data_vector() + + +def test_statistic_get_data_vector_after_read(sacc_data_for_trivial_stat): + s = stat.TrivialStatistic() + s.read(sacc_data_for_trivial_stat) + assert_allclose(s.get_data_vector(), [1.0, 4.0, -3.0]) + + +def test_statistic_get_theory_vector_before_compute(): + s = stat.TrivialStatistic() + with pytest.raises( + RuntimeError, match="The theory for statistic .* has not been computed yet\\." + ): + _ = s.get_theory_vector() + + +def test_statistic_get_theory_vector_after_compute(sacc_data_for_trivial_stat): + s = stat.TrivialStatistic() + s.read(sacc_data_for_trivial_stat) + s.update(ParamsMap(mean=10.5)) + s.compute_theory_vector(ModelingTools()) + assert_allclose(s.get_theory_vector(), [10.5, 10.5, 10.5]) + + +def test_statistic_get_theory_vector_after_reset(sacc_data_for_trivial_stat): + s = stat.TrivialStatistic() + s.read(sacc_data_for_trivial_stat) + s.update(ParamsMap(mean=10.5)) + s.compute_theory_vector(ModelingTools()) + s.reset() + with pytest.raises( + RuntimeError, match="The theory for statistic .* has not been computed yet\\." + ): + _ = s.get_theory_vector() + + +def test_statistic_compute_theory_vector_before_read(): + s = stat.TrivialStatistic() + with pytest.raises( + RuntimeError, + match="The statistic .* has not been updated with parameters\\.", + ): + s.compute_theory_vector(ModelingTools()) + + +def test_statistic_compute_theory_vector_before_update(sacc_data_for_trivial_stat): + s = stat.TrivialStatistic() + s.read(sacc_data_for_trivial_stat) + with pytest.raises( + RuntimeError, + match="The statistic .* has not been updated with parameters\\.", + ): + s.compute_theory_vector(ModelingTools()) + + +def test_statistic_compute_theory_vector_after_update(sacc_data_for_trivial_stat): + s = stat.TrivialStatistic() + s.read(sacc_data_for_trivial_stat) + s.update(ParamsMap(mean=10.5)) + + assert_allclose(s.compute_theory_vector(ModelingTools()), [10.5, 10.5, 10.5]) diff --git a/tests/likelihood/gauss_family/test_const_gaussian.py b/tests/likelihood/gauss_family/test_const_gaussian.py index c0f9f9974..98cce98fc 100644 --- a/tests/likelihood/gauss_family/test_const_gaussian.py +++ b/tests/likelihood/gauss_family/test_const_gaussian.py @@ -1,14 +1,21 @@ """Unit testsing for ConstGaussian """ +import re +from typing import Tuple + import pytest import numpy as np +import numpy.typing as npt +from numpy.testing import assert_allclose +from scipy.stats import chi2 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.modeling_tools import ModelingTools from firecrown.parameters import ( RequiredParameters, @@ -16,6 +23,17 @@ ) +class StatisticWithoutIndices(TrivialStatistic): + """This is a statistic that has no indices when read. It is only for + testing. + """ + + def read(self, sacc_data: sacc.Sacc) -> None: + """Read the TrivialStatistic data, then nullify the sacc_indices.""" + super().read(sacc_data) + self.sacc_indices = None + + def test_require_nonempty_statistics(): with pytest.raises(ValueError): _ = ConstGaussian(statistics=[]) @@ -27,18 +45,161 @@ def test_get_cov_fails_before_read(trivial_stats): _ = likelihood.get_cov() +def test_read_with_wrong_statistic_fails(sacc_data_for_trivial_stat): + # Make the first statistic defective. + defective_stat = StatisticWithoutIndices() + likelihood = ConstGaussian(statistics=[defective_stat]) + with pytest.raises( + RuntimeError, + match="The statistic .* has no sacc_indices", + ): + likelihood.read(sacc_data_for_trivial_stat) + + def test_get_cov_works_after_read(trivial_stats, sacc_data_for_trivial_stat): likelihood = ConstGaussian(statistics=trivial_stats) likelihood.read(sacc_data_for_trivial_stat) assert np.all(likelihood.get_cov() == np.diag([4.0, 9.0, 16.0])) +def test_get_cov_with_one_statistic(trivial_stats, sacc_data_for_trivial_stat): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + + assert np.all(likelihood.get_cov(trivial_stats[0]) == np.diag([4.0, 9.0, 16.0])) + + +def test_get_cov_with_list_of_statistics(trivial_stats, sacc_data_for_trivial_stat): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + + assert isinstance(trivial_stats, list) + cov = likelihood.get_cov(trivial_stats) + assert np.all(cov == np.diag([4.0, 9.0, 16.0])) + + +def test_get_data_vector_fails_before_read(trivial_stats): + likelihood = ConstGaussian(statistics=trivial_stats) + with pytest.raises( + AssertionError, + match=re.escape("read() must be called before get_data_vector()"), + ): + _ = likelihood.get_data_vector() + + +def test_get_data_vector_works_after_read(trivial_stats, sacc_data_for_trivial_stat): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + assert np.all(likelihood.get_data_vector() == np.array([1.0, 4.0, -3.0])) + + +def test_compute_theory_vector_fails_before_read(trivial_stats): + likelihood = ConstGaussian(statistics=trivial_stats) + with pytest.raises( + AssertionError, + match=re.escape("update() must be called before compute_theory_vector()"), + ): + _ = likelihood.compute_theory_vector(ModelingTools()) + + +def test_compute_theory_vector_fails_before_update( + trivial_stats, sacc_data_for_trivial_stat +): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + with pytest.raises( + AssertionError, + match=re.escape("update() must be called before compute_theory_vector()"), + ): + _ = likelihood.compute_theory_vector(ModelingTools()) + + +def test_compute_theory_vector_works_after_read_and_update( + trivial_stats, sacc_data_for_trivial_stat +): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + likelihood.update(firecrown.parameters.ParamsMap(mean=10.5)) + assert np.all( + likelihood.compute_theory_vector(ModelingTools()) + == np.array([10.5, 10.5, 10.5]) + ) + + +def test_get_theory_vector_fails_before_read(trivial_stats): + likelihood = ConstGaussian(statistics=trivial_stats) + with pytest.raises( + AssertionError, + match=re.escape( + "compute_theory_vector() must be called before get_theory_vector()" + ), + ): + _ = likelihood.get_theory_vector() + + +def test_get_theory_vector_fails_before_update( + trivial_stats, sacc_data_for_trivial_stat +): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + with pytest.raises( + AssertionError, + match=re.escape( + "compute_theory_vector() must be called before get_theory_vector()" + ), + ): + _ = likelihood.get_theory_vector() + + +def test_get_theory_vector_fails_before_compute_theory_vector( + trivial_stats, sacc_data_for_trivial_stat +): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + likelihood.update(firecrown.parameters.ParamsMap(mean=10.5)) + with pytest.raises( + AssertionError, + match=re.escape( + "compute_theory_vector() must be called before get_theory_vector()", + ), + ): + _ = likelihood.get_theory_vector() + + +def test_get_theory_vector_works_after_read_update_and_compute_theory_vector( + trivial_stats, sacc_data_for_trivial_stat +): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + likelihood.update(firecrown.parameters.ParamsMap(mean=10.5)) + likelihood.compute_theory_vector(ModelingTools()) + assert np.all(likelihood.get_theory_vector() == np.array([10.5, 10.5, 10.5])) + + +def test_get_theory_vector_fails_after_read_update_compute_theory_vector_and_reset( + trivial_stats, sacc_data_for_trivial_stat +): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + likelihood.update(firecrown.parameters.ParamsMap(mean=10.5)) + likelihood.compute_theory_vector(ModelingTools()) + likelihood.reset() + with pytest.raises( + AssertionError, + match=re.escape( + "compute_theory_vector() must be called before get_theory_vector()" + ), + ): + _ = likelihood.get_theory_vector() + + def test_compute_chisq_fails_before_read(trivial_stats): """Note that the error message from the direct call to compute_chisq notes that update() must be called; this can only be called after read().""" likelihood = ConstGaussian(statistics=trivial_stats) with pytest.raises( - AssertionError, match=r"update\(\) must be called before compute_chisq\(\)" + AssertionError, + match=re.escape("update() must be called before compute_chisq()"), ): _ = likelihood.compute_chisq(ModelingTools()) @@ -50,6 +211,37 @@ def test_compute_chisq(trivial_stats, sacc_data_for_trivial_stat, trivial_params assert likelihood.compute_chisq(ModelingTools()) == 2.0 +def test_deprecated_compute(trivial_stats, sacc_data_for_trivial_stat, trivial_params): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + likelihood.update(trivial_params) + with pytest.warns(DeprecationWarning): + data_vector, theory_vector = likelihood.compute(ModelingTools()) + assert np.all(data_vector == np.array([1.0, 4.0, -3.0])) + assert np.all(theory_vector == np.array([1.0, 1.0, 1.0])) + + +def test_chisquared_compute_vector_not_implemented( + trivial_stats, sacc_data_for_trivial_stat, trivial_params +): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + likelihood.update(trivial_params) + + def compute_theory_vector(_tools: ModelingTools) -> npt.NDArray[np.float64]: + raise NotImplementedError() + + def compute( + _tools: ModelingTools, + ) -> Tuple[npt.NDArray[np.float64], npt.NDArray[np.float64]]: + return np.array([1.0, 4.0, -3.0]), np.array([1.0, 1.0, 1.0]) + + likelihood.compute_theory_vector = compute_theory_vector # type: ignore + likelihood.compute = compute # type: ignore + + assert likelihood.compute_chisq(ModelingTools()) == 2.0 + + def test_required_parameters(trivial_stats, sacc_data_for_trivial_stat, trivial_params): likelihood = ConstGaussian(statistics=trivial_stats) likelihood.read(sacc_data_for_trivial_stat) @@ -87,14 +279,6 @@ def test_missing_covariance(trivial_stats, sacc_with_data_points: sacc.Sacc): likelihood.read(sacc_with_data_points) -def test_get_data_vector_fails_before_read(trivial_stats): - likelihood = ConstGaussian(statistics=trivial_stats) - with pytest.raises( - AssertionError, match=r"read\(\) must be called before get_data_vector\(\)" - ): - _ = likelihood.get_data_vector() - - def test_using_good_sacc( trivial_stats, sacc_data_for_trivial_stat: sacc.Sacc, @@ -117,3 +301,123 @@ def test_after_read_all_statistics_are_ready( for gs in likelihood.statistics: stat: Statistic = gs.statistic assert stat.ready + + +def test_make_realization_chisq( + trivial_stats, + sacc_data_for_trivial_stat: sacc.Sacc, + tools_with_vanilla_cosmology: ModelingTools, +): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + params = firecrown.parameters.ParamsMap(mean=10.5) + likelihood.update(params) + likelihood.compute_chisq(tools_with_vanilla_cosmology) + + new_sacc = likelihood.make_realization(sacc_data_for_trivial_stat) + + new_likelihood = ConstGaussian(statistics=[TrivialStatistic()]) + new_likelihood.read(new_sacc) + params = firecrown.parameters.ParamsMap(mean=10.5) + new_likelihood.update(params) + chisq = new_likelihood.compute_chisq(tools_with_vanilla_cosmology) + + # The new likelihood chisq is distributed as a chi-squared with 3 degrees of + # freedom. We want to check that the new chisq is within the 1-10^-6 quantile + # of the chi-squared distribution. This is equivalent to checking that the + # new chisq is less than the 1-10^-6 quantile. This is expected to fail + # 1 in 10^6 times. + assert chisq < chi2.ppf(1.0 - 1.0e-6, df=3) + + +def test_make_realization_chisq_mean( + trivial_stats, + sacc_data_for_trivial_stat: sacc.Sacc, + tools_with_vanilla_cosmology: ModelingTools, +): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + params = firecrown.parameters.ParamsMap(mean=10.5) + likelihood.update(params) + likelihood.compute_chisq(tools_with_vanilla_cosmology) + + chisq_list = [] + for _ in range(1000): + new_sacc = likelihood.make_realization(sacc_data_for_trivial_stat) + + new_likelihood = ConstGaussian(statistics=[TrivialStatistic()]) + new_likelihood.read(new_sacc) + params = firecrown.parameters.ParamsMap(mean=10.5) + new_likelihood.update(params) + chisq = new_likelihood.compute_chisq(tools_with_vanilla_cosmology) + chisq_list.append(chisq) + + # The new likelihood chisq is distributed as a chi-squared with 3 degrees of + # freedom, so the mean is 3.0 and the variance is 6.0. Since we are computing + # the mean of 1000 realizations, the variance of the mean is 6.0 / 1000.0. + # We want to check that the new chisq is within 5 sigma of the mean. + assert_allclose(np.mean(chisq_list), 3.0, atol=5.0 * np.sqrt(6.0 / 1000.0)) + + +def test_make_realization_data_vector( + trivial_stats, + sacc_data_for_trivial_stat: sacc.Sacc, + tools_with_vanilla_cosmology: ModelingTools, +): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + params = firecrown.parameters.ParamsMap(mean=10.5) + likelihood.update(params) + likelihood.compute_chisq(tools_with_vanilla_cosmology) + + data_vector_list = [] + for _ in range(1000): + new_sacc = likelihood.make_realization(sacc_data_for_trivial_stat) + + new_likelihood = ConstGaussian(statistics=[TrivialStatistic()]) + new_likelihood.read(new_sacc) + params = firecrown.parameters.ParamsMap(mean=10.5) + new_likelihood.update(params) + data_vector = new_likelihood.get_data_vector() + data_vector_list.append(data_vector) + + # The new likelihood data vector is distributed as a Gaussian with mean + # equal to the theory vector and covariance equal to the covariance matrix. + # We want to check that the new data vector is within 5 sigma of the mean. + var_exact = np.array([4.0, 9.0, 16.0]) + assert_allclose( + (np.mean(data_vector_list, axis=0) - np.array([10.5, 10.5, 10.5])) + / np.sqrt(var_exact), + 0.0, + atol=5.0 / np.sqrt(1000.0), + ) + + # The covariance can be computed as the covariance of the data vectors + # minus the covariance of the theory vectors. + covariance = np.cov(np.array(data_vector_list).T) + assert_allclose( + (covariance.diagonal() - var_exact) / np.sqrt(2.0 * var_exact**2 / 999.0), + 1.0, + atol=5, + ) + + +def test_make_realization_no_noise( + trivial_stats, + sacc_data_for_trivial_stat: sacc.Sacc, + tools_with_vanilla_cosmology: ModelingTools, +): + likelihood = ConstGaussian(statistics=trivial_stats) + likelihood.read(sacc_data_for_trivial_stat) + params = firecrown.parameters.ParamsMap(mean=10.5) + likelihood.update(params) + likelihood.compute_chisq(tools_with_vanilla_cosmology) + + new_sacc = likelihood.make_realization(sacc_data_for_trivial_stat, add_noise=False) + + new_likelihood = ConstGaussian(statistics=[TrivialStatistic()]) + new_likelihood.read(new_sacc) + params = firecrown.parameters.ParamsMap(mean=10.5) + new_likelihood.update(params) + + assert_allclose(new_likelihood.get_data_vector(), likelihood.get_theory_vector()) diff --git a/tests/test_modeling_tools.py b/tests/test_modeling_tools.py index 036585f89..4ec3959e6 100644 --- a/tests/test_modeling_tools.py +++ b/tests/test_modeling_tools.py @@ -51,7 +51,7 @@ def test_no_adding_pk_twice(dummy_powerspectrum: pyccl.Pk2D): def test_modeling_tool_prepare_without_update(): tools = ModelingTools() cosmo = pyccl.CosmologyVanillaLCDM() - with pytest.raises(RuntimeError, match="ModelingTools has not been prepared"): + with pytest.raises(RuntimeError, match="ModelingTools has not been updated."): tools.prepare(cosmo) diff --git a/tests/test_pt_systematics.py b/tests/test_pt_systematics.py index 22ae01cb2..31dc7d7e0 100644 --- a/tests/test_pt_systematics.py +++ b/tests/test_pt_systematics.py @@ -143,6 +143,7 @@ def test_pt_systematics(weak_lensing_source, number_counts_source, sacc_data): # TODO: We need to have a way to test systematics without requiring # digging into the innards of a likelihood object. s0 = likelihood.statistics[0].statistic + assert isinstance(s0, TwoPoint) ells = s0.ells cells_GG = s0.cells[("shear", "shear")] cells_GI = s0.cells[("shear", "intrinsic_pt")] @@ -151,12 +152,14 @@ def test_pt_systematics(weak_lensing_source, number_counts_source, sacc_data): # print(list(likelihood.statistics[2].cells.keys())) s2 = likelihood.statistics[2].statistic + assert isinstance(s2, TwoPoint) cells_gG = s2.cells[("galaxies", "shear")] cells_gI = s2.cells[("galaxies", "intrinsic_pt")] cells_mI = s2.cells[("magnification+rsd", "intrinsic_pt")] # print(list(likelihood.statistics[3].cells.keys())) s3 = likelihood.statistics[3].statistic + assert isinstance(s3, TwoPoint) cells_gg = s3.cells[("galaxies", "galaxies")] cells_gm = s3.cells[("galaxies", "magnification+rsd")] cells_mm = s3.cells[("magnification+rsd", "magnification+rsd")] @@ -303,6 +306,7 @@ def test_pt_mixed_systematics(sacc_data): # pylint: disable=no-member s0 = likelihood.statistics[0].statistic + assert isinstance(s0, TwoPoint) ells = s0.ells # print(list(likelihood.statistics[2].cells.keys())) diff --git a/tests/test_updatable.py b/tests/test_updatable.py index a28da1ea3..2efcca0aa 100644 --- a/tests/test_updatable.py +++ b/tests/test_updatable.py @@ -60,15 +60,17 @@ def test_simple_updatable(): assert obj.required_parameters() == expected_params assert obj.x is None assert obj.y is None + assert not obj.is_updated() new_params = ParamsMap({"x": -1.0, "y": 5.5}) obj.update(new_params) assert obj.x == -1.0 assert obj.y == 5.5 + assert obj.is_updated() # pylint: disable-msg=E1101 def test_updatable_collection_appends(): - coll = UpdatableCollection() + coll: UpdatableCollection = UpdatableCollection() assert len(coll) == 0 coll.append(SimpleUpdatable()) @@ -84,7 +86,7 @@ def test_updatable_collection_appends(): def test_updatable_collection_updates(): - coll = UpdatableCollection() + coll: UpdatableCollection = UpdatableCollection() assert len(coll) == 0 coll.append(SimpleUpdatable()) @@ -100,11 +102,11 @@ def test_updatable_collection_updates(): def test_updatable_collection_rejects_nonupdatables(): - coll = UpdatableCollection() + coll: UpdatableCollection = UpdatableCollection() assert len(coll) == 0 with pytest.raises(TypeError): - coll.append(3) # type: ignore + coll.append(3) assert len(coll) == 0 @@ -231,6 +233,21 @@ def test_update_rejects_internal_parameters(): assert my_updatable.the_meaning_of_life == 42.0 +def test_updatable_collection_is_updated(): + obj: UpdatableCollection = UpdatableCollection([SimpleUpdatable()]) + new_params = {"x": -1.0, "y": 5.5} + + assert not obj.is_updated() + obj.update(ParamsMap(new_params)) + assert obj.is_updated() + + +def test_updatablecollection_without_derived_parameters(): + obj: UpdatableCollection = UpdatableCollection() + + assert obj.get_derived_parameters() is None + + @pytest.fixture(name="nested_updatables", params=permutations(range(3))) def fixture_nested_updatables(request): updatables = np.array(