Speed up the computation of the models using the compiled interp function...

Speed up the computation of the models using the compiled interp function rather than the regular numpy one so we bypass checks on the size and type of the arrays.

CHANGELOG.md

@@ -18,6 +18,7 @@
 
 ### Optimised
 - By default the MKL library created many threads for each for the parallel processes. Not only was this not necessary as a high-level parallelisation already exists, but it generated a strong oversubscription on the CPU and on the RAM. The slowdown was over a factor of ~2 in some cases. Now we mandate KML to use only 1 thread fo each process. (Médéric Boquien & Yannick Roehlly)
+- The generic numpy interpolation function was used. As we are in a well-controlled environment, this generated unnecessary verifications on the type and shape of the arrays. The compiled numpy interpolation function is now used, bypassing those checks. This generates a gain of 5-10% in computing speed for the generation of the models. (Médéric Boquien)
 
 ## 0.10.0 (2016-09-15)
 ### Added

pcigale/sed/__init__.py

@@ -30,6 +30,7 @@ Such SED is characterised by:
 """
 
 import numpy as np
+from numpy.core.multiarray import interp # Compiled version
 from scipy.constants import c, parsec
 
 from . import utils
@@ -308,13 +309,13 @@ class SED(object):
             w = np.where((wavelength >= lambda_min) &
                          (wavelength <= lambda_max))
             wavelength_r = utils.best_grid(wavelength[w], trans_table[0], key)
-            transmission_r = np.interp(wavelength_r, trans_table[0],
-                                       trans_table[1])
+            transmission_r = interp(wavelength_r, trans_table[0],
+                                    trans_table[1])
 
             self.cache_filters[key] = (wavelength_r, transmission_r,
                                        lambda_eff)
 
-        l_lambda_r = np.interp(wavelength_r, wavelength, self.luminosity)
+        l_lambda_r = interp(wavelength_r, wavelength, self.luminosity)
 
         f_lambda = utils.luminosity_to_flux(
             utils.flux_trapz(transmission_r * l_lambda_r, wavelength_r, key),

pcigale/sed/utils.py

@@ -4,6 +4,7 @@
 # Authors: Yannick Roehlly, Médéric Boquien
 
 import numpy as np
+from numpy.core.multiarray import interp # Compiled version
 from scipy.constants import c, pi
 
 # Cache of dx for integrate(y,dx) done by flux_trapz
@@ -357,8 +358,8 @@ def interpolate_lumin(wl, lumin, wl_new, lumin_new):
         # already defined.
         w = np.where((wl_unique > wl[0]) & (wl_unique < wl[-1]))
         for i in range(lumin.shape[0]):
-            lumin_out[i, lumin.shape[1]+w[0]] = np.interp(wl_unique[w], wl,
-                                                          lumin[i, :])
+            lumin_out[i, lumin.shape[1]+w[0]] = interp(wl_unique[w], wl,
+                                                       lumin[i, :])
         wl_best = np.concatenate((wl, wl_unique))
         s = argsort_wl(wl_best)
         wl_best = wl_best[s]
@@ -368,7 +369,7 @@ def interpolate_lumin(wl, lumin, wl_new, lumin_new):
         lumin_out[:-1, :] = lumin
 
     # We interpolate the new component on the new merged wavelength grid.
-    lumin_out[-1, :] = np.interp(wl_best, wl_new, lumin_new, left=0., right=0.)
+    lumin_out[-1, :] = interp(wl_best, wl_new, lumin_new, left=0., right=0.)
 
     return (wl_best, lumin_out)