Add ExternalNonLinearRatio NonLinearModel subclass

camb/_compilers.py

@@ -4,8 +4,6 @@
 import struct
 import subprocess
 
-from camb.baseconfig import CAMBError
-
 is_windows = platform.system() == "Windows"
 
 is_32_bit = struct.calcsize("P") == 4
@@ -86,7 +84,7 @@ def check_gfortran(version=gfortran_min, msg=False, retry=False):
         version_str = str(subprocess.check_output("gfortran -dumpmachine", shell=True, env=compiler_environ))
         ok = gfortran_bits in version_str
     if not ok and msg:
-        raise CAMBError(
+        raise RuntimeError(
             f"You need ifort or gfortran {version} or higher to compile (found: {gfortran_version}).\n"
             "See https://camb.readthedocs.io/en/latest/fortran_compilers.htm\n"
             "or install from pypi using pip ('pip install camb') to just use pre-built binary wheels."

camb/baseconfig.py

@@ -200,7 +200,7 @@ def AllocatableObject(cls=None):
 
 
 class _AllocatableArray(FortranAllocatable):  # member corresponding to allocatable :: d(:) member in fortran
-    _fields_ = ("allocatable", ctypes.c_void_p * (_f_allocatable_array_size // ctypes.sizeof(ctypes.c_void_p)))
+    _fields_ = (("allocatable", ctypes.c_void_p * (_f_allocatable_array_size // ctypes.sizeof(ctypes.c_void_p))),)
 
     def get_allocatable(self):
         size = self._get_allocatable_1D_array(byref(self), byref(_reuse_pointer))

camb/initialpower.py

@@ -190,8 +190,11 @@ def set_params(
             # set from inflationary consistency
             if ntrun:
                 raise CAMBError("ntrun set but using inflation consistency (nt=None)")
-            if self.tensor_parameterization != tensor_param_rpivot:
-                raise CAMBError("tensor parameterization not tensor_param_rpivot with inflation consistency")
+            if self.tensor_parameterization != "tensor_param_rpivot":
+                raise CAMBError(
+                    "tensor parameterization not tensor_param_rpivot with inflation consistency: "
+                    f"{self.tensor_parameterization}"
+                )
             self.nt = -r / 8.0 * (2.0 - ns - r / 8.0)
             self.ntrun = r / 8.0 * (r / 8.0 + ns - 1)
         else:

camb/nonlinear.py

@@ -1,6 +1,9 @@
-from ctypes import c_double, c_int
+from ctypes import POINTER, byref, c_double, c_int
 
-from .baseconfig import F2003Class, fortran_class
+import numpy as np
+from numpy.ctypeslib import ndpointer
+
+from .baseconfig import F2003Class, fortran_class, numpy_1d
 
 
 class NonLinearModel(F2003Class):
@@ -107,3 +110,57 @@ class SecondOrderPK(NonLinearModel):
 
     def set_params(self):
         pass
+
+
+@fortran_class
+class ExternalNonLinearRatio(NonLinearModel):
+    """
+    Non-linear model that applies a user-supplied ratio sqrt(P_NL/P_L)
+    from an external source (e.g. CCL, axionHMcode).
+
+    Use :meth:`set_ratio` to provide the ratio grid, then assign this
+    as ``params.NonLinearModel`` before calling :func:`camb.get_results`.
+    Can be called after only time transfers computed to then get
+    lensed C_l with a consistent non-linear model. If linear P_k sampling
+    points in k_h are used, no interpolation from provided grid is required.
+    """
+
+    _fortran_class_module_ = "ExternalNonLinearRatio"
+    _fortran_class_name_ = "TExternalNonLinearRatio"
+
+    _methods_ = (
+        (
+            "SetRatio",
+            [
+                POINTER(c_int),
+                POINTER(c_int),
+                numpy_1d,
+                numpy_1d,
+                ndpointer(c_double, flags="F_CONTIGUOUS", ndim=2),
+            ],
+        ),
+        ("ClearRatio", []),
+    )
+
+    def set_ratio(self, k_h, z, ratio):
+        """
+        Set the non-linear ratio grid sqrt(P_NL/P_L).
+
+        :param k_h: 1D array of k values in h/Mpc units (ascending)
+        :param z: 1D array of redshift values (ascending)
+        :param ratio: 2D array of sqrt(P_NL/P_L), shape (len(z), len(k_h)),
+                      matching the convention of CAMB's get_matter_power_spectrum
+        """
+        k_h = np.ascontiguousarray(k_h, dtype=np.float64)
+        z = np.ascontiguousarray(z, dtype=np.float64)
+        if ratio.shape != (len(z), len(k_h)):
+            raise ValueError(f"ratio shape {ratio.shape} must be (len(z), len(k_h)) = ({len(z)}, {len(k_h)})")
+        # Fortran expects (nk, nz) column-major; C-order (nz, nk) has the same memory layout
+        ratio_f = np.asfortranarray(ratio.T, dtype=np.float64)
+        self.f_SetRatio(byref(c_int(len(k_h))), byref(c_int(len(z))), k_h, z, ratio_f)
+
+    def clear_ratio(self):
+        """
+        Clear the stored ratio, releasing interpolation data.
+        """
+        self.f_ClearRatio()

fortran/ExternalNonLinearRatio.f90

@@ -0,0 +1,100 @@
+
+    module ExternalNonLinearRatio
+    ! Non-linear model that applies a user-supplied ratio sqrt(P_NL/P_L)
+    ! from an external source (e.g. CCL, axionHMcode).
+    ! The ratio is stored as a 2D interpolation grid over (k/h, z).
+    use results
+    use transfer
+    use Interpolation, only : TInterpGrid2D
+    implicit none
+    private
+
+    type, extends(TNonLinearModel) :: TExternalNonLinearRatio
+        Type(TInterpGrid2D) :: Ratio
+        logical :: ratio_set = .false.
+    contains
+    procedure :: GetNonLinRatios => TExternalNonLinearRatio_GetNonLinRatios
+    procedure :: GetNonLinRatios_All => TExternalNonLinearRatio_GetNonLinRatios_All
+    procedure :: SetRatio => TExternalNonLinearRatio_SetRatio
+    procedure :: ClearRatio => TExternalNonLinearRatio_ClearRatio
+    procedure, nopass :: SelfPointer => TExternalNonLinearRatio_SelfPointer
+    end type TExternalNonLinearRatio
+
+    public TExternalNonLinearRatio
+    contains
+
+    subroutine TExternalNonLinearRatio_SelfPointer(cptr, P)
+    use iso_c_binding
+    Type(c_ptr) :: cptr
+    Type(TExternalNonLinearRatio), pointer :: PType
+    class(TPythonInterfacedClass), pointer :: P
+
+    call c_f_pointer(cptr, PType)
+    P => PType
+
+    end subroutine TExternalNonLinearRatio_SelfPointer
+
+    subroutine TExternalNonLinearRatio_SetRatio(this, nk, nz, k_h_arr, z_arr, ratio_arr)
+    ! Set the non-linear ratio grid sqrt(P_NL/P_L)
+    class(TExternalNonLinearRatio) :: this
+    integer, intent(in) :: nk, nz
+    real(dl), intent(in) :: k_h_arr(nk)
+    real(dl), intent(in) :: z_arr(nz)
+    real(dl), intent(in) :: ratio_arr(nk, nz)
+
+    call this%Ratio%Init(k_h_arr, z_arr, ratio_arr)
+    this%ratio_set = .true.
+
+    end subroutine TExternalNonLinearRatio_SetRatio
+
+    subroutine TExternalNonLinearRatio_ClearRatio(this)
+    ! Clear the stored ratio, releasing interpolation data
+    class(TExternalNonLinearRatio) :: this
+
+    call this%Ratio%Clear()
+    this%ratio_set = .false.
+
+    end subroutine TExternalNonLinearRatio_ClearRatio
+
+    subroutine TExternalNonLinearRatio_GetNonLinRatios(this, State, CAMB_Pk)
+    ! Fill CAMB_Pk%nonlin_ratio from the stored interpolation grid.
+    ! Clamps k and z to the grid range to prevent unreliable polynomial extrapolation.
+    class(TExternalNonLinearRatio) :: this
+    class(TCAMBdata) :: State
+    type(MatterPowerData), target :: CAMB_Pk
+    integer :: ik, iz
+    real(dl) :: kh, z, kh_clamped, z_clamped
+    real(dl) :: k_min, k_max, z_min, z_max
+
+    if (.not. this%ratio_set) &
+        error stop 'ExternalNonLinearRatio: ratio not set. Call SetRatio first.'
+
+    ! Get bounds of the ratio grid to clamp values
+    k_min = this%Ratio%x(1)
+    k_max = this%Ratio%x(this%Ratio%nx)
+    z_min = this%Ratio%y(1)
+    z_max = this%Ratio%y(this%Ratio%ny)
+
+    CAMB_Pk%nonlin_ratio = 1.0_dl
+    do iz = 1, CAMB_Pk%num_z
+        z = CAMB_Pk%redshifts(iz)
+        z_clamped = max(z_min, min(z_max, z))
+        do ik = 1, CAMB_Pk%num_k
+            kh = exp(CAMB_Pk%log_kh(ik))
+            kh_clamped = max(k_min, min(k_max, kh))
+            CAMB_Pk%nonlin_ratio(ik, iz) = this%Ratio%Value(kh_clamped, z_clamped)
+        end do
+    end do
+
+    end subroutine TExternalNonLinearRatio_GetNonLinRatios
+
+    subroutine TExternalNonLinearRatio_GetNonLinRatios_All(this, State, CAMB_Pk)
+    class(TExternalNonLinearRatio) :: this
+    class(TCAMBdata) :: State
+    type(MatterPowerData), target :: CAMB_Pk
+
+    error stop 'ExternalNonLinearRatio: GetNonLinRatios_All not supported (no velocity corrections)'
+
+    end subroutine TExternalNonLinearRatio_GetNonLinRatios_All
+
+    end module ExternalNonLinearRatio

fortran/Makefile_main

@@ -5,7 +5,7 @@
 POWERSPECTRUM_FILES ?= InitialPower
 REIONIZATION_FILES ?= reionization
 RECOMBINATION_FILES ?= recfast
-NONLINEAR_FILES     ?= halofit SecondOrderPK
+NONLINEAR_FILES     ?= halofit SecondOrderPK ExternalNonLinearRatio
 DARKENERGY_FILES  ?= DarkEnergyFluid DarkEnergyPPF PowellMinimize DarkEnergyQuintessence
 
 BISPECTRUM ?= SeparableBispectrum