From ecd2319433eca65df079e129b8d621f6f7f8324a Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Tue, 31 Mar 2026 15:31:14 +0100
Subject: [PATCH 01/23] Classy, thanks to Alex Nemtinova

---
 cosmopower/spectra_classy.py | 56 ++++++++++++++++++++++++++++++++++++
 1 file changed, 56 insertions(+)
 create mode 100644 cosmopower/spectra_classy.py

diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
new file mode 100644
index 0000000..9f571d5
--- /dev/null
+++ b/cosmopower/spectra_classy.py
@@ -0,0 +1,56 @@
+from .parser import YAMLParser
+import numpy as np
+from classy import Class
+
+
+def initialize(parser, extra_args={}):
+
+    cosmo = Class()
+
+    state = {}
+
+    for quantity in parser.quantities:
+        state[quantity + ".modes"] = parser.modes(quantity)
+
+    state["code"] = "classy"
+    state["module"] = cosmo
+    state["params"] = extra_args
+    state["derived"] = {}
+
+    return state
+
+
+def get_spectra(parser, state, args={}, quantities=[], extra_args={}):
+
+    cosmo = state["module"]
+
+    params = state["params"].copy()
+    params.update(args)
+
+    if "ln10^{10}A_s" in params:
+
+        params["A_s"] = 1.e-10 * np.exp(params["ln10^{10}A_s"])
+
+        del params["ln10^{10}A_s"]
+
+    cosmo.set(params)
+
+    cosmo.compute()
+
+    cls = cosmo.raw_cl()
+    ell = cls["ell"]
+
+    for quantity in quantities:
+
+        spec = quantity.split("/")[1]
+
+        Cl = cls[spec]
+
+        Dl = Cl * (ell * (ell + 1) / (2 * np.pi))
+
+        state[quantity] = Dl[state[quantity + ".modes"]]
+
+    cosmo.struct_cleanup()
+    cosmo.empty()
+
+    return True

From 13a3d84141cce6d3a27af9685042632658159f59 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Tue, 31 Mar 2026 16:20:23 +0100
Subject: [PATCH 02/23] Update deprecation of np.in1d for numpy 2.2

---
 cosmopower/dataset.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/cosmopower/dataset.py b/cosmopower/dataset.py
index 76b75f6..f86cd0d 100644
--- a/cosmopower/dataset.py
+++ b/cosmopower/dataset.py
@@ -170,7 +170,7 @@ def read_parameters(self, indices: Union[int, np.ndarray]) -> dict:
         cast_1d = (type(indices) is int)
         indices = np.atleast_1d(indices)
 
-        entries = np.where(np.in1d(self.file["data"]["indices"][:],
+        entries = np.where(np.isin(self.file["data"]["indices"][:],
                                    indices))[0]
 
         parameters = {
@@ -179,7 +179,7 @@ def read_parameters(self, indices: Union[int, np.ndarray]) -> dict:
         }
 
         if cast_1d:
-            parameters = {k: float(v) for k, v in parameters.items()}
+            parameters = {k: float(v[0]) for k, v in parameters.items()}
 
         return parameters
 
@@ -193,7 +193,7 @@ def read_spectra(self, indices: Union[int, np.ndarray]) -> np.ndarray:
         cast_1d = (type(indices) is int)
         indices = np.atleast_1d(indices)
 
-        entries = np.where(np.in1d(self.file["data"]["indices"][:],
+        entries = np.where(np.isin(self.file["data"]["indices"][:],
                                    indices))[0]
         spec = self.file["data"]["spectra"][entries, :]
         if cast_1d:

From 134bd4250c3a718cc094e4ab302e86208a339c25 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Tue, 31 Mar 2026 16:23:48 +0100
Subject: [PATCH 03/23] Basic management of classy module akin to camb
 interface.

---
 cosmopower/spectra_classy.py | 55 ++++++++++++++++++------------------
 1 file changed, 27 insertions(+), 28 deletions(-)

diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index 9f571d5..56ea839 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -1,54 +1,53 @@
 from .parser import YAMLParser
 import numpy as np
-from classy import Class
 
 
-def initialize(parser, extra_args={}):
-
-    cosmo = Class()
+def initialize(parser: YAMLParser, extra_args: dict = {}) -> dict:
+    import classy
 
     state = {}
 
     for quantity in parser.quantities:
         state[quantity + ".modes"] = parser.modes(quantity)
 
+    extra_args["l_max_scalars"] = extra_args.get("l_max_scalars", 2)
+    for quantity in parser.quantities:
+        if parser.modes_label(quantity) == "l" and \
+           extra_args["l_max_scalars"] < parser.modes(quantity).max():
+           extra_args["l_max_scalars"] = parser.modes(quantity).max()
+
     state["code"] = "classy"
-    state["module"] = cosmo
+    state["module"] = classy
     state["params"] = extra_args
     state["derived"] = {}
 
     return state
 
 
-def get_spectra(parser, state, args={}, quantities=[], extra_args={}):
-
-    cosmo = state["module"]
-
-    params = state["params"].copy()
-    params.update(args)
-
-    if "ln10^{10}A_s" in params:
+def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
+                quantities: list = [], extra_args: dict = {}) -> bool:
 
-        params["A_s"] = 1.e-10 * np.exp(params["ln10^{10}A_s"])
+    classy = state["module"]
 
-        del params["ln10^{10}A_s"]
-
-    cosmo.set(params)
-
-    cosmo.compute()
-
-    cls = cosmo.raw_cl()
-    ell = cls["ell"]
+    try:
+        cosmo = classy.Class()
+        cosmo.set(state["params"] | args)
+        cosmo.compute()
+    except classy.CosmoError as e:
+        print(e)
+        return False
 
     for quantity in quantities:
+        qpath = quantity.split("/")
+        if qpath[0] == "Cl":
+            spec = qpath[1]
 
-        spec = quantity.split("/")[1]
-
-        Cl = cls[spec]
-
-        Dl = Cl * (ell * (ell + 1) / (2 * np.pi))
+            cls = cosmo.raw_cl()
+            ell = cls["ell"]
+            Cl = cls[spec]
+            Dl = Cl * (ell * (ell + 1) / (2 * np.pi))
 
-        state[quantity] = Dl[state[quantity + ".modes"]]
+            state[quantity] = Dl[state[quantity + ".modes"]]
 
     cosmo.struct_cleanup()
     cosmo.empty()

From 0662d7b0312ce417b689305b3cd65c15f65b9a47 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Tue, 31 Mar 2026 16:24:50 +0100
Subject: [PATCH 04/23] Cl now calls lensed_cl, moved raw_cl to Cl_unlensed.

---
 cosmopower/spectra_classy.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index 56ea839..f19347f 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -42,6 +42,15 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
         if qpath[0] == "Cl":
             spec = qpath[1]
 
+            cls = cosmo.lensed_cl()
+            ell = cls["ell"]
+            Cl = cls[spec]
+            Dl = Cl * (ell * (ell + 1) / (2 * np.pi))
+
+            state[quantity] = Dl[state[quantity + ".modes"]]
+        elif qpath[0] == "Cl_unlensed":
+            spec = qpath[1]
+
             cls = cosmo.raw_cl()
             ell = cls["ell"]
             Cl = cls[spec]

From fe8946fd91e5710c01d5a64c769100ebb606a1c1 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Tue, 31 Mar 2026 17:50:32 +0100
Subject: [PATCH 05/23] Unlensed_Cl not Cl_unlensed (for cobaya)

---
 cosmopower/spectra_classy.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index f19347f..ec93f8d 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -48,7 +48,7 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
             Dl = Cl * (ell * (ell + 1) / (2 * np.pi))
 
             state[quantity] = Dl[state[quantity + ".modes"]]
-        elif qpath[0] == "Cl_unlensed":
+        elif qpath[0] == "unlensed_Cl":
             spec = qpath[1]
 
             cls = cosmo.raw_cl()

From af7e007429cf15dd7050be17275228ededa2fd8c Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Wed, 1 Apr 2026 09:42:45 +0100
Subject: [PATCH 06/23] Using the correct error type.

---
 cosmopower/spectra_classy.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index ec93f8d..a57cda8 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -33,8 +33,7 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
         cosmo = classy.Class()
         cosmo.set(state["params"] | args)
         cosmo.compute()
-    except classy.CosmoError as e:
-        print(e)
+    except (classy.CosmoComputationError, classy.CosmoSevereError) as e:
         return False
 
     for quantity in quantities:

From 5cc31af7a2b729a7a0ea0b569dced239feb291d5 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Wed, 1 Apr 2026 09:56:24 +0100
Subject: [PATCH 07/23] Scalar vs Tensor lmax parsing.

---
 cosmopower/spectra_camb.py   | 12 ++++++------
 cosmopower/spectra_classy.py | 31 +++++++++++++++++++++++++++----
 2 files changed, 33 insertions(+), 10 deletions(-)

diff --git a/cosmopower/spectra_camb.py b/cosmopower/spectra_camb.py
index 7411cd0..311c4d1 100644
--- a/cosmopower/spectra_camb.py
+++ b/cosmopower/spectra_camb.py
@@ -48,13 +48,13 @@ def initialize(parser: YAMLParser, extra_args: dict = {}) -> dict:
 
         # Check that the maximum values for mode ranges matches
         # the requested outputs.
-        if parser.modes_label(quantity) == "l" and \
-           extra_args["lmax"] < parser.modes(quantity).max():
-            extra_args["lmax"] = parser.modes(quantity).max()
+        if parser.modes_label(quantity) == "l":
+           extra_args["lmax"] = max(extra_args["lmax"],
+                                    parser.modes(quantity).max())
 
-        if parser.modes_label(quantity) == "k" and \
-           extra_args["kmax"] < parser.modes(quantity).max():
-            extra_args["kmax"] = parser.modes(quantity).max()
+        if parser.modes_label(quantity) == "k":
+           extra_args["kmax"] = max(extra_args["kmax"],
+                                    parser.modes(quantity).max())
 
         if parser.modes_label(quantity) == "z":
             z1 = extra_args["redshifts"]
diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index a57cda8..dfd051a 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -10,11 +10,24 @@ def initialize(parser: YAMLParser, extra_args: dict = {}) -> dict:
     for quantity in parser.quantities:
         state[quantity + ".modes"] = parser.modes(quantity)
 
-    extra_args["l_max_scalars"] = extra_args.get("l_max_scalars", 2)
+    extra_args["l_max_scalars"] = extra_args.get("l_max_scalars", 0)
+    extra_args["l_max_tensors"] = extra_args.get("l_max_tensors", 0)
+
+    cl_modes = []
+
     for quantity in parser.quantities:
-        if parser.modes_label(quantity) == "l" and \
-           extra_args["l_max_scalars"] < parser.modes(quantity).max():
-           extra_args["l_max_scalars"] = parser.modes(quantity).max()
+        if quantity.split("/")[0] == "Cl":
+            cl_modes.append("s")
+            extra_args["l_max_scalars"] = max(extra_args["l_max_scalars"],
+                                              parser.modes(quantity).max())
+        if quantity.split("/")[0] == "tensor_Cl":
+            cl_modes.append("t")
+            extra_args["l_max_tensors"] = max(extra_args["l_max_tensors"],
+                                              parser.modes(quantity).max())
+
+    if "s" not in cl_modes: extra_args.pop("l_max_scalars")
+    if "t" not in cl_modes: extra_args.pop("l_max_tensors")
+    extra_args["modes"] = " ".join(cl_modes)
 
     state["code"] = "classy"
     state["module"] = classy
@@ -34,6 +47,7 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
         cosmo.set(state["params"] | args)
         cosmo.compute()
     except (classy.CosmoComputationError, classy.CosmoSevereError) as e:
+        print(str(e))
         return False
 
     for quantity in quantities:
@@ -55,6 +69,15 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
             Cl = cls[spec]
             Dl = Cl * (ell * (ell + 1) / (2 * np.pi))
 
+            state[quantity] = Dl[state[quantity + ".modes"]]
+        elif qpath[0] == "tensor_Cl":
+            spec = qpath[1]
+
+            cls = cosmo.raw_cl()
+            ell = cls["ell"]
+            Cl = cls[spec]
+            Dl = Cl * (ell * (ell + 1) / (2 * np.pi))
+
             state[quantity] = Dl[state[quantity + ".modes"]]
 
     cosmo.struct_cleanup()

From 930da1a0ca93aea20617695815210d6369d14249 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Wed, 1 Apr 2026 10:53:57 +0100
Subject: [PATCH 08/23] Optional automatic class parameters.

---
 cosmopower/spectra.py                    |  4 +-
 cosmopower/spectra_classy.py             | 54 +++++++++++++++++-------
 cosmopower/wrappers/cobaya/cosmopower.py |  1 -
 3 files changed, 41 insertions(+), 18 deletions(-)

diff --git a/cosmopower/spectra.py b/cosmopower/spectra.py
index 8947abd..1f6a154 100644
--- a/cosmopower/spectra.py
+++ b/cosmopower/spectra.py
@@ -59,7 +59,7 @@ def init_boltzmann_code(parser: YAMLParser) -> dict:
     res = {}
 
     code = parser.boltzmann_code
-    extra_args = parser.boltzmann_extra_args
+    extra_args = parser.boltzmann_extra_args or {}
 
     # Try to import the correct python module for spectrum generation
     try:
@@ -256,7 +256,7 @@ def generate_spectra(args: list = None) -> None:
           f"files {first_file}--{last_file}.")
 
     state = init_boltzmann_code(parser)
-    extra_args = parser.boltzmann_extra_args
+    extra_args = parser.boltzmann_extra_args or {}
 
     accepted = 0
     tbar = tqdm.tqdm(np.arange(first, last + 1))
diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index dfd051a..27229a3 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -10,32 +10,57 @@ def initialize(parser: YAMLParser, extra_args: dict = {}) -> dict:
     for quantity in parser.quantities:
         state[quantity + ".modes"] = parser.modes(quantity)
 
+    if extra_args.get("auto", False):
+        extra_args.pop("auto")
+        extra_args = classy_args_interpret(parser, extra_args)
+        print(f"Auto-generated class settings: {extra_args}")
+
+    state["code"] = "classy"
+    state["module"] = classy
+    state["params"] = extra_args
+    state["derived"] = {}
+
+    return state
+
+
+def classy_args_interpret(parser: YAMLParser, extra_args: dict = {}) -> dict:
     extra_args["l_max_scalars"] = extra_args.get("l_max_scalars", 0)
     extra_args["l_max_tensors"] = extra_args.get("l_max_tensors", 0)
 
-    cl_modes = []
+    extra_args["lensing"] = "no"
+    cl_spec = set([])
+    cl_modes = set([])
 
     for quantity in parser.quantities:
-        if quantity.split("/")[0] == "Cl":
-            cl_modes.append("s")
+        qpath = quantity.split("/")
+        if qpath[0] == "Cl" or qpath[0] == "unlensed_Cl":
+            if qpath[1] == "tt" or qpath[1] == "te":
+                cl_spec.add("tCl")
+            if qpath[1] == "te" or qpath[1] == "ee":
+                cl_spec.add("pCl")
+            if qpath[0] == "Cl" or qpath[1] == "pp":
+                extra_args["lensing"] = "yes"
+                cl_spec.add("lCl")
+
+            cl_modes.add("s")
             extra_args["l_max_scalars"] = max(extra_args["l_max_scalars"],
                                               parser.modes(quantity).max())
-        if quantity.split("/")[0] == "tensor_Cl":
-            cl_modes.append("t")
+        elif qpath[0] == "tensor_Cl":
+            if qpath[1] == "tt" or qpath[1] == "te":
+                cl_spec.add("tCl")
+            if qpath[1] == "te" or qpath[1] == "ee":
+                cl_spec.add("pCl")
+
+            cl_modes.add("t")
             extra_args["l_max_tensors"] = max(extra_args["l_max_tensors"],
                                               parser.modes(quantity).max())
 
     if "s" not in cl_modes: extra_args.pop("l_max_scalars")
     if "t" not in cl_modes: extra_args.pop("l_max_tensors")
-    extra_args["modes"] = " ".join(cl_modes)
-
-    state["code"] = "classy"
-    state["module"] = classy
-    state["params"] = extra_args
-    state["derived"] = {}
-
-    return state
+    extra_args["output"] = " ".join(list(cl_spec))
+    extra_args["modes"] = " ".join(list(cl_modes))
 
+    return extra_args
 
 def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
                 quantities: list = [], extra_args: dict = {}) -> bool:
@@ -47,8 +72,7 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
         cosmo.set(state["params"] | args)
         cosmo.compute()
     except (classy.CosmoComputationError, classy.CosmoSevereError) as e:
-        print(str(e))
-        return False
+        raise e
 
     for quantity in quantities:
         qpath = quantity.split("/")
diff --git a/cosmopower/wrappers/cobaya/cosmopower.py b/cosmopower/wrappers/cobaya/cosmopower.py
index 645af7a..77c0beb 100644
--- a/cosmopower/wrappers/cobaya/cosmopower.py
+++ b/cosmopower/wrappers/cobaya/cosmopower.py
@@ -284,7 +284,6 @@ def get_Cl(self, ell_factor: bool = False,
 
 
     def get_requirements(self):
-
         return self._params_required_by_networks
 
     def get_in_dict(self, dct: dict, path: str) -> Any:

From 15b90ba75f22e8dd3e46353baaf46a1cb47aa762 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Wed, 1 Apr 2026 11:18:00 +0100
Subject: [PATCH 09/23] Renamed tensor_Cl -> unlensed_Cl for cobaya parity.

---
 cosmopower/spectra_classy.py | 33 +++++++++------------------------
 1 file changed, 9 insertions(+), 24 deletions(-)

diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index 27229a3..8e34b46 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -28,8 +28,8 @@ def classy_args_interpret(parser: YAMLParser, extra_args: dict = {}) -> dict:
     extra_args["l_max_tensors"] = extra_args.get("l_max_tensors", 0)
 
     extra_args["lensing"] = "no"
-    cl_spec = set([])
-    cl_modes = set([])
+    cl_spec = set(extra_args.get("output", "").split(" "))
+    cl_modes = set(extra_args.get("modes", "s").split(" "))
 
     for quantity in parser.quantities:
         qpath = quantity.split("/")
@@ -42,18 +42,12 @@ def classy_args_interpret(parser: YAMLParser, extra_args: dict = {}) -> dict:
                 extra_args["lensing"] = "yes"
                 cl_spec.add("lCl")
 
-            cl_modes.add("s")
-            extra_args["l_max_scalars"] = max(extra_args["l_max_scalars"],
-                                              parser.modes(quantity).max())
-        elif qpath[0] == "tensor_Cl":
-            if qpath[1] == "tt" or qpath[1] == "te":
-                cl_spec.add("tCl")
-            if qpath[1] == "te" or qpath[1] == "ee":
-                cl_spec.add("pCl")
-
-            cl_modes.add("t")
-            extra_args["l_max_tensors"] = max(extra_args["l_max_tensors"],
-                                              parser.modes(quantity).max())
+            if "s" in cl_modes:
+                extra_args["l_max_scalars"] = max(extra_args["l_max_scalars"],
+                                                  parser.modes(quantity).max())
+            if "t" in cl_modes:
+                extra_args["l_max_tensors"] = max(extra_args["l_max_tensors"],
+                                                  parser.modes(quantity).max())
 
     if "s" not in cl_modes: extra_args.pop("l_max_scalars")
     if "t" not in cl_modes: extra_args.pop("l_max_tensors")
@@ -85,16 +79,7 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
             Dl = Cl * (ell * (ell + 1) / (2 * np.pi))
 
             state[quantity] = Dl[state[quantity + ".modes"]]
-        elif qpath[0] == "unlensed_Cl":
-            spec = qpath[1]
-
-            cls = cosmo.raw_cl()
-            ell = cls["ell"]
-            Cl = cls[spec]
-            Dl = Cl * (ell * (ell + 1) / (2 * np.pi))
-
-            state[quantity] = Dl[state[quantity + ".modes"]]
-        elif qpath[0] == "tensor_Cl":
+        elif qpath[0] == "unlensed_Cl" or qpath[0] == "tensor_Cl":
             spec = qpath[1]
 
             cls = cosmo.raw_cl()

From b438a0c29096d63bfb63b2079df478c2836c7e16 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Wed, 1 Apr 2026 11:19:58 +0100
Subject: [PATCH 10/23] cleanup.

---
 cosmopower/spectra_classy.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index 8e34b46..2d75834 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -79,7 +79,7 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
             Dl = Cl * (ell * (ell + 1) / (2 * np.pi))
 
             state[quantity] = Dl[state[quantity + ".modes"]]
-        elif qpath[0] == "unlensed_Cl" or qpath[0] == "tensor_Cl":
+        elif qpath[0] == "unlensed_Cl":
             spec = qpath[1]
 
             cls = cosmo.raw_cl()

From b16c03a5c1281d36db72b394f450566f8b0e26cf Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Wed, 1 Apr 2026 15:20:35 +0100
Subject: [PATCH 11/23] Derived parameters from class.

---
 cosmopower/examples/planck_cmb.yaml |  4 ---
 cosmopower/spectra_camb.py          | 14 +++++------
 cosmopower/spectra_classy.py        | 39 ++++++++++++++++++++++++-----
 3 files changed, 40 insertions(+), 17 deletions(-)

diff --git a/cosmopower/examples/planck_cmb.yaml b/cosmopower/examples/planck_cmb.yaml
index d0f30a4..b1493b0 100644
--- a/cosmopower/examples/planck_cmb.yaml
+++ b/cosmopower/examples/planck_cmb.yaml
@@ -40,10 +40,6 @@ networks:
   - quantity: "derived"
     type: NN
     inputs: [ombh2, omch2, cosmomc_theta, ns, logA, tau]
-    log: True
-    modes:
-      label: l
-      range: [2, 2508]
     n_traits:
       n_hidden: [512,512,512,512]
     training:
diff --git a/cosmopower/spectra_camb.py b/cosmopower/spectra_camb.py
index 311c4d1..8262714 100644
--- a/cosmopower/spectra_camb.py
+++ b/cosmopower/spectra_camb.py
@@ -49,12 +49,12 @@ def initialize(parser: YAMLParser, extra_args: dict = {}) -> dict:
         # Check that the maximum values for mode ranges matches
         # the requested outputs.
         if parser.modes_label(quantity) == "l":
-           extra_args["lmax"] = max(extra_args["lmax"],
-                                    parser.modes(quantity).max())
+            extra_args["lmax"] = max(extra_args["lmax"],
+                                     parser.modes(quantity).max())
 
         if parser.modes_label(quantity) == "k":
-           extra_args["kmax"] = max(extra_args["kmax"],
-                                    parser.modes(quantity).max())
+            extra_args["kmax"] = max(extra_args["kmax"],
+                                     parser.modes(quantity).max())
 
         if parser.modes_label(quantity) == "z":
             z1 = extra_args["redshifts"]
@@ -223,12 +223,12 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
     except (camb.CAMBError, camb.CAMBFortranError):
         return False
 
+    state["derived"] = get_camb_derived(state["params"], state["results"],
+                                        list(state["derived"].keys()))
+
     for quantity in quantities:
         qpath = quantity.split("/")
 
-        state["derived"] = get_camb_derived(state["params"], state["results"],
-                                            list(state["derived"].keys()))
-
         if qpath[0] == "derived":
             continue
         elif qpath[0] == "Cl":
diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index 2d75834..63772bb 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -1,10 +1,11 @@
+from __future__ import annotations
 from .parser import YAMLParser
 import numpy as np
+import classy
+from typing import Sequence
 
 
 def initialize(parser: YAMLParser, extra_args: dict = {}) -> dict:
-    import classy
-
     state = {}
 
     for quantity in parser.quantities:
@@ -18,7 +19,7 @@ def initialize(parser: YAMLParser, extra_args: dict = {}) -> dict:
     state["code"] = "classy"
     state["module"] = classy
     state["params"] = extra_args
-    state["derived"] = {}
+    state["derived"] = {k: np.nan for k in parser.computed_parameters}
 
     return state
 
@@ -49,13 +50,33 @@ def classy_args_interpret(parser: YAMLParser, extra_args: dict = {}) -> dict:
                 extra_args["l_max_tensors"] = max(extra_args["l_max_tensors"],
                                                   parser.modes(quantity).max())
 
-    if "s" not in cl_modes: extra_args.pop("l_max_scalars")
-    if "t" not in cl_modes: extra_args.pop("l_max_tensors")
+    if "s" not in cl_modes:
+        extra_args.pop("l_max_scalars")
+    if "t" not in cl_modes:
+        extra_args.pop("l_max_tensors")
     extra_args["output"] = " ".join(list(cl_spec))
     extra_args["modes"] = " ".join(list(cl_modes))
 
     return extra_args
 
+
+def get_classy_derived(params: dict, cosmo: classy.Class,
+                       parameters: Sequence[str]) -> dict:
+    res = {}
+
+    for p in parameters:
+        if p == "rs_drag":
+            res[p] = cosmo.rs_drag()
+        elif p == "Omega_nu":
+            res[p] = cosmo.Omega_nu
+        elif p == "T_cmb":
+            res[p] = cosmo.T_cmb()
+        else:
+            res[p] = cosmo.get_current_derived_parameters([p])[p]
+
+    return res
+
+
 def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
                 quantities: list = [], extra_args: dict = {}) -> bool:
 
@@ -68,9 +89,15 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
     except (classy.CosmoComputationError, classy.CosmoSevereError) as e:
         raise e
 
+    state["derived"] = get_classy_derived(state["params"], cosmo,
+                                          list(state["derived"].keys()))
+
     for quantity in quantities:
         qpath = quantity.split("/")
-        if qpath[0] == "Cl":
+
+        if qpath[0] == "derived":
+            continue
+        elif qpath[0] == "Cl":
             spec = qpath[1]
 
             cls = cosmo.lensed_cl()

From 24eba47367b3a575a45e01f30e537e44ed2182de Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Wed, 1 Apr 2026 15:22:54 +0100
Subject: [PATCH 12/23] Auto inference of mPk requested for sigma8 derived.

---
 cosmopower/spectra_classy.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index 63772bb..d57a43d 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -33,6 +33,10 @@ def classy_args_interpret(parser: YAMLParser, extra_args: dict = {}) -> dict:
     cl_modes = set(extra_args.get("modes", "s").split(" "))
 
     for quantity in parser.quantities:
+        if quantity == "derived":
+            if "sigma8" in parser.computed_parameters:
+                cl_spec.add("mPk")
+
         qpath = quantity.split("/")
         if qpath[0] == "Cl" or qpath[0] == "unlensed_Cl":
             if qpath[1] == "tt" or qpath[1] == "te":

From 5fc7c4a975d08d02b58a4d6ca23820e8c65f7e80 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Thu, 9 Apr 2026 10:35:36 +0100
Subject: [PATCH 13/23] Updated weight inclusion for tf>2.15

---
 cosmopower/cosmopower_NN.py                   | 41 +++++++++++--------
 cosmopower/cosmopower_PCAplusNN.py            | 39 +++++++++++-------
 cosmopower/examples/1_create_training_data.py | 11 ++---
 cosmopower/examples/2b_train_PCA_emulator.py  |  2 +-
 cosmopower/train.py                           |  8 ++++
 5 files changed, 64 insertions(+), 37 deletions(-)

diff --git a/cosmopower/cosmopower_NN.py b/cosmopower/cosmopower_NN.py
index 94e089d..b74dfc0 100644
--- a/cosmopower/cosmopower_NN.py
+++ b/cosmopower/cosmopower_NN.py
@@ -118,22 +118,32 @@ def __init__(self, parameters: Optional[Sequence[str]] = None,
         self.alphas = []
         self.betas = []
         for i in range(self.n_layers):
-            self.W.append(tf.Variable(
-                tf.random.normal([
-                    self.architecture[i], self.architecture[i + 1]
-                ], 0.0, 1e-3),
-                name="W_" + str(i), trainable=trainable))
-            self.b.append(tf.Variable(tf.zeros([self.architecture[i + 1]]),
-                                      name="b_" + str(i),
-                                      trainable=trainable))
+            self.W.append(self.add_weight(
+                initializer=tf.keras.initializers.RandomNormal(stddev=1e-3),
+                shape=(self.architecture[i], self.architecture[i + 1]),
+                name="W_" + str(i),
+                trainable=trainable
+            ))
+            self.b.append(self.add_weight(
+                initializer="zeros",
+                shape=(self.architecture[i + 1],),
+                name="b_" + str(i),
+                trainable=trainable
+            ))
 
         for i in range(self.n_layers - 1):
-            self.alphas.append(tf.Variable(
-                tf.random.normal([self.architecture[i + 1]]),
-                name="alphas_" + str(i), trainable=trainable))
-            self.betas.append(tf.Variable(
-                tf.random.normal([self.architecture[i + 1]]),
-                name="betas_" + str(i), trainable=trainable))
+            self.alphas.append(self.add_weight(
+                initializer="random_normal",
+                shape=(self.architecture[i + 1],),
+                name="alphas_" + str(i),
+                trainable=trainable
+            ))
+            self.betas.append(self.add_weight(
+                initializer="random_normal",
+                shape=(self.architecture[i + 1],),
+                name="betas_" + str(i),
+                trainable=trainable
+            ))
 
         # restore weights if restoring
         if restore_filename is not None:
@@ -546,7 +556,6 @@ def ten_to_rescaled_predictions_np(self, parameters_dict: dict
         return 10.**self.rescaled_predictions_np(parameters_dict)
 
     # Infrastructure for network training
-
     @tf.function
     def compute_loss(self, training_parameters: tf.Tensor,
                      training_features: tf.Tensor) -> tf.Tensor:
@@ -595,7 +604,6 @@ def compute_loss_and_gradients(self,
         """
         # compute loss on the tape
         with tf.GradientTape() as tape:
-
             # loss
             loss = tf.sqrt(
                 tf.reduce_mean(
@@ -888,7 +896,6 @@ def train(self, training_data: Sequence,
                 for epoch in t:
                     # loop over batches
                     for theta, feats in training_data:
-
                         # training step: check whether to accumulate gradients
                         # or not (only worth doing this for very large batch
                         # sizes)
diff --git a/cosmopower/cosmopower_PCAplusNN.py b/cosmopower/cosmopower_PCAplusNN.py
index abc8c4b..cabe407 100644
--- a/cosmopower/cosmopower_PCAplusNN.py
+++ b/cosmopower/cosmopower_PCAplusNN.py
@@ -128,21 +128,32 @@ def __init__(self, cp_pca: Optional[object] = None,
         self.alphas = []
         self.betas = []
         for i in range(self.n_layers):
-            self.W.append(tf.Variable(
-                tf.random.normal([
-                    self.architecture[i], self.architecture[i + 1]
-                ], 0.0, np.sqrt(2.0 / self.n_parameters)),
-                name="W_" + str(i), trainable=trainable))
-            self.b.append(tf.Variable(tf.zeros([self.architecture[i + 1]]),
-                                      name="b_" + str(i),
-                                      trainable=trainable))
+            self.W.append(self.add_weight(
+                initializer=tf.keras.initializers.RandomNormal(stddev=1e-3),
+                shape=(self.architecture[i], self.architecture[i + 1]),
+                name="W_" + str(i),
+                trainable=trainable
+            ))
+            self.b.append(self.add_weight(
+                initializer="zeros",
+                shape=(self.architecture[i + 1],),
+                name="b_" + str(i),
+                trainable=trainable
+            ))
+
         for i in range(self.n_layers - 1):
-            self.alphas.append(tf.Variable(
-                tf.random.normal([self.architecture[i + 1]]),
-                name="alphas_" + str(i), trainable=trainable))
-            self.betas.append(tf.Variable(
-                tf.random.normal([self.architecture[i + 1]]),
-                name="betas_" + str(i), trainable=trainable))
+            self.alphas.append(self.add_weight(
+                initializer="random_normal",
+                shape=(self.architecture[i + 1],),
+                name="alphas_" + str(i),
+                trainable=trainable
+            ))
+            self.betas.append(self.add_weight(
+                initializer="random_normal",
+                shape=(self.architecture[i + 1],),
+                name="betas_" + str(i),
+                trainable=trainable
+            ))
 
         # restore weights if restoring
         if restore_filename is not None:
diff --git a/cosmopower/examples/1_create_training_data.py b/cosmopower/examples/1_create_training_data.py
index ce695c0..981f5b3 100644
--- a/cosmopower/examples/1_create_training_data.py
+++ b/cosmopower/examples/1_create_training_data.py
@@ -2,7 +2,7 @@
 
 # Author: A. Spurio Mancini, H. T. Jense
 
-import tqdm
+from tqdm import tqdm
 import numpy as np
 from cosmopower.parser import YAMLParser
 from cosmopower.spectra import init_boltzmann_code, get_boltzmann_spectra
@@ -51,7 +51,8 @@
 samples, validation_samples, testing_samples = \
     parser.get_parameter_samples()
 
-print(samples)
+for k in samples:
+    print(k, samples[k][:10])
 
 parser.save_samples_to_file(samples, validation_samples)
 
@@ -128,7 +129,7 @@
             # This function returns False if the sample is invalid, so it's a
             # simple check to discard invalid datapoints.
             network_params = np.array([
-                params[k]
+                samples[k][n]
                 for k in parser.network_input_parameters("Pk/lin")
             ])
             dataset.write_data(n, network_params,
@@ -147,7 +148,7 @@
             # This function returns False if the sample is invalid, so it's a
             # simple check to discard invalid datapoints.
             network_params = np.array([
-                params[k]
+                samples[k][n]
                 for k in parser.network_input_parameters("Pk/lin")
             ])
             dataset.write_data(n, network_params,
@@ -163,7 +164,7 @@
             # This function returns False if the sample is invalid, so it's a
             # simple check to discard invalid datapoints.
             network_params = np.array([
-                params[k]
+                samples[k][n]
                 for k in parser.network_input_parameters("Pk/lin")
             ])
             dataset.write_data(n, network_params,
diff --git a/cosmopower/examples/2b_train_PCA_emulator.py b/cosmopower/examples/2b_train_PCA_emulator.py
index 1cf09ce..f3b71ea 100644
--- a/cosmopower/examples/2b_train_PCA_emulator.py
+++ b/cosmopower/examples/2b_train_PCA_emulator.py
@@ -42,7 +42,7 @@
 network = cosmopower_PCAplusNN(cp_pca=pca, verbose=True, trainable=True,
                                **settings.get("n_traits", {}))
 
-with tf.device("/device:CPU:0"):
+with tf.device(None):
     network.train(training_data=datasets,
                   filename_saved_model=output_file,
                   validation=validation,
diff --git a/cosmopower/train.py b/cosmopower/train.py
index 6c4127a..68859c0 100644
--- a/cosmopower/train.py
+++ b/cosmopower/train.py
@@ -46,6 +46,14 @@ def train_network_NN(parser: YAMLParser, quantity: str, device: str = "",
                             modes=parser.modes(quantity), verbose=True,
                             **settings.get("n_traits", {}))
 
+    print(network.trainable)
+    print("_layers=", network._layers)
+    print("_non_trainable_variables=", network._non_trainable_variables)
+    print("_trainable_variables=", network._trainable_variables)
+    print("layers=", network.layers)
+    print("variables=", network.variables)
+    print("trainable_variables=", network.trainable_variables)
+    
     datasets = [Dataset(parser, quantity, os.path.basename(filename))
                 for filename in filenames]
     training_params = parser.network_training_parameters(quantity)

From 307a41c4a97e4f4c372fbdd97f39e01dd7e18768 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Thu, 9 Apr 2026 10:37:33 +0100
Subject: [PATCH 14/23] Remove debug prints.

---
 cosmopower/train.py | 8 --------
 1 file changed, 8 deletions(-)

diff --git a/cosmopower/train.py b/cosmopower/train.py
index 68859c0..6c4127a 100644
--- a/cosmopower/train.py
+++ b/cosmopower/train.py
@@ -46,14 +46,6 @@ def train_network_NN(parser: YAMLParser, quantity: str, device: str = "",
                             modes=parser.modes(quantity), verbose=True,
                             **settings.get("n_traits", {}))
 
-    print(network.trainable)
-    print("_layers=", network._layers)
-    print("_non_trainable_variables=", network._non_trainable_variables)
-    print("_trainable_variables=", network._trainable_variables)
-    print("layers=", network.layers)
-    print("variables=", network.variables)
-    print("trainable_variables=", network.trainable_variables)
-    
     datasets = [Dataset(parser, quantity, os.path.basename(filename))
                 for filename in filenames]
     training_params = parser.network_training_parameters(quantity)

From 918332839c248159c09cce8dd20caaf9143c6e97 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Thu, 9 Apr 2026 11:45:19 +0100
Subject: [PATCH 15/23] Added linear matter power spectrum.

---
 cosmopower/examples/example.yaml |  2 +-
 cosmopower/spectra_classy.py     | 34 ++++++++++++++++++++++++++++----
 2 files changed, 31 insertions(+), 5 deletions(-)

diff --git a/cosmopower/examples/example.yaml b/cosmopower/examples/example.yaml
index 4ba920b..9dd3e26 100644
--- a/cosmopower/examples/example.yaml
+++ b/cosmopower/examples/example.yaml
@@ -4,7 +4,7 @@ path: example
 emulated_code:
   name: camb
   # which parameters are passed to camb
-  inputs: [ombh2, omch2, H0, ns, As]
+  inputs: [ombh2, omch2, H0, ns, As, z]
   extra_args:
     # Some extra accuracy and speed parameters
     kmax: 20.0
diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index d57a43d..b60c8b0 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -31,6 +31,7 @@ def classy_args_interpret(parser: YAMLParser, extra_args: dict = {}) -> dict:
     extra_args["lensing"] = "no"
     cl_spec = set(extra_args.get("output", "").split(" "))
     cl_modes = set(extra_args.get("modes", "s").split(" "))
+    want_cl = False
 
     for quantity in parser.quantities:
         if quantity == "derived":
@@ -39,6 +40,7 @@ def classy_args_interpret(parser: YAMLParser, extra_args: dict = {}) -> dict:
 
         qpath = quantity.split("/")
         if qpath[0] == "Cl" or qpath[0] == "unlensed_Cl":
+            want_cl = True
             if qpath[1] == "tt" or qpath[1] == "te":
                 cl_spec.add("tCl")
             if qpath[1] == "te" or qpath[1] == "ee":
@@ -53,10 +55,23 @@ def classy_args_interpret(parser: YAMLParser, extra_args: dict = {}) -> dict:
             if "t" in cl_modes:
                 extra_args["l_max_tensors"] = max(extra_args["l_max_tensors"],
                                                   parser.modes(quantity).max())
-
-    if "s" not in cl_modes:
+        elif qpath[0] == "Pk":
+            cl_spec.add("mPk")
+            zval = parser.parameter_value("z")
+            if type(zval) == float:
+                z_max_pk = zval
+            elif type(zval) == tuple:
+                z_max_pk = max(zval)
+            extra_args["z_max_pk"] = max(extra_args.get("z_max_pk", 0.0), z_max_pk)
+            k_max_pk = parser.modes(quantity).max()
+            extra_args["P_k_max_1/Mpc"] = max(extra_args.get("P_k_max_1/Mpc", 1.0), k_max_pk)
+
+            if qpath[1] == "nonlin":
+                extra_args["non_linear"] = "yes"
+
+    if "s" not in cl_modes or not want_cl:
         extra_args.pop("l_max_scalars")
-    if "t" not in cl_modes:
+    if "t" not in cl_modes or not want_cl:
         extra_args.pop("l_max_tensors")
     extra_args["output"] = " ".join(list(cl_spec))
     extra_args["modes"] = " ".join(list(cl_modes))
@@ -85,6 +100,9 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
                 quantities: list = [], extra_args: dict = {}) -> bool:
 
     classy = state["module"]
+    
+    if "z" in args:
+        z_pk = args.pop("z")
 
     try:
         cosmo = classy.Class()
@@ -117,8 +135,16 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
             ell = cls["ell"]
             Cl = cls[spec]
             Dl = Cl * (ell * (ell + 1) / (2 * np.pi))
+        elif qpath[0] == "Pk":
+            k = state[quantity + ".modes"]
+            spec = qpath[1]
+            
+            if spec == "lin":
+                Pk = np.zeros_like(k)
+                for i in range(len(k)):
+                    Pk[i] = cosmo.pk_lin(k[i], z_pk)
 
-            state[quantity] = Dl[state[quantity + ".modes"]]
+            state[quantity] = Pk
 
     cosmo.struct_cleanup()
     cosmo.empty()

From 75d4adb9e013adfcb794a074c76f72c6aad29742 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Thu, 9 Apr 2026 11:49:46 +0100
Subject: [PATCH 16/23] Added non-linear matter power.

---
 cosmopower/spectra_classy.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index b60c8b0..2e6ca34 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -67,7 +67,7 @@ def classy_args_interpret(parser: YAMLParser, extra_args: dict = {}) -> dict:
             extra_args["P_k_max_1/Mpc"] = max(extra_args.get("P_k_max_1/Mpc", 1.0), k_max_pk)
 
             if qpath[1] == "nonlin":
-                extra_args["non_linear"] = "yes"
+                extra_args["non_linear"] = "hmcode"
 
     if "s" not in cl_modes or not want_cl:
         extra_args.pop("l_max_scalars")
@@ -76,6 +76,8 @@ def classy_args_interpret(parser: YAMLParser, extra_args: dict = {}) -> dict:
     extra_args["output"] = " ".join(list(cl_spec))
     extra_args["modes"] = " ".join(list(cl_modes))
 
+    print(extra_args)
+
     return extra_args
 
 
@@ -143,6 +145,10 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
                 Pk = np.zeros_like(k)
                 for i in range(len(k)):
                     Pk[i] = cosmo.pk_lin(k[i], z_pk)
+            if spec == "nonlin":
+                Pk = np.zeros_like(k)
+                for i in range(len(k)):
+                    Pk[i] = cosmo.pk(k[i], z_pk)
 
             state[quantity] = Pk
 

From d525c2a22ce2b97c7a3d4d8b578d04f0c0f60334 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Thu, 9 Apr 2026 12:01:02 +0100
Subject: [PATCH 17/23] Ignore class errors and continue generating.

---
 cosmopower/spectra_classy.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index 2e6ca34..8f51d6f 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -111,7 +111,10 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
         cosmo.set(state["params"] | args)
         cosmo.compute()
     except (classy.CosmoComputationError, classy.CosmoSevereError) as e:
-        raise e
+        print(str(e))
+        cosmo.struct_cleanup()
+        cosmo.empty()
+        return False
 
     state["derived"] = get_classy_derived(state["params"], cosmo,
                                           list(state["derived"].keys()))

From 070af987f799f957853906f809251ad9056d19d9 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Thu, 9 Apr 2026 12:06:15 +0100
Subject: [PATCH 18/23] Flake8

---
 cosmopower/spectra.py        |  2 +-
 cosmopower/spectra_classy.py | 15 +++++++++------
 2 files changed, 10 insertions(+), 7 deletions(-)

diff --git a/cosmopower/spectra.py b/cosmopower/spectra.py
index 1f6a154..ad06a88 100644
--- a/cosmopower/spectra.py
+++ b/cosmopower/spectra.py
@@ -265,7 +265,7 @@ def generate_spectra(args: list = None) -> None:
 
     for n in tbar:
         tbar.set_description(("" if MPI is None else f"[{rank}] ")
-                             + f"{accepted/max(n,1):.1%} success rate")
+                             + f"{accepted/max(n, 1):.1%} success rate")
 
         boltzmann_params = {k: training[k][n] for k in parser.boltzmann_inputs}
 
diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index 8f51d6f..39d700c 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -58,13 +58,16 @@ def classy_args_interpret(parser: YAMLParser, extra_args: dict = {}) -> dict:
         elif qpath[0] == "Pk":
             cl_spec.add("mPk")
             zval = parser.parameter_value("z")
-            if type(zval) == float:
+            if type(zval) is float:
                 z_max_pk = zval
-            elif type(zval) == tuple:
+            elif type(zval) is tuple:
                 z_max_pk = max(zval)
-            extra_args["z_max_pk"] = max(extra_args.get("z_max_pk", 0.0), z_max_pk)
+            extra_args["z_max_pk"] = max(extra_args.get("z_max_pk", 0.0),
+                                         z_max_pk)
             k_max_pk = parser.modes(quantity).max()
-            extra_args["P_k_max_1/Mpc"] = max(extra_args.get("P_k_max_1/Mpc", 1.0), k_max_pk)
+            extra_args["P_k_max_1/Mpc"] = max(extra_args.get("P_k_max_1/Mpc",
+                                                             1.0),
+                                              k_max_pk)
 
             if qpath[1] == "nonlin":
                 extra_args["non_linear"] = "hmcode"
@@ -102,7 +105,7 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
                 quantities: list = [], extra_args: dict = {}) -> bool:
 
     classy = state["module"]
-    
+
     if "z" in args:
         z_pk = args.pop("z")
 
@@ -143,7 +146,7 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
         elif qpath[0] == "Pk":
             k = state[quantity + ".modes"]
             spec = qpath[1]
-            
+
             if spec == "lin":
                 Pk = np.zeros_like(k)
                 for i in range(len(k)):

From 560d2a97e8aa13f960808a9444b9b0457d815174 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Thu, 9 Apr 2026 13:54:03 +0100
Subject: [PATCH 19/23] numpy inf deprecation.

---
 cosmopower/cosmopower_PCAplusNN.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/cosmopower/cosmopower_PCAplusNN.py b/cosmopower/cosmopower_PCAplusNN.py
index cabe407..e831282 100644
--- a/cosmopower/cosmopower_PCAplusNN.py
+++ b/cosmopower/cosmopower_PCAplusNN.py
@@ -934,8 +934,8 @@ def train(self, training_data: Sequence,
                 ).shuffle(n_training).batch(batch_sizes[i])
 
             # set up training loss
-            validation_loss = [np.infty]
-            best_loss = np.infty
+            validation_loss = [np.inf]
+            best_loss = np.inf
             early_stopping_counter = 0
 
             # loop over epochs

From 918ac73f1280b38c7ecc9f152a69321f674133f0 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Mon, 13 Apr 2026 14:25:48 +0100
Subject: [PATCH 20/23] Missing line.

---
 cosmopower/spectra_classy.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index 39d700c..5a956e0 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -143,6 +143,8 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
             ell = cls["ell"]
             Cl = cls[spec]
             Dl = Cl * (ell * (ell + 1) / (2 * np.pi))
+
+            state[quantity] = Dl[state[quantity + ".modes"]]
         elif qpath[0] == "Pk":
             k = state[quantity + ".modes"]
             spec = qpath[1]

From cc689169331a3f6e17529d97f9eb0225b4cd95bc Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Mon, 13 Apr 2026 16:14:04 +0100
Subject: [PATCH 21/23] Updated glob matching for >10 files.

---
 cosmopower/train.py    | 14 +++++++-------
 cosmopower/validate.py |  2 +-
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/cosmopower/train.py b/cosmopower/train.py
index 6c4127a..0647f1d 100644
--- a/cosmopower/train.py
+++ b/cosmopower/train.py
@@ -25,7 +25,7 @@ def train_network_NN(parser: YAMLParser, quantity: str, device: str = "",
     """
     filenames = glob.glob(os.path.join(parser.path, "spectra",
                                        quantity.replace("/", "_")
-                                       + ".[0-9].hdf5"))
+                                       + ".[0-9]*.hdf5"))
 
     if len(filenames) == 0:
         raise IOError(f"No files found to train quantity {quantity} with.")
@@ -55,7 +55,7 @@ def train_network_NN(parser: YAMLParser, quantity: str, device: str = "",
     else:
         filenames = glob.glob(os.path.join(parser.path, "spectra",
                                            quantity.replace("/", "_") +
-                                           ".validation.[0-9].hdf5"))
+                                           ".validation.[0-9]*.hdf5"))
         validation = [Dataset(parser, quantity, os.path.basename(filename))
                       for filename in filenames]
 
@@ -85,7 +85,7 @@ def train_network_PCAplusNN(parser: YAMLParser, quantity: str,
     """
     filenames = glob.glob(os.path.join(parser.path, "spectra",
                                        quantity.replace("/", "_")
-                                       + ".[0-9].hdf5"))
+                                       + ".[0-9]*.hdf5"))
 
     saved_filename = os.path.join(parser.path, "networks",
                                   parser.network_filename(quantity))
@@ -111,7 +111,7 @@ def train_network_PCAplusNN(parser: YAMLParser, quantity: str,
     else:
         filenames = glob.glob(os.path.join(parser.path, "spectra",
                                            quantity.replace("/", "_") +
-                                           ".validation.[0-9].hdf5"))
+                                           ".validation.[0-9]*.hdf5"))
         validation = [Dataset(parser, quantity, os.path.basename(filename))
                       for filename in filenames]
 
@@ -171,7 +171,7 @@ def show_training(args: Optional[list] = None) -> None:
     for n, quantity in enumerate(plot_quantities):
         filenames = glob.glob(os.path.join(parser.path, "spectra",
                                            quantity.replace("/", "_")
-                                           + ".[0-9].hdf5"))
+                                           + ".[0-9]*.hdf5"))
 
         if len(filenames) == 0:
             raise IOError(f"No files found for {quantity}.")
@@ -273,7 +273,7 @@ def show_validation(args: Optional[list] = None) -> None:
 
         filenames = glob.glob(os.path.join(parser.path, "spectra",
                                            quantity.replace("/", "_")
-                                           + ".validation.[0-9].hdf5"))
+                                           + ".validation.[0-9]*.hdf5"))
 
         if len(filenames) == 0:
             raise IOError(f"No files found for {quantity}.")
@@ -463,7 +463,7 @@ def show_validation(args: Optional[list] = None) -> None:
         ny = (len(parser.computed_parameters) // nx)
 
         filenames = glob.glob(os.path.join(parser.path, "spectra",
-                                           "derived.[0-9].hdf5"))
+                                           "derived.[0-9]*.hdf5"))
 
         if len(filenames) == 0:
             raise IOError("No files found for derived parameters.")
diff --git a/cosmopower/validate.py b/cosmopower/validate.py
index 567196e..923fbee 100644
--- a/cosmopower/validate.py
+++ b/cosmopower/validate.py
@@ -15,7 +15,7 @@ def find_files(parser: YAMLParser, quantity: str) -> np.ndarray:
     import glob
 
     fn = os.path.join(parser.path, "spectra",
-                      quantity.replace("/", "_") + ".*.hdf5")
+                      quantity.replace("/", "_") + ".[0-9]*.hdf5")
 
     return glob.glob(fn)
 

From 39414806aa296b55428bc72ae309a7511b419a1b Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Tue, 21 Apr 2026 14:39:59 +0100
Subject: [PATCH 22/23] Added DA, Hubble, sigma8

---
 cosmopower/spectra_camb.py   | 6 +++---
 cosmopower/spectra_classy.py | 9 +++++++++
 2 files changed, 12 insertions(+), 3 deletions(-)

diff --git a/cosmopower/spectra_camb.py b/cosmopower/spectra_camb.py
index 8262714..d0865b2 100644
--- a/cosmopower/spectra_camb.py
+++ b/cosmopower/spectra_camb.py
@@ -275,14 +275,14 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
                 state["params"], state["results"], parser.modes(qpath[0])
             )
         elif qpath[0] == "DA" or qpath[0] == "angular_diameter_distance":
-            # angular diameter distance
-            state[quantity] = get_camb_sigma8(
+            # angular diameter distance DA(z)
+            state[quantity] = get_camb_angular_diameter_distance(
                 state["params"], state["results"], parser.modes(qpath[0])
             )
         elif qpath[0] in ["Omega_b", "Omega_cdm"]:
             # background evolution of densities
             var = {"Omega_b": "baryon", "Omega_cdm": "cdm"}.get(qpath[0])
-            state[quantity] = get_camb_sigma8(
+            state[quantity] = get_camb_Omega(
                 state["params"], state["results"], parser.modes(qpath[0]), var
             )
         else:
diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index 5a956e0..0a0ec23 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -159,6 +159,15 @@ def get_spectra(parser: YAMLParser, state: dict, args: dict = {},
                     Pk[i] = cosmo.pk(k[i], z_pk)
 
             state[quantity] = Pk
+        elif qpath[0] == "DA" or qpath[0] == "angular_diameter_distance":
+            z = state[quantity + ".modes"]
+            state[quantity] = cosmo.angular_distance(z)
+        elif qpath[0] == "Hubble":
+            z = state[quantity + ".modes"]
+            state[quantity] = cosmo.Hubble(z)
+        elif qpath[0] == "sigma8":
+            z = state[quantity + ".modes"]
+            state[quantity] = cosmo.sigma(8, z)
 
     cosmo.struct_cleanup()
     cosmo.empty()

From 8a4ae3a86ed53dbafb8fd0a80d711064f7ba5905 Mon Sep 17 00:00:00 2001
From: Hidde Jense <JenseH@cardiff.ac.uk>
Date: Tue, 21 Apr 2026 15:20:33 +0100
Subject: [PATCH 23/23] Fixes

---
 cosmopower/spectra_classy.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/cosmopower/spectra_classy.py b/cosmopower/spectra_classy.py
index 0a0ec23..66fcc48 100644
--- a/cosmopower/spectra_classy.py
+++ b/cosmopower/spectra_classy.py
@@ -71,6 +71,15 @@ def classy_args_interpret(parser: YAMLParser, extra_args: dict = {}) -> dict:
 
             if qpath[1] == "nonlin":
                 extra_args["non_linear"] = "hmcode"
+        elif qpath[0] == "sigma8":
+            cl_spec.add("mPk")
+            z = parser.modes(quantity)
+            extra_args["z_max_pk"] = max(extra_args.get("z_max_pk", 0.0),
+                                         z.max())
+            k_max_pk = 1e-1
+            extra_args["P_k_max_1/Mpc"] = max(extra_args.get("P_k_max_1/Mpc",
+                                                             1.0),
+                                              k_max_pk)
 
     if "s" not in cl_modes or not want_cl:
         extra_args.pop("l_max_scalars")