Version used in the KiDS-1000 beyond flat LCDM paper

cosmosis/cosmosis_reaction_module.py

@@ -1,6 +1,4 @@
-import os
 import numpy as np
-import scipy.interpolate
 
 from cosmosis.datablock import option_section, names
 
@@ -12,8 +10,14 @@ def setup(options):
     config = {}
 
     config["module"] = pyreact.ReACT()
+    config["mode"] = options.get_string(option_section, "mode", "f(R)").lower()
+    if not config["mode"] in ("f(r)", "dgp", "gr"):
+        raise ValueError(f"ReACT mode {config['mode']} not supported.")
+
+    config["log10_fR0"] = options.get_bool(option_section, "log10_fR0", True)
     config["verbose"] = options.get_int(option_section, "verbose", 1)
     config["massloop"] = options.get_int(option_section, "massloop", 30)
+    config["z_max"] = options.get_double(option_section, "z_max", 2.5)
     config["reaction_output_section"] = options.get_string(option_section, "reaction_output_section", "reaction")
     config["linear_matter_power_output_section"] = options.get_string(option_section, "linear_matter_power_output_section", names.matter_power_lin)
 
@@ -25,23 +29,50 @@ def execute(block, config):
     omega_b = block[names.cosmological_parameters, "omega_b"]
     sigma_8 = block[names.cosmological_parameters, "sigma_8"]
     n_s = block[names.cosmological_parameters, "n_s"]
-    mg1 = block[names.cosmological_parameters, "mg1"]
 
-    Pk = block[names.matter_power_lin, "p_k"][0]
+    fR0 = None
+    Omega_rc = None
+
+    if config["mode"] == "f(r)":
+        if config["log10_fR0"]:
+            fR0 = 10**block[names.cosmological_parameters, "log10_fR0"]
+        else:
+            fR0 = block[names.cosmological_parameters, "fR0"]
+    elif config["mode"] == "dgp":
+        Omega_rc = block[names.cosmological_parameters, "Omega_rc"]
+
+    Pk = block[names.matter_power_lin, "p_k"]
     k_h = block[names.matter_power_lin, "k_h"]
     z = block[names.matter_power_lin, "z"]
 
-    # np.savetxt("pofk.txt", Pk)
-    # np.savetxt("k.txt", k_h)
-    # np.savetxt("z.txt", z)
+    z_react = z[z < config["z_max"]]
 
-    reaction, pofk_lin = config["module"].compute_reaction(
-                                h, n_s, omega_m, omega_b, sigma_8, mg1, 
-                                z, k_h, Pk, is_transfer=False, mass_loop=config["massloop"],
+    try:
+        reaction, pofk_lin, sigma_8_MG = config["module"].compute_reaction(
+                                h, n_s, omega_m, omega_b, sigma_8, 
+                                z_react, k_h, Pk[0], 
+                                model=config["mode"], fR0=fR0, Omega_rc=Omega_rc, 
+                                is_transfer=False, mass_loop=config["massloop"],
                                 verbose=config["verbose"])
+    except:
+        return 1
+
+    # Replace linear power spectrum below z_max with MG linear power spectrum
+    Pk[z < config["z_max"]] = pofk_lin
+    # Pad the reaction with 1 above z_max
+    reaction = np.concatenate((reaction, np.ones((np.count_nonzero(z >= config["z_max"]), len(k_h)))), axis=0)
 
     block.put_grid(config["reaction_output_section"], "z", z, "k_h", k_h, "reaction", reaction)
-    block.replace_grid(config["linear_matter_power_output_section"], "z", z, "k_h", k_h, "p_k", pofk_lin)
+    block.replace_grid(config["linear_matter_power_output_section"], "z", z, "k_h", k_h, "p_k", Pk)
+
+    # Replace sigma8 with MG sigma8
+    S8_LCDM = block[names.cosmological_parameters, "S_8"]
+
+    block[names.cosmological_parameters, "sigma_8"] = sigma_8_MG
+    block[names.cosmological_parameters, "S_8"] = sigma_8_MG*np.sqrt(omega_m/0.3)
+
+    block[names.cosmological_parameters, "sigma_8_LCDM"] = sigma_8
+    block[names.cosmological_parameters, "S_8_LCDM"] = S8_LCDM
 
     return 0
 

pyreact/react.py

@@ -36,16 +36,33 @@ def test_func(self, a):
         r = f(*array_arg(a))
         return r
 
-    def compute_reaction(self, h, n_s, omega_m, omega_b, sigma_8, mg1,
-                               z, k, Pk, is_transfer=False, mass_loop=30,
+    def compute_reaction(self, h, n_s, omega_m, omega_b, sigma_8,
+                               z, k, Pk, 
+                               model="f(R)", fR0=None, Omega_rc=None, 
+                               is_transfer=False, mass_loop=30,
                                verbose=True):
 
         if max(z) > 2.5:
             raise ValueError("ReACT is unstable above z=2.5, try limiting the range of z values.")
 
         if len(z) > 1 and z[0] > z[-1]:
             raise ValueError("The z array needs to be ordered from low to high redshift.")
-
+
+        if model.lower() == "f(r)":
+            reaction_model = 2
+            modified_gravity_param = fR0
+        elif model.lower() == "dgp":
+            reaction_model = 3
+            modified_gravity_param = Omega_rc
+        elif model.lower() == "gr":
+            reaction_model = 1
+            modified_gravity_param = 0.0
+        else:
+            raise ValueError(f"model '{model}' not supported.")
+
+        if modified_gravity_param is None:
+            raise ValueError("fR0 or Omega_rc need to be specified.")
+
         f = self.get_function("compute_reaction")
         f.restype = np.int
         f.argtypes = [*array_ctype(ndim=1, dtype=np.float64), # P(k, z=0)
@@ -59,27 +76,33 @@ def compute_reaction(self, h, n_s, omega_m, omega_b, sigma_8, mg1,
                       ct.POINTER(ct.c_double),     # sigma_8
                       ct.POINTER(ct.c_double),     # modified gravity param
                       ct.POINTER(ct.c_int),        # mass_loop
+                      ct.POINTER(ct.c_int),        # model (1: GR, 2: f(R), 3; DGP)
                       *array_ctype(ndim=2, dtype=np.float64), # reaction (output)
                       *array_ctype(ndim=2, dtype=np.float64), # linear MG power spectrum (output)
+                      np.ctypeslib.ndpointer(ndim=1, dtype=np.float64, flags="C"),     # modified sigma_8 storage variable
                       ct.POINTER(ct.c_int),        # verbose
                      ]
         reaction = np.zeros((len(k), len(z)), dtype=Pk.dtype, order="C")
         p_lin = np.zeros((len(k), len(z)), dtype=Pk.dtype, order="C")
+        sigma8 = np.zeros(1, dtype=Pk.dtype, order="C")
 
         r =   f(*array_arg(np.ascontiguousarray(Pk, dtype=np.float64)),
                 *array_arg(np.ascontiguousarray(k, dtype=np.float64)),
                 *array_arg(np.ascontiguousarray(z[::-1], dtype=np.float64)), # ReACT expect z order from high to low
                 ct.c_bool(is_transfer),
                 ct.c_double(h), ct.c_double(n_s), ct.c_double(omega_m), ct.c_double(omega_b), ct.c_double(sigma_8),
-                ct.c_double(mg1),
+                ct.c_double(modified_gravity_param),
                 ct.c_int(mass_loop),
+                ct.c_int(reaction_model),
                 *array_arg(reaction),
                 *array_arg(p_lin),
+                sigma8,
                 ct.c_int(verbose),
                 )
         if r != 0:
             string_type = ct.c_char * ct.c_int.in_dll(self.lib, "ERROR_MESSAGE_LEN").value
             error_message = string_type.in_dll(self.lib, "error_message").value.decode()
             raise RuntimeError(f"ReACT code terminated with an error: {error_message}")
         # Get into CAMB ordering (z, k), with increasing z
-        return reaction[:,::-1].T, p_lin[:,::-1].T
+        #Tilman: to check output of modsig8 
+        return reaction[:,::-1].T, p_lin[:,::-1].T, sigma8[0]

pyreact/react_wrapper.cpp

@@ -31,7 +31,7 @@ extern "C" {
         #pragma omp critical
         {
             snprintf(error_message, truncate, "Reaction error: %s", msg);
-            printf("%s\n", error_message);
+            fprintf(stderr, "%s\n", error_message);
         }
     }
 
@@ -54,11 +54,13 @@ extern "C" {
                          int* N_z, double* zvals,
                          bool* is_transfer,
                          double* h, double* n_s, double* Omega_m, double* Omega_b, double* sigma_8,
-                         double* mg1, int* mass_loop,
+                         double* mg1, int* mass_loop, int* model,
                          int* N_k_react, int* N_z_react, double* output_react,
                          int* N_k_pl, int* N_z_pl, double* output_pl,
+                         double* modsig8,
                          int* verbose)
     {
+        int status = 0;
         if(*N_k != *N_k_pk || *N_k != *N_k_react || *N_k != *N_k_pl){
             react_error("Inconsistency in k array sizes");
             return 1;
@@ -92,6 +94,8 @@ extern "C" {
             std::cout<<"sigma_8: " << *sigma_8 << "\n";
             std::cout<<"mg1: " << *mg1 << "\n";
             std::cout<<"mass loop: " << *mass_loop << "\n";
+            std::cout<<"model: "  << *model << " (1:GR, 2:f(R), 3:DGP)\n";
+
         }
 
 
@@ -136,13 +140,15 @@ extern "C" {
         /* declare splines */
         Spline mysr,mysp,myreact;
         /*declare variable array*/
-        double vars[6];
+        double vars[7];
         vars[0] = 1.;
         vars[1] = *Omega_m;
         vars[2] = *mg1;//pow(10.0,*mg1); // edit this for new mg param
         vars[3] = 1.;
         vars[4] = 1.;
         vars[5] = *mass_loop;
+
+        int mod = *model;
         // initialise power spectrum normalisation before running 1-loop computations
         iow.initnorm(vars);
 
@@ -154,7 +160,7 @@ extern "C" {
         }
 
         // 1-loop computations at all redshifts @ k0
-        spt.ploop_init(ploopr,ploopp, zvals , *N_z, vars, k0);
+        spt.ploop_init(ploopr,ploopp, zvals , *N_z, vars, mod, k0);
 
         double myscalef[*N_z];
         for(int i = 0; i<*N_z ; i++){
@@ -173,8 +179,18 @@ extern "C" {
             iow.initnorm(vars);
             // Spherical collapse stuff
             /// initialise delta_c(M), a_vir(M), delta_avir(M) and v(M)
-            halo.scol_init(vars);
-            halo.scol_initp(vars);
+            status = halo.scol_init(vars,mod);
+            status |= halo.scol_initp(vars,mod);
+
+            // store modified sigma 8 at z=0
+            if(j == *N_z-1) {
+                modsig8[0] = vars[6];
+            }
+
+            if(status != 0) {
+                react_error("Failed to compute spherical collapse.");
+                return 1;
+            }
 
               // initialise k_star and mathcal{E}
             halo.react_init2(vars,mysr,mysp);
@@ -196,7 +212,7 @@ extern "C" {
 
             for(int i =0; i < *N_k;  i ++) {
                 output_react[i*(*N_z)+j] =  myreact(kvals[i]);
-                output_pl[i*(*N_z)+j] = halo.plinear_cosmosis(kvals[i]);
+                output_pl[i*(*N_z)+j] = pow2(halo.Lin_Grow(kvals[i]))*powerspectrum[i];//halo.plinear_cosmosis(kvals[i]);
                 if(*verbose > 1) {
                     printf(" %e %e %e \n",zvals[j], kvals[i],halo.plinear_cosmosis(kvals[i]));
                 }

reactions/.gitignore

@@ -4,6 +4,7 @@ __pycache__
 Makefile
 config.log
 config.status
+config.h
 libtool
 
 *.o
@@ -13,8 +14,9 @@ libtool
 *.la
 *.lo
 
+include
 lib
 src/.deps
 src/.libs
 
-cosmosis_output
+cosmosis_output