Merge branch 'master' into cluster_example_fixingSaccDemo_2024retreat

LSSTDESC · Sep 24, 2024 · 6916bef · 6916bef
2 parents 6598365 + 097769c
commit 6916bef
Show file tree

Hide file tree

Showing 79 changed files with 6,800 additions and 2,099 deletions.
diff --git a/examples/des_y1_3x2pt/des_y1_3x2pt.py b/examples/des_y1_3x2pt/des_y1_3x2pt.py
@@ -7,6 +7,8 @@
 import firecrown.likelihood.number_counts as nc
 from firecrown.likelihood.two_point import TwoPoint
 from firecrown.likelihood.gaussian import ConstGaussian
+from firecrown.modeling_tools import ModelingTools
+from firecrown.ccl_factory import CCLFactory
 
 
 # The likelihood used for DES Y1 3x2pt analysis is a Gaussian likelihood, which
@@ -111,6 +113,8 @@ def build_likelihood(_):
     # file and the sources their respective dndz.
     likelihood.read(sacc_data)
 
+    modeling_tools = ModelingTools(ccl_factory=CCLFactory(require_nonlinear_pk=True))
+
     # This script will be loaded by the appropriated connector. The framework
     # will call the factory function that should return a Likelihood instance.
-    return likelihood
+    return likelihood, modeling_tools
diff --git a/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py b/examples/des_y1_3x2pt/des_y1_3x2pt_PT.py
@@ -20,6 +20,8 @@
 from firecrown.parameters import ParamsMap
 from firecrown.modeling_tools import ModelingTools
 from firecrown.likelihood.likelihood import Likelihood
+from firecrown.ccl_factory import CCLFactory
+from firecrown.updatable import get_default_params_map
 
 
 saccfile = os.path.expanduser(
@@ -137,7 +139,9 @@ def build_likelihood(_) -> tuple[Likelihood, ModelingTools]:
         nk_per_decade=20,
     )
 
-    modeling_tools = ModelingTools(pt_calculator=pt_calculator)
+    modeling_tools = ModelingTools(
+        pt_calculator=pt_calculator, ccl_factory=CCLFactory(require_nonlinear_pk=True)
+    )
     # Note that the ordering of the statistics is relevant, because the data
     # vector, theory vector, and covariance matrix will be organized to
     # follow the order used here.
@@ -173,9 +177,11 @@ def run_likelihood() -> None:
     z, nz = src0_tracer.z, src0_tracer.nz
     lens_z, lens_nz = lens0_tracer.z, lens0_tracer.nz
 
-    # Define a ccl.Cosmology object using default parameters
-    ccl_cosmo = ccl.CosmologyVanillaLCDM()
-    ccl_cosmo.compute_nonlin_power()
+    # Setting cosmological parameters
+    cosmo_params = get_default_params_map(tools)
+    tools.update(cosmo_params)
+    tools.prepare()
+    ccl_cosmo = tools.get_ccl_cosmology()
 
     cs = CclSetup()
     c_1, c_d, c_2 = pyccl.nl_pt.translate_IA_norm(
@@ -224,9 +230,6 @@ def run_likelihood() -> None:
     # Apply the systematics parameters
     likelihood.update(systematics_params)
 
-    # Prepare the cosmology object
-    tools.prepare(ccl_cosmo)
-
     # Compute the log-likelihood, using the ccl.Cosmology object as the input
     log_like = likelihood.compute_loglike(tools)
 
@@ -296,7 +299,7 @@ def plot_predicted_and_measured_statistics(
     # pylint: disable=import-outside-toplevel
     import matplotlib.pyplot as plt
 
-    ells = stat0.ells
+    ells = stat0.ells_for_xi
     cells = CElls(stat0, stat2, stat3)
 
     # Code that computes effect from IA using that Pk2D object
@@ -339,6 +342,9 @@ def plot_predicted_and_measured_statistics(
     fig, ax = plt.subplots(2, 1, sharex=True, figsize=(6, 6))
     fig.subplots_adjust(hspace=0)
     # ax[0].plot(x, y_theory, label="Total")
+    assert isinstance(ax, np.ndarray)
+    assert isinstance(ax[0], plt.Axes)
+    assert isinstance(ax[1], plt.Axes)
     ax[0].plot(ells, cells.GG, label="GG firecrown")
     ax[0].plot(ells, cl_GG, ls="--", label="GG CCL")
     ax[0].plot(ells, -cells.GI, label="-GI firecrown")

diff --git a/examples/des_y1_3x2pt/des_y1_3x2pt_default_factory.ini b/examples/des_y1_3x2pt/des_y1_3x2pt_default_factory.ini
@@ -0,0 +1,60 @@
+[runtime]
+sampler = test
+root = ${PWD}
+
+[DEFAULT]
+fatal_errors = T
+
+[output]
+filename = output/des_y1_3x2pt_samples.txt
+format = text
+verbosity = 0
+
+[pipeline]
+modules = consistency camb firecrown_likelihood
+values = ${FIRECROWN_DIR}/examples/des_y1_3x2pt/des_y1_3x2pt_values.ini
+likelihoods = firecrown
+quiet = T
+debug = T
+timing = T
+extra_output = TwoPoint/NumberCountsScale_lens0 TwoPoint/NumberCountsScale_lens1 TwoPoint/NumberCountsScale_lens2 TwoPoint/NumberCountsScale_lens3 TwoPoint/NumberCountsScale_lens4
+
+[consistency]
+file = ${CSL_DIR}/utility/consistency/consistency_interface.py
+
+[camb]
+file = ${CSL_DIR}/boltzmann/camb/camb_interface.py
+
+mode = all
+lmax = 2500
+feedback = 0
+zmin = 0.0
+zmax = 4.0
+nz = 100
+kmin = 1e-4
+kmax = 50.0
+nk = 1000
+
+[firecrown_likelihood]
+;; Fix this to use an environment variable to find the files.
+;; Set FIRECROWN_DIR to the base of the firecrown installation (or build, if you haven't
+;; installed it)
+file = ${FIRECROWN_DIR}/firecrown/connector/cosmosis/likelihood.py
+likelihood_source = firecrown.likelihood.factories.build_two_point_likelihood
+likelihood_config = ${FIRECROWN_DIR}/examples/des_y1_3x2pt/des_y1_3x2pt_experiment.yaml
+;; Connector settings
+require_nonlinear_pk = True
+sampling_parameters_sections = firecrown_two_point
+
+[test]
+fatal_errors = T
+save_dir = des_y1_3x2pt_output
+
+[metropolis]
+samples = 1000
+nsteps = 1
+
+[emcee]
+walkers = 64
+samples = 400
+nsteps = 10
diff --git a/examples/des_y1_3x2pt/des_y1_3x2pt_experiment.yaml b/examples/des_y1_3x2pt/des_y1_3x2pt_experiment.yaml
@@ -0,0 +1,29 @@
+---
+
+# This file contains all information about the DES Y1 3x2pt experiment.
+
+# The 'data_source:' field points to the SACC file containing the DES Y1 3x2pt data vector 
+# and covariance matrix, which is used for analysis.
+
+# The 'two_point_factory:' field points to the factory that will create the TwoPoint 
+# objects. These objects represent the chosen theoretical models for each data point 
+# found in the SACC file.
+
+data_source:
+  sacc_data_file: des_y1_3x2pt_sacc_data.fits
+
+# The two point statistics are defined by the TwoPoint objects. The TwoPoint statistics
+# are created using the factories defined in this file.
+two_point_factory:
+  correlation_space: real
+  number_counts_factory:
+    global_systematics: []
+    per_bin_systematics:
+      - type: PhotoZShiftFactory
+  weak_lensing_factory:
+    global_systematics:
+      - alphag: 1
+        type: LinearAlignmentSystematicFactory
+    per_bin_systematics:
+      - type: MultiplicativeShearBiasFactory
+      - type: PhotoZShiftFactory
diff --git a/examples/des_y1_3x2pt/des_y1_cosmic_shear_TATT.py b/examples/des_y1_3x2pt/des_y1_cosmic_shear_TATT.py
@@ -11,7 +11,9 @@
 from firecrown.parameters import ParamsMap
 from firecrown.modeling_tools import ModelingTools
 from firecrown.likelihood.likelihood import Likelihood
-
+from firecrown.ccl_factory import CCLFactory
+from firecrown.updatable import get_default_params_map
+from firecrown.metadata_types import TracerNames, TRACER_NAMES_TOTAL
 
 SACCFILE = os.path.expanduser(
     os.path.expandvars(
@@ -38,7 +40,9 @@ def build_likelihood(_) -> tuple[Likelihood, ModelingTools]:
         nk_per_decade=20,
     )
 
-    modeling_tools = ModelingTools(pt_calculator=pt_calculator)
+    modeling_tools = ModelingTools(
+        pt_calculator=pt_calculator, ccl_factory=CCLFactory(require_nonlinear_pk=True)
+    )
     likelihood = ConstGaussian(statistics=list(stats.values()))
 
     # Read the two-point data from the sacc file
@@ -102,14 +106,33 @@ def run_likelihood() -> None:
     src0_tracer = sacc_data.get_tracer("src0")
     z, nz = src0_tracer.z, src0_tracer.nz
 
-    # Define a ccl.Cosmology object using default parameters
-    ccl_cosmo = ccl.CosmologyVanillaLCDM()
-    ccl_cosmo.compute_nonlin_power()
-
     # Bare CCL setup
     a_1 = 1.0
     a_2 = 0.5
     a_d = 0.5
+    # Set the parameters for our systematics
+    systematics_params = ParamsMap(
+        {
+            "ia_a_1": a_1,
+            "ia_a_2": a_2,
+            "ia_a_d": a_d,
+            "src0_delta_z": 0.000,
+            "src1_delta_z": 0.003,
+            "src2_delta_z": -0.001,
+            "src3_delta_z": 0.002,
+        }
+    )
+    # Prepare the cosmology object
+    params = ParamsMap(get_default_params_map(tools) | systematics_params)
+
+    # Apply the systematics parameters
+    likelihood.update(params)
+
+    # Prepare the cosmology object
+    tools.update(params)
+    tools.prepare()
+    ccl_cosmo = tools.get_ccl_cosmology()
+
     c_1, c_d, c_2 = pyccl.nl_pt.translate_IA_norm(
         ccl_cosmo, z=z, a1=a_1, a1delta=a_d, a2=a_2, Om_m2_for_c2=False
     )
@@ -131,25 +154,6 @@ def run_likelihood() -> None:
     pk_im = ptc.get_biased_pk2d(tracer1=ptt_i, tracer2=ptt_m)
     pk_ii = ptc.get_biased_pk2d(tracer1=ptt_i, tracer2=ptt_i)
 
-    # Set the parameters for our systematics
-    systematics_params = ParamsMap(
-        {
-            "ia_a_1": a_1,
-            "ia_a_2": a_2,
-            "ia_a_d": a_d,
-            "src0_delta_z": 0.000,
-            "src1_delta_z": 0.003,
-            "src2_delta_z": -0.001,
-            "src3_delta_z": 0.002,
-        }
-    )
-
-    # Apply the systematics parameters
-    likelihood.update(systematics_params)
-
-    # Prepare the cosmology object
-    tools.prepare(ccl_cosmo)
-
     # Compute the log-likelihood, using the ccl.Cosmology object as the input
     log_like = likelihood.compute_loglike(tools)
 
@@ -180,11 +184,11 @@ def make_plot(ccl_cosmo, nz, pk_ii, pk_im, two_point_0, z):
     import numpy as np  # pylint: disable-msg=import-outside-toplevel
     import matplotlib.pyplot as plt  # pylint: disable-msg=import-outside-toplevel
 
-    ells = two_point_0.ells
-    cells_gg = two_point_0.cells[("shear", "shear")]
-    cells_gi = two_point_0.cells[("shear", "intrinsic_pt")]
-    cells_ii = two_point_0.cells[("intrinsic_pt", "intrinsic_pt")]
-    cells_total = two_point_0.cells["total"]
+    ells = two_point_0.ells_for_xi
+    cells_gg = two_point_0.cells[TracerNames("shear", "shear")]
+    cells_gi = two_point_0.cells[TracerNames("shear", "intrinsic_pt")]
+    cells_ii = two_point_0.cells[TracerNames("intrinsic_pt", "intrinsic_pt")]
+    cells_total = two_point_0.cells[TRACER_NAMES_TOTAL]
     # pylint: enable=no-member
     # Code that computes effect from IA using that Pk2D object
     t_lens = ccl.WeakLensingTracer(ccl_cosmo, dndz=(z, nz))

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
@@ -12,10 +12,12 @@
 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.parameters import ParamsMap, register_new_updatable_parameter
 from firecrown.modeling_tools import ModelingTools, PowerspectrumModifier
 from firecrown.likelihood.likelihood import Likelihood
-
+from firecrown.ccl_factory import CCLFactory
+from firecrown.updatable import get_default_params_map
+from firecrown.metadata_types import TracerNames
 
 SACCFILE = os.path.expanduser(
     os.path.expandvars(
@@ -36,7 +38,7 @@ def __init__(self, pk_to_modify: str = "delta_matter:delta_matter"):
         super().__init__()
         self.pk_to_modify = pk_to_modify
         self.vD19 = pyccl.baryons.BaryonsvanDaalen19()
-        self.f_bar = create()
+        self.f_bar = register_new_updatable_parameter(default_value=0.5)
 
     def compute_p_of_k_z(self, tools: ModelingTools) -> pyccl.Pk2D:
         """Compute the 3D power spectrum P(k, z)."""
@@ -56,7 +58,9 @@ def build_likelihood(_) -> tuple[Likelihood, ModelingTools]:
 
     # Define the power spectrum modification and add it to the ModelingTools
     pk_modifier = vanDaalen19Baryonfication(pk_to_modify="delta_matter:delta_matter")
-    modeling_tools = ModelingTools(pk_modifiers=[pk_modifier])
+    modeling_tools = ModelingTools(
+        pk_modifiers=[pk_modifier], ccl_factory=CCLFactory(require_nonlinear_pk=True)
+    )
 
     # Create the likelihood from the statistics
     likelihood = ConstGaussian(statistics=list(stats.values()))
@@ -123,18 +127,7 @@ def run_likelihood() -> None:
     src0_tracer = sacc_data.get_tracer("src0")
     z, nz = src0_tracer.z, src0_tracer.nz
 
-    # Define a ccl.Cosmology object using default parameters
-    ccl_cosmo = ccl.CosmologyVanillaLCDM()
-    ccl_cosmo.compute_nonlin_power()
-
     f_bar = 0.5
-
-    # Calculate the barynic effects directly with CCL
-    vD19 = pyccl.BaryonsvanDaalen19(fbar=f_bar)
-    pk_baryons = vD19.include_baryonic_effects(
-        cosmo=ccl_cosmo, pk=ccl_cosmo.get_nonlin_power()
-    )
-
     # Set the parameters for our systematics
     systematics_params = ParamsMap(
         {
@@ -145,13 +138,22 @@ def run_likelihood() -> None:
             "src3_delta_z": 0.002,
         }
     )
+    # Prepare the cosmology object
+    params = ParamsMap(get_default_params_map(tools) | systematics_params)
 
-    # Apply the systematics parameters
-    likelihood.update(systematics_params)
+    tools.update(params)
+    tools.prepare()
 
-    # Prepare the cosmology object
-    tools.update(systematics_params)
-    tools.prepare(ccl_cosmo)
+    ccl_cosmo = tools.get_ccl_cosmology()
+
+    # Calculate the barynic effects directly with CCL
+    vD19 = pyccl.BaryonsvanDaalen19(fbar=f_bar)
+    pk_baryons = vD19.include_baryonic_effects(
+        cosmo=ccl_cosmo, pk=ccl_cosmo.get_nonlin_power()
+    )
+
+    # Apply the systematics parameters
+    likelihood.update(params)
 
     # Compute the log-likelihood, using the ccl.Cosmology object as the input
     log_like = likelihood.compute_loglike(tools)
@@ -166,7 +168,7 @@ def run_likelihood() -> None:
 
     # Predict CCL Cl
     wl_tracer = ccl.WeakLensingTracer(ccl_cosmo, dndz=(z, nz))
-    ell = two_point_0.ells
+    ell = two_point_0.ells_for_xi
     cl_dm = ccl.angular_cl(
         cosmo=ccl_cosmo,
         tracer1=wl_tracer,
@@ -191,7 +193,7 @@ def make_plot(ell, cl_dm, cl_baryons, two_point_0):
     """Create and show a diagnostic plot."""
     import matplotlib.pyplot as plt  # pylint: disable-msg=import-outside-toplevel
 
-    cl_firecrown = two_point_0.cells[("shear", "shear")]
+    cl_firecrown = two_point_0.cells[TracerNames("shear", "shear")]
 
     plt.plot(ell, cl_firecrown / cl_dm, label="firecrown w/ baryons")
     plt.plot(ell, cl_baryons / cl_dm, ls="--", label="CCL w/ baryons")