more consistent error propagation from Fortran; effective n_s for spl…

…ined initial power

camb/_config.py

@@ -1,4 +1,5 @@
-from .baseconfig import import_property, c_int, c_bool, c_double
+from .baseconfig import import_property, CAMBError
+from ctypes import c_char, c_int, c_bool, c_double
 
 lensing_method_curv_corr = 1
 lensing_method_flat_corr = 2
@@ -30,11 +31,31 @@ class _config(object):
 
     transfer_power_var = import_property(c_int, "transfer", "transfer_power_var")
 
+    _global_error_message = import_property(c_char * 1024, "config", "global_error_message")
+
+    def global_error_message(self):
+        return bytearray(self._global_error_message).decode('ascii').strip()
+
+    def check_global_error(self, reference=''):
+        code = self.global_error_flag
+        if code:
+            err = config.global_error_message()
+            self.global_error_flag = 0
+            if reference:
+                reference = 'Error in Fortran called from %s:\n' % reference
+            else:
+                reference = ''
+            if err:
+                raise CAMBError(reference + '%s' % err)
+            else:
+                raise CAMBError(reference + 'Error code: %s' % code)
+
     def __repr__(self):
         s = ''
         for x in dir(self):
             if x[0] != '_':
                 s += '%s = %s\n' % (x, getattr(self, x))
         return s
 
+
 config = _config()

camb/camb.py

@@ -217,12 +217,13 @@ def validate_ini_file(filename):
     import subprocess
     import sys
     try:
+        err = ''
         command = '"%s" "%s" "%s" --validate' % (
             sys.executable, os.path.join(os.path.dirname(__file__), '_command_line.py'), filename)
         subprocess.check_output(command, stderr=subprocess.STDOUT, shell=True)
     except subprocess.CalledProcessError as E:
         err = E.output.decode().replace('ERROR STOP', '').strip()
-        raise CAMBValueError(err + ' (%s)' % filename)
+    if err: raise CAMBValueError(err + ' (%s)' % filename)
     return True
 
 
@@ -239,12 +240,11 @@ def run_ini(ini_filename, no_validate=False):
         raise CAMBValueError('File not found: %s' % ini_filename)
     if not no_validate: validate_ini_file(ini_filename)
     runIni = camblib.__camb_MOD_camb_runinifile
-    runIni.argtypes = [ctypes.c_char_p, POINTER(ctypes.c_long), ctypes.c_char_p, POINTER(ctypes.c_long)]
+    runIni.argtypes = [ctypes.c_char_p, POINTER(ctypes.c_long)]
     runIni.restype = c_bool
     s = ctypes.create_string_buffer(six.b(ini_filename))
-    err = ctypes.create_string_buffer(1025)
-    if not runIni(s, ctypes.c_long(len(ini_filename)), err, ctypes.c_long(1024)):
-        raise CAMBError(err.value.decode('ascii'))
+    if not runIni(s, ctypes.c_long(len(ini_filename))):
+        config.check_global_error('run_ini')
 
 
 def read_ini(ini_filename, no_validate=False):
@@ -259,13 +259,11 @@ def read_ini(ini_filename, no_validate=False):
     if not no_validate: validate_ini_file(ini_filename)
     cp = model.CAMBparams()
     readIni = camblib.__camb_MOD_camb_readparamfile
-    readIni.argtypes = [POINTER(CAMBparams), ctypes.c_char_p, POINTER(ctypes.c_long), ctypes.c_char_p,
-                        POINTER(ctypes.c_long)]
+    readIni.argtypes = [POINTER(CAMBparams), ctypes.c_char_p, POINTER(ctypes.c_long)]
     readIni.restype = ctypes.c_bool
     s = ctypes.create_string_buffer(six.b(ini_filename))
-    err = ctypes.create_string_buffer(1025)
-    if not readIni(cp, s, ctypes.c_long(len(ini_filename)), err, ctypes.c_long(1024)):
-        raise CAMBValueError(err.value.decode('ascii'))
+    if not readIni(cp, s, ctypes.c_long(len(ini_filename))):
+        config.check_global_error('read_ini')
     return cp
 
 

camb/initialpower.py

@@ -26,6 +26,9 @@ class SplinedInitialPower(InitialPower):
     """
     _fortran_class_name_ = 'TSplinedInitialPower'
 
+    _fields_ = [
+        ('effective_ns_for_nonlinear', c_double, "Effective n_s to use for approximate non-linear correction models")]
+
     _methods_ = [('HasTensors', [], c_int),
                  ('SetScalarTable', [POINTER(c_int), numpy_1d, numpy_1d]),
                  ('SetTensorTable', [POINTER(c_int), numpy_1d, numpy_1d]),

camb/model.py

@@ -86,7 +86,8 @@ class AccuracyParams(CAMB_Structure):
     _fields_ = [
         ("AccuracyBoost", c_double,
          "general accuracy setting effecting everything related to step sizes etc. (including separate settings below except the next two)"),
-        ("lSampleBoost", c_double, "accuracy for sampling in ell for interpolation for the C_l (if >=50, all ell are calculated)"),
+        ("lSampleBoost", c_double,
+         "accuracy for sampling in ell for interpolation for the C_l (if >=50, all ell are calculated)"),
         ("lAccuracyBoost", c_double, "Boosts number of multipoles integrated in Boltzman heirarchy"),
         ("AccuratePolarization", c_bool, "Do you care about the accuracy of the polarization Cls?"),
         ("AccurateBB", c_bool, "Do you care about BB accuracy (e.g. in lensing)"),
@@ -229,7 +230,8 @@ class in model.py to add the new parameter in the corresponding location of the
         ("transfer_21cm_cl", c_bool, "Get 21cm C_L at a given fixed redshift"),
         ("Log_lvalues", c_bool, "Use log spacing for sampling in L"),
         (
-        "use_cl_spline_template", c_bool, "When interpolating use a fiducial spectrum shape to define ratio to spline"),
+            "use_cl_spline_template", c_bool,
+            "When interpolating use a fiducial spectrum shape to define ratio to spline"),
 
         ("SourceWindows", AllocatableObjectArray(SourceWindow)),
         ("CustomSources", CustomSources)
@@ -274,7 +276,8 @@ def set_accuracy(self, AccuracyBoost=1., lSampleBoost=1., lAccuracyBoost=1., DoL
         self.DoLateRadTruncation = DoLateRadTruncation
         return self
 
-    def set_initial_power_function(self, P_scalar, P_tensor=None, kmin=1e-6, kmax=100., N_min=200, rtol=5e-5, args=()):
+    def set_initial_power_function(self, P_scalar, P_tensor=None, kmin=1e-6, kmax=100., N_min=200, rtol=5e-5,
+                                   effective_ns_for_nonlinear=None, args=()):
         r"""
         Set the initial power spectrum from a function P_scalar(k, \*args), and optionally also the tensor spectrum.
         The function is called to make a pre-computed array which is then interpolated inside CAMB. The sampling in k
@@ -286,6 +289,7 @@ def set_initial_power_function(self, P_scalar, P_tensor=None, kmin=1e-6, kmax=10
         :param kmax: maximum wavenumber to compute
         :param N_min: minimum number of spline points for the pre-computation
         :param rtol: relative tolerance for deciding how many points are enough
+        :param effective_ns_for_nonlinear: an effective n_s for use with approximate non-linear corrections
         :param args: optional list of arguments passed to P_scalar (and P_tensor)
         :return: self
         """
@@ -307,20 +311,23 @@ def set_initial_power_function(self, P_scalar, P_tensor=None, kmin=1e-6, kmax=10
             PK_test = PK
             ktest = ks
         PK_t = None if P_tensor is None else P_tensor(ks, *args)
-        self.set_initial_power_table(ks, PK, PK_t)
+        self.set_initial_power_table(ks, PK, PK_t, effective_ns_for_nonlinear)
         return self
 
-    def set_initial_power_table(self, k, pk=None, pk_tensor=None):
+    def set_initial_power_table(self, k, pk=None, pk_tensor=None, effective_ns_for_nonlinear=None):
         """
         Set a general intial power spectrum from tabulated values. It's up to you to ensure the sampling
         of the k values is high enough that it can be interpolated accurately.
 
         :param k: array of k values (Mpc^{-1})
         :param pk: array of primordial curvature perturbation power spectrum values P(k_i)
         :param pk_tensor: array of tensor spectrum values
+        :param effective_ns_for_nonlinear: an effective n_s for use with approximate non-linear corrections
         """
         self.InitPower = SplinedInitialPower()
         initpower = self.InitPower
+        if effective_ns_for_nonlinear is not None:
+            initpower.effective_ns_for_nonlinear = effective_ns_for_nonlinear
         if pk is None:
             pk = np.asarray([])
         elif len(k) != len(pk):

camb/results.py

@@ -244,8 +244,7 @@ def calc_background_no_thermo(self, params):
         """
         self._check_params(params)
         self.f_SetParams(byref(params), None, None, None)
-        if config.global_error_flag != 0:
-            raise CAMBError('Error %s settings parameters' % config.global_error_flag)
+        config.check_global_error('calc_background_no_thermo')
 
     def calc_background(self, params):
         """
@@ -254,9 +253,8 @@ def calc_background(self, params):
         :param params:  :class:`~.model.CAMBparams` instance to use
         """
         self._check_params(params)
-        res = CAMBdata_CalcBackgroundTheory(byref(self), byref(params))
-        if res:
-            raise CAMBError('Error %s in calc_background' % res)
+        if CAMBdata_CalcBackgroundTheory(byref(self), byref(params)):
+            config.check_global_error('calc_background')
 
     def calc_transfers(self, params, only_transfers=True):
         """
@@ -268,8 +266,8 @@ def calc_transfers(self, params, only_transfers=True):
         """
         self._check_params(params)
         if not only_transfers: self._check_powers(params)
-        res = CAMBdata_gettransfers(byref(self), byref(params), byref(c_int(1 if only_transfers else 0)))
-        return res
+        if CAMBdata_gettransfers(byref(self), byref(params), byref(c_int(1 if only_transfers else 0))):
+            config.check_global_error('calc_transfer')
 
     def _check_powers(self, params=None):
         if params is None: params = self.Params
@@ -286,12 +284,11 @@ def calc_power_spectra(self, params=None):
 
         """
         if params is not None:
-            result = self.calc_transfers(params, only_transfers=False)
-            if result != 0:
-                raise CAMBError('Error getting transfer functions: %u' % result)
+            self.calc_transfers(params, only_transfers=False)
         else:
             self._check_powers()
             CAMBdata_transferstopowers(byref(self))
+            config.check_global_error()
 
     def power_spectra_from_transfer(self, initial_power_params, silent=False):
         """
@@ -307,12 +304,14 @@ def power_spectra_from_transfer(self, initial_power_params, silent=False):
         :param silent: suppress warnings about non-linear corrections not being recalculated
         """
         if not silent and self.Params.NonLinear in [model.NonLinear_lens, model.NonLinear_both] and \
-                self.Params.WantScalars and self.Params.WantCls:
+                self.Params.WantScalars and self.Params.WantCls and not getattr(self, '_suppress_power_warn', False):
             logging.warning(
                 'power_spectra_from_transfer with non-linear lensing does not recalculate the non-linear correction')
+            self._suppress_power_warn = True
         self.Params.set_initial_power(initial_power_params)
         self._check_powers()
         CAMBdata_transferstopowers(byref(self))
+        config.check_global_error()
 
     def _CMB_unit(self, CMB_unit):
         if isinstance(CMB_unit, six.string_types):
@@ -495,10 +494,9 @@ def get_time_evolution(self, q, eta, vars=model.evolve_names, lAccuracyBoost=4,
             nvars = num_standard_names + ncustom
             outputs = np.empty((k.shape[0], times.shape[0], nvars))
             if CAMB_TimeEvolution(byref(self), byref(c_int(k.shape[0])), k, byref(c_int(times.shape[0])),
-                                  times[indices],
-                                  byref(c_int(nvars)), outputs,
+                                  times[indices], byref(c_int(nvars)), outputs,
                                   byref(c_int(ncustom)), byref(custom_source_func)):
-                raise CAMBError('Error in evolution')
+                config.check_global_error('get_time_evolution')
             i_rev = np.zeros(times.shape, dtype=int)
             i_rev[indices] = np.arange(times.shape[0])
             outputs = outputs[:, i_rev, :]

camb_tests/camb_test.py

@@ -494,6 +494,21 @@ def testSources(self):
         ls = np.arange(0, PP.shape[0])
         self.assertTrue(np.allclose(PP / 4 * (ls * (ls + 1)), cls['W1xW1'], rtol=1e-3))
         self.assertTrue(np.allclose(PP / 2 * np.sqrt(ls * (ls + 1)), cls['PxW1'], rtol=1e-3))
+        # test something sharp with redshift distortions (tricky..)
+        from scipy import signal
+        zs = np.arange(1.9689, 2.1057, (2.1057 - 1.9689) / 2000)
+        W = signal.tukey(len(zs), alpha=0.1)
+        pars = camb.CAMBparams()
+        pars.set_cosmology(H0=67.5, ombh2=0.022, omch2=0.122)
+        pars.InitPower.set_params(As=2e-9, ns=0.965)
+        pars.set_for_lmax(4000)
+        pars.SourceWindows = [SplinedSourceWindow(z=zs, W=W, source_type='counts')]
+        pars.SourceTerms.counts_redshift = True
+        results = camb.get_results(pars)
+        cls = results.get_source_cls_dict()
+        self.assertAlmostEqual(np.sum(cls['PxW1'][10:3000:20]), 0.00020001, places=5)
+        self.assertAlmostEqual(np.sum(cls['W1xW1'][10:3000:20]), 2.26348, places=3)
+        self.assertAlmostEqual(np.sum(cls['W1xW1'][10]), 0.0001097, places=6)
 
     def testSymbolic(self):
         if fast: return

fortran/InitialPower.f90

@@ -65,6 +65,7 @@ module InitialPower
     end Type TInitialPowerLaw
 
     Type, extends(TInitialPower) :: TSplinedInitialPower
+        real(dl) :: effective_ns_for_nonlinear = -1._dl !used for halofit
         class(TSpline1D), allocatable :: Pscalar, Ptensor
     contains
     procedure :: SetScalarTable => TSplinedInitialPower_SetScalarTable
@@ -74,6 +75,7 @@ module InitialPower
     procedure :: ScalarPower => TSplinedInitialPower_ScalarPower
     procedure :: TensorPower => TSplinedInitialPower_TensorPower
     procedure :: HasTensors => TSplinedInitialPower_HasTensors
+    procedure :: Effective_ns => TSplinedInitialPower_Effective_ns
     procedure, nopass :: PythonClass => TSplinedInitialPower_PythonClass
     procedure, nopass :: SelfPointer => TSplinedInitialPower_SelfPointer
     end Type TSplinedInitialPower
@@ -340,6 +342,18 @@ subroutine TSplinedInitialPower_SetTensorLogRegular(this, kmin, kmax, n, PK)
 
     end subroutine TSplinedInitialPower_SetTensorLogRegular
 
+    function TSplinedInitialPower_Effective_ns(this)
+    use config
+    class(TSplinedInitialPower) :: this
+    real(dl) :: TSplinedInitialPower_Effective_ns
+
+    if (this%effective_ns_for_nonlinear==-1._dl) then
+        call GlobalError('TSplinedInitialPower: effective_ns_for_nonlinear not set',error_inital_power)
+    else
+        TSplinedInitialPower_Effective_ns = this%effective_ns_for_nonlinear
+    end if
+    end function TSplinedInitialPower_Effective_ns
+
 
 
     end module InitialPower

fortran/camb.cbp

@@ -34,52 +34,36 @@
 				</Linker>
 			</Target>
 		</Build>
+		<Unit filename="DarkAge21cm.f90" />
+		<Unit filename="DarkEnergyFluid.f90" />
+		<Unit filename="DarkEnergyInterface.f90" />
+		<Unit filename="DarkEnergyPPF.f90" />
+		<Unit filename="InitialPower.f90" />
 		<Unit filename="Makefile">
 			<Option virtualFolder="makefile/" />
 		</Unit>
 		<Unit filename="Makefile_main">
 			<Option virtualFolder="makefile/" />
 		</Unit>
-		<Unit filename="bessels.f90">
-			<Option weight="3" />
-		</Unit>
-		<Unit filename="camb.f90">
-			<Option weight="6" />
-		</Unit>
-		<Unit filename="cmbmain.f90">
-			<Option weight="5" />
-		</Unit>
-		<Unit filename="constants.f90">
-			<Option weight="0" />
-		</Unit>
-		<Unit filename="equations.f90">
-			<Option weight="3" />
-		</Unit>
-		<Unit filename="halofit_ppf.f90">
-			<Option weight="4" />
-		</Unit>
-		<Unit filename="inidriver.F90" />
-		<Unit filename="inifile.f90">
-			<Option weight="0" />
-		</Unit>
-		<Unit filename="lensing.f90">
-			<Option weight="3" />
-		</Unit>
-		<Unit filename="modules.f90">
-			<Option weight="2" />
-		</Unit>
-		<Unit filename="power_tilt.f90">
-			<Option weight="1" />
-		</Unit>
-		<Unit filename="recfast.f90">
-			<Option weight="1" />
-		</Unit>
-		<Unit filename="reionization.f90">
-			<Option weight="1" />
-		</Unit>
-		<Unit filename="subroutines.f90">
-			<Option weight="1" />
-		</Unit>
+		<Unit filename="MathUtils.f90" />
+		<Unit filename="SecondOrderPK.f90" />
+		<Unit filename="SeparableBispectrum.f90" />
+		<Unit filename="SourceWindows.f90" />
+		<Unit filename="bessels.f90" />
+		<Unit filename="camb.f90" />
+		<Unit filename="camb_python.f90" />
+		<Unit filename="classes.f90" />
+		<Unit filename="cmbmain.f90" />
+		<Unit filename="config.f90" />
+		<Unit filename="constants.f90" />
+		<Unit filename="equations.f90" />
+		<Unit filename="halofit.f90" />
+		<Unit filename="lensing.f90" />
+		<Unit filename="massive_neutrinos.f90" />
+		<Unit filename="recfast.f90" />
+		<Unit filename="reionization.f90" />
+		<Unit filename="results.f90" />
+		<Unit filename="subroutines.f90" />
 		<Extensions>
 			<envvars />
 			<code_completion />