Fixed RealSpaceInterface and recent shooting bug

.github/workflows/test_nightly.yml

@@ -93,7 +93,7 @@ jobs:
       run: |
         source venv/bin/activate
         cd main_class/python
-        TEST_LEVEL=3 CLASS_VERBOSE=1 python "$(which nosetests)" -a dump_ini_files test_class.py
+        TEST_LEVEL=3 VALGRIND_MODE=1 CLASS_VERBOSE=1 python "$(which nosetests)" -a dump_ini_files test_class.py
         deactivate
     - name: Valgrind 🤖
       run: |

.github/workflows/test_on_pull_request.yml

@@ -38,7 +38,7 @@ jobs:
       run: |
         source venv/bin/activate
         cd main_class/python
-        OMP_NUM_THREADS=16 COMPARE_OUTPUT_REF=1 TEST_LEVEL=2 nosetests -v -a test_scenario test_class.py --nologcapture --nocapture
+        OMP_NUM_THREADS=16 COMPARE_OUTPUT_REF=1 TEST_LEVEL=3 nosetests -v -a test_scenario test_class.py --nologcapture --nocapture
         deactivate
     - name: Upload plots 📤
       if: success() || failure()

explanatory.ini

@@ -151,6 +151,9 @@ gauge = synchronous
 #      (default: set to 'no')
 #nbody_gauge_transfer_functions = yes
 
+# 4.c) Do you want the source functions for total non-relativistic matter, delta_m and theta_m, and baryon+cdm, delta_cb and theta_cb, to be outputed in the current gauge (the one selected in 4.a), instead of being automatically expressed as a gauge-invariant (GI) quantity, likeL: delta_m^GI=delta_m + 3 a H theta_m/k2, theta_m^GI=theta_m+alpha*k2 (default: no, that is, convert to GI)
+#matter_source_in_current_gauge= no
+
 # 5) Hubble parameter : either 'H0' in km/s/Mpc or 'h' (or 'theta_s_100'), where
 #    the latter is the peak scale parameter defined exactly as 100(ds_dec/da_dec)
 #    with a decoupling time given by maximum of visibility function (quite different

external/RealSpaceInterface/Calc2D/CalculationClass.py

@@ -71,7 +71,7 @@ def getData(self, redshiftindex):
         Valuenew = dict()
         FValue_abs = np.abs(self.FValue)
         _min, _max = FValue_abs.min(), FValue_abs.max()
-        dimensions = (self.endshape[0] / 2, self.endshape[1])
+        dimensions = (self.endshape[0] // 2, self.endshape[1])
         for quantity, FT in FValuenew.items():
             FT_abs = np.abs(FT)
             FT_normalized = cv2.resize(FT_abs, dimensions).ravel()
@@ -98,7 +98,7 @@ def getInitialData(self):
         return Value.ravel(), cv2.resize(
             (np.abs(self.FValue) - np.abs(self.FValue).min()) /
             (np.abs(self.FValue).max() - np.abs(self.FValue).min()),
-            (self.endshape[0] / 2, self.endshape[1])).ravel(), (minimum,
+            (self.endshape[0] // 2, self.endshape[1])).ravel(), (minimum,
                                                                 maximum)
 
     def getTransferData(self, redshiftindex):

external/RealSpaceInterface/Calc2D/Database.py

@@ -14,22 +14,22 @@ def __init__(self, directory, db_file="database.dat"):
         self.db_path = os.path.join(directory, db_file)
         if not os.path.exists(self.db_path):
             logging.info("No database found; Creating one at {}.".format(self.db_path))
-            with open(self.db_path, "w") as f:
+            with open(self.db_path, "wb") as f:
                 pickle.dump(dict(), f)
 
         self.db = self.__read_database()
 
     def __read_database(self):
-        with open(self.db_path) as f:
-            return pickle.load(f)
+        with open(self.db_path, "rb") as f:
+            return pickle.load(f, encoding='latin1')
 
     def __write_database(self):
-        with open(self.db_path, "w") as f:
-            pickle.dump(self.db, f)
+        with open(self.db_path, "wb") as f:
+            pickle.dump(self.db, f, protocol=pickle.HIGHEST_PROTOCOL)
 
     def __create_file(self, data):
         filename = str(uuid.uuid4())
-        with open(os.path.join(self.directory, filename), "w") as f:
+        with open(os.path.join(self.directory, filename), "wb") as f:
             pickle.dump(data, f)
         return filename
 
@@ -40,8 +40,8 @@ def __getitem__(self, key):
         frozen_key = self.__get_frozen_key(key)
         if frozen_key in self.db:
             filename = self.db[frozen_key]
-            with open(os.path.join(self.directory, filename)) as f:
-                return pickle.load(f)
+            with open(os.path.join(self.directory, filename), "rb") as f:
+                return pickle.load(f, encoding='latin1')
         else:
             raise KeyError("No data for key: {}".format(key))
 
@@ -55,4 +55,4 @@ def __contains__(self, key):
         Return whether `self` contains a record
         for the given `key`.
         """
-        return self.__get_frozen_key(key) in self.db
+        return self.__get_frozen_key(key) in self.db

external/RealSpaceInterface/Calc2D/TransferFunction.py

@@ -42,6 +42,7 @@ def ComputeTransferFunctionList(cosmologicalParameters, redshift, kperdecade=200
     class_settings.update({
         "output": "mTk",
         "gauge": "newtonian",
+        "matter_source_in_current_gauge": "yes",
         "evolver": "1",
         "P_k_max_h/Mpc": P_k_max,
         "k_per_decade_for_pk": kperdecade,

external/RealSpaceInterface/README

@@ -1,8 +1,14 @@
+CLASS real space interface
+07.2015: started by Max Beutelspacher
+09.2018: improved by Georgios Samaras and released in class v2.7.0
+10.2025: for v3.2.5: made compatible with python 3, credits Timothy Morton and Mowen Zhao
+
+
 For installation of python packages, run
 
     pip install -r requirements.txt
 
-Launch the application with
+To see the real space interface, go to this directory (external/RealSpaceInterface/) and launch the application with
 
     python tornadoserver.py
 
@@ -15,4 +21,3 @@ Then in any browser open the URL
 Cache files are located in cache/, so to clear the cache, run
 
     rm cache/*
-

external/RealSpaceInterface/tornadoserver.py

@@ -143,8 +143,8 @@ def send_frame(self, redindex):
         self.write_message(json.dumps({'type': 'extrema', 'extrema': extrema}))
         progress = float(redindex) / len(self.calc.redshift)
 
-        real = {quantity: base64.b64encode(data.astype(np.float32)) for quantity, data in Valuenew.iteritems()}
-        transfer = {quantity: base64.b64encode(data.astype(np.float32)) for quantity, data in TransferData.iteritems()}
+        real = {quantity: base64.b64encode(data.astype(np.float32)).decode() for quantity, data in Valuenew.items()}
+        transfer = {quantity: base64.b64encode(data.astype(np.float32)).decode() for quantity, data in TransferData.items()}
         self.write_message(
             json.dumps({
                 'type': 'data',
@@ -166,9 +166,9 @@ def send_initial_state(self):
         extremastring = json.dumps({'type': 'extrema', 'extrema': extrema})
         datastring = json.dumps({
             'type': 'data',
-            'real': base64.b64encode(Value.astype(np.float32)),
+            'real': base64.b64encode(Value.astype(np.float32)).decode(),
             'fourier': [],
-            'transfer': base64.b64encode(TransferData.astype(np.float32)),
+            'transfer': base64.b64encode(TransferData.astype(np.float32)).decode(),
             'k': krange.tolist()
             })
         self.write_message(extremastring)
@@ -212,6 +212,7 @@ def set_cosmological_parameters(self, cosmologicalParameters):
 
             z_of_decoupling = self.calc.z_dec
             frame_of_decoupling = np.argmin(np.abs(z_of_decoupling - self.calc.redshift))
+            frame_of_decoupling = int(frame_of_decoupling) # this is np.int64 object, which json doesn't like
             if self.calc.redshift[frame_of_decoupling] > z_of_decoupling:
                 frame_of_decoupling -= 1
             messages.append({

include/common.h

@@ -15,7 +15,7 @@
 #ifndef __COMMON__
 #define __COMMON__
 
-#define _VERSION_ "v3.2.4"
+#define _VERSION_ "v3.2.5"
 
 /* @cond INCLUDE_WITH_DOXYGEN */
 

include/fourier.h

@@ -356,6 +356,7 @@ extern "C" {
 
   int fourier_get_k_list(
                          struct precision *ppr,
+                         struct primordial *ppm,
                          struct perturbations * ppt,
                          struct fourier * pfo
                          );

include/perturbations.h

@@ -198,6 +198,8 @@ struct perturbations
 
   enum possible_gauges gauge; /**< gauge in which to perform this calculation */
 
+  short has_matter_source_in_current_gauge; /**< whether to keep matter and baryon+CDM sources in current gauge, instead of automatic conversion to gauge-invariant variables */
+
   //@}
 
   /** @name - indices running on modes (scalar, vector, tensor) */

python/test_class.py

@@ -93,6 +93,8 @@
 COMPARE_OUTPUT_GAUGE = bool(int(os.getenv('COMPARE_OUTPUT_GAUGE', '0'))) # Compare synchronous and Newtonian gauge outputs?
 COMPARE_OUTPUT_REF = bool(int(os.getenv('COMPARE_OUTPUT_REF', '0'))) # Compare classy with classyref?
 POWER_ALL = bool(int(os.getenv('POWER_ALL', '0'))) # Combine every extension with each other? (Very slow!)
+VALGRIND_MODE = bool(int(os.getenv('VALGRIND_MODE', '0')))
+
 TEST_LEVEL = int(os.getenv('TEST_LEVEL', '0')) # 0 <= TEST_LEVEL <= 3
 
 if COMPARE_OUTPUT_REF:
@@ -108,6 +110,7 @@
 COMPARE_PK_RELATIVE_ERROR = 1e-2
 COMPARE_PK_RELATIVE_ERROR_GAUGE = 5*1e-2
 COMPARE_PK_ABSOLUTE_ERROR = 1e-20
+COMPARE_SHOOTING_RELATIVE_ERROR = 1e-5
 
 # Dictionary of models to test the wrapper against. Each of these scenario will
 # be run against all the possible output choices (nothing, tCl, mPk, etc...),
@@ -116,6 +119,10 @@
 # against. Indeed, when not specifying a field, CLASS takes the default input.
 CLASS_INPUT = {}
 
+# power = compare ALL settings with each of the entries
+# normal = compare only the power settings with this entry
+# onlyfull = compare only for 'full' output settings (see below)
+
 CLASS_INPUT['Output_spectra'] = (
     [{'output': 'mPk', 'P_k_max_1/Mpc': 2},
      {'output': 'tCl'},
@@ -126,7 +133,8 @@
     'power')
 
 CLASS_INPUT['Nonlinear'] = (
-    [{'non linear': 'halofit'}],
+    [{'non_linear': 'halofit'},
+    {'non_linear': 'hmcode'}],
     'power')
 
 CLASS_INPUT['Lensing'] = (
@@ -148,44 +156,65 @@
     CLASS_INPUT['modes'] = (
         [{'modes': 't'},
         {'modes': 's, t'}],
-        'power')
+        'normal')
 
     CLASS_INPUT['Tensor_method'] = (
         [{'tensor method': 'exact'},
         {'tensor method': 'photons'}],
-        'power')
+        'onlyfull')
 
 if TEST_LEVEL > 2:
     CLASS_INPUT['Isocurvature_modes'] = (
         [{'ic': 'ad,nid,cdi', 'c_ad_cdi': -0.5}],
-        'normal')
+        'onlyfull')
 
     CLASS_INPUT['Scalar_field'] = (
         [{'Omega_scf': 0.1, 'attractor_ic_scf': 'yes',
         'scf_parameters': '10, 0, 0, 0'}],
-        'normal')
+        'onlyfull')
 
     CLASS_INPUT['Inflation'] = (
-        [{'P_k_ini type': 'inflation_V'},
-        {'P_k_ini type': 'inflation_H'}],
-        'normal')
-    # CLASS_INPUT['Inflation'] = (
-    #     [{'P_k_ini type': 'inflation_V'},
-    #     {'P_k_ini type': 'inflation_H'},
-    #     {'P_k_ini type': 'inflation_V_end'}],
-    #     'normal')
+        [{'Pk_ini_type': 'inflation_V'},
+        {'Pk_ini_type': 'inflation_H'}],
+        'onlyfull')
+
+    CLASS_INPUT['Shooting'] = (
+        [{'sigma8': 0.8},
+        {'S8': 0.8},
+        {'theta_s_100': 1.04},
+        {'Neff': 4.0},
+        {'Neff': 2.0},
+        {'Neff': 2.0,'N_ncdm':1,'m_ncdm':0.06},
+        {'Neff': 3.0,'N_ncdm':1,'m_ncdm':0.06},
+        {'Neff': 4.0,'N_ncdm':1,'m_ncdm':0.06},
+        {'Neff': 4.0,'N_ncdm':1,'m_ncdm':0.06, 'theta_s_100':1.04,'sigma8':0.8},
+        {'Neff': 4.0,'N_ncdm':1,'m_ncdm':0.06, 'theta_s_100':1.04,'S8':0.8},
+        {'tau_reio':0.05}
+        ],
+        'onlyfull')
+# TODO :: possibly another one with nonlinear=hmcode in the future, but for now this is more stable
+FULL_SETTINGS = {'output':'mPk mTk vTk tCl pCl lCl sCl nCl','P_k_max_1/Mpc': 2,'lensing':'yes','modes':'s,t'}
 
 if POWER_ALL:
     for k, v in iteritems(CLASS_INPUT):
         models, state = v
         CLASS_INPUT[k] = (models, 'power')
+if VALGRIND_MODE:
+    for k, v in iteritems(CLASS_INPUT):
+        if k == 'Output_spectra':
+            continue
+        models, state = v
+        CLASS_INPUT[k] = (models, 'normal')
 
 INPUTPOWER = []
 INPUTNORMAL = [{}]
+INPUTONLYFULL = [{}]
 for key, value in list(CLASS_INPUT.items()):
     models, state = value
     if state == 'power':
         INPUTPOWER.append([{}]+models)
+    elif state == 'onlyfull':
+        INPUTONLYFULL.extend(models)
     else:
         INPUTNORMAL.extend(models)
 
@@ -199,6 +228,10 @@
                 temp_dict.update(elem)
             DICTARRAY.append(temp_dict)
 
+    for onlyfull in INPUTONLYFULL:
+        temp_dict = onlyfull.copy()
+        temp_dict.update(FULL_SETTINGS)
+        DICTARRAY.append(temp_dict)
 
 TUPLE_ARRAY = []
 for e in DICTARRAY:
@@ -212,7 +245,7 @@ def powerset(iterable):
         itertools.combinations(xs, n) for n in range(1, len(xs)+1))
 
 def custom_name_func(testcase_func, param_num, param):
-    special_keys = ['N_ncdm']
+    special_keys = ['N_ncdm','Pk_ini_type','ic','Neff','sigma8','S8','theta_s_100','tau_reio']
     somekeys = []
     for key in param.args[0]:
         if key in special_keys:
@@ -435,9 +468,9 @@ def has_incompatible_input(self):
             if not has_tensor(self.scenario):
                 should_fail = True
 
-        # If we have specified non linear, we must have some form of
+        # If we have specified non_linear, we must have some form of
         # perturbations output.
-        if 'non linear' in self.scenario:
+        if 'non_linear' in self.scenario:
             if 'output' not in self.scenario:
                 should_fail = True
 
@@ -460,7 +493,7 @@ def has_incompatible_input(self):
 
         # If we use inflation module, we must have scalar modes,
         # tensor modes, no vector modes and we should only have adiabatic IC:
-        if 'P_k_ini type' in self.scenario and self.scenario['P_k_ini type'].find('inflation') != -1:
+        if 'Pk_ini_type' in self.scenario and self.scenario['Pk_ini_type'].find('inflation') != -1:
             if 'modes' not in self.scenario:
                 should_fail = True
             else:
@@ -472,6 +505,8 @@ def has_incompatible_input(self):
                     should_fail = True
             if 'ic' in self.scenario and self.scenario['ic'].find('i') != -1:
                 should_fail = True
+            if 'non_linear' in self.scenario and self.scenario['non_linear'].find('hmcode') != -1:
+                should_fail = True
 
 
         return should_fail
@@ -545,6 +580,36 @@ def compare_output(self, reference, reference_name, candidate, candidate_name, r
                     self.pk_faulty_plot(k, reference_pk, reference_name, candidate_pk, candidate_name, rtol_pk)
                     status_pass = False
 
+        if 'output' in self.scenario and self.scenario['output'].find('mPk') != -1:
+            ref_sigma8 = reference.sigma8()
+            can_sigma8 = candidate.sigma8()
+            ref_sigma = reference.sigma(8, 0, h_units=True)
+            can_sigma = candidate.sigma(8, 0, h_units=True)
+            try:
+                np.testing.assert_allclose([ref_sigma8, ref_sigma, can_sigma8],
+                                           [can_sigma8, can_sigma, can_sigma],
+                                           rtol=rtol_pk,
+                                           atol=COMPARE_PK_ABSOLUTE_ERROR)
+            except AssertionError:
+                status_pass = False
+
+        if 'sigma8' in self.scenario:
+            if not np.isclose(self.scenario['sigma8'], candidate.sigma8(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+            if not np.isclose(candidate.sigma8(), reference.sigma8(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+
+        if 'Neff' in self.scenario:
+            if not np.isclose(self.scenario['Neff'], candidate.Neff(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+            if not np.isclose(candidate.Neff(), reference.Neff(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+
+        if 'theta_s_100' in self.scenario:
+            if not np.isclose(self.scenario['theta_s_100'], candidate.theta_s_100(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
+            if not np.isclose(candidate.theta_s_100(), reference.theta_s_100(), rtol=COMPARE_SHOOTING_RELATIVE_ERROR, atol=0):
+              status_pass = False
         return status_pass
 
     def store_ini_file(self, path):