In rare circumstances requiring a specific distribution of redshifts the...

In rare circumstances requiring a specific distribution of redshifts the integration of the spectrum in some filters was not done correctly, inducing relative errors of ~10¯⁵-10¯⁶.

CHANGELOG.md

@@ -30,6 +30,7 @@
 - In sfhfromfile, if the SFR is made of integers cigale crashed. Now we systematically convert it to float. (Médéric Boquien)
 - The order of the parameters for the analysis modules would change each time a new pcigale.ini was generated. Now the order is fixed. (Médéric Boquien)
 - In the output the sfh.age parameter would correspond to the input value minus 1. Now both values are consistent with one another. (Laure Ciesla & Médéric Boquien) 
+- In rare circumstances requiring a specific distribution of redshifts the integration of the spectrum in some filters was not done correctly, inducing relative errors of ~10¯⁵-10¯⁶. (Médéric Boquien)
 
 ### Optimised
 - Prior to version 0.7.0, we needed to maintain the list of redshifts for all the computed models. Past 0.7.0 we just infer the redshift from a list unique redshifts. This means that we can now discard the list of redshifts for all the models and only keep the list of unique redshifts. This saves ~8 MB of memory for every 10⁶ models. the models should be computed slightly faster but it is in the measurement noise. (Médéric Boquien)

pcigale/sed/__init__.py

@@ -307,7 +307,7 @@ class SED(object):
             # filter one.
             w = np.where((wavelength >= lambda_min) &
                          (wavelength <= lambda_max))
-            wavelength_r = utils.best_grid(wavelength[w], trans_table[0])
+            wavelength_r = utils.best_grid(wavelength[w], trans_table[0], key)
             transmission_r = np.interp(wavelength_r, trans_table[0],
                                        trans_table[1])
 

pcigale/sed/utils.py

@@ -8,6 +8,7 @@ from scipy.constants import c, pi
 
 # Cache of dx for integrate(y,dx) done by flux_trapz
 dx_cache = {}
+best_grid_cache = {}
 
 
 def lambda_to_nu(wavelength):
@@ -233,8 +234,7 @@ def memoise_2var(f):
     return memoise_helper
 
 
-@memoise_2var
-def best_grid(wavelengths1, wavelengths2):
+def best_grid(wavelengths1, wavelengths2, key):
     """
     Return the best wavelength grid to regrid to arrays
 
@@ -247,6 +247,8 @@ def best_grid(wavelengths1, wavelengths2):
     ----------
     wavelengths1, wavelengths2: array of floats
         The wavelength grids to be 'regridded'.
+    key: tuple
+        Key to key the results in cache.
 
     Returns
     -------
@@ -255,10 +257,13 @@ def best_grid(wavelengths1, wavelengths2):
 
     """
 
+    if key in best_grid_cache:
+        return best_grid_cache[key]
     wl = np.concatenate((wavelengths1, wavelengths2))
     wl.sort(kind='mergesort')
     flag = np.ones(len(wl), dtype=bool)
     np.not_equal(wl[1:], wl[:-1], out=flag[1:])
+    best_grid_cache[key] = wl[flag]
     return wl[flag]