test better

examples/run_lightcone.py

@@ -21,7 +21,7 @@
 
 if preprocess is True:
     z_bins, z_centers, k_bins = p21c.define_grid_modes_redshifts(6., 8 * units.MHz, z_max = 22, k_min = 0.1 / units.Mpc, k_max = 1 / units.Mpc)
-    p21c.Run(output_dir, "Lightcone_rs" + str(lightcone.random_seed) + "_" + run_id + ".h5", z_bins, z_centers, k_bins, False, lightcone = lightcone, verbose = False)
+    p21c.Run(output_dir, "Lightcone_rs" + str(lightcone.random_seed) + "_" + run_id + ".h5", z_bins, z_centers, k_bins, False, lightcone = lightcone, load = False, save = True, verbose = False)
 else:
     lightcone.save(fname = "Lightcone_rs" + str(lightcone.random_seed) + "_" + run_id + ".h5", direc = output_dir)
 

pyproject.toml

@@ -6,7 +6,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "py21cmcast"
-version = "1.0.2"
+version = "1.0.4"
 description = "Fisher forecasts with 21cm experiments"
 readme = "README.md"
 authors = [{ name = "Gaétan Facchinetti", email = "gaetanfacc@gmail.com" }]

src/py21cmcast.egg-info/PKG-INFO

@@ -1,6 +1,6 @@
 Metadata-Version: 2.1
 Name: py21cmcast
-Version: 1.0.2
+Version: 1.0.4
 Summary: Fisher forecasts with 21cm experiments
 Author-email: Gaétan Facchinetti <gaetanfacc@gmail.com>
 License:                     GNU GENERAL PUBLIC LICENSE

src/py21cmcast/__init__.py

@@ -8,6 +8,7 @@
 from .runs import (
     init_runs,
     init_grid_runs,
+    init_random_runs,
     make_config_one_varying_param,
     run_lightcone_from_config,
 )

src/py21cmcast/core.py

@@ -221,7 +221,10 @@ def __init__(self, dir_path : str, lightcone_name : str,  z_bins : list = None,
             # Compute the optical depth to reionization
             if PY21CMFAST:
                 try: 
-                    self._tau_ion = p21f.compute_tau(redshifts=self._z_glob, global_xHI = self.xH_box, user_params=self._user_params, cosmo_params=self._cosmo_params)
+                    if np.all(self._xHIIdb == 0):
+                        self._tau_ion = p21f.compute_tau(redshifts=self._z_glob, global_xHI = self.xH_box, user_params=self._user_params, cosmo_params=self._cosmo_params)
+                    else:
+                        self._tau_ion = p21f.compute_tau(redshifts=self._z_glob, global_xHI = 1-self._xHIIdb, user_params=self._user_params, cosmo_params=self._cosmo_params)
                 except:
                     self._tau_ion = 0.0
                     if self._verbose:
@@ -264,6 +267,12 @@ def _get_power_spectra(self):
             self._k_array = _k_array[0]
 
         except Exception as e:
+
+            self._power_spectrum = np.array([])
+            self._ps_poisson_noise = np.array([])
+            self._k_array = np.array([])
+            self._chunk_indices = np.array([])
+
             if self._verbose:
                 warnings.warn("Impossible to compute the power spectra: \n" + str(e))
 
@@ -276,17 +285,24 @@ def _get_global_quantities(self):
 
         _global_signal = self._lightcone.global_quantities.get('brightness_temp', np.zeros(len(_lc_glob_redshifts), dtype=np.float64))
         _xH_box        = self._lightcone.global_quantities.get('xH_box', np.zeros(len(_lc_glob_redshifts), dtype=np.float64))
+        _xHIIdb        = self._lightcone.global_quantities.get('xHIIdb', np.zeros(len(_lc_glob_redshifts), dtype=np.float64))
         _x_e_box       = self._lightcone.global_quantities.get('x_e_box', np.zeros(len(_lc_glob_redshifts), dtype=np.float64))
         _Ts_box        = self._lightcone.global_quantities.get('Ts_box', np.zeros(len(_lc_glob_redshifts), dtype=np.float64))
         _Tk_box        = self._lightcone.global_quantities.get('Tk_box', np.zeros(len(_lc_glob_redshifts), dtype=np.float64))
 
+        # constrain the value of the ionization fraction between 0 and 1
+        _xHIIdb[ _xHIIdb > 1.0 ] = 1.0
+        _xHIIdb[ _xHIIdb < 0.0 ] = 0.0
+
         self._global_signal = interpolate.interp1d(_lc_glob_redshifts, _global_signal)(self._z_glob)
         self._xH_box        = interpolate.interp1d(_lc_glob_redshifts, _xH_box)(self._z_glob)
+        self._xHIIdb        = interpolate.interp1d(_lc_glob_redshifts, _xHIIdb)(self._z_glob)
         self._x_e_box       = interpolate.interp1d(_lc_glob_redshifts, _x_e_box)(self._z_glob)
         self._Ts_box        = interpolate.interp1d(_lc_glob_redshifts, _Ts_box)(self._z_glob)
         self._Tk_box        = interpolate.interp1d(_lc_glob_redshifts, _Tk_box)(self._z_glob)
 
 
+
     def _save(self):
         """
         ##  Saves all the attributes of the class to be easily reload later if necessary
@@ -301,6 +317,7 @@ def _save(self):
                             global_signal     = self.global_signal,
                             chunk_indices     = self.chunk_indices,
                             xH_box            = self.xH_box,
+                            xHIIdb            = self.xHIIdb,
                             x_e_box           = self.x_e_box,
                             Ts_box            = self.Ts_box,
                             Tk_box            = self.Tk_box,
@@ -344,6 +361,7 @@ def _load(self):
                 self._k_bins            = data['k_bins']  / units.Mpc
                 self._global_signal     = data['global_signal']
                 self._xH_box            = data['xH_box']
+                self._xHIIdb            = data['xHIIdb']
                 self._x_e_box           = data.get('x_e_box', None)
                 self._Ts_box            = data.get('Ts_box', None)
                 self._Tk_box            = data.get('Tk_box', None)
@@ -410,6 +428,10 @@ def xH_box(self):
     def x_e_box(self):
         return self._x_e_box
 
+    @property
+    def xHIIdb(self):
+        return self._xHIIdb
+
     @property
     def Ts_box(self):
         return self._Ts_box

src/py21cmcast/runs.py

@@ -50,6 +50,7 @@
 from py21cmcast import tools as p21c_tools
 import warnings
 import numpy as np
+import random
 
 PY21CMFAST = True
 try:
@@ -328,6 +329,118 @@ def indices_combinations(*arrays, index_combinations=None):
 
 
 
+def init_random_runs(config_file: str, *, n_sample = 1000, random_seed = None, clean_existing_dir: bool = False, verbose:bool = False) -> None:
+
+    config = configparser.ConfigParser(delimiters=':')
+    config.optionxform = str
+
+    config.read(config_file)
+
+    name            = config.get('run', 'name')
+    output_dir      = config.get('run', 'output_dir')
+    cache_dir       = config.get('run', 'cache_dir')
+
+    output_run_dir = output_dir + "/" + name.upper() + "/"
+    existing_dir = p21c_tools.make_directory(output_run_dir, clean_existing_dir = clean_existing_dir)
+
+    if existing_dir is True:
+        print('WARNING: Cannot create a clean new folder because clean_existing_dir is False')
+        return 
+
+    try:
+        lightcone_q       = p21c_tools.read_config_params(config.items('lightcone_quantities'))
+    except configparser.NoSectionError:
+        if verbose: 
+            warnings.warn("No section lightcone_quantities provided")
+        lightcone_q = {}
+
+    try:
+        global_q       = p21c_tools.read_config_params(config.items('global_quantities'))
+    except configparser.NoSectionError:
+        if verbose: 
+            warnings.warn("No section global_quantities provided")
+        global_q = {}    
+
+
+    extra_params  = p21c_tools.read_config_params(config.items('extra_params'), int_type = False)
+    user_params   = p21c_tools.read_config_params(config.items('user_params'))
+    flag_options  = p21c_tools.read_config_params(config.items('flag_options'))
+    astro_params  = p21c_tools.read_config_params(config.items('astro_params'), int_type = False, allow_lists=True)
+    cosmo_params  = p21c_tools.read_config_params(config.items('cosmo_params'), int_type = False, allow_lists=True)
+
+
+    params_astro = []
+    params_cosmo = []
+    draws  = []
+
+    # set the random seed
+    random.seed(random_seed)
+
+    for i in range(0, n_sample): 
+
+        # Initialise the parameters
+        params_astro.append({ k : 0 for k in astro_params.keys()})
+        params_cosmo.append({ k : 0 for k in cosmo_params.keys()})
+        draws.append({ k : -1 for k in (astro_params | cosmo_params).keys()})
+
+        for param_key, param_values in (astro_params | cosmo_params).items():
+
+            # if the user enters a wrong input table
+            assert(len(param_values) < 3), "For random generator, can only define min and max values."
+
+            # parameter to randomly draw between min and max values
+            if len(param_values) == 2:
+                u = random.random()
+                if param_key in astro_params:
+                    params_astro[i][param_key] = (param_values[1] - param_values[0])*u  + param_values[0]
+                if param_key in cosmo_params:
+                    params_cosmo[i][param_key] = (param_values[1] - param_values[0])*u  + param_values[0]
+                draws[i][param_key] = u
+
+            # parameter fixed at a given value
+            if len(param_values) == 1:
+                if param_key in astro_params:
+                    params_astro[i][param_key] = param_values[0]
+                if param_key in cosmo_params:
+                    params_cosmo[i][param_key] = param_values[0]
+
+                # remove the parameter in the dictionnary of draws as nothing is drawn for that parameter
+                draws[i].pop(param_key)
+
+
+    with open(output_run_dir + "/Database.txt", 'a') as file:
+
+        print("# random seed = " + str(random_seed), file=file)
+
+        for i in range(0, n_sample): 
+            p21c_tools.write_config_params(output_run_dir + '/Config_' + str(i) + ".config", name, output_run_dir, cache_dir, 
+                                    lightcone_q, global_q, extra_params, user_params, flag_options, params_astro[i], params_cosmo[i], str(i))
+
+            print(str(i) + ' : ' + str(params_astro[i] | params_cosmo[i]), file=file)
+
+        file.close()
+
+    with open(output_run_dir + "Database.npz", 'wb') as file:
+        np.savez(file, params_astro = params_astro, params_cosmo = params_cosmo, random_seed = random_seed)
+        file.close()
+
+    # Create a file containing the uniform distributions
+    with open(output_run_dir + "/Uniform.txt", 'a') as file:
+
+        print("# random seed = " + str(random_seed), file=file)
+
+        for i in range(0, n_sample): 
+            print(str(i) + ' : ' + str(draws[i]), file=file)
+
+        file.close()
+
+    with open(output_run_dir + "Uniform.npz", 'wb') as file:
+        np.savez(file, draws = draws, random_seed = random_seed)
+        file.close()
+
+    print(str(n_sample) + ' config files generated with random seed : ' +  str(random_seed))
+
+