Euclid requisites -- CAMB & CLASS

CHANGELOG.md

@@ -2,12 +2,20 @@
 
 - Documented uses of `Model` class in general contexts (previously only cosmo)
 - `Model` methods to compute log-probabilities and derived parameters now have an `as_dict` keyword (default `False`), for more informative return value.
+- Added ``--minimize`` flag to ``cobaya-run`` for quick minimization (replaces sampler, uses previous output).
 
 ### Cosmological likelihoods and theory codes
 
+- `Pk_interpolator`: added extrapolation up to `extrap_kmin` and improved robustness
+
 #### CAMB
 
 - Removed problematic `zrei: zre` alias (fixes #199, thanks @pcampeti)
+- Added `Omega_b|cdm|nu_massive(z)` and `angular_diameter_distance_2`
+
+#### CLASS
+
+- Updated to v3.1.1
 
 #### BAO
 

cobaya/likelihoods/base_classes/des.py

@@ -567,10 +567,10 @@ def chi_squared(self, theory, return_theory_vector=False):
 
     def logp(self, **params_values):
         PKdelta = self.provider.get_Pk_interpolator(("delta_tot", "delta_tot"),
-                                                    extrap_kmax=500 * self.acc)
+                                                    extrap_kmax=3000 * self.acc)
         if self.use_Weyl:
             PKWeyl = self.provider.get_Pk_interpolator(("Weyl", "Weyl"),
-                                                       extrap_kmax=500 * self.acc)
+                                                       extrap_kmax=3000 * self.acc)
         else:
             PKWeyl = None
 

cobaya/run.py

@@ -36,6 +36,7 @@ def run(info_or_yaml_or_file: Union[InputDict, str, os.PathLike],
         debug: Union[bool, int, None] = None,
         stop_at_error: Optional[bool] = None,
         resume: bool = False, force: bool = False,
+        minimize: Optional[bool] = None,
         no_mpi: bool = False, test: bool = False,
         override: Optional[InputDict] = None,
         ) -> Union[InfoSamplerTuple, PostTuple]:
@@ -77,6 +78,9 @@ def run(info_or_yaml_or_file: Union[InputDict, str, os.PathLike],
             info["resume"] = bool(resume)
             info["force"] = bool(force)
         if info.get("post"):
+            if minimize:
+                raise ValueError(
+                    "``minimize`` option is incompatible with post-processing.")
             if isinstance(output, str) or output is False:
                 info["post"]["output"] = output or None
             return post(info)
@@ -94,6 +98,11 @@ def run(info_or_yaml_or_file: Union[InputDict, str, os.PathLike],
         # GetDist needs to know the original sampler, so don't overwrite if minimizer
         try:
             which_sampler = list(info["sampler"])[0]
+            if minimize:
+                # Preserve options if "minimize" was already the sampler
+                if which_sampler.lower() != "minimize":
+                    info["sampler"] = {"minimize": None}
+                    which_sampler = "minimize"
         except (KeyError, TypeError):
             raise LoggedError(
                 logger_run, "You need to specify a sampler using the 'sampler' key "
@@ -192,6 +201,9 @@ def run_script(args=None):
                                    "(use with care!)")
     parser.add_argument("--%s" % "test", action="store_true",
                         help="Initialize model and sampler, and exit.")
+    parser.add_argument("--minimize", action="store_true",
+                        help=("Replaces the sampler in the input and runs a minimization "
+                              "process (incompatible with post-processing)."))
     parser.add_argument("--version", action="version", version=get_version())
     parser.add_argument("--no-mpi", action='store_true',
                         help="disable MPI when mpi4py installed but MPI does "

cobaya/samplers/minimize/minimize.py

@@ -6,12 +6,10 @@
 
 This is a **maximizer** for posteriors or likelihoods, based on
 `scipy.optimize.Minimize <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html>`_
-and `Py-BOBYQA <https://numericalalgorithmsgroup.github.io/pybobyqa/build/html/index.html>`_
-(added in 2.0).
+and `Py-BOBYQA <https://numericalalgorithmsgroup.github.io/pybobyqa/build/html/index.html>`_.
 
-.. note::
-
-   BOBYQA tends to work better on Cosmological problems with the default settings.
+The default is BOBYQA, which tends to work better on Cosmological problems with default
+settings.
 
 .. |br| raw:: html
 
@@ -36,24 +34,30 @@
    **If you use scipy**, you can find `the appropriate references here
    <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html>`_.
 
-It works more effectively when run on top of a Monte Carlo sample: just change the sampler
-for ``Minimize`` with the desired options, and it will use as a starting point the
-*maximum a posteriori* (MAP) or best fit (maximum likelihood, or minimal :math:`\chi^2`)
-found so far, as well as the covariance matrix of the sample for rescaling of the
-parameter jumps.
+It works more effectively when run on top of a Monte Carlo sample: it will use the maximum
+a posteriori as a starting point (or the best fit, depending on whether the prior is
+ignored, :ref:`see below <minimize_like>`), and the recovered covariance matrix of the
+posterior to rescale the variables.
+
+To take advantage of a previous run with a Monte Carlo sampler, either:
+
+- change the ``sampler`` to ``minimize`` in the input file,
+
+- or, if running from the shell, repeat the ``cobaya-run`` command used for the original
+  run, adding the ``--minimize`` flag.
 
-As text output, it produces two different files:
+When called from a Python script, Cobaya's ``run`` function returns the updated info
+and the products described below in the method
+:func:`samplers.minimize.Minimize.products` (see below).
+
+If text output is requested, it produces two different files:
 
 - ``[output prefix].minimum.txt``, in
   :ref:`the same format as Cobaya samples <output_format>`,
   but containing a single line.
 
 - ``[output prefix].minimum``, the equivalent **GetDist-formatted** file.
 
-If ``ignore_prior: True``, those files are named ``.bestfit[.txt]`` instead of ``minimum``,
-and contain the best-fit (maximum of the likelihood) instead of the MAP
-(maximum of the posterior).
-
 .. warning::
 
    For historical reasons, in the first two lines of the GetDist-formatted output file
@@ -62,10 +66,6 @@
    :math:`-2` times the sum of the individual :math:`\chi^2` (``chi2__``, with double
    underscore) in the table that follows these first lines.
 
-When called from a Python script, Cobaya's ``run`` function returns the updated info
-and the products described below in the method
-:func:`products <samplers.Minimize.Minimize.products>`.
-
 It is recommended to run a couple of parallel MPI processes:
 it will finally pick the best among the results.
 
@@ -77,6 +77,18 @@
    want to refine the convergence parameters (see ``override`` options in the ``yaml``
    below).
 
+
+.. _minimize_like:
+
+Maximizing the likelihood instead of the posterior
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+To maximize the likelihood, add ``ignore_prior: True`` in the ``minimize`` input block.
+
+When producing text output, the generated files are named ``.bestfit[.txt]`` instead of
+``minimum``, and contain the best-fit (maximum of the likelihood) instead of the MAP
+(maximum of the posterior).
+
 """
 
 # Global

cobaya/theories/camb/camb.py

@@ -180,7 +180,7 @@
 from cobaya.log import LoggedError, get_logger
 from cobaya.install import download_github_release, check_gcc_version, NotInstalledError
 from cobaya.tools import getfullargspec, get_class_methods, get_properties, load_module, \
-    VersionCheckError, str_to_list
+    VersionCheckError, str_to_list, Pool1D, Pool2D, PoolND, check_2d
 from cobaya.theory import HelperTheory
 from cobaya.typing import InfoDict
 
@@ -190,6 +190,7 @@ class Collector(NamedTuple):
     method: Callable
     args: list = []
     kwargs: dict = {}
+    z_pool: Optional[PoolND] = None
 
 
 class CAMBOutputs(NamedTuple):
@@ -343,13 +344,17 @@ def must_provide(self, **requirements):
                     method=CAMBdata.get_cmb_power_spectra,
                     kwargs={"spectra": ["unlensed_total"], "raw_cl": False})
             elif k == "Hubble":
-                self.collectors[k] = Collector(
-                    method=CAMBdata.h_of_z,
-                    kwargs={"z": self._combine_z(k, v)})
+                self.set_collector_with_z_pool(k, v["z"], CAMBdata.h_of_z)
+            elif k in ["Omega_b", "Omega_cdm", "Omega_nu_massive"]:
+                varnames = {
+                    "Omega_b": "baryon", "Omega_cdm": "cdm", "Omega_nu_massive": "nu"}
+                self.set_collector_with_z_pool(
+                    k, v["z"], CAMBdata.get_Omega, kwargs={"var": varnames[k]})
             elif k in ("angular_diameter_distance", "comoving_radial_distance"):
-                self.collectors[k] = Collector(
-                    method=getattr(CAMBdata, k),
-                    kwargs={"z": self._combine_z(k, v)})
+                self.set_collector_with_z_pool(k, v["z"], getattr(CAMBdata, k))
+            elif k == "angular_diameter_distance_2":
+                self.set_collector_with_z_pool(
+                    k, v["z_pairs"], CAMBdata.angular_diameter_distance2, d=2)
             elif k == "sigma8_z":
                 self.add_to_redshifts(v["z"])
                 self.collectors[k] = Collector(
@@ -470,16 +475,33 @@ def get_sigmaR(results, **tmp):
         return must_provide
 
     def add_to_redshifts(self, z):
-        self.extra_args["redshifts"] = np.sort(np.unique(np.concatenate(
-            (np.atleast_1d(z), self.extra_args.get("redshifts", [])))))[::-1]
-
-    def _combine_z(self, k, v):
-        c = self.collectors.get(k, None)
-        if c:
-            return np.sort(
-                np.unique(np.concatenate((c.kwargs['z'], np.atleast_1d(v['z'])))))
+        """
+        Adds redshifts to the list of them for which CAMB computes perturbations.
+        """
+        if not hasattr(self, "z_pool_for_perturbations"):
+            self.z_pool_for_perturbations = Pool1D(z)
+        else:
+            self.z_pool_for_perturbations.update(z)
+        self.extra_args["redshifts"] = np.flip(self.z_pool_for_perturbations.values)
+
+    def set_collector_with_z_pool(self, k, zs, method, args=[], kwargs={}, d=1):
+        """
+        Creates a collector for a z-dependent quantity, keeping track of the pool of z's.
+        """
+        if k in self.collectors:
+            z_pool = self.collectors[k].z_pool
+            z_pool.update(zs)
         else:
-            return np.sort(np.atleast_1d(v['z']))
+            Pool = {1: Pool1D, 2: Pool2D}[d]
+            z_pool = Pool(zs)
+        if d == 1:
+            kwargs_with_z = {"z": z_pool.values}
+        else:
+            kwargs_with_z = {"z1": np.array(z_pool.values[:, 0]),
+                             "z2": np.array(z_pool.values[:, 1])}
+        kwargs_with_z.update(kwargs)
+        self.collectors[k] = Collector(
+            method=method, z_pool=z_pool, kwargs=kwargs_with_z, args=args)
 
     def calculate(self, state, want_derived=True, **params_values_dict):
         try:
@@ -630,17 +652,47 @@ def get_unlensed_Cl(self, ell_factor=False, units="FIRASmuK2"):
         return self._get_Cl(ell_factor=ell_factor, units=units, lensed=False)
 
     def _get_z_dependent(self, quantity, z):
-        if quantity in ["sigma8_z", "fsigma8"]:
-            computed_redshifts = self.extra_args["redshifts"]
-            i_kwarg_z = np.concatenate(
-                [np.where(computed_redshifts == zi)[0] for zi in np.atleast_1d(z)])
-        else:
-            computed_redshifts = self.collectors[quantity].kwargs["z"]
-            i_kwarg_z = np.searchsorted(computed_redshifts, np.atleast_1d(z))
+        try:
+            if quantity in ["sigma8_z", "fsigma8"]:
+                pool = self.z_pool_for_perturbations
+                # Mind that CAMB stores redshifts for perturbations in inv order
+                i_kwarg_z_inv = pool.find_indices(z)
+                i_kwarg_z = len(pool) - 1 - i_kwarg_z_inv
+            else:
+                pool = self.collectors[quantity].z_pool
+                i_kwarg_z = pool.find_indices(z)
+        except ValueError:
+            raise LoggedError(self.log, f"{quantity} not computed for all z requested. "
+                                        f"Requested z are {z}, but computed ones are "
+                                        f"{pool.values}.")
         return np.array(self.current_state[quantity], copy=True)[i_kwarg_z]
 
-    def get_sigma8_z(self, z):
-        return self._get_z_dependent("sigma8_z", z)
+    def _get_z_pair_dependent(self, quantity, z_pairs, inv_value=0):
+        """
+        ``inv_value`` (default=0) is assigned to pairs for which ``z1 > z2``.
+        """
+        try:
+            check_2d(z_pairs)
+        except ValueError:
+            raise LoggedError(self.log, f"{z_pairs=} not correctly formatted for "
+                                        f"{quantity}. It should be a list of pairs.")
+        # Only recover for correctly sorted pairs
+        z_pairs_arr = np.array(z_pairs)
+        i_right = (z_pairs_arr[:, 0] <= z_pairs_arr[:, 1])
+        pool = self.collectors[quantity].z_pool
+        i_z_pair = pool.find_indices(z_pairs_arr[i_right])
+        result = np.full(len(z_pairs), inv_value, dtype=float)
+        result[i_right] = np.array(self.current_state[quantity], copy=True)[i_z_pair]
+        return result
+
+    def get_Omega_b(self, z):
+        return self._get_z_dependent("Omega_b", z)
+
+    def get_Omega_cdm(self, z):
+        return self._get_z_dependent("Omega_cdm", z)
+
+    def get_Omega_nu_massive(self, z):
+        return self._get_z_dependent("Omega_nu_massive", z)
 
     def get_fsigma8(self, z):
         return self._get_z_dependent("fsigma8", z)
@@ -812,7 +864,7 @@ def get_path(cls, path):
     def is_installed(cls, **kwargs):
         if not kwargs.get("code", True):
             return True
-        log = get_logger(cls.__name__)
+        log = get_logger(cls.__name__, color=cls._logger_color)
         import platform
         check = kwargs.get("check", True)
         func = log.info if check else log.error