Maximise the CPU cache hit rate with computing the χ² by storing contiguously...

Maximise the CPU cache hit rate with computing the χ² by storing contiguously all the models at a given redshift.

CHANGELOG.md

@@ -27,6 +27,7 @@
 ### Optimised
 - Slight speedup of the computation of the likelihood from the χ² (Médéric Boquien)
 - Speedup of the computation of the χ² by ~10% (Médéric Boquien)
+- Thanks to a change in the layout of the models storage in RAM, the computation of the χ² is now massively faster when the run contains multiple redshifts. (Médéric Boquien)
 - The the fritz2006 module should now run faster thanks to an optimisation of the computation of the luminosity of the various AGN components (Médéric Boquien & Guang Yang)
 
 ## 2018.0 (2018-11-06)

pcigale/analysis_modules/pdf_analysis/workers.py

@@ -103,6 +103,9 @@ def sed(idx, midx):
     sed = gbl_warehouse.get_sed(gbl_models.params.modules,
                                 gbl_models.params.from_index(midx))
 
+    # The redshift is the fastest varying variable but we want to store it
+    # as the slowest one so that models at a given redshift are contiguous
+    idx = (idx % gbl_models.nz) * gbl_models.nm + idx // gbl_models.nz
     if 'sfh.age' in sed.info and sed.info['sfh.age'] > sed.info['universe.age']:
         for band in gbl_models.flux:
             gbl_models.flux[band][idx] = np.nan
@@ -117,10 +120,11 @@ def sed(idx, midx):
             gbl_models.extprop[prop][idx] = sed.info[prop]
         for prop in gbl_models.intprop:
             gbl_models.intprop[prop][idx] = sed.info[prop]
+    gbl_models.index[idx] = midx
 
     gbl_counter.inc()
 
-
+@profile
 def analysis(idx, obs):
     """Worker process to analyse the PDF and estimate parameters values and
     store them in an instance of ResultsManager.
@@ -141,8 +145,10 @@ def analysis(idx, obs):
         # work on views of the arrays and not on copies to save on RAM.
         z = np.array(
             gbl_models.conf['sed_modules_params']['redshifting']['redshift'])
-        wz = slice(np.abs(obs.redshift-z).argmin(), None, z.size)
-        corr_dz = compute_corr_dz(z[wz.start], obs)
+        length = gbl_models.nm
+        zidx = np.abs(obs.redshift-z).argmin()
+        wz = slice(zidx * length, (zidx + 1) * length, 1)
+        corr_dz = compute_corr_dz(z[zidx], obs)
     else:  # We do not know the redshift so we use the full grid
         wz = slice(0, None, 1)
         corr_dz = 1.
@@ -205,8 +211,7 @@ def analysis(idx, obs):
         best_idx_z = np.nanargmin(chi2)
         gbl_results.best.chi2[idx] = chi2[best_idx_z]
         gbl_results.best.scaling[idx] = scaling[best_idx_z]
-        gbl_results.best.index[idx] = (wz.start + best_idx_z*wz.step +
-                                       gbl_models.block.start)
+        gbl_results.best.index[idx] = gbl_models.index[wz][best_idx_z]
 
     gbl_counter.inc()
 

pcigale/managers/models.py

@@ -41,11 +41,16 @@ class ModelsManager(object):
             bandnames = set(obs.bands+conf['analysis_params']['bands'])
         else:
             bandnames = obs.bands
+
         size = len(params.blocks[iblock])
+        if conf['parameters_file'] is "":
+            self.nz = len(conf['sed_modules_params']['redshifting']['redshift'])
+            self.nm = size // self.nz
 
         self.flux = {band: SharedArray(size) for band in bandnames}
         self.intprop = {prop: SharedArray(size) for prop in self.intpropnames}
         self.extprop = {prop: SharedArray(size) for prop in self.extpropnames}
+        self.index = SharedArray(size)
 
     def save(self, filename):
         """Save the fluxes and properties of all the models into a table.