More testing of CosmoSIS module (#313)

* Bugfix for and more testing of CosmoSIS module
* Remove extraneous import
* FIx oddly-formatted concatenated strings
* Specify "utf-8" encoding for all open calls
* Test CosmoSIS module with derived parameters
* Do not insert derived_param0 in data from sampler
* Fix documentation and comments
* Add len support for RequiredParameters
* Improve error message on CosmoSIS misconfiguration
* Improve error report due to a missing parameter
* If sacc_tracer is not found, raise an error
* Make get_data_vector, compute_theory_vector abstract
* Introduce StatisticUnreadError
* Unpack tuple at call site
* Add GuardedStatistic
* Move TrivialStatistic
* More error checking in GaussFamily.read
* Update use of mypy
* Relax pylint rules, and apply everywhere
* Fix pylint complaints in srd_sn examples
* Fix example so that it runs
* Refactoring guided by pylint
* Attempt to diagnose some statistic read failures
* Fix malformed SACC object (add tag to data points)

.github/workflows/ci.yml

@@ -86,19 +86,15 @@ jobs:
       - name: Running mypy
         shell: bash -l {0}
         run: |
-          mypy -p firecrown --ignore-missing-imports
-          mypy -p examples --ignore-missing-imports
-          mypy -p tests --ignore-missing-imports
+          mypy -p firecrown
+          mypy -p examples
+          mypy -p tests
       - name: Running pylint
         shell: bash -l {0}
         run: |
           pylint --rcfile tests/pylintrc  tests
           pylint --rcfile firecrown/models/pylintrc firecrown/models
-          pylint firecrown/connector
-          pylint firecrown/*.py
-          pylint firecrown/likelihood/*.py
-          pylint firecrown/likelihood/gauss_family/*.py
-
+          pylint firecrown
       - name: Running pytest
         shell: bash -l {0}
         run: python -m pytest -vv

docs/contrib.rst

@@ -26,7 +26,7 @@ Please run:
 
 .. code:: bash
 
-    mypy firecrown/ --ignore-missing-imports 
+    mypy -p firecrown -p examples -p tests
 
 and fix any errors reported before pushing commits to the GitHub repository.
 

examples/des_y1_3x2pt/des_y1_3x2pt_PT.py

@@ -1,11 +1,12 @@
 #!/usr/bin/env python
 
 """Example factory function for DES Y1 3x2pt likelihood."""
-
+from dataclasses import dataclass
 import os
 
 from typing import Dict, Union, Tuple
 
+import numpy as np
 import sacc
 import pyccl as ccl
 import pyccl.nl_pt
@@ -26,6 +27,48 @@
 )
 
 
+@dataclass
+class CclSetup:
+    """A package of related CCL parameters, to reduce the number of variables
+    used in the :python:`run_likelihood` method."""
+
+    a_1: float = 1.0
+    a_2: float = 0.5
+    a_d: float = 0.5
+    b_1: float = 2.0
+    b_2: float = 1.0
+    b_s: float = 1.0
+    mag_bias: float = 1.0
+
+
+@dataclass
+class CElls:
+    GG: np.ndarray
+    GI: np.ndarray
+    II: np.ndarray
+    cs_total: np.ndarray
+    gG: np.ndarray
+    gI: np.ndarray
+    mI: np.ndarray
+    gg: np.ndarray
+    gm: np.ndarray
+    gg_total: np.ndarray
+
+    def __init__(self, stat0: TwoPoint, stat2: TwoPoint, stat3: TwoPoint):
+        self.GG = stat0.cells[("shear", "shear")]
+        self.GI = stat0.cells[("shear", "intrinsic_pt")]
+        self.II = stat0.cells[("intrinsic_pt", "intrinsic_pt")]
+        self.cs_total = stat0.cells["total"]
+
+        self.gG = stat2.cells[("galaxies", "shear")]
+        self.gI = stat2.cells[("galaxies", "intrinsic_pt")]
+        self.mI = stat2.cells[("magnification+rsd", "intrinsic_pt")]
+
+        self.gg = stat3.cells[("galaxies", "galaxies")]
+        self.gm = stat3.cells[("galaxies", "magnification+rsd")]
+        self.gg_total = stat3.cells["total"]
+
+
 def build_likelihood(_) -> Tuple[Likelihood, ModelingTools]:
     """Likelihood factory function for DES Y1 3x2pt analysis."""
 
@@ -115,11 +158,6 @@ def run_likelihood() -> None:
     """Produce some plots using the likelihood function built by
     :python:`build_likelihood`.
     """
-    # We do imports here to save a bit of time when importing this module but
-    # not using the run_likelihood function.
-    # pylint: disable=import-outside-toplevel
-    import numpy as np
-    import matplotlib.pyplot as plt
 
     # pylint: enable=import-outside-toplevel
 
@@ -137,19 +175,9 @@ def run_likelihood() -> None:
     ccl_cosmo = ccl.CosmologyVanillaLCDM()
     ccl_cosmo.compute_nonlin_power()
 
-    # Bare CCL setup
-    a_1 = 1.0
-    a_2 = 0.5
-    a_d = 0.5
-
-    b_1 = 2.0
-    b_2 = 1.0
-    b_s = 1.0
-
-    mag_bias = 1.0
-
+    cs = CclSetup()
     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
+        ccl_cosmo, z=z, a1=cs.a_1, a1delta=cs.a_d, a2=cs.a_2, Om_m2_for_c2=False
     )
 
     # Code that creates a Pk2D object:
@@ -165,27 +193,27 @@ def run_likelihood() -> None:
         c1=(z, c_1), c2=(z, c_2), cdelta=(z, c_d)
     )
     ptt_m = pyccl.nl_pt.PTMatterTracer()
-    ptt_g = pyccl.nl_pt.PTNumberCountsTracer(b1=b_1, b2=b_2, bs=b_s)
+    ptt_g = pyccl.nl_pt.PTNumberCountsTracer(b1=cs.b_1, b2=cs.b_2, bs=cs.b_s)
     # IA
-    pk_im = ptc.get_biased_pk2d(ccl_cosmo, ptt_i, tracer2=ptt_m)
-    pk_ii = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_i, tracer2=ptt_i)
-    pk_gi = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_g, tracer2=ptt_i)
+    pk_im = ptc.get_biased_pk2d(ptt_i, tracer2=ptt_m)
+    pk_ii = ptc.get_biased_pk2d(ptt_i, tracer2=ptt_i)
+    pk_gi = ptc.get_biased_pk2d(ptt_g, tracer2=ptt_i)
     # Galaxies
-    pk_gm = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_g, tracer2=ptt_m)
-    pk_gg = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_g, tracer2=ptt_g)
+    pk_gm = ptc.get_biased_pk2d(ptt_g, tracer2=ptt_m)
+    pk_gg = ptc.get_biased_pk2d(ptt_g, tracer2=ptt_g)
     # Magnification: just a matter-matter P(k)
-    pk_mm = ptc.get_biased_pk2d(ccl_cosmo, tracer1=ptt_m, tracer2=ptt_m)
+    pk_mm = ptc.get_biased_pk2d(ptt_m, tracer2=ptt_m)
 
     # Set the parameters for our systematics
     systematics_params = ParamsMap(
         {
-            "ia_a_1": a_1,
-            "ia_a_2": a_2,
-            "ia_a_d": a_d,
-            "lens0_bias": b_1,
-            "lens0_b_2": b_2,
-            "lens0_b_s": b_s,
-            "lens0_mag_bias": mag_bias,
+            "ia_a_1": cs.a_1,
+            "ia_a_2": cs.a_2,
+            "ia_a_d": cs.a_d,
+            "lens0_bias": cs.b_1,
+            "lens0_b_2": cs.b_2,
+            "lens0_b_s": cs.b_s,
+            "lens0_mag_bias": cs.mag_bias,
             "src0_delta_z": 0.000,
             "lens0_delta_z": 0.000,
         }
@@ -202,33 +230,70 @@ def run_likelihood() -> None:
 
     print(f"Log-like = {log_like:.1f}")
 
-    # Plot the predicted and measured statistic
-    # x = likelihood.statistics[0].ell_or_theta_
-    # y_data = likelihood.statistics[0].measured_statistic_
-
     assert isinstance(likelihood, ConstGaussian)
     assert likelihood.cov is not None
 
+    stat0 = likelihood.statistics[0].statistic
+
+    # x = likelihood.statistics[0].ell_or_theta_
+    # y_data = likelihood.statistics[0].measured_statistic_
+
     # y_err = np.sqrt(np.diag(likelihood.cov))[: len(x)]
     # y_theory = likelihood.statistics[0].predicted_statistic_
 
-    print(list(likelihood.statistics[0].cells.keys()))
+    print(list(stat0.cells.keys()))
 
-    ells = likelihood.statistics[0].ells
-    cells_GG = likelihood.statistics[0].cells[("shear", "shear")]
-    cells_GI = likelihood.statistics[0].cells[("shear", "intrinsic_pt")]
-    cells_II = likelihood.statistics[0].cells[("intrinsic_pt", "intrinsic_pt")]
-    cells_cs_total = likelihood.statistics[0].cells["total"]
+    stat2 = likelihood.statistics[2].statistic
+    assert isinstance(stat2, TwoPoint)
+    print(list(stat2.cells.keys()))
 
-    print(list(likelihood.statistics[2].cells.keys()))
-    cells_gG = likelihood.statistics[2].cells[("galaxies", "shear")]
-    cells_gI = likelihood.statistics[2].cells[("galaxies", "intrinsic_pt")]
-    cells_mI = likelihood.statistics[2].cells[("magnification+rsd", "intrinsic_pt")]
+    stat3 = likelihood.statistics[3].statistic
+    print(list(stat3.cells.keys()))
 
-    print(list(likelihood.statistics[3].cells.keys()))
-    cells_gg = likelihood.statistics[3].cells[("galaxies", "galaxies")]
-    cells_gm = likelihood.statistics[3].cells[("galaxies", "magnification+rsd")]
-    cells_gg_total = likelihood.statistics[3].cells["total"]
+    plot_predicted_and_measured_statistics(
+        ccl_cosmo,
+        cs,
+        lens_nz,
+        lens_z,
+        nz,
+        pk_gg,
+        pk_gi,
+        pk_gm,
+        pk_ii,
+        pk_im,
+        pk_mm,
+        stat0,
+        stat2,
+        stat3,
+        z,
+    )
+
+
+def plot_predicted_and_measured_statistics(
+    ccl_cosmo,
+    cs,
+    lens_nz,
+    lens_z,
+    nz,
+    pk_gg,
+    pk_gi,
+    pk_gm,
+    pk_ii,
+    pk_im,
+    pk_mm,
+    stat0,
+    stat2,
+    stat3,
+    z,
+):
+    """Plot the predictions and measurements."""
+    # We do imports here to save a bit of time when importing this module but
+    # not using the run_likelihood function.
+    # pylint: disable=import-outside-toplevel
+    import matplotlib.pyplot as plt
+
+    ells = stat0.ells
+    cells = CElls(stat0, stat2, stat3)
 
     # Code that computes effect from IA using that Pk2D object
     t_lens = ccl.WeakLensingTracer(ccl_cosmo, dndz=(z, nz))
@@ -250,13 +315,12 @@ def run_likelihood() -> None:
         has_rsd=True,
         dndz=(lens_z, lens_nz),
         bias=None,
-        mag_bias=(lens_z, mag_bias * np.ones_like(lens_z)),
+        mag_bias=(lens_z, cs.mag_bias * np.ones_like(lens_z)),
     )
     cl_GI = ccl.angular_cl(ccl_cosmo, t_lens, t_ia, ells, p_of_k_a=pk_im)
     cl_II = ccl.angular_cl(ccl_cosmo, t_ia, t_ia, ells, p_of_k_a=pk_ii)
     # The weak gravitational lensing power spectrum
     cl_GG = ccl.angular_cl(ccl_cosmo, t_lens, t_lens, ells)
-
     # Galaxies
     cl_gG = ccl.angular_cl(ccl_cosmo, t_g, t_lens, ells, p_of_k_a=pk_gm)
     cl_gI = ccl.angular_cl(ccl_cosmo, t_g, t_ia, ells, p_of_k_a=pk_gi)
@@ -265,36 +329,32 @@ def run_likelihood() -> None:
     cl_mI = ccl.angular_cl(ccl_cosmo, t_m, t_ia, ells, p_of_k_a=pk_im)
     cl_gm = ccl.angular_cl(ccl_cosmo, t_g, t_m, ells, p_of_k_a=pk_gm)
     cl_mm = ccl.angular_cl(ccl_cosmo, t_m, t_m, ells, p_of_k_a=pk_mm)
-
     # The observed angular power spectrum is the sum of the two.
     cl_cs_theory = cl_GG + 2 * cl_GI + cl_II
     cl_gg_theory = cl_gg + 2 * cl_gm + cl_mm
-
     fig, ax = plt.subplots(2, 1, sharex=True, figsize=(6, 6))
     fig.subplots_adjust(hspace=0)
     # ax[0].plot(x, y_theory, label="Total")
-    ax[0].plot(ells, cells_GG, label="GG firecrown")
+    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")
+    ax[0].plot(ells, -cells.GI, label="-GI firecrown")
     ax[0].plot(ells, -cl_GI, ls="--", label="-GI CCL")
-    ax[0].plot(ells, cells_II, label="II firecrown")
+    ax[0].plot(ells, cells.II, label="II firecrown")
     ax[0].plot(ells, cl_II, ls="--", label="II CCL")
-    ax[0].plot(ells, -cells_gI, label="-Ig firecrown")
+    ax[0].plot(ells, -cells.gI, label="-Ig firecrown")
     ax[0].plot(ells, -cl_gI, ls="--", label="-Ig CCL")
-    ax[0].plot(ells, cells_cs_total, label="total CS firecrown")
+    ax[0].plot(ells, cells.cs_total, label="total CS firecrown")
     ax[0].plot(ells, cl_cs_theory, ls="--", label="total CS CCL")
-
-    ax[1].plot(ells, cells_gG, label="Gg firecrown")
+    ax[1].plot(ells, cells.gG, label="Gg firecrown")
     ax[1].plot(ells, cl_gG, ls="--", label="Gg CCL")
-    ax[1].plot(ells, cells_gg, label="gg firecrown")
+    ax[1].plot(ells, cells.gg, label="gg firecrown")
     ax[1].plot(ells, cl_gg, ls="--", label="gg CCL")
-    ax[1].plot(ells, -cells_mI, label="-mI firecrown")
+    ax[1].plot(ells, -cells.mI, label="-mI firecrown")
     ax[1].plot(ells, -cl_mI, ls="--", label="-mI CCL")
-    ax[1].plot(ells, cells_gm, label="gm firecrown")
+    ax[1].plot(ells, cells.gm, label="gm firecrown")
     ax[1].plot(ells, cl_gm, ls="--", label="gm CCL")
-    ax[1].plot(ells, cells_gg_total, label="total gg firecrown")
+    ax[1].plot(ells, cells.gg_total, label="total gg firecrown")
     ax[1].plot(ells, cl_gg_theory, ls="--", label="total gg CCL")
-
     # ax[0].errorbar(x, y_data, y_err, ls="none", marker="o")
     ax[0].set_xscale("log")
     ax[1].set_xlabel(r"$\ell$")
@@ -303,10 +363,8 @@ def run_likelihood() -> None:
         a.set_yscale("log")
         a.set_ylabel(r"$C_\ell$")
         a.legend(fontsize="small")
-
     fig.suptitle("PT Cls, including IA, galaxy bias, magnification")
     fig.savefig("pt_cls.png", facecolor="white", dpi=300)
-
     plt.show()