Error handling for unallocated transfers; setuptools requirements

Author: cmbant

.github/workflows/tests.yml

@@ -13,9 +13,9 @@ jobs:
         - name: "Latest Python 3.x"
           os: ubuntu-latest
           python-version: 3.x
-        - name: "OS X Python 3.9"
+        - name: "OS X Python 3.10"
           os: macos-latest
-          python-version: 3.9
+          python-version: 3.10
         - name: "Windows Python 3.12"
           os: windows-latest
           python-version: 3.12

camb/baseconfig.py

@@ -318,7 +318,7 @@ class AllocatableArrayDouble(_AllocatableArray):
 AllocatableArrayDouble._set_allocatable_1D_array.argtypes = [POINTER(AllocatableArrayDouble), numpy_1d, POINTER(c_int)]
 
 
-def fortran_array(c_pointer, shape, dtype=np.float64, order='F', own_data=True):
+def fortran_array(c_pointer, shape, dtype: type =np.float64, order='F', own_data=True):
     if not hasattr(shape, '__len__'):
         shape = np.atleast_1d(shape)
     arr_size = np.prod(shape[:]) * np.dtype(dtype).itemsize
@@ -640,7 +640,7 @@ def __reduce__(self):
 
     def __getstate__(self) -> dict:
         state = {}
-        for field_name, field_type in self.get_all_fields():
+        for field_name, _ in self.get_all_fields():
             obj = getattr(self, field_name)
             if isinstance(obj, (ctypes.Array, np.ndarray, _ArrayOfAllocatable)):
                 state[field_name] = list(obj[:len(obj)])
@@ -685,6 +685,12 @@ class F2003Class(CAMB_Structure):
 
     __slots__ = ()
 
+    # _fields_ store ctypes fields, but here overridden  
+    _fields_: list[tuple[str, type] | tuple[str, type, str]]
+
+    # _methods_ for calling fortran functions name, arguments and optional return value  
+    _methods_: list[tuple[str, list, type] | tuple[str, list]]
+
     # pointer to fortran class; generated once per instance using _fortran_selfpointer_function then replaced with
     # the actual value
     fortran_self = _FortranSelf()

camb/results.py

@@ -90,8 +90,8 @@ class MatterTransferData:
 
     nq: int
     q: np.ndarray
-    sigma_8: np.ndarray
-    sigma2_vdelta_8: np.ndarray
+    sigma_8: np.ndarray | None
+    sigma2_vdelta_8: np.ndarray | None
     transfer_data: np.ndarray
 
     def transfer_z(self, name, z_index=0):
@@ -372,8 +372,7 @@ def _CMB_unit(self, CMB_unit):
 
     def _scale_cls(self, cls, CMB_unit=None, raw_cl=False, lens_potential=False):
         if raw_cl:
-            ls = np.arange(1, cls.shape[0])[..., np.newaxis]
-            ls = np.float64(ls * (ls + 1))
+            ls = np.arange(1, cls.shape[0], dtype=np.float64)[..., np.newaxis]
             if lens_potential:
                 cls[1:, 0] /= ls[:, 0] ** 2 / (2 * np.pi)
                 cls[1:, 1:] /= ls ** (3. / 2) / (2 * np.pi)
@@ -739,6 +738,8 @@ def get_linear_matter_power_spectrum(self, var1=None, var2=None, hubble_units=Tr
         num_k = c_int(0)
         CAMBdata_mattertransferks(byref(self), byref(num_k), np.array([]))
         nk = num_k.value
+        if not nk:
+            raise CAMBError('No k values available for matter power spectrum')
 
         ks = np.empty(nk, dtype=np.float64)
         CAMBdata_mattertransferks(byref(self), byref(num_k), ks)
@@ -754,9 +755,9 @@ def get_linear_matter_power_spectrum(self, var1=None, var2=None, hubble_units=Tr
         PK = np.empty((nz, nk))
         if nonlinear:
             CAMBdata_GetNonLinearMatterPower(byref(self), PK, byref(var1), byref(var2), byref(hubble_units))
-            config.check_global_error('get_[non]linear_matter_power_spectrum')
         else:
             CAMBdata_GetLinearMatterPower(byref(self), PK, byref(var1), byref(var2), byref(hubble_units))
+        config.check_global_error('get_[non]linear_matter_power_spectrum')
 
         z = self.Params.Transfer.PK_redshifts[:nz]
         z.reverse()
@@ -945,6 +946,10 @@ def get_matter_power_interpolator(self, nonlinear=True, var1=None, var2=None, hu
         class PKInterpolator(RectBivariateSpline):
             islog: bool
             logsign: int
+            kmin : float
+            kmax : float
+            zmin : float
+            zmax : float
 
             def P(self, z, kh, grid=None):
                 if grid is None:
@@ -1395,7 +1400,7 @@ def angular_diameter_distance2(self, z1, z2):
         else:
             return self.f_AngularDiameterDistance2(byref(c_double(z1)), byref(c_double(z2)))
 
-    def comoving_radial_distance(self, z, tol=1e-4):
+    def comoving_radial_distance(self, z: float, tol=1e-4):
         """
         Get comoving radial distance from us to redshift z in Mpc. This is efficient for arrays.
 

camb/tests/camb_test.py

@@ -11,7 +11,7 @@
     sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..', '..')))
     import camb
 from camb import model, correlations, bbn, dark_energy, initialpower
-from camb.baseconfig import CAMBParamRangeError, CAMBValueError
+from camb.baseconfig import CAMBParamRangeError, CAMBValueError, CAMBError
 
 fast = 'ci fast' in os.getenv("GITHUB_ACTIONS", '')
 
@@ -250,6 +250,15 @@ def testBackground(self):
         self.assertAlmostEqual(d, data.conformal_time_a1_a2(0, 1))
         self.assertAlmostEqual(d, sum(data.conformal_time_a1_a2([0, 0.5], [0.5, 1])))
 
+    def testErrors(self):
+        redshifts = np.logspace(-1, np.log10(1089))
+        pars = camb.set_params(H0=67.5, ombh2=0.022, omch2=0.122, As=2e-9, ns=0.95,
+                               redshifts=redshifts, kmax=0.1)
+
+        results = camb.get_background(pars)
+        with self.assertRaises(CAMBError):
+            results.get_matter_power_interpolator()
+
     def testEvolution(self):
         redshifts = [0.4, 31.5]
         pars = camb.set_params(H0=67.5, ombh2=0.022, omch2=0.122, As=2e-9, ns=0.95,

fortran/camb_python.f90

@@ -307,6 +307,11 @@ subroutine CAMBdata_GetMatterTransferks(State, nk, ks)
     real(dl) ks(nk)
 
     if (nk/=0) then
+        if (.not. allocated(State%MT%TransferData)) then
+            global_error_flag = error_allocation
+            global_error_message = 'TransferData not allocated: make sure transfer functions computed'
+            return
+        end if
         ks =  State%MT%TransferData(Transfer_kh,:,1) * (State%CP%H0/100)
     end if
     nk = State%MT%num_q_trans

fortran/config.f90

@@ -60,6 +60,7 @@ module config
     integer, parameter :: error_darkenergy=6
     integer, parameter :: error_ini=7
     integer, parameter :: error_nonlinear=8
+    integer, parameter :: error_allocation=9
 
     contains
 

pyproject.toml

@@ -1,5 +1,5 @@
 [build-system]
-requires = ["setuptools", "packaging", "wheel"]
+requires = ["setuptools>=77", "packaging>=24", "wheel"]
 build-backend = "setuptools.build_meta"
 
 [project]
@@ -13,21 +13,20 @@ readme = "docs/README_pypi.rst"
 license = "LicenseRef-LGPLv3-arxiv"
 license-files = ["LICENCE.txt"]
 dynamic = ["version"]
-requires-python = ">=3.9.0"
+requires-python = ">=3.10.0"
 classifiers = [
     "Development Status :: 5 - Production/Stable",
     "Operating System :: OS Independent",
     "Environment :: Console",
     "Intended Audience :: Science/Research",
     "Topic :: Scientific/Engineering :: Astronomy",
-    "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
     "Programming Language :: Python :: 3.13"
 ]
 dependencies = [
-    "numpy", "scipy>=1.0", "sympy>=1.0", "packaging"
+    "numpy", "scipy>=1.0", "sympy>=1.0", "packaging>=24"
 ]
 
 [project.optional-dependencies]

requirements.txt

@@ -1,4 +1,4 @@
 numpy
 sympy>=1.0
 scipy>=1.0
-packaging
+packaging>=24